model.f90 Source File

! This file is part of s-dftd3.
! SPDX-Identifier: LGPL-3.0-or-later
!
! s-dftd3 is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! s-dftd3 is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with s-dftd3.  If not, see <https://www.gnu.org/licenses/>.

module dftd3_model
   use, intrinsic :: ieee_arithmetic, only : ieee_is_nan
   use dftd3_citation, only : citation_type, get_citation, doi_fourier_d3
   use dftd3_data, only : get_r4r2_val, get_vdw_rad
   use dftd3_fourier_decomposition, only : d3_lowrank_c6, d3_lowrank_config, new_lowrank_c6
   use dftd3_partition, only : work_partition, owns_pair
   use dftd3_reference, only : number_of_references, reference_cn, get_c6, init_reference_c6
   use mctc_data, only : get_covalent_rad
   use mctc_env, only : wp
   use mctc_io, only : structure_type
   implicit none
   private

   public :: d3_model, new_d3_model
   public :: d3_lowrank_config


   !> Base D3 dispersion model to evaluate C6 coefficients
   type :: d3_model

      !> Covalent radii for coordination number
      real(wp), allocatable :: rcov(:)

      !> Van-der-Waals radii for damping function
      real(wp), allocatable :: rvdw(:, :)

      !> Expectation values for C8 extrapolation
      real(wp), allocatable :: r4r2(:)

      !> Weighting factor for CN interpolation
      real(wp), allocatable :: wf

      !> Number of reference systems
      integer, allocatable :: ref(:)

      !> Reference coordination numbers
      real(wp), allocatable :: cn(:, :)

      !> Reference C6 coefficients
      real(wp), allocatable :: c6(:, :, :, :)

      !> Masked
      logical, allocatable :: ghost(:)

      !> Separable representation of the reference C6 coefficients, only allocated
      !> if the model was set up with a low-rank approximation
      type(d3_lowrank_c6), allocatable :: lowrank

   contains

      !> Generate weights for all reference systems
      procedure :: weight_references

      !> Set up a separable low-rank representation of the C6 coefficients
      procedure :: set_lowrank

      !> Evaluate C6 coefficient
      procedure :: get_atomic_c6

   end type d3_model


   !> Default weighting factor for coordination number interpolation
   real(wp), parameter :: wf_default = 4.0_wp


contains


