From a5c79956019cfe59af6c004ce9c224d96378dc9a Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 27 Jan 2026 00:31:01 +0800 Subject: [PATCH 1/2] Feature: support ELF with non-collinear spin (nspin = 4) --- docs/advanced/input_files/input-main.md | 2 + source/source_io/write_elf.cpp | 51 +++++++++++++++++++++++-- 2 files changed, 50 insertions(+), 3 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index c022f0520d..c0426da0eb 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2106,6 +2106,8 @@ These variables are used to control the output of properties. - nspin = 2: - ELF_SPIN1.cube, ELF_SPIN2.cube: ${\rm{ELF}}_\sigma = \frac{1}{1+\chi_\sigma^2}$, $\chi_\sigma = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i,\sigma}|^2} - \frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}$; - ELF.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i,\sigma}{f_i |\nabla\psi_{i,\sigma}|^2} - \sum_{\sigma}{\frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}}{\sum_{\sigma}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}}$; + - nspin = 4 (noncollinear): + - ELF0.cube, ELF1.cube, ELF2.cube, ELF3.cube: ELF for each Pauli matrix component. Component 0 represents the total charge density, while components 1-3 represent the magnetization in x, y, and z directions respectively. Each component is calculated as ${\rm{ELF}}_i = \frac{1}{1+\chi_i^2}$, where $\chi_i = \frac{\tau_i - \tau_{vW,i}}{\tau_{TF,i}}$, with $\tau_i$ being the kinetic energy density, $\tau_{vW,i} = \frac{1}{2}|\nabla\sqrt{\rho_i}|^2$ the von Weizsäcker kinetic energy density, and $\tau_{TF,i} = \frac{3}{10}(3\pi^2)^{2/3}\rho_i^{5/3}$ the Thomas-Fermi kinetic energy density. The second integer controls the precision of the kinetic energy density output, if not given, will use `3` as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use `10`. diff --git a/source/source_io/write_elf.cpp b/source/source_io/write_elf.cpp index 6c34a08773..bcb694609c 100644 --- a/source/source_io/write_elf.cpp +++ b/source/source_io/write_elf.cpp @@ -75,6 +75,24 @@ void write_elf( } } } + else if (nspin == 4) + { + for (int is = 0; is < nspin; ++is) + { + for (int ir = 0; ir < rho_basis->nrxx; ++ir) + { + // Handle negative densities for numerical stability + if (rho[is][ir] > 0.0) + { + tau_TF[is][ir] = c_tf * std::pow(rho[is][ir], 5.0 / 3.0); + } + else + { + tau_TF[is][ir] = 0.0; + } + } + } + } // 3) calculate the enhancement factor F = (tau_KS - tau_vw) / tau_TF, and then ELF = 1 / (1 + F^2) double eps = 1.0e-5; // suppress the numerical instability in LCAO (Ref: Acta Phys. -Chim. Sin. 2011, 27(12), 2786-2792. doi: 10.3866/PKU.WHXB20112786) @@ -82,8 +100,15 @@ void write_elf( { for (int ir = 0; ir < rho_basis->nrxx; ++ir) { - elf[is][ir] = (tau[is][ir] - tau_vw[is][ir] + eps) / tau_TF[is][ir]; - elf[is][ir] = 1. / (1. + elf[is][ir] * elf[is][ir]); + if (tau_TF[is][ir] > 1.0e-12) + { + elf[is][ir] = (tau[is][ir] - tau_vw[is][ir] + eps) / tau_TF[is][ir]; + elf[is][ir] = 1. / (1. + elf[is][ir] * elf[is][ir]); + } + else + { + elf[is][ir] = 0.0; + } } } @@ -147,7 +172,27 @@ void write_elf( ef_tmp, ucell_, precision, - out_fermi); + out_fermi); + } + else if (nspin == 4) + { + for (int is = 0; is < nspin; ++is) + { + std::string fn = out_dir + "/elf" + std::to_string(is) + ".cube"; + + int ispin = is; + + ModuleIO::write_vdata_palgrid(pgrid, + elf[is].data(), + ispin, + nspin, + istep_in, + fn, + ef_tmp, + ucell_, + precision, + out_fermi); + } } } } From dca22ca67661f8777cf9f17791cec149c01a3684 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 27 Jan 2026 10:45:34 +0800 Subject: [PATCH 2/2] fix: UT for write_elf --- source/source_io/test/CMakeLists.txt | 6 + .../source_io/test/write_elf_logic_test.cpp | 177 ++++++++++++++++++ 2 files changed, 183 insertions(+) create mode 100644 source/source_io/test/write_elf_logic_test.cpp diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index 407eadecec..c4e01f0a0e 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -278,3 +278,9 @@ AddTest( ../../source_basis/module_nao/radial_collection.cpp ) endif() + +AddTest( + TARGET MODULE_IO_write_elf_logic_test + SOURCES write_elf_logic_test.cpp +) + diff --git a/source/source_io/test/write_elf_logic_test.cpp b/source/source_io/test/write_elf_logic_test.cpp new file mode 100644 index 0000000000..25269172f7 --- /dev/null +++ b/source/source_io/test/write_elf_logic_test.cpp @@ -0,0 +1,177 @@ +#include "gtest/gtest.h" +#include +#include + +/************************************************ + * unit test of write_elf logic + ***********************************************/ + +/** + * This test verifies the ELF calculation logic for nspin=4 + * by testing the key formulas directly without file I/O. + */ + +class ElfLogicTest : public ::testing::Test +{ +protected: + // Test the Thomas-Fermi kinetic energy density calculation + double calculate_tau_TF(double rho) { + const double c_tf = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0) * 2.0; + if (rho > 0.0) { + return c_tf * std::pow(rho, 5.0 / 3.0); + } else { + return 0.0; + } + } + + // Test the ELF calculation + double calculate_elf(double tau, double tau_vw, double tau_TF) { + const double eps = 1.0e-5; + if (tau_TF > 1.0e-12) { + double chi = (tau - tau_vw + eps) / tau_TF; + return 1.0 / (1.0 + chi * chi); + } else { + return 0.0; + } + } +}; + +TEST_F(ElfLogicTest, ThomasFermiPositiveDensity) +{ + // Test with positive density + double rho = 0.1; + double tau_TF = calculate_tau_TF(rho); + + EXPECT_GT(tau_TF, 0.0); + EXPECT_LT(tau_TF, 1.0); // Should be reasonable value +} + +TEST_F(ElfLogicTest, ThomasFermiNegativeDensity) +{ + // Test with negative density (for magnetization components) + double rho = -0.02; + double tau_TF = calculate_tau_TF(rho); + + EXPECT_EQ(tau_TF, 0.0); // Should return 0 for negative density +} + +TEST_F(ElfLogicTest, ThomasFermiZeroDensity) +{ + // Test with zero density + double rho = 0.0; + double tau_TF = calculate_tau_TF(rho); + + EXPECT_EQ(tau_TF, 0.0); +} + +TEST_F(ElfLogicTest, ElfCalculationNormal) +{ + // Test ELF calculation with normal values + double tau = 0.05; + double tau_vw = 0.02; + double tau_TF = 0.03; + + double elf = calculate_elf(tau, tau_vw, tau_TF); + + EXPECT_GE(elf, 0.0); + EXPECT_LE(elf, 1.0); +} + +TEST_F(ElfLogicTest, ElfCalculationSmallTauTF) +{ + // Test ELF calculation with very small tau_TF + double tau = 0.05; + double tau_vw = 0.02; + double tau_TF = 1.0e-15; // Very small + + double elf = calculate_elf(tau, tau_vw, tau_TF); + + EXPECT_EQ(elf, 0.0); // Should return 0 for very small tau_TF +} + +TEST_F(ElfLogicTest, ElfCalculationZeroTauTF) +{ + // Test ELF calculation with zero tau_TF + double tau = 0.05; + double tau_vw = 0.02; + double tau_TF = 0.0; + + double elf = calculate_elf(tau, tau_vw, tau_TF); + + EXPECT_EQ(elf, 0.0); // Should return 0 for zero tau_TF +} + +TEST_F(ElfLogicTest, ElfValueRange) +{ + // Test that ELF is always in [0, 1] for various inputs + std::vector tau_values = {0.01, 0.05, 0.1, 0.5, 1.0}; + std::vector tau_vw_values = {0.005, 0.02, 0.05, 0.2, 0.5}; + std::vector tau_TF_values = {0.01, 0.03, 0.08, 0.3, 0.8}; + + for (double tau : tau_values) { + for (double tau_vw : tau_vw_values) { + for (double tau_TF : tau_TF_values) { + double elf = calculate_elf(tau, tau_vw, tau_TF); + EXPECT_GE(elf, 0.0) << "ELF should be >= 0"; + EXPECT_LE(elf, 1.0) << "ELF should be <= 1"; + } + } + } +} + +TEST_F(ElfLogicTest, Nspin4ComponentHandling) +{ + // Test that we can handle 4 components independently + int nspin = 4; + std::vector rho(nspin); + std::vector tau_TF(nspin); + + // Component 0: total charge (positive) + rho[0] = 0.1; + tau_TF[0] = calculate_tau_TF(rho[0]); + EXPECT_GT(tau_TF[0], 0.0); + + // Components 1-3: magnetization (can be negative) + for (int i = 1; i < nspin; ++i) { + rho[i] = -0.01 * i; // Negative + tau_TF[i] = calculate_tau_TF(rho[i]); + EXPECT_EQ(tau_TF[i], 0.0); // Should be 0 for negative density + } +} + +TEST_F(ElfLogicTest, Nspin4AllPositive) +{ + // Test with all positive densities + int nspin = 4; + std::vector rho(nspin); + std::vector tau_TF(nspin); + + for (int i = 0; i < nspin; ++i) { + rho[i] = 0.05 + 0.01 * i; + tau_TF[i] = calculate_tau_TF(rho[i]); + EXPECT_GT(tau_TF[i], 0.0); + } +} + +TEST_F(ElfLogicTest, Nspin4MixedSigns) +{ + // Test with mixed positive and negative densities + int nspin = 4; + std::vector rho = {0.1, -0.02, 0.03, -0.01}; + std::vector tau_TF(nspin); + + for (int i = 0; i < nspin; ++i) { + tau_TF[i] = calculate_tau_TF(rho[i]); + if (rho[i] > 0.0) { + EXPECT_GT(tau_TF[i], 0.0); + } else { + EXPECT_EQ(tau_TF[i], 0.0); + } + } +} + +int main(int argc, char** argv) +{ + testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +}