Merge branch 'develop' of github.com:deepmodeling/abacus-develop into HSolver

Author: dyzheng

README.md

@@ -17,12 +17,14 @@ ABACUS is an electronic structure package based on density functional theory(DFT
 Please refer to our [GitHub repository](https://github.com/deepmodeling/abacus-develop) for more information and support.
 # Table of contents
 
+- [Table of contents](#table-of-contents)
 - [Features](#features)
 - [Download and install](#download-and-install)
 - [Quickstart guide](#quickstart-guide)
   - [Input files](#input-files)
+  - [Run ABACUS](#run-abacus)
   - [Output files](#output-files)
-- [Features](#features)
+- [Features](#features-1)
 - [Functionalities](#functionalities)
 - [Examples](#examples)
 - [For developers](#for-developers)
@@ -92,6 +94,22 @@ The following files are the central input files for ABACUS. Before executing the
 
 [back to top](#readme-top)
 
+## Run ABACUS
+
+After putting all required input files under one folder, enter this folder.
+```bash
+cd input_folder
+```
+
+Perform calculation by:
+```bash
+mpirun -np 4 abacus
+```
+
+You can replace `4` with your desired number of process, typically the result of the command `nproc`.
+
+[back to top](#readme-top)
+
 ## Output files
 
 When the calculation finishes, the program will create an output directory (default: OUT.ABACUS/),

docs/input-main.md

@@ -28,7 +28,7 @@
 
 - [Electronic structure (SDFT)](#electronic-structure-sdft)
 
-    [nbands_sto](#nbands_sto) | [nche_sto](#nche_sto) | [emin_sto](#emin_sto) | [emax_sto](#emax_sto) | [seed_sto](#seed_sto)
+    [method_sto](#method_sto) | [nbands_sto](#nbands_sto) | [nche_sto](#nche_sto) | [emin_sto](#emin_sto) | [emax_sto](#emax_sto) | [seed_sto](#seed_sto) | [initsto_freq](#initsto_freq)
 
 - [Geometry relaxation](#geometry-relaxation)
 
@@ -469,7 +469,7 @@ calculations.
 #### nbands
 
 - **Type**: Integer
-- **Description**: Number of bands to calculate. It is recommended you setup this value, especially when you use smearing techniques, more bands should be included.
+- **Description**: Number of Kohn-Sham orbitals to calculate. It is recommended you setup this value, especially when you use smearing techniques, more bands should be included.
 - **Default**:
   - nspin=1: 1.2\*occupied_bands, occupied_bands+10)
   - nspin=2: max(1.2\*nelec, nelec+20)
@@ -577,14 +577,25 @@ calculations.
 
 ### Electronic structure (SDFT)
 
-This part of variables are used to control the parameters of stochastic DFT (SDFT),  mix stochastic-deterministic DFT (MDFT), or complete-basis Chebyshev method (CT).
+This part of variables are used to control the parameters of stochastic DFT (SDFT),  mix stochastic-deterministic DFT (MDFT), or complete-basis Chebyshev method (CT). We now support calculation of "sto-scf" and "sto-md".
+
+#### method_sto
+
+- **Type**: Integer
+- **Description**: 
+  - Different method to do SDFT.
+  - 1: SDFT calculates $T_n(\hat{h})\ket{\chi}$ twice, where $T_n(x)$ is the n-th order Chebyshev polynomial and $\hat{h}=\frac{\hat{H}-\bar{E}}{\Delta E}$ owning eigen-value $\in(-1,1)$. This method cost less memory but slow.
+  - 2: SDFT calculates $T_n(\hat{h})\ket{\chi}$ once but need much more memory. This method is fast but when memory is not enough. Only method 1 can be used.
+  - other: use 1
+- **Default**: 1
 
 #### nbands_sto
 
 - **Type**: Integer
 - **Description**: 
   - nbands_sto>0: Number of stochastic orbitals to calculate in SDFT and 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.
+  - If you want to do MDFT. [nbands](#nbands) which represents the number of KS orbitals should be set.
 - **Default**: 0
 
 #### nche_sto
@@ -609,8 +620,15 @@ This part of variables are used to control the parameters of stochastic DFT (SDF
 
 - **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
+  - 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
+
+#### initsto_freq
+
+- **Type**: Integer
+- **Description**: Frequency (once each initsto_freq steps) to generate new stochastic orbitals when running md.
+- **Default**:1000
 
 ### Geometry relaxation
 

source/Makefile.Objects

@@ -219,7 +219,6 @@ parallel_kpoints.o\
 parallel_common.o\
 parallel_reduce.o\
 parallel_pw.o\
-ft.o\
 parallel_grid.o\
 
 OBJS_ESOLVER=esolver.o\
@@ -273,11 +272,13 @@ diago_elpa.o\
 diago_iter_assist.o\
 hsolver_lcao.o\
 hsolver_pw.o\
+hsolver_pw_sdft.o
 
 OBJ_ELECSTATES=elecstate.o\
 dm2d_to_grid.o\
 elecstate_lcao.o\
 elecstate_pw.o\
+elecstate_pw_sdft.o
 
 OBJ_PSI=psi.o
 

source/input.cpp

@@ -143,6 +143,8 @@ void Input::Default(void)
     seed_sto = 0;
     bndpar = 1;
     kpar = 1;
+    initsto_freq = 1000;
+    method_sto = 1;
     berry_phase = false;
     gdir = 3;
     towannier90 = false;
@@ -421,6 +423,10 @@ void Input::Default(void)
     sigma_k = 0.6;
     nc_k = 0.00037;
 
+    //==========================================================
+    //    compensating charge        donghs added on 2022-06-23
+    //==========================================================
+    comp_chg = 0;
     comp_q = 0.0;
     comp_l = 1.0;
     comp_center = 0.0;
@@ -569,6 +575,14 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, emin_sto);
         }
+        else if (strcmp("initsto_freq", word) == 0)
+        {
+            read_value(ifs, initsto_freq);
+        }
+        else if (strcmp("method_sto", word) == 0)
+        {
+            read_value(ifs, method_sto);
+        }
         else if (strcmp("bndpar", word) == 0)
         {
             read_value(ifs, bndpar);
@@ -1503,6 +1517,10 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, nc_k);
         }
+        else if (strcmp("comp_chg", word) == 0)
+        {
+            read_value(ifs, comp_chg);
+        }
         else if (strcmp("comp_q", word) == 0)
         {
             read_value(ifs, comp_q);
@@ -1899,6 +1917,8 @@ void Input::Bcast()
     Parallel_Common::bcast_int(pw_seed);
     Parallel_Common::bcast_double(emax_sto);
     Parallel_Common::bcast_double(emin_sto);
+    Parallel_Common::bcast_int(initsto_freq);
+    Parallel_Common::bcast_int(method_sto);
     Parallel_Common::bcast_int(bndpar);
     Parallel_Common::bcast_int(kpar);
     Parallel_Common::bcast_bool(berry_phase);
@@ -2189,6 +2209,12 @@ void Input::Bcast()
     Parallel_Common::bcast_double(sigma_k);
     Parallel_Common::bcast_double(nc_k);
 
+    Parallel_Common::bcast_bool(comp_chg);
+    Parallel_Common::bcast_double(comp_q);
+    Parallel_Common::bcast_double(comp_l);
+    Parallel_Common::bcast_double(comp_center);
+    Parallel_Common::bcast_int(comp_dim);
+
     return;
 }
 #endif

source/input.h

@@ -70,6 +70,8 @@ class Input
     double emax_sto; // Emax & Emin to normalize H
     double emin_sto;
     int bndpar; //parallel for stochastic/deterministic bands
+    int initsto_freq; //frequency to init stochastic orbitals when running md
+    int method_sto; //different methods for sdft, 1: slow, less memory  2: fast, more memory
 
     //==========================================================
     // electrons / spin
@@ -400,7 +402,8 @@ class Input
     double tau;
     double sigma_k;
     double nc_k;
-    // gauss charge
+    // compensating charge
+    bool comp_chg;
     double comp_q;
     double comp_l;
     double comp_center;

source/input_conv.cpp

@@ -494,6 +494,10 @@ void Input_Conv::Convert(void)
     GlobalV::sigma_k = INPUT.sigma_k;
     GlobalV::nc_k = INPUT.nc_k;
 
+    //-----------------------------------------------
+    // compensating charge
+    //-----------------------------------------------
+    GlobalV::comp_chg = INPUT.comp_chg;
     GlobalC::solvent_model.comp_q = INPUT.comp_q;
     GlobalC::solvent_model.comp_l = INPUT.comp_l;
     GlobalC::solvent_model.comp_center = INPUT.comp_center;

source/module_base/global_variable.cpp

@@ -205,4 +205,6 @@ double eb_k = 80.0;
 double tau = 1.0798 * 1e-5;
 double sigma_k = 0.6;
 double nc_k = 0.00037;
+
+bool comp_chg = false; // compensating charge
 } // namespace GlobalV

source/module_base/global_variable.h

@@ -232,5 +232,8 @@ extern double eb_k;
 extern double tau;
 extern double sigma_k;
 extern double nc_k;
+
+// compensating charge
+extern bool comp_chg;
 } // namespace GlobalV
 #endif

source/module_base/math_chebyshev.h

@@ -104,7 +104,7 @@ class Chebyshev
     void calfinalvec(T *ptr, 
 		void (T::*funA)(REAL *in, REAL *out, const int), 
 		REAL *wavein, REAL *waveout, 
-		const int N, const int LDA = 1,  const int m = 1);
+		const int N, const int LDA = 1,  const int m = 1); //do not define yet
 
 	// calfinalvec_real uses C_n[f], where f is a function of real number and A is a complex Operator.
 	template<class T>
@@ -129,8 +129,25 @@ class Chebyshev
 		T *ptr, void (T::*funA)(std::complex<REAL> *in, std::complex<REAL> *out, const int), 
 		std::complex<REAL> *wavein, 
 		const int N, const int LDA = 1,  const int m = 1);
-	//get T_n(x)
+	
+	// get T_n(x)
 	void getpolyval(REAL x, REAL* polyval, const int N);
+
+	// get each order of vector: {T_0(A)v, T_1(A)v, ..., T_n(A)v}
+	// Note: use it carefully, it will cost a lot of memory!
+	// calpolyvec_real: f(x) = \sum_n C_n*T_n(x), f is a real function
+	template<class T>
+	void calpolyvec_real(T *ptr, 
+		void (T::*funA)(REAL *in, REAL *out, const int), 
+		REAL *wavein, REAL *waveout, 
+		const int N, const int LDA = 1,  const int m = 1);//do not define yet
+	// calpolyvec_complex: f(x) = \sum_n C_n*T_n(x), f is a complex function
+	template<class T>
+	void calpolyvec_complex(T *ptr, 
+		void (T::*funA)(std::complex<REAL> *in, std::complex<REAL> *out, const int), 
+		std::complex<REAL> *wavein, std::complex<REAL> *waveout, 
+		const int N, const int LDA = 1,  const int m = 1);
+
 
 	// IV.
 	// recurs fomula: v_{n+1} = 2Av_n - v_{n-1}

source/module_base/math_chebyshev_def.h

@@ -339,6 +339,40 @@ void Chebyshev<REAL>::calfinalvec_complex(T *ptr,
     return;
 }
 
+template<typename REAL>
+template<class T>
+void Chebyshev<REAL>::calpolyvec_complex(T *ptr, 
+	void (T::*funA)(std::complex<REAL> *in, std::complex<REAL> *out, const int), 
+	std::complex<REAL> *wavein, 
+	std::complex<REAL> *polywaveout, 
+	const int N, const int LDA, const int m)
+{
+
+    assert(N>=0 && LDA >= N);
+    int ndmxt;
+    if(m == 1) ndmxt = N * m;
+    else       ndmxt = LDA * m; 
+
+    std::complex<REAL> *arraynp1 = polywaveout + 2 * ndmxt;
+	std::complex<REAL> *arrayn = polywaveout + ndmxt;
+	std::complex<REAL> *arrayn_1 = polywaveout;
+  
+    ModuleBase::GlobalFunc::DCOPY(wavein, arrayn_1, ndmxt);
+    
+    //1-st order
+    (ptr->*funA)(arrayn_1, arrayn, m);
+
+    //more than 1-st orders
+    for(int ior = 2; ior < norder; ++ior)
+    {
+        recurs_complex(ptr, funA, arraynp1, arrayn, arrayn_1, N, LDA, m);
+        arrayn_1 += ndmxt;
+        arrayn += ndmxt;
+        arraynp1 += ndmxt; 
+    }
+    return;
+}
+
 template<typename REAL>
 template<class T>
 void Chebyshev<REAL>::tracepolyA(