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

Author: dyzheng

CMakeLists.txt

@@ -13,7 +13,8 @@ project(ABACUS
 
 option(ENABLE_DEEPKS "Enable DeePKS functionality" OFF)
 option(ENABLE_LIBXC "Enable LibXC functionality" OFF)
-option(USE_CUDA "Enable support to CUDA." OFF)
+option(USE_CUDA "Enable support to CUDA for PW." OFF)
+option(USE_CUSOLVER_LCAO "Enable support to CUSOLVER for LCAO." OFF)
 option(USE_ROCM "Enable support to ROCm." OFF)
 option(USE_OPENMP " Enable OpenMP in abacus." ON)
 option(ENABLE_ASAN "Enable AddressSanitizer" OFF)
@@ -58,21 +59,32 @@ set(CMAKE_CXX_STANDARD 11)
 include(CheckLanguage)
 check_language(CUDA)
 if(CMAKE_CUDA_COMPILER)
-  if(NOT DEFINED USE_CUDA)
-    message("CUDA components detected. \nWill build the CUDA version of ABACUS.")
-    set(USE_CUDA ON)
+  if(NOT DEFINED USE_CUDA OR NOT DEFINED USE_CUSOLVER_LCAO)
+    if (NOT DEFINED USE_CUDA AND NOT DEFINED USE_CUSOLVER_LCAO)
+      message("CUDA components detected. \nWill build the CUDA for PW version of ABACUS by default.")
+      set(USE_CUDA ON)
+      set(USE_CUSOLVER_LCAO OFF)
+    elseif (NOT DEFINED USE_CUDA)
+      set(USE_CUDA OFF)
+    else()
+      set(USE_CUSOLVER_LCAO OFF)
+    endif()
   else()
-    if(NOT USE_CUDA)
-      message(WARNING "CUDA components detected, but USE_CUDA set to OFF. \nNOT building CUDA version of ABACUS.")
+    if(NOT USE_CUDA AND NOT USE_CUSOLVER_LCAO)
+      message(WARNING "CUDA components detected, but both USE_CUDA and USE_CUSOLVER_LCAO set to OFF. \nNOT building CUDA version of ABACUS.")
+    elseif (USE_CUDA AND USE_CUSOLVER_LCAO)
+      message(FATAL_ERROR "USE_CUDA and USE_CUSOLVER_LCAO set, but now they not allowed to coexist.")
     endif()
   endif()
 else() # CUDA not found
-  if (USE_CUDA)
-    message(FATAL_ERROR "USE_CUDA set but no CUDA components found.")
+  if (USE_CUDA OR USE_CUSOLVER_LCAO)
+    message(FATAL_ERROR "USE_CUDA or USE_CUSOLVER_LCAO set but no CUDA components found.")
     set(USE_CUDA OFF)
+    set(USE_CUSOLVER_LCAO OFF)
   endif()
 endif()
-if(USE_CUDA)
+
+if(USE_CUDA OR USE_CUSOLVER_LCAO)
   set(CMAKE_CXX_STANDARD 14)
   set(CMAKE_CXX_EXTENSIONS ON)
   set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -92,11 +104,18 @@ if(USE_CUDA)
   60 # P100
   70 # V100
   75 # T4
+  80 # A100
   )
   include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
-  add_compile_definitions(__CUDA)
+  if (USE_CUDA)
+    add_compile_definitions(__CUDA)
+  endif()
+  if (USE_CUSOLVER_LCAO)
+    add_compile_definitions(__CUSOLVER_LCAO)
+  endif()
 endif()
 
+
 # Warning: CMake add support to HIP in version 3.21. This is rather a new version.
 # Use cmake with AMD-ROCm: https://rocmdocs.amd.com/en/latest/Installation_Guide/Using-CMake-with-AMD-ROCm.html
 if(USE_ROCM)

docs/features.md

@@ -74,7 +74,7 @@ equation. For PW basis, there are CG and Blocked Davidson methods for solving th
 equation for each basis.
 
 - PW: ks_solver = ‘cg’ or ‘dav’
