Merge branch 'develop' into prbranch

BariumOxide13716 · web-flow · commit 10c69d78cb81 · 2024-11-11T16:24:35.000+08:00
diff --git a/source/module_base/blas_connector.cpp b/source/module_base/blas_connector.cpp
@@ -69,17 +69,17 @@ float BlasConnector::dot( const int n, const float *X, const int incX, const flo
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		return sdot_(&n, X, &incX, Y, &incY);
+	}
 	return sdot_(&n, X, &incX, Y, &incY);
 }
-}
 
 double BlasConnector::dot( const int n, const double *X, const int incX, const double *Y, const int incY, base_device::AbacusDevice_t device_type)
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		return ddot_(&n, X, &incX, Y, &incY);
+	}
 	return ddot_(&n, X, &incX, Y, &incY);
 }
-}
 
 // C = a * A.? * B.? + b * C
 void BlasConnector::gemm(const char transa, const char transb, const int m, const int n, const int k,
@@ -196,39 +196,39 @@ float BlasConnector::nrm2( const int n, const float *X, const int incX, base_dev
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		return snrm2_( &n, X, &incX );
+	}
 	return snrm2_( &n, X, &incX );
 }
-}
 
 
 double BlasConnector::nrm2( const int n, const double *X, const int incX, base_device::AbacusDevice_t device_type )
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		return dnrm2_( &n, X, &incX );
+	}
 	return dnrm2_( &n, X, &incX );
 }
-}
 
 
 double BlasConnector::nrm2( const int n, const std::complex<double> *X, const int incX, base_device::AbacusDevice_t device_type )
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		return dznrm2_( &n, X, &incX );
+	}
 	return dznrm2_( &n, X, &incX );
 }
-}
 
 // copies a into b
 void BlasConnector::copy(const long n, const double *a, const int incx, double *b, const int incy, base_device::AbacusDevice_t device_type)
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		dcopy_(&n, a, &incx, b, &incy);
-}
+	}
 }
 
 void BlasConnector::copy(const long n, const std::complex<double> *a, const int incx, std::complex<double> *b, const int incy, base_device::AbacusDevice_t device_type)
 {
 	if (device_type == base_device::AbacusDevice_t::CpuDevice) {
 		zcopy_(&n, a, &incx, b, &incy);
-}
+	}
 }
