compute_sr_forces_ions Subroutine

subroutine compute_sr_forces_ions(atoms, nuclei, &
                                  at_st, at_end)

    !> list of classical atoms
    type(atom_type), dimension(:), intent(inout) :: atoms
    !> list of ions to compute interactions with
    class(nucleus_type), dimension(:), intent(inout) :: nuclei
    !> index of the atom to start electrostatic treatment
    integer, intent(in) :: at_st
    !> index of the atom to finish electrostatic treatment
    integer, intent(in) :: at_end

    integer :: i, j
    integer :: num_atoms
    real(dp), dimension(3) :: r_ij
    real(dp), dimension(3) :: force_temp
    real(dp) :: rij
    real(dp) :: rij_sq
    real(dp) :: rij_4
    real(dp) :: radius_4
    real(dp) :: radius_5
    real(dp) :: nominator
    real(dp) :: denominator
    real(dp) :: inv_denominator
    real(dp), dimension(:), allocatable :: temp_force

    call timer_start("compute_sr_forces_ions")

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

    num_atoms = size(atoms)
    allocate(temp_force(3 * size(nuclei)))
    temp_force = 0.0_dp

    !$OMP PARALLEL DO private(i, radius_4, radius_5, j, r_ij, rij_sq, rij, &
    !$OMP& rij_4, nominator, denominator, inv_denominator, force_temp) &
    !$OMP& reduction(-:temp_force)
    do i = at_st, at_end
        if (abs(atoms(i)%charge) <= 1.0e-12_dp) cycle
        radius_4 = atoms(i)%radius**ORDER
        radius_5 = radius_4 * atoms(i)%radius
        do j = 1, size(nuclei)
                r_ij = atoms(i)%coordinate - nuclei(j)%coordinate
                rij_sq = r_ij(1)*r_ij(1) + r_ij(2)*r_ij(2) + r_ij(3)*r_ij(3)
                rij = sqrt(rij_sq)
                rij_4 = rij_sq * rij_sq

                nominator = radius_4 - rij_4
                denominator = radius_5 - rij_4 * rij

                inv_denominator = 1.0_dp / denominator

                force_temp = r_ij * atoms(i)%charge * nuclei(j)%charge * rij_sq * &
                             (4.0_dp * denominator - 5.0_dp * rij * nominator) * &
                             (inv_denominator * inv_denominator)
                atoms(i)%force = atoms(i)%force + force_temp
                temp_force((j-1)*3 + 1 : j*3) = temp_force((j-1)*3 + 1 : j*3) - force_temp
        enddo
    enddo
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO private(i)
    do i = 1, size(nuclei)
        nuclei(i)%force = nuclei(i)%force + temp_force((i-1)*3 + 1 : i*3)
    end do
    !$OMP END PARALLEL DO

    call timer_stop

end subroutine compute_sr_forces_ions