compute_sr_energy_ions Subroutine

subroutine compute_sr_energy_ions(atoms, nuclei, &
                                  energy, 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
    !> resulting energy of interaction
    real(dp), intent(out) :: energy
    !> 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) :: 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

    call timer_start("compute_sr_energy_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

    energy = 0.0_dp

    num_atoms = size(atoms)

    !$OMP PARALLEL DO private(i, radius_4, radius_5, j, r_ij, rij_sq, rij, &
    !$OMP& rij_4, nominator, denominator, inv_denominator) reduction(+:energy)
    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

            energy = energy + atoms(i)%charge * nuclei(j)%charge * &
                              nominator * inv_denominator
        enddo
    enddo
    !$OMP END PARALLEL DO

    call timer_stop

end subroutine compute_sr_energy_ions