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

Author: dyzheng

.gitignore

@@ -12,3 +12,7 @@ result.out
 html
 *.log
 STRU_READIN_ADJUST.cif
+*.egg
+*.egg-info
+build
+dist

docs/input-main.md

@@ -1051,13 +1051,26 @@ This part of variables are relevant when using hybrid functionals
 
 #### exx_hse_omega
 
-- **Type**:
+- **Type**: Real
 - **Description**: range-separation parameter in HSE functional, such that $1/r=erfc(\omega r)/r+erf(\omega r)/r$.
 - **Default**: 0.11
 
-adial integration for pseudopotentials, in Bohr.
-@@ -214,6 +279,13 @@ This part of variables are used to control general system para
+#### exx_separate_loop
+
+- **Type**: Boolean
+- **Description**: There are two types of iterative approach provided by ABACUS to evaluate Fock exchange. If this parameter is set to 0, it will start with a GGA-Loop, and then Hybrid-Loop, in which EXX Hamiltonian $H_{exx}$ is updated with electronic iterations. If this parameter is set to 1, a two-step method is employed, i.e. in the inner iterations, density matrix is updated, while in the outer iterations, $H_{exx}$ is calculated based on density matrix that converges in the inner iteration. 
+- **Default**: 1
+
+#### exx_hybrid_step
 
+- **Type**: Integer
+- **Description**: This variable indicates the maximal electronic iteration number in the evaluation of Fock exchange.
+- **Default**: 100
+
+#### exx_lambda
+
+- **Type**: Real
+- **Description**: It is used to compensate for divergence points at G=0 in the evaluation of Fock exchange using *lcao_in_pw* method.
 - **Default**: 0.3
 
 #### exx_pca_threshold

docs/install.md

@@ -100,7 +100,7 @@ If environment variable `MKLROOT` exists, `cmake` will take MKL as a preference,
 You can also choose to build with which components, e.g.:
 
 ```bash
-cmake -B build -DUSE_LIBXC=1 -DUSE_CUDA=1
+cmake -B build -DENABLE_LIBXC=1 -DUSE_CUDA=1
 ```
 
 If Libxc is not installed in standard path (i.e. installed with a custom prefix path), you can set `Libxc_DIR` to the corresponding directory.

source/module_base/blas_connector.h

@@ -36,9 +36,10 @@ extern "C"
 	double dznrm2_( const int *n, const std::complex<double> *X, const int *incX );
 
 	// level 2: matrix-std::vector operations, O(n^2) data and O(n^2) work.
-	void dgemv_(const char *transa, const int *m, const int *n, const double *alpha,  const double *a,
-		const int *lda, const double *x, const int *incx, const double *beta, double *y, const int *incy);
-
+	void dgemv_(const char*const transa, const int*const m, const int*const n,
+		const double*const alpha, const double*const a, const int*const lda, const double*const x, const int*const incx,
+		const double*const beta, double*const y, const int*const incy);
+		
 	void zgemv_(const char *trans, const int *m, const int *n, const std::complex<double> *alpha,
 			const std::complex<double> *a, const int *lda, const std::complex<double> *x, const int *incx,
 			const std::complex<double> *beta, std::complex<double> *y, const int *incy);

source/module_base/lapack_connector.h

@@ -60,8 +60,8 @@ extern "C"
     // Peize Lin add dsptrf and dsptri 2016-06-21, to compute inverse real symmetry indefinit matrix.
     // dpotrf computes the Cholesky factorization of a real symmetric positive definite matrix
     // while dpotri taks its output to perform matrix inversion
-    void dpotrf_(const char* uplo,const int* n, double* A, const int* lda, int *info);
-    void dpotri_(const char* uplo,const int* n, double* A, const int* lda, int *info);
+    void dpotrf_(const char*const uplo, const int*const n, double*const A, const int*const lda, int*const info);
+    void dpotri_(const char*const uplo, const int*const n, double*const A, const int*const lda, int*const info);
 
     // zgetrf computes the LU factorization of a general matrix
     // while zgetri takes its output to perform matrix inversion
@@ -345,20 +345,34 @@ class LapackConnector
 
 	// Peize Lin add 2016-07-09
 	static inline
-	void dpotrf( char uplo, const int n, ModuleBase::matrix &a, const int lda, int *info )
+	void dpotrf( const char &uplo, const int &n, double*const A, const int &lda, int &info )
 	{
 		const char uplo_changed = change_uplo(uplo);
-		dpotrf_( &uplo_changed, &n, a.c, &lda, info );
-	}
-
+		dpotrf_( &uplo_changed, &n, A, &lda, &info );
+	}	
+	
 	// Peize Lin add 2016-07-09
 	static inline
