get_counterpoise_hessian_api Subroutine

subroutine get_counterpoise_hessian_api(verror, vmol, vgcp, &
      & energy, c_hessian) &
      & bind(C, name=namespace//"get_counterpoise_hessian")
   type(c_ptr), value :: verror
   type(vp_error), pointer :: error
   type(c_ptr), value :: vmol
   type(vp_structure), pointer :: mol
   type(c_ptr), value :: vgcp
   type(vp_gcp), pointer :: gcp
   real(c_double), intent(out) :: energy
   real(c_double), intent(out) :: c_hessian(*)
   real(wp), allocatable :: hessian(:, :)
   type(realspace_cutoff) :: cutoff
   integer :: ndim

   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(vgcp)) then
      call fatal_error(error%ptr, "Counter-poise parameters are missing")
      return
   end if
   call c_f_pointer(vgcp, gcp)

   ndim = 3*mol%ptr%nat
   allocate(hessian(ndim, ndim))

   cutoff = realspace_cutoff()
   if (allocated(gcp%cutoff)) then
      cutoff = gcp%cutoff
   end if
   energy = 0.0_wp
   if (allocated(gcp%comm)) then
      call get_counterpoise_mpi(error%ptr, mol%ptr, gcp%ptr, cutoff, gcp%comm, &
         & energy, hessian=hessian)
      if (allocated(error%ptr)) return
   else
      call get_geometric_counterpoise(mol%ptr, gcp%ptr, cutoff, energy, &
         & partition=gcp%partition)
      call get_geometric_counterpoise_hessian(mol%ptr, gcp%ptr, cutoff, hessian, &
         & gcp%partition)
   end if

   c_hessian(:ndim*ndim) = reshape(hessian, [ndim*ndim])

end subroutine get_counterpoise_hessian_api