refactor: add SDFT for ESolver_SDFT_PW

	move SDFT program in sdf-version ABACUS to ESolver_SDFT_PW
	Save on schedule and no tests for SDFT are set up.

Author: Qianruipku

CMakeLists.txt

@@ -17,7 +17,7 @@ option(USE_CUDA "Enable support to CUDA." OFF)
 option(USE_ROCM "Enable support to ROCm." OFF)
 option(USE_OPENMP " Enable OpenMP in abacus." ON)
 option(ENABLE_ASAN "Enable AddressSanitizer" OFF)
-option(BUILD_TESTING "Build ABACUS unit tests" OFF)
+option(BUILD_TESTING "Build ABACUS unit tests" ON)
 option(GENERATE_TEST_REPORTS "Enable test report generation" OFF)
 
 set(ABACUS_BIN_NAME abacus)

source/Makefile.Objects

@@ -57,9 +57,8 @@ soc.o\
 to_wannier90.o \
 unk_overlap_pw.o \
 berryphase.o \
-sto_elec.o\
-sto_wf.o\
 sto_iter.o\
+sto_wf.o\
 sto_hchi.o\
 sto_che.o\
 
@@ -226,6 +225,7 @@ parallel_reduce.o\
 parallel_pw.o\
 ft.o\
 parallel_grid.o\
+parallel_stochi.o\
 
 OBJS_ESOLVER=esolver.o\
     esolver_ks.o\

source/Makefile.vars

@@ -4,12 +4,12 @@ CPLUSPLUS_MPI = mpiicpc
 
 LAPACK_DIR    = $(MKLROOT)
 
-FFTW_DIR      = /public/software/fftw_3.3.8
+FFTW_DIR      = /home/qianrui/intelcompile/fftw_3.3.8
 
-ELPA_DIR      = /public/software/elpa_21.05.002
+ELPA_DIR      = /home/qianrui/intelcompile/elpa_21.05.002
 ELPA_INCLUDE_DIR = ${ELPA_DIR}/include/elpa-2021.05.002
 
-CEREAL_DIR    = /public/software/cereal
+CEREAL_DIR    = /home/qianrui/headfile/cereal
 
 # LIBXC_DIR     = /public/software/libxc-5.0.0
 

source/input.h

@@ -64,11 +64,11 @@ class Input
     // Stochastic DFT
     //==========================================================
     int nche_sto; // number of orders for Chebyshev expansion in stochastic DFT //qinarui 2021-2-5
+    int nbands_sto;			// number of stochastic bands //qianrui 2021-2-5
     int seed_sto; // random seed for sDFT
     double emax_sto; // Emax & Emin to normalize H
     double emin_sto;
     std::string stotype;
-    int nbands_sto; // number of stochastic bands //qianrui 2021-2-5
 
     //==========================================================
     // electrons / spin

source/module_base/global_variable.cpp

@@ -101,10 +101,14 @@ int out_mul = 0; // qifeng add 2019/9/10
 //----------------------------------------------------------
 int NPROC = 1;
 int KPAR = 1;
+int NSTOGROUP = 1;
 int MY_RANK = 0;
 int MY_POOL = 0;
+int MY_STOGROUP = 0;
 int NPROC_IN_POOL = 1;
+int NPROC_IN_STOGROUP = 1;
 int RANK_IN_POOL = 0;
+int RANK_IN_STOGROUP = 0;
 int DRANK = -1; // mohan add 2012-01-13, must be -1, so we can recognize who didn't in DIAG_WORLD
 int DSIZE = KPAR;
 int DCOLOR = -1;

source/module_base/global_variable.h

@@ -116,10 +116,14 @@ extern int out_mul; // qifeng add 2019/9/10
 //========================================================================
 extern int NPROC;
 extern int KPAR;
+extern int NSTOGROUP;
 extern int MY_RANK;
 extern int MY_POOL;
+extern int MY_STOGROUP;
 extern int NPROC_IN_POOL;
+extern int NPROC_IN_STOGROUP;
 extern int RANK_IN_POOL;
+extern int RANK_IN_STOGROUP;
 extern int DRANK;
 extern int DSIZE;
 extern int DCOLOR;

source/module_cell/unitcell_pseudo.h

@@ -31,7 +31,7 @@ class UnitCell_pseudo : public UnitCell
 	int lmax_ppwf;
 	int lmaxmax; // liuyu 2021-07-04
 	bool init_vel; // liuyu 2021-07-15
