Add PySCF B3LYP variant names and DM21 functionals

python/dftd3/pyscf.py

@@ -116,7 +116,7 @@ def __init__(self, mol, xc="hf", version="d3bj", atm=False):
 
     def dump_flags(self, verbose=None):
         """
-        Show options used for the DFT-D4 dispersion correction.
+        Show options used for the DFT-D3 dispersion correction.
         """
         lib.logger.info(self, "** DFTD3 parameter **")
         lib.logger.info(self, "func %s", self.xc)
@@ -127,7 +127,7 @@ def dump_flags(self, verbose=None):
 
     def kernel(self) -> Tuple[float, np.ndarray]:
         """
-        Compute the DFT-D4 dispersion correction.
+        Compute the DFT-D3 dispersion correction.
 
         The dispersion model as well as the parameters are created locally and
         not part of the state of the instance.

python/dftd3/test_pyscf.py

@@ -56,6 +56,37 @@ def test_energy_r2scan_d3():
     assert d3.kernel()[0] == approx(-0.00578401192369041, abs=1.0e-7)
 
 
+@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
+@pytest.mark.parametrize("xc", ["b3lyp", "b3lypg", "b3lyp5", "b3lyp3"])
+def test_energy_b3lyp_d3(xc: str):
+
+    mol = gto.M(
+        atom="""
+             C   -0.755422531  -0.796459123  -1.023590391
+             C    0.634274834  -0.880017014  -1.075233285
+             C    1.406955202   0.199695367  -0.653144334
+             C    0.798863737   1.361204515  -0.180597909
+             C   -0.593166787   1.434312023  -0.133597923
+             C   -1.376239198   0.359205222  -0.553258516
+             I   -1.514344238   3.173268101   0.573601106
+             H    1.110906949  -1.778801728  -1.440619836
+             H    1.399172302   2.197767355   0.147412751
+             H    2.486417780   0.142466525  -0.689380574
+             H   -2.454252250   0.422581120  -0.512807958
+             H   -1.362353593  -1.630564523  -1.348743149
+             S   -3.112683203   6.289227834   1.226984439
+             H   -4.328789697   5.797771251   0.973373089
+             C   -2.689135032   6.703163830  -0.489062886
+             H   -1.684433029   7.115457372  -0.460265708
+             H   -2.683867206   5.816530502  -1.115183775
+             H   -3.365330613   7.451201412  -0.890098894
+             """
+    )
+
+    d3 = disp.DFTD3Dispersion(mol, xc=xc)
+    assert d3.kernel()[0] == approx(-0.034889411100056264, abs=1.0e-7)
+
+
 @pytest.mark.skipif(pyscf is None, reason="requires pyscf")
 def test_gradient_b97m_d3():
     mol = gto.M(
@@ -147,4 +178,4 @@ def test_gradient_hf():
         ]
     )
     grad = disp.energy(scf.RHF(mol)).run().nuc_grad_method()
-    assert grad.kernel() == approx(ref, abs=1.0e-8)
+    assert grad.kernel() == approx(ref, abs=1.0e-7)

src/dftd3/param.f90

@@ -40,7 +40,7 @@ module dftd3_param
    enum, bind(C)
       enumerator :: p_invalid, &
          & p_bp_df, p_blyp_df, p_revpbe_df, p_rpbe_df, p_b97d_df, p_pbe_df, &
-         & p_rpw86pbe_df, p_b3lyp_df, p_tpss_df, p_hf_df, p_tpss0_df, &
+         & p_rpw86pbe_df, p_b3lyp_df, p_b3lyp_g_df, p_tpss_df, p_hf_df, p_tpss0_df, &
          & p_pbe0_df, p_hse06_df, p_revpbe38_df, p_pw6b95_df, p_b2plyp_df, &
          & p_dsdblyp_df, p_dsdblypfc_df, p_bop_df, p_mpwlyp_df, p_olyp_df, &
          & p_pbesol_df, p_bpbe_df, p_opbe_df, p_ssb_df, p_revssb_df, p_otpss_df, &
@@ -60,7 +60,7 @@ module dftd3_param
          & p_mpw2plyp_df, p_m11_df, p_sogga11x_df, p_n12sx_df, p_mn12sx_df, &
          & p_ms2_df, p_ms2h_df, p_mpw1lyp_df, p_mpwkcis1k_df, p_pkzb_df, p_n12_df, &
          & p_m08hx_df, p_m11l_df, p_mn15l_df, p_pwp_df, p_r2scanh_df, p_r2scan0_df, &
-         & p_r2scan50_df, p_b973c_df
+         & p_r2scan50_df, p_b973c_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df
    end enum
 
 contains
@@ -90,7 +90,8 @@ function get_method_id(method) result(id)
    case("b1p", "b1p86"); id = p_b1p_df
    case("b2gpplyp"); id = p_b2gpplyp_df
    case("b2plyp"); id = p_b2plyp_df
-   case("b3lyp"); id = p_b3lyp_df
+   case("b3lyp", "b3lyp5"); id = p_b3lyp_df
+   case("b3lypg", "b3lyp3"); id = p_b3lyp_g_df
    case("b3lyp/631gd"); id = p_b3lyp_631gd_df
    case("b3p", "b3p86"); id = p_b3p_df
    case("b3pw91"); id = p_b3pw91_df
@@ -108,6 +109,10 @@ function get_method_id(method) result(id)
    case("bpbe"); id = p_bpbe_df
    case("camb3lyp"); id = p_camb3lyp_df
    case("dftb3"); id = p_dftb3_df
+   case("dm21"); id = p_dm21_df
+   case("dm21m"); id = p_dm21m_df
+   case("dm21mc"); id = p_dm21mc_df
+   case("dm21mu"); id = p_dm21mu_df
    case("dsdblyp"); id = p_dsdblyp_df
    case("dsdblypfc"); id = p_dsdblypfc_df
    case("hcth120"); id = p_hcth120_df
@@ -236,7 +241,7 @@ subroutine get_rational_damping(param, method, error, s9)
       param = d3_param(a1=0.4289_wp, s8=0.7875_wp, a2=4.4407_wp)
    case(p_rpw86pbe_df)
       param = d3_param(a1=0.4613_wp, s8=1.3845_wp, a2=4.5062_wp)
-   case(p_b3lyp_df)
+   case(p_b3lyp_df, p_b3lyp_g_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df)
       param = d3_param(a1=0.3981_wp, s8=1.9889_wp, a2=4.4211_wp)
    case(p_tpss_df)
       param = d3_param(a1=0.4535_wp, s8=1.9435_wp, a2=4.4752_wp)
@@ -457,7 +462,7 @@ subroutine get_zero_damping(param, method, error, s9)
       param = d3_param(rs6=0.872_wp, s8=0.514_wp)
    case(p_tpss_df)
       param = d3_param(rs6=1.166_wp, s8=1.105_wp)
-   case(p_b3lyp_df)
+   case(p_b3lyp_df, p_b3lyp_g_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df)
       param = d3_param(rs6=1.261_wp, s8=1.703_wp)
    case(p_pbe0_df)
       param = d3_param(rs6=1.287_wp, s8=0.928_wp)
@@ -634,7 +639,7 @@ subroutine get_mrational_damping(param, method, error, s9)
       return
    case(p_b2plyp_df)
       param = d3_param(a1=0.486434_wp, s8=0.672820_wp, a2=3.656466_wp, s6=0.640000_wp)
-   case(p_b3lyp_df)
+   case(p_b3lyp_df, p_b3lyp_g_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df)
       param = d3_param(a1=0.278672_wp, s8=1.466677_wp, a2=4.606311_wp)
    case(p_b97d_df)
       param = d3_param(a1=0.240184_wp, s8=1.206988_wp, a2=3.864426_wp)
@@ -677,7 +682,7 @@ subroutine get_mzero_damping(param, method, error, s9)
       return
    case(p_b2plyp_df)
       param = d3_param(rs6=1.313134_wp, s8=0.717543_wp, bet=0.016035_wp, s6=0.640000_wp)
-   case(p_b3lyp_df)
+   case(p_b3lyp_df, p_b3lyp_g_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df)
       param = d3_param(rs6=1.338153_wp, s8=1.532981_wp, bet=0.013988_wp)
    case(p_b97d_df)
       param = d3_param(rs6=1.151808_wp, s8=1.020078_wp, bet=0.035964_wp)
@@ -728,7 +733,7 @@ subroutine get_optimizedpower_damping(param, method, error, s9)
       param = d3_param(s6=1.0_wp, s8=0.12467_wp, a1=0.575_wp, a2=3.000_wp, bet=10.0_wp)
    case(p_blyp_df)
       param = d3_param(s6=1.0_wp, s8=1.31867_wp, a1=0.425_wp, a2=3.50_wp, bet=2.0_wp)
-   case(p_b3lyp_df)
+   case(p_b3lyp_df, p_b3lyp_g_df, p_dm21_df, p_dm21m_df, p_dm21mc_df, p_dm21mu_df)
       param = d3_param(s6=1.0_wp, s8=0.78311_wp, a1=0.300_wp, a2=4.25_wp, bet=4.0_wp)
    case(p_b97d_df)
       param = d3_param(s6=1.0_wp, s8=1.46861_wp, a1=0.600_wp, a2=2.50_wp, bet=0.0_wp)

test/unit/test_param.f90

@@ -110,7 +110,8 @@ subroutine test_d3bj_mb01(error)
       & "b1lyp", "mpw1pw", "mpw1kcis", "pbeh1pbe", "pbe1kcis", "x3lyp", "o3lyp", &
       & "b971", "b972", "b98", "hiss", "hse03", "revtpssh", "revtpss0", "tpss1kcis", &
       & "tauhcthhyb", "mn15", "lcwhpbe", "mpw2plyp", "m11", "sogga11x", "n12sx", "mn12sx", &
-      & "r2scanh", "r2scan0", "r2scan50", "b973c"]
+      & "r2scanh", "r2scan0", "r2scan50", "b973c", "b3lyp3", "b3lyp5", "b3lypg", &
+      & "dm21", "dm21m", "dm21mc", "dm21mu"]
    real(wp), parameter :: ref(*) = [&
       &-2.9551694676908012E-2_wp,-1.6638703086788331E-2_wp,-1.6725877716130381E-2_wp, &
       &-3.3014429592265318E-2_wp,-2.2051435219996540E-2_wp,-3.3481565825316001E-2_wp, &
@@ -142,7 +143,10 @@ subroutine test_d3bj_mb01(error)
       &-2.6879433016433214E-3_wp,-4.7068039876219803E-5_wp,-2.2830467677890126E-2_wp, &
       &-8.9046741008620024E-3_wp,-4.2882619531553470E-3_wp,-3.4182032740020921E-2_wp, &
       &-1.7851923991011053E-2_wp,-7.2749942776428288E-3_wp,-5.9043780290722733E-3_wp, &
-      &-6.5370091760150479E-3_wp,-7.3564457249637944E-3_wp,-4.2958524322550894E-2_wp]
+      &-6.5370091760150479E-3_wp,-7.3564457249637944E-3_wp,-4.2958524322550894E-2_wp, &
+      &-3.3014429592265318E-2_wp,-3.3014429592265318E-2_wp,-3.3014429592265318E-2_wp, &
+      &-3.3014429592265318E-2_wp,-3.3014429592265318E-2_wp,-3.3014429592265318E-2_wp, &
+      &-3.3014429592265318E-2_wp]
 
    call get_structure(mol, "MB16-43", "01")
    do ii = 1, size(func)