diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 6acbd732b8..38dc562bc2 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -243,8 +243,8 @@ - [exx\_opt\_orb\_ecut](#exx_opt_orb_ecut) - [exx\_opt\_orb\_tolerence](#exx_opt_orb_tolerence) - [exx\_real\_number](#exx_real_number) - - [exx\_symmetry\_realspace](#exx_symmetry_realspace) - [rpa\_ccp\_rmesh\_times](#rpa_ccp_rmesh_times) + - [exx\_symmetry\_realspace](#exx_symmetry_realspace) - [out\_ri\_cv](#out_ri_cv) - [Molecular dynamics](#molecular-dynamics) - [md\_type](#md_type) @@ -273,6 +273,9 @@ - [lj\_epsilon](#lj_epsilon) - [lj\_sigma](#lj_sigma) - [pot\_file](#pot_file) + - [dp\_rescaling](#dp_rescaling) + - [dp\_fparam](#dp_fparam) + - [dp\_aparam](#dp_aparam) - [msst\_direction](#msst_direction) - [msst\_vel](#msst_vel) - [msst\_vis](#msst_vis) @@ -422,11 +425,12 @@ - [nocc](#nocc) - [nvirt](#nvirt) - [lr\_nstates](#lr_nstates) + - [lr\_unrestricted](#lr_unrestricted) - [abs\_wavelen\_range](#abs_wavelen_range) - [out\_wfc\_lr](#out_wfc_lr) - [abs\_broadening](#abs_broadening) - [ri\_hartree\_benchmark](#ri_hartree_benchmark) - - [aims_nbasis](#aims_nbasis) + - [aims\_nbasis](#aims_nbasis) [back to top](#full-list-of-input-keywords) ## System variables @@ -2908,46 +2912,38 @@ These variables are used to control vdW-corrected related parameters. - **Type**: String - **Description**: Specifies the method used for Van der Waals (VdW) correction. Available options are: - `d2`: [Grimme's D2](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.20495) dispersion correction method - - `d3_0`: [Grimme's DFT-D3(0)](https://aip.scitation.org/doi/10.1063/1.3382344) dispersion correction method - - `d3_bj`: [Grimme's DFTD3(BJ)](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21759) dispersion correction method + - `d3_0`: [Grimme's DFT-D3(0)](https://aip.scitation.org/doi/10.1063/1.3382344) dispersion correction method (zero-damping) + - `d3_bj`: [Grimme's DFTD3(BJ)](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21759) dispersion correction method (BJ-damping) - `none`: no vdW correction - **Default**: none +- **Note**: ABACUS supports automatic setting on DFT-D3 parameters for common functionals after version 3.8.3 (and several develop versions earlier). To benefit from this feature, please specify the parameter `dft_functional` explicitly (for more details on this parameter, please see [dft_functional](#dft_functional)), otherwise the autoset procedure will crash with error message like `cannot find DFT-D3 parameter for XC(***)`. If not satisfied with those in-built parameters, any manually setting on `vdw_s6`, `vdw_s8`, `vdw_a1` and `vdw_a2` will overwrite. +- **Special**: There are special cases for functional family wB97 (Omega-B97): if want to use the functional wB97X-D3BJ, one needs to specify the `dft_functional` as `HYB_GGA_WB97X_V` and `vdw_method` as `d3_bj`. If want to use the functional wB97X-D3, specify `dft_functional` as `HYB_GGA_WB97X_D3` and `vdw_method` as `d3_0`. ### vdw_s6 - **Type**: Real - **Availability**: `vdw_method` is set to `d2`, `d3_0`, or `d3_bj` -- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE. +- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). If not set, will use values of PBE functional. For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value. - **Default**: - 0.75: if `vdw_method` is set to `d2` - - 1.0: if `vdw_method` is set to `d3_0` or `d3_bj` ### vdw_s8 - **Type**: Real - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj` -- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE. -- **Default**: - - 0.722: if `vdw_method` is set to `d3_0` - - 0.7875: if `vdw_method` is set to `d3_bj` +- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value. ### vdw_a1 - **Type**: Real - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj` -- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE. -- **Default**: - - 1.217: if `vdw_method` is set to `d3_0` - - 0.4289: if `vdw_method` is set to `d3_bj` +- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value. ### vdw_a2 - **Type**: Real - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj` -- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). The default value of this parameter in ABACUS is set to be the recommended value for PBE. -- **Default**: - - 1.0: if `vdw_method` is set to `d3_0` - - 4.4407: if `vdw_method` is set to `d3_bj` +- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value. ### vdw_d diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 4933766cfe..a245a44358 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -210,7 +210,7 @@ void ESolver_KS_PW::before_scf(const int istep) //---------------------------------------------------------- // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); + auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 3fe23490a5..b50722c883 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -198,7 +198,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) //---------------------------------------------------------- // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); + auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); diff --git a/source/module_hamilt_general/module_vdw/dftd3_xc_name.h b/source/module_hamilt_general/module_vdw/dftd3_xc_name.h new file mode 100644 index 0000000000..f2f1280bc8 --- /dev/null +++ b/source/module_hamilt_general/module_vdw/dftd3_xc_name.h @@ -0,0 +1,610 @@ +#ifndef DFTD3_XC_NAME_H +#define DFTD3_XC_NAME_H +/** + * Intro + * ----- + * This file stores the mapping from LibXC xcname to the "conventional" + * + * XCNotSupportedError + * ------------------- + * GGA_X_REVSSB_D + * GGA_X_SSB_D + * + * in J. Chem. Phys. 131, 094103 2009, a simplified version of PBC (the + * correlation part of PBE XC) is used as the correlation part, but libXC + * does not directly support one named as + * GGA_C_SPBEC. + * + * Certainly, those XC with dispersion correction in form of non-local + * correlation are not supported. Such as: + * + * vdw-DF family nonlocal dispersion correction included are not supported: + * GGA_X_OPTB86B_VDW + * GGA_X_OPTB88_VDW + * GGA_X_OPTPBE_VDW + * GGA_X_PBEK1_VDW + * + * VV09, VV10 and rVV10 nonlocal correlation included are not supported: + * GGA_XC_VV10 + * HYB_GGA_XC_LC_VV10 + * HYB_MGGA_XC_WB97M_V + * HYB_GGA_XC_WB97X_V + * MGGA_X_VCML + * MGGA_C_REVSCAN_VV10 + * MGGA_C_SCAN_VV10 + * MGGA_C_SCANL_VV10 + * MGGA_XC_B97M_V + * MGGA_XC_VCML_RVV10 + * + * There is also one quite special, the wB97X-D3BJ functional uses the + * wB97X-V functionals excluding the VV10 part, then use its own DFT-D3(BJ) + * parameters. This seems not recorded in simple-dftd3, so it is not supported + * temporarily: + * HYB_GGA_XC_WB97X_D3BJ + * HYB_GGA_XC_WB97X_V + */ +#include +#include +#include +#include "module_base/formatter.h" +#include +#include +#include "module_base/tool_quit.h" + +namespace DFTD3 { + const std::map xcname_libxc_xc_ = { + {"XC_LDA_XC_TETER93", "teter93"}, + {"XC_LDA_XC_ZLP", "zlp"}, + {"XC_MGGA_XC_OTPSS_D", "otpss_d"}, // DFT-D2 + {"XC_GGA_XC_OPBE_D", "opbe_d"}, // DFT-D2 + {"XC_GGA_XC_OPWLYP_D", "opwlyp_d"}, // DFT-D2 + {"XC_GGA_XC_OBLYP_D", "oblyp_d"}, // DFT-D2 + {"XC_GGA_XC_HCTH_407P", "hcth_407p"}, + {"XC_GGA_XC_HCTH_P76", "hcth_p76"}, + {"XC_GGA_XC_HCTH_P14", "hcth_p14"}, + {"XC_GGA_XC_B97_GGA1", "b97_gga1"}, + {"XC_GGA_XC_KT2", "kt2"}, + {"XC_GGA_XC_TH1", "th1"}, + {"XC_GGA_XC_TH2", "th2"}, + {"XC_GGA_XC_TH3", "th3"}, + {"XC_GGA_XC_TH4", "th4"}, + {"XC_GGA_XC_HCTH_93", "hcth_93"}, + {"XC_GGA_XC_HCTH_120", "hcth_120"}, + {"XC_GGA_XC_HCTH_147", "hcth_147"}, + {"XC_GGA_XC_HCTH_407", "hcth_407"}, + {"XC_GGA_XC_EDF1", "edf1"}, + {"XC_GGA_XC_XLYP", "xlyp"}, + {"XC_GGA_XC_KT1", "kt1"}, + {"XC_GGA_XC_B97_D", "b97_d"}, // DFT-D2? + {"XC_GGA_XC_PBE1W", "pbe1w"}, + {"XC_GGA_XC_MPWLYP1W", "mpwlyp1w"}, + {"XC_GGA_XC_PBELYP1W", "pbelyp1w"}, + {"XC_HYB_LDA_XC_LDA0", "lda0"}, + {"XC_HYB_LDA_XC_CAM_LDA0", "cam_lda0"}, + {"XC_GGA_XC_NCAP", "ncap"}, + {"XC_GGA_XC_MOHLYP", "mohlyp"}, + {"XC_GGA_XC_MOHLYP2", "mohlyp2"}, + {"XC_GGA_XC_TH_FL", "th_fl"}, + {"XC_GGA_XC_TH_FC", "th_fc"}, + {"XC_GGA_XC_TH_FCFO", "th_fcfo"}, + {"XC_GGA_XC_TH_FCO", "th_fco"}, + {"XC_MGGA_XC_CC06", "cc06"}, + {"XC_MGGA_XC_TPSSLYP1W", "tpsslyp1w"}, + {"XC_MGGA_XC_B97M_V", "b97m_v"}, + {"XC_GGA_XC_VV10", "vv10"}, + {"XC_LDA_XC_KSDT", "ksdt"}, + {"XC_HYB_GGA_XC_B97_1P", "b97_1p"}, + {"XC_HYB_GGA_XC_PBE_MOL0", "pbe_mol0"}, + {"XC_HYB_GGA_XC_PBE_SOL0", "pbe_sol0"}, + {"XC_HYB_GGA_XC_PBEB0", "pbeb0"}, + {"XC_HYB_GGA_XC_PBE_MOLB0", "pbe_molb0"}, + {"XC_GGA_XC_BEEFVDW", "beefvdw"}, + {"XC_MGGA_XC_HLE17", "hle17"}, + {"XC_HYB_GGA_XC_PBE50", "pbe50"}, + {"XC_HYB_GGA_XC_HFLYP", "hflyp"}, + {"XC_HYB_GGA_XC_B3P86_NWCHEM", "b3p86_nwchem"}, + {"XC_LDA_XC_CORRKSDT", "corrksdt"}, + {"XC_HYB_GGA_XC_RELPBE0", "relpbe0"}, + {"XC_GGA_XC_B97_3C", "b97_3c"}, + {"XC_HYB_MGGA_XC_BR3P86", "br3p86"}, + {"XC_HYB_GGA_XC_CASE21", "case21"}, + {"XC_HYB_GGA_XC_PBE_2X", "pbe_2x"}, + {"XC_HYB_GGA_XC_PBE38", "pbe38"}, + {"XC_HYB_GGA_XC_B3LYP3", "b3lyp3"}, + {"XC_HYB_GGA_XC_CAM_O3LYP", "cam_o3lyp"}, + {"XC_HYB_MGGA_XC_TPSS0", "tpss0"}, + {"XC_HYB_MGGA_XC_B94_HYB", "b94_hyb"}, + {"XC_HYB_GGA_XC_WB97X_D3", "wb97x_d3"}, // DFT-D3(0) + {"XC_HYB_GGA_XC_LC_BLYP", "lc_blyp"}, + {"XC_HYB_GGA_XC_B3PW91", "b3pw91"}, + {"XC_HYB_GGA_XC_B3LYP", "b3lyp"}, + {"XC_HYB_GGA_XC_B3P86", "b3p86"}, + {"XC_HYB_GGA_XC_O3LYP", "o3lyp"}, + {"XC_HYB_GGA_XC_MPW1K", "mpw1k"}, + {"XC_HYB_GGA_XC_PBEH", "pbeh"}, + {"XC_HYB_GGA_XC_B97", "b97"}, + {"XC_HYB_GGA_XC_B97_1", "b97_1"}, + {"XC_HYB_GGA_XC_APF", "apf"}, + {"XC_HYB_GGA_XC_B97_2", "b97_2"}, + {"XC_HYB_GGA_XC_X3LYP", "x3lyp"}, + {"XC_HYB_GGA_XC_B1WC", "b1wc"}, + {"XC_HYB_GGA_XC_B97_K", "b97_k"}, + {"XC_HYB_GGA_XC_B97_3", "b97_3"}, + {"XC_HYB_GGA_XC_MPW3PW", "mpw3pw"}, + {"XC_HYB_GGA_XC_B1LYP", "b1lyp"}, + {"XC_HYB_GGA_XC_B1PW91", "b1pw91"}, + {"XC_HYB_GGA_XC_MPW1PW", "mpw1pw"}, + {"XC_HYB_GGA_XC_MPW3LYP", "mpw3lyp"}, + {"XC_HYB_GGA_XC_SB98_1A", "sb98_1a"}, + {"XC_HYB_GGA_XC_SB98_1B", "sb98_1b"}, + {"XC_HYB_GGA_XC_SB98_1C", "sb98_1c"}, + {"XC_HYB_GGA_XC_SB98_2A", "sb98_2a"}, + {"XC_HYB_GGA_XC_SB98_2B", "sb98_2b"}, + {"XC_HYB_GGA_XC_SB98_2C", "sb98_2c"}, + {"XC_HYB_GGA_XC_HSE03", "hse03"}, + {"XC_HYB_GGA_XC_HSE06", "hse06"}, + {"XC_HYB_GGA_XC_HJS_PBE", "hjs_pbe"}, + {"XC_HYB_GGA_XC_HJS_PBE_SOL", "hjs_pbe_sol"}, + {"XC_HYB_GGA_XC_HJS_B88", "hjs_b88"}, + {"XC_HYB_GGA_XC_HJS_B97X", "hjs_b97x"}, + {"XC_HYB_GGA_XC_CAM_B3LYP", "cam_b3lyp"}, + {"XC_HYB_GGA_XC_TUNED_CAM_B3LYP", "tuned_cam_b3lyp"}, + {"XC_HYB_GGA_XC_BHANDH", "bhandh"}, + {"XC_HYB_GGA_XC_BHANDHLYP", "bhandhlyp"}, + {"XC_HYB_GGA_XC_MB3LYP_RC04", "mb3lyp_rc04"}, + {"XC_HYB_MGGA_XC_B88B95", "b88b95"}, + {"XC_HYB_MGGA_XC_B86B95", "b86b95"}, + {"XC_HYB_MGGA_XC_PW86B95", "pw86b95"}, + {"XC_HYB_MGGA_XC_BB1K", "bb1k"}, + {"XC_HYB_MGGA_XC_MPW1B95", "mpw1b95"}, + {"XC_HYB_MGGA_XC_MPWB1K", "mpwb1k"}, + {"XC_HYB_MGGA_XC_X1B95", "x1b95"}, + {"XC_HYB_MGGA_XC_XB1K", "xb1k"}, + {"XC_HYB_MGGA_XC_PW6B95", "pw6b95"}, + {"XC_HYB_MGGA_XC_PWB6K", "pwb6k"}, + {"XC_HYB_GGA_XC_MPWLYP1M", "mpwlyp1m"}, + {"XC_HYB_GGA_XC_REVB3LYP", "revb3lyp"}, + {"XC_HYB_GGA_XC_CAMY_BLYP", "camy_blyp"}, + {"XC_HYB_GGA_XC_PBE0_13", "pbe0_13"}, + {"XC_HYB_MGGA_XC_TPSSH", "tpssh"}, + {"XC_HYB_MGGA_XC_REVTPSSH", "revtpssh"}, + {"XC_HYB_GGA_XC_B3LYPS", "b3lyps"}, + {"XC_HYB_GGA_XC_QTP17", "qtp17"}, + {"XC_HYB_GGA_XC_B3LYP_MCM1", "b3lyp_mcm1"}, + {"XC_HYB_GGA_XC_B3LYP_MCM2", "b3lyp_mcm2"}, + {"XC_HYB_GGA_XC_WB97", "wb97"}, + {"XC_HYB_GGA_XC_WB97X", "wb97x"}, + {"XC_HYB_GGA_XC_LRC_WPBEH", "lrc_wpbeh"}, + {"XC_HYB_GGA_XC_WB97X_V", "wb97x_v"}, + {"XC_HYB_GGA_XC_LCY_PBE", "lcy_pbe"}, + {"XC_HYB_GGA_XC_LCY_BLYP", "lcy_blyp"}, + {"XC_HYB_GGA_XC_LC_VV10", "lc_vv10"}, + {"XC_HYB_GGA_XC_CAMY_B3LYP", "camy_b3lyp"}, + {"XC_HYB_GGA_XC_WB97X_D", "wb97x_d"}, // DFT-D2 + {"XC_HYB_GGA_XC_HPBEINT", "hpbeint"}, + {"XC_HYB_GGA_XC_LRC_WPBE", "lrc_wpbe"}, + {"XC_HYB_GGA_XC_B3LYP5", "b3lyp5"}, + {"XC_HYB_GGA_XC_EDF2", "edf2"}, + {"XC_HYB_GGA_XC_CAP0", "cap0"}, + {"XC_HYB_GGA_XC_LC_WPBE", "lc_wpbe"}, + {"XC_HYB_GGA_XC_HSE12", "hse12"}, + {"XC_HYB_GGA_XC_HSE12S", "hse12s"}, + {"XC_HYB_GGA_XC_HSE_SOL", "hse_sol"}, + {"XC_HYB_GGA_XC_CAM_QTP_01", "cam_qtp_01"}, + {"XC_HYB_GGA_XC_MPW1LYP", "mpw1lyp"}, + {"XC_HYB_GGA_XC_MPW1PBE", "mpw1pbe"}, + {"XC_HYB_GGA_XC_KMLYP", "kmlyp"}, + {"XC_HYB_GGA_XC_LC_WPBE_WHS", "lc_wpbe_whs"}, + {"XC_HYB_GGA_XC_LC_WPBEH_WHS", "lc_wpbeh_whs"}, + {"XC_HYB_GGA_XC_LC_WPBE08_WHS", "lc_wpbe08_whs"}, + {"XC_HYB_GGA_XC_LC_WPBESOL_WHS", "lc_wpbesol_whs"}, + {"XC_HYB_GGA_XC_CAM_QTP_00", "cam_qtp_00"}, + {"XC_HYB_GGA_XC_CAM_QTP_02", "cam_qtp_02"}, + {"XC_HYB_GGA_XC_LC_QTP", "lc_qtp"}, + {"XC_HYB_GGA_XC_BLYP35", "blyp35"}, + {"XC_HYB_MGGA_XC_WB97M_V", "wb97m_v"}, + {"XC_LDA_XC_1D_EHWLRG_1", "1d_ehwlrg_1"}, + {"XC_LDA_XC_1D_EHWLRG_2", "1d_ehwlrg_2"}, + {"XC_LDA_XC_1D_EHWLRG_3", "1d_ehwlrg_3"}, + {"XC_GGA_XC_HLE16", "hle16"}, + {"XC_LDA_XC_LP_A", "lp_a"}, + {"XC_LDA_XC_LP_B", "lp_b"}, + {"XC_HYB_MGGA_XC_B0KCIS", "b0kcis"}, + {"XC_MGGA_XC_LP90", "lp90"}, + {"XC_HYB_MGGA_XC_MPW1KCIS", "mpw1kcis"}, + {"XC_HYB_MGGA_XC_MPWKCIS1K", "mpwkcis1k"}, + {"XC_HYB_MGGA_XC_PBE1KCIS", "pbe1kcis"}, + {"XC_HYB_MGGA_XC_TPSS1KCIS", "tpss1kcis"}, + {"XC_HYB_GGA_XC_B5050LYP", "b5050lyp"}, + {"XC_LDA_XC_GDSMFB", "gdsmfb"}, + {"XC_GGA_XC_KT3", "kt3"}, + {"XC_HYB_LDA_XC_BN05", "bn05"}, + {"XC_HYB_GGA_XC_LB07", "lb07"}, + {"XC_HYB_MGGA_XC_B98", "b98"}, + {"XC_LDA_XC_TIH", "tih"}, + {"XC_HYB_GGA_XC_APBE0", "apbe0"}, + {"XC_HYB_GGA_XC_HAPBE", "hapbe"}, + {"XC_HYB_GGA_XC_RCAM_B3LYP", "rcam_b3lyp"}, + {"XC_HYB_GGA_XC_WC04", "wc04"}, + {"XC_HYB_GGA_XC_WP04", "wp04"}, + {"XC_HYB_GGA_XC_CAMH_B3LYP", "camh_b3lyp"}, + {"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_BLYPR", "lc_blypr"}, + {"XC_HYB_GGA_XC_MCAM_B3LYP", "mcam_b3lyp"}, + {"XC_MGGA_XC_VCML_RVV10", "vcml_rvv10"}, + {"XC_HYB_MGGA_XC_GAS22", "gas22"}, + {"XC_HYB_MGGA_XC_R2SCANH", "r2scanh"}, + {"XC_HYB_MGGA_XC_R2SCAN0", "r2scan0"}, + {"XC_HYB_MGGA_XC_R2SCAN50", "r2scan50"}, + {"XC_HYB_GGA_XC_CAM_PBEH", "cam_pbeh"}, + {"XC_HYB_GGA_XC_CAMY_PBEH", "camy_pbeh"}, + {"XC_HYB_MGGA_XC_EDMGGAH", "edmggah"}, + {"XC_HYB_MGGA_XC_LC_TMLYP", "lc_tmlyp"}, + }; + const std::map xcname_libxc_xplusc_ = { + {"XC_GGA_X_GAM+XC_GGA_C_GAM", "gam"}, + {"XC_GGA_X_HCTH_A+XC_GGA_C_HCTH_A", "hcth_a"}, + {"XC_HYB_MGGA_X_DLDF+XC_MGGA_C_DLDF", "dldf"}, + {"XC_GGA_X_Q2D+XC_GGA_C_Q2D", "q2d"}, + {"XC_GGA_X_PBE_MOL+XC_GGA_C_PBE_MOL", "pbe_mol"}, + {"XC_GGA_X_PBEINT+XC_GGA_C_PBEINT", "pbeint"}, + {"XC_HYB_GGA_X_N12_SX+XC_GGA_C_N12_SX", "n12_sx"}, + {"XC_GGA_X_N12+XC_GGA_C_N12", "n12"}, + {"XC_GGA_X_PBE+XC_GGA_C_PBE", "pbe"}, + {"XC_GGA_X_B88+XC_MGGA_C_B88", "b88"}, + {"XC_GGA_X_PW91+XC_GGA_C_PW91", "pw91"}, + {"XC_GGA_X_PBE_SOL+XC_GGA_C_PBE_SOL", "pbe_sol"}, + {"XC_GGA_X_AM05+XC_GGA_C_AM05", "am05"}, + {"XC_GGA_X_XPBE+XC_GGA_C_XPBE", "xpbe"}, + {"XC_GGA_X_RGE2+XC_GGA_C_RGE2", "rge2"}, + {"XC_GGA_X_SOGGA11+XC_GGA_C_SOGGA11", "sogga11"}, + {"XC_GGA_X_APBE+XC_GGA_C_APBE", "apbe"}, + {"XC_MGGA_X_TPSS+XC_MGGA_C_TPSS", "tpss"}, + {"XC_MGGA_X_M06_L+XC_MGGA_C_M06_L", "m06_l"}, + {"XC_HYB_MGGA_X_TAU_HCTH+XC_GGA_C_TAU_HCTH", "tau_hcth"}, + {"XC_MGGA_X_REVTPSS+XC_MGGA_C_REVTPSS", "revtpss"}, + {"XC_MGGA_X_PKZB+XC_MGGA_C_PKZB", "pkzb"}, + {"XC_MGGA_X_M11_L+XC_MGGA_C_M11_L", "m11_l"}, + {"XC_MGGA_X_MN12_L+XC_MGGA_C_MN12_L", "mn12_l"}, + {"XC_HYB_MGGA_X_MN12_SX+XC_MGGA_C_MN12_SX", "mn12_sx"}, + {"XC_MGGA_X_MN15_L+XC_MGGA_C_MN15_L", "mn15_l"}, + {"XC_MGGA_X_SCAN+XC_MGGA_C_SCAN", "scan"}, + {"XC_GGA_X_PBEFE+XC_GGA_C_PBEFE", "pbefe"}, + {"XC_HYB_MGGA_X_MN15+XC_MGGA_C_MN15", "mn15"}, + {"XC_HYB_MGGA_X_BMK+XC_GGA_C_BMK", "bmk"}, + {"XC_MGGA_X_REVM06_L+XC_MGGA_C_REVM06_L", "revm06_l"}, + {"XC_HYB_MGGA_X_M08_HX+XC_MGGA_C_M08_HX", "m08_hx"}, + {"XC_HYB_MGGA_X_M08_SO+XC_MGGA_C_M08_SO", "m08_so"}, + {"XC_HYB_MGGA_X_M11+XC_MGGA_C_M11", "m11"}, + {"XC_GGA_X_CHACHIYO+XC_GGA_C_CHACHIYO", "chachiyo"}, + {"XC_HYB_MGGA_X_REVM11+XC_MGGA_C_REVM11", "revm11"}, + {"XC_HYB_MGGA_X_REVM06+XC_MGGA_C_REVM06", "revm06"}, + {"XC_HYB_MGGA_X_M06_SX+XC_MGGA_C_M06_SX", "m06_sx"}, + {"XC_GGA_X_PBE_GAUSSIAN+XC_GGA_C_PBE_GAUSSIAN", "pbe_gaussian"}, + {"XC_HYB_GGA_X_SOGGA11_X+XC_GGA_C_SOGGA11_X", "sogga11_x"}, + {"XC_HYB_MGGA_X_M05+XC_MGGA_C_M05", "m05"}, + {"XC_HYB_MGGA_X_M05_2X+XC_MGGA_C_M05_2X", "m05_2x"}, + {"XC_HYB_MGGA_X_M06_HF+XC_MGGA_C_M06_HF", "m06_hf"}, + {"XC_HYB_MGGA_X_M06+XC_MGGA_C_M06", "m06"}, + {"XC_HYB_MGGA_X_M06_2X+XC_MGGA_C_M06_2X", "m06_2x"}, + {"XC_MGGA_X_RSCAN+XC_MGGA_C_RSCAN", "rscan"}, + {"XC_MGGA_X_R2SCAN+XC_MGGA_C_R2SCAN", "r2scan"}, + {"XC_GGA_X_SG4+XC_GGA_C_SG4", "sg4"}, + {"XC_MGGA_X_TM+XC_MGGA_C_TM", "tm"}, + {"XC_MGGA_X_REVSCAN+XC_MGGA_C_REVSCAN", "revscan"}, + {"XC_MGGA_X_REGTPSS+XC_GGA_C_REGTPSS", "regtpss"}, + {"XC_MGGA_X_R2SCAN01+XC_MGGA_C_R2SCAN01", "r2scan01"}, + {"XC_MGGA_X_RPPSCAN+XC_MGGA_C_RPPSCAN", "rppscan"}, + {"XC_MGGA_X_REVTM+XC_MGGA_C_REVTM", "revtm"}, + {"XC_MGGA_X_SCANL+XC_MGGA_C_SCANL", "scanl"}, + {"XC_MGGA_X_MGGAC+XC_GGA_C_MGGAC", "mggac"}, + {"XC_MGGA_X_R2SCANL+XC_MGGA_C_R2SCANL", "r2scanl"}, + {"XC_GGA_X_B88+XC_GGA_C_LYP", "blyp"}, + {"XC_GGA_X_B88+XC_GGA_C_P86", "bp86"}, + {"XC_GGA_X_PW91+XC_GGA_C_PW91", "pw91"}, + {"XC_GGA_X_PBE+XC_GGA_C_PBE", "pbe"}, + {"XC_GGA_X_PBE_SOL+XC_GGA_C_PBE_SOL", "pbesol"}, + {"XC_MGGA_X_PKZB+XC_MGGA_C_PKZB", "pkzb"}, + {"XC_MGGA_X_TPSS+XC_MGGA_C_TPSS", "tpss"}, + {"XC_MGGA_X_REVTPSS+XC_MGGA_C_REVTPSS", "revtpss"}, + {"XC_MGGA_X_SCAN+XC_MGGA_C_SCAN", "scan"}, + {"XC_GGA_X_SOGGA+XC_GGA_C_PBE", "sogga"}, + {"XC_MGGA_X_BLOC+XC_MGGA_C_TPSSLOC", "bloc"}, + {"XC_GGA_X_OPTX+XC_GGA_C_LYP", "olyp"}, + {"XC_GGA_X_RPBE+XC_GGA_C_PBE", "rpbe"}, + {"XC_GGA_X_B88+XC_GGA_C_PBE", "bpbe"}, + {"XC_GGA_X_MPW91+XC_GGA_C_PW91", "mpw91"}, + {"XC_MGGA_X_MS0+XC_GGA_C_REGTPSS", "ms0"}, + {"XC_MGGA_X_MS1+XC_GGA_C_REGTPSS", "ms1"}, + {"XC_MGGA_X_MS2+XC_GGA_C_REGTPSS", "ms2"}, + {"XC_HYB_MGGA_X_MS2H+XC_GGA_C_REGTPSS", "ms2h"}, + {"XC_MGGA_X_MVS+XC_GGA_C_REGTPSS", "mvs"}, + {"XC_HYB_MGGA_X_MVSH+XC_GGA_C_REGTPSS", "mvsh"}, + {"XC_GGA_X_SOGGA11+XC_GGA_C_SOGGA11", "sogga11"}, + {"XC_HYB_GGA_X_SOGGA11_X+XC_GGA_C_SOGGA11_X", "sogga11-x"}, + {"XC_HYB_MGGA_X_DLDF+XC_MGGA_C_DLDF", "dldf"}, + {"XC_GGA_X_GAM+XC_GGA_C_GAM", "gam"}, + {"XC_MGGA_X_M06_L+XC_MGGA_C_M06_L", "m06-l"}, + {"XC_MGGA_X_M11_L+XC_MGGA_C_M11_L", "m11-l"}, + {"XC_MGGA_X_MN12_L+XC_MGGA_C_MN12_L", "mn12-l"}, + {"XC_MGGA_X_MN15_L+XC_MGGA_C_MN15_L", "mn15-l"}, + {"XC_GGA_X_N12+XC_GGA_C_N12", "n12"}, + {"XC_HYB_GGA_X_N12_SX+XC_GGA_C_N12_SX", "n12-sx"}, + {"XC_HYB_MGGA_X_MN12_SX+XC_MGGA_C_MN12_SX", "mn12-sx"}, + {"XC_HYB_MGGA_X_MN15+XC_MGGA_C_MN15", "mn15"}, + {"XC_MGGA_X_MBEEF+XC_GGA_C_PBE_SOL", "mbeef"}, + {"XC_HYB_MGGA_X_SCAN0+XC_MGGA_C_SCAN", "scan0"}, + {"XC_GGA_X_PBE+XC_GGA_C_OP_PBE", "pbeop"}, + {"XC_GGA_X_B88+XC_GGA_C_OP_B88", "bop"} + }; + + void _xcname_libxc_xplusc(const std::string& xcpattern, std::string& xname) + { + std::vector xc_words = FmtCore::split(xcpattern, "+"); + std::for_each(xc_words.begin(), xc_words.end(), [](std::string& s) { + s = (FmtCore::startswith(s, "XC_")? s: "XC_" + s); }); // add XC_ if not present + assert(xc_words.size() == 2); + + std::vector words = FmtCore::split(xc_words[0], "_"); + const std::string key = (words[2] == "X")? + xc_words[0] + "+" + xc_words[1]: xc_words[1] + "+" + xc_words[0]; + + if (xcname_libxc_xplusc_.find(key) != xcname_libxc_xplusc_.end()) { + xname = xcname_libxc_xplusc_.at(key); + } else { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::xcname_libxc_xplusc", + "XC's LibXC-notation on `" + xcpattern + "` not recognized"); + } + } + + void _xcname_libxc_xc(const std::string& xcpattern, std::string& xname) + { + // add XC_ if not present + const std::string key = FmtCore::startswith(xcpattern, "XC_")? xcpattern: "XC_" + xcpattern; + + if (xcname_libxc_xc_.find(key) != xcname_libxc_xc_.end()) { + xname = xcname_libxc_xc_.at(key); + } else { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::xcname_libxc_xc", + "XC's LibXC-notation on `" + xcpattern + "` not recognized"); + } + } + + void _xcname_libxc(const std::string& xcpattern, std::string& xname) + { + if (xcpattern.find("+") != std::string::npos) { + _xcname_libxc_xplusc(xcpattern, xname); + } else { + _xcname_libxc_xc(xcpattern, xname); + } + } + + std::string _xcname(const std::string& xcpattern) + { + std::string xcname = xcpattern; + const std::regex pattern("(LDA|GGA|MGGA|HYB|HYB_LDA|HYB_GGA|HYB_MGGA)_(X|C|XC|K)_(.*)"); + // as long as there is piece in xcpattern that can match, we can search for the corresponding name + if (std::regex_search(xcpattern, pattern)) { + _xcname_libxc(xcpattern, xcname); + } + return xcname; + } +} + +/** +import os +import re +def read_xc_func_h(fn): + with open(fn) as f: + lines = f.readlines() + out = {} + for line in lines: + words = line.strip().split() + xc, xcid = words[1], int(words[2]) + xc_annos = ' '.join(words[4:-1]) + out[xc] = {'id': xcid, 'annos': xc_annos} + return out + +def sort_xc(xc_data): + '''Sort the xc functionals into x, c, xc, k functionals. + + Parameters + ---------- + xc_data : dict + from function read_xc_func_h + + Returns + ------- + dict, dict, dict, dict + The dictionaries of x, c, xc, k functionals, whose keys are the + like LDA, GGA, MGGA, HYB, HYB_LDA, HYB_GGA, HYB_MGGA, values are + the dictionaries of the functionals, whose keys are the conventional + xc name, values include approx, annos, id, full. + ''' + x, c, xc, k = {}, {}, {}, {} + dictmap = {'X': x, 'C': c, 'XC': xc, 'K': k} + xcpat = r'XC_(LDA|GGA|MGGA|HYB|HYB_LDA|HYB_GGA|HYB_MGGA)_(X|C|XC|K)_(.*)' + for xc_name, data in xc_data.items(): + m = re.match(xcpat, xc_name) + if m is None: + print('Warning: cannot match', xc_name) + continue + approx, type_, name = m.groups() + dictmap[type_][name] = {'approx': approx, 'annos': data['annos'], + 'id': data['id'], 'full': xc_name} + return x, c, xc, k + +def pair_xc(x, c): + ''' + Pair the x and c functionals. + + Parameters + ---------- + x : dict + The dictionary of x functionals, whose keys are the conventional + xc name, values include approx, annos, id, full. + + c : dict + the same as x + + Returns + ------- + dict, dict + The dictionary of paired and unpaired x and c functionals, whose keys are the + conventional xc name, values are the dictionary of x and c functionals. + ''' + paired, unpaired = {}, {} + for xc_name, data in x.items(): + if xc_name in c: + paired[xc_name] = {'x': data, 'c': c[xc_name]} + else: + unpaired[xc_name] = data + return paired, unpaired + +def xc_to_stdmap(xc, conventional_lower=True): + '''print the xc in the way of c++ std::map. + + Parameters + ---------- + xc : dict + The dictionary of xc functionals, whose keys are the conventional + xc name, values include approx, annos, id, full. + conventional_lower : bool + Whether to convert the conventional name to lower case. + + Returns + ------- + str + The string of c++ code, std::map mapping + the full name of xc to its conventional name. + ''' + out = 'const std::map xcname_libxc_xc_ = {\n' + for name, data in xc.items(): + name = name.lower() if conventional_lower else name + out += ' {"%s", "%s"},\n' % (data['full'], name) + out += '};\n' + return out + +def paired_xc_to_stdmap(pairs, conventional_lower=True): + '''print the xc in the way of c++ std::map. + + Parameters + ---------- + pairs : dict + The dictionary of xc functionals, whose keys are the conventional + xc name, values include approx, annos, id, full. + conventional_lower : bool + Whether to convert the conventional name to lower case. + + Returns + ------- + str + The string of c++ code, std::map mapping + the full name of xc to its conventional name. + ''' + out = 'const std::map xcname_libxc_xplusc_ = {\n' + for name, data in pairs.items(): + name = name.lower() if conventional_lower else name + plus = f'{data["x"]["full"]}+{data["c"]["full"]}' + out += ' {"%s", "%s"},\n' % (plus, name) + # sulp = f'{data["c"]["full"]}+{data["x"]["full"]}' + # out += ' {"%s", "%s"},\n' % (sulp, name) + out += '};\n' + return out + +def special_x_and_c(x, c): + '''Special pairings of x and c functionals. The following data sheet is + from Pyscf: + https://github.com/pyscf/pyscf/blob/master/pyscf/dft/xcfun.py + Thanks for pointing out the bug by @QuantumMiska and the help from wsr (@hebrewsnabla) + + + Parameters + ---------- + x : dict + The dictionary of x functionals, whose keys are the conventional + xc name, values include approx, annos, id, full. + + c : dict + the same as x + + Returns + ------- + dict + The dictionary of special pairings of x and c functionals. + ''' + DATA = { + 'BLYP' : 'B88,LYP', + 'BP86' : 'B88,P86', + 'PW91' : 'PW91,PW91', + 'PBE' : 'PBE,PBE', + 'REVPBE' : 'REVPBE,PBE', + 'PBESOL' : 'PBE_SOL,PBE_SOL', + 'PKZB' : 'PKZB,PKZB', + 'TPSS' : 'TPSS,TPSS', + 'REVTPSS' : 'REVTPSS,REVTPSS', + 'SCAN' : 'SCAN,SCAN', + 'SOGGA' : 'SOGGA,PBE', + 'BLOC' : 'BLOC,TPSSLOC', + 'OLYP' : 'OPTX,LYP', + 'RPBE' : 'RPBE,PBE', + 'BPBE' : 'B88,PBE', + 'MPW91' : 'MPW91,PW91', + 'HFLYP' : 'HF,LYP', + 'HFPW92' : 'HF,PWMOD', + 'SPW92' : 'SLATER,PWMOD', + 'SVWN' : 'SLATER,VWN', + 'MS0' : 'MS0,REGTPSS', + 'MS1' : 'MS1,REGTPSS', + 'MS2' : 'MS2,REGTPSS', + 'MS2H' : 'MS2H,REGTPSS', + 'MVS' : 'MVS,REGTPSS', + 'MVSH' : 'MVSH,REGTPSS', + 'SOGGA11' : 'SOGGA11,SOGGA11', + 'SOGGA11-X': 'SOGGA11_X,SOGGA11_X', + 'KT1' : 'KT1X,VWN', + 'DLDF' : 'DLDF,DLDF', + 'GAM' : 'GAM,GAM', + 'M06-L' : 'M06_L,M06_L', + 'M11-L' : 'M11_L,M11_L', + 'MN12-L' : 'MN12_L,MN12_L', + 'MN15-L' : 'MN15_L,MN15_L', + 'N12' : 'N12,N12', + 'N12-SX' : 'N12_SX,N12_SX', + 'MN12-SX' : 'MN12_SX,MN12_SX', + 'MN15' : 'MN15,MN15', + 'MBEEF' : 'MBEEF,PBE_SOL', + 'SCAN0' : 'SCAN0,SCAN', + 'PBEOP' : 'PBE,OP_PBE', + 'BOP' : 'B88,OP_B88', + } + paired = {} + for name, data in DATA.items(): + xname, cname = data.split(',') + if xname in x and cname in c: + paired[name] = {'x': x[xname], 'c': c[cname]} + else: + print(f'Warning: {name} not found in x or c: {xname}, {cname}') + return paired + +def print_xc(xc): + print(f'{"Name":20s} {"Full":30s} {"Appr":10s} {"Annos"}') + for name, data in xc.items(): + print(f'{name:20s} {data["full"]:30s} {data["approx"]:10s} {data["annos"]}') + +if __name__ == '__main__': + libxc = '/root/soft/libxc/libxc-6.2.2' + f = 'src/xc_funcs.h' + xc_data = read_xc_func_h(os.path.join(libxc, f)) + x, c, xc, k = sort_xc(xc_data) + pairs, others = pair_xc(x, c) + special = special_x_and_c(x, c) + # print(xc_to_stdmap(xc)) + # print(paired_xc_to_stdmap(pairs)) + # print_xc(others) + print(paired_xc_to_stdmap(special)) + */ +#endif // DFTD3_XCNAME_H_ diff --git a/source/module_hamilt_general/module_vdw/dftd3_xc_param.h b/source/module_hamilt_general/module_vdw/dftd3_xc_param.h new file mode 100644 index 0000000000..0090c6f308 --- /dev/null +++ b/source/module_hamilt_general/module_vdw/dftd3_xc_param.h @@ -0,0 +1,595 @@ +#ifndef DFTD3_XC_PARAM_H_ +#define DFTD3_XC_PARAM_H_ +/** + * Intro + * ----- + * This file stores XC dependent DFT-D3 parameters for Grimme-D3 + * dispersion correction. + * + * Supported forms: + * + * DFT-D3(0): zero-damping + * DFT-D3(BJ): Becke-Johnson damping + * DFT-D3M(0): zero-damping with modified damping function + * DFT-D3M(BJ): Becke-Johnson damping with modified damping function + * + * A detailed introduction of undamped, and BJ damping, the modified + * damping can be found in DFT-D3 software manual, see: + * https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/man.pdf + * + * Other excellent learning materials (where you can find expression + * of both DFT-D2 and DFT-D3): + * DFT-D2: https://www.vasp.at/wiki/index.php/DFT-D2 + * DFT-D3: https://www.vasp.at/wiki/index.php/DFT-D3 + * + * Usage + * ----- + * call function DFTD3::search(xc, method, param) to get the DFT-D3 parameters + * for the given XC functional. The obtained param should be a std::vector, + * in which the first 9 elements are the DFT-D3 parameters: + * 's6', 'sr6', 'a1', 's8', 'sr8', 'a2', 's9', 'alp', 'bet' + * + * ParamNotFoundError + * ------------------ + * If the requested D3 parameters of XC are not found, then the ABACUS will + * WARNING_QUIT with the message "DFT-D3 parameters for XC not found". + * + * Other dispersion correction + * --------------------------- + * there are other kinds of dispersion correction, such as the xc VV09, VV10, + * and rVV10, and the vdw-DF family nonlocal dispersion correction. They will + * be mixed directly with the correlation and exchange part, which act + * differently from the DFT-D2 and D3 methods. + * + * Special: Omega-B97 family + * ------------------------- + * (thanks for help and discussion with @hhebrewsnabla and @moolawooda) + * wB97 XC family is special, their DFT-D3 supports are quite complicated. + * + * wB97 long-range exx with B97 + * wB97X wB97 with additional short-range exx + * wB97X-D wB97X_D from libXC with DFTD2, not in DFTD3 framework + * wB97X-D3 wB97X_D3 from libXC with DFTD3(0) + * wB97X-D3(BJ) wB97X_V from libXC with DFTD3(BJ) + * wB97X-V with VV10, not in DFTD3 framework + * wB97M-V with VV10, not in DFTD3 framework + * + * Recommended: http://bbs.keinsci.com/thread-19076-1-1.html + * Related information from Pyscf Github repo: + * https://github.com/pyscf/pyscf/issues/2069 + * + */ +#include +#include +#include +#include +#include +#include +#include "dftd3_xc_name.h" +// 's6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet' +namespace DFTD3 { + // DFT-D3(BJ) + std::map> bj = { + {"bp", {1.0, 0.3946, 0.3946, 3.2822, 4.8516, 4.8516, 1.0, 14.0, 0.0}}, + {"blyp", {1.0, 0.4298, 0.4298, 2.6996, 4.2359, 4.2359, 1.0, 14.0, 0.0}}, + {"revpbe", {1.0, 0.5238, 0.5238, 2.355, 3.5016, 3.5016, 1.0, 14.0, 0.0}}, + {"rpbe", {1.0, 0.182, 0.182, 0.8318, 4.0094, 4.0094, 1.0, 14.0, 0.0}}, + {"b97_d", {1.0, 0.5545, 0.5545, 2.2609, 3.2297, 3.2297, 1.0, 14.0, 0.0}}, + {"b973c", {1.0, 0.37, 0.37, 1.5, 4.1, 4.1, 1.0, 14.0, 0.0}}, + {"pbe", {1.0, 0.4289, 0.4289, 0.7875, 4.4407, 4.4407, 1.0, 14.0, 0.0}}, + {"rpw86pbe", {1.0, 0.4613, 0.4613, 1.3845, 4.5062, 4.5062, 1.0, 14.0, 0.0}}, + {"b3lyp", {1.0, 0.3981, 0.3981, 1.9889, 4.4211, 4.4211, 1.0, 14.0, 0.0}}, + {"tpss", {1.0, 0.4535, 0.4535, 1.9435, 4.4752, 4.4752, 1.0, 14.0, 0.0}}, + {"hf", {1.0, 0.3385, 0.3385, 0.9171, 2.883, 2.883, 1.0, 14.0, 0.0}}, + {"tpss0", {1.0, 0.3768, 0.3768, 1.2576, 4.5865, 4.5865, 1.0, 14.0, 0.0}}, + {"pbe0", {1.0, 0.4145, 0.4145, 1.2177, 4.8593, 4.8593, 1.0, 14.0, 0.0}}, + {"hse06", {1.0, 0.383, 0.383, 2.31, 5.685, 5.685, 1.0, 14.0, 0.0}}, + {"hse", {1.0, 0.383, 0.383, 2.31, 5.685, 5.685, 1.0, 14.0, 0.0}}, // ABACUS implements HSE06 as HSE + {"revpbe38", {1.0, 0.4309, 0.4309, 1.476, 3.9446, 3.9446, 1.0, 14.0, 0.0}}, + {"pw6b95", {1.0, 0.2076, 0.2076, 0.7257, 6.375, 6.375, 1.0, 14.0, 0.0}}, + {"b2plyp", {0.64, 0.3065, 0.3065, 0.9147, 5.057, 5.057, 1.0, 14.0, 0.0}}, + {"dsdblyp", {0.5, 0.0, 0.0, 0.213, 6.0519, 6.0519, 1.0, 14.0, 0.0}}, + {"dsdblypfc", {0.5, 0.0009, 0.0009, 0.2112, 5.9807, 5.9807, 1.0, 14.0, 0.0}}, + {"dodscan66", {0.3152, 0.0, 0.0, 0.0, 5.75, 5.75, 1.0, 14.0, 0.0}}, + {"revdsdblyp", {0.5451, 0.0, 0.0, 0.0, 5.2, 5.2, 1.0, 14.0, 0.0}}, + {"revdsdpbep86", {0.4377, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"revdsdpbeb95", {0.3686, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"revdsdpbe", {0.5746, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"revdodblyp", {0.6145, 0.0, 0.0, 0.0, 5.2, 5.2, 1.0, 14.0, 0.0}}, + {"revdodpbep86", {0.477, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"revdodpbeb95", {0.4107, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"revdodpbe", {0.6067, 0.0, 0.0, 0.0, 5.5, 5.5, 1.0, 14.0, 0.0}}, + {"bop", {1.0, 0.487, 0.487, 3.295, 3.5043, 3.5043, 1.0, 14.0, 0.0}}, + {"mpwlyp", {1.0, 0.4831, 0.4831, 2.0077, 4.5323, 4.5323, 1.0, 14.0, 0.0}}, + {"olyp", {1.0, 0.5299, 0.5299, 2.6205, 2.8065, 2.8065, 1.0, 14.0, 0.0}}, + {"pbesol", {1.0, 0.4466, 0.4466, 2.9491, 6.1742, 6.1742, 1.0, 14.0, 0.0}}, + {"bpbe", {1.0, 0.4567, 0.4567, 4.0728, 4.3908, 4.3908, 1.0, 14.0, 0.0}}, + {"opbe", {1.0, 0.5512, 0.5512, 3.3816, 2.9444, 2.9444, 1.0, 14.0, 0.0}}, + {"ssb", {1.0, -0.0952, -0.0952, -0.1744, 5.217, 5.217, 1.0, 14.0, 0.0}}, + {"revssb", {1.0, 0.472, 0.472, 0.4389, 4.0986, 4.0986, 1.0, 14.0, 0.0}}, + {"otpss", {1.0, 0.4634, 0.4634, 2.7495, 4.3153, 4.3153, 1.0, 14.0, 0.0}}, + {"b3pw91", {1.0, 0.4312, 0.4312, 2.8524, 4.4693, 4.4693, 1.0, 14.0, 0.0}}, + {"bhlyp", {1.0, 0.2793, 0.2793, 1.0354, 4.9615, 4.9615, 1.0, 14.0, 0.0}}, + {"revpbe0", {1.0, 0.4679, 0.4679, 1.7588, 3.7619, 3.7619, 1.0, 14.0, 0.0}}, + {"tpssh", {1.0, 0.4529, 0.4529, 2.2382, 4.655, 4.655, 1.0, 14.0, 0.0}}, + {"mpw1b95", {1.0, 0.1955, 0.1955, 1.0508, 6.4177, 6.4177, 1.0, 14.0, 0.0}}, + {"pwb6k", {1.0, 0.1805, 0.1805, 0.9383, 7.7627, 7.7627, 1.0, 14.0, 0.0}}, + {"b1b95", {1.0, 0.2092, 0.2092, 1.4507, 5.5545, 5.5545, 1.0, 14.0, 0.0}}, + {"bmk", {1.0, 0.194, 0.194, 2.086, 5.9197, 5.9197, 1.0, 14.0, 0.0}}, + {"camb3lyp", {1.0, 0.3708, 0.3708, 2.0674, 5.4743, 5.4743, 1.0, 14.0, 0.0}}, + {"lcwpbe", {1.0, 0.3919, 0.3919, 1.8541, 5.0897, 5.0897, 1.0, 14.0, 0.0}}, + {"b2gpplyp", {0.56, 0.0, 0.0, 0.2597, 6.3332, 6.3332, 1.0, 14.0, 0.0}}, + {"ptpss", {0.75, 0.0, 0.0, 0.2804, 6.5745, 6.5745, 1.0, 14.0, 0.0}}, + {"pwpb95", {0.82, 0.0, 0.0, 0.2904, 7.3141, 7.3141, 1.0, 14.0, 0.0}}, + {"hf_mixed", {1.0, 0.5607, 0.5607, 3.9027, 4.5622, 4.5622, 1.0, 14.0, 0.0}}, + {"hf_sv", {1.0, 0.4249, 0.4249, 2.1849, 4.2783, 4.2783, 1.0, 14.0, 0.0}}, + {"hf_minis", {1.0, 0.1702, 0.1702, 0.9841, 3.8506, 3.8506, 1.0, 14.0, 0.0}}, + {"b3lyp_631gd", {1.0, 0.5014, 0.5014, 4.0672, 4.8409, 4.8409, 1.0, 14.0, 0.0}}, + {"hcth120", {1.0, 0.3563, 0.3563, 1.0821, 4.3359, 4.3359, 1.0, 14.0, 0.0}}, + {"dftb3", {1.0, 0.5719, 0.5719, 0.5883, 3.6017, 3.6017, 1.0, 14.0, 0.0}}, + {"pw1pw", {1.0, 0.3807, 0.3807, 2.3363, 5.8844, 5.8844, 1.0, 14.0, 0.0}}, + {"pwgga", {1.0, 0.2211, 0.2211, 2.691, 6.7278, 6.7278, 1.0, 14.0, 0.0}}, + {"hsesol", {1.0, 0.465, 0.465, 2.9215, 6.2003, 6.2003, 1.0, 14.0, 0.0}}, + {"hf3c", {1.0, 0.4171, 0.4171, 0.8777, 2.9149, 2.9149, 1.0, 14.0, 0.0}}, + {"hf3cv", {1.0, 0.3063, 0.3063, 0.5022, 3.9856, 3.9856, 1.0, 14.0, 0.0}}, + {"pbeh3c", {1.0, 0.486, 0.486, 0.0, 4.5, 4.5, 1.0, 14.0, 0.0}}, + {"scan", {1.0, 0.538, 0.538, 0.0, 5.42, 5.42, 1.0, 14.0, 0.0}}, + {"rscan", {1.0, 0.47023427, 0.47023427, 1.08859014, 5.73408312, 5.73408312, 1.0, 14.0, 0.0}}, + {"r2scan", {1.0, 0.49484001, 0.49484001, 0.78981345, 5.73083694, 5.73083694, 1.0, 14.0, 0.0}}, + {"r2scanh", {1.0, 0.4709, 0.4709, 1.1236, 5.9157, 5.9157, 1.0, 14.0, 0.0}}, + {"r2scan0", {1.0, 0.4534, 0.4534, 1.1846, 5.8972, 5.8972, 1.0, 14.0, 0.0}}, + {"r2scan50", {1.0, 0.4311, 0.4311, 1.3294, 5.924, 5.924, 1.0, 14.0, 0.0}}, + {"wb97x_v", {1.0, 0.0, 0.0, 0.2641, 5.4959, 5.4959, 1.0, 14.0, 0.0}}, + // NOTE: the key `wb97x_v` directly corresonding to HYB_GGA_XC_WB97X_V, which can be further + // employed to construct either wB97X-V with VV10, or wB97X-D3BJ with D3BJ. Here it is the D3BJ + // parameter of wB97X-D3BJ, instead of those of wB97X-V. + {"wb97m", {1.0, 0.566, 0.566, 0.3908, 3.128, 3.128, 1.0, 14.0, 0.0}}, + {"b97m", {1.0, -0.078, -0.078, 0.1384, 5.5946, 5.5946, 1.0, 14.0, 0.0}}, + {"pbehpbe", {1.0, 0.0, 0.0, 1.1152, 6.7184, 6.7184, 1.0, 14.0, 0.0}}, + {"xlyp", {1.0, 0.0809, 0.0809, 1.5669, 5.3166, 5.3166, 1.0, 14.0, 0.0}}, + {"mpwpw", {1.0, 0.3168, 0.3168, 1.7974, 4.7732, 4.7732, 1.0, 14.0, 0.0}}, + {"hcth407", {1.0, 0.0, 0.0, 0.649, 4.8162, 4.8162, 1.0, 14.0, 0.0}}, + {"revtpss", {1.0, 0.4326, 0.4326, 1.4023, 4.4723, 4.4723, 1.0, 14.0, 0.0}}, + {"tauhcth", {1.0, 0.0, 0.0, 1.2626, 5.6162, 5.6162, 1.0, 14.0, 0.0}}, + {"b3p", {1.0, 0.4601, 0.4601, 3.3211, 4.9858, 4.9858, 1.0, 14.0, 0.0}}, + {"b1p", {1.0, 0.4724, 0.4724, 3.5681, 4.9858, 4.9858, 1.0, 14.0, 0.0}}, + {"b1lyp", {1.0, 0.1986, 0.1986, 2.1167, 5.3875, 5.3875, 1.0, 14.0, 0.0}}, + {"mpwb1k", {1.0, 0.1474, 0.1474, 0.9499, 6.6223, 6.6223, 1.0, 14.0, 0.0}}, + {"mpw1pw", {1.0, 0.3342, 0.3342, 1.8744, 4.9819, 4.9819, 1.0, 14.0, 0.0}}, + {"mpw1kcis", {1.0, 0.0576, 0.0576, 1.0893, 5.5314, 5.5314, 1.0, 14.0, 0.0}}, + {"pbeh1pbe", {1.0, 0.0, 0.0, 1.4877, 7.0385, 7.0385, 1.0, 14.0, 0.0}}, + {"pbe1kcis", {1.0, 0.0, 0.0, 0.7688, 6.2794, 6.2794, 1.0, 14.0, 0.0}}, + {"x3lyp", {1.0, 0.2022, 0.2022, 1.5744, 5.4184, 5.4184, 1.0, 14.0, 0.0}}, + {"o3lyp", {1.0, 0.0963, 0.0963, 1.8171, 5.994, 5.994, 1.0, 14.0, 0.0}}, + {"b97_1", {1.0, 0.0, 0.0, 0.4814, 6.2279, 6.2279, 1.0, 14.0, 0.0}}, + {"b97_2", {1.0, 0.0, 0.0, 0.9448, 5.994, 5.994, 1.0, 14.0, 0.0}}, + {"b98", {1.0, 0.0, 0.0, 0.7086, 6.0672, 6.0672, 1.0, 14.0, 0.0}}, + {"hiss", {1.0, 0.0, 0.0, 1.6112, 7.3539, 7.3539, 1.0, 14.0, 0.0}}, + {"hse03", {1.0, 0.0, 0.0, 1.1243, 6.8889, 6.8889, 1.0, 14.0, 0.0}}, + {"revtpssh", {1.0, 0.266, 0.266, 1.4076, 5.3761, 5.3761, 1.0, 14.0, 0.0}}, + {"revtpss0", {1.0, 0.2218, 0.2218, 1.6151, 5.7985, 5.7985, 1.0, 14.0, 0.0}}, + {"tpss1kcis", {1.0, 0.0, 0.0, 1.0542, 6.0201, 6.0201, 1.0, 14.0, 0.0}}, + {"tauhcthhyb", {1.0, 0.0, 0.0, 0.9585, 10.1389, 10.1389, 1.0, 14.0, 0.0}}, + {"m11", {1.0, 0.0, 0.0, 2.8112, 10.1389, 10.1389, 1.0, 14.0, 0.0}}, + {"sogga11x", {1.0, 0.133, 0.133, 1.1426, 5.7381, 5.7381, 1.0, 14.0, 0.0}}, + {"n12sx", {1.0, 0.3283, 0.3283, 2.49, 5.7898, 5.7898, 1.0, 14.0, 0.0}}, + {"mn12sx", {1.0, 0.0983, 0.0983, 1.1674, 8.0259, 8.0259, 1.0, 14.0, 0.0}}, + {"mn12l", {1.0, 0.0, 0.0, 2.2674, 9.1494, 9.1494, 1.0, 14.0, 0.0}}, + {"mn15", {1.0, 2.0971, 2.0971, 0.7862, 7.5923, 7.5923, 1.0, 14.0, 0.0}}, + {"lc_whpbe", {1.0, 0.2746, 0.2746, 1.1908, 5.3157, 5.3157, 1.0, 14.0, 0.0}}, + {"mpw2plyp", {0.66, 0.4105, 0.4105, 0.6223, 5.0136, 5.0136, 1.0, 14.0, 0.0}}, + {"pw91", {1.0, 0.6319, 0.6319, 1.9598, 4.5718, 4.5718, 1.0, 14.0, 0.0}}, + {"drpa75", {0.3754, 0.0, 0.0, 0.0, 4.5048, 4.5048, 1.0, 14.0, 0.0}}, + {"scsdrpa75", {0.2528, 0.0, 0.0, 0.0, 4.505, 4.505, 1.0, 14.0, 0.0}}, + {"optscsdrpa75", {0.2546, 0.0, 0.0, 0.0, 4.505, 4.505, 1.0, 14.0, 0.0}}, + {"dsdpbedrpa75", {0.3223, 0.0, 0.0, 0.0, 4.505, 4.505, 1.0, 14.0, 0.0}}, + {"dsdpbep86drpa75", {0.3012, 0.0, 0.0, 0.0, 4.505, 4.505, 1.0, 14.0, 0.0}}, + {"dsdpbep86_2011", {0.418, 0.0, 0.0, 0.0, 5.65, 5.65, 1.0, 14.0, 0.0}}, + {"dsdsvwn5", {0.46, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dsdsp86", {0.3, 0.0, 0.0, 0.0, 5.8, 5.8, 1.0, 14.0, 0.0}}, + {"dsdslyp", {0.3, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dsdspbe", {0.4, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + {"dsdbvwn5", {0.61, 0.0, 0.0, 0.0, 5.2, 5.2, 1.0, 14.0, 0.0}}, + {"dsdblyp_2013", {0.57, 0.0, 0.0, 0.0, 5.4, 5.4, 1.0, 14.0, 0.0}}, + {"dsdbpbe", {1.22, 0.0, 0.0, 0.0, 6.6, 6.6, 1.0, 14.0, 0.0}}, + {"dsdbp86", {0.76, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + {"dsdbpw91", {1.14, 0.0, 0.0, 0.0, 6.5, 6.5, 1.0, 14.0, 0.0}}, + {"dsdbb95", {1.02, 0.0, 0.0, 0.0, 6.8, 6.8, 1.0, 14.0, 0.0}}, + {"dsdpbevwn5", {0.54, 0.0, 0.0, 0.0, 5.1, 5.1, 1.0, 14.0, 0.0}}, + {"dsdpbelyp", {0.43, 0.0, 0.0, 0.0, 5.2, 5.2, 1.0, 14.0, 0.0}}, + {"dsdpbe", {0.78, 0.0, 0.0, 0.0, 6.1, 6.1, 1.0, 14.0, 0.0}}, + {"dsdpbep86", {0.48, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dsdpbepw91", {0.73, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + {"dsdpbeb95", {0.61, 0.0, 0.0, 0.0, 6.2, 6.2, 1.0, 14.0, 0.0}}, + {"dsdpbehb95", {0.58, 0.0, 0.0, 0.0, 6.2, 6.2, 1.0, 14.0, 0.0}}, + {"dsdpbehp86", {0.46, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dsdmpwlyp", {0.48, 0.0, 0.0, 0.0, 5.3, 5.3, 1.0, 14.0, 0.0}}, + {"dsdmpwpw91", {0.9, 0.0, 0.0, 0.0, 6.2, 6.2, 1.0, 14.0, 0.0}}, + {"dsdmpwp86", {0.59, 0.0, 0.0, 0.0, 5.8, 5.8, 1.0, 14.0, 0.0}}, + {"dsdmpwpbe", {0.96, 0.0, 0.0, 0.0, 6.3, 6.3, 1.0, 14.0, 0.0}}, + {"dsdmpwb95", {0.82, 0.0, 0.0, 0.0, 6.6, 6.6, 1.0, 14.0, 0.0}}, + {"dsdhsepbe", {0.79, 0.0, 0.0, 0.0, 6.1, 6.1, 1.0, 14.0, 0.0}}, + {"dsdhsepw91", {0.74, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + {"dsdhsep86", {0.46, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dsdhselyp", {0.4, 0.0, 0.0, 0.0, 5.2, 5.2, 1.0, 14.0, 0.0}}, + {"dsdtpss", {0.72, 0.0, 0.0, 0.0, 6.5, 6.5, 1.0, 14.0, 0.0}}, + {"dsdtpssb95", {0.91, 0.0, 0.0, 0.0, 7.9, 7.9, 1.0, 14.0, 0.0}}, + {"dsdolyp", {0.93, 0.0, 0.0, 0.0, 5.8, 5.8, 1.0, 14.0, 0.0}}, + {"dsdxlyp", {0.51, 0.0, 0.0, 0.0, 5.3, 5.3, 1.0, 14.0, 0.0}}, + {"dsdxb95", {0.92, 0.0, 0.0, 0.0, 6.7, 6.7, 1.0, 14.0, 0.0}}, + {"dsdb98", {0.07, 0.0, 0.0, 0.0, 3.7, 3.7, 1.0, 14.0, 0.0}}, + {"dsdbmk", {0.17, 0.0, 0.0, 0.0, 3.9, 3.9, 1.0, 14.0, 0.0}}, + {"dsdthcth", {0.39, 0.0, 0.0, 0.0, 4.8, 4.8, 1.0, 14.0, 0.0}}, + {"dsdhcth407", {0.53, 0.0, 0.0, 0.0, 5.0, 5.0, 1.0, 14.0, 0.0}}, + {"dodsvwn5", {0.57, 0.0, 0.0, 0.0, 5.6, 5.6, 1.0, 14.0, 0.0}}, + {"dodblyp", {0.96, 0.0, 0.0, 0.0, 5.1, 5.1, 1.0, 14.0, 0.0}}, + {"dodpbe", {0.91, 0.0, 0.0, 0.0, 5.9, 5.9, 1.0, 14.0, 0.0}}, + {"dodpbep86", {0.72, 0.0, 0.0, 0.0, 5.4, 5.4, 1.0, 14.0, 0.0}}, + {"dodpbeb95", {0.71, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + {"dodhsep86", {0.69, 0.0, 0.0, 0.0, 5.4, 5.4, 1.0, 14.0, 0.0}}, + {"dodpbehb95", {0.67, 0.0, 0.0, 0.0, 6.0, 6.0, 1.0, 14.0, 0.0}}, + }; + // DFT-D3(0) + std::map> zero = { + {"slaterdirac", {1.0, 0.999, 0.999, -1.957, 0.697, 0.697, 1.0, 14.0, 0.0}}, + {"bp", {1.0, 1.139, 1.139, 1.683, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"blyp", {1.0, 1.094, 1.094, 1.682, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"revpbe", {1.0, 0.923, 0.923, 1.01, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"rpbe", {1.0, 0.872, 0.872, 0.514, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b97_d", {1.0, 0.892, 0.892, 0.909, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b973c", {1.0, 1.06, 1.06, 1.5, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pbe", {1.0, 1.217, 1.217, 0.722, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"rpw86pbe", {1.0, 1.224, 1.224, 0.901, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b3lyp", {1.0, 1.261, 1.261, 1.703, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tpss", {1.0, 1.166, 1.166, 1.105, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hf", {1.0, 1.158, 1.158, 1.746, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tpss0", {1.0, 1.252, 1.252, 1.242, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pbe0", {1.0, 1.287, 1.287, 0.928, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hse06", {1.0, 1.129, 1.129, 0.109, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hse", {1.0, 1.129, 1.129, 0.109, 1.0, 1.0, 1.0, 14.0, 0.0}}, // ABACUS implements HSE06 as HSE + {"revpbe38", {1.0, 1.021, 1.021, 0.862, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pw6b95", {1.0, 1.532, 1.532, 0.862, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b2plyp", {0.64, 1.427, 1.427, 1.022, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"dsdblyp", {0.5, 1.569, 1.569, 0.705, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpwlyp", {1.0, 1.239, 1.239, 1.098, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"olyp", {1.0, 0.806, 0.806, 1.764, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"bpbe", {1.0, 1.087, 1.087, 2.033, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"opbe", {1.0, 0.837, 0.837, 2.033, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"ssb", {1.0, 1.215, 1.215, 0.663, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"revssb", {1.0, 1.221, 1.221, 0.56, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"otpss", {1.0, 1.128, 1.128, 1.494, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b3pw91", {1.0, 1.176, 1.176, 1.775, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"bhlyp", {1.0, 1.37, 1.37, 1.442, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tpssh", {1.0, 1.223, 1.223, 1.219, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpw1b95", {1.0, 1.605, 1.605, 1.118, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pwb6k", {1.0, 1.66, 1.66, 0.55, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b1b95", {1.0, 1.613, 1.613, 1.868, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"bmk", {1.0, 1.931, 1.931, 2.168, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"camb3lyp", {1.0, 1.378, 1.378, 1.217, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"lcwpbe", {1.0, 1.355, 1.355, 1.279, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b2gpplyp", {0.56, 1.586, 1.586, 0.76, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"ptpss", {0.75, 1.541, 1.541, 0.879, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pwpb95", {0.82, 1.557, 1.557, 0.705, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pw1pw", {1.0, 1.4968, 1.4968, 1.1786, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"scan", {1.0, 1.324, 1.324, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"wb97x_d3", {1.0, 1.281, 1.281, 1.0, 1.094, 1.094, 1.0, 14.0, 0.0}}, + // NOTE: simple-dftd3 assign the D3(0) parameters of functional wB97X-D3 + // to a key `wb97x`, but the functional wB97X itself does not own these params. + // instead, there is a XC in libxc really names HYB_GGA_WB97X_D3 + {"pbehpbe", {1.0, 1.5703, 1.5703, 1.401, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"xlyp", {1.0, 0.9384, 0.9384, 0.7447, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpwpw", {1.0, 1.3725, 1.3725, 1.9467, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hcth407", {1.0, 4.0426, 4.0426, 2.7694, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"revtpss", {1.0, 1.3491, 1.3491, 1.3666, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tauhcth", {1.0, 0.932, 0.932, 0.5662, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b3p", {1.0, 1.1897, 1.1897, 1.1961, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b1p", {1.0, 1.1815, 1.1815, 1.1209, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b1lyp", {1.0, 1.3725, 1.3725, 1.9467, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpwb1k", {1.0, 1.671, 1.671, 1.061, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpw1lyp", {1.0, 2.0512, 2.0512, 1.9529, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpw1pw", {1.0, 1.2892, 1.2892, 1.4758, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpw1kcis", {1.0, 1.7231, 1.7231, 2.2917, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpwkcis1k", {1.0, 1.4853, 1.4853, 1.7553, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pbeh1pbe", {1.0, 1.3719, 1.3719, 1.043, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pbe1kcis", {1.0, 3.6355, 3.6355, 1.7934, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"x3lyp", {1.0, 1.0, 1.0, 0.299, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"o3lyp", {1.0, 1.406, 1.406, 1.8058, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b97_1", {1.0, 3.7924, 3.7924, 1.6418, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b97_2", {1.0, 1.7066, 1.7066, 1.6418, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"b98", {1.0, 2.6895, 2.6895, 1.9078, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hiss", {1.0, 1.3338, 1.3338, 0.7615, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"hse03", {1.0, 1.3944, 1.3944, 1.0156, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"revtpssh", {1.0, 1.3224, 1.3224, 1.2504, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"revtpss0", {1.0, 1.2881, 1.2881, 1.0649, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tpss1kcis", {1.0, 1.7729, 1.7729, 2.0902, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"tauhcthhyb", {1.0, 1.5001, 1.5001, 1.6302, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pkzb", {1.0, 0.6327, 0.6327, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"n12", {1.0, 1.3493, 1.3493, 2.3916, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mpw2plyp", {0.66, 1.5527, 1.5527, 0.7529, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m05", {1.0, 1.373, 1.373, 0.595, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m052x", {1.0, 1.417, 1.417, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m06l", {1.0, 1.581, 1.581, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m06", {1.0, 1.325, 1.325, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m062x", {1.0, 1.619, 1.619, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m08hx", {1.0, 1.6247, 1.6247, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"m11l", {1.0, 2.3933, 2.3933, 1.1129, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"mn15l", {1.0, 3.3388, 3.3388, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"pwp", {1.0, 2.104, 2.104, 0.8747, 1.0, 1.0, 1.0, 14.0, 0.0}}, + }; + // DFT-D3M(BJ): not implemented for beta + std::map> bjm = { + {"bp", {1.0, 0.82185, 0.82185, 3.140281, 2.728151, 2.728151, 1.0, 14.0, 0.0}}, + {"blyp", {1.0, 0.448486, 0.448486, 1.875007, 3.610679, 3.610679, 1.0, 14.0, 0.0}}, + {"b97_d", {1.0, 0.240184, 0.240184, 1.206988, 3.864426, 3.864426, 1.0, 14.0, 0.0}}, + {"pbe", {1.0, 0.012092, 0.012092, 0.35894, 5.938951, 5.938951, 1.0, 14.0, 0.0}}, + {"b3lyp", {1.0, 0.278672, 0.278672, 1.466677, 4.606311, 4.606311, 1.0, 14.0, 0.0}}, + {"pbe0", {1.0, 0.007912, 0.007912, 0.528823, 6.162326, 6.162326, 1.0, 14.0, 0.0}}, + {"b2plyp", {0.64, 0.486434, 0.486434, 0.67282, 3.656466, 3.656466, 1.0, 14.0, 0.0}}, + {"lcwpbe", {1.0, 0.563761, 0.563761, 0.906564, 3.59368, 3.59368, 1.0, 14.0, 0.0}}, + }; + // DFT-D3M(0): not implemented for beta + std::map> zerom = { + {"bp", {1.0, 1.23346, 1.23346, 1.945174, 1.0, 1.0, 1.0, 14.0, 0.0}}, + {"blyp", {1.0, 1.279637, 1.279637, 1.841686, 1.0, 1.0, 1.0, 14.0, 0.01437}}, + {"b97_d", {1.0, 1.151808, 1.151808, 1.020078, 1.0, 1.0, 1.0, 14.0, 0.035964}}, + {"pbe", {1.0, 2.340218, 2.340218, 0.0, 1.0, 1.0, 1.0, 14.0, 0.129434}}, + {"b3lyp", {1.0, 1.338153, 1.338153, 1.532981, 1.0, 1.0, 1.0, 14.0, 0.013988}}, + {"pbe0", {1.0, 2.077949, 2.077949, 8.1e-05, 1.0, 1.0, 1.0, 14.0, 0.116755}}, + {"b2plyp", {0.64, 1.313134, 1.313134, 0.717543, 1.0, 1.0, 1.0, 14.0, 0.016035}}, + {"lcwpbe", {1.0, 1.366361, 1.366361, 1.280619, 1.0, 1.0, 1.0, 14.0, 0.00316}}, + }; + // DFT-D3(OptimizedPower) + std::map> op = { + {"blyp", {1.0, 0.425, 0.425, 1.31867, 3.5, 3.5, 1.0, 14.0, 2.0}}, + {"revpbe", {1.0, 0.6, 0.6, 1.44765, 2.5, 2.5, 1.0, 14.0, 0.0}}, + {"b97_d", {1.0, 0.6, 0.6, 1.46861, 2.5, 2.5, 1.0, 14.0, 0.0}}, + {"pbe", {0.91826, 0.2, 0.2, 0.0, 4.75, 4.75, 1.0, 14.0, 6.0}}, + {"b3lyp", {1.0, 0.3, 0.3, 0.78311, 4.25, 4.25, 1.0, 14.0, 4.0}}, + {"tpss", {1.0, 0.575, 0.575, 0.51581, 3.0, 3.0, 1.0, 14.0, 8.0}}, + {"pbe0", {0.8829, 0.15, 0.15, 0.0, 4.75, 4.75, 1.0, 14.0, 6.0}}, + {"revpbe0", {1.0, 0.725, 0.725, 1.25684, 2.25, 2.25, 1.0, 14.0, 0.0}}, + {"tpssh", {1.0, 0.575, 0.575, 0.43185, 3.0, 3.0, 1.0, 14.0, 8.0}}, + {"revtpss", {1.0, 0.7, 0.7, 0.27632, 2.5, 2.5, 1.0, 14.0, 8.0}}, + {"b97_1", {0.97388, 0.15, 0.15, 0.0, 4.25, 4.25, 1.0, 14.0, 6.0}}, + {"revtpssh", {1.0, 0.575, 0.575, 0.12467, 3.0, 3.0, 1.0, 14.0, 10.0}}, + {"ms2", {1.0, 0.7, 0.7, 0.90743, 4.0, 4.0, 1.0, 14.0, 2.0}}, + {"ms2h", {1.0, 0.65, 0.65, 1.69464, 4.75, 4.75, 1.0, 14.0, 0.0}}, + }; + // a declaration + std::string _lowercase(const std::string& s); + /** + * @brief Get the dftd3 params object. + * dftd3 method fall back: xc-bjm -> xc-bj -> pbe-bj + * xc-zerom -> xc-zero -> pbe-zero + * + * @param xc the functional name + * @param d3method the d3 method, can be "bj", "zero-damping", "bj-modified", "zero-damping-modified", "op" + * @param param the dftd3 parameters, ALL_KEYS = {'s6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet'} + */ + void _search(const std::string& xc, + const std::string& method, + std::vector& param) + { + const std::string xc_lowercase = _lowercase(xc); + const std::vector allowed_ = { "bj", "zero", "bjm", "zerom", "op" }; + assert(std::find(allowed_.begin(), allowed_.end(), method) != allowed_.end()); + if (method == "op") + { + if (op.find(xc_lowercase) != op.end()) + { + param = op[xc_lowercase]; + } + else + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "XC (`" + xc + "`)'s DFT-D3(OP) parameters not found"); + } + } + else if (method == "bjm") + { + if (bjm.find(xc_lowercase) != bjm.end()) + { + param = bjm[xc_lowercase]; + } + else + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "XC (`" + xc + "`)'s DFT-D3M(BJ) parameters not found"); + } + } + else if (method == "bj") + { + if (bj.find(xc_lowercase) != bj.end()) + { + param = bj[xc_lowercase]; + } + else + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "XC (`" + xc + "`)'s DFT-D3(BJ) parameters not found"); + } + } + else if (method == "zerom") + { + if (zerom.find(xc_lowercase) != zerom.end()) + { + param = zerom[xc_lowercase]; + } + else + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "XC (`" + xc + "`)'s DFT-D3M(0) parameters not found"); + } + } + else if (method == "zero") + { + if (zero.find(xc_lowercase) != zero.end()) + { + param = zero[xc_lowercase]; + } + else + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "XC (`" + xc + "`)'s DFT-D3(0) parameters not found"); + } + } + else // should not reach here + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "Unknown DFT-D3 method: " + method); + } + } + + std::string _lowercase(const std::string& s) + { + std::string result = s; + std::transform(s.begin(), s.end(), result.begin(), ::tolower); + return result; + } + + std::string _uppercase(const std::string& s) + { + std::string result = s; + std::transform(s.begin(), s.end(), result.begin(), ::toupper); + return result; + } + + /** + * @brief Get DFT-D3 parameters. If if there are parameters defined, + * then it will overwrite the search result. If all parameters are + * defined already by user, then search will not performed. + * + * @param xc XC functional name + * @param d3method can be "d3_0" or "d3_bj" + * @param s6_in user defined s6, default is "default" + * @param s8_in user defined s8, default is "default" + * @param a1_in user defined a1, default is "default" + * @param a2_in user defined a2, default is "default" + * @param s6 [out] s6 parameter + * @param s8 [out] s8 parameter + * @param a1 [out] a1 parameter + * @param a2 [out] a2 parameter + */ + void dftd3_params(const std::string& xc_in, + const std::string& d3method, + const std::string& s6_in, + const std::string& s8_in, + const std::string& a1_in, + const std::string& a2_in, + double& s6, + double& s8, + double& a1, + double& a2, + std::ofstream* plog = nullptr) + { + const std::map param_map = { + {"d3_bj", "bj"}, {"d3_0", "zero"}, {"d3_bjm", "bjm"}, {"d3_0m", "zerom"}, + {"op", "op"}}; + + const std::vector flag = {s6_in, s8_in, a1_in, a2_in}; + const bool autoset = std::any_of(flag.begin(), flag.end(), [](const std::string& s) { return s == "default"; }); + if (!autoset) // all parameters are defined + { + s6 = std::stod(s6_in); + s8 = std::stod(s8_in); + a1 = std::stod(a1_in); + a2 = std::stod(a2_in); + } + else + { + std::vector param; + const std::string xc = DFTD3::_xcname(xc_in); + _search(xc, param_map.at(d3method), param); + s6 = (s6_in == "default") ? param[0] : std::stod(s6_in); + s8 = (s8_in == "default") ? param[3] : std::stod(s8_in); + a1 = (a1_in == "default") ? param[2] : std::stod(a1_in); + a2 = (a2_in == "default") ? param[5] : std::stod(a2_in); + if (plog != nullptr) // logging the autoset + { + param = {s6, s8, a1, a2}; + FmtTable vdwd3tab({"Parameters", "Original", "Autoset"}, 4, {"%10s", "%10s", "%10.4f"}); + const std::vector items = {"s6", "s8", "a1", "a2"}; + vdwd3tab << items << flag << param; + (*plog) << "\nDFT-D3 Dispersion correction parameters autoset\n" << vdwd3tab.str() + << "XC functional: " << xc_in << std::endl; + } + + } + } +} + +/* +''' +dftd3 parameters from +https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml + +this script is to convert the toml file to c++ map +''' + +import toml + +def load(fn): + with open(fn, 'r') as f: + data = toml.load(f) + return data + +def xc_indexing(data): + out = {'bj': {}, 'zero': {}, 'bjm': {}, 'zerom': {}, 'op': {}} + for xc, param in data['parameter'].items(): + for vdw_method, value in param['d3'].items(): + out[vdw_method][xc] = {k: v for k, v in value.items() if k != 'doi'} + return out + +def complete(vdw_method, value): + ''' + for each functional, the zero damping version must be provided + for each vdw method, all parameters including + s6, rs6/a1, s8, rs8/a2, s9, alp, bet must be provided, otherwise + use the default value + ''' + DEFAULT = { + 'bj': {'s6': 1.0, 's9': 1.0, 'alp': 14.0}, + 'zero': {'s6': 1.0, 's9': 1.0, 'rs8': 1.0, 'alp': 14.0}, + 'bjm': {'s6': 1.0, 's9': 1.0, 'alp': 14.0}, + 'zerom': {'s6': 1.0, 's9': 1.0, 'rs8': 1.0, 'alp': 14.0}, + 'op': {'s9': 1.0, 'alp': 14.0} + } + ALL_KEYS = {'s6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet'} + EQUIVALENT = {'rs6': 'a1', 'a1': 'rs6', 'rs8': 'a2', 'a2': 'rs8'} + out = value.copy() + for k in ALL_KEYS: + equilk = EQUIVALENT.get(k, k) + val = [out.get(k), out.get(equilk), + DEFAULT[vdw_method].get(k), DEFAULT[vdw_method].get(equilk)] + val = [v for v in val if v is not None] + val = [0.0] if not val else val + out[k] = val[0] + out[equilk] = out[k] + # equivalent? + # according to + # abacus-develop/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp + # https://abacus.deepmodeling.com/en/latest/advanced/input_files/input-main.html + + return out + +def make_stdmap(data): + for vdw_method, param in data.items(): + print(f'std::map> {vdw_method} = {{') + for xc, value in param.items(): + print(f' {{\"{xc}\", {{', end='') + print(', '.join([f'{v}' for v in value.values()]), end='') + print('}},') + print('};') + +if __name__ == '__main__': + fn = 'dftd3.toml' + data = load(fn) + data = xc_indexing(data) + for vdw_method, param in data.items(): + for xc, value in param.items(): + raw = complete(vdw_method, value) + data[vdw_method][xc] = {k: raw[k] + for k in ['s6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet']} + make_stdmap(data) +*/ +#endif // DFTD3_XC_PARAM_H diff --git a/source/module_hamilt_general/module_vdw/test/CMakeLists.txt b/source/module_hamilt_general/module_vdw/test/CMakeLists.txt index 7084800956..8833b3cf08 100644 --- a/source/module_hamilt_general/module_vdw/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_vdw/test/CMakeLists.txt @@ -9,4 +9,10 @@ AddTest( TARGET vdwTest LIBS parameter ${math_libs} base device vdw SOURCES vdw_test.cpp +) + +AddTest( + TARGET dftd3_xc + LIBS parameter ${math_libs} base device + SOURCES dftd3_xc_test.cpp ) \ No newline at end of file diff --git a/source/module_hamilt_general/module_vdw/test/dftd3_xc_test.cpp b/source/module_hamilt_general/module_vdw/test/dftd3_xc_test.cpp new file mode 100644 index 0000000000..3ed3a42146 --- /dev/null +++ b/source/module_hamilt_general/module_vdw/test/dftd3_xc_test.cpp @@ -0,0 +1,104 @@ +#include +#include "module_hamilt_general/module_vdw/dftd3_xc_name.h" +#include "module_hamilt_general/module_vdw/dftd3_xc_param.h" + +TEST(DFTD3XCTest, SearchXcnameLibXCXplusC) +{ + std::string xname; + DFTD3::_xcname_libxc_xplusc("XC_GGA_X_PBE+XC_GGA_C_OP_PBE", xname); + EXPECT_EQ(xname, "pbeop"); + // then test the case with out prefix XC_ + DFTD3::_xcname_libxc_xplusc("GGA_X_PBE+GGA_C_OP_PBE", xname); + EXPECT_EQ(xname, "pbeop"); +} + +TEST(DFTD3XCTest, SearchXcnameLibXCXC) +{ + std::string xname; + DFTD3::_xcname_libxc_xc("XC_LDA_XC_TETER93", xname); + EXPECT_EQ(xname, "teter93"); + // then test the case with out prefix XC_ + DFTD3::_xcname_libxc_xc("LDA_XC_TETER93", xname); + EXPECT_EQ(xname, "teter93"); +} + +TEST(DFTD3XCTest, SearchXcnameLibXC) +{ + std::string xname; + DFTD3::_xcname_libxc("XC_GGA_X_PBE+XC_GGA_C_OP_PBE", xname); + EXPECT_EQ(xname, "pbeop"); + // then test the case with out prefix XC_ + DFTD3::_xcname_libxc("GGA_X_PBE+GGA_C_OP_PBE", xname); + EXPECT_EQ(xname, "pbeop"); +} + +TEST(DFTD3XCTest, SearchXcname) +{ + std::string xcpattern = "XC_GGA_X_PBE+XC_GGA_C_OP_PBE"; + std::string xname = DFTD3::_xcname(xcpattern); + EXPECT_EQ(xname, "pbeop"); + + xcpattern = "XC_LDA_XC_TETER93"; + xname = DFTD3::_xcname(xcpattern); + EXPECT_EQ(xname, "teter93"); + + xcpattern = "default"; + xname = DFTD3::_xcname(xcpattern); + EXPECT_EQ(xname, "default"); + + xcpattern = "pbe"; + xname = DFTD3::_xcname(xcpattern); + EXPECT_EQ(xname, "pbe"); +} + +TEST(DFTD3XCTest, SuccessfulSearch) +{ + std::string xc = "pbe"; + std::string d3method = "d3_0"; + std::string s6_in = "default"; + std::string s8_in = "default"; + std::string a1_in = "default"; + std::string a2_in = "default"; + double s6, s8, a1, a2; + DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2); + EXPECT_DOUBLE_EQ(s6, 1.0); + EXPECT_DOUBLE_EQ(s8, 0.722); + EXPECT_DOUBLE_EQ(a1, 1.217); + EXPECT_DOUBLE_EQ(a2, 1.0); + + // a more complicated case: MGGA_X_SCAN+MGGA_C_SCAN + xc = "XC_MGGA_X_SCAN+XC_MGGA_C_SCAN"; + DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2); + EXPECT_DOUBLE_EQ(s6, 1.0); + EXPECT_DOUBLE_EQ(s8, 0.0); + EXPECT_DOUBLE_EQ(a1, 1.324); + EXPECT_DOUBLE_EQ(a2, 1.0); + + // user defines all parameters + s6_in = "1.1"; + s8_in = "0.1"; + a1_in = "1.325"; + a2_in = "1.1"; + DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2); + EXPECT_DOUBLE_EQ(s6, 1.1); + EXPECT_DOUBLE_EQ(s8, 0.1); + EXPECT_DOUBLE_EQ(a1, 1.325); + EXPECT_DOUBLE_EQ(a2, 1.1); + + // user defines one parameter + s6_in = "1.1"; + s8_in = "default"; + a1_in = "default"; + a2_in = "default"; + DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2); + EXPECT_DOUBLE_EQ(s6, 1.1); + EXPECT_DOUBLE_EQ(s8, 0.0); + EXPECT_DOUBLE_EQ(a1, 1.324); + EXPECT_DOUBLE_EQ(a2, 1.0); +} + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} \ No newline at end of file diff --git a/source/module_hamilt_general/module_vdw/vdw.cpp b/source/module_hamilt_general/module_vdw/vdw.cpp index 4644ab7035..b070d527ae 100644 --- a/source/module_hamilt_general/module_vdw/vdw.cpp +++ b/source/module_hamilt_general/module_vdw/vdw.cpp @@ -6,30 +6,35 @@ namespace vdw { -std::unique_ptr make_vdw(const UnitCell &ucell, const Input_para &input) +std::unique_ptr make_vdw(const UnitCell &ucell, + const Input_para &input, + std::ofstream* plog) { - if (ucell.nat < 2 && input.vdw_method != "none") - { - ModuleBase::WARNING("VDW", "Only one atom in this system, and will not do the calculation of VDW"); - return nullptr; - } - else if (input.vdw_method == "d2") + // if (ucell.nat < 2 && input.vdw_method != "none") + // { + // ModuleBase::WARNING("VDW", "Only one atom in this system, and will not do the calculation of VDW"); + // return nullptr; + // } + if (input.vdw_method == "d2") { std::unique_ptr vdw_ptr = make_unique(ucell); - vdw_ptr->parameter().initial_parameters(input); + vdw_ptr->parameter().initial_parameters(input, plog); vdw_ptr->parameter().initset(ucell); return vdw_ptr; } else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj") { std::unique_ptr vdw_ptr = make_unique(ucell); - vdw_ptr->parameter().initial_parameters(input); + vdw_ptr->parameter().initial_parameters(input, plog); return vdw_ptr; } - else + else if (input.vdw_method != "none") { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", + "Unrecognized Van der Waals correction method: " + input.vdw_method); return nullptr; } + return nullptr; // "none" method } } // namespace vdw \ No newline at end of file diff --git a/source/module_hamilt_general/module_vdw/vdw.h b/source/module_hamilt_general/module_vdw/vdw.h index 94b54239cb..257b80cd98 100644 --- a/source/module_hamilt_general/module_vdw/vdw.h +++ b/source/module_hamilt_general/module_vdw/vdw.h @@ -48,7 +48,17 @@ class Vdw virtual void cal_stress() = 0; }; -std::unique_ptr make_vdw(const UnitCell &ucell, const Input_para &input); +/** + * @brief make vdw correction object + * + * @param ucell UnitCell instance + * @param input Parameter instance + * @param plog optional, for logging the parameter setting process + * @return std::unique_ptr + */ +std::unique_ptr make_vdw(const UnitCell &ucell, + const Input_para &input, + std::ofstream* plog = nullptr); } // namespace vdw diff --git a/source/module_hamilt_general/module_vdw/vdwd2_parameters.cpp b/source/module_hamilt_general/module_vdw/vdwd2_parameters.cpp index de8bda7bcc..a562244be3 100644 --- a/source/module_hamilt_general/module_vdw/vdwd2_parameters.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd2_parameters.cpp @@ -9,7 +9,7 @@ namespace vdw { -void Vdwd2Parameters::initial_parameters(const Input_para &input) +void Vdwd2Parameters::initial_parameters(const Input_para &input, std::ofstream* plog) { scaling_ = std::stod(input.vdw_s6); damping_ = input.vdw_d; diff --git a/source/module_hamilt_general/module_vdw/vdwd2_parameters.h b/source/module_hamilt_general/module_vdw/vdwd2_parameters.h index ef981b2c0a..06d4a52793 100644 --- a/source/module_hamilt_general/module_vdw/vdwd2_parameters.h +++ b/source/module_hamilt_general/module_vdw/vdwd2_parameters.h @@ -31,7 +31,15 @@ class Vdwd2Parameters : public VdwParameters void R0_input(const std::string &file, const std::string &unit); void initset(const UnitCell &ucell); // init sets of vdwd2 once this correction is called - void initial_parameters(const Input_para &input); // initial parameters of Vdwd2 with INPUT file + + /** + * @brief initial parameters of Vdwd2 with INPUT file + * + * @param input Parameter instance + * @param plog optional, for logging the parameter setting process (not implemented) + */ + void initial_parameters(const Input_para &input, + std::ofstream* plog = nullptr); inline const std::map C6() const { return C6_; } inline const std::map R0() const { return R0_; } diff --git a/source/module_hamilt_general/module_vdw/vdwd3.cpp b/source/module_hamilt_general/module_vdw/vdwd3.cpp index be68c13651..b5b0798c63 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd3.cpp @@ -24,29 +24,33 @@ void Vdwd3::init() std::vector at_kind = atom_kind(); iz_.reserve(ucell_.nat); xyz_.reserve(ucell_.nat); - for (size_t it = 0; it != ucell_.ntype; it++) + for (size_t it = 0; it != ucell_.ntype; it++) { for (size_t ia = 0; ia != ucell_.atoms[it].na; ia++) { iz_.emplace_back(at_kind[it]); xyz_.emplace_back(ucell_.atoms[it].tau[ia] * ucell_.lat0); } +} std::vector tau_max(3); if (para_.model() == "radius") { rep_vdw_.resize(3); set_criteria(para_.rthr2(), lat_, tau_max); - for (size_t i = 0; i < 3; i++) + for (size_t i = 0; i < 3; i++) { rep_vdw_[i] = std::ceil(tau_max[i]); +} } - else if (para_.model() == "period") + else if (para_.model() == "period") { rep_vdw_ = {para_.period().x, para_.period().y, para_.period().z}; +} rep_cn_.resize(3); set_criteria(para_.cn_thr2(), lat_, tau_max); - for (size_t i = 0; i < 3; i++) + for (size_t i = 0; i < 3; i++) { rep_cn_[i] = ceil(tau_max[i]); } +} void Vdwd3::set_criteria(double rthr, const std::vector> &lat, std::vector &tau_max) { @@ -66,13 +70,15 @@ void Vdwd3::set_criteria(double rthr, const std::vector Vdwd3::atom_kind() { std::vector atom_kind(ucell_.ntype); - for (size_t i = 0; i != ucell_.ntype; i++) - for (int j = 0; j != ModuleBase::element_name.size(); j++) + for (size_t i = 0; i != ucell_.ntype; i++) { + for (int j = 0; j != ModuleBase::element_name.size(); j++) { if (ucell_.atoms[i].ncpp.psd == ModuleBase::element_name[j]) { atom_kind[i] = j; break; } +} +} return atom_kind; } @@ -91,24 +97,25 @@ void Vdwd3::cal_energy() if (para_.version() == "d3_0") // DFT-D3(zero-damping) { double tmp; - for (int iat = 0; iat != ucell_.nat - 1; iat++) + for (int iat = 0; iat != ucell_.nat - 1; iat++) { for (int jat = iat + 1; jat != ucell_.nat; jat++) { get_c6(iz_[iat], iz_[jat], cn[iat], cn[jat], c6); - if (para_.abc()) + if (para_.abc()) // three-body term { ij = lin(iat, jat); cc6ab[ij] = std::sqrt(c6); } - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; - r2 = (xyz_[iat] - xyz_[jat] + tau).norm2(); - if (r2 > para_.rthr2()) + r2 = (xyz_[iat] - xyz_[jat] + tau).norm2(); // |r+T|^2 + if (r2 > para_.rthr2()) { // neglect the distance larger than rthr2 continue; +} rr = para_.r0ab()[iz_[iat]][iz_[jat]] / std::sqrt(r2); // zero-damping function tmp = para_.rs6() * rr; @@ -123,7 +130,10 @@ void Vdwd3::cal_energy() r8 = r6 * r2; e8 += c8 * damp8 / r8; } // end tau +} +} } // end jat +} for (int iat = 0; iat != ucell_.nat; iat++) { @@ -134,17 +144,19 @@ void Vdwd3::cal_energy() ij = lin(iat, jat); cc6ab[ij] = std::sqrt(c6); } - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { - if (taux == 0 && tauy == 0 && tauz == 0) + if (taux == 0 && tauy == 0 && tauz == 0) { continue; +} tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = tau.norm2(); - if (r2 > para_.rthr2()) + if (r2 > para_.rthr2()) { continue; +} rr = para_.r0ab()[iz_[iat]][iz_[jat]] / std::sqrt(r2); // zero-damping function @@ -160,6 +172,8 @@ void Vdwd3::cal_energy() r8 = r6 * r2; e8 += c8 * damp8 / r8 * 0.5; } // end tau +} +} } // end iat } // end d3_0 else if (para_.version() == "d3_bj") // DFT-D3(BJ-damping) @@ -175,15 +189,16 @@ void Vdwd3::cal_energy() ij = lin(iat, jat); cc6ab[ij] = std::sqrt(c6); } - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = (xyz_[iat] - xyz_[jat] + tau).norm2(); - if (r2 > para_.rthr2()) + if (r2 > para_.rthr2()) { continue; +} rr = para_.r0ab()[iz_[iat]][iz_[jat]] / std::sqrt(r2); // BJ-damping function @@ -198,6 +213,8 @@ void Vdwd3::cal_energy() r8 = r6 * r2; e8 += c8 / (r8 + damp8); } // end tau +} +} } // end jat int jat = iat; get_c6(iz_[iat], iz_[jat], cn[iat], cn[jat], c6); @@ -209,17 +226,19 @@ void Vdwd3::cal_energy() ij = lin(iat, jat); cc6ab[ij] = std::sqrt(c6); } - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { - if (taux == 0 && tauy == 0 && tauz == 0) + if (taux == 0 && tauy == 0 && tauz == 0) { continue; +} tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = tau.norm2(); - if (r2 > para_.rthr2()) + if (r2 > para_.rthr2()) { continue; +} rr = para_.r0ab()[iz_[iat]][iz_[jat]] / std::sqrt(r2); r6 = std::pow(r2, 3); @@ -229,6 +248,8 @@ void Vdwd3::cal_energy() r8 = r6 * r2; e8 += c8 / (r8 + damp8) * 0.5; } // end tau +} +} } // end iat } // end d3_bj @@ -256,8 +277,9 @@ void Vdwd3::cal_force() pbc_gdisp(g, smearing_sigma); - for (size_t iat = 0; iat != ucell_.nat; iat++) + for (size_t iat = 0; iat != ucell_.nat; iat++) { force_[iat] = -2.0 * g[iat]; +} ModuleBase::timer::tick("Vdwd3", "cal_force"); } @@ -286,7 +308,7 @@ void Vdwd3::get_c6(int iat, int jat, double nci, double ncj, double &c6) { double c6mem = -1e99, rsum = 0.0, csum = 0.0, r_save = 1e99; double cn1, cn2, r; - for (size_t i = 0; i != para_.mxc()[iat]; i++) + for (size_t i = 0; i != para_.mxc()[iat]; i++) { for (size_t j = 0; j != para_.mxc()[jat]; j++) { c6 = para_.c6ab()[0][j][i][jat][iat]; @@ -305,6 +327,7 @@ void Vdwd3::get_c6(int iat, int jat, double nci, double ncj, double &c6) csum += tmp1 * c6; } } +} c6 = (rsum > 1e-99) ? csum / rsum : c6mem; } @@ -315,21 +338,26 @@ void Vdwd3::pbc_ncoord(std::vector &cn) double xn = 0.0; ModuleBase::Vector3 tau; double r2, rr; - for (size_t iat = 0; iat != ucell_.nat; iat++) - for (int taux = -rep_cn_[0]; taux <= rep_cn_[0]; taux++) - for (int tauy = -rep_cn_[1]; tauy <= rep_cn_[1]; tauy++) + for (size_t iat = 0; iat != ucell_.nat; iat++) { + for (int taux = -rep_cn_[0]; taux <= rep_cn_[0]; taux++) { + for (int tauy = -rep_cn_[1]; tauy <= rep_cn_[1]; tauy++) { for (int tauz = -rep_cn_[2]; tauz <= rep_cn_[2]; tauz++) { - if (iat == i && taux == 0 && tauy == 0 && tauz == 0) + if (iat == i && taux == 0 && tauy == 0 && tauz == 0) { continue; +} tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = (xyz_[iat] - xyz_[i] + tau).norm2(); - if (r2 > para_.cn_thr2()) + if (r2 > para_.cn_thr2()) { continue; +} rr = (para_.rcov()[iz_[i]] + para_.rcov()[iz_[iat]]) / std::sqrt(r2); xn += 1.0 / (1.0 + exp(-para_.k1() * (rr - 1.0))); } +} +} +} cn[i] = xn; } } @@ -346,7 +374,7 @@ void Vdwd3::pbc_three_body(const std::vector &iz, double r0ij, r0ik, r0jk, c9, rij2, rik2, rjk2, rr0ij, rr0ik, rr0jk, geomean, fdamp, tmp1, tmp2, tmp3, tmp4, ang; ModuleBase::Vector3 ijvec, ikvec, jkvec, jtau, ktau; std::vector repmin(3), repmax(3); - for (int iat = 2; iat != ucell_.nat; iat++) + for (int iat = 2; iat != ucell_.nat; iat++) { for (int jat = 1; jat != iat; jat++) { ijvec = xyz_[jat] - xyz_[iat]; @@ -378,25 +406,28 @@ void Vdwd3::pbc_three_body(const std::vector &iz, jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = (ijvec + jtau).norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / r0ij; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / r0ik; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / r0jk; geomean = std::pow(rr0ij * rr0ik * rr0jk, 1.0 / 3.0); @@ -410,11 +441,14 @@ void Vdwd3::pbc_three_body(const std::vector &iz, eabc += ang * c9 * fdamp; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux } // end kat } // end jat +} // end iat for (int iat = 1; iat != ucell_.nat; iat++) @@ -445,29 +479,33 @@ void Vdwd3::pbc_three_body(const std::vector &iz, { repmin[2] = std::max(-rep_cn_[2], jtauz - rep_cn_[2]); repmax[2] = std::min(rep_cn_[2], jtauz + rep_cn_[2]); - if (jtaux == 0 && jtauy == 0 && jtauz == 0) + if (jtaux == 0 && jtauy == 0 && jtauz == 0) { continue; +} jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = (ijvec + jtau).norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / r0ij; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / r0ik; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / r0jk; geomean = std::pow(rr0ij * rr0ik * rr0jk, 1.0 / 3.0); @@ -481,13 +519,15 @@ void Vdwd3::pbc_three_body(const std::vector &iz, eabc += ang * c9 * fdamp / 2.0; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux } // end kat } // end iat - for (int iat = 1; iat != ucell_.nat; iat++) + for (int iat = 1; iat != ucell_.nat; iat++) { for (int jat = 0; jat != iat; jat++) { int kat = jat; @@ -518,26 +558,30 @@ void Vdwd3::pbc_three_body(const std::vector &iz, jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = (ijvec + jtau).norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / r0ij; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { - if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) + if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) { continue; +} ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / r0ik; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / r0jk; geomean = std::pow(rr0ij * rr0ik * rr0jk, 1.0 / 3.0); @@ -551,10 +595,13 @@ void Vdwd3::pbc_three_body(const std::vector &iz, eabc += ang * c9 * fdamp / 2.0; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux } // end jat +} // end iat for (int iat = 0; iat != ucell_.nat; iat++) @@ -587,31 +634,37 @@ void Vdwd3::pbc_three_body(const std::vector &iz, repmax[2] = std::min(rep_cn_[2], jtauz + rep_cn_[2]); jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; - if (jtaux == 0 && jtauy == 0 && jtauz == 0) + if (jtaux == 0 && jtauy == 0 && jtauz == 0) { continue; +} rij2 = jtau.norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / r0ij; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { - if (ktaux == 0 && ktauy == 0 && ktauz == 0) + if (ktaux == 0 && ktauy == 0 && ktauz == 0) { continue; - if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) +} + if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) { continue; +} ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = ktau.norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / r0ik; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / r0jk; geomean = std::pow(rr0ij * rr0ik * rr0jk, 1.0 / 3.0); @@ -625,6 +678,8 @@ void Vdwd3::pbc_three_body(const std::vector &iz, eabc += ang * c9 * fdamp / 6.0; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux @@ -637,7 +692,7 @@ void Vdwd3::get_dc6_dcnij(int mxci, int mxcj, double cni, double cnj, int izi, i double r_save = 9999.0, c6mem = -1e99, zaehler = 0.0, nenner = 0.0; double dzaehler_i = 0.0, dnenner_i = 0.0, dzaehler_j = 0.0, dnenner_j = 0.0; double c6ref = 0.0, cn_refi = 0.0, cn_refj = 0.0, r = 0.0, expterm = 0.0, term = 0.0; - for (size_t a = 0; a != mxci; a++) + for (size_t a = 0; a != mxci; a++) { for (size_t b = 0; b != mxcj; b++) { c6ref = para_.c6ab()[0][b][a][izj][izi]; @@ -664,6 +719,7 @@ void Vdwd3::get_dc6_dcnij(int mxci, int mxcj, double cni, double cnj, int izi, i dnenner_j += term; } } +} if (nenner > 1e-99) { @@ -710,8 +766,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m r42 = para_.r2r4()[iz_[iat]] * para_.r2r4()[iz_[iat]]; rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[iat]]; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] @@ -748,6 +804,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc6_rest_sum[linii] += dc6_rest; } } // end tau +} +} for (int jat = 0; jat != iat; jat++) { get_dc6_dcnij(para_.mxc()[iz_[iat]], para_.mxc()[iz_[jat]], cn[iat], cn[jat], @@ -760,15 +818,16 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[jat]]; dc6ij[jat][iat] = dc6iji; dc6ij[iat][jat] = dc6ijj; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = (xyz_[jat] - xyz_[iat] + tau).norm2(); - if (r2 > para_.rthr2()) + if (r2 > para_.rthr2()) { continue; +} r = std::sqrt(r2); r6 = std::pow(r2, 3); @@ -794,6 +853,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc6i[jat] += dc6_rest * dc6ijj; dc6_rest_sum[linij] += dc6_rest; } // end tau +} +} } // end jat } // end iat } // end d3_0 @@ -812,8 +873,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m r0 = para_.rs6() * std::sqrt(3.0 * r42) + para_.rs18(); rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[iat]]; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] @@ -843,6 +904,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc6_rest_sum[linii] += dc6_rest; } } // end tau +} +} for (int jat = 0; jat != iat; jat++) { get_dc6_dcnij(para_.mxc()[iz_[iat]], para_.mxc()[iz_[jat]], cn[iat], cn[jat], @@ -855,15 +918,16 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[jat]]; dc6ij[jat][iat] = dc6iji; dc6ij[iat][jat] = dc6ijj; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = (xyz_[jat] - xyz_[iat] + tau).norm2(); - if (r2 > para_.rthr2()) + if (r2 > para_.rthr2()) { continue; +} r = std::sqrt(r2); r4 = r2 * r2; @@ -884,6 +948,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc6i[jat] += dc6_rest * dc6ijj; dc6_rest_sum[linij] += dc6_rest; } // end tau +} +} } // end jat } // end iat } // end d3_bj @@ -929,24 +995,27 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = (ijvec + jtau).norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / para_.r0ab()[iz_[jat]][iz_[iat]]; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / para_.r0ab()[iz_[kat]][iz_[iat]]; rr0jk = std::sqrt(rjk2) / para_.r0ab()[iz_[kat]][iz_[jat]]; @@ -1004,6 +1073,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc9 = (dc6ij[iat][kat] / c6ik + dc6ij[jat][kat] / c6jk) * c9 * 0.5; dc6i[kat] += dc6_rest * dc9; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux @@ -1040,29 +1111,33 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m { repmin[2] = std::max(-rep_cn_[2], jtauz - rep_cn_[2]); repmax[2] = std::min(rep_cn_[2], jtauz + rep_cn_[2]); - if (jtaux == 0 && jtauy == 0 && jtauz == 0) + if (jtaux == 0 && jtauy == 0 && jtauz == 0) { continue; +} jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = jtau.norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / para_.r0ab()[iz_[jat]][iz_[iat]]; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / para_.r0ab()[iz_[kat]][iz_[iat]]; rr0jk = std::sqrt(rjk2) / para_.r0ab()[iz_[kat]][iz_[jat]]; @@ -1120,13 +1195,15 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc9 = (dc6ij[iat][kat] / c6ik + dc6ij[jat][kat] / c6jk) * c9 * 0.5; dc6i[kat] += dc6_rest * dc9; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux } // end kat } // end iat - for (int iat = 1; iat != ucell_.nat; iat++) + for (int iat = 1; iat != ucell_.nat; iat++) { for (int jat = 0; jat != iat; jat++) { int kat = jat; @@ -1157,27 +1234,31 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = (ijvec + jtau).norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / para_.r0ab()[iz_[jat]][iz_[iat]]; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { - if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) + if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) { continue; +} ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = (ikvec + ktau).norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / para_.r0ab()[iz_[kat]][iz_[iat]]; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / para_.r0ab()[iz_[kat]][iz_[jat]]; geomean2 = rij2 * rjk2 * rik2; @@ -1234,10 +1315,13 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc9 = (dc6ij[iat][kat] / c6ik + dc6ij[jat][kat] / c6jk) * c9 * 0.5; dc6i[kat] += dc6_rest * dc9; } // end ktau +} +} } // end jtauz } // end jtauy } // end jtaux } // end jat +} // end iat for (int iat = 0; iat != ucell_.nat; iat++) @@ -1267,33 +1351,39 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m { repmin[2] = std::max(-rep_cn_[2], jtauz - rep_cn_[2]); repmax[2] = std::min(rep_cn_[2], jtauz + rep_cn_[2]); - if (jtaux == 0 && jtauy == 0 && jtauz == 0) + if (jtaux == 0 && jtauy == 0 && jtauz == 0) { continue; +} jtau = static_cast(jtaux) * lat_[0] + static_cast(jtauy) * lat_[1] + static_cast(jtauz) * lat_[2]; rij2 = jtau.norm2(); - if (rij2 > para_.cn_thr2()) + if (rij2 > para_.cn_thr2()) { continue; +} rr0ij = std::sqrt(rij2) / para_.r0ab()[iz_[jat]][iz_[iat]]; - for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) - for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) + for (int ktaux = repmin[0]; ktaux <= repmax[0]; ktaux++) { + for (int ktauy = repmin[1]; ktauy <= repmax[1]; ktauy++) { for (int ktauz = repmin[2]; ktauz <= repmax[2]; ktauz++) { - if (ktaux == 0 && ktauy == 0 && ktauz == 0) + if (ktaux == 0 && ktauy == 0 && ktauz == 0) { continue; - if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) +} + if (jtaux == ktaux && jtauy == ktauy && jtauz == ktauz) { continue; +} ktau = static_cast(ktaux) * lat_[0] + static_cast(ktauy) * lat_[1] + static_cast(ktauz) * lat_[2]; rik2 = ktau.norm2(); - if (rik2 > para_.cn_thr2()) + if (rik2 > para_.cn_thr2()) { continue; +} rr0ik = std::sqrt(rik2) / para_.r0ab()[iz_[kat]][iz_[iat]]; rjk2 = (jkvec + ktau - jtau).norm2(); - if (rjk2 > para_.cn_thr2()) + if (rjk2 > para_.cn_thr2()) { continue; +} rr0jk = std::sqrt(rjk2) / para_.r0ab()[iz_[kat]][iz_[jat]]; geomean2 = rij2 * rjk2 * rik2; @@ -1350,6 +1440,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m dc9 = (dc6ij[iat][kat] / c6ik + dc6ij[jat][kat] / c6jk) * c9 * 0.5; dc6i[kat] += dc6_rest * dc9; } // end ktau +} +} } // end jtauz } // end jtauy } // jtaux @@ -1359,29 +1451,31 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m // dE/dr_ij * dr_ij/dxyz_i double expterm, dcnn, x1; ModuleBase::Vector3 rij, vec3; - for (int iat = 1; iat != ucell_.nat; iat++) + for (int iat = 1; iat != ucell_.nat; iat++) { for (int jat = 0; jat != iat; jat++) { linij = lin(iat, jat); rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[jat]]; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; rij = xyz_[jat] - xyz_[iat] + tau; r2 = rij.norm2(); - if (r2 > para_.rthr2() || r2 < 0.5) + if (r2 > para_.rthr2() || r2 < 0.5) { continue; +} r = std::sqrt(r2); if (r2 < para_.cn_thr2()) { expterm = exp(-para_.k1() * (rcovij / r - 1.0)); dcnn = -para_.k1() * rcovij * expterm / (r2 * (expterm + 1.0) * (expterm + 1.0)); } - else + else { dcnn = 0.0; +} x1 = drij[linij][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] + dcnn * (dc6i[iat] + dc6i[jat]); vec3 = x1 * rij / r; @@ -1390,23 +1484,28 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m std::vector vec = {vec3.x, vec3.y, vec3.z}; std::vector rij_vec = {rij.x, rij.y, rij.z}; - for (size_t i = 0; i != 3; i++) + for (size_t i = 0; i != 3; i++) { for (size_t j = 0; j != 3; j++) { smearing_sigma(i, j) += vec[j] * rij_vec[i]; } +} } // end tau +} +} } // end iat, jat +} for (int iat = 0; iat != ucell_.nat; iat++) { linii = lin(iat, iat); rcovij = para_.rcov()[iz_[iat]] + para_.rcov()[iz_[iat]]; - for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) - for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) + for (int taux = -rep_vdw_[0]; taux <= rep_vdw_[0]; taux++) { + for (int tauy = -rep_vdw_[1]; tauy <= rep_vdw_[1]; tauy++) { for (int tauz = -rep_vdw_[2]; tauz <= rep_vdw_[2]; tauz++) { - if (taux == 0 && tauy == 0 && tauz == 0) + if (taux == 0 && tauy == 0 && tauz == 0) { continue; +} tau = static_cast(taux) * lat_[0] + static_cast(tauy) * lat_[1] + static_cast(tauz) * lat_[2]; r2 = tau.norm2(); @@ -1416,19 +1515,23 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m expterm = exp(-para_.k1() * (rcovij / r - 1.0)); dcnn = -para_.k1() * rcovij * expterm / (r2 * (expterm + 1.0) * (expterm + 1.0)); } - else + else { dcnn = 0.0; +} x1 = drij[linii][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] + dcnn * dc6i[iat]; vec3 = x1 * tau / r; std::vector vec = {vec3.x, vec3.y, vec3.z}; std::vector tau_vec = {tau.x, tau.y, tau.z}; - for (size_t i = 0; i != 3; i++) + for (size_t i = 0; i != 3; i++) { for (size_t j = 0; j != 3; j++) { smearing_sigma(i, j) += vec[j] * tau_vec[i]; } +} } // end tau +} +} } // end iat } diff --git a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp index e0336c3828..4e1bed1ae0 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp @@ -5,14 +5,15 @@ //========================================================== #include "vdwd3_parameters.h" - #include "module_base/constants.h" - +#include +#include "dftd3_xc_param.h" namespace vdw { -void Vdwd3Parameters::initial_parameters(const Input_para &input) +void Vdwd3Parameters::initial_parameters(const Input_para &input, std::ofstream* plog) { + // initialize the dftd3 parameters mxc_.resize(max_elem_, 1); r0ab_.resize(max_elem_, std::vector(max_elem_, 0.0)); @@ -22,11 +23,11 @@ void Vdwd3Parameters::initial_parameters(const Input_para &input) std::vector>>( 5, std::vector>(max_elem_, std::vector(max_elem_, 0.0))))); - - s6_ = std::stod(input.vdw_s6); - s18_ = std::stod(input.vdw_s8); - rs6_ = std::stod(input.vdw_a1); - rs18_ = std::stod(input.vdw_a2); + + DFTD3::dftd3_params(input.dft_functional, input.vdw_method, + input.vdw_s6, input.vdw_s8, input.vdw_a1, input.vdw_a2, + s6_, s18_, rs6_, rs18_, /* rs6: a1, rs18: a2 */ + plog); abc_ = input.vdw_abc; version_ = input.vdw_method; model_ = input.vdw_cutoff_type; diff --git a/source/module_hamilt_general/module_vdw/vdwd3_parameters.h b/source/module_hamilt_general/module_vdw/vdwd3_parameters.h index c5f005a7b7..95d1a03143 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_parameters.h +++ b/source/module_hamilt_general/module_vdw/vdwd3_parameters.h @@ -21,7 +21,14 @@ class Vdwd3Parameters : public VdwParameters ~Vdwd3Parameters() = default; - void initial_parameters(const Input_para &input); + /** + * @brief initialize the parameter by either input (from user setting) or autoset by dft XC + * + * @param input Parameter instance + * @param plog optional, for logging the parameter setting process + */ + void initial_parameters(const Input_para &input, + std::ofstream* plog = nullptr); // for logging the parameter autoset inline const std::string &version() const { return version_; } diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.h b/source/module_hamilt_general/module_xc/xc_functional_libxc.h index 75b2693d0e..26044c1b90 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -11,7 +11,7 @@ #include #include - class Charge; +class Charge; namespace XC_Functional_Libxc { @@ -19,158 +19,158 @@ namespace XC_Functional_Libxc // xc_functional_libxc.cpp //------------------- - // sets functional type, which allows combination of LIBXC keyword connected by "+" - // for example, "XC_LDA_X+XC_LDA_C_PZ" - extern std::pair> set_xc_type_libxc(const std::string xc_func_in); + // sets functional type, which allows combination of LIBXC keyword connected by "+" + // for example, "XC_LDA_X+XC_LDA_C_PZ" + extern std::pair> set_xc_type_libxc(const std::string xc_func_in); - // converts func_id into corresponding xc_func_type vector - extern std::vector init_func(const std::vector &func_id, const int xc_polarized); + // converts func_id into corresponding xc_func_type vector + extern std::vector init_func(const std::vector &func_id, const int xc_polarized); - extern void finish_func(std::vector &funcs); + extern void finish_func(std::vector &funcs); //------------------- // xc_functional_libxc_vxc.cpp //------------------- - extern std::tuple v_xc_libxc( - const std::vector &func_id, - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr); // charge density + extern std::tuple v_xc_libxc( + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge* const chr); // charge density - // for mGGA functional - extern std::tuple v_xc_meta( - const std::vector &func_id, - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr); + // for mGGA functional + extern std::tuple v_xc_meta( + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge* const chr); //------------------- // xc_functional_libxc_tools.cpp //------------------- - // converting rho (abacus=>libxc) - extern std::vector convert_rho( - const int nspin, - const std::size_t nrxx, - const Charge* const chr); - - // converting rho (abacus=>libxc) - extern std::tuple, std::vector> convert_rho_amag_nspin4( - const int nspin, - const std::size_t nrxx, - const Charge* const chr); - - // calculating grho - extern std::vector>> cal_gdr( - const int nspin, - const std::size_t nrxx, - const std::vector &rho, - const double tpiba, - const Charge* const chr); - - // converting grho (abacus=>libxc) - extern std::vector convert_sigma( - const std::vector>> &gdr); - - // sgn for threshold mask - extern std::vector cal_sgn( - const double rho_threshold, - const double grho_threshold, - const xc_func_type &func, - const int nspin, - const std::size_t nrxx, - const std::vector &rho, - const std::vector &sigma); - - // converting etxc from exc (libxc=>abacus) - extern double convert_etxc( - const int nspin, - const std::size_t nrxx, - const std::vector &sgn, - const std::vector &rho, - std::vector exc); - - // converting vtxc and v from vrho and vsigma (libxc=>abacus) - extern std::pair convert_vtxc_v( - const xc_func_type &func, - const int nspin, - const std::size_t nrxx, - const std::vector &sgn, - const std::vector &rho, - const std::vector>> &gdr, - const std::vector &vrho, - const std::vector &vsigma, - const double tpiba, - const Charge* const chr); - - // dh for gga v - extern std::vector> cal_dh( - const int nspin, - const std::size_t nrxx, - const std::vector &sgn, - const std::vector>> &gdr, - const std::vector &vsigma, - const double tpiba, - const Charge* const chr); - - // convert v for NSPIN=4 - extern ModuleBase::matrix convert_v_nspin4( - const std::size_t nrxx, - const Charge* const chr, - const std::vector &amag, - const ModuleBase::matrix &v); + // converting rho (abacus=>libxc) + extern std::vector convert_rho( + const int nspin, + const std::size_t nrxx, + const Charge* const chr); + + // converting rho (abacus=>libxc) + extern std::tuple, std::vector> convert_rho_amag_nspin4( + const int nspin, + const std::size_t nrxx, + const Charge* const chr); + + // calculating grho + extern std::vector>> cal_gdr( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba, + const Charge* const chr); + + // converting grho (abacus=>libxc) + extern std::vector convert_sigma( + const std::vector>> &gdr); + + // sgn for threshold mask + extern std::vector cal_sgn( + const double rho_threshold, + const double grho_threshold, + const xc_func_type &func, + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &sigma); + + // converting etxc from exc (libxc=>abacus) + extern double convert_etxc( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + std::vector exc); + + // converting vtxc and v from vrho and vsigma (libxc=>abacus) + extern std::pair convert_vtxc_v( + const xc_func_type &func, + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + const std::vector>> &gdr, + const std::vector &vrho, + const std::vector &vsigma, + const double tpiba, + const Charge* const chr); + + // dh for gga v + extern std::vector> cal_dh( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const double tpiba, + const Charge* const chr); + + // convert v for NSPIN=4 + extern ModuleBase::matrix convert_v_nspin4( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &amag, + const ModuleBase::matrix &v); //------------------- // xc_functional_libxc_wrapper_xc.cpp //------------------- - extern void xc_spin_libxc( - const std::vector &func_id, - const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw); + extern void xc_spin_libxc( + const std::vector &func_id, + const double &rhoup, const double &rhodw, + double &exc, double &vxcup, double &vxcdw); //------------------- // xc_functional_libxc_wrapper_gcxc.cpp //------------------- - // the entire GGA functional, for nspin=1 case - extern void gcxc_libxc( - const std::vector &func_id, - const double &rho, const double &grho, - double &sxc, double &v1xc, double &v2xc); + // the entire GGA functional, for nspin=1 case + extern void gcxc_libxc( + const std::vector &func_id, + const double &rho, const double &grho, + double &sxc, double &v1xc, double &v2xc); - // the entire GGA functional, for nspin=2 case - extern void gcxc_spin_libxc( - const std::vector &func_id, - const double rhoup, const double rhodw, - const ModuleBase::Vector3 gdr1, const ModuleBase::Vector3 gdr2, - double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); + // the entire GGA functional, for nspin=2 case + extern void gcxc_spin_libxc( + const std::vector &func_id, + const double rhoup, const double rhodw, + const ModuleBase::Vector3 gdr1, const ModuleBase::Vector3 gdr2, + double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); //------------------- // xc_functional_libxc_wrapper_tauxc.cpp //------------------- - // wrapper for the mGGA functionals - extern void tau_xc( - const std::vector &func_id, - const double &rho, const double &grho, const double &atau, double &sxc, - double &v1xc, double &v2xc, double &v3xc); - - extern void tau_xc_spin( - const std::vector &func_id, - double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, - double tauup, double taudw, - double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, - double &v3xcup, double &v3xcdw); + // wrapper for the mGGA functionals + extern void tau_xc( + const std::vector &func_id, + const double &rho, const double &grho, const double &atau, double &sxc, + double &v1xc, double &v2xc, double &v3xc); + + extern void tau_xc_spin( + const std::vector &func_id, + double rhoup, double rhodw, + ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double tauup, double taudw, + double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, + double &v3xcup, double &v3xcdw); } // namespace XC_Functional_Libxc #endif // USE_LIBXC diff --git a/source/module_io/read_input_item_model.cpp b/source/module_io/read_input_item_model.cpp index 964e0f2c52..30549e5d6f 100644 --- a/source/module_io/read_input_item_model.cpp +++ b/source/module_io/read_input_item_model.cpp @@ -147,10 +147,10 @@ void ReadInput::item_model() { para.input.vdw_s6 = "0.75"; } - else if (para.input.vdw_method == "d3_0" || para.input.vdw_method == "d3_bj") - { - para.input.vdw_s6 = "1.0"; - } + // else if (para.input.vdw_method == "d3_0" || para.input.vdw_method == "d3_bj") + // { + // para.input.vdw_s6 = "1.0"; + // } } }; read_sync_string(input.vdw_s6); @@ -160,17 +160,17 @@ void ReadInput::item_model() Input_Item item("vdw_s8"); item.annotation = "scale parameter of d3_0/d3_bj"; item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.vdw_s8 == "default") - { - if (para.input.vdw_method == "d3_0") - { - para.input.vdw_s8 = "0.722"; - } - else if (para.input.vdw_method == "d3_bj") - { - para.input.vdw_s8 = "0.7875"; - } - } + // if (para.input.vdw_s8 == "default") + // { + // if (para.input.vdw_method == "d3_0") + // { + // para.input.vdw_s8 = "0.722"; + // } + // else if (para.input.vdw_method == "d3_bj") + // { + // para.input.vdw_s8 = "0.7875"; + // } + // } }; read_sync_string(input.vdw_s8); this->add_item(item); @@ -179,17 +179,17 @@ void ReadInput::item_model() Input_Item item("vdw_a1"); item.annotation = "damping parameter of d3_0/d3_bj"; item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.vdw_a1 == "default") - { - if (para.input.vdw_method == "d3_0") - { - para.input.vdw_a1 = "1.217"; - } - else if (para.input.vdw_method == "d3_bj") - { - para.input.vdw_a1 = "0.4289"; - } - } + // if (para.input.vdw_a1 == "default") + // { + // if (para.input.vdw_method == "d3_0") + // { + // para.input.vdw_a1 = "1.217"; + // } + // else if (para.input.vdw_method == "d3_bj") + // { + // para.input.vdw_a1 = "0.4289"; + // } + // } }; read_sync_string(input.vdw_a1); this->add_item(item); @@ -198,17 +198,17 @@ void ReadInput::item_model() Input_Item item("vdw_a2"); item.annotation = "damping parameter of d3_bj"; item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.vdw_a2 == "default") - { - if (para.input.vdw_method == "d3_0") - { - para.input.vdw_a2 = "1.0"; - } - else if (para.input.vdw_method == "d3_bj") - { - para.input.vdw_a2 = "4.4407"; - } - } + // if (para.input.vdw_a2 == "default") + // { + // if (para.input.vdw_method == "d3_0") + // { + // para.input.vdw_a2 = "1.0"; + // } + // else if (para.input.vdw_method == "d3_bj") + // { + // para.input.vdw_a2 = "4.4407"; + // } + // } }; read_sync_string(input.vdw_a2); this->add_item(item); diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp index 013462fe3f..4a5a21218d 100644 --- a/source/module_io/test_serial/read_input_item_test.cpp +++ b/source/module_io/test_serial/read_input_item_test.cpp @@ -1049,46 +1049,54 @@ TEST_F(InputTest, Item_test2) it->second.reset_value(it->second, param); EXPECT_EQ(param.input.vdw_s6, "0.75"); + // dftd3 parameter will not get its value here param.input.vdw_s6 = "default"; param.input.vdw_method = "d3_0"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_s6, "1.0"); + // EXPECT_EQ(param.input.vdw_s6, "1.0"); + EXPECT_EQ(param.input.vdw_s6, "default"); } { // vdw_s8 auto it = find_label("vdw_s8", readinput.input_lists); param.input.vdw_s8 = "default"; param.input.vdw_method = "d3_0"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_s8, "0.722"); + // EXPECT_EQ(param.input.vdw_s8, "0.722"); + EXPECT_EQ(param.input.vdw_s8, "default"); param.input.vdw_s8 = "default"; param.input.vdw_method = "d3_bj"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_s8, "0.7875"); + // EXPECT_EQ(param.input.vdw_s8, "0.7875"); + EXPECT_EQ(param.input.vdw_s8, "default"); } { // vdw_a1 auto it = find_label("vdw_a1", readinput.input_lists); param.input.vdw_a1 = "default"; param.input.vdw_method = "d3_0"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_a1, "1.217"); + // EXPECT_EQ(param.input.vdw_a1, "1.217"); + EXPECT_EQ(param.input.vdw_a1, "default"); param.input.vdw_a1 = "default"; param.input.vdw_method = "d3_bj"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_a1, "0.4289"); + // EXPECT_EQ(param.input.vdw_a1, "0.4289"); + EXPECT_EQ(param.input.vdw_a1, "default"); } { // vdw_a2 auto it = find_label("vdw_a2", readinput.input_lists); param.input.vdw_a2 = "default"; param.input.vdw_method = "d3_0"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_a2, "1.0"); + // EXPECT_EQ(param.input.vdw_a2, "1.0"); + EXPECT_EQ(param.input.vdw_a2, "default"); param.input.vdw_a2 = "default"; param.input.vdw_method = "d3_bj"; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.vdw_a2, "4.4407"); + // EXPECT_EQ(param.input.vdw_a2, "4.4407"); + EXPECT_EQ(param.input.vdw_a2, "default"); } { // vdw_c6_unit auto it = find_label("vdw_c6_unit", readinput.input_lists); diff --git a/tests/integrate/150_PW_15_CR_VDW3/INPUT b/tests/integrate/150_PW_15_CR_VDW3/INPUT index 75926221a8..fd89158b41 100644 --- a/tests/integrate/150_PW_15_CR_VDW3/INPUT +++ b/tests/integrate/150_PW_15_CR_VDW3/INPUT @@ -20,3 +20,4 @@ cal_force 1 cal_stress 1 pseudo_dir ../../PP_ORB orbital_dir ../../PP_ORB +dft_functional pbe diff --git a/tests/integrate/250_NO_KP_CR_VDW3/INPUT b/tests/integrate/250_NO_KP_CR_VDW3/INPUT index 2ee242a2dc..62b007f191 100644 --- a/tests/integrate/250_NO_KP_CR_VDW3/INPUT +++ b/tests/integrate/250_NO_KP_CR_VDW3/INPUT @@ -20,3 +20,4 @@ cal_force 1 cal_stress 1 pseudo_dir ../../PP_ORB orbital_dir ../../PP_ORB +dft_functional pbe diff --git a/tests/integrate/250_NO_KP_CR_VDW3ABC/INPUT b/tests/integrate/250_NO_KP_CR_VDW3ABC/INPUT index 7a1ecc4762..83222fd70c 100644 --- a/tests/integrate/250_NO_KP_CR_VDW3ABC/INPUT +++ b/tests/integrate/250_NO_KP_CR_VDW3ABC/INPUT @@ -17,3 +17,4 @@ cal_force 1 cal_stress 1 pseudo_dir ../../PP_ORB orbital_dir ../../PP_ORB +dft_functional pbe diff --git a/tests/integrate/250_NO_KP_CR_VDW3BJ/INPUT b/tests/integrate/250_NO_KP_CR_VDW3BJ/INPUT index e3eda14906..a0753cdbc6 100644 --- a/tests/integrate/250_NO_KP_CR_VDW3BJ/INPUT +++ b/tests/integrate/250_NO_KP_CR_VDW3BJ/INPUT @@ -20,3 +20,4 @@ cal_force 1 cal_stress 1 pseudo_dir ../../PP_ORB orbital_dir ../../PP_ORB +dft_functional pbe