Fix: Enhance the computational stability of ELF.

Update corresponding integrate tests.

Author: sunliang98

source/module_io/write_elf.cpp

@@ -22,27 +22,34 @@ void write_elf(
     std::vector<std::vector<double>> elf(nspin, 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<double> phi(rho_basis->nrxx, 0.); // phi = sqrt(rho)
+
     for (int is = 0; is < nspin; ++is)
     {
-        std::vector<std::vector<double>> gradient_rho(3, std::vector<double>(rho_basis->nrxx, 0.));
-
-        std::vector<std::complex<double>> recip_rho(rho_basis->npw, 0.0);
-        std::vector<std::complex<double>> recip_gradient_rho(rho_basis->npw, 0.0);
-        rho_basis->real2recip(rho[is], recip_rho.data());
+        for (int ir = 0; ir < rho_basis->nrxx; ++ir)
+        {
+            phi[ir] = std::sqrt(std::abs(rho[is][ir]));
+        }
+        
+        std::vector<std::vector<double>> gradient_phi(3, std::vector<double>(rho_basis->nrxx, 0.));
+        std::vector<std::complex<double>> recip_phi(rho_basis->npw, 0.0);
+        std::vector<std::complex<double>> recip_gradient_phi(rho_basis->npw, 0.0);
+        
+        rho_basis->real2recip(phi.data(), recip_phi.data());
 
         std::complex<double> img(0.0, 1.0);
         for (int j = 0; j < 3; ++j)
         {
             for (int ip = 0; ip < rho_basis->npw; ++ip)
             {
-                recip_gradient_rho[ip] = img * rho_basis->gcar[ip][j] * recip_rho[ip] * rho_basis->tpiba;
+                recip_gradient_phi[ip] = img * rho_basis->gcar[ip][j] * recip_phi[ip] * rho_basis->tpiba;
             }
 
-            rho_basis->recip2real(recip_gradient_rho.data(), gradient_rho[j].data());
+            rho_basis->recip2real(recip_gradient_phi.data(), gradient_phi[j].data());
 
             for (int ir = 0; ir < rho_basis->nrxx; ++ir)
             {
-                tau_vw[is][ir] += gradient_rho[j][ir] * gradient_rho[j][ir] / (8. * rho[is][ir]) * 2.0; // convert Ha to Ry.
+                tau_vw[is][ir] += gradient_phi[j][ir] * gradient_phi[j][ir] / 2. * 2.; // convert Ha to Ry.
             }
         }
     }

tests/integrate/119_PW_out_elf/INPUT

@@ -8,7 +8,7 @@ symmetry		1
 pseudo_dir	../../PP_ORB
 
 #Parameters (2.Iteration)
-ecutwfc			5
+ecutwfc			10
 scf_thr				1e-9
 scf_nmax			100
 

tests/integrate/119_PW_out_elf/refELF.cube

@@ -1,7 +1,7 @@
-etotref -194.6234112615314
-etotperatomref -97.3117056308
+etotref -197.140564442101
+etotperatomref -98.5702822211
 ComparePot1_pass 0
 pointgroupref T_d
 spacegroupref O_h
 nksibzref 1
-totaltimeref 0.25
+totaltimeref 0.34

tests/integrate/119_PW_out_elf/result.ref

@@ -10,7 +10,7 @@ orbital_dir	../../PP_ORB
 out_elf			1
 
 #Parameters (2.Iteration)
-ecutwfc			5
+ecutwfc			20
 scf_thr				1e-6
 scf_nmax			100
 

tests/integrate/219_NO_out_elf/INPUT

@@ -1,7 +1,7 @@
-etotref -196.6254053459516
-etotperatomref -98.3127026730
+etotref -196.6219104988249
+etotperatomref -98.3109552494
 ComparePot1_pass 0
 pointgroupref T_d
 spacegroupref O_h
 nksibzref 1
-totaltimeref 1.14
+totaltimeref 2.58