-- LCAO: ks_solver = ‘hpseps’ , ‘genelpa’ or ‘lapack’
+- LCAO: ks_solver = ‘hpseps’ , ‘genelpa’ , ‘scalapack_gvx’ or 'cusolver'
 - LCAO_in_PW: ks_solver = ‘lapack’
 
 If you set ks_solver=‘hpseps’ for basis_type=‘pw’, the program will be stopped with an error

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
@@ -435,6 +441,7 @@ calculations.
   - genelpa: This method should be used if you choose localized orbitals.
   - hpseps: old method, still used.
   - lapack: lapack can be used for localized orbitals, but is only used for single processor.
+  - cusolver: this method needs building with the cusolver component for lcao and at least one gpu is available.
 
    If you set ks_solver=`hpseps` for basis_type=`pw`, the program will be stopped with an error message:
 
@@ -453,6 +460,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 +566,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.

docs/install.md

@@ -75,6 +75,9 @@ You can also choose to build with which components.
 ```bash
 cmake -B build -DUSE_LIBXC=1 -DUSE_CUDA=1
 ```
+```bash
+cmake -B build -DUSE_CUSOLVER_LCAO=1
+```
 
 If Libxc is not installed in standard path (i.e. installed with a custom prefix path), you may add the installation prefix of `FindLibxc.cmake` to `CMAKE_MODULE_PATH` environment variable, or set `Libxc_DIR` to the directory containing the file.
 

examples/performance/clean.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 
 for i in P*;do
-    rm -rf $i/result.log $i/OUT.* $i/result.out
+    rm -rf $i/result.log $i/time.log $i/OUT.* $i/result.out
 done
 
 rm -rf *cpu  *kpar *bxyz sum.dat* *.old log

examples/performance/run.sh

