compute_sr_potential_electrons Subroutine

subroutine compute_sr_potential_electrons(atoms, potential, cell, &
                                          x_start, x_end, id_start, id_end)
    USE, INTRINSIC :: ISO_C_BINDING
    !> list of SR atoms
    type(atom_type), dimension(:), intent(in) :: atoms
    !> potential on the electronic grid
    real(dp), dimension(:,:,:), contiguous, intent(inout), target :: potential
    !> lattice data
    type(cell_type), intent(in) :: cell
    !> index of starting X-plane (CPMD parallelization stuff)
    integer, intent(in) :: x_start
    !> index of final X-plane (CPMD parallelization stuff)
    integer, intent(in) :: x_end
    !> Index of the atom to start electrostatic treatment
    integer, intent(in) :: id_start
    !> Index of the atom to finish electrostatic treatment
    integer, intent(in) :: id_end

    ! cell elementary vectors
    real(dp), dimension(3) :: da, db, dc
    ! coordinate of the j-th atom
    real(dp), dimension(3) :: r_j
    real(dp) :: dV
    real(dp) :: rij, x_i, y_i, z_i, x_ij, y_ij, z_ij
    real(dp) :: rij_sq
    real(dp) :: charge
    real(dp) :: inv_denominator
    real(dp) :: point_potential
    real(dp) :: cov_r_4, cov_r_5

    integer :: n_atom
    integer :: nz_step, n
    integer :: ny_step
    integer :: nx_step
    integer :: n_points
    integer :: nx, ny, nz, nx_real
    integer, dimension(:), allocatable :: x_ind, y_ind
    real(dp), dimension(:), allocatable :: grid_x, grid_y, grid_z
    real(dp), dimension(:,:), contiguous, pointer :: potential_p

    call timer_start("compute_sr_potential_electrons")

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

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

    n_points = cell%num_points(1) * cell%num_points(2) * cell%num_points(3)
    dV = cell%volume / real(n_points, dp)
    da = cell%lattice(:,1) / real(cell%num_points(1), dp)
    db = cell%lattice(:,2) / real(cell%num_points(2), dp)
    dc = cell%lattice(:,3) / real(cell%num_points(3), dp)
    nx = cell%num_points_r(1)
    ny = cell%num_points_r(2)
    nz = cell%num_points(3)
   
    call c_f_pointer(c_loc(potential), potential_p, [nx * ny, nz])

    allocate(x_ind(nx * ny))
    allocate(y_ind(nx * ny))
    allocate(grid_x(nx * ny))
    allocate(grid_y(nx * ny))
    allocate(grid_z(nx * ny))

    !$OMP PARALLEL DO private(nx_step)
    do ny_step = 1, ny
        do nx_step = 1, nx
            grid_x((ny_step - 1) * nx + nx_step) = db(1) * real(ny_step - 1, dp) + da(1) * (nx_step + x_start - 2)
            grid_y((ny_step - 1) * nx + nx_step) = db(2) * real(ny_step - 1, dp) + da(2) * (nx_step + x_start - 2)
            grid_z((ny_step - 1) * nx + nx_step) = db(3) * real(ny_step - 1, dp) + da(3) * (nx_step + x_start - 2)
        enddo
    enddo
    !$OMP END PARALLEL DO

    ! first calculate action of classical atoms on the electronic grid
    do n_atom = id_start, id_end
        if (abs(atoms(n_atom)%charge) <= 1.0e-12_dp) cycle
        cov_r_4 = atoms(n_atom)%radius**ORDER
        cov_r_5 = cov_r_4 * atoms(n_atom)%radius
        r_j = atoms(n_atom)%coordinate
        charge = atoms(n_atom)%charge
        ! z-loop over grid
        !$OMP PARALLEL DO SCHEDULE(DYNAMIC) DEFAULT(SHARED) &
        !$OMP& private(rij_sq, x_i, y_i, z_i, x_ij, y_ij, z_ij, &
        !$OMP& rij, inv_denominator, point_potential, n)
        do nz_step = 1,nz
            x_i = cell%origin(1) + dc(1) * real(nz_step - 1, dp)
            y_i = cell%origin(2) + dc(2) * real(nz_step - 1, dp)
            z_i = cell%origin(3) + dc(3) * real(nz_step - 1, dp)

            ! y-loop over grid
            !dir$ vector aligned
            !$OMP SIMD
            do n = 1, ny * nx
                ! x-loop over grid
                ! do nx_step = 1, x_end - x_start + 1
                x_ij = r_j(1) - (x_i + grid_x(n))
                y_ij = r_j(2) - (y_i + grid_y(n))
                z_ij = r_j(3) - (z_i + grid_z(n))

                rij_sq = x_ij * x_ij + y_ij * y_ij + z_ij * z_ij
                rij = sqrt(rij_sq)

                inv_denominator = 1.0_dp / (cov_r_5 - rij_sq * rij_sq * rij)

                ! external potential on a grid point
                point_potential = (cov_r_4 - rij_sq * rij_sq) * inv_denominator

                potential_p(n, nz_step) =  &
                        potential_p(n, nz_step) - point_potential * charge
            enddo !end ny_step
        enddo !end nz_step
        !$OMP END PARALLEL DO
    enddo !end n_atom

    nx_real = (x_end - x_start + 1)
    ! z-loop over grid
    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) DEFAULT(SHARED) &
    !$OMP& private(nz_step, ny_step)
    do nz_step = 1,nz
        potential(:, ny, nz_step) = 0.0_dp
        if (nx /= nx_real) then
            do ny_step = 1, ny
                potential(nx, ny_step, nz_step) = 0.0_dp
            end do
        end if
    enddo !end nz_step
    !$OMP END PARALLEL DO

    call timer_stop

end subroutine compute_sr_potential_electrons