-	void dpotri( char uplo, const int n, ModuleBase::matrix &a, const int lda, int *info )
+	void dpotri( const char &uplo, const int &n, double*const A, const int &lda, int &info )
 	{
 		const char uplo_changed = change_uplo(uplo);
-		dpotri_( &uplo_changed, &n, a.c, &lda, info);
-	}
+		dpotri_( &uplo_changed, &n, A, &lda, &info);		
+	}	
 
+	// Peize Lin add 2016-07-09
+	static inline
+	void dpotrf( const char &uplo, const int &n, ModuleBase::matrix &A, const int &lda, int &info )
+	{
+		dpotrf( uplo, n, A.c, lda, info );
+	}	
+	
+	// Peize Lin add 2016-07-09
+	static inline
+	void dpotri( const char &uplo, const int &n, ModuleBase::matrix &A, const int &lda, int &info )
+	{
+		dpotri( uplo, n, A.c, lda, info);		
+	}	
+	
 	// Peize Lin add 2019-04-14
 	// if trans=='N':	C = a * A * A.H + b * C
 	// if trans=='C':	C = a * A.H * A + b * C

source/module_base/tool_quit.cpp

@@ -129,24 +129,41 @@ void WARNING_QUIT(const std::string &file,const std::string &description)
 }
 
 
-//Input judgement and communicate , if any judgement is true, do WARNING_QUIT
-void CHECK_WARNING_QUIT(bool error_in, const std::string &file,const std::string &description)
+//Check and print warning information for all cores.
+//Maybe in the future warning.log should be replaced by error.log.
+void CHECK_WARNING_QUIT(const bool error_in, const std::string &file,const std::string &description)
 {
-	int error = (int)error_in;
 #ifdef __NORMAL
 // only for UT, do nothing here
 #else
-#ifdef __MPI
-	int error_max = error;
-	MPI_Reduce(&error, &error_max, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD);
-	MPI_Bcast(&error_max, 1, MPI_INT, 0, MPI_COMM_WORLD);
-	error = error_max;
-#endif
-#endif
-	if(error)
+	if(error_in)
 	{
-		WARNING_QUIT(file, description);
+		//All cores will print inforamtion
+		std::cout.clear();
+		if(!GlobalV::ofs_running.is_open()) 
+		{
+			std::string logfile = GlobalV::global_out_dir + "running_" + GlobalV::CALCULATION + ".log";
+			GlobalV::ofs_running.open( logfile.c_str(), std::ios::app );
+		}
+		if(!GlobalV::ofs_warning.is_open()) 
+		{
+			std::string warningfile = GlobalV::global_out_dir + "warning.log";
+			GlobalV::ofs_warning.open( warningfile.c_str(), std::ios::app );
+		}
+
+		//print error information
+		std::cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
+		std::cout << " ERROR! " << description << std::endl;
+		std::cout << " CHECK IN FILE : " << GlobalV::global_out_dir << "warning.log" << std::endl;
+		std::cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
+		GlobalV::ofs_running << " ERROR! CHECK IN FILE : " << GlobalV::global_out_dir << "warning.log" << std::endl;
+		GlobalV::ofs_warning << std::endl;
+		GlobalV::ofs_warning << " ERROR! " << file << ", core " << GlobalV::MY_RANK+1 << ": " << description << std::endl;
+		GlobalV::ofs_warning << std::endl;
+		exit(0);
 	}
+#endif
+	return;
 }
 
 }

source/module_base/tool_quit.h

@@ -49,7 +49,7 @@ void WARNING_QUIT(const std::string &file, const std::string &description);
  * @param file The file where warning happens
  * @param description The warning information
  */
-void CHECK_WARNING_QUIT(bool error, const std::string &file,const std::string &description);
+void CHECK_WARNING_QUIT(const bool error, const std::string &file,const std::string &description);
 
 } // namespace ModuleBase
 

source/module_hamilt/ks_pw/meta_pw.cpp

@@ -61,11 +61,13 @@ void Meta<OperatorPW>::act
             }
 
             wfcpw->recip2real(porter, porter, ik);
-
-            const double* pvk = &(this->vk[0](current_spin, 0));
-            for (int ir = 0; ir < this->vk->nc; ir++)
+            if(this->vk->nc != 0)
             {
-                porter[ir] *= pvk[ir];
+                const double* pvk = &(this->vk[0](current_spin, 0));
+                for (int ir = 0; ir < this->vk->nc; ir++)
+                {
+                    porter[ir] *= pvk[ir];
+                }
             }
             wfcpw->real2recip(porter, porter, ik);
 

source/module_orbital/ORB_control.cpp