-	//double nelec;
+	double nelec;
 
 public: // member functions
 	UnitCell_pseudo();

source/module_esolver/esolver_ks_pw.cpp

@@ -211,7 +211,10 @@ void ESolver_KS_PW:: eachiterinit(const int iter)
 
     //(2) save change density as previous charge,
     // prepared fox mixing.
-    GlobalC::CHR.save_rho_before_sum_band();
+    if(GlobalV::MY_STOGROUP == 0)
+	{
+        GlobalC::CHR.save_rho_before_sum_band();
+    }
 }
 
 //Temporary, it should be replaced by hsolver later.

source/module_esolver/esolver_sdft_pw.cpp

@@ -1,4 +1,14 @@
 #include "./esolver_sdft_pw.h"
+#include "time.h"
+#include <fstream>
+#include <algorithm>
+#include "../module_base/timer.h"
+
+//-------------------Temporary------------------
+#include "../module_base/global_variable.h"
+#include "../src_pw/global.h"
+#include "../src_pw/symmetry_rho.h"
+//----------------------------------------------
 
 namespace ModuleESolver
 {
@@ -9,18 +19,22 @@ ESolver_SDFT_PW::ESolver_SDFT_PW()
     basisname = "PW";
 }
 
+ESolver_SDFT_PW::~ESolver_SDFT_PW()
+{
+}
+
 void ESolver_SDFT_PW::Init(Input &inp, UnitCell_pseudo &cell)
 {
-    // ESolver_KS_PW::Init(inp, cell);
-    // STO_WF.alloc();
-	// stoiter.alloc( wf.npwx );
+    ESolver_KS_PW::Init(inp, cell);
+	stowf.init(GlobalC::kv.nks);
+	if(INPUT.nbands_sto != 0)	Init_Sto_Orbitals(this->stowf, INPUT.seed_sto);
+	else						Init_Com_Orbitals(this->stowf, GlobalC::kv);
+	stoiter.init(GlobalC::wf.npwx, this->stowf.nchip);
 }
 
 void ESolver_SDFT_PW::beforescf()
 {
     ESolver_KS_PW::beforescf();
-    // STO_WF.init();
-	// stoiter.init( wf.npwx );
 	// if(NITER==0)
 	// {
 	// 	int iter = 1;
@@ -62,11 +76,11 @@ void ESolver_SDFT_PW::beforescf()
     // 	    if(kv.nks > 1) hm.hpw.init_k(ik);
     // 		stoiter.stoche.ndmin = wf.npw;
 	// 		complex<double> * out, *pchi;
-	// 		out = STO_WF.shchi[ik].c;
+	// 		out = stowf.shchi[ik].c;
     //     	if(NBANDS > 0)
-    //     	    pchi = STO_WF.chiortho[ik].c;
+    //     	    pchi = stowf.chiortho[ik].c;
     //     	else
-    //     	    pchi = STO_WF.chi0[ik].c;
+    //     	    pchi = stowf.chi0[ik].c;
 
     //     	stoiter.stoche.calfinalvec(stoiter.stohchi.hchi_reciprocal, pchi, out, stoiter.nchip[ik]);
 	// 	}
@@ -76,21 +90,190 @@ void ESolver_SDFT_PW::beforescf()
 
 void ESolver_SDFT_PW::eachiterfinish(int iter, bool conv_elec)
 {
-
+	//print_eigenvalue(GlobalV::ofs_running);
+    GlobalC::en.calculate_etot();
 }
 void ESolver_SDFT_PW::afterscf(bool conv_elec)
 {
-
+	for(int is=0; is<GlobalV::NSPIN; is++)
+    {
+        std::stringstream ssc;
+        std::stringstream ss1;
+        ssc << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG";
+		ss1 << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG.cube";
+        GlobalC::CHR.write_rho(GlobalC::CHR.rho_save[is], is, 0, ssc.str() );//mohan add 2007-10-17
+	    GlobalC::CHR.write_rho_cube(GlobalC::CHR.rho_save[is], is, ss1.str(), 3);
+    }
+    if(conv_elec)
+    {
+        //GlobalV::ofs_running << " convergence is achieved" << std::endl;			
+        //GlobalV::ofs_running << " !FINAL_ETOT_IS " << GlobalC::en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; 
+        GlobalV::ofs_running << "\n charge density convergence is achieved" << std::endl;
+        GlobalV::ofs_running << " final etot is " << GlobalC::en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
+    }
+    else
+    {
+        GlobalV::ofs_running << " convergence has NOT been achieved!" << std::endl;
+    }
 }
 
 void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr)
 {
-    
+    double *h_diag = new double[GlobalC::wf.npwx * GlobalV::NPOL];
+	GlobalV::ofs_running << " "  <<setw(8) << "K-point" << setw(15) << "CG iter num" << setw(15) << "Time(Sec)"<< std::endl;
+	GlobalV::ofs_running << setprecision(6) << setiosflags(ios::fixed) << setiosflags(ios::showpoint);
+	for (int ik = 0;ik < GlobalC::kv.nks;ik++)
+	{
+		if(GlobalV::NBANDS > 0 && GlobalV::MY_STOGROUP == 0)
+		{
+			this->c_bands_k(ik,h_diag,istep+1,iter);
+		}
+		else
+		{
+			GlobalC::hm.hpw.init_k(ik); 
+			//In fact, hm.hpw.init_k has been done in wf.wfcinit();
+		}
+		
+#ifdef __MPI
+			if(GlobalV::NBANDS > 0)
+			{
+				MPI_Bcast(GlobalC::wf.evc[ik].c, GlobalC::wf.npwx*GlobalV::NBANDS*2, MPI_DOUBLE , 0, PARAPW_WORLD);
+				MPI_Bcast(GlobalC::wf.ekb[ik], GlobalV::NBANDS, MPI_DOUBLE, 0, PARAPW_WORLD);
+			}
+#endif
+			stoiter.stoche.ndmin = GlobalC::wf.npw;
+			stoiter.orthog(ik,this->stowf);
+			stoiter.checkemm(ik,iter,this->stowf);	//check and reset emax & emin
+		}
+		for (int ik = 0;ik < GlobalC::kv.nks;ik++)
+		{
+			//init k
+			if(GlobalC::kv.nks > 1) GlobalC::hm.hpw.init_k(ik);
+			stoiter.stoche.ndmin = GlobalC::wf.npw;
+
+			stoiter.sumpolyval_k(ik, this->stowf);
+		}
+		delete [] h_diag;
+		GlobalC::en.eband  = 0.0;
+        GlobalC::en.demet  = 0.0;
+        GlobalC::en.ef     = 0.0;
+        GlobalC::en.ef_up  = 0.0;
+        GlobalC::en.ef_dw  = 0.0;
+		stoiter.itermu(iter);
+		//(5) calculate new charge density 
+		// calculate KS rho.
+		if(GlobalV::NBANDS > 0)
+		{
+			if(GlobalV::MY_STOGROUP == 0)
+			{
+				GlobalC::CHR.sum_band();
+			}
+			else
+			{
+				for(int is=0; is < GlobalV::NSPIN; is++)
+				{
+					ModuleBase::GlobalFunc::ZEROS(GlobalC::CHR.rho[is], GlobalC::pw.nrxx);
+				}
+			}
+			MPI_Bcast(&GlobalC::en.eband,1, MPI_DOUBLE, 0,PARAPW_WORLD);
+		}
+		else
+		{
+			for(int is=0; is < GlobalV::NSPIN; is++)
+			{
+				ModuleBase::GlobalFunc::ZEROS(GlobalC::CHR.rho[is], GlobalC::pw.nrxx);
+			}
+		}
+	// calculate stochastic rho
+		stoiter.sum_stoband(this->stowf);
+		
+
+		//(6) calculate the delta_harris energy 
+		// according to new charge density.
+		// mohan add 2009-01-23
+		//en.calculate_harris(2);
+
+		if(GlobalV::MY_STOGROUP==0)
+		{
+			Symmetry_rho srho;
+			for(int is=0; is < GlobalV::NSPIN; is++)
+			{
+				srho.begin(is, GlobalC::CHR,GlobalC::pw, GlobalC::Pgrid, GlobalC::symm);
+			}
+		}
+		else
+		{
+			if(ModuleSymmetry::Symmetry::symm_flag)	MPI_Barrier(MPI_COMM_WORLD);
+		}
 }
 
