get_realspace_cutoff Function

public function get_realspace_cutoff(mol, disp, accuracy) result(cutoff)

Realspace cutoffs reproducing the two-body dispersion energy within a requested accuracy per atom.

Beyond the damping region the pair sum is dominated by the bare C6 tail. Integrating it over the d periodic dimensions gives the truncation error

dE/atom = rho * * S(d) / (2 * (6 - d) * R**(6 - d))

with the number density rho of the periodic cell, the mean pair coefficient and the surface S(d) of the d-dimensional unit sphere. Inverting this relation yields the cutoff. The estimate assumes a scaling factor s6 of one and neglects the faster decaying C8 tail, both covered by a safety margin.

Only the two-body cutoff is derived, the coordination number cutoff is not an accuracy parameter because the counting function has a non-vanishing limit, and the three-body cutoff is left at its default.

Arguments

Type IntentOptional Attributes Name
class(structure_type), intent(in) :: mol

Molecular structure data
class(d3_model), intent(in) :: disp

Dispersion model
real(kind=wp), intent(in) :: accuracy

Requested accuracy of the two-body dispersion energy in Hartree per atom

Return Value type(realspace_cutoff)

Realspace cutoffs


Source Code

function get_realspace_cutoff(mol, disp, accuracy) result(cutoff)

   !> Molecular structure data
   class(structure_type), intent(in) :: mol

   !> Dispersion model
   class(d3_model), intent(in) :: disp

   !> Requested accuracy of the two-body dispersion energy in Hartree per atom
   real(wp), intent(in) :: accuracy

   !> Realspace cutoffs
   type(realspace_cutoff) :: cutoff

   integer :: mref, ndim, iat, izp, isp, jsp, iref, jref, jat
   real(wp) :: measure, c6mean, amp, r2
   real(wp), allocatable :: cn(:), gwvec(:, :), wsum(:, :), lattr(:, :)

   cutoff = realspace_cutoff()
   if (accuracy <= 0.0_wp) return

   mref = maxval(disp%ref)
   allocate(cn(mol%nat), gwvec(mref, mol%nat))
   call get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
   call get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, cn)
   call disp%weight_references(mol, cn, gwvec)

   ! contract the reference coefficients with the accumulated weights per species
   allocate(wsum(mref, mol%nid), source=0.0_wp)
   do iat = 1, mol%nat
      if (disp%ghost(iat)) cycle
      izp = mol%id(iat)
      wsum(:, izp) = wsum(:, izp) + gwvec(:, iat)
   end do
   c6mean = 0.0_wp
   do isp = 1, mol%nid
      do jsp = 1, mol%nid
         do iref = 1, disp%ref(isp)
            do jref = 1, disp%ref(jsp)
               c6mean = c6mean + disp%c6(iref, jref, isp, jsp) &
                  & * wsum(iref, isp) * wsum(jref, jsp)
            end do
         end do
      end do
   end do
   c6mean = c6mean / mol%nat**2

   ndim = count(mol%periodic)
   measure = get_periodic_measure(mol, ndim)

   select case(ndim)
   case(3)
      amp = 2.0_wp*pi/3.0_wp * mol%nat/measure * c6mean
      cutoff%disp2 = (cutoff_safety * amp / accuracy)**(1.0_wp/3.0_wp)
   case(2)
      amp = 0.25_wp*pi * mol%nat/measure * c6mean
      cutoff%disp2 = (cutoff_safety * amp / accuracy)**(0.25_wp)
   case(1)
      amp = 0.2_wp * mol%nat/measure * c6mean
      cutoff%disp2 = (cutoff_safety * amp / accuracy)**(0.2_wp)
   case default
      ! a finite system is summed exactly once all pairs are covered
      r2 = 0.0_wp
      do iat = 1, mol%nat
         do jat = 1, iat - 1
            r2 = max(r2, sum((mol%xyz(:, iat) - mol%xyz(:, jat))**2))
         end do
      end do
      cutoff%disp2 = sqrt(r2)
   end select

   cutoff%disp2 = max(cutoff%disp2, cutoff_minimum)
   cutoff%disp3 = min(cutoff%disp3, cutoff%disp2)

end function get_realspace_cutoff