compute_lr_potential Subroutine

subroutine compute_lr_potential(quantum_fragment, &
                                a_start, a_end, tensor_sums)

    !> quantum fragment
    type(quantum_fragment_type), intent(in) :: quantum_fragment
    !> 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
    !> Array stroing tensor sums
    real(dp), dimension(:), intent(in) :: tensor_sums

    integer :: i, j, k
    integer :: x, y, z
    integer :: a, b, c
    integer :: a_index
    integer :: comp
    integer :: order
    integer :: multipole_order
    real(dp) :: temp
    real(dp), dimension(3) :: ra, rb, rc, r
    real(dp), dimension(3) :: multipole_origin
    real(dp), dimension(3,3) :: lattice_steps

    call timer_start("compute_lr_potential")

    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

    multipole_origin = quantum_fragment%nuclear_multipoles%origin
    multipole_order = quantum_fragment%nuclear_multipoles%order

    associate (cell => quantum_fragment%cell)

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

    associate (potential => quantum_fragment%potential%field)

    !$OMP PARALLEL DO private(c, rc, b, rb, a, a_index, ra, r, &
    !$OMP& comp, order, i, j, k, temp)
    do c = 1, cell%num_points(3)
        rc = cell%origin + real(c - 1, dp) * lattice_steps(:,3)
        do b = 1, 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 - multipole_origin
                comp = 1
                do order = 0, multipole_order
                    do i = order, 0, - 1
                        do j = order - i, 0, - 1
                            k = order - i - j
                            temp = 1.0_dp
                            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
                            potential(a_index,b,c) = potential(a_index,b,c) + temp * tensor_sums(comp)
                            comp = comp + 1
                        end do
                    end do
                end do
            end do
        end do
    end do
    !$OMP END PARALLEL DO

    end associate

    end associate

    call timer_stop

end subroutine compute_lr_potential