Merge pull request #917 from deepmodeling/SDFT

Merge SDFT to develop

Author: mohanchen

docs/input-main.md

@@ -8,7 +8,7 @@
 
 - [System variables](#system-variables)
 
-    [suffix](#suffix) | [ntype](#ntype) | [calculation](#calculation) | [symmetry](#symmetry) | [kpar](#kpar) | [latname](#latname) | [init_wfc](#init_wfc) | [init_chg](#init_chg) | [init_vel](#init_vel) | [nelec](#nelec) | [tot_magnetization](#tot-magnetization) | [dft_functional](#dft-functional) | [pseudo_type](#pseudo-type) |  [pseudo_rcut](#pseudo-rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem-saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline)
+    [suffix](#suffix) | [ntype](#ntype) | [calculation](#calculation) | [symmetry](#symmetry) | [kpar](#kpar) | [bndpar](#bndpar) | [latname](#latname) | [init_wfc](#init_wfc) | [init_chg](#init_chg) | [init_vel](#init_vel) | [nelec](#nelec) | [tot_magnetization](#tot-magnetization) | [dft_functional](#dft-functional) | [pseudo_type](#pseudo-type) |  [pseudo_rcut](#pseudo-rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem-saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline)
 
 - [Variables related to input files](#variables-related-to-input-files)
 
@@ -24,7 +24,7 @@
 
 - [Electronic structure](#electronic-structure)
 
-    [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing_method](#smearing_method) | [smearing_sigma](#smearing_sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [scf_nmax](#scf_nmax) | [scf_thr](#scf_thr) | [chg_extrap](#chg_extrap)
+    [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_sto](#nbands_sto) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing_method](#smearing_method) | [smearing_sigma](#smearing_sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [scf_nmax](#scf_nmax) | [scf_thr](#scf_thr) | [chg_extrap](#chg_extrap) | [nche_sto](#nche_sto) | [emin_sto](#emin_sto) | [emax_sto](#emax_sto) | [seed_sto](#seed_sto)
 
 - [Geometry relaxation](#geometry-relaxation)
 
@@ -150,6 +150,12 @@ This part of variables are used to control general system parameters.
 - **Description**: devide all processors into kpar groups, and k points will be distributed among each group. The value taken should be less than or equal to the number of k points as well as the number of MPI threads.
 - **Default**: 1
 
+#### bndpar
+
+- **Type**: Integer
+- **Description**: devide all processors into bndpar groups, and bands (only stochastic orbitals now) will be distributed among each group. It should be larger than 0.
+- **Default**: 1
+
 #### latname
 
 - **Type**: String
@@ -453,6 +459,14 @@ calculations.
   - nspin=1: 1.2\*occupied_bands, occupied_bands+10)
   - nspin=2: max(1.2\*nelec, nelec+20)
 
+#### nbands_sto
+
+- **Type**: Integer
+- **Description**: 
+  - nbands_sto>0: Number of stochastic orbitals to calculate in stochastic DFT (SDFT) or mix stochastic-deterministic (MDFT). More bands obtain more precise results or smaller stochastic errors ($ \propto 1/\sqrt{N_{\chi}}$); 
+  - nbands_sto=0: Complete basis will be used to replace stochastic orbitals with the Chebyshev method (CT) and it will get the results the same as KSDFT without stochastic errors.
+- **Default**: 0
+
 #### nbands_istate
 
 - **Type**: Integer
@@ -551,6 +565,32 @@ calculations.
   - second-order: second-order extrapolation
 - **Default**:atomic
 
+#### nche_sto
+
+- **Type**: Integer
+- **Description**: Chebyshev expansion orders for SDFT, MDFT, CT methods.
+- **Default**:5
+
+#### emin_sto
+
+- **Type**: Real
+- **Description**: Trial energy to guess the lower bound of eigen energies of the Hamitonian Operator $\hat{H}$. The unit is Ry.
+- **Default**:0.0
+
+#### emax_sto
+
+- **Type**: Real
+- **Description**: Trial energy to guess the upper bound of eigen energies of the Hamitonian Operator $\hat{H}$. The unit is Ry.
+- **Default**:0.0
+
+#### seed_sto
+
+- **Type**: Integer
+- **Description**: The random seed to generate stochastic orbitals.
+  - seed_sto>=0: Stochastic orbitals have the form of $\exp(i2\pi\theta(G))$, where $\theta$ is a uniform distribution in $(0,1)$. If seed_sto=0, the seed is decided by time(NULL).
+  - seed_sto<=-1: Stochastic orbitals have the form of $\pm1$ with the equal probability. If seed_sto=-1, the seed is decided by time(NULL).
+- **Default**:0
+
 ### Geometry relaxation
 
 This part of variables are used to control the geometry relaxation.

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\
 

source/driver.cpp

@@ -96,7 +96,8 @@ void Driver::atomic_world(void)
 	ModuleESolver::ESolver *p_esolver;
 	if(GlobalV::BASIS_TYPE=="pw" || GlobalV::BASIS_TYPE=="lcao_in_pw")
 	{
-		use_ensol = "ksdft_pw"; 
+		if(GlobalV::CALCULATION.substr(0,3) == "sto")	use_ensol = "sdft_pw";
+		else											use_ensol = "ksdft_pw"; 
 		//We set it temporarily
 		//Finally, we have ksdft_pw, ksdft_lcao, sdft_pw, ofdft, lj, eam, etc.
 		ModuleESolver::init_esolver(p_esolver, use_ensol);

source/input.cpp

@@ -138,9 +138,9 @@ void Input::Default(void)
     nbands_sto = 0;
     nbands_istate = 5;
     pw_seed = 1;
-    nche_sto = 0;
+    nche_sto = 5;
     seed_sto = 0;
-    stotype = "pw";
+    bndpar = 1;
     kpar = 1;
     berry_phase = false;
     gdir = 3;
@@ -519,8 +519,8 @@ bool Input::Read(const std::string &fn)
         else if (strcmp("nbands", word) == 0) // number of atom bands
         {
             read_value(ifs, nbands);
-            if (nbands <= 0)
-                ModuleBase::WARNING_QUIT("Input", "NBANDS must > 0");
+            if (nbands < 0)
+                ModuleBase::WARNING_QUIT("Input", "NBANDS must >= 0");
         }
         else if (strcmp("nbands_sto", word) == 0) // number of stochastic bands
         {
@@ -553,9 +553,9 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, emin_sto);
         }
-        else if (strcmp("stotype", word) == 0)
+        else if (strcmp("bndpar", word) == 0)
         {
-            read_value(ifs, stotype);
+            read_value(ifs, bndpar);
         }
         else if (strcmp("kpar", word) == 0) // number of pools
         {
@@ -1803,6 +1803,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
             vdw_radius = "95";
         }
     }
+    if(calculation.substr(0,3) != "sto")    bndpar = 1;
 }
 #ifdef __MPI
 void Input::Bcast()
@@ -1833,7 +1834,7 @@ void Input::Bcast()
     Parallel_Common::bcast_int(pw_seed);
     Parallel_Common::bcast_double(emax_sto);
     Parallel_Common::bcast_double(emin_sto);
-    Parallel_Common::bcast_string(stotype);
+    Parallel_Common::bcast_int(bndpar);
     Parallel_Common::bcast_int(kpar);
     Parallel_Common::bcast_bool(berry_phase);
     Parallel_Common::bcast_int(gdir);
@@ -2124,7 +2125,7 @@ void Input::Check(void)
     ModuleBase::TITLE("Input", "Check");
 
     if (nbands < 0)
-        ModuleBase::WARNING_QUIT("Input", "NBANDS must > 0");
+        ModuleBase::WARNING_QUIT("Input", "NBANDS must >= 0");
     //	if(nbands_istate < 0) ModuleBase::WARNING_QUIT("Input","NBANDS_ISTATE must > 0");
     if (nb2d < 0)
         ModuleBase::WARNING_QUIT("Input", "nb2d must > 0");
@@ -2180,7 +2181,7 @@ void Input::Check(void)
         */
         this->relax_nmax = 1;
     }
-    else if (calculation == "scf-sto") // qianrui 2021-2-20
+    else if (calculation == "sto-scf") // qianrui 2021-2-20
     {
         if (mem_saver == 1)
         {

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
+    int bndpar; //parallel for stochastic/deterministic bands
 
     //==========================================================
     // electrons / spin
@@ -179,9 +179,9 @@ class Input
     //==========================================================
     // charge mixing
     //==========================================================
-    std::string mixing_mode; // "plain","broden",...
+    std::string mixing_mode; // "plain","broyden",...
     double mixing_beta; // 0 : no_mixing
-    int mixing_ndim; // used in Broden method
+    int mixing_ndim; // used in Broyden method
     double mixing_gg0; // used in kerker method. mohan add 2014-09-27
 
     //==========================================================

source/input_conv.cpp

@@ -61,6 +61,7 @@ void Input_Conv::Convert(void)
     GlobalV::KPAR = temp_nproc;
 #else
     GlobalV::KPAR = INPUT.kpar;
+    GlobalV::NSTOGROUP = INPUT.bndpar;
 #endif
     GlobalV::CALCULATION = INPUT.calculation;
 

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_elecstate/elecstate_pw.cpp

@@ -41,15 +41,19 @@ void ElecStatePW::parallelK()
 {
 #ifdef __MPI
     charge->rho_mpi();
-    if (GlobalV::CALCULATION != "scf-sto" && GlobalV::CALCULATION != "relax-sto"
-        && GlobalV::CALCULATION != "md-sto") // qinarui add it temporarily.
-    {
-        //==================================
-        // Reduce all the Energy in each cpu
-        //==================================
-        this->eband /= GlobalV::NPROC_IN_POOL;
-        Parallel_Reduce::reduce_double_all(this->eband);
-    }
+    if(GlobalV::CALCULATION.substr(0,3) == "sto") //qinarui add it 2021-7-21
+	{
+		GlobalC::en.eband /= GlobalV::NPROC_IN_POOL;
+		MPI_Allreduce(MPI_IN_PLACE, &GlobalC::en.eband, 1, MPI_DOUBLE, MPI_SUM , STO_WORLD);
+	}
+	else
+	{
+    	//==================================
+    	// Reduce all the Energy in each cpu
+    	//==================================
+		GlobalC::en.eband /= GlobalV::NPROC_IN_POOL;
+		Parallel_Reduce::reduce_double_all( GlobalC::en.eband );
+	}
 #endif
     return;
 }