compute_nuclear_multipoles Subroutine

subroutine compute_nuclear_multipoles(this, origin, order)

    class(quantum_fragment_type), intent(inout) :: this
    !> origin of the multipole expansion
    real(dp), dimension(3), intent(in) :: origin
    !> order of the multipole expansion
    integer, intent(in) :: order

    integer :: nuc, ord, comp
    integer :: i, j, k
    integer :: x, y, z
    real(dp) :: temp
    real(dp), dimension(3) :: r
    real(dp), dimension(:), allocatable :: multipoles

    call timer_start("compute_nuclear_multipoles")

    allocate(multipoles((order+3)*(order+2)*(order+1)/6))

    multipoles = 0.0_dp

    do nuc = 1, this%num_nuclei
        r = this%nuclei(nuc)%coordinate - origin
        comp = 1
        do ord = 0, order
            do i = ord, 0, -1
                do j = ord - i, 0, -1
                    k = ord - i - j
                    temp = this%nuclei(nuc)%charge
                    do x = 1, i
                        temp = temp * r(1)
                    end do
                    do y = 1, j
                        temp = temp * r(2)
                    end do
                    do z = 1, k
                        temp = temp * r(3)
                    end do
                    multipoles(comp) = multipoles(comp) + temp
                    comp = comp + 1
                end do
            end do
        end do
    end do

    call this%nuclear_multipoles%init(id=0, origin=origin, order=order, values=multipoles)

    deallocate(multipoles)

    call timer_stop

end subroutine compute_nuclear_multipoles