Merge pull request #1177 from dyzheng/v2.2.3

Merge PR172 PR174 PR175 PR178 from abacusmodeling/develop to deepmodeling/develop

Author: hongriTianqi

.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_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"

source/module_orbital/parallel_orbitals.cpp

@@ -20,8 +20,6 @@ Parallel_Orbitals::Parallel_Orbitals()
     // default value of nb is 1,
     // but can change to larger value from input.
     nb = 1;
-    MatrixInfo.row_set = nullptr;
-    MatrixInfo.col_set = nullptr;
 
     // in multi-k, 2D-block-division variables for FT (R<->k)
     nnr = 1;
@@ -34,18 +32,16 @@ Parallel_Orbitals::~Parallel_Orbitals()
     delete[] trace_loc_row;
     delete[] trace_loc_col;
     delete[] loc_sizes;
-
-    if (alloc_Z_LOC) // xiaohui add 2014-12-22
-    {
-        for (int is = 0; is < this->nspin; is++)
-        {
-            delete[] Z_LOC[is];
-        }
-        delete[] Z_LOC;
-    }
-    delete[] MatrixInfo.row_set;
-    delete[] MatrixInfo.col_set;
-
+    
+    if (alloc_Z_LOC)//xiaohui add 2014-12-22
+	{
+		for(int is=0; is<this->nspin; is++)
+		{
+			delete[] Z_LOC[is];
+		}
+		delete[] Z_LOC;
+	}
+    
     delete[] nlocdim;
     delete[] nlocstart;
 }
@@ -287,24 +283,22 @@ void ORB_control::divide_HS_2d(
     pv->nloc = nlocal * nlocal;
     this->set_parameters(ofs_running, ofs_warning);
     pv->MatrixInfo.row_b = 1;
-    pv->MatrixInfo.row_num = nlocal;
-    delete[] pv->MatrixInfo.row_set;
-    pv->MatrixInfo.row_set = new int[nlocal];
-    for (int i = 0; i < nlocal; i++)
-    {
-        pv->MatrixInfo.row_set[i] = i;
-    }
-    pv->MatrixInfo.row_pos = 0;
-
-    pv->MatrixInfo.col_b = 1;
-    pv->MatrixInfo.col_num = nlocal;
-    delete[] pv->MatrixInfo.col_set;
-    pv->MatrixInfo.col_set = new int[nlocal];
-    for (int i = 0; i < nlocal; i++)
-    {
-        pv->MatrixInfo.col_set[i] = i;
-    }
-    pv->MatrixInfo.col_pos = 0;
+	pv->MatrixInfo.row_num = nlocal;
+	pv->MatrixInfo.row_set.resize(nlocal);
+	for(int i=0; i<nlocal; i++)
+	{
+		pv->MatrixInfo.row_set[i]=i;
+	}
+	pv->MatrixInfo.row_pos=0;
+
+	pv->MatrixInfo.col_b = 1;
+	pv->MatrixInfo.col_num = nlocal;
+	pv->MatrixInfo.col_set.resize(nlocal);
+	for(int i=0; i<nlocal; i++)
+	{
+		pv->MatrixInfo.col_set[i]=i;
+	}
+	pv->MatrixInfo.col_pos=0;
 #endif
 
     assert(pv->nloc > 0);

source/src_pdiag/pdiag_common.h

@@ -6,11 +6,12 @@
 #endif
 #include "../module_base/blas_connector.h"
 #include "../module_base/lapack_connector.h"			// Peize Lin add 2016-08-04
+#include <vector>
 
 struct LocalMatrix
 {
-    int *row_set;
-    int *col_set;
+    std::vector<int> row_set;				// Peize Lin change int* to vector 2022.08.03
+    std::vector<int> col_set;
 
 	int col_num;
     int row_num;

source/src_ri/exx_abfs-inverse_matrix_double.cpp

@@ -35,7 +35,7 @@ void Exx_Abfs::Inverse_Matrix_Double::cal_inverse( const Method &method )
 
 void Exx_Abfs::Inverse_Matrix_Double::using_dpotrf()
 {
-	LapackConnector::dpotrf('U',dim,A,dim,&info);
+	LapackConnector::dpotrf('U',dim,A,dim,info);
 
 	if(info!=0)
 	{
@@ -45,7 +45,7 @@ void Exx_Abfs::Inverse_Matrix_Double::using_dpotrf()
 		ModuleBase::QUIT();
 	}
 
-	LapackConnector::dpotri('U',dim,A,dim,&info);
+	LapackConnector::dpotri('U',dim,A,dim,info);
 
 	if(info!=0)
 	{