compute_sr_forces_electrons Subroutine

subroutine compute_sr_forces_electrons(atoms, density, cell, &
                                       x_start, x_end, id_start, id_end)
    USE, INTRINSIC :: ISO_C_BINDING

    !> list of atoms
    type(atom_type), dimension(:), intent(in) :: atoms
    !> electronic density
    real(dp), dimension(:,:,:), contiguous, intent(in), target :: density
    !> 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) :: force_x, force_y, force_z, temp_force
    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
    integer, dimension(:), allocatable :: x_ind, y_ind
    real(dp), dimension(:), allocatable :: grid_x, grid_y, grid_z
    real(dp), dimension(:,:), contiguous, pointer :: density_p

    call timer_start("compute_sr_forces_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(density), density_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
        force_x = 0.0_dp
        force_y = 0.0_dp
        force_z = 0.0_dp
        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 REDUCTION(-:force_x, force_y, force_z) 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, temp_force, 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 REDUCTION(-:force_x, force_y, force_z)
            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

                temp_force = density_p(n, nz_step) * &
                        rij_sq * (4.0_dp - 5.0_dp * rij * point_potential) &
                        * inv_denominator

                ! this is not a force yet
                force_x = force_x - x_ij * temp_force

                force_y = force_y - y_ij * temp_force

                force_z = force_z - z_ij * temp_force
            enddo !end n
        enddo !end nz_step
        !$OMP END PARALLEL DO
        atoms(n_atom)%force(1) = atoms(n_atom)%force(1) + force_x * atoms(n_atom)%charge * dV
        atoms(n_atom)%force(2) = atoms(n_atom)%force(2) + force_y * atoms(n_atom)%charge * dV
        atoms(n_atom)%force(3) = atoms(n_atom)%force(3) + force_z * atoms(n_atom)%charge * dV
    enddo !end n_atom

    call timer_stop

end subroutine compute_sr_forces_electrons