get_atm_dispersion Subroutine

subroutine get_atm_dispersion(mol, trans, cutoff, s9, a1, a2, alp, r4r2, &
      & c6, dc6dcn, dc6dq, energy, dEdcn, dEdq, gradient, sigma)

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

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

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

   !> Scaling for dispersion coefficients
   real(wp), intent(in) :: s9

   !> Scaling parameter for critical radius
   real(wp), intent(in) :: a1

   !> Offset parameter for critical radius
   real(wp), intent(in) :: a2

   !> Exponent of zero damping function
   real(wp), intent(in) :: alp

   !> Expectation values for r4 over r2 operator
   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), optional :: dc6dcn(:, :)

   !> Derivative of the C6 w.r.t. the partial charges
   real(wp), intent(in), optional :: dc6dq(:, :)

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

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

   !> Derivative of the energy w.r.t. the partial charges
   real(wp), intent(inout), optional :: dEdq(:)

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

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

   logical :: grad

   if (abs(s9) < epsilon(1.0_wp)) return
   grad = present(dc6dcn) .and. present(dEdcn) .and. present(dc6dq) &
      & .and. present(dEdq) .and. present(gradient) .and. present(sigma)

   if (grad) then
      call get_atm_dispersion_derivs(mol, trans, cutoff, s9, a1, a2, alp, r4r2, &
         & c6, dc6dcn, dc6dq, energy, dEdcn, dEdq, gradient, sigma)
   else
      call get_atm_dispersion_energy(mol, trans, cutoff, s9, a1, a2, alp, r4r2, &
         & c6, energy)
   end if

end subroutine get_atm_dispersion