Update printout to include damping parameters

- include r⁴/r² expectation values for actinides / superheavies
- add rational damping parameters for r²SCAN hybrids
- allow changing of realspace cutoff in the API

app/driver.f90

@@ -239,7 +239,7 @@ subroutine run_driver(config, error)
       if (config%json) then
          open(file=config%json_output, newunit=unit)
          call json_results(unit, "  ", energy=energy, gradient=gradient, sigma=sigma, &
-            & pairwise_energy2=pair_disp2, pairwise_energy3=pair_disp3)
+            & pairwise_energy2=pair_disp2, pairwise_energy3=pair_disp3, param=param)
          close(unit)
          if (config%verbosity > 0) then
             write(output_unit, '(a)') &

assets/parameters.toml

@@ -225,6 +225,15 @@ d3.bj = {a1=0.47023427, s8=1.08859014, a2=5.73408312}
 [parameter.r2scan]
 d3.bj = {a1=0.49484001, s8=0.78981345, a2=5.73083694}
 
+[parameter.r2scanh]
+d3.bj = {s8=1.1236, a1=0.4709, a2=5.9157}
+
+[parameter.r2scan0]
+d3.bj = {s8=1.1846, a1=0.4534, a2=5.8972}
+
+[parameter.r2scan50]
+d3.bj = {s8=1.3294, a1=0.4311, a2=5.9240}
+
 [parameter.wb97x]
 d3.bj = {a1=0.0000, s8=0.2641, a2=5.4959}
 d3.zero = {rs6=1.281, s8=1.0, rs8=1.094}

include/s-dftd3.h

@@ -103,6 +103,14 @@ SDFTD3_API_ENTRY dftd3_model SDFTD3_API_CALL
 dftd3_new_d3_model(dftd3_error /* error */,
                    dftd3_structure /* mol */) SDFTD3_API_SUFFIX__V_0_2;
 
+/// Set realspace cutoffs (quantities in Bohr)
+SDFTD3_API_ENTRY void SDFTD3_API_CALL
+dftd3_set_model_realspace_cutoff(dftd3_error /* error */,
+                                 dftd3_model /* model */,
+                                 double /* disp2 */,
+                                 double /* disp3 */,
+                                 double /* cn */) SDFTD3_API_SUFFIX__V_0_5;
+
 /// Delete dispersion model
 SDFTD3_API_ENTRY void SDFTD3_API_CALL
 dftd3_delete_model(dftd3_model* /* disp */) SDFTD3_API_SUFFIX__V_0_2;

python/dftd3/interface.py

@@ -245,6 +245,7 @@ class ZeroDampingParam(DampingParam):
     of the dispersion energy a repulsive contribution to the gradient can arise, which
     is considered artificial.\ :footcite:`grimme2011`
     """
+
     def __init__(self, **kwargs):
         DampingParam.__init__(self, **kwargs)
 
@@ -278,6 +279,7 @@ class ModifiedRationalDampingParam(DampingParam):
     from the library rather than the original one. Providing a full parameter set is
     functionally equivalent to using the `RationalDampingParam` constructor.
     """
+
     def __init__(self, **kwargs):
         DampingParam.__init__(self, **kwargs)
 
@@ -311,6 +313,7 @@ class ModifiedZeroDampingParam(DampingParam):
 
        This damping function is identical to zero damping for ``bet=0.0``.
     """
+
     def __init__(self, **kwargs):
         DampingParam.__init__(self, **kwargs)
 
@@ -345,6 +348,7 @@ class OptimizedPowerDampingParam(DampingParam):
     from the library rather than the original one. Providing the parameter `bet=0` is
     equivalent to using rational the `RationalDampingParam` constructor.
     """
+
     def __init__(self, **kwargs):
         DampingParam.__init__(self, **kwargs)
 