@@ -80,7 +80,7 @@ run_abacus(){
         lastword=`tail -1 result.log | awk '{print $1}'`
     fi
     if [[ $lastword != "SEE" ]];then
-        mpirun -n $1 -env OMP_NUM_THREADS=$2 $abacus > result.log
+        /usr/bin/time -v mpirun -n $1 -env OMP_NUM_THREADS=$2 $abacus > result.log 2>time.log
     else
         printf "**result.log is normal end, skip this job** "
     fi

examples/performance/sumdat.sh

@@ -13,9 +13,9 @@ fi
 test -f $outf && rm $outf
 
 #title
-printf "%20s %7s %8s %8s %6s %6s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s\n" \
+printf "%20s %7s %8s %8s %6s %6s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s\n" \
 "example" "Natoms" "EneCut" "k-points" "NProc" "Niter" "TotTime" "1stSCF" "SCF/iter" \
-"s_c%" "c_bands%" "s_bands%" "h_psi%" "vloc%" "vnl%" "FFT%" "stress%" "force%"> $outf
+"Run%" "c_bands%" "s_bands%" "h_psi%" "vloc%" "vnl%" "FFT%" "stress%" "force%" "MaxResSize" > $outf
 
 for i in `cat $allcase`;do
     if [[ ! -f $i/result.log ]];then
@@ -34,13 +34,13 @@ for i in `cat $allcase`;do
         niter=`sed -n '/ITER   ETOT(eV)/,/><><><><>/'p ${i}/result.log | wc -l|awk '{print $1-2}'`
         tottime=`awk '$1=="total"{printf"%.2f", $2}' ${i}/result.log`
         scf1=`grep -A 1 "ITER   ETOT(eV)" ${i}/result.log | awk 'END{printf"%.2f", $NF}'`
-        totalscf=`awk '$2=="self_consistent"{print $3}' ${i}/result.log`
+        totalscf=`awk '$2=="Run"{print $3}' ${i}/result.log`
         scfpiter=`awk -v a=$totalscf -v b=$scf1 -v c=$niter 'BEGIN{printf"%.2f",(a-b)/(c-1)}'`
         fft=`awk '$2=="FFT3D"{printf"%.1f",$6}' ${i}/result.log`
         hpsi=`awk '$2=="h_psi"{printf"%.1f",$6}' ${i}/result.log`
         vloc=`awk '$2=="vloc"{printf"%.1f",$6}' ${i}/result.log`
         vnl=`awk '$2=="vnl"{printf"%.1f",$6}' ${i}/result.log`
-        sc=`awk '$2=="self_consistent"{printf"%.1f",$6}' ${i}/result.log`
+        sc=`awk '$2=="Run"{printf"%.1f",$6}' ${i}/result.log`
         cbands=`awk '$2=="c_bands"{printf"%.1f",$6}' ${i}/result.log`
         sbands=`awk '$2=="sum_band"{printf"%.1f",$6}' ${i}/result.log`
         stress=`awk '$2=="cal_stress"{printf"%.1f",$6}' ${i}/result.log`
@@ -53,13 +53,13 @@ for i in `cat $allcase`;do
         niter=`sed -n '/ITER   ETOT(eV)/,/><><><><>/'p ${i}/result.log | wc -l|awk '{print $1-2}'`
         tottime=`awk '$1=="total"{printf"%.2f", $2}' ${i}/result.log`
         scf1=`grep -A 1 "ITER   ETOT(eV)" ${i}/result.log | awk 'END{printf"%.2f", $NF}'`
-        totalscf=`awk '$1=="ELEC_scf"{print $3}' ${i}/result.log`
+        totalscf=`awk '$1=="Run"{print $3}' ${i}/result.log`
         scfpiter=`awk -v a=$totalscf -v b=$scf1 -v c=$niter 'BEGIN{printf"%.2f",(a-b)/(c-1)}'`
         fft="-"
         hpsi="-"
         vloc=`awk '$2=="vlocal"{printf"%.1f",$6}' ${i}/result.log`
         vnl="-"
-        sc=`awk '$1=="ELEC_scf"{printf"%.1f",$6}' ${i}/result.log`
+        sc=`awk '$2=="Run"{printf"%.1f",$6}' ${i}/result.log`
         cbands=`awk '$2=="cal_bands"{printf"%.1f",$6}' ${i}/result.log`
         sbands=`awk '$2=="sum_bands"{printf"%.1f",$6}' ${i}/result.log`
         stress=`awk '$2=="evaluate_vl_stress"{printf"%.1f",$6}' ${i}/result.log`
@@ -68,9 +68,10 @@ for i in `cat $allcase`;do
         echo "ERROR: UNKNOW basis type $basis"
         continue
     fi
+    maxres=`grep "Maximum resident set size" ${i}/time.log | awk '{print $NF}'`
 
-    printf "%20s %7s %8s %8s %6s %6s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s\n" \
+    printf "%20s %7s %8s %8s %6s %6s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %s\n" \
     $i $natoms $encut $kpt $nproc $niter $tottime $scf1 $scfpiter $sc $cbands $sbands $hpsi $vloc $vnl $fft \
-    $stress $force >> $outf
+    $stress $force $maxres >> $outf
 
 done

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)
         {
@@ -2484,6 +2485,12 @@ void Input::Check(void)
         {
             ModuleBase::WARNING_QUIT("Input", "not ready for linear_scaling method in lcao .");
         }
+        else if (ks_solver == "cusolver")
+        {
+#ifndef __MPI
+            ModuleBase::WARNING_QUIT("Input","Cusolver can not be used for series version.");
+#endif
+        }
         else
         {
             ModuleBase::WARNING_QUIT("Input", "please check the ks_solver parameter!");
@@ -2665,7 +2672,7 @@ void Input::Check(void)
             if (!(calculation == "nscf"))
                 ModuleBase::WARNING_QUIT("Input", "calculate berry phase, please set calculation = nscf");
         }
-        else if (basis_type == "lcao" && (ks_solver == "genelpa" || ks_solver == "scalapack_gvx"))
+        else if (basis_type == "lcao" && ks_solver == "genelpa" || ks_solver == "scalapack_gvx" || ks_solver == "cusolver")
         {
             if (!(calculation == "nscf"))
                 ModuleBase::WARNING_QUIT("Input", "calculate berry phase, please set calculation = nscf");