refactor: finish forces and stress for sdft

md is avilable for sdft with "calculation=sto-md"

Author: Qianruipku

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

@@ -2262,9 +2262,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,6 +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(istep > 0 && INPUT.nbands_sto != 0) Update_Sto_Orbitals(this->stowf, INPUT.seed_sto);
 	// if(NITER==0)
 	// {
 	// 	int iter = 1;
@@ -292,11 +298,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