change ucell in pwdft/stress_func_mgga.cpp

A-006 · A-006 · commit a9626e7b926c · 2024-12-01T11:23:23.000+08:00
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -1008,7 +1008,7 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
     // local pseudopotential stress:
     // use charge density; plane wave; local pseudopotential;
     //--------------------------------------------------------
-    sc_pw.stress_loc(sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr);
+    sc_pw.stress_loc(ucell,sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr);
 
     //--------------------------------------------------------
     // hartree term
diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp
@@ -76,7 +76,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
     stress_gga(ucell,sigmaxc, this->rhopw, pelec->charge);
 
     // local contribution
-    stress_loc(sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge);
+    stress_loc(ucell,sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge);
 
     // nlcc
     stress_cc(sigmaxcc, this->rhopw, p_sf, true, locpp.numeric, pelec->charge);
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h
@@ -84,7 +84,8 @@ class Stress_Func
                     const bool is_pw); // ewald part in PW or LCAO basis
 
     // 4) the stress from the local pseudopotentials
-    void stress_loc(ModuleBase::matrix& sigma,
+    void stress_loc(const UnitCell& ucell,
+                    ModuleBase::matrix& sigma,
                     ModulePW::PW_Basis* rho_basis,
                     const ModuleBase::matrix& vloc,
                     const Structure_Factor* p_sf,
@@ -105,7 +106,8 @@ class Stress_Func
      *        D V(g^2) / D g^2 = 4pi e^2/omegai /G^4
      *
      */
-    void dvloc_coulomb(const FPTYPE& zp,
+    void dvloc_coulomb(const UnitCell& ucell,
+                       const FPTYPE& zp,
                        FPTYPE* dvloc,
                        ModulePW::PW_Basis* rho_basis); // used in local pseudopotential stress
 
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp
@@ -8,7 +8,8 @@
 
 //calculate local pseudopotential stress in PW or VL_dVL stress in LCAO
 template <typename FPTYPE, typename Device>
-void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
+void Stress_Func<FPTYPE, Device>::stress_loc(const UnitCell& ucell,
+											 ModuleBase::matrix& sigma,
                                              ModulePW::PW_Basis* rho_basis,
 											 const ModuleBase::matrix& vloc,
                                              const Structure_Factor* p_sf,
@@ -69,7 +70,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
 	if(is_pw)
 	{
 #pragma omp parallel for collapse(2) reduction(+:evloc)
-		for (int it=0; it<GlobalC::ucell.ntype; it++)
+		for (int it=0; it<ucell.ntype; it++)
 		{
 			for (int ig=0; ig<rho_basis->npw; ig++)
 			{
@@ -83,15 +84,15 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
             }
 		}
     }
-    for (int it = 0; it < GlobalC::ucell.ntype; ++it)
+    for (int it = 0; it < ucell.ntype; ++it)
     {
-        const Atom* atom = &GlobalC::ucell.atoms[it];
+        const Atom* atom = &ucell.atoms[it];
 		if(atom->coulomb_potential)
 		{
 		//
 		// special case: pseudopotential is coulomb 1/r potential
 		//
-			this->dvloc_coulomb (atom->ncpp.zv, dvloc.data(), rho_basis);
+			this->dvloc_coulomb (ucell,atom->ncpp.zv, dvloc.data(), rho_basis);
 		//
 		}
 		else
@@ -100,7 +101,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
 		// normal case: dvloc contains dV_loc(G)/dG
 		//
 			this->dvloc_of_g ( atom->ncpp.msh, atom->ncpp.rab.data(), atom->ncpp.r.data(),
-					atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, GlobalC::ucell);
+					atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, ucell);
 		//
 		}
 #ifndef _OPENMP
@@ -119,7 +120,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
 				{
                     local_sigma(l, m) = local_sigma(l, m)
                                         + (conj(aux[ig]) * p_sf->strucFac(it, ig)).real() * 2.0
-                                              * dvloc[rho_basis->ig2igg[ig]] * GlobalC::ucell.tpiba2
+                                              * dvloc[rho_basis->ig2igg[ig]] * ucell.tpiba2
                                               * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] * fact;
                 }
 			}
@@ -277,7 +278,10 @@ const UnitCell& ucell_in
 }
 
 template <typename FPTYPE, typename Device>
-void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc, ModulePW::PW_Basis* rho_basis)
+void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const UnitCell& ucell, 
+												const FPTYPE& zp, 
+												FPTYPE* dvloc, 
+												ModulePW::PW_Basis* rho_basis)
 {
     int igl0;
 	// start from |G|=0 or not.
@@ -295,8 +299,8 @@ void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc,
 #endif
     for (int i = igl0; i < rho_basis->ngg; i++)
     {
-        dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / GlobalC::ucell.omega
-                   / pow((GlobalC::ucell.tpiba2 * rho_basis->gg_uniq[i]), 2);
+        dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / ucell.omega
+                   / pow((ucell.tpiba2 * rho_basis->gg_uniq[i]), 2);
     }
 
     return;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp
@@ -89,7 +89,7 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
     }
 
     // local contribution
-    this->stress_loc(sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge);
+    this->stress_loc(ucell,sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge);
 
     // nlcc
     this->stress_cc(sigmaxcc, rho_basis, p_sf, 1, nlpp.numeric, pelec->charge);
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
@@ -52,7 +52,7 @@ void Sto_Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
     this->stress_gga(ucell_in,sigmaxc, rho_basis, chr);
 
     // local contribution
