get_dispersion2_hessian Subroutine

   subroutine get_dispersion2_hessian(self, mol, trans, cutoff, width, rvdw, r4r2, &
         & c6, dc6dcn, d2c6dcn2, d2c6dcnij, hessian, dEdcn, dEdcndr, dEdcndcn, partition)

      !> Damping parameters
      class(damping_param), intent(in) :: self

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

      !> Lattice points
      real(wp), intent(in) :: trans(:, :)

      !> Real space cutoff
      real(wp), intent(in) :: cutoff

      !> Width of smooth cutoff
      real(wp), intent(in) :: width

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

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

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

      !> Derivatives of the C6 w.r.t. the coordination number
      real(wp), intent(in) :: dc6dcn(:, :), d2c6dcn2(:, :), d2c6dcnij(:, :)

      !> Second derivative of the energy w.r.t. the Cartesian coordinates
      real(wp), intent(inout) :: hessian(:, :)

      !> Derivative of the energy w.r.t. the coordination number
      real(wp), intent(inout) :: dEdcn(:)

      !> Mixed derivative w.r.t. coordination number and Cartesian coordinates
      real(wp), intent(inout) :: dEdcndr(:, :)

      !> Second derivative w.r.t. the coordination numbers
      real(wp), intent(inout) :: dEdcndcn(:, :)

      !> Work partition of the atom pairs
      type(work_partition), intent(in), optional :: partition

      integer :: iat, jat, izp, jzp, jtr, ic, jc, ii, jj, ia, ib, nc
      integer :: cnat(2)
      real(wp) :: vec(3), r2, cutoff2, c6ij, fac
      real(wp) :: e0, e0r, e0c, e0rr, e0rc, e0cc
      real(wp) :: sw, swr, swrr, fr, frr, fc, frc, fcc
      real(wp) :: block(3, 3), cnd(2), cnd2(2, 2), dr2(3, 2)

      ! Kept serial on purpose: the loop is O(N^2) but bound by the scattered
      ! writes into the O(N^2) hessian, so thread-private copies only add traffic
      cutoff2 = cutoff*cutoff

      do iat = 1, mol%nat
         izp = mol%id(iat)
         do jat = 1, iat
            if (.not.owns_pair(partition, iat, jat)) cycle
            jzp = mol%id(jat)
            c6ij = c6(jat, iat)
            if (iat /= jat) then
               fac = 1.0_wp
               nc = 2
               cnat(1) = iat
               cnat(2) = jat
               cnd(1) = dc6dcn(iat, jat)
               cnd(2) = dc6dcn(jat, iat)
               cnd2(1, 1) = d2c6dcn2(iat, jat)
               cnd2(2, 2) = d2c6dcn2(jat, iat)
               cnd2(1, 2) = d2c6dcnij(iat, jat)
               cnd2(2, 1) = d2c6dcnij(iat, jat)
            else
               fac = 0.5_wp
               nc = 1
               cnat(1) = iat
               cnd(1) = 2.0_wp*dc6dcn(iat, iat)
               cnd2(1, 1) = 2.0_wp*d2c6dcn2(iat, iat) + 2.0_wp*d2c6dcnij(iat, iat)
            end if
            do jtr = 1, size(trans, 2)
               vec(:) = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, jtr))
               r2 = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
               if (r2 > cutoff2 .or. r2 < epsilon(1.0_wp)) cycle

               call self%get_damping_kernel(izp, jzp, rvdw, r4r2, r2, c6ij, &
                  & e0, e0r, e0c, e0rr, e0rc, e0cc)
               call smooth_cutoff_r2(r2, cutoff, width, sw, swr, swrr)

               fr = fac*(swr*e0 + sw*e0r)
               frr = fac*(swrr*e0 + 2.0_wp*swr*e0r + sw*e0rr)
               fc = fac*sw*e0c
               frc = fac*(swr*e0c + sw*e0rc)
               fcc = fac*sw*e0cc

               ! Cartesian second derivatives at fixed coordination number
               if (iat /= jat) then
                  do ic = 1, 3
                     do jc = 1, 3
                        block(ic, jc) = 4.0_wp*frr*vec(ic)*vec(jc)
                     end do
                     block(ic, ic) = block(ic, ic) + 2.0_wp*fr
                  end do
                  do ic = 1, 3
                     do jc = 1, 3
                        hessian(3*(iat-1)+ic, 3*(iat-1)+jc) = &
                           & hessian(3*(iat-1)+ic, 3*(iat-1)+jc) + block(ic, jc)
                        hessian(3*(jat-1)+ic, 3*(jat-1)+jc) = &
                           & hessian(3*(jat-1)+ic, 3*(jat-1)+jc) + block(ic, jc)
                        hessian(3*(iat-1)+ic, 3*(jat-1)+jc) = &
                           & hessian(3*(iat-1)+ic, 3*(jat-1)+jc) - block(ic, jc)
                        hessian(3*(jat-1)+ic, 3*(iat-1)+jc) = &
                           & hessian(3*(jat-1)+ic, 3*(iat-1)+jc) - block(ic, jc)
                     end do
                  end do
                  dr2(:, 1) = 2.0_wp*vec
                  dr2(:, 2) = -2.0_wp*vec
               else
                  dr2(:, :) = 0.0_wp
               end if

               ! Coordination number dependence
               do ia = 1, nc
                  dEdcn(cnat(ia)) = dEdcn(cnat(ia)) + fc*cnd(ia)
                  do ib = 1, nc
                     dEdcndcn(cnat(ia), cnat(ib)) = dEdcndcn(cnat(ia), cnat(ib)) &
                        & + fcc*cnd(ia)*cnd(ib) + fc*cnd2(ia, ib)
                  end do
               end do

               if (iat /= jat) then
                  do ia = 1, nc
                     do ic = 1, 3
                        ii = 3*(iat-1) + ic
                        jj = 3*(jat-1) + ic
                        dEdcndr(ii, cnat(ia)) = dEdcndr(ii, cnat(ia)) &
                           & + frc*dr2(ic, 1)*cnd(ia)
                        dEdcndr(jj, cnat(ia)) = dEdcndr(jj, cnat(ia)) &
                           & + frc*dr2(ic, 2)*cnd(ia)
                     end do
                  end do
               end if
            end do
         end do
      end do

   end subroutine get_dispersion2_hessian