get_dispersion_hessian Subroutine

subroutine get_dispersion_hessian(mol, disp, param, cutoff, hessian)
   !DEC$ ATTRIBUTES DLLEXPORT :: get_dispersion_hessian

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

   !> Dispersion model
   class(dispersion_model), intent(in) :: disp

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

   !> Realspace cutoffs
   type(realspace_cutoff), intent(in) :: cutoff

   !> Dispersion hessian
   real(wp), intent(out) :: hessian(:, :, :, :)

   integer :: iat, ix
   real(wp), parameter :: step = 1.0e-4_wp
   type(structure_type) :: displ
   real(wp) :: el, er
   real(wp), allocatable :: gl(:, :), gr(:, :), sl(:, :), sr(:, :)

   hessian(:, :, :, :) = 0.0_wp
   !$omp parallel default(none) &
   !$omp private(iat, ix, displ, er, el, gr, gl, sr, sl) &
   !$omp shared(mol, disp, param, cutoff, hessian)
   displ = mol
   allocate(gl(3, mol%nat), gr(3, mol%nat), sl(3, 3), sr(3, 3))
   !$omp do schedule(dynamic) collapse(2)
   do iat = 1, mol%nat
      do ix = 1, 3
         displ%xyz(ix, iat) = mol%xyz(ix, iat) + step
         call get_dispersion(displ, disp, param, cutoff, el, gl, sl)

         displ%xyz(ix, iat) = mol%xyz(ix, iat) - step
         call get_dispersion(displ, disp, param, cutoff, er, gr, sr)

         displ%xyz(ix, iat) = mol%xyz(ix, iat)
         hessian(:, :, ix, iat) = (gl - gr) / (2 * step)
      end do
   end do
   !$omp end parallel
end subroutine get_dispersion_hessian