Add a new para esolver_type to specify energy solver (#1475)

* refactor: add a para esolver_type

* docs: update docs and examples

* test: update test INPUT

* refactor: update para esolver_type

* refactor: convert Hartree to Rydberg in lj/dp

docs/advanced/input_files/input-main.md

@@ -2,7 +2,7 @@
 
 - [System variables](#system-variables)
 
-  [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) | [xc_temperature](#xc_temperature) | [pseudo_rcut](#pseudo_rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem_saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline) | [kspacing](#kspacing)  | [min_dist_coef](#min_dist_coef) | [symmetry_prec](#symmetry_prec)
+  [suffix](#suffix) | [ntype](#ntype) | [calculation](#calculation) | [esolver_type](#esolver_type) | [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) | [xc_temperature](#xc_temperature) | [pseudo_rcut](#pseudo_rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem_saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline) | [kspacing](#kspacing)  | [min_dist_coef](#min_dist_coef) | [symmetry_prec](#symmetry_prec)
 - [Variables related to input files](#variables-related-to-input-files)
 
   [stru_file](#stru_file) | [kpoint_file](#kpoint_file) | [pseudo_dir](#pseudo_dir) | [orbital_dir](#orbital_dir) | [read_file_dir](#read_file_dir) | [wannier_card](#wannier_card)
@@ -33,7 +33,7 @@
 
 - [Molecular dynamics](#molecular-dynamics)
 
-  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_ensolver](#md_ensolver) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_tfreq](#md_tfreq) | [md_tchain](#md_tchain) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
+  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_tfreq](#md_tfreq) | [md_tchain](#md_tchain) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
 
 - [vdW correction](#vdw-correction)
 
@@ -43,7 +43,7 @@
   [berry_phase](#berry_phase) | [gdir](#gdir) | [towannier90](#towannier90) | [nnkpfile](#nnkpfile) | [wannier_spin](#wannier_spin)
 - [TDDFT: time dependent density functional theory](#tddft-time-dependent-density-functional-theory) (Under tests)
 
-  [tddft](#tddft) | [td_scf_thr](#td_scf_thr) | [td_dt](#td_dt) | [td_force_dt](#td_force_dt) | [td_vext](#td_vext) | [td_vext_dire](#td_vext_dire) | [td_timescale](#td_timescale) | [td_vexttype](#td_vexttype) | [td_vextout](#td_vextout) | [td_dipoleout](#td_dipoleout) | [ocp](#ocp) | [ocp_set](#ocp_set) | [td_val_elec_01](#td_val_elec_01) | [td_val_elec_02](#td_val_elec_02) |[td_val_elec_03](#td_val_elec_03)
+  [td_scf_thr](#td_scf_thr) | [td_dt](#td_dt) | [td_force_dt](#td_force_dt) | [td_vext](#td_vext) | [td_vext_dire](#td_vext_dire) | [td_timescale](#td_timescale) | [td_vexttype](#td_vexttype) | [td_vextout](#td_vextout) | [td_dipoleout](#td_dipoleout) | [ocp](#ocp) | [ocp_set](#ocp_set) | [td_val_elec_01](#td_val_elec_01) | [td_val_elec_02](#td_val_elec_02) |[td_val_elec_03](#td_val_elec_03)
 - [DFT+*U* correction](#dftu-correction) (Under development)
 
   [dft_plus_u](#dft_plus_u) | [orbital_corr](#orbital_corr) | [hubbard_u](#hubbard_u) | [hund_j](#hund_j) | [yukawa_potential](#yukawa_potential) | [yukawa_lambda](#yukawa_lambda) | [omc](#omc)
@@ -99,16 +99,24 @@ This part of variables are used to control general system parameters.
   - *istate*: For LCAO basis. Please see the explanation for variable `nbands_istate`.
   - *ienvelope*: Envelope function for LCAO basis. Please see the explanation for variable `nbands_istate`.
   - *md*: molecular dynamics
-  - *ofdft*: do orbital-free dft calculation.
-  - *of-md*: do molecular dynamics with OFDFT.
-  - *sto-scf*: do self-consistent electronic structure calculation with [stochastic DFT](#electronic-structure-sdft)
-  - *sto-md*: molecular dynamics with [stochastic DFT](#electronic-structure-sdft)
   - *test_memory* : checks memory required for the calculation. The number is not quite reliable, please use with care
   - *test_neighbour* : only performs neighbouring atom search
   - *gen_bessel* : generates projectors (a series of bessel functions) for DeePKS; see also keywords bessel_lmax, bessel_rcut and bessel_tol. A file named `jle.orb` will be generated which contains the projectors. An example is provided in examples/H2O-deepks-pw.
   - *get_S* : only works for multi-k calculation with LCAO basis. Generates and writes the overlap matrix to a file names `SR.csr` in the working directory. The format of the file will be the same as that generated by [out_mat_hs2](#out_mat_hs2).
 - **Default**: scf
 
+### esolver_type
+
+- **Type**: String
+- **Description**: choose the energy solver.
+  - ksdft: Kokn-Sham density functional theory;
+  - ofdft: orbital-free density functional theory;
+  - sdft: [stochastic density functional theory](#electronic-structure-sdft);
+  - tddft: real time time dependent density functional theory (TDDFT);
+  - lj: Leonard Jones potential;
+  - dp: DeeP potential;
+- **Default**: ksdft
+
 ### symmetry
 
 - **Type**: Integer
@@ -575,7 +583,7 @@ 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). To use it,  [calculation](#calculation) need to be set to "sto-scf" and "sto-md". We suggest using SDFT to calculate high-temperature systems and we only support [smearing_method](#smearing_method) "fd".
+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 suggest using SDFT to calculate high-temperature systems and we only support [smearing_method](#smearing_method) "fd".
 
 ### method_sto
 
@@ -1417,15 +1425,6 @@ This part of variables are used to control the molecular dynamics calculations.
 - **Description**: the total number of md steps.
 - **Default**: 10
 
-### md_ensolver
-
-- **Type**: String
-- **Description**: choose the energy solver for MD.
-  - FP: First-Principles MD;
-  - LJ: Leonard Jones potential;
-  - DP: DeeP potential;
-- **Default**: FP
-
 ### md_restart
 
 - **Type**: Boolean
@@ -1781,14 +1780,6 @@ This part of variables are used to control berry phase and wannier90 interfacae
 
 ## TDDFT: time dependent density functional theory
 
-### tddft
-
-- **Type**: Integer
-- **Description**:
-  - 1: calculate the real time time dependent density functional theory (TDDFT).
-  - 0: do not calculate TDDFT.
-- **Default**: 0
-
 ### td_scf_thr
 
 - **Type**: Real

docs/advanced/md.md

@@ -2,7 +2,7 @@
 
 Molecular dynamics (MD) is a computer simulation method for analyzing the physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic "evolution" of the system. In the most common version, the trajectories of atoms and molecules are determined by numerically solving Newton's equations of motion for a system of interacting particles, where forces between the particles and their potential energies are calculated using first-principles calculations (first-principles molecular dynamics, FPMD), or interatomic potentials and molecular mechanics force fields (classical molecular dynamics, CMD).
 
-By setting `calculation` to be `md` or `sto-md`, ABACUS currently provides six different MD evolution methods, which is specified by keyword `md_type` in the `INPUT` file:
+By setting `calculation` to be `md`, ABACUS currently provides six different MD evolution methods, which is specified by keyword `md_type` in the `INPUT` file:
 
   - -1: FIRE method
   - 0: velocity Verlet algorithm (default: NVE ensemble)
@@ -20,7 +20,7 @@ When `md_type` is set to 0, `md_thermostat` is used to specify the thermostat ba
 
 Furthermore, ABACUS also provides a [list of keywords](./input_files/input-main.md#molecular-dynamics) to control relevant parmeters used in MD simulations.
 
-To employ CMD calculations, `md_ensolver` should be set to be `LJ` or `DP`.
+To employ CMD calculations, `esolver_type` should be set to be `lj` or `dp`.
 If DP model is selected, the filename of DP model is specified by keyword `pot_file`.
 
 [Examples](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/md/lcao_gammaonly_Sn64) of MD simulations are also provided.

docs/advanced/scf/hsolver.md

@@ -9,7 +9,7 @@ When "basis_type = pw", `ks_solver` can be `cg` or `dav`. The default setting `c
 When "basis_type = lcao", `ks_solver` can be `genelpa` or `scalapack_gvx`. The default setting `genelpa` is recommended, which is based on ELPA (EIGENVALUE SOLVERS FOR PETAFLOP APPLICATIONS) (https://elpa.mpcdf.mpg.de/) and the kernel is auto choosed by GENELPA(https://github.com/pplab/GenELPA), usually faster than the setting of "scalapack_gvx", which is based on ScaLAPACK(Scalable Linear Algebra PACKage)  
 
 ## Stochasic DFT
-We support stochastic DFT calculation (SDFT) or mixed stochastic-deterministic DFT (MDFT) with plane-wave basis [[Phys. Rev. B 106, 125132 (2022)](https://doi.org/10.1103/PhysRevB.106.125132)]. Different from traditional KSDFT with the explicit diagonalization method, SDFT and MDFT calculate physical quantities with trace of the corresponding operators. The advantages of SDFT and MDFT compared to the traditional KSDFT are the ability to simulate larger sizes and higher temperatures. In our package, SDFT and MDFT can be used by setting the `calculation` parameter to `sto-scf` or `sto-md` for SCF calculations or MD calculations. To start with, you can refer to two [examples](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/stochastic) and an explanation of the [input variables](../input_files/input-main.md#electronic-structure-sdft).
+We support stochastic DFT calculation (SDFT) or mixed stochastic-deterministic DFT (MDFT) with plane-wave basis [[Phys. Rev. B 106, 125132 (2022)](https://doi.org/10.1103/PhysRevB.106.125132)]. Different from traditional KSDFT with the explicit diagonalization method, SDFT and MDFT calculate physical quantities with trace of the corresponding operators. The advantages of SDFT and MDFT compared to the traditional KSDFT are the ability to simulate larger sizes and higher temperatures. In our package, SDFT and MDFT can be used by setting the `esolver_type` parameter to `sdft` for SCF calculations or MD calculations. To start with, you can refer to two [examples](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/stochastic) and an explanation of the [input variables](../input_files/input-main.md#electronic-structure-sdft).
 
 When we have a hamiltonian, the electronic density can be calculated with:
 

examples/md/README

@@ -5,7 +5,7 @@
 /*******************************************************************************/
 
 These examples show how to do a md calculation.
-One can set 'calculation' to be 'md' or 'sto-md' to do the MD calculation.
+One can set 'calculation' to be 'md' to do the MD calculation.
 
 More detail introduction can be found in ../../docs/advanced/md.md
 

examples/md/lcao_gammaonly_Sn64/README

@@ -4,7 +4,7 @@
 /
 /*******************************************************************************/
 
-By setting `calculation` to be `md` or `sto-md`, ABACUS currently provides five different algorithms to integrate the equation of motion, which is specified by keyword `md_type` in the `INPUT` file:
+By setting `calculation` to be `md`, ABACUS currently provides five different algorithms to integrate the equation of motion, which is specified by keyword `md_type` in the `INPUT` file:
 
   - -1: FIRE method
   - 0: velocity Verlet algorithm (default: NVE ensemble)

examples/stochastic/pw_Si2/INPUT

@@ -1,6 +1,7 @@
 INPUT_PARAMETERS
 #Parameters	(General)
-calculation     sto-scf
+calculation     scf
+esolver_type    sdft
 pseudo_dir              ../../../tests/PP_ORB
 ntype			1	
 nbands			4

examples/stochastic/pw_md_Al/INPUT

@@ -1,6 +1,7 @@
 INPUT_PARAMETERS
 #Parameters	(General)
-calculation     sto-md
+calculation     md
+esolver_type    sdft
 pseudo_dir      ../../../tests/PP_ORB
 ntype			1	
 nbands			0

source/driver_run.cpp

@@ -49,23 +49,27 @@ void Driver::driver_run()
 
     //------------------------------------------------------------
     // This part onward needs to be refactored.
-    //---------------------------MD/Relax------------------
-    if(GlobalV::CALCULATION == "md" && GlobalV::BASIS_TYPE=="lcao")
+    //---------------------------MD/Relax-------------------------
+    if(GlobalV::CALCULATION == "md")
     {
+        // In the future, I will universalize pw, lcao and classic line
+        if(GlobalV::ESOLVER_TYPE == "lj" || GlobalV::ESOLVER_TYPE == "dp")
+        {
+            Run_MD_CLASSIC run_md_classic;
+            run_md_classic.classic_md_line(GlobalC::ucell, p_esolver);
+        }
+        else if(GlobalV::BASIS_TYPE == "lcao")
+        {
 #ifdef __LCAO
-        Run_MD_LCAO run_md_lcao;
-        run_md_lcao.opt_ions(p_esolver);
+            Run_MD_LCAO run_md_lcao;
+            run_md_lcao.opt_ions(p_esolver);
 #endif
-    }
-    else if(INPUT.mdp.md_ensolver == "LJ" || INPUT.mdp.md_ensolver == "DP")
-    {
-        Run_MD_CLASSIC run_md_classic;
-        run_md_classic.classic_md_line(GlobalC::ucell, p_esolver);
-    }
-    else if(GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "sto-md" || GlobalV::CALCULATION == "of-md")
-    {
-        Run_MD_PW run_md_pw;
-        run_md_pw.md_ions_pw(p_esolver);
+        }
+        else
+        {
+            Run_MD_PW run_md_pw;
+            run_md_pw.md_ions_pw(p_esolver);
+        }
     }
     else // scf; cell relaxation; nscf; etc
     {

source/input.cpp

@@ -134,6 +134,7 @@ void Input::Default(void)
     // xiaohui modify 2015-09-15, relax -> scf
     // calculation = "relax";
     calculation = "scf";
+    esolver_type = "ksdft";
     pseudo_rcut = 15.0; // qianrui add this parameter 2021-5
     pseudo_mesh = false; // qianrui add this pararmeter
     ntype = 0;
@@ -383,7 +384,6 @@ void Input::Default(void)
     //----------------------------------------------------------			//Fuxiang He add 2016-10-26
     // tddft
     //----------------------------------------------------------
-    tddft = 0;
     td_scf_thr = 1e-9;
     td_dt = 0.02;
     td_force_dt = 0.02;
@@ -573,6 +573,10 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, calculation);
         }
+        else if (strcmp("esolver_type", word) == 0)
+        {
+            read_value(ifs, esolver_type);
+        }
         else if (strcmp("ntype", word) == 0) // number of atom types
         {
             read_value(ifs, ntype);
@@ -1257,10 +1261,6 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, mdp.lj_sigma);
         }
-        else if (strcmp("md_ensolver", word) == 0)
-        {
-            read_value(ifs, mdp.md_ensolver);
-        }
         else if (strcmp("msst_direction", word) == 0)
         {
             read_value(ifs, mdp.msst_direction);
@@ -1357,10 +1357,6 @@ bool Input::Read(const std::string &fn)
         // tddft
         // Fuxiang He add 2016-10-26
         //----------------------------------------------------------
-        else if (strcmp("tddft", word) == 0)
-        {
-            read_value(ifs, tddft);
-        }
         else if (strcmp("td_scf_thr", word) == 0)
         {
             read_value(ifs, td_scf_thr);
@@ -2086,7 +2082,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
             vdw_cutoff_radius = "95";
         }
     }
-    if(calculation.substr(0,3) != "sto")    bndpar = 1;
+    if(esolver_type != "sdft")    bndpar = 1;
     if(bndpar > GlobalV::NPROC) bndpar = GlobalV::NPROC;
     if(method_sto != 1 && method_sto != 2) 
     {
@@ -2114,6 +2110,7 @@ void Input::Bcast()
     Parallel_Common::bcast_string(wannier_card);
     Parallel_Common::bcast_string(latname);
     Parallel_Common::bcast_string(calculation);
+    Parallel_Common::bcast_string(esolver_type);
     Parallel_Common::bcast_double(pseudo_rcut);
     Parallel_Common::bcast_bool(pseudo_mesh);
     Parallel_Common::bcast_int(ntype);
@@ -2300,7 +2297,6 @@ void Input::Bcast()
     Parallel_Common::bcast_double(mdp.lj_rcut);
     Parallel_Common::bcast_double(mdp.lj_epsilon);
     Parallel_Common::bcast_double(mdp.lj_sigma);
-    Parallel_Common::bcast_string(mdp.md_ensolver);
     Parallel_Common::bcast_int(mdp.msst_direction);
     Parallel_Common::bcast_double(mdp.msst_vel);
     Parallel_Common::bcast_double(mdp.msst_vis);
@@ -2360,7 +2356,6 @@ void Input::Bcast()
     Parallel_Common::bcast_int(vdw_cutoff_period.y);
     Parallel_Common::bcast_int(vdw_cutoff_period.z);
     // Fuxiang He add 2016-10-26
-    Parallel_Common::bcast_int(tddft);
     Parallel_Common::bcast_int(td_val_elec_01);
     Parallel_Common::bcast_int(td_val_elec_02);
     Parallel_Common::bcast_int(td_val_elec_03);
@@ -2522,7 +2517,7 @@ void Input::Check(void)
     //----------------------------------------------------------
     // main parameters / electrons / spin ( 1/16 )
     //----------------------------------------------------------
-    if (calculation == "scf" || calculation == "ofdft")
+    if (calculation == "scf")
     {
         if (mem_saver == 1)
         {
@@ -2543,15 +2538,6 @@ void Input::Check(void)
         */
         this->relax_nmax = 1;
     }
-    else if (calculation == "sto-scf") // qianrui 2021-2-20
-    {
-        if (mem_saver == 1)
-        {
-            mem_saver = 0;
-            ModuleBase::GlobalFunc::AUTO_SET("mem_savre", "0");
-        }
-        this->relax_nmax = 1;
-    }
     else if (calculation == "relax") // pengfei 2014-10-13
     {
         if (mem_saver == 1)
@@ -2631,9 +2617,9 @@ void Input::Check(void)
             ModuleBase::WARNING_QUIT("Input::Check", "calculate = istate is only availble for LCAO.");
         }
     }
-    else if (calculation == "md" || calculation == "sto-md" || calculation == "of-md") // mohan add 2011-11-04
+    else if (calculation == "md") // mohan add 2011-11-04
     {
-        GlobalV::CALCULATION = calculation;
+        GlobalV::CALCULATION = "md";
         symmetry = 0;
         cal_force = 1;
         if (mdp.md_nstep == 0)
@@ -2647,7 +2633,7 @@ void Input::Check(void)
         // deal with input parameters , 2019-04-30
         if (mdp.md_dt < 0)
             ModuleBase::WARNING_QUIT("Input::Check", "time interval of MD calculation should be set!");
-        if (mdp.md_tfirst < 0 && tddft==0)
+        if (mdp.md_tfirst < 0 && esolver_type != "tddft")
             ModuleBase::WARNING_QUIT("Input::Check", "temperature of MD calculation should be set!");
         if (mdp.md_tlast < 0.0)
             mdp.md_tlast = mdp.md_tfirst;
@@ -2668,7 +2654,7 @@ void Input::Check(void)
         {
             init_vel = 1;
         }
-        if(mdp.md_ensolver == "LJ" || mdp.md_ensolver == "DP" || mdp.md_type == 4)
+        if(esolver_type == "lj" || esolver_type == "dp" || mdp.md_type == 4)
         {
             cal_stress = 1;
         }
@@ -2679,7 +2665,7 @@ void Input::Check(void)
                 ModuleBase::WARNING_QUIT("Input::Check", "msst_qmass must be greater than 0!");
             }
         }
-        if(mdp.md_ensolver == "DP")
+        if(esolver_type == "dp")
         {
             if (access(mdp.pot_file.c_str(), 0) == -1)
             {