compute_folded_tensors Subroutine

subroutine compute_folded_tensors(tensors, tensor_sums, &
                                  multipole_origin, multipole_order, &
                                  atoms)

    !> Array storing tensors
    real(dp), dimension(:), intent(out) :: tensors
    !> Array stroing tensor sums
    real(dp), dimension(:), intent(out) :: tensor_sums
    !> Coordinates of the origin of multipolar expansion
    real(dp), dimension(3), intent(in) :: multipole_origin
    !> Maximal order of multipoles
    integer, intent(in) :: multipole_order
    !> long-range atoms
    type(atom_type), dimension(:), intent(in) :: atoms

    integer :: i, j, k, order
    real(dp) :: taylor
    real(dp), dimension(3) :: r

    call timer_start("compute_folded_tensors")

    tensors = 0.0_dp
    tensor_sums = 0.0_dp

    !$OMP PARALLEL DO private(i, j, r, k, order, taylor, tensors) &
    !$OMP& reduction(+:tensor_sums)
    do i = 1, size(atoms)
        r = multipole_origin - atoms(i)%coordinate
        k = 0
        do order = 0, multipole_order
            j = k + 1
            k = k + tensor_size(order)
            call folded_interaction_tensor(r, tensors(j:k), order, 0)
            taylor = merge(1.0_dp, - 1.0_dp, (mod(order, 2) == 0)) / factorial(order)
            tensors(j:k) = taylor * tensors(j:k)
        end do
        tensor_sums = tensor_sums + atoms(i)%charge * tensors
    end do
    !$OMP END PARALLEL DO

    call timer_stop

end subroutine compute_folded_tensors