get_dispersion_api Subroutine

subroutine get_dispersion_api(verror, vmol, vdisp, vparam, &
      & energy, c_gradient, c_sigma) &
      & bind(C, name=namespace//"get_dispersion")
   !DEC$ ATTRIBUTES DLLEXPORT :: get_dispersion_api
   type(c_ptr), value :: verror
   type(vp_error), pointer :: error
   type(c_ptr), value :: vmol
   type(vp_structure), pointer :: mol
   type(c_ptr), value :: vdisp
   type(vp_model), pointer :: disp
   type(c_ptr), value :: vparam
   type(vp_param), pointer :: param
   real(c_double), intent(out) :: energy
   real(c_double), intent(out), optional :: c_gradient(3, *)
   real(wp), allocatable :: gradient(:, :)
   real(c_double), intent(out), optional :: c_sigma(3, 3)
   real(wp), allocatable :: sigma(:, :)
   logical :: has_grad, has_sigma


   if (debug) print'("[Info]",1x, a)', "get_dispersion"

   if (.not.c_associated(verror)) return
   call c_f_pointer(verror, error)

   if (.not.c_associated(vmol)) then
      call fatal_error(error%ptr, "Molecular structure data is missing")
      return
   end if
   call c_f_pointer(vmol, mol)

   if (.not.c_associated(vdisp)) then
      call fatal_error(error%ptr, "Dispersion model is missing")
      return
   end if
   call c_f_pointer(vdisp, disp)

   if (.not.c_associated(vparam)) then
      call fatal_error(error%ptr, "Damping parameters are missing")
      return
   end if
   call c_f_pointer(vparam, param)

   if (.not.allocated(param%ptr)) then
      call fatal_error(error%ptr, "Damping parameters are not initialized")
      return
   end if

   has_grad = present(c_gradient) 
   if (has_grad) then
      gradient = c_gradient(:3, :mol%ptr%nat)
   endif

   has_sigma = present(c_sigma) 
   if (has_sigma) then
      sigma = c_sigma(:3, :3)
   ! Still needs to be passed into dispersion subroutines,
   ! just won't be returned through the API. 
   ! Would need to refactor disperision
   ! subroutines to make sigma truly optional. 
   else if (has_grad) then
      allocate(sigma(3,3)) 
   endif

   ! Evaluate energy, gradient (optional), and 
   ! sigma (optional) analytically
   call get_dispersion(mol%ptr, disp%ptr, param%ptr, realspace_cutoff(), &
      & energy, gradient, sigma)

   if (has_grad) then
      c_gradient(:3, :mol%ptr%nat) = gradient
   endif

   if (has_sigma) then
      c_sigma(:3, :3) = sigma
   endif

end subroutine get_dispersion_api