-void ESolver_SDFT_PW::cal_Energy(energy &en)
+void ESolver_SDFT_PW:: c_bands_k(const int ik, double* h_diag, const int istep, const int iter)
 {
+	ModuleBase::timer::tick(this->classname,"c_bands_k");
+	int precondition_type = 2;
+	GlobalC::hm.hpw.init_k(ik);
+	
+    //===========================================
+    // Conjugate-Gradient diagonalization
+    // h_diag is the precondition matrix
+    // h_diag(1:npw) = MAX( 1.0, g2kin(1:npw) );
+    //===========================================
+    if (precondition_type==1)
+    {
+        for (int ig = 0;ig < GlobalC::wf.npw; ++ig)
+		{
+			h_diag[ig] = std::max(1.0, GlobalC::wf.g2kin[ig]);
+			if(GlobalV::NPOL==2) h_diag[ig+GlobalC::wf.npwx] = h_diag[ig];
+		}
+    }
+    else if (precondition_type==2)
+    {
+        for (int ig = 0;ig < GlobalC::wf.npw; ig++)
+		{
+			h_diag[ig] = 1 + GlobalC::wf.g2kin[ig] + sqrt( 1 + (GlobalC::wf.g2kin[ig] - 1) * (GlobalC::wf.g2kin[ig] - 1));
+			if(GlobalV::NPOL==2) h_diag[ig+GlobalC::wf.npwx] = h_diag[ig];
+		}
+    }
+	//h_diag can't be zero!  //zhengdy-soc
+	if(GlobalV::NPOL==2)
+	{
+		for(int ig = GlobalC::wf.npw;ig < GlobalC::wf.npwx; ig++)
+		{
+			h_diag[ig] = 1.0;
+			h_diag[ig+ GlobalC::wf.npwx] = 1.0;
+		}
+	}
+	clock_t start=clock();
+
+	//============================================================
+	// diago the hamiltonian!!
+	// In plane wave method, firstly using cinitcgg to diagnolize,
+	// then using cg method.
+	//
+	// In localized orbital presented in plane wave case,
+	// only using cinitcgg.
+	//
+	// In linear scaling method, using sparse matrix and
+	// adjacent searching code and cg method to calculate the
+	// eigenstates.
+	//=============================================================
+	double avg_iter_k = 0.0;
+	GlobalC::hm.diagH_pw(istep, iter, ik, h_diag, avg_iter_k);
+
+	GlobalC::en.print_band(ik);
+	clock_t finish=clock();
+	const double duration = static_cast<double>(finish - start) / CLOCKS_PER_SEC;
+	GlobalV::ofs_running << " " << setw(8) 
+		<< ik+1 << setw(15) 
+		<< avg_iter_k << setw(15) << duration << endl;
+
+	ModuleBase::timer::tick(this->classname,"c_bands_k");
+}
 