-    this->stress_loc(sigmaloc, rho_basis, nlpp->vloc, p_sf, true, chr);
+    this->stress_loc(ucell_in,sigmaloc, rho_basis, nlpp->vloc, p_sf, true, chr);
 
     // nlcc
     this->stress_cc(sigmaxcc, rho_basis, p_sf, true, nlpp->numeric, chr);

Original file line number	Diff line number	Diff line change
`@@ -8,7 +8,8 @@`
`8`	`8`
`9`	`9`	`//calculate local pseudopotential stress in PW or VL_dVL stress in LCAO`
`10`	`10`	`template <typename FPTYPE, typename Device>`
`11`		`-void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,`
	`11`	`+void Stress_Func<FPTYPE, Device>::stress_loc(const UnitCell& ucell,`
	`12`	`+ ModuleBase::matrix& sigma,`
`12`	`13`	`ModulePW::PW_Basis* rho_basis,`
`13`	`14`	`const ModuleBase::matrix& vloc,`
`14`	`15`	`const Structure_Factor* p_sf,`
`@@ -69,7 +70,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,`
`69`	`70`	`if(is_pw)`
`70`	`71`	`{`
`71`	`72`	`#pragma omp parallel for collapse(2) reduction(+:evloc)`
`72`		`- for (int it=0; it<GlobalC::ucell.ntype; it++)`
	`73`	`+ for (int it=0; it<ucell.ntype; it++)`
`73`	`74`	`{`
`74`	`75`	`for (int ig=0; ig<rho_basis->npw; ig++)`
`75`	`76`	`{`
`@@ -83,15 +84,15 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,`
`83`	`84`	`}`
`84`	`85`	`}`
`85`	`86`	`}`
`86`		`- for (int it = 0; it < GlobalC::ucell.ntype; ++it)`
	`87`	`+ for (int it = 0; it < ucell.ntype; ++it)`
`87`	`88`	`{`
`88`		`- const Atom* atom = &GlobalC::ucell.atoms[it];`
	`89`	`+ const Atom* atom = &ucell.atoms[it];`
`89`	`90`	`if(atom->coulomb_potential)`
`90`	`91`	`{`
`91`	`92`	`//`
`92`	`93`	`// special case: pseudopotential is coulomb 1/r potential`
`93`	`94`	`//`
`94`		`- this->dvloc_coulomb (atom->ncpp.zv, dvloc.data(), rho_basis);`
	`95`	`+ this->dvloc_coulomb (ucell,atom->ncpp.zv, dvloc.data(), rho_basis);`
`95`	`96`	`//`
`96`	`97`	`}`
`97`	`98`	`else`
`@@ -100,7 +101,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,`
`100`	`101`	`// normal case: dvloc contains dV_loc(G)/dG`
`101`	`102`	`//`
`102`	`103`	`this->dvloc_of_g ( atom->ncpp.msh, atom->ncpp.rab.data(), atom->ncpp.r.data(),`
`103`		`- atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, GlobalC::ucell);`
	`104`	`+ atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, ucell);`
`104`	`105`	`//`
`105`	`106`	`}`
`106`	`107`	`#ifndef _OPENMP`
`@@ -119,7 +120,7 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,`
`119`	`120`	`{`
`120`	`121`	`local_sigma(l, m) = local_sigma(l, m)`
`121`	`122`	`+ (conj(aux[ig]) * p_sf->strucFac(it, ig)).real() * 2.0`
`122`		`- * dvloc[rho_basis->ig2igg[ig]] * GlobalC::ucell.tpiba2`
	`123`	`+ * dvloc[rho_basis->ig2igg[ig]] * ucell.tpiba2`
`123`	`124`	`* rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] * fact;`
`124`	`125`	`}`
`125`	`126`	`}`
`@@ -277,7 +278,10 @@ const UnitCell& ucell_in`
`277`	`278`	`}`
`278`	`279`
`279`	`280`	`template <typename FPTYPE, typename Device>`
`280`		`-void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc, ModulePW::PW_Basis* rho_basis)`
	`281`	`+void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const UnitCell& ucell,`
	`282`	`+ const FPTYPE& zp,`
	`283`	`+ FPTYPE* dvloc,`
	`284`	`+ ModulePW::PW_Basis* rho_basis)`
`281`	`285`	`{`
`282`	`286`	`int igl0;`
`283`	`287`	`// start from \|G\|=0 or not.`
`@@ -295,8 +299,8 @@ void Stress_Func<FPTYPE, Device>::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc,`
`295`	`299`	`#endif`
`296`	`300`	`for (int i = igl0; i < rho_basis->ngg; i++)`
`297`	`301`	`{`
`298`		`- dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / GlobalC::ucell.omega`
`299`		`- / pow((GlobalC::ucell.tpiba2 * rho_basis->gg_uniq[i]), 2);`
	`302`	`+ dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / ucell.omega`
	`303`	`+ / pow((ucell.tpiba2 * rho_basis->gg_uniq[i]), 2);`
`300`	`304`	`}`
`301`	`305`
`302`	`306`	`return;`
Original file line number	Diff line number	Diff line change
`@@ -89,7 +89,7 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,`
`89`	`89`	`}`
`90`	`90`
`91`	`91`	`// local contribution`
`92`		`- this->stress_loc(sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge);`
	`92`	`+ this->stress_loc(ucell,sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge);`
`93`	`93`
`94`	`94`	`// nlcc`
`95`	`95`	`this->stress_cc(sigmaxcc, rho_basis, p_sf, 1, nlpp.numeric, pelec->charge);`