Fix: SDFT bug after merge

Author: dyzheng

source/src_pw/stress_func_kin.cpp

@@ -51,7 +51,15 @@ void Stress_Func::stress_kin(ModuleBase::matrix& sigma)
 			{
 				for(int ibnd=0;ibnd<GlobalV::NBANDS;ibnd++)
 				{
-					const std::complex<double>* ppsi = &(GlobalC::wf.psi[0](ik, ibnd, 0));
+					const std::complex<double>* ppsi=nullptr;
+					if(GlobalC::wf.psi!=nullptr)
+					{
+						ppsi = &(GlobalC::wf.psi[0](ik, ibnd, 0));
+					}
+					else
+					{
+						ppsi = &(GlobalC::wf.evc[ik](ibnd, 0));
+					}
 					for(int i=0;i<npw;i++)
 					{
 						s_kin[l][m] +=

source/src_pw/stress_func_mgga.cpp

@@ -44,9 +44,18 @@ void Stress_Func::stress_mgga(ModuleBase::matrix& sigma)
 		for (int ibnd = 0; ibnd < GlobalV::NBANDS; ibnd++)
 		{
 			const double w1 = GlobalC::wf.wg(ik, ibnd) / GlobalC::ucell.omega;
+			const std::complex<double>* ppsi=nullptr;
+			if(GlobalC::wf.psi!=nullptr)
+			{
+				ppsi = &(GlobalC::wf.psi[0](ik, ibnd, 0));
+			}
+			else
+			{
+				ppsi = &(GlobalC::wf.evc[ik](ibnd, 0));
+			}
 			for(int ig = 0; ig<npw; ig++)
 			{
-				psi[ig]=GlobalC::wf.psi[0](ik, ibnd, ig);
+				psi[ig] = ppsi[ig];
 			}
 			XC_Functional::grad_wfc(psi, ik, gradwfc, npw);
 

source/src_pw/stress_func_nl.cpp

@@ -60,7 +60,15 @@ void Stress_Func::stress_nl(ModuleBase::matrix& sigma)
 			if(GlobalC::wf.wg(ik, ib) < ModuleBase::threshold_wg) continue;
 			for (int i = 0; i < nkb; i++) 
 			{
-				const std::complex<double>* ppsi = &(GlobalC::wf.psi[0](ik, ib, 0));
+				const std::complex<double>* ppsi=nullptr;
+				if(GlobalC::wf.psi!=nullptr)
+				{
+					ppsi = &(GlobalC::wf.psi[0](ik, ib, 0));
+				}
+				else
+				{
+					ppsi = &(GlobalC::wf.evc[ik](ib, 0));
+				}
 				const std::complex<double>* pvkb = &(GlobalC::ppcell.vkb(i, 0));
                 for (int ig = 0; ig < GlobalC::wf.npw; ig++) 
 				{
@@ -146,7 +154,15 @@ void Stress_Func::stress_nl(ModuleBase::matrix& sigma)
 					if(GlobalC::wf.wg(ik, ib) < ModuleBase::threshold_wg) continue;
 					for (int i=0; i<nkb; i++)
 					{
-						const std::complex<double>* ppsi = &(GlobalC::wf.psi[0](ik, ib, 0));
+						const std::complex<double>* ppsi=nullptr;
+						if(GlobalC::wf.psi!=nullptr)
+						{
+							ppsi = &(GlobalC::wf.psi[0](ik, ib, 0));
+						}
+						else
+						{
+							ppsi = &(GlobalC::wf.evc[ik](ib, 0));
+						}
 						const std::complex<double>* pdbecp_noevc = &(dbecp_noevc(i, 0));
 						for (int ig=0; ig<GlobalC::wf.npw; ig++) 
 						{

source/src_pw/wavefunc.cpp

@@ -128,6 +128,14 @@ void wavefunc::allocate(const int nks)
 		std::cout << " MEMORY FOR PSI (MB)  : " << 
 		ModuleBase::Memory::record("wavefunc","evc",nks2*GlobalV::NBANDS*(prefactor*npwx),"complexmatrix") << std::endl;
 	}
+	else if(GlobalV::CALCULATION.substr(0,3) == "sto")
+	{
+		this->evc = new ModuleBase::ComplexMatrix [nks2];
+		for (int ik = 0; ik < nks2; ik++)
+		{
+			this->evc[ik].create(GlobalV::NBANDS, npwx * GlobalV::NPOL);//added by zhengdy-soc
+		}
+	}
 	else
 	{
 		//initial psi rather than evc
@@ -359,7 +367,14 @@ void wavefunc::wfcinit_k(void)
 			// get the wave functions
 			// by first diagolize PAO
 			// wave functions.
-			this->diago_PAO_in_pw_k(ik, this->psi[0]);
+			if(GlobalV::CALCULATION.substr(0,3) == "sto")
+			{
+				this->diago_PAO_in_pw_k(ik, this->evc[ik]);
+			}
+			else
+			{
+				this->diago_PAO_in_pw_k(ik, this->psi[0]);
+			}
 		}
 #ifdef __LCAO
 		else if(GlobalV::BASIS_TYPE=="lcao_in_pw")
@@ -823,6 +838,15 @@ void wavefunc::init_after_vc(const int nks)
 }
 
 //temporary function for nscf
+void wavefunc::diago_PAO_in_pw_k(const int &ik, ModuleBase::ComplexMatrix &wvf)
+{
+	ModuleBase::TITLE("wavefunc","diago_PAO_in_pw_k");
+
+	GlobalC::hm.hpw.init_k(ik);
+    this->diago_PAO_in_pw_k2(ik, wvf);
+
+	return;
+}
 void wavefunc::diago_PAO_in_pw_k2(const int &ik, ModuleBase::ComplexMatrix &wvf)
 {
 	ModuleBase::TITLE("wavefunc","diago_PAO_in_pw_k2");

source/src_pw/wavefunc.h

@@ -45,6 +45,7 @@ class wavefunc : public WF_atomic
 	// evc: get the initial wave functions from diagnalized the PAO
 	// orbitals first.
 	void diago_PAO_in_pw_k(const int &ik, psi::Psi<std::complex<double>> &wvf);
+    void diago_PAO_in_pw_k(const int &ik, ModuleBase::ComplexMatrix &wvf);
 
 	// used if k dependent staff is ready.
 	void prepare_k(void);