diff --git a/source/module_base/lapack_connector.h b/source/module_base/lapack_connector.h
@@ -61,10 +61,10 @@ extern "C"
 
 
     void dsygvx_(const int* itype, const char* jobz, const char* range, const char* uplo,
-                const int* n, double* A, const int* lda, double* B, const int* ldb,
-                const double* vl, const double* vu, const int* il, const int* iu,
-                const double* abstol, const int* m, double* w, double* Z, const int* ldz,
-                double* work, int* lwork, int*iwork, int* ifail, int* info);
+        const int* n, double* A, const int* lda, double* B, const int* ldb,
+        const double* vl, const double* vu, const int* il, const int* iu,
+        const double* abstol, const int* m, double* w, double* Z, const int* ldz,
+        double* work, const int* lwork, int* iwork, int* ifail, int* info);
 
     void chegvx_(const int* itype,const char* jobz,const char* range,const char* uplo,
              const int* n,std::complex<float> *a,const int* lda,std::complex<float> *b,
diff --git a/source/module_base/lapack_wrapper.h b/source/module_base/lapack_wrapper.h
@@ -1,6 +1,6 @@
 #ifndef LAPACK_HPP
 #define LAPACK_HPP
-
+#include <iostream>
 extern "C"
 {
     // =================================================================================
@@ -49,10 +49,10 @@ extern "C"
     // =================================================================================
     // gvx
     void dsygvx_(const int* itype, const char* jobz, const char* range, const char* uplo,
-                const int* n, double* A, const int* lda, double* B, const int* ldb,
-                const double* vl, const double* vu, const int* il, const int* iu,
-                const double* abstol, const int* m, double* w, double* Z, const int* ldz,
-                double* work, int* lwork, int*iwork, int* ifail, int* info);
+        const int* n, double* A, const int* lda, double* B, const int* ldb,
+        const double* vl, const double* vu, const int* il, const int* iu,
+        const double* abstol, const int* m, double* w, double* Z, const int* ldz,
+        double* work, const int* lwork, int* iwork, int* ifail, int* info);
 
     void chegvx_(const int* itype,const char* jobz,const char* range,const char* uplo,
              const int* n,std::complex<float> *a,const int* lda,std::complex<float> *b,
@@ -424,8 +424,8 @@ class LapackWrapper
                               int* ifail,
                               int& info)
     {
-        // dsygvx_(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl,
-        //         &vu, &il,&iu, &abstol, &m, w, z, &ldz, work, &lwork, rwork, iwork, ifail, &info);
+        dsygvx_(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl,
+            &vu, &il, &iu, &abstol, &m, w, z, &ldz, work, &lwork, iwork, ifail, &info);
     }
 
     // wrap function of fortran lapack routine xhegvx ( pointer version ).
diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -1708,11 +1708,6 @@ bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho, con
     ModuleBase::TITLE("Charge_Mixing", "if_scf_oscillate");
     ModuleBase::timer::tick("Charge_Mixing", "if_scf_oscillate");
 
-    if(threshold >= 0) // close the function
-    {
-        return false;
-    }
-
     if(this->_drho_history.size() == 0)
     {
         this->_drho_history.resize(PARAM.inp.scf_nmax);
@@ -1721,8 +1716,13 @@ bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho, con
     // add drho into history
     this->_drho_history[iteration - 1] = drho;
 
+    if(threshold >= 0) // close the function
+    {
+        return false;
+    }
+
     // check if the history is long enough
-    if(iteration < iternum_used)
+    if(iteration < iternum_used + this->mixing_restart_last)
     {
         return false;
     }
diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h
@@ -102,8 +102,9 @@ class Charge_Mixing
     Base_Mixing::Mixing* get_mixing() const {return mixing;}
 
     // for mixing restart
-    int mixing_restart_step = 0; //which step to restart mixing during SCF
+    int mixing_restart_step = 0; //which step to restart mixing during SCF, always equal to scf_namx except for the mixing restart
     int mixing_restart_count = 0; // the number of restart mixing during SCF. Do not set mixing_restart_count as bool since I want to keep some flexibility in the future
+    int mixing_restart_last = 0; // the label of mixing restart step, store the step number of the last mixing restart
 
     // to calculate the slope of drho curve during SCF, which is used to determine if SCF oscillate
     bool if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold);
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
@@ -681,6 +681,7 @@ void ESolver_KS<T, Device>::iter_finish(const int istep, int& iter)
     // notice for restart
     if (PARAM.inp.mixing_restart > 0 && iter == this->p_chgmix->mixing_restart_step - 1 && iter != PARAM.inp.scf_nmax)
     {
+        this->p_chgmix->mixing_restart_last = iter;
         std::cout << " SCF restart after this step!" << std::endl;
     }
 }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp
