Compute folded multipole interaction tensor containing derivatives of , i.e. the da'th derivative wrt and db'th derivative wrt . Tensor is returned in a one-dimensional array with elements in anticanonical ordering, e.g. xx, xy, xz, yy, yz, zz.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp), | intent(in), | dimension(:) | :: | r_ab | coordinate |
|
real(kind=dp), | intent(out), | dimension(:) | :: | T | the interaction tensor |
|
integer, | intent(in) | :: | da | order of derivative wrt |
||
integer, | intent(in) | :: | db | order of derivative wrt |
pure subroutine folded_interaction_tensor(r_ab, T, da, db)
!> coordinate \(\mathbf{r}_{ab} = \mathbf{r}_{b}-\mathbf{r}_{a}\)
real(dp), dimension(:), intent(in) :: r_ab
!> the interaction tensor
real(dp), dimension(:), intent(out) :: T
!> order of derivative wrt \(\mathbf{r}_{a}\)
integer, intent(in) :: da
!> order of derivative wrt \(\mathbf{r}_{b}\)
integer, intent(in) :: db
integer :: i, j, k, l
integer :: error_code
real(dp) :: sgn
error_code = 0
if (size(T) /= (da + db + 1) * (da + db + 2) / 2) then
! 'MiMiC ERROR in tensor: wrong size of T_ba'
error_code = error_code + 1
end if
! tensor_elements gives derivatives of 1 / |r_b - r_a| wrt r_b, so therefore the sign
! alternates for each order derivative, and thus we here change sign for equal order
! derivatives with respect to r_a, because it should not change sign
sgn = merge(-1.0_dp, 1.0_dp, (mod(da, 2) /= 0))
l = 1
do i = da + db, 0, - 1
do j = da + db - i, 0, - 1
k = da + db - i - j
T(l) = sgn * trinomial(i, j, k) * tensor_element(i, j, k, r_ab)
l = l + 1
end do
end do
end subroutine folded_interaction_tensor