get_atomic_c6_hessian Subroutine

subroutine get_atomic_c6_hessian(self, mol, gwvec, gwdcn, gwd2cn, c6, &
      & dc6dcn, d2c6dcn2, d2c6dcnij)

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

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

   !> Reference weights and their CN derivatives
   real(wp), intent(in) :: gwvec(:, :, :), gwdcn(:, :, :), gwd2cn(:, :, :)

   !> C6 coefficients and their CN derivatives
   real(wp), intent(out) :: c6(:, :), dc6dcn(:, :)
   real(wp), intent(out) :: d2c6dcn2(:, :), d2c6dcnij(:, :)

   integer :: ipair, npair, iat, jat, izp, jzp, iref, jref
   real(wp) :: refc6, c6ij, dc6dcni, dc6dcnj, d2c6dcni, d2c6dcnj, d2c6mix

   ! Flatten the lower triangle so OpenMP can schedule O(N^2) pair tasks.
   npair = mol%nat*(mol%nat + 1)/2
   !$omp parallel do default(none) schedule(runtime) &
   !$omp shared(c6, dc6dcn, d2c6dcn2, d2c6dcnij, mol, self, gwvec, gwdcn, gwd2cn, npair) &
   !$omp private(ipair, iat, jat, izp, jzp, iref, jref, refc6, c6ij, dc6dcni, dc6dcnj, &
   !$omp& d2c6dcni, d2c6dcnj, d2c6mix)
   do ipair = 1, npair
      iat = int(0.5_wp*(sqrt(8.0_wp*real(ipair, wp) + 1.0_wp) - 1.0_wp))
      if (iat*(iat + 1)/2 < ipair) iat = iat + 1
      jat = ipair - iat*(iat - 1)/2
      izp = mol%id(iat)
      jzp = mol%id(jat)
      c6ij = 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)
            c6ij = c6ij + gwvec(iref, iat, 1)*gwvec(jref, jat, 1)*refc6
            dc6dcni = dc6dcni + gwdcn(iref, iat, 1)*gwvec(jref, jat, 1)*refc6
            dc6dcnj = dc6dcnj + gwvec(iref, iat, 1)*gwdcn(jref, jat, 1)*refc6
            d2c6dcni = d2c6dcni + gwd2cn(iref, iat, 1)*gwvec(jref, jat, 1)*refc6
            d2c6dcnj = d2c6dcnj + gwvec(iref, iat, 1)*gwd2cn(jref, jat, 1)*refc6
            d2c6mix = d2c6mix + gwdcn(iref, iat, 1)*gwdcn(jref, jat, 1)*refc6
         end do
      end do
      c6(iat, jat) = c6ij
      c6(jat, iat) = c6ij
      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 subroutine get_atomic_c6_hessian