@@ -416,8 +416,7 @@ int ORB_control::mat_2d(MPI_Comm vu,
     // (5) row_set, it's a global index :
     // save explicitly : every row in this processor
     // belongs to which row in the global matrix.
-    delete[] LM.row_set;
-    LM.row_set = new int[LM.row_num];
+    LM.row_set.resize(LM.row_num);
     j = 0;
     for (i = 0; i < LM.row_b; i++)
     {
@@ -475,8 +474,7 @@ int ORB_control::mat_2d(MPI_Comm vu,
 
     if (pv->testpb)ModuleBase::GlobalFunc::OUT(ofs_running, "Local columns (including nb)", LM.row_num);
 
-    delete[] LM.col_set;
-    LM.col_set = new int[LM.col_num];
+    LM.col_set.resize(LM.col_num);
 
     j = 0;
     for (i = 0; i < LM.col_b; i++)

source/module_orbital/ORB_table_phi.cpp

@@ -5,6 +5,10 @@
 #include "../module_base/constants.h"
 #include "../module_base/timer.h"
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 double ORB_table_phi::dr = -1.0;
 
 ORB_table_phi::ORB_table_phi()
@@ -128,9 +132,7 @@ void ORB_table_phi::cal_ST_Phi12_R
 {
 	ModuleBase::timer::tick("ORB_table_phi", "cal_ST_Phi12_R");
 
-	double* k1_dot_k2 = new double[kmesh];
-	double* k1_dot_k2_dot_kpoint = new double[kmesh];
-
+	std::vector<double> k1_dot_k2(kmesh);
 	// Peize Lin change 2017-12-12
 	switch(job)
 	{
@@ -166,6 +168,7 @@ void ORB_table_phi::cal_ST_Phi12_R
 			break;
 	}
 
+	std::vector<double> k1_dot_k2_dot_kpoint(kmesh);
 	for (int ik = 0; ik < kmesh; ik++)
 	{
 		k1_dot_k2_dot_kpoint[ik] = k1_dot_k2[ik] * this->kpoint[ik];
@@ -176,7 +179,7 @@ void ORB_table_phi::cal_ST_Phi12_R
 
 	//previous version
 
-	double* integrated_func = new double[kmesh];
+	//double* integrated_func = new double[kmesh];
 
 	int ll;
 	if(l==0) ll=0;
@@ -186,8 +189,12 @@ void ORB_table_phi::cal_ST_Phi12_R
 	const std::vector<std::vector<double>> &jl = pSB->get_jlx()[l];
 	const std::vector<std::vector<double>> &jlp1 = pSB->get_jlx()[l+1];
 
+#ifdef _OPENMP
+	#pragma omp parallel for schedule(static)
+#endif
 	for (int ir = 0; ir < rmesh; ir++)
 	{
+		std::vector<double> integrated_func(kmesh);
 		const std::vector<double> &jl_r = jl[ir];
 		for (int ik=0; ik<kmesh; ++ik)
 		{
@@ -196,7 +203,7 @@ void ORB_table_phi::cal_ST_Phi12_R
 		// Call simpson integration
 		double temp = 0.0;
 
-		ModuleBase::Integral::Simpson_Integral(kmesh,integrated_func,dk,temp);
+		ModuleBase::Integral::Simpson_Integral(kmesh, integrated_func.data(), dk, temp);
 		rs[ir] = temp * ModuleBase::FOUR_PI ;
 
 		// Peize Lin accelerate 2017-10-02
@@ -218,7 +225,7 @@ void ORB_table_phi::cal_ST_Phi12_R
 			}
 		}
 
-		ModuleBase::Integral::Simpson_Integral(kmesh,integrated_func,dk,temp);
+		ModuleBase::Integral::Simpson_Integral(kmesh, integrated_func.data(), dk, temp);
 		drs[ir] = -ModuleBase::FOUR_PI*(l+1)/(2.0*l+1) * temp;
 	}
 
@@ -228,25 +235,19 @@ void ORB_table_phi::cal_ST_Phi12_R
 
 	if (l > 0)
 	{
-		ModuleBase::GlobalFunc::ZEROS(integrated_func,kmesh);
+		std::vector<double> integrated_func(kmesh);
 		double temp = 0.0;
 
 		for (int ik = 0; ik < kmesh; ik++)
 		{
-			integrated_func[ik] = k1_dot_k2[ik] * pow (kpoint[ik], l);
+			integrated_func[ik] = k1_dot_k2[ik] * std::pow (kpoint[ik], l);
 		}
 
-		ModuleBase::Integral::Simpson_Integral(kmesh,integrated_func,kab,temp);
+		ModuleBase::Integral::Simpson_Integral(kmesh, integrated_func.data(), kab, temp);
 		rs[0] = ModuleBase::FOUR_PI / ModuleBase::Mathzone_Add1::dualfac (2*l+1) * temp;
 	}
 
-	delete [] integrated_func;
-	delete [] k1_dot_k2;
-	delete [] k1_dot_k2_dot_kpoint;	
-
 	ModuleBase::timer::tick("ORB_table_phi", "cal_ST_Phi12_R");
-	
-	return;
 }
 
 #include "../module_base/constants.h"