Merge pull request #920 from Qianruipku/SDFT

refactor: finish forces and stress for SDFT

Author: mohanchen

source/Makefile.Objects

@@ -61,6 +61,8 @@ sto_iter.o\
 sto_wf.o\
 sto_hchi.o\
 sto_che.o\
+sto_forces.o\
+sto_stress_pw.o
 
 OBJS_TOOLS=complexarray.o\
 complexmatrix.o \

source/input.cpp

@@ -1804,6 +1804,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         }
     }
     if(calculation.substr(0,3) != "sto")    bndpar = 1;
+    if(bndpar > GlobalV::NPROC) bndpar = GlobalV::NPROC;
 }
 #ifdef __MPI
 void Input::Bcast()
@@ -2270,9 +2271,9 @@ void Input::Check(void)
             ModuleBase::WARNING_QUIT("Input::Check", "calculate = istate is only availble for LCAO.");
         }
     }
-    else if (calculation == "md") // mohan add 2011-11-04
+    else if (calculation == "md" || calculation == "sto-md") // mohan add 2011-11-04
     {
-        GlobalV::CALCULATION = "md";
+        GlobalV::CALCULATION = calculation;
         symmetry = false;
         cal_force = 1;
         if (mdp.md_nstep == 0)

source/module_esolver/esolver_sdft_pw.cpp

@@ -9,6 +9,11 @@
 #include "../src_pw/global.h"
 #include "../src_pw/symmetry_rho.h"
 //----------------------------------------------
+//-----force-------------------
+#include "../src_pw/sto_forces.h"
+//-----stress------------------
+#include "../src_pw/sto_stress_pw.h"
+//---------------------------------------------------
 
 namespace ModuleESolver
 {
@@ -43,57 +48,7 @@ void ESolver_SDFT_PW::Init(Input &inp, UnitCell_pseudo &cell)
 void ESolver_SDFT_PW::beforescf(const int istep)
 {
     ESolver_KS_PW::beforescf(istep);
-	// if(NITER==0)
-	// {
-	// 	int iter = 1;
-	// 	ETHR = 0.1*DRHO2/ std::max(1.0, CHR.nelec);
-	// 	double *h_diag = new double[wf.npwx * NPOL];
-	// 	for (int ik = 0;ik < kv.nks;ik++)
-	// 	{
-	// 		if(NBANDS > 0 && MY_STOGROUP == 0)
-	// 		{
-	// 			this->c_bands_k(ik,h_diag,istep+1);
-	// 		}
-	// 		else
-	// 		{
-	// 			hm.hpw.init_k(ik); 
-	// 		}
-	// 		if(NBANDS > 0)
-	// 		{
-	// 			MPI_Bcast(wf.evc[ik].c, wf.npwx*NBANDS*2, MPI_DOUBLE , 0, PARAPW_WORLD);
-	// 			MPI_Bcast(wf.ekb[ik], NBANDS, MPI_DOUBLE, 0, PARAPW_WORLD);
-	// 		}
-	// 		stoiter.stoche.ndmin = wf.npw;
-	// 		stoiter.orthog(ik);
-	// 		stoiter.checkemm(ik,iter);
-	// 	}
-	// 	for (int ik = 0;ik < kv.nks;ik++)
-	// 	{
-	// 		//init k
-	// 		if(kv.nks > 1) hm.hpw.init_k(ik);
-	// 		stoiter.stoche.ndmin = wf.npw;
-
-	// 		stoiter.sumpolyval_k(ik);
-	// 	}
-	// 	delete [] h_diag;
-	// 	stoiter.itermu(iter);
-	// 	stoiter.stoche.calcoef(Stochastic_Iter::nroot_fd);
-	// 	for(int ik = 0; ik < kv.nks; ++ik)
-    // 	{
-    // 	    //init k
-    // 	    if(kv.nks > 1) hm.hpw.init_k(ik);
-    // 		stoiter.stoche.ndmin = wf.npw;
-	// 		complex<double> * out, *pchi;
-	// 		out = stowf.shchi[ik].c;
-    //     	if(NBANDS > 0)
-    //     	    pchi = stowf.chiortho[ik].c;
-    //     	else
-    //     	    pchi = stowf.chi0[ik].c;
-	
-    //     	stoiter.stoche.calfinalvec(stoiter.stohchi.hchi_reciprocal, pchi, out, stoiter.nchip[ik]);
-	// 	}
-
-	// }
+	if(istep > 0 && INPUT.nbands_sto != 0) Update_Sto_Orbitals(this->stowf, INPUT.seed_sto);
 }
 
 void ESolver_SDFT_PW::eachiterfinish(int iter)
@@ -292,11 +247,13 @@ void ESolver_SDFT_PW::cal_Energy(energy &en)
 
 void ESolver_SDFT_PW::cal_Force(ModuleBase::matrix &force)
 {
-	
+	Sto_Forces ff;
+    ff.init(force, this->stowf);
 }
 void ESolver_SDFT_PW::cal_Stress(ModuleBase::matrix &stress)
 {
-
+	Sto_Stress_PW ss;
+    ss.cal_stress(stress, this->stowf);
 }
 
 

source/run_pw.cpp

@@ -41,7 +41,7 @@ void Run_pw::plane_wave_line(ModuleESolver::ESolver *p_esolver)
 
 
 
-    if(GlobalV::CALCULATION == "md")
+    if(GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "sto-md")
     {
         Run_MD_PW run_md_pw;
         run_md_pw.md_cells_pw(p_esolver);

source/src_ions/ions.cpp

@@ -141,7 +141,7 @@ void Ions::opt_ions_pw(ModuleESolver::ESolver *p_esolver)
 		time_t fstart = time(NULL);
 
 
-        if (GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax")
+        if (GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax" || GlobalV::CALCULATION.substr(0,3)=="sto")
         {
 			stop = this->after_scf(p_esolver, istep, force_step, stress_step);    // pengfei Li 2018-05-14
 		}
@@ -171,7 +171,7 @@ void Ions::opt_ions_pw(ModuleESolver::ESolver *p_esolver)
 
     }
 
-    if(GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax")
+    if(GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax" || GlobalV::CALCULATION.substr(0,3)=="sto")
     {
         GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
         GlobalV::ofs_running << std::setprecision(16);

source/src_pw/CMakeLists.txt

@@ -26,6 +26,8 @@ list(APPEND objects
     sto_iter.cpp
     sto_hchi.cpp
     sto_wf.cpp
+    sto_forces.cpp 
+    sto_stress_pw.cpp
     stress_func_cc.cpp
     stress_func_ewa.cpp
     stress_func_gga.cpp

source/src_pw/forces.h

@@ -23,7 +23,7 @@ class Forces
 
     void init(ModuleBase::matrix& force);
 
-private:
+protected:
 
     int nat;
 	static double output_acc;

source/src_pw/run_md_pw.cpp

@@ -206,7 +206,7 @@ void Run_MD_PW::md_force_virial(
 
     // mohan added eiter to count for the electron iteration number, 2021-01-28
     int eiter = 0;
-    if (GlobalV::CALCULATION == "md")
+    if (GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "sto-md")
     {
         Electrons elec;
 #ifdef __LCAO