compute tensor coefficients needed to compute open-ended derivatives of (from C. E. Dykstra, J. Comput. Chem., 9 (1988), 476, where C^(n)_ij are stored as tensor_coefficient(j,i,n))
subroutine compute_tensor_coefficients()
integer :: i, j, k, n
tensor_coefficient(:,:,:) = 0.0_dp
tensor_coefficient(0,0,:) = 1.0_dp
do n = 1, 2*MAX_ORDER + 1
if (mod(n, 2) == 0) cycle
do i = 1, MAX_ORDER
k = merge(i, i - 1, (mod(i, 2) == 0))
do j = 0, i
if (mod(i + j, 2) /= 0) cycle
if (j == 0) then
tensor_coefficient(j,i,n) = tensor_coefficient(j+1,i-1,n)
else if (j /= i) then
tensor_coefficient(j,i,n) = (real(j, dp) + 1) * tensor_coefficient(j+1,i-1,n)
tensor_coefficient(j,i,n) = tensor_coefficient(j,i,n) &
- (real(n, dp) + real(k, dp)) * tensor_coefficient(j-1,i-1,n)
k = k + 2
else if (j == i) then
tensor_coefficient(j,i,n) = - (real(n, dp) + real(k, dp)) * tensor_coefficient(j-1,i-1,n)
end if
end do
end do
end do
end subroutine compute_tensor_coefficients