Feature: long-range corrected hybrid functionals implementation with Spencer-Alavi truncation Coulomb potential (#6351)

* Use template to reconstruct parse_expression

* Feature: output R matrix at each MD step

* Modify'matrix_HS' to 'matrix' for R matrix output

* Merge branches 'develop' and 'develop' of https://github.com/1041176461/abacus-develop into develop

* Fix: modify index in parse_expression

* Fix: add regfree for parse_expression

* Doc: update phonopy doc

* Doc: update phonopy doc

* fix tdos plot for nspin=2

* optimize dosplot for nspin=2

* fix legend for dosplot

* Add files via upload

* Update cal_edm_tddft.cpp

* Refactor: modify exx loop for relax/md

* Update result.ref

* Fix wrong parameters in integrate test

* Update result.ref

* Update result.ref

* Update result.ref

* Update result.ref

* Update result.ref

* Update result.ref

* Update result.ref

* Update result.ref

* add exx when istep>0

* Update dos.py

* Update esolver_sdft_pw.cpp

* Update lcao_before_scf.cpp

* Update Exx_LRI_interface.h

* Update Exx_LRI_interface.hpp

* Fix: compile error

* Fix: compile error

* Fix: change HSE relax/md result.ref for new framework

* Fix: compile error

* compatible with exx_iter_finish

* Add files via upload

* Update esolver_ks_lcao_tddft.h

* Update esolver_ks_lcao_tddft.cpp

* Fix: support negative value in parse_expression

* [pre-commit.ci lite] apply automatic fixes

* Fix: parse_expression for scientific notation

* Update input_conv.h

* add complex erf function

* refactor ccp framework

* fix compile

* fix compile

* fix compile

* fix compile

* fix compile

* fix compile

* fix compile

* fix empty error

* add revised spencer

* add cal_psi_erfc_spencer

* add exx_singularity_correction

* add exx_singularity_correction test

* add lc-corrected hybrid functionals

* fix compile

* fix compile

* fix compile error

* fix libxc error

* fix libxc error

* fix libxc error

* add klist.h in RI_Util.h

* fix include

* fix ptest for hse

* fix lc_pbe name for vdw autoset

* fix test error

* fix cam_pbeh test

* fix fock settings

* support rsh functionals

* fix libxc param

* add input

* fix input

---------

Co-authored-by: jiyuang <[email protected]>
Co-authored-by: Qianrui <[email protected]>
Co-authored-by: HTZhao <[email protected]>
Co-authored-by: maki49 <[email protected]>
Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

Author: 6 people

docs/advanced/input_files/input-main.md

@@ -2769,37 +2769,41 @@ Support for hybrid functionals in the *pw [basis_type](#basis_type)* is under ac
 
 The following parameters apply to *[basis_type](#basis_type)==lcao/lcao_in_pw/pw*. For basis specific parameters, see the sections *[Exact Exchange (LCAO/LCAO in PW)](#exact-exchange-lcaolcao-in-pw)* and *[Exact Exchange (PW)](#exact-exchange-pw)*.
 
-**Availablity**: *[dft_functional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*. 
+#### Hybrid Functional Parameters {#hybrid_func_params}
+|[dft_functional](#dft_functional)|[exx_fock_alpha](#exx_fock_alpha)|[exx_erfc_alpha](#exx_erfc_alpha)|[exx_erfc_omega](#exx_erfc_omega)|
+|:-:|:-:|:-:|:-:|
+|hf|1|0|0|
+|lc_pbe|1|-1|0.33|
+|lc_wpbe|1|-1|0.4|
+|lrc_wpbe|1|-1|0.3|
+|lrc_wpbeh|1|-0.8|0.2|
+|pbe0|0.25|0|0|
+|hse|0.0|0.25|0.11|
+|scan0|0.25|0|0|
+|cam_pbeh|0.2|0.8|0.7|
+|b3lyp|0.2|0|0|
+|muller|1|0|0|
+|power|1|0|0|
+|wp22|1|-1|0.11|
+|cwp22|0|1|0.11|
 
 ### exx_fock_alpha
 
 - **Type**: Real \[Real...\](optional)
-- **Description**: Fraction of Fock exchange 1/r in hybrid functionals, so that $E_{X} = \alpha E_{X} + (1-\alpha)E_{X,\text{LDA/GGA}}$
-- **Default**:
-  - 1: if *[dft_functional](#dft_functional)==hf*
-  - 0.25: if *[dft_functional](#dft_functional)==pbe0*
-  - 0.2: if *[dft_functional](#dft_functional)==b3lyp*
-  - 0.25: if *[dft_functional](#dft_functional)==scan0*
-  - 1: if *[dft_functional](#dft_functional)==muller*
-  - 1: if *[dft_functional](#dft_functional)==power*
-  - 1: if *[dft_functional](#dft_functional)==wp22*
-  - 0: else
+- **Description**: Fraction of full-ranged Fock exchange 1/r ($\alpha$) in range-separated hybrid funtionals, so that $E_{X} = \alpha E_{X}^\text{HF-LR}+(\alpha+\beta) E_{X}^\text{HF-SR}+(1-\alpha)E_{X}^\text{KS-LR}+[1-(\alpha+\beta)]E_{X}^\text{KS-SR}$.
+- **Default**: see [hybrid_func_params](#hybrid_func_params)
 
 ### exx_erfc_alpha
 
 - **Type**: Real \[Real...\](optional)
-- **Description**: Fraction of exchange erfc(wr)/r in hybrid functionals, so that $E_{X} = \alpha E_{X}^{\text{SR}} + (1-\alpha)E_{X,\text{LDA/GGA}}^{\text{SR}} + E_{X,\text{LDA/GGA}}^{\text{LR}}$
-- **Default**:
-  - 0.25: if *[dft_functional](#dft_functional)==hse*
-  - 1: if *[dft_functional](#dft_functional)==cwp22*
-  - -1: if *[dft_functional](#dft_functional)==wp22*
-  - 0: else
+- **Description**: Fraction of short-ranged Fock exchange erfc(wr)/r ($\beta$) in range-separated hybrid funtionals, so that $E_{X} = \alpha E_{X}^\text{HF-LR}+(\alpha+\beta) E_{X}^\text{HF-SR}+(1-\alpha)E_{X}^\text{KS-LR}+[1-(\alpha+\beta)]E_{X}^\text{KS-SR}$.
+- **Default**: see [hybrid_func_params](#hybrid_func_params)
 
 ### exx_erfc_omega
 
 - **Type**: Real \[Real...\](optional)
 - **Description**: Range-separation parameter in exchange, such that $1/r=\text{erfc}(\omega r)/r+\text{erf}(\omega r)/r$
-- **Default**: 0.11
+- **Default**: see [hybrid_func_params](#hybrid_func_params)
 
 ### exx_separate_loop
 
@@ -2934,6 +2938,15 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
   - True: if gamma_only
   - False: else
 
+### exx_singularity_correction
+
+- **Type**: String
+- **Description**:
+  - spencer: see Phys. Rev. B 77, 193110 (2008).
+  - revised_spencer: see Phys. Rev. Mater. 5, 013807 (2021).
+  Set the scheme of Coulomb singularity correction.
+- **Default**: default
+
 ### rpa_ccp_rmesh_times
 
 - **Type**: Real

source/module_parameter/input_parameter.h

@@ -523,6 +523,7 @@ struct Input_para
     std::vector<std::string> exx_erfc_omega = {"default"};      ///< range-separation parameter in HSE functional
     bool exx_separate_loop = true;               ///< if 1, a two-step method is employed, else it will start
                                                  ///< with a GGA-Loop, and then Hybrid-Loop
+    std::string exx_singularity_correction = "default";    ///< set the scheme of Coulomb singularity correction
     int exx_hybrid_step = 100;                   ///< the maximal electronic iteration number in
                                                  ///< the evaluation of Fock exchange
     double exx_mixing_beta = 1.0;                ///< mixing_beta for outer-loop when exx_separate_loop=1

source/source_hamilt/module_vdw/vdwd3_autoset_xcname.cpp

@@ -229,7 +229,7 @@ const std::map<std::string, std::string> xcname_libxc_xc_ = {
     {"XC_HYB_GGA_XC_WHPBE0", "whpbe0"},
     {"XC_HYB_GGA_XC_LC_BLYP_EA", "lc_blyp_ea"},
     {"XC_HYB_GGA_XC_LC_BOP", "lc_bop"},
-    {"XC_HYB_GGA_XC_LC_PBEOP", "lc_pbeop"},
+    {"XC_HYB_GGA_XC_LC_PBEOP", "lc_pbe"},
     {"XC_HYB_GGA_XC_LC_BLYPR", "lc_blypr"},
     {"XC_HYB_GGA_XC_MCAM_B3LYP", "mcam_b3lyp"},
     {"XC_MGGA_XC_VCML_RVV10", "vcml_rvv10"},

source/source_hamilt/module_xc/exx_info.h

@@ -18,14 +18,14 @@ struct Exx_Info
 
 		// Fock:
 		//		"alpha":		"0"
-		//		"Rcut_type":	"limits" / "spencer"
+		//		"singularity_correction":	"limits" / "spencer" / "revised_spencer"
 		//		"lambda":		"0.3"
-		//		//"Rcut"
+        //      "Rcut"
 		// Erfc:
 		//		"alpha":		"0"
 		//		"omega":		"0.11"
-		//		"Rcut_type":	"limits"
-		//		//"Rcut"
+		//		"singularity_correction":	"limits" / "spencer" / "revised_spencer"
+        //      "Rcut"
 
         Conv_Coulomb_Pot_K::Ccp_Type ccp_type;
         double hybrid_alpha = 0.25;
@@ -51,10 +51,7 @@ struct Exx_Info
 
     struct Exx_Info_RI
     {
-        std::map<Conv_Coulomb_Pot_K::Coulomb_Method, 
-            std::pair<bool, 
-                std::map<Conv_Coulomb_Pot_K::Coulomb_Type, 
-                    std::vector<std::map<std::string,std::string>>>>> coulomb_settings;
+        const std::map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> &coulomb_param;
 
         bool real_number = false;
 
@@ -72,10 +69,15 @@ struct Exx_Info
         double kmesh_times = 4;
 
         int abfs_Lmax = 0; // tmp
+
+        Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
+            : coulomb_param(info_global.coulomb_param)
+        {
+        }
     };
     Exx_Info_RI info_ri;
 
-    Exx_Info() : info_lip(this->info_global)
+    Exx_Info() : info_lip(this->info_global), info_ri(this->info_global)
     {
     }
 };

source/source_hamilt/module_xc/xc_functional.cpp

@@ -29,9 +29,10 @@ void XC_Functional::set_xc_first_loop(const UnitCell& ucell)
 the first scf iteration only calculate the functional without exact
 exchange. but in "nscf" calculation, there is no need of "two-level"
 method. */
-    if (ucell.atoms[0].ncpp.xc_func == "HF"
-        || ucell.atoms[0].ncpp.xc_func == "PBE0"
-        || ucell.atoms[0].ncpp.xc_func == "HSE") {
+    if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "HSE" 
+        || ucell.atoms[0].ncpp.xc_func == "PBE0"|| ucell.atoms[0].ncpp.xc_func == "LC_PBE" 
+        || ucell.atoms[0].ncpp.xc_func == "LC_WPBE" || ucell.atoms[0].ncpp.xc_func == "LRC_WPBEH" 
+        || ucell.atoms[0].ncpp.xc_func == "CAM_PBEH") {
         XC_Functional::set_xc_type("pbe");
     }
     else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") {
@@ -149,6 +150,36 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         func_type = 5;
         use_libxc = true;
 	}
+    else if( xc_func == "LC_PBE")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP);
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "LC_WPBE")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_LC_WPBE);
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "LRC_WPBE")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_LRC_WPBE);
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "LRC_WPBEH")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_LRC_WPBEH);
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "CAM_PBEH")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_CAM_PBEH);
+        func_type = 4;
+        use_libxc = true;
+    }
 #endif
     else if( xc_func == "HF")
     {
@@ -271,9 +302,11 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
 
 #ifndef USE_LIBXC
     if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" 
-        || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22")
+        || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22" ||
+        xc_func == "LC_PBE" || xc_func == "LC_WPBE" || xc_func == "LRC_WPBE" ||
+        xc_func == "LRC_PBEH" || xc_func == "CAM_PBEH")
     {
-        ModuleBase::WARNING_QUIT("set_xc_type","to use SCAN, SCAN0, or HSE, LIBXC is required");
+        ModuleBase::WARNING_QUIT("set_xc_type","to use SCAN, SCAN0, HSE, long-range corrected (LC_PBE, LC_WPBE...) or CAM_PBEH LIBXC is required");
     }
     use_libxc = false;
 #endif

source/source_hamilt/module_xc/xc_functional_libxc.cpp

@@ -175,6 +175,50 @@ const std::vector<double> in_built_xc_func_ext_params(const int id)
 		case XC_GGA_C_LYPR:
 			return {0.04918, 0.132, 0.2533, 0.349, 
 					0.35/2.29, 2.0/2.29, GlobalC::exx_info.info_global.hse_omega};
+		// Long-range corrected functionals:
+        case XC_HYB_GGA_XC_LC_PBEOP:  // LC version of PBE
+		{
+            // This is a range-separated hybrid functional with range-separation constant  0.330,
+            // and  0.0% short-range and 100.0% long-range exact exchange,
+            // using the error function kernel.
+			return { GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.33
+		}
+        case XC_HYB_GGA_XC_LC_WPBE:  // Long-range corrected PBE (LC-wPBE) by Vydrov and Scuseria
+		{
+            // This is a range-separated hybrid functional with range-separation constant  0.400,
+            // and  0.0% short-range and 100.0% long-range exact exchange,
+            // using the error function kernel.
+			return { std::stod(PARAM.inp.exx_fock_alpha[0]),  //Fraction of Hartree-Fock exchange: 1.0
+				std::stod(PARAM.inp.exx_erfc_alpha[0]),  //Fraction of short-range exact exchange: -1.0
+				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.4
+		}
+        case XC_HYB_GGA_XC_LRC_WPBE:  // Long-range corrected PBE (LRC-wPBE) by by Rohrdanz, Martins and Herbert
+		{
+            // This is a range-separated hybrid functional with range-separation constant  0.300,
+            // and  0.0% short-range and 100.0% long-range exact exchange,
+            // using the error function kernel.
+			return { std::stod(PARAM.inp.exx_fock_alpha[0]),  //Fraction of Hartree-Fock exchange: 1.0
+				std::stod(PARAM.inp.exx_erfc_alpha[0]),  //Fraction of short-range exact exchange: -1.0
+				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.3
+		}
+        case XC_HYB_GGA_XC_LRC_WPBEH:  // Long-range corrected short-range hybrid PBE (LRC-wPBEh) by Rohrdanz, Martins and Herbert
+		{
+            // This is a range-separated hybrid functional with range-separation constant  0.200,
+            // and 20.0% short-range and 100.0% long-range exact exchange,
+            // using the error function kernel.	
+			return { std::stod(PARAM.inp.exx_fock_alpha[0]),  //Fraction of Hartree-Fock exchange: 1.0
+				std::stod(PARAM.inp.exx_erfc_alpha[0]),  //Fraction of short-range exact exchange: -0.8
+				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.2
+		}
+        case XC_HYB_GGA_XC_CAM_PBEH:  // CAM hybrid screened exchange PBE version
+		{
+            // This is a range-separated hybrid functional with range-separation constant  0.700,
+            // and 100.0% short-range and 20.0% long-range exact exchange,
+            // using the error function kernel.
+			return { std::stod(PARAM.inp.exx_fock_alpha[0]),  //Fraction of Hartree-Fock exchange: 0.2
+				std::stod(PARAM.inp.exx_erfc_alpha[0]),  //Fraction of short-range exact exchange: 0.8
+				GlobalC::exx_info.info_global.hse_omega }; //Range separation constant: 0.7
+		}
 #endif
 		default:
 			return std::vector<double>{};

source/source_io/input_conv.cpp

@@ -297,7 +297,12 @@ void Input_Conv::Convert()
         if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0"
             || dft_functional_lower == "hse"
             || dft_functional_lower == "opt_orb"
-            || dft_functional_lower == "scan0") {
+            || dft_functional_lower == "scan0"
+            || dft_functional_lower == "lc_pbe"
+            || dft_functional_lower == "lc_wpbe" 
+            || dft_functional_lower == "lrc_wpbe"
+            || dft_functional_lower == "lrc_wpbeh"
+            || dft_functional_lower == "cam_pbeh") {
             GlobalC::restart.info_load.load_charge = true;
             GlobalC::restart.info_load.load_H = true;
         }
@@ -322,10 +327,15 @@ void Input_Conv::Convert()
                    dft_functional_lower.begin(),
                    tolower);
     if (dft_functional_lower == "hf"
-     || dft_functional_lower == "pbe0" || dft_functional_lower == "b3lyp" || dft_functional_lower == "hse"
-     || dft_functional_lower == "scan0"
-     || dft_functional_lower == "muller" || dft_functional_lower == "power"
-     || dft_functional_lower == "cwp22" || dft_functional_lower == "wp22")
+    || dft_functional_lower == "pbe0" || dft_functional_lower == "b3lyp" || dft_functional_lower == "hse"
+    || dft_functional_lower == "scan0"
+    || dft_functional_lower == "muller" || dft_functional_lower == "power"
+    || dft_functional_lower == "cwp22" || dft_functional_lower == "wp22" 
+    || dft_functional_lower == "lc_pbe"
+    || dft_functional_lower == "lc_wpbe" 
+    || dft_functional_lower == "lrc_wpbe"
+    || dft_functional_lower == "lrc_wpbeh"
+    || dft_functional_lower == "cam_pbeh")
     {
         GlobalC::exx_info.info_global.cal_exx = true;
 
@@ -352,15 +362,13 @@ void Input_Conv::Convert()
         {
             if(PARAM.inp.basis_type == "lcao")
             {
-                std::map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
-                coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock].resize(fock_alpha.size());
+                GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock].resize(fock_alpha.size());
                 for(std::size_t i=0; i<fock_alpha.size(); ++i)
                 {
-                    coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock][i] = {{
+                    GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
                         {"alpha", ModuleBase::GlobalFunc::TO_STRING(fock_alpha[i])},
-                        {"Rcut_type", "spencer"} }};
+                        {"singularity_correction", PARAM.inp.exx_singularity_correction} }};
                 }
-                GlobalC::exx_info.info_ri.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
             }
             else if(PARAM.inp.basis_type == "lcao_in_pw")
             {
@@ -392,16 +400,14 @@ void Input_Conv::Convert()
             assert(erfc_alpha.size() == PARAM.inp.exx_erfc_omega.size());
             if(PARAM.inp.basis_type == "lcao")
             {
-                std::map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
-                coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc].resize(erfc_alpha.size());
+                GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc].resize(erfc_alpha.size());
                 for(std::size_t i=0; i<erfc_alpha.size(); ++i)
                 {
-                    coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
+                    GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
                         {"alpha", ModuleBase::GlobalFunc::TO_STRING(erfc_alpha[i])},
                         {"omega", ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.exx_erfc_omega[i])},
-                        {"Rcut_type", "limits"} }};
+                        {"singularity_correction", PARAM.inp.exx_singularity_correction} }};
                 }
-                GlobalC::exx_info.info_ri.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
             }
         }
     }