+
+void ESolver_SDFT_PW::cal_Energy(energy &en)
+{
+	
 }
 
 void ESolver_SDFT_PW::cal_Force(ModuleBase::matrix &force)

source/module_esolver/esolver_sdft_pw.h

@@ -1,22 +1,36 @@
 #include "./esolver_ks_pw.h"
+#include "../src_pw/sto_wf.h"
+#include "../src_pw/sto_iter.h"
+#include "../src_pw/sto_che.h"
+#include "../src_pw/sto_hchi.h"
+
 namespace ModuleESolver
 {
 
 class ESolver_SDFT_PW: public ESolver_KS_PW
 {
 public:
     ESolver_SDFT_PW();
+    ~ESolver_SDFT_PW();
     void Init(Input &inp, UnitCell_pseudo &cell) override;
     void cal_Energy(energy& en) override;
-    void cal_Force(ModuleBase::matrix &force) override;
-    void cal_Stress(ModuleBase::matrix &stress) override;
+    void cal_Force(ModuleBase::matrix& force) override;
+    void cal_Stress(ModuleBase::matrix& stress) override;
+public:
+    Stochastic_WF stowf;
+    Stochastic_Iter stoiter;
+    // Stochastic_Chebychev stoche;
+    // Stochastic_hchi stohchi;
 
 protected:
     virtual void beforescf() override; 
     // virtual void eachiterinit(int iter) override; 
     virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
     virtual void eachiterfinish(const int iter, const bool conv) override; 
     virtual void afterscf(const bool) override;
+private:
+    void c_bands_k(const int ik, double* h_diag, const int istep, const int iter);
+
 };
 
 }