!> Create new dispersion model from molecular structure input
subroutine new_d3_model(self, mol, wf, ghost, lowrank, citation)

   !> Instance of the dispersion model
   type(d3_model), intent(out) :: self

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

   !> Weighting factor for coordination number interpolation
   real(wp), intent(in), optional :: wf

   !> Atom indices to ignore for dispersion calculation
   integer, intent(in), optional :: ghost(:)

   !> Setup for a separable low-rank representation of the C6 coefficients
   type(d3_lowrank_config), intent(in), optional :: lowrank

   !> Reference for the separable low-rank representation
   type(citation_type), intent(out), optional :: citation

   integer :: isp, izp, iref, jsp, jzp, jref
   integer :: mref

   call init_reference_c6

   if (present(wf)) then
      self%wf = wf
   else
      self%wf = wf_default
   end if

   allocate(self%ref(mol%nid))
   self%ref(:) = number_of_references(mol%num)
   mref = maxval(self%ref)

   allocate(self%rcov(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      self%rcov(isp) = get_covalent_rad(izp)
   end do

   allocate(self%r4r2(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      self%r4r2(isp) = get_r4r2_val(izp)
   end do

   allocate(self%rvdw(mol%nid, mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      do jsp = 1, isp
         jzp = mol%num(jsp)
         self%rvdw(jsp, isp) = get_vdw_rad(jzp, izp)
         self%rvdw(isp, jsp) = self%rvdw(jsp, isp)
      end do
   end do

   allocate(self%cn(mref, mol%nid), source=0.0_wp)
   do isp = 1, mol%nid
      izp = mol%num(isp)
      do iref = 1, self%ref(isp)
         self%cn(iref, isp) = reference_cn(iref, izp)
      end do
   end do

   allocate(self%c6(mref, mref, mol%nid, mol%nid), source=0.0_wp)
   do isp = 1, mol%nid
      izp = mol%num(isp)
      do jsp = 1, isp
         jzp = mol%num(jsp)
         do iref = 1, self%ref(isp)
            do jref = 1, self%ref(jsp)
               self%c6(jref, iref, jsp, isp) = get_c6(jref, iref, jzp, izp)
               self%c6(iref, jref, isp, jsp) = self%c6(jref, iref, jsp, isp)
            end do
         end do
      end do
   end do

   allocate(self%ghost(mol%nat))
   self%ghost(:) = .false.
   if (present(ghost)) then
      do iref = 1, size(ghost)
         if (ghost(iref) > 0 .and. ghost(iref) <= mol%nat) then
            self%ghost(ghost(iref)) = .true.
         end if
      end do
   end if

   if (present(lowrank)) then
      call self%set_lowrank(lowrank)
      if (present(citation)) citation = get_citation(doi_fourier_d3)
   end if

end subroutine new_d3_model


!> Set up a separable low-rank representation of the reference C6 coefficients
subroutine set_lowrank(self, config)

   !> Instance of the dispersion model
   class(d3_model), intent(inout) :: self

   !> Setup of the separable representation
   type(d3_lowrank_config), intent(in) :: config

   if (.not.allocated(self%lowrank)) allocate(self%lowrank)
   call new_lowrank_c6(self%lowrank, self%ref, self%c6, config)

end subroutine set_lowrank


!> Calculate the weights of the reference system and the derivatives w.r.t.
!> coordination number for later use.
subroutine weight_references(self, mol, cn, gwvec, gwdcn, gwd2cn)

   !> Instance of the dispersion model
   class(d3_model), intent(in) :: self

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

   !> Coordination number of every atom.
   real(wp), intent(in) :: cn(:)

   !> weighting for the atomic reference systems
   real(wp), intent(out) :: gwvec(:, :)

   !> derivative of the weighting function w.r.t. the coordination number
   real(wp), intent(out), optional :: gwdcn(:, :)

   !> second derivative of the weighting function w.r.t. the coordination number
   real(wp), intent(out), optional :: gwd2cn(:, :)

   integer :: iat, izp, iref, mref
   real(wp) :: norm, dnorm, d2norm, gw, expw, expd, expdd, gwk, dgwk, d2gwk
   real(wp) :: wf2, dcn
   real(wp), allocatable :: gwt(:)

   mref = maxval(self%ref)

   if (present(gwdcn)) then
      gwvec(:, :) = 0.0_wp
      gwdcn(:, :) = 0.0_wp
      if (present(gwd2cn)) gwd2cn(:, :) = 0.0_wp
      wf2 = 2 * self%wf

      !$omp parallel default(none) &
      !$omp shared(gwvec, gwdcn, gwd2cn, mol, self, cn, wf2, mref) &
      !$omp private(iat, izp, iref, norm, dnorm, d2norm, gw, expw, expd, expdd, &
      !$omp& gwk, dgwk, d2gwk, dcn, gwt)
      allocate(gwt(mref))
      !$omp do schedule(runtime)
      do iat = 1, mol%nat
         izp = mol%id(iat)
         norm = 0.0_wp
         dnorm = 0.0_wp
         d2norm = 0.0_wp
         do iref = 1, self%ref(izp)
            gw = weight_cn(self%wf, cn(iat), self%cn(iref, izp))
            gwt(iref) = gw
            dcn = self%cn(iref, izp) - cn(iat)
            norm = norm + gw
            dnorm = dnorm + wf2 * dcn * gw
            d2norm = d2norm + (wf2*wf2*dcn*dcn - wf2) * gw
         end do
         norm = 1.0_wp / norm
         do iref = 1, self%ref(izp)
            expw = gwt(iref)
            dcn = self%cn(iref, izp) - cn(iat)
            expd = wf2 * dcn * expw
            expdd = (wf2*wf2*dcn*dcn - wf2) * expw
            gwk = expw * norm
            if (is_exceptional(gwk)) then
               if (maxval(self%cn(:self%ref(izp), izp)) == self%cn(iref, izp)) then
                  gwk = 1.0_wp
               else
                  gwk = 0.0_wp
               end if
            end if
            gwvec(iref, iat) = gwk

            dgwk = expd * norm - expw * dnorm * norm**2
            if (is_exceptional(dgwk)) then
               dgwk = 0.0_wp
            end if
            gwdcn(iref, iat) = dgwk

            if (present(gwd2cn)) then
               d2gwk = expdd * norm - 2.0_wp * expd * dnorm * norm**2 &
                  & - expw * d2norm * norm**2 + 2.0_wp * expw * dnorm**2 * norm**3
               if (is_exceptional(d2gwk)) then
                  d2gwk = 0.0_wp
               end if
               gwd2cn(iref, iat) = d2gwk
            end if
         end do
      end do
      !$omp end do
      deallocate(gwt)
      !$omp end parallel

   else

      gwvec(:, :) = 0.0_wp

      !$omp parallel default(none) &
      !$omp shared(gwvec, mol, self, cn, mref) &
      !$omp private(iat, izp, iref, norm, gw, expw, gwk, gwt)
      allocate(gwt(mref))
      !$omp do schedule(runtime)
      do iat = 1, mol%nat
         izp = mol%id(iat)
         norm = 0.0_wp
         do iref = 1, self%ref(izp)
            gw = weight_cn(self%wf, cn(iat), self%cn(iref, izp))
            gwt(iref) = gw
            norm = norm + gw
         end do
         norm = 1.0_wp / norm
         do iref = 1, self%ref(izp)
            gwk = gwt(iref) * norm
            if (is_exceptional(gwk)) then
               if (maxval(self%cn(:self%ref(izp), izp)) == self%cn(iref, izp)) then
                  gwk = 1.0_wp
               else
                  gwk = 0.0_wp
               end if
            end if
            gwvec(iref, iat) = gwk
         end do
      end do
      !$omp end do
      deallocate(gwt)
      !$omp end parallel
   end if

end subroutine weight_references


!> Check whether we are dealing with an exceptional value, NaN or Inf
elemental function is_exceptional(val)
   real(wp), intent(in) :: val
   logical :: is_exceptional
   is_exceptional = ieee_is_nan(val) .or. abs(val) > huge(val)
end function is_exceptional


!> Calculate atomic dispersion coefficients and their derivatives w.r.t.
!> the coordination number.
!>
!> The second derivatives require a model without a low-rank expansion, which
!> only provides the coefficients and their first derivatives.
subroutine get_atomic_c6(self, mol, gwvec, gwdcn, c6, dc6dcn, gwd2cn, d2c6dcn2, &
      & d2c6dcnij, partition)

   !> Instance of the dispersion model
   class(d3_model), intent(in) :: self

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

   !> Weighting function for the atomic reference systems
   real(wp), intent(in) :: gwvec(:, :)

   !> Derivative of the weighting function w.r.t. the coordination number
   real(wp), intent(in), optional :: gwdcn(:, :)

   !> C6 coefficients for all atom pairs.
   real(wp), intent(out) :: c6(:, :)

   !> Derivative of the C6 w.r.t. the coordination number
   real(wp), intent(out), optional :: dc6dcn(:, :)

   !> Second derivative of the weighting function w.r.t. the coordination number
   real(wp), intent(in), optional :: gwd2cn(:, :)

   !> Second derivative of the C6 w.r.t. the coordination number of one atom
   real(wp), intent(out), optional :: d2c6dcn2(:, :)

   !> Mixed second derivative of the C6 w.r.t. both coordination numbers
   real(wp), intent(out), optional :: d2c6dcnij(:, :)

   !> Only evaluate the pairs of this part, the caller has to know that the
   !> missing coefficients are not consumed, see dftd3_partition
   type(work_partition), intent(in), optional :: partition

   integer :: iat, jat, izp, jzp, iref, jref
   real(wp) :: refc6, dc6, dc6dcni, dc6dcnj
   real(wp) :: d2c6dcni, d2c6dcnj, d2c6mix
   logical :: second

   if (allocated(self%lowrank)) then
      call self%lowrank%get_atomic_c6(mol, self%ghost, gwvec, gwdcn, c6, dc6dcn)
      return
   end if

   second = present(gwd2cn) .and. present(d2c6dcn2) .and. present(d2c6dcnij)

   if (second) then
      c6(:, :) = 0.0_wp
      dc6dcn(:, :) = 0.0_wp
      d2c6dcn2(:, :) = 0.0_wp
      d2c6dcnij(:, :) = 0.0_wp

      !$omp parallel do schedule(runtime) default(none) &
      !$omp shared(c6, dc6dcn, d2c6dcn2, d2c6dcnij, mol, self, gwvec, gwdcn, gwd2cn) &
      !$omp shared(partition) &
      !$omp private(iat, jat, izp, jzp, iref, jref, refc6, dc6, dc6dcni, dc6dcnj, &
      !$omp& d2c6dcni, d2c6dcnj, d2c6mix)
      do iat = 1, mol%nat
         if (self%ghost(iat)) cycle
         izp = mol%id(iat)
         do jat = 1, iat
            if (self%ghost(jat)) cycle
            if (.not.owns_pair(partition, iat, jat)) cycle
            jzp = mol%id(jat)
            dc6 = 0.0_wp
            dc6dcni = 0.0_wp
            dc6dcnj = 0.0_wp
            d2c6dcni = 0.0_wp
            d2c6dcnj = 0.0_wp
            d2c6mix = 0.0_wp
            do iref = 1, self%ref(izp)
               do jref = 1, self%ref(jzp)
                  refc6 = self%c6(iref, jref, izp, jzp)
                  dc6 = dc6 + gwvec(iref, iat) * gwvec(jref, jat) * refc6
                  dc6dcni = dc6dcni + gwdcn(iref, iat) * gwvec(jref, jat) * refc6
                  dc6dcnj = dc6dcnj + gwvec(iref, iat) * gwdcn(jref, jat) * refc6
                  d2c6dcni = d2c6dcni + gwd2cn(iref, iat) * gwvec(jref, jat) * refc6
                  d2c6dcnj = d2c6dcnj + gwvec(iref, iat) * gwd2cn(jref, jat) * refc6
                  d2c6mix = d2c6mix + gwdcn(iref, iat) * gwdcn(jref, jat) * refc6
               end do
            end do
            c6(iat, jat) = dc6
            c6(jat, iat) = dc6
            dc6dcn(iat, jat) = dc6dcni
            dc6dcn(jat, iat) = dc6dcnj
            d2c6dcn2(iat, jat) = d2c6dcni
            d2c6dcn2(jat, iat) = d2c6dcnj
            d2c6dcnij(iat, jat) = d2c6mix
            d2c6dcnij(jat, iat) = d2c6mix
         end do
      end do

   else if (present(gwdcn).and.present(dc6dcn)) then
      c6(:, :) = 0.0_wp
      dc6dcn(:, :) = 0.0_wp

      !$omp parallel do schedule(runtime) default(none) &
      !$omp shared(c6, dc6dcn, mol, self, gwvec, gwdcn, partition) &
      !$omp private(iat, jat, izp, jzp, iref, jref, refc6, dc6, dc6dcni, dc6dcnj)
      do iat = 1, mol%nat
         if (self%ghost(iat)) cycle
         izp = mol%id(iat)
         do jat = 1, iat
            if (self%ghost(jat)) cycle
            if (.not.owns_pair(partition, iat, jat)) cycle
            jzp = mol%id(jat)
            dc6 = 0.0_wp
            dc6dcni = 0.0_wp
            dc6dcnj = 0.0_wp
            do iref = 1, self%ref(izp)
               do jref = 1, self%ref(jzp)
                  refc6 = self%c6(iref, jref, izp, jzp)
                  dc6 = dc6 + gwvec(iref, iat) * gwvec(jref, jat) * refc6
                  dc6dcni = dc6dcni + gwdcn(iref, iat) * gwvec(jref, jat) * refc6
                  dc6dcnj = dc6dcnj + gwvec(iref, iat) * gwdcn(jref, jat) * refc6
               end do
            end do
            c6(iat, jat) = dc6
            c6(jat, iat) = dc6
            dc6dcn(iat, jat) = dc6dcni
            dc6dcn(jat, iat) = dc6dcnj
         end do
      end do

   else

      c6(:, :) = 0.0_wp

      !$omp parallel do schedule(runtime) default(none) &
      !$omp shared(c6, mol, self, gwvec, partition) &
      !$omp private(iat, jat, izp, jzp, iref, jref, refc6, dc6)
      do iat = 1, mol%nat
         if (self%ghost(iat)) cycle
         izp = mol%id(iat)
         do jat = 1, iat
            if (self%ghost(jat)) cycle
            if (.not.owns_pair(partition, iat, jat)) cycle
            jzp = mol%id(jat)
            dc6 = 0.0_wp
            do iref = 1, self%ref(izp)
               do jref = 1, self%ref(jzp)
                  refc6 = self%c6(iref, jref, izp, jzp)
                  dc6 = dc6 + gwvec(iref, iat) * gwvec(jref, jat) * refc6
               end do
            end do
            c6(iat, jat) = dc6
            c6(jat, iat) = dc6
         end do
      end do
   end if

end subroutine get_atomic_c6


elemental function weight_cn(wf,cn,cnref) result(cngw)
   real(wp),intent(in) :: wf, cn, cnref
   real(wp) :: cngw
   intrinsic :: exp
   cngw = exp ( -wf * ( cn - cnref )**2 )
end function weight_cn


end module dftd3_model