Feature : implement stress for mgga nspin = 2 (deepmodeling#3263)

* Fix : remove obsolete error message * Feature : mgga stress for nspin = 2 * add use_libxc directory * fix bug --------- Co-authored-by: wenfei-li <[email protected]> Co-authored-by: Qianrui Liu <[email protected]>
maki49 · Nov 26, 2023 · 7d237b5 · 7d237b5
1 parent 23c7361
commit 7d237b5
Show file tree

Hide file tree

Showing 8 changed files with 118 additions and 17 deletions.
diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp
@@ -200,10 +200,10 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
     {
         ModuleBase::WARNING_QUIT("set_xc_type","meta-GGA has not been implemented for nspin = 4 yet");
     }
-    if((func_type == 3 || func_type == 5) && GlobalV::CAL_STRESS == 1 && GlobalV::NSPIN!=1)
-    {
-        ModuleBase::WARNING_QUIT("set_xc_type","mgga stress not implemented for polarized case yet");
-    }
+    //if((func_type == 3 || func_type == 5) && GlobalV::CAL_STRESS == 1 && GlobalV::NSPIN!=1)
+    //{
+    //    ModuleBase::WARNING_QUIT("set_xc_type","mgga stress not implemented for polarized case yet");
+    //}
 
 #ifndef __EXX
     if(func_type == 4 || func_type == 5)

diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h
@@ -171,13 +171,20 @@ class XC_Functional
 // This file contains wrapper for the mGGA functionals
 // it includes 1 subroutine:
 // 1. tau_xc
+// 2. tau_xc_spin
 
 // NOTE : mGGA is realized through LIBXC
 
 #ifdef USE_LIBXC
 	// mGGA
 	static void tau_xc(const double &rho, const double &grho, const double &atau, double &sxc,
           double &v1xc, double &v2xc, double &v3xc);
+
+	static void tau_xc_spin(double rhoup, double rhodw, 
+			ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+			double tauup, double taudw,
+			double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud,
+			double &v3xcup, double &v3xcdw);
 #endif 
 
 //-------------------

diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp
@@ -290,9 +290,21 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		{
 			if(use_libxc)
 			{
+#ifdef USE_LIBXC
 				double sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud;
-				XC_Functional::gcxc_spin_libxc(rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], 
-					sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud);
+				if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA
+				{
+					double v3xcup, v3xcdw;
+					double atau1 = chr->kin_r[0][ir]/2.0;
+					double atau2 = chr->kin_r[1][ir]/2.0;
+					XC_Functional::tau_xc_spin( rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], 
+						atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw);
+				}
+				else
+				{
+					XC_Functional::gcxc_spin_libxc(rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], 
+						sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud);
+				}
 				if(is_stress)
 				{
 					double tt1[3],tt2[3];
@@ -330,6 +342,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 					local_vtxcgc = local_vtxcgc + ModuleBase::e2 * v1xcdw * ( rhotmp2[ir] - chr->rho_core[ir] * fac );
 					local_etxcgc = local_etxcgc + ModuleBase::e2 * sxc;
 				}