@@ -47,9 +47,6 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(Grid_Driver* Grid
     ModuleBase::TITLE("DFTU", "initialize_HR");
     ModuleBase::timer::tick("DFTU", "initialize_HR");
 
-    auto* paraV = this->hR->get_paraV();// get parallel orbitals from HR
-    // TODO: if paraV is nullptr, AtomPair can not use paraV for constructor, I will repair it in the future.
-
     this->adjs_all.clear();
     this->adjs_all.reserve(this->ucell->nat);
     for (int iat0 = 0; iat0 < ucell->nat; iat0++)
@@ -58,8 +55,9 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(Grid_Driver* Grid
         int T0, I0;
         ucell->iat2iait(iat0, &I0, &T0);
         const int target_L = this->dftu->orbital_corr[T0];
-        if (target_L == -1)
+        if (target_L == -1) {
             continue;
+}
 
         AdjacentAtomInfo adjs;
         GridD->Find_atom(*ucell, tau0, T0, I0, &adjs);
@@ -92,8 +90,9 @@ template <typename TK, typename TR>
 void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_nlm_all(const Parallel_Orbitals* paraV)
 {
     ModuleBase::TITLE("DFTU", "cal_nlm_all");
-    if (this->precal_nlm_done)
+    if (this->precal_nlm_done) {
         return;
+}
     ModuleBase::timer::tick("DFTU", "cal_nlm_all");
     nlm_tot.resize(this->ucell->nat);
     const int npol = this->ucell->get_npol();
@@ -104,8 +103,9 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_nlm_all(const Parallel_Orbi
         int T0, I0;
         ucell->iat2iait(iat0, &I0, &T0);
         const int target_L = this->dftu->orbital_corr[T0];
-        if (target_L == -1)
+        if (target_L == -1) {
             continue;
+}
         const int tlp1 = 2 * target_L + 1;
         AdjacentAtomInfo& adjs = this->adjs_all[atom_index++];
 
@@ -144,7 +144,7 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_nlm_all(const Parallel_Orbi
                 const int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2;
 
                 ModuleBase::Vector3<double> dtau = tau0 - tau1;
-                intor_->snap(T1, L1, N1, M1, T0, dtau * this->ucell->lat0, 0 /*cal_deri*/, nlm);
+                intor_->snap(T1, L1, N1, M1, T0, dtau * this->ucell->lat0, false /*cal_deri*/, nlm);
                 // select the elements of nlm with target_L
                 for (int iw = 0; iw < this->ucell->atoms[T0].nw; iw++)
                 {
@@ -178,8 +178,9 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
     else
     {
         // will update this->dftu->locale and this->dftu->EU
-        if (this->current_spin == 0)
+        if (this->current_spin == 0) {
             this->dftu->EU = 0.0;
+}
     }
     ModuleBase::timer::tick("DFTU", "contributeHR");
 
@@ -196,8 +197,9 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
         int T0, I0;
         ucell->iat2iait(iat0, &I0, &T0);
         const int target_L = this->dftu->orbital_corr[T0];
-        if (target_L == -1)
+        if (target_L == -1) {
             continue;
+}
         const int tlp1 = 2 * target_L + 1;
         AdjacentAtomInfo& adjs = this->adjs_all[atom_index++];
 
@@ -241,8 +243,9 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
             // save occ to dftu
             for (int i = 0; i < occ.size(); i++)
             {
-                if (this->nspin == 1)
+                if (this->nspin == 1) {
                     occ[i] *= 0.5;
+}
                 this->dftu->locale[iat0][target_L][0][this->current_spin].c[i] = occ[i];
             }
         }
@@ -297,17 +300,20 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
     }
 
     // energy correction for NSPIN=1
-    if (this->nspin == 1)
+    if (this->nspin == 1) {
         this->dftu->EU *= 2.0;
+}
     // for readin onsite_dm, set initialed_locale to false to avoid using readin locale in next iteration
-    if (this->current_spin == this->nspin - 1 || this->nspin == 4)
+    if (this->current_spin == this->nspin - 1 || this->nspin == 4) {
         this->dftu->initialed_locale = false;
+}
 
     // update this->current_spin: only nspin=2 iterate change it between 0 and 1
     // the key point is only nspin=2 calculate spin-up and spin-down separately,
     // and won't calculate spin-up twice without spin-down
-    if (this->nspin == 2)
+    if (this->nspin == 2) {
         this->current_spin = 1 - this->current_spin;
+}
 
     ModuleBase::timer::tick("DFTU", "contributeHR");
 }
@@ -488,7 +494,7 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_v_of_u(const std::vector<do
 {
     // calculate the local matrix
     int spin_fold = occ.size() / m_size / m_size;
-    if (spin_fold < 4)
+    if (spin_fold < 4) {
         for (int is = 0; is < spin_fold; ++is)
         {
             int start = is * m_size * m_size;
@@ -501,7 +507,7 @@ void hamilt::DFTU<hamilt::OperatorLCAO<TK, TR>>::cal_v_of_u(const std::vector<do
                 }
             }
         }
-    else
+    } else
     {
         for (int m1 = 0; m1 < m_size; m1++)
         {
diff --git a/source/module_hsolver/kernels/cuda/math_kernel_op.cu b/source/module_hsolver/kernels/cuda/math_kernel_op.cu
@@ -12,7 +12,7 @@
 namespace hsolver
 {
 const int warp_size = 32;
-const unsigned int full_mask = 0xffffffff;
+// const unsigned int full_mask = 0xffffffff;
 const int thread_per_block = 256;
 }
 
@@ -65,11 +65,11 @@ void destoryBLAShandle(){
     }
 }
 
-template <typename FPTYPE>
-__forceinline__ __device__ void warp_reduce(FPTYPE& val) {
-    for (int offset = 16; offset > 0; offset >>= 1)
-        val += __shfl_down_sync(full_mask, val, offset);
-}
+// template <typename FPTYPE>
+// __forceinline__ __device__ void warp_reduce(FPTYPE& val) {
+//     for (int offset = 16; offset > 0; offset >>= 1)
+//         val += __shfl_down_sync(full_mask, val, offset);
+// }
 
 template <typename Real>
 __global__ void line_minimize_with_block(

Original file line number	Diff line number	Diff line change
`@@ -69,17 +69,17 @@ float BlasConnector::dot( const int n, const float *X, const int incX, const flo`
`69`	`69`	`{`
`70`	`70`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`71`	`71`	`return sdot_(&n, X, &incX, Y, &incY);`
	`72`	`+ }`
`72`	`73`	`return sdot_(&n, X, &incX, Y, &incY);`
`73`	`74`	`}`
`74`		`-}`
`75`	`75`
`76`	`76`	`double BlasConnector::dot( const int n, const double X, const int incX, const double Y, const int incY, base_device::AbacusDevice_t device_type)`
`77`	`77`	`{`
`78`	`78`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`79`	`79`	`return ddot_(&n, X, &incX, Y, &incY);`
	`80`	`+ }`
`80`	`81`	`return ddot_(&n, X, &incX, Y, &incY);`
`81`	`82`	`}`
`82`		`-}`
`83`	`83`
`84`	`84`	`// C = a * A.? * B.? + b * C`
`85`	`85`	`void BlasConnector::gemm(const char transa, const char transb, const int m, const int n, const int k,`
`@@ -196,39 +196,39 @@ float BlasConnector::nrm2( const int n, const float *X, const int incX, base_dev`
`196`	`196`	`{`
`197`	`197`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`198`	`198`	`return snrm2_( &n, X, &incX );`
	`199`	`+ }`
`199`	`200`	`return snrm2_( &n, X, &incX );`
`200`	`201`	`}`
`201`		`-}`
`202`	`202`
`203`	`203`
`204`	`204`	`double BlasConnector::nrm2( const int n, const double *X, const int incX, base_device::AbacusDevice_t device_type )`
`205`	`205`	`{`
`206`	`206`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`207`	`207`	`return dnrm2_( &n, X, &incX );`
	`208`	`+ }`
`208`	`209`	`return dnrm2_( &n, X, &incX );`
`209`	`210`	`}`
`210`		`-}`
`211`	`211`
`212`	`212`
`213`	`213`	`double BlasConnector::nrm2( const int n, const std::complex<double> *X, const int incX, base_device::AbacusDevice_t device_type )`
`214`	`214`	`{`
`215`	`215`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`216`	`216`	`return dznrm2_( &n, X, &incX );`
	`217`	`+ }`
`217`	`218`	`return dznrm2_( &n, X, &incX );`
`218`	`219`	`}`
`219`		`-}`
`220`	`220`
`221`	`221`	`// copies a into b`
`222`	`222`	`void BlasConnector::copy(const long n, const double a, const int incx, double b, const int incy, base_device::AbacusDevice_t device_type)`
`223`	`223`	`{`
`224`	`224`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`225`	`225`	`dcopy_(&n, a, &incx, b, &incy);`
`226`		`-}`
	`226`	`+ }`
`227`	`227`	`}`
`228`	`228`
`229`	`229`	`void BlasConnector::copy(const long n, const std::complex<double> a, const int incx, std::complex<double> b, const int incy, base_device::AbacusDevice_t device_type)`
`230`	`230`	`{`
`231`	`231`	`if (device_type == base_device::AbacusDevice_t::CpuDevice) {`
`232`	`232`	`zcopy_(&n, a, &incx, b, &incy);`
`233`		`-}`
	`233`	`+ }`
`234`	`234`	`}`
Original file line number	Diff line number	Diff line change
`@@ -681,6 +681,7 @@ void ESolver_KS<T, Device>::iter_finish(const int istep, int& iter)`
`681`	`681`	`// notice for restart`
`682`	`682`	`if (PARAM.inp.mixing_restart > 0 && iter == this->p_chgmix->mixing_restart_step - 1 && iter != PARAM.inp.scf_nmax)`
`683`	`683`	`{`
	`684`	`+ this->p_chgmix->mixing_restart_last = iter;`
`684`	`685`	`std::cout << " SCF restart after this step!" << std::endl;`
`685`	`686`	`}`
`686`	`687`	`}`
Original file line number	Diff line number	Diff line change
`@@ -12,7 +12,7 @@`
`12`	`12`	`namespace hsolver`
`13`	`13`	`{`
`14`	`14`	`const int warp_size = 32;`
`15`		`-const unsigned int full_mask = 0xffffffff;`
	`15`	`+// const unsigned int full_mask = 0xffffffff;`
`16`	`16`	`const int thread_per_block = 256;`
`17`	`17`	`}`
`18`	`18`
`@@ -65,11 +65,11 @@ void destoryBLAShandle(){`
`65`	`65`	`}`
`66`	`66`	`}`
`67`	`67`
`68`		`-template <typename FPTYPE>`
`69`		`-__forceinline__ __device__ void warp_reduce(FPTYPE& val) {`
`70`		`- for (int offset = 16; offset > 0; offset >>= 1)`
`71`		`- val += __shfl_down_sync(full_mask, val, offset);`
`72`		`-}`
	`68`	`+// template <typename FPTYPE>`
	`69`	`+// __forceinline__ __device__ void warp_reduce(FPTYPE& val) {`
	`70`	`+// for (int offset = 16; offset > 0; offset >>= 1)`
	`71`	`+// val += __shfl_down_sync(full_mask, val, offset);`
	`72`	`+// }`
`73`	`73`
`74`	`74`	`template <typename Real>`
`75`	`75`	`__global__ void line_minimize_with_block(`