Refactor: Refine ELF Output for Noncollinear Spin (nspin=4) (deepmodeling#6914)

* Fix: Remove the warning for cal_elf in NSPIN=4 case

* Refactor: Update the calculation of ELF for NSPIN=4 case.

* Doc: Update the documentation of out_elf

* Test: Add an integrate test `tests/03_NAO_multik/63_NO_KP_out_elf` for out_elf with NSPIN=4

* Fix: Update the annotation of write_elf

* Fix: Enhance the stability of write_elf for spin=2 case

* Test: update 63_NO_KP_out_elf

* Test: Reduce the cell length of 63_NO_KP_out_elf

Author: sunliang98

docs/advanced/input_files/input-main.md

@@ -2133,12 +2133,12 @@ These variables are used to control the output of properties.
 - **Availability**: Only for Kohn-Sham DFT and Orbital Free DFT.
 - **Description**: Whether to output the electron localization function (ELF) in the folder `OUT.${suffix}`. The files are named as
     - nspin = 1:
-      - ELF.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$;
+      - elf.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$;
     - 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}}}$;
+      - elf1.cube, elf2.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.
+      - elf.cube: ELF for total charge density, ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$
 
   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/read_input_item_output.cpp

@@ -564,10 +564,6 @@ void ReadInput::item_output()
             {
                 ModuleBase::WARNING_QUIT("ReadInput", "ELF is only aviailable for ksdft and ofdft");
             }
-            if (para.input.out_elf[0] > 0 && para.input.nspin == 4)
-            {
-                ModuleBase::WARNING_QUIT("ReadInput", "ELF is not aviailable for nspin = 4");
-            }
         };
         sync_intvec(input.out_elf, 2, 0);
         this->add_item(item);

source/source_io/write_elf.cpp

@@ -18,12 +18,20 @@ void write_elf(
     const UnitCell* ucell_,
     const int& precision)
 {
-    std::vector<std::vector<double>> elf(nspin, std::vector<double>(rho_basis->nrxx, 0.));
+    // For nspin = 4, we only calculate the total ELF using the rho_total and tau_total,
+    // containing in the first channel of rho and tau.
+    // What's more, we have not introduced the U(1) and SU(2) gauge invariance corrections
+    // proposed by Desmarais J K, Vignale G, Bencheikh K, et al. Physical Review Letters, 2024, 133(13): 136401,
+    // where the current density is also included in the ELF calculation.
+
+    int nspin_eff = (nspin == 4) ? 1 : nspin;
+
+    std::vector<std::vector<double>> elf(nspin_eff, std::vector<double>(rho_basis->nrxx, 0.));
     // 1) calculate the kinetic energy density of vW KEDF
-    std::vector<std::vector<double>> tau_vw(nspin, std::vector<double>(rho_basis->nrxx, 0.));
+    std::vector<std::vector<double>> tau_vw(nspin_eff, std::vector<double>(rho_basis->nrxx, 0.));
     std::vector<double> phi(rho_basis->nrxx, 0.); // phi = sqrt(rho)
 
-    for (int is = 0; is < nspin; ++is)
+    for (int is = 0; is < nspin_eff; ++is)
     {
         for (int ir = 0; ir < rho_basis->nrxx; ++ir)
         {
@@ -54,37 +62,34 @@ void write_elf(
     }
 
     // 2) calculate the kinetic energy density of TF KEDF
-    std::vector<std::vector<double>> tau_TF(nspin, std::vector<double>(rho_basis->nrxx, 0.));
+    std::vector<std::vector<double>> tau_TF(nspin_eff, std::vector<double>(rho_basis->nrxx, 0.));
     const double c_tf
         = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0)
           * 2.0; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2)
-    if (nspin == 1)
+    if (nspin == 1 || nspin == 4)
     {
         for (int ir = 0; ir < rho_basis->nrxx; ++ir)
         {
-            tau_TF[0][ir] = c_tf * std::pow(rho[0][ir], 5.0 / 3.0);
-        }
-    }
-    else if (nspin == 2)
-    {
-        for (int is = 0; is < nspin; ++is)
-        {
-            for (int ir = 0; ir < rho_basis->nrxx; ++ir)
+            if (rho[0][ir] > 0.0)
             {
-                tau_TF[is][ir] = 0.5 * c_tf * std::pow(2.0 * rho[is][ir], 5.0 / 3.0);
+                tau_TF[0][ir] = c_tf * std::pow(rho[0][ir], 5.0 / 3.0);
+            }
+            else
+            {
+                tau_TF[0][ir] = 0.0;
             }
         }
     }
-    else if (nspin == 4)
+    else if (nspin == 2)
     {
+        // the spin-scaling law: tau_TF[rho_up, rho_dn] = 1/2 * (tau_TF[2*rho_up] + tau_TF[2*rho_dn])
         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);
+                    tau_TF[is][ir] = 0.5 * c_tf * std::pow(2.0 * rho[is][ir], 5.0 / 3.0);
                 }
                 else
                 {
@@ -96,7 +101,7 @@ void write_elf(
 
     // 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 is = 0; is < nspin_eff; ++is)
     {
         for (int ir = 0; ir < rho_basis->nrxx; ++ir)
         {
@@ -116,7 +121,7 @@ void write_elf(
     double ef_tmp = 0.0;
     int out_fermi = 0;
 
-    if (nspin == 1)
+    if (nspin == 1 || nspin == 4)
     {
         std::string fn = out_dir + "/elf.cube";
 
@@ -174,25 +179,5 @@ void write_elf(
             precision,
             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);
-        }
-    }
 }
 }

tests/03_NAO_multik/63_NO_KP_out_elf/INPUT

@@ -0,0 +1,28 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix			autotest
+calculation     scf
+
+nbands			20
+pseudo_dir	../../PP_ORB
+orbital_dir	../../PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc			5
+scf_thr				1e-7
+scf_nmax			100
+
+nspin           4
+#Parameters (3.Basis)
+basis_type		lcao
+
+#Parameters (4.Smearing)
+smearing_method		gauss
+smearing_sigma			0.002
+
+#Parameters (5.Mixing)
+mixing_type		plain
+mixing_beta		0.7
+mixing_gg0      0.0
+
+out_elf	1

tests/03_NAO_multik/63_NO_KP_out_elf/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 1 1 0 0 0

tests/03_NAO_multik/63_NO_KP_out_elf/README

@@ -0,0 +1 @@
+test out_elf for nspin=4

tests/03_NAO_multik/63_NO_KP_out_elf/STRU

@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+C  12.0  C_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+C_gga_8au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+1.89035917
+
+LATTICE_VECTORS
+4.0 0.0 0.0 #latvec3
+0.0 4.0 0.0
+0.0 0.0 4.0
+
+ATOMIC_POSITIONS
+Direct
+
+C #label
+0 #magnetism
+1 #number of atoms
+0.0 0.0000000000000000  0.0000000000000000 1 1 1
+