get_dispersion_derivs Subroutine

subroutine get_dispersion_derivs(self, mol, trans, cutoff, rvdw, r4r2, c6, dc6dcn, &
      & energy, dEdcn, gradient, sigma)

   !> Damping parameters
   class(mzero_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

   !> 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(:, :)

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

   !> Dispersion energy
   real(wp), intent(inout) :: energy(:)

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

   !> Dispersion gradient
   real(wp), intent(inout) :: gradient(:, :)

   !> Dispersion virial
   real(wp), intent(inout) :: sigma(:, :)

   integer :: iat, jat, izp, jzp, jtr
   real(wp) :: vec(3), r2, r1, r6, r8, t6, t8, d6, d8, f6, f8, alp6, alp8
   real(wp) :: edisp, gdisp, cutoff2, r0ij, rrij, c6ij, dE, dG(3), dS(3, 3)

   cutoff2 = cutoff*cutoff
   alp6 = self%alp
   alp8 = self%alp + 2.0_wp

   !$omp parallel do schedule(runtime) default(none) &
   !$omp reduction(+:energy, gradient, sigma, dEdcn) &
   !$omp shared(mol, self, c6, dc6dcn, trans, cutoff2, alp6, alp8, rvdw, r4r2) &
   !$omp private(iat, jat, izp, jzp, jtr, vec, r2, r1, r6, r8, t6, t8, d6, &
   !$omp& d8, f6, f8, edisp, gdisp, r0ij, rrij, c6ij, dE, dG, dS)
   do iat = 1, mol%nat
      izp = mol%id(iat)
      do jat = 1, iat
         jzp = mol%id(jat)
         rrij = 3*r4r2(izp)*r4r2(jzp)
         r0ij = rvdw(jzp, izp)
         c6ij = c6(jat, iat)
         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
            r1 = sqrt(r2)

            r6 = r2*r2*r2
            r8 = r6*r2

            t6 = (r1/(self%rs6*r0ij)+self%bet*r0ij)**(-alp6)
            t8 = (r1/(self%rs8*r0ij)+self%bet*r0ij)**(-alp8)

            f6 = 1.0_wp / (1.0_wp + 6.0_wp*t6)
            f8 = 1.0_wp / (1.0_wp + 6.0_wp*t8)

            d6 = -6.0_wp * f6 / r2 &
               & + 6.0_wp*alp6*t6*f6**2 / (r2+self%rs6*self%bet*r0ij**2*r1)
            d8 = -8.0_wp * f8 / r2 &
               & + 6.0_wp*alp8*t8*f8**2 / (r2+self%rs8*self%bet*r0ij**2*r1)

            edisp = self%s6 * f6 / r6 + self%s8 * rrij * f8 / r8
            gdisp = self%s6 * d6 / r6 + self%s8 * rrij * d8 / r8

            dE = -c6ij*edisp * 0.5_wp
            dG(:) = -c6ij*gdisp*vec
            dS(:, :) = spread(dG, 1, 3) * spread(vec, 2, 3) * 0.5_wp

            energy(iat) = energy(iat) + dE
            dEdcn(iat) = dEdcn(iat) - dc6dcn(iat, jat) * edisp
            sigma(:, :) = sigma + dS
            if (iat /= jat) then
               energy(jat) = energy(jat) + dE
               dEdcn(jat) = dEdcn(jat) - dc6dcn(jat, iat) * edisp
               gradient(:, iat) = gradient(:, iat) + dG
               gradient(:, jat) = gradient(:, jat) - dG
               sigma(:, :) = sigma + dS
            end if
         end do
      end do
   end do

end subroutine get_dispersion_derivs