+#endif
 			}
 			else
 			{

diff --git a/source/module_hamilt_general/module_xc/xc_functional_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_wrapper_tauxc.cpp
@@ -40,4 +40,77 @@ void XC_Functional::tau_xc(const double &rho, const double &grho, const double &
 	return;
 }
 
+
+void XC_Functional::tau_xc_spin(double rhoup, double rhodw, 
+        ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+        double tauup, double taudw,
+        double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud,
+        double &v3xcup, double &v3xcdw)
+{
+
+	std::vector<xc_func_type> funcs = init_func(XC_POLARIZED);
+
+    double *rho, *grho, *tau, *v1xc, *v2xc, *v3xc, *sgn, s;
+    sxc = v1xcup = v1xcdw = 0.0;
+    v2xcup = v2xcdw = v2xcud = 0.0;
+    v3xcup = v3xcdw = 0.0;
+    rho = new double[2];
+    grho= new double[3];
+    tau = new double[2];
+    v1xc= new double[2];
+    v2xc= new double[3];
+    v3xc= new double[2];
+    sgn = new double[2];
+
+    rho[0] = rhoup;
+    rho[1] = rhodw;
+    grho[0] = gdr1.norm2();
+    grho[1] = gdr1 * gdr2;
+    grho[2] = gdr2.norm2();
+    tau[0] = tauup;
+    tau[1] = taudw;
+
+    double lapl[2] = {0.0,0.0};
+    double vlapl[2];
+
+    const double rho_threshold = 1E-6;
+    const double grho_threshold = 1E-10;
+
+    for(xc_func_type &func : funcs)
+    {
+        if( func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA)
+        {
+            sgn[0] = sgn[1] = 1.0;
+            // call Libxc function: xc_mgga_exc_vxc
+            xc_mgga_exc_vxc( &func, 1, rho, grho, lapl, tau, &s, v1xc, v2xc, vlapl, v3xc);
+            if(func.info->kind==XC_CORRELATION)
+            {
+                if ( rho[0]<rho_threshold || sqrt(std::abs(grho[0]))<grho_threshold )
+                    sgn[0] = 0.0;
+                if ( rho[1]<rho_threshold || sqrt(std::abs(grho[2]))<grho_threshold )
+                    sgn[1] = 0.0;
+            }
+            sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]);
+            v1xcup += v1xc[0] * sgn[0];
+            v1xcdw += v1xc[1] * sgn[1];
+            v2xcup += 2.0 * v2xc[0] * sgn[0];
+            v2xcud += v2xc[1] * sgn[0] * sgn[1];
+            v2xcdw += 2.0 * v2xc[2] * sgn[1];
+            v3xcup += v3xc[0] * sgn[0];
+            v3xcdw += v3xc[1] * sgn[1];
+        }
+    }
+
+    delete[] grho;
+    delete[] rho;
+    delete[] tau;
+    delete[] v1xc;
+    delete[] v2xc;
+    delete[] v3xc;
+    delete[] sgn;
+    finish_func(funcs);
+
+	return;
+}
+
 #endif
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp
@@ -85,6 +85,14 @@ void Stress_Func<FPTYPE, Device>::stress_mgga(ModuleBase::matrix& sigma,
 
     delete[] gradwfc;
 
+    for(int ix = 0; ix < 3; ix++)
+    {
+        for(int iy = 0; iy < 3; iy++)
+        {
+            sigma_mgga[ix][iy] = 0.0;
+        }
+    }
+
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
         for (int ix = 0; ix < 3; ix++)
@@ -94,7 +102,6 @@ void Stress_Func<FPTYPE, Device>::stress_mgga(ModuleBase::matrix& sigma,
                 FPTYPE delta = 0.0;
                 if (ix == iy)
                     delta = 1.0;
-                sigma_mgga[ix][iy] = 0.0;
                 for (int ir = 0; ir < nrxx; ir++)
                 {
                     FPTYPE x

diff --git a/tests/integrate/116_PW_scan_Si2/result.ref b/tests/integrate/116_PW_scan_Si2/result.ref
@@ -1,5 +1,5 @@
-etotref -204.1114225159662396
-etotperatomref -102.0557112580
-totalforceref 16.621252
-totalstressref 1353.880229
-totaltimeref +0.88323
+etotref -204.0730437254555625
+etotperatomref -102.0365218627
+totalforceref 16.672700
+totalstressref 1356.938302
+totaltimeref 1.57
diff --git a/tests/integrate/116_PW_scan_Si2_nspin2/INPUT b/tests/integrate/116_PW_scan_Si2_nspin2/INPUT
@@ -26,3 +26,4 @@ mixing_type		broyden
 mixing_beta		0.7
 
 cal_force 1
+cal_stress 1
diff --git a/tests/integrate/116_PW_scan_Si2_nspin2/result.ref b/tests/integrate/116_PW_scan_Si2_nspin2/result.ref
@@ -1,5 +1,5 @@
-etotref -204.1131747676821249
-etotperatomref -102.0565873838
-totalforceref 16.565288
-totaltimeref 
-3.07882
+etotref -204.0730158607786393
+etotperatomref -102.0365079304
+totalforceref 16.660110
+totalstressref 1357.594842
+totaltimeref 2.39