compute_electronic_multipoles Subroutine

subroutine compute_electronic_multipoles(this, expansion_origin, expansion_order, &
                                         a_start, a_end)

    class(quantum_fragment_type), intent(inout) :: this
    !> origin of the multipole expansion
    real(dp), dimension(3), intent(in) :: expansion_origin
    !> order of the multipole expansion
    integer, intent(in) :: expansion_order
    !> start index of real space mesh of lattice vector a
    integer, intent(in) :: a_start
    !> end index of real space mesh of lattice vector a
    integer, intent(in) :: a_end

    integer :: i, j, k
    integer :: x, y, z
    integer :: a, b, c
    integer :: a_index
    integer :: mi
    integer :: order
    real(dp) :: dV
    real(dp) :: temp
    real(dp), dimension(3) :: ra, rb, rc, r
    real(dp), dimension(3) :: cell_origin
    real(dp), dimension(3,3) :: lattice_steps
    real(dp), dimension(:), allocatable :: multipoles

    call timer_start("compute_electronic_multipoles")

    if (a_start > a_end) then
        call handle_error(SEVERITY_FATAL, &
                          TYPE_INCORRECT_ARG, &
                          "a_start should be smaller than a_end", &
                          __FILE__, __LINE__)
    endif

    allocate(multipoles(polytensor_size(expansion_order)))

    multipoles = 0.0_dp

    cell_origin = this%cell%origin

    do i = 1, 3
        lattice_steps(:,i) = this%cell%lattice(:,i) / real(this%cell%num_points(i), dp)
    end do

    dV = this%cell%volume / real(product(this%cell%num_points), dp)

    associate (density => this%density%field)
    !$OMP PARALLEL DO PRIVATE(order, i, j, k, c, rc, b, rb, a, a_index, ra, r, temp, &
    !$OMP& x, y, z, mi) REDUCTION(+:multipoles)
    do c = 1, this%cell%num_points(3)
        rc = cell_origin + real(c - 1, dp) * lattice_steps(:,3)
        do b = 1, this%cell%num_points(2)
            rb = rc + real(b - 1, dp) * lattice_steps(:,2)
            do a = a_start, a_end
                a_index = a - a_start + 1
                ra = rb + real(a - 1, dp) * lattice_steps(:,1)
                r = ra - expansion_origin
                mi = 1
                do order = 0, expansion_order
                    do i = order, 0, -1
                        do j = order - i, 0, -1
                            k = order - i - j
                            temp = - density(a_index,b,c) * dV
                            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(mi) = multipoles(mi) + temp
                            mi = mi + 1
                        end do
                    end do
                end do
            end do
        end do
    end do
    !$OMP END PARALLEL DO
    end associate

    call this%electronic_multipoles%init(id=0, origin=expansion_origin, order=expansion_order, values=multipoles)

    deallocate(multipoles)

    call timer_stop

end subroutine compute_electronic_multipoles