Feature: support ELF with non-collinear spin (nspin = 4)

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`.
 

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
+)
+

source/source_io/test/write_elf_logic_test.cpp

@@ -0,0 +1,177 @@
+#include "gtest/gtest.h"
+#include <vector>
+#include <cmath>
+
+/************************************************
+ *  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<double> tau_values = {0.01, 0.05, 0.1, 0.5, 1.0};
+    std::vector<double> tau_vw_values = {0.005, 0.02, 0.05, 0.2, 0.5};
+    std::vector<double> 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<double> rho(nspin);
+    std::vector<double> 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<double> rho(nspin);
+    std::vector<double> 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<double> rho = {0.1, -0.02, 0.03, -0.01};
+    std::vector<double> 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();
+}

source/source_io/write_elf.cpp

@@ -75,15 +75,40 @@ 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)
     for (int is = 0; is < nspin; ++is)
     {
         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);
+        }
     }
 }
 }