subroutine add_ecpder(basis, coord, denab, de)
real(real64), contiguous, intent(in) :: coord(:,:)
type(basis_set), intent(in) :: basis
REAL(kind=dp), INTENT(INOUT) :: denab(:)
REAL(kind=dp), intent(INOUT) :: de(:,:)
type(ecp_result) :: result_ptr
type(c_ptr) :: integrator
real(c_double), pointer :: libecp_res(:)
real(real64), allocatable :: res_x(:), res_y(:), res_z(:)
real(real64), allocatable :: deloc(:,:)
integer :: i, j, c, n, natm, prim
integer :: tri_size, full_size
integer(c_int) :: driv_order
if (.not.(basis%ecp_params%is_ecp)) then
return
end if
driv_order = 1
tri_size = basis%nbf * (basis%nbf + 1) / 2
full_size = basis%nbf * basis%nbf
allocate(res_x(full_size))
allocate(res_y(full_size))
allocate(res_z(full_size))
natm = size(coord, dim=2)
allocate(deloc(3, natm))
deloc = 0
call set_integrator(integrator, basis, coord, driv_order)
result_ptr = compute_first_derivs(integrator)
call c_f_pointer(result_ptr%data, libecp_res, [result_ptr%size])
do n = 1, natm
res_x = libecp_res(full_size * (3 * n - 3)+1:full_size * (3 * n - 2))
res_y = libecp_res(full_size * (3 * n - 2)+1:full_size * (3 * n - 1))
res_z = libecp_res(full_size * (3 * n - 1)+1:full_size * (3 * n))
call transform_matrix(basis, res_x)
call transform_matrix(basis, res_y)
call transform_matrix(basis, res_z)
do j = 1, basis%nbf
do i = 1, j
c = j * (j - 1) / 2 + i
if (i == j) then
prim = 1
else
prim = 2
end if
deloc(1, n) = deloc(1, n) + prim * res_x((i - 1) * basis%nbf + j) * denab(c)
deloc(2, n) = deloc(2, n) + prim * res_y((i - 1) * basis%nbf + j) * denab(c)
deloc(3, n) = deloc(3, n) + prim * res_z((i - 1) * basis%nbf + j) * denab(c)
end do
end do
end do
de(:, 1:natm) = de(:, 1:natm) + deloc(:, 1:natm)
result_ptr%data = c_null_ptr
result_ptr%size = 0
nullify(libecp_res)
call free_integrator(integrator)
call free_result(result_ptr)
end subroutine add_ecpder