@@ -392,6 +396,11 @@ def __init__(
 
         self._disp = library.new_d3_model(self._mol)
 
+    def set_realspace_cutoff(self, disp2: float, disp3: float, cn: float):
+        """Set realspace cutoff for evaluation of interactions"""
+
+        library.set_model_realspace_cutoff(self._disp, disp2, disp3, cn)
+
     def get_dispersion(self, param: DampingParam, grad: bool) -> dict:
         """Perform actual evaluation of the dispersion correction"""
 

python/dftd3/library.py

@@ -105,6 +105,9 @@ def new_d3_model(mol):
     return ffi.gc(error_check(lib.dftd3_new_d3_model)(mol), _delete_model)
 
 
+set_model_realspace_cutoff = error_check(lib.dftd3_set_model_realspace_cutoff)
+
+
 def _delete_param(param):
     """Delete a DFT-D3 damping parameteter object"""
     ptr = ffi.new("dftd3_param *")

src/dftd3/api.f90

@@ -75,6 +75,8 @@ module dftd3_api
    type :: vp_model
       !> Actual payload
       type(d3_model) :: ptr
+      !> Additional real space cutoff
+      type(realspace_cutoff), allocatable :: cutoff
    end type vp_model
 
    !> Void pointer to damping parameters
@@ -304,6 +306,29 @@ function new_d3_model_api(verror, vmol) &
 end function new_d3_model_api
 
 
+subroutine set_model_realspace_cutoff(verror, vdisp, disp2, disp3, cn) &
+      & bind(C, name=namespace//"set_model_realspace_cutoff")
+   type(c_ptr), value :: verror
+   type(vp_error), pointer :: error
+   type(c_ptr), value :: vdisp
+   type(vp_model), pointer :: disp
+   real(c_double), value, intent(in) :: disp2
+   real(c_double), value, intent(in) :: disp3
+   real(c_double), value, intent(in) :: cn
+
+   if (.not.c_associated(verror)) return
+   call c_f_pointer(verror, error)
+
+   if (.not.c_associated(vdisp)) then
+      call fatal_error(error%ptr, "D3 dispersion model is missing")
+      return
+   end if
+   call c_f_pointer(vdisp, disp)
+
+   disp%cutoff = realspace_cutoff(disp2=disp2, disp3=disp3, cn=cn)
+end subroutine set_model_realspace_cutoff
+
+
 !> Delete dispersion model
 subroutine delete_model_api(vdisp) &
       & bind(C, name=namespace//"delete_model")
@@ -681,6 +706,7 @@ subroutine get_dispersion_api(verror, vmol, vdisp, vparam, &
    real(wp), allocatable :: gradient(:, :)
    real(c_double), intent(out), optional :: c_sigma(3, 3)
    real(wp), allocatable :: sigma(:, :)
+   type(realspace_cutoff) :: cutoff
 
    if (.not.c_associated(verror)) return
    call c_f_pointer(verror, error)
@@ -716,7 +742,11 @@ subroutine get_dispersion_api(verror, vmol, vdisp, vparam, &
       sigma = c_sigma(:3, :3)
    endif
 
-   call get_dispersion(mol%ptr, disp%ptr, param%ptr, realspace_cutoff(), &
+   cutoff = realspace_cutoff()
+   if (allocated(disp%cutoff)) then
+      cutoff = disp%cutoff
+   end if
+   call get_dispersion(mol%ptr, disp%ptr, param%ptr, cutoff, &
       & energy, gradient, sigma)
 
    if (present(c_gradient)) then
@@ -746,6 +776,7 @@ subroutine get_pairwise_dispersion_api(verror, vmol, vdisp, vparam, &
    real(wp), pointer :: pair_energy2(:, :)
    type(c_ptr), value, intent(in) :: c_pair_energy3
    real(wp), pointer :: pair_energy3(:, :)
+   type(realspace_cutoff) :: cutoff
 
    if (.not.c_associated(verror)) return
    call c_f_pointer(verror, error)
@@ -776,7 +807,11 @@ subroutine get_pairwise_dispersion_api(verror, vmol, vdisp, vparam, &
    call c_f_pointer(c_pair_energy2, pair_energy2, [mol%ptr%nat, mol%ptr%nat])
    call c_f_pointer(c_pair_energy3, pair_energy3, [mol%ptr%nat, mol%ptr%nat])
 
-   call get_pairwise_dispersion(mol%ptr, disp%ptr, param%ptr, realspace_cutoff(), &
+   cutoff = realspace_cutoff()
+   if (allocated(disp%cutoff)) then
+      cutoff = disp%cutoff
+   end if
+   call get_pairwise_dispersion(mol%ptr, disp%ptr, param%ptr, cutoff, &
       & pair_energy2, pair_energy3)
 
 end subroutine get_pairwise_dispersion_api

src/dftd3/data/r4r2.f90

@@ -40,6 +40,7 @@ module dftd3_data_r4r2
    !  not replaced but recalculated (PBE0/cc-pVQZ) were
    !   H: 8.0589 ->10.9359, Li:29.0974 ->39.7226, Be:14.8517 ->17.7460
    !  also new super heavies Cn,Nh,Fl,Lv,Og
+   !  Am-Rg calculated at 4c-PBE/Dyall-AE4Z (Dirac 2022)
    real(wp), parameter :: r4_over_r2(max_elem) = [  &
       &  8.0589_wp, 3.4698_wp, & ! H,He
       & 29.0974_wp,14.8517_wp,11.8799_wp, 7.8715_wp, 5.5588_wp, 4.7566_wp, 3.8025_wp, 3.1036_wp, & ! Li-Ne
@@ -59,12 +60,11 @@ module dftd3_data_r4r2
       &            10.2671_wp, 8.3549_wp, 7.8496_wp, 7.3278_wp, 7.4820_wp, & ! -Hg
       &                       13.5124_wp,11.6554_wp,10.0959_wp, 9.7340_wp, 8.8584_wp, 8.0125_wp, & ! Tl-Rn
       & 29.8135_wp,26.3157_wp, & ! Fr,Ra
-      &            19.1885_wp,15.8542_wp,16.1305_wp,15.6161_wp,15.1226_wp,16.1576_wp, 0.0000_wp, & ! Ac-Am
-      &             0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, & ! Cm-No
-      &             0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, & ! Lr-
-      &             0.0000_wp, 0.0000_wp, 0.0000_wp, 0.0000_wp, 5.4929_wp, & ! -Cn
+      &            19.1885_wp,15.8542_wp,16.1305_wp,15.6161_wp,15.1226_wp,16.1576_wp,14.6510_wp, & ! Ac-Am
+      &            14.7178_wp,13.9108_wp,13.5623_wp,13.2326_wp,12.9189_wp,12.6133_wp,12.3142_wp, & ! Cm-No
+      &            14.8326_wp,12.3771_wp,10.6378_wp, 9.3638_wp, 8.2297_wp, & ! Lr-
+      &             7.5667_wp, 6.9456_wp, 6.3946_wp, 5.9159_wp, 5.4929_wp, & ! -Cn
       &                        6.7286_wp, 6.5144_wp,10.9169_wp,10.3600_wp, 9.4723_wp, 8.6641_wp ] ! Nh-Og
-
    integer :: idum
    real(wp), parameter :: sqrt_z_r4_over_r2(max_elem) = &
       &  sqrt(0.5_wp*(r4_over_r2*[(sqrt(real(idum,wp)),idum=1,max_elem)]))

src/dftd3/data/vdwrad.f90

@@ -18,6 +18,7 @@ module dftd3_data_vdwrad
    use mctc_env, only : wp
    use mctc_io_convert, only : aatoau
    use mctc_io_symbols, only : to_number
+   use dftd3_data_covrad, only : get_covalent_rad
    implicit none
    private
 
@@ -594,14 +595,27 @@ elemental function get_vdw_rad_pair_num(num1, num2) result(rad)
    !> Van-der-Waals radius
    real(wp) :: rad
 
+   ! Magic number from ratio of van-der-Waals and covalent radii over element 1-94.
+   !
+   ! L. Trombach, S. Ehlert, S. Grimme, P Schwerdtfeger, J.-M. Mewes, PCCP, 2019.
+   ! DOI: 10.1039/c9cp02455g
+   real(wp), parameter :: cov2vdw = 1.45_wp
+
    if (num1 > 0 .and. num1 <= max_elem .and. num2 > 0 .and. num2 <= max_elem) then
       if (num1 > num2) then
          rad = vdwrad(num2 + num1*(num1-1)/2)
       else
          rad = vdwrad(num1 + num2*(num2-1)/2)
-      endif
+      end if
    else
-      rad = 0.0_wp
+      ! best estimate
+      if (num1 > 0 .and. num1 <= max_elem) then
+         rad = 0.5_wp * vdwrad(num1*(num1+1)/2) + cov2vdw * get_covalent_rad(num2)
+      else if (num2 > 0 .and. num2 <= max_elem) then
+         rad = 0.5_wp * vdwrad(num2*(num2+1)/2) + cov2vdw * get_covalent_rad(num1)
+      else
+         rad = cov2vdw * (get_covalent_rad(num1) + get_covalent_rad(num2))
+      end if
    end if
 
 end function get_vdw_rad_pair_num

src/dftd3/output.f90

@@ -528,7 +528,7 @@ end subroutine getline
 
 
 subroutine json_results(unit, indentation, energy, gradient, sigma, cn, c6, &
-      & pairwise_energy2, pairwise_energy3)
+      & pairwise_energy2, pairwise_energy3, param)
    integer, intent(in) :: unit
    character(len=*), intent(in), optional :: indentation
    real(wp), intent(in), optional :: energy
@@ -538,6 +538,7 @@ subroutine json_results(unit, indentation, energy, gradient, sigma, cn, c6, &
    real(wp), intent(in), optional :: c6(:, :)
    real(wp), intent(in), optional :: pairwise_energy2(:, :)
    real(wp), intent(in), optional :: pairwise_energy3(:, :)
+   class(damping_param), intent(in), optional :: param
    character(len=:), allocatable :: indent, version_string
    character(len=*), parameter :: jsonkey = "('""',a,'"":',1x)"
    real(wp), allocatable :: array(:)
@@ -599,11 +600,97 @@ subroutine json_results(unit, indentation, energy, gradient, sigma, cn, c6, &
       array = reshape(pairwise_energy3, [size(pairwise_energy3)])
       call write_json_array(unit, array, indent)
    end if
+   if (present(param)) then
+      write(unit, '(",")', advance='no')
+      if (allocated(indent)) write(unit, '(/,a)', advance='no') repeat(indent, 1)
+      write(unit, jsonkey, advance='no') 'damping parameters'
+      select type(param)
+      type is(rational_damping_param)
+         call write_json_param(unit, "rational", s6=param%s6, s8=param%s8, &
+            & s9=param%s9, a1=param%a1, a2=param%a2, alp=param%alp, indent=indent)
+      type is(zero_damping_param)
+         call write_json_param(unit, "zero", s6=param%s6, s8=param%s8, &
+            & s9=param%s9, rs6=param%rs6, rs8=param%rs8, alp=param%alp, indent=indent)
+      type is(mzero_damping_param)
+         call write_json_param(unit, "mzero", s6=param%s6, s8=param%s8, s9=param%s9, &
+            & rs6=param%rs6, rs8=param%rs8, alp=param%alp, bet=param%bet, indent=indent)
+      type is(optimizedpower_damping_param)
+         call write_json_param(unit, "optimizedpower", s6=param%s6, s8=param%s8, &
+            & s9=param%s9, a1=param%a1, a2=param%a2, alp=param%alp, bet=param%bet, &
+            & indent=indent)
+      class default
+         call write_json_param(unit, "unknown", indent=indent)
+      end select
+   end if
    if (allocated(indent)) write(unit, '(/)', advance='no')
    write(unit, '("}")')
 
 end subroutine json_results
 
+subroutine write_json_param(unit, damping, s6, s8, s9, a1, a2, rs6, rs8, alp, bet, indent)
+   integer, intent(in) :: unit
+   character(len=*), intent(in) :: damping
+   real(wp), intent(in), optional :: s6, s8, s9, a1, a2, rs6, rs8, alp, bet
+   character(len=:), allocatable, intent(in) :: indent
+   character(len=*), parameter :: jsonkeyval = "('""',a,'"":',1x,'""',a,'""')"
+
+   write(unit, '("{")', advance='no')
+
+   if (allocated(indent)) write(unit, '(/,a)', advance='no') repeat(indent, 2)
+   write(unit, jsonkeyval, advance='no') 'damping', damping
+
+   if (present(s6)) then
+      call write_json_keyval(unit, 's6', s6, indent)
+   end if
+
+   if (present(s8)) then
+      call write_json_keyval(unit, 's8', s8, indent)
+   end if
+
+   if (present(s9)) then
+      call write_json_keyval(unit, 's9', s9, indent)
+   end if
+
+   if (present(a1)) then
+      call write_json_keyval(unit, 'a1', a1, indent)
+   end if
+
+   if (present(a2)) then
+      call write_json_keyval(unit, 'a2', a2, indent)
+   end if
+
+   if (present(rs6)) then
+      call write_json_keyval(unit, 'rs6', rs6, indent)
+   end if
+
+   if (present(rs8)) then
+      call write_json_keyval(unit, 'rs8', rs8, indent)
+   end if
+
+   if (present(alp)) then
+      call write_json_keyval(unit, 'alp', alp, indent)
+   end if
+
+   if (present(bet)) then
+      call write_json_keyval(unit, 'bet', bet, indent)
+   end if
+
+   if (allocated(indent)) write(unit, '(/,a)', advance='no') repeat(indent, 1)
+   write(unit, '("}")', advance='no')
+end subroutine write_json_param
+
+subroutine write_json_keyval(unit, key, val, indent)
+   integer, intent(in) :: unit
+   character(len=*), intent(in) :: key
+   real(wp), intent(in) :: val
+   character(len=:), allocatable, intent(in) :: indent
+   character(len=*), parameter :: jsonkeyval = "('""',a,'"":',1x,es23.16)"
+
+   write(unit, '(",")', advance='no')
+   if (allocated(indent)) write(unit, '(/,a)', advance='no') repeat(indent, 2)
+   write(unit, jsonkeyval, advance='no') key, val
+
+end subroutine write_json_keyval
 
 subroutine write_json_array(unit, array, indent)
    integer, intent(in) :: unit

src/dftd3/param.f90

@@ -59,7 +59,8 @@ module dftd3_param
          & p_revtpss0_df, p_tpss1kcis_df, p_tauhcthhyb_df, p_mn15_df, p_lc_whpbe_df, &
          & 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_m08hx_df, p_m11l_df, p_mn15l_df, p_pwp_df, p_r2scanh_df, p_r2scan0_df, &
+         & p_r2scan50_df
    end enum
 
 contains
@@ -168,6 +169,9 @@ function get_method_id(method) result(id)
    case("pwgga"); id = p_pwgga_df
    case("pwpb95"); id = p_pwpb95_df
    case("r2scan"); id = p_r2scan_df
+   case("r2scanh"); id = p_r2scanh_df
+   case("r2scan0"); id = p_r2scan0_df
+   case("r2scan50"); id = p_r2scan50_df
    case("revpbe"); id = p_revpbe_df
    case("revpbe0"); id = p_revpbe0_df
    case("revpbe38"); id = p_revpbe38_df
@@ -325,6 +329,12 @@ subroutine get_rational_damping(param, method, error, s9)
       param = d3_param(a1=0.47023427_wp, s8=1.08859014_wp, a2=5.73408312_wp)
    case(p_r2scan_df)
       param = d3_param(a1=0.49484001_wp, s8=0.78981345_wp, a2=5.73083694_wp)
+   case(p_r2scanh_df)
+      param = d3_param(s8=1.1236_wp, a1=0.4709_wp, a2=5.9157_wp)
+   case(p_r2scan0_df)
+      param = d3_param(s8=1.1846_wp, a1=0.4534_wp, a2=5.8972_wp)
+   case(p_r2scan50_df)
+      param = d3_param(s8=1.3294_wp, a1=0.4311_wp, a2=5.9240_wp)
    case(p_wb97x_df)
       param = d3_param(a1=0.0000_wp, s8=0.2641_wp, a2=5.4959_wp)
    case(p_wb97m_df)

test/api/api-test.c

@@ -150,6 +150,9 @@ test (void) {
    if (dftd3_check_error(error)) {return 1;}
    dftd3_delete_param(&param);
 
+   dftd3_set_model_realspace_cutoff(error, disp, 50.0, 30.0, 25.0);
+   if (dftd3_check_error(error)) {return 1;}
+
    // DSD-BLYP-D3(BJ)-ATM
    param = dftd3_load_rational_damping(error, "dsdblyp", true);
    if (dftd3_check_error(error)) {return 1;}