ion_interaction.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module ion_interaction_oct_m
  use atom_oct_m
  use comm_oct_m
  use debug_oct_m
  use global_oct_m
  use jellium_oct_m
  use distributed_oct_m
  use kpoints_oct_m
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_oct_m
  use species_factory_oct_m
 
  implicit none
 
  private
  public ::                           &
    ion_interaction_t,                &
    ion_interaction_init,             &
    ion_interaction_end,              &
    ion_interaction_calculate,        &
    ion_interaction_init_parallelization, &
    ion_interaction_test,             &
    ion_interaction_stress
 
  type ion_interaction_t
    real(real64) :: alpha
    type(distributed_t) :: dist
  end type ion_interaction_t
 
  integer, parameter ::            &
    ION_COMPONENT_REAL     = 1,    &
    ion_component_self     = 2,    &
    ion_component_fourier  = 3,    &
    ion_num_components     = 3
 
contains
 
  subroutine ion_interaction_init(this, namespace, space, natoms)
    type(ion_interaction_t),      intent(out)   :: this
    type(namespace_t),            intent(in)    :: namespace
    class(space_t),               intent(in)    :: space
    integer,                      intent(in)    :: natoms
 
    push_sub(ion_interaction_init)
 
    !%Variable EwaldAlpha
    !%Type float
    !%Default 0.21
    !%Section Hamiltonian
    !%Description
    !% The value 'Alpha' that controls the splitting of the Coulomb
    !% interaction in the Ewald sum used to calculation the ion-ion
    !% interaction for periodic systems. This value affects the speed
    !% of the calculation, normally users do not need to modify it.
    !%End
    call parse_variable(namespace, 'EwaldAlpha', 0.21_real64, this%alpha)
 
    call distributed_init_serial(this%dist, natoms)
 
    if (space%periodic_dim == 1) then
      call messages_experimental('Ewald summation for 1D-periodic systems', namespace=namespace)
    end if
 
    pop_sub(ion_interaction_init)
  end subroutine ion_interaction_init
 
  subroutine ion_interaction_init_parallelization(this, natoms, mc)
    type(ion_interaction_t),      intent(inout) :: this
    integer,                      intent(in)    :: natoms
    type(multicomm_t),            intent(in)    :: mc
 
    push_sub(ion_interaction_init_parallelization)
 
    !As the code below is not parallelized with any of k-point, states nor domain
    !we can safely parallelize it over atoms
    if (debug%info) then
      call distributed_init(this%dist, natoms, mc%master_comm, "Ions")
    else
      call distributed_init(this%dist, natoms, mc%master_comm)
    end if
 
    pop_sub(ion_interaction_init_parallelization)
  end subroutine ion_interaction_init_parallelization
 
  subroutine ion_interaction_end(this)
    type(ion_interaction_t), intent(inout) :: this
 
    push_sub(ion_interaction_end)
 
    this%alpha = -m_one
 
    call distributed_end(this%dist)
 
    pop_sub(ion_interaction_end)
  end subroutine ion_interaction_end
 
  subroutine ion_interaction_calculate(this, space, latt, atom, natoms, pos, lsize, energy, force, &
    energy_components, force_components)
    type(ion_interaction_t), intent(inout) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), contiguous, intent(in) :: lsize(:)
    real(real64), intent(out) :: energy
    real(real64), contiguous, intent(out) :: force(:, :)
    real(real64), contiguous, optional, intent(out) :: energy_components(:)
    real(real64), contiguous, optional, intent(out) :: force_components(:, :, :)
 
 
    push_sub(ion_interaction_calculate)
    call profiling_in("ION_ION_INTERACTION")
 
    if (present(energy_components)) then
      assert(ubound(energy_components, dim = 1) == ion_num_components)
      energy_components = m_zero
    end if
 
    if (present(force_components)) then
      assert(all(ubound(force_components) == (/space%dim, natoms, ion_num_components/)))
      force_components = m_zero
    end if
 
    if (space%is_periodic() .and. any_species_is_jellium_sphere(atom)) then
      call messages_not_implemented('No periodic implementation of ion-ion energy for the jellium sphere')
    end if
 
    if (space%is_periodic()) then
      if (all_species_are_jellium_slab(atom)) then
        energy = jellium_slab_energy_periodic(space, atom, lsize)
        force = 0._real64
      else
        call ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, energy_components, force_components)
      end if
    else
      call ion_interaction_finite(this%dist, space, atom, pos, lsize, energy, force)
      energy = energy + jellium_self_energy_finite(this%dist, latt, atom, lsize)
    end if
 
    call profiling_out("ION_ION_INTERACTION")
    pop_sub(ion_interaction_calculate)
 
  end subroutine ion_interaction_calculate
 
  function jellium_slab_energy_periodic(space, atom, lsize) result(energy)
    class(space_t), intent(in) :: space
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), contiguous, intent(in) :: lsize(:)
    real(real64) :: energy
 
    real(real64) :: area
 
    ! Implementation assumes a single atom
    assert(size(atom) == 1)
    ! This is only allowed if periodic dim = 2. In that case the lattice volume is in fact an area.
    assert(space%periodic_dim == 2)
 
    select type(spec => atom(1)%species)
    type is (jellium_slab_t)
      area = lsize(1) * lsize(2) * m_four
      energy = m_pi * spec%get_density(lsize) **2 * area * spec%thickness()**3 / m_three
    class default
      assert(.false.)
    end select
 
  end function jellium_slab_energy_periodic
 
  function jellium_self_energy_finite(dist, latt, atom, lsize) result(energy)
    type(distributed_t), intent(in) :: dist
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), contiguous, intent(in) :: lsize(:)
    real(real64) :: energy
 
    real(real64)              :: zi
    integer                   :: iatom
    logical                   :: lattice_is_orthogonal
    class(species_t), pointer :: spec
 
    push_sub(jellium_self_energy_finite)
 
    energy = 0._real64
    lattice_is_orthogonal = .not. latt%nonorthogonal
 
    do iatom = dist%start, dist%end
      spec => atom(iatom)%species
      zi = spec%get_zval()
 
      select type(spec)
      type is (jellium_sphere_t)
        energy = energy + (m_three / m_five) * zi**2 / spec%radius()
        ! The part depending on the simulation sphere is neglected
 
      type is (jellium_slab_t)
        ! Jellium slab energy only implemented for orthogonal cells.
        ! One would need to replace (lsize(1) * lsize(2)) * spec%thickness()) with the triple product
        assert(lattice_is_orthogonal)
        energy = energy + m_pi * zi**2 / (m_four * lsize(1)*lsize(2)) * spec%thickness() / m_three
        ! The part depending on the simulation box transverse dimension is neglected
      end select
      nullify(spec)
    enddo
 
    call comm_allreduce(dist%mpi_grp, energy)
 
    pop_sub(jellium_self_energy_finite)
 
  end function jellium_self_energy_finite
 
  subroutine ion_interaction_finite(dist, space, atom, pos, lsize, energy, force)
    type(distributed_t), intent(in) :: dist
    class(space_t), intent(in) :: space
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), contiguous, intent(in) :: pos(:,:)
    real(real64), contiguous, intent(in) :: lsize(:)
    real(real64), intent(out) :: energy
    real(real64), contiguous, intent(out) :: force(:, :)
 
    class(species_t), pointer :: species_i, species_j
    real(real64)              :: r(space%dim), f(space%dim)
    real(real64)              :: r_mag
    real(real64)              :: u_e
    real(real64)              :: zi, zj
    integer                   :: iatom, jatom, natoms
 
    push_sub(ion_interaction_finite)
 
    natoms = size(atom)
    energy = m_zero
    force(1:space%dim, 1:natoms) = m_zero
 
    do iatom = dist%start, dist%end
      species_i => atom(iatom)%species
      zi = species_i%get_zval()
 
      do jatom = iatom + 1, natoms
        species_j => atom(jatom)%species
        zj = species_j%get_zval()
 
        r = pos(:, iatom) - pos(:, jatom)
        r_mag = norm2(r)
        u_e = zi * zj / r_mag
 
        energy = energy + u_e
        f(1:space%dim) = (u_e / r_mag**2) * r(1:space%dim)
        force(1:space%dim, iatom) = force(1:space%dim, iatom) + f(1:space%dim)
        force(1:space%dim, jatom) = force(1:space%dim, jatom) - f(1:space%dim)
      end do
    end do
 
    call comm_allreduce(dist%mpi_grp, energy)
    call comm_allreduce(dist%mpi_grp, force)
 
    nullify(species_i, species_j)
 
    pop_sub(ion_interaction_finite)
 
  end subroutine ion_interaction_finite
 
  subroutine ion_interaction_periodic(this, space, latt, atom, natoms, pos, energy, force, &
    energy_components, force_components)
    type(ion_interaction_t), intent(in) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(out) :: energy
    real(real64), contiguous, intent(out) :: force(:, :)
    real(real64), contiguous, optional, intent(out) :: energy_components(:)
    real(real64), contiguous, optional, intent(out) :: force_components(:, :, :)
 
    real(real64) :: ereal, efourier, epseudo, eself
    real(real64) :: charge
    real(real64), allocatable :: force_real(:,:)
 
    push_sub(ion_interaction_periodic)
 
    energy = m_zero
    force(1:space%dim, 1:natoms) = m_zero
 
    call ewald_short(this%dist, space, latt, atom, pos, this%alpha, ereal, force)
    if (present(force_components)) then
      force_components(1:space%dim, 1:natoms, ion_component_real) = force(1:space%dim, 1:natoms)
    end if
 
    call ewald_self_interaction(this%dist, atom, this%alpha, eself, charge)
 
    call profiling_in("EWALD_LONG")
    select case (space%periodic_dim)
    case (1)
      ! The 1D-periodic Ewald sum uses a soft-Coulomb range separation (see ewald_1D),
      ! different from the shared erfc Ewald_short/Ewald_self_interaction computed above.
      ! Recompute all components and discard the erfc real-space contribution.
      ! The long range part of the pseudopotential is not accounted for here (same situation as the 2D case). See Issue #950.
      epseudo = m_zero
      force(1:space%dim, 1:natoms) = m_zero
      safe_allocate(force_real(1:space%dim, 1:natoms))
      call ewald_1d(this, space, latt, atom, natoms, pos, ereal, eself, efourier, force, force_real)
      if (present(force_components)) then
        force_components(1:space%dim, 1:natoms, ion_component_real) = force_real(1:space%dim, 1:natoms)
      end if
      safe_deallocate_a(force_real)
    case (2)
      ! The energy contribution of the long range part of the pseudo is
      ! not correctly accounted for in systems periodic in 1D or 2D, however
      ! this term should not appear here anyway. See Issue #950.
      epseudo = m_zero
      call ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
    case (3)
      call ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
      !TODO(Alex/Nicolas) Issue #950. Refactor: Move G=0 correction from ion-ion energy to pseudopotential energy
      call pseudopotential_correction_3d(this%dist, latt, atom, charge, epseudo)
    end select
    call profiling_out("EWALD_LONG")
 
    if (present(energy_components)) then
      energy_components(ion_component_real) = ereal
      energy_components(ion_component_self) = eself
      energy_components(ion_component_fourier) = efourier
    end if
 
    if (present(force_components)) then
      ! This is dependent on the order in which the force terms are computed
      force_components(1:space%dim, 1:natoms, ion_component_fourier) = &
        force(1:space%dim, 1:natoms) - force_components(1:space%dim, 1:natoms, ion_component_real)
    end if
 
    energy = ereal + efourier + eself + epseudo
 
    pop_sub(ion_interaction_periodic)
  end subroutine ion_interaction_periodic
 
  subroutine ewald_short(dist, space, latt, atom, pos, alpha, ereal, force)
    type(distributed_t), intent(in) :: dist
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), contiguous, intent(in) :: pos(:, :)
 
    real(real64), intent(in) :: alpha
    real(real64), intent(out) :: ereal
    real(real64), contiguous, intent(inout) :: force(:, :)
    !                                                            Intent(inout) allows force contributions to be summed
    integer                  :: iatom, jatom, icopy, natoms
    real(real64)             :: rnorm, xi(space%dim)
    real(real64)             :: force_real(space%dim)
    real(real64)             :: zi, zj
    real(real64)             :: erfc_rnorm
    real(real64)             :: rcut
    type(lattice_iterator_t) :: latt_iter
    real(real64)             :: charge, coeff
 
    push_sub_with_profile(ewald_short)
 
    ereal = m_zero
    ! Origin of this value is not documented
    rcut = 6.0_real64 / alpha
    latt_iter = lattice_iterator_t(latt, rcut)
    natoms = size(atom)
 
    charge = m_zero
    do iatom = dist%start, dist%end
      if (.not. atom(iatom)%species%represents_real_atom()) cycle
      zi = atom(iatom)%species%get_zval()
      charge = charge + zi**2
    end do
 
    ! Diagonal terms iatom == jatom for all cells, except T=(0,0,0)
    ! Note: Only half of the copies are needed, by symmetries
    do icopy = 1, latt_iter%n_cells
      rnorm = norm2(latt_iter%get(icopy))
      if (rnorm < r_min_atom_dist) cycle
      if (rnorm > rcut) cycle
      ereal = ereal + m_half * charge * erfc(alpha * rnorm) /rnorm
    end do
 
    coeff = m_two * alpha / sqrt(m_pi)
 
    !$omp parallel default(shared) &
    !$omp& private(iatom, jatom, zi, zj, icopy, xi, rnorm, erfc_rnorm, force_real, charge) &
    !$omp& reduction(+: ereal, force)
    do iatom = dist%start, dist%end
      if (.not. atom(iatom)%species%represents_real_atom()) cycle
      zi = atom(iatom)%species%get_zval()
 
      ! Upper triangle, for all replica cells
      do jatom = iatom + 1, natoms
        zj = atom(jatom)%species%get_zval()
 
        charge = zi*zj
 
        !$omp do
        do icopy = 1, latt_iter%n_cells
          xi = pos(:, iatom) + latt_iter%get(icopy)
          rnorm = norm2(xi - pos(:, jatom))
          if (rnorm > rcut) cycle
 
          erfc_rnorm = erfc(alpha * rnorm) / rnorm
 
          ! Factor 1/2 omitted as one is only summing over upper triangle
          ereal = ereal + charge * erfc_rnorm
 
          force_real(:) = charge * (xi - pos(:, jatom)) * &
            (erfc_rnorm + coeff *exp(-(alpha*rnorm)**2)) / rnorm**2
 
          ! Upper trianglar contribution
          force(1:space%dim, jatom) = force(1:space%dim, jatom) - force_real
 
          ! Lower triangular contribution
          force(1:space%dim, iatom) = force(1:space%dim, iatom) + force_real
        end do
        !$omp end do
 
      end do
    end do
    !$omp end parallel
 
    call comm_allreduce(dist%mpi_grp, ereal)
    call comm_allreduce(dist%mpi_grp, force)
 
    pop_sub_with_profile(ewald_short)
  end subroutine ewald_short
 
  subroutine ewald_self_interaction(dist, atom, alpha, eself, charge)
    type(distributed_t), intent(in) :: dist
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), intent(in) :: alpha
    real(real64), intent(out) :: eself
    real(real64), intent(out) :: charge
 
    integer                  :: iatom
    real(real64)             :: zi
 
    push_sub(ewald_self_interaction)
 
    eself = m_zero
    charge = m_zero
 
    do iatom = dist%start, dist%end
      zi = atom(iatom)%species%get_zval()
      charge = charge + zi
      eself = eself - alpha / sqrt(m_pi) * zi**2
    end do
 
    call comm_allreduce(dist%mpi_grp, eself)
    call comm_allreduce(dist%mpi_grp, charge)
 
    pop_sub(ewald_self_interaction)
  end subroutine ewald_self_interaction
 
  subroutine ewald_long_3d(this, space, latt, atom, natoms, pos, efourier, force, charge)
    type(ion_interaction_t), intent(in) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), contiguous, intent(in) :: pos(:,:)
    real(real64), intent(inout) :: efourier
    real(real64), contiguous, intent(inout) :: force(:, :)
    real(real64), intent(in) :: charge
 
    real(real64) :: rcut, gmax_squared
    integer :: iatom
    integer :: ix, iy, iz, isph
    real(real64)   :: gvec(3), gred(3), gg2, gx
    real(real64)   :: factor
    complex(real64)   :: sumatoms, tmp(3), aa
 
    complex(real64), allocatable :: phase(:)
 
    push_sub(ewald_long_3d)
 
    assert(space%dim == 3)
    assert(space%periodic_dim == 3)
 
    ! And the long-range part, using an Ewald sum
    safe_allocate(phase(1:natoms))
 
    ! get a converged value for the cutoff in g
    rcut = sqrt(minval(sum(latt%klattice**2, dim=1)))
 
    ! 9.5 is a constant that controls the range separation
    isph = ceiling(9.5_real64*this%alpha/rcut)
 
    ! First the G = 0 term (charge was calculated previously)
    efourier = -m_pi*charge**2/(m_two*this%alpha**2*latt%rcell_volume)
 
    ! Getting the G-shell cutoff
    gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
 
    do ix = -isph, isph
      do iy = -isph, isph
        do iz = -isph, isph
 
          ! Exploit k <-> -k symmetry
          ! Only process one half of reciprocal space.
          ! g=0 must also be removed from the sum
          if (ix < 0) cycle
          if (ix == 0 .and. iy < 0) cycle
          if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
 
          gred = [ix, iy, iz]
          call kpoints_to_absolute(latt, gred, gvec)
          gg2 = dot_product(gvec, gvec)
 
          if (gg2 > gmax_squared*1.001_real64) cycle
 
          gx = -0.25_real64*gg2/this%alpha**2
 
          if (gx < -36.0_real64) cycle
 
          ! We have used the k-> -k symmetry, hence the factor 4
          factor = m_four*m_pi/latt%rcell_volume*exp(gx)/gg2
 
          if (factor < epsilon(factor)) cycle
 
          sumatoms = m_z0
          !$omp parallel do private(iatom, gx, aa) reduction(+:sumatoms)
          do iatom = 1, natoms
            gx = sum(gvec*pos(:,iatom))
            aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
            phase(iatom) = aa
            sumatoms = sumatoms + aa
          end do
 
          efourier = efourier + factor * real(sumatoms*conjg(sumatoms), real64)
 
          do iatom = 1, natoms
            tmp = m_zi*gvec*phase(iatom)
            force(1:space%dim, iatom) = force(1:space%dim, iatom) - factor*real(conjg(tmp)*sumatoms + tmp*conjg(sumatoms), real64)
 
          end do
 
        end do
      end do
    end do
 
    safe_deallocate_a(phase)
 
    pop_sub(ewald_long_3d)
 
  end subroutine ewald_long_3d
 
  subroutine ewald_long_2d(this, space, latt, atom, natoms, pos, efourier, force)
    type(ion_interaction_t), intent(in) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(inout) :: efourier
    real(real64), contiguous, intent(inout) :: force(:, :)
 
    real(real64) :: rcut, gmax_squared
    integer :: iatom, jatom
    integer :: ix, iy, ix_max, iy_max
    real(real64)   :: gvec(space%dim), gg2, gx, gg_abs
    real(real64)   :: factor,factor1,factor2, coeff
    real(real64)   :: dz_max, dz_ij, erfc1, erfc2, tmp_erf
    real(real64), allocatable :: force_tmp(:,:)
    real(real64), parameter :: tol = 1e-10_real64
 
    push_sub(ewald_long_2d)
 
    assert(space%periodic_dim == 2)
    assert(space%dim == 2 .or. space%dim == 3)
 
    ! And the long-range part, using an Ewald sum
 
    ! Searching maximum distance
    if (space%dim == 3) then
      dz_max = m_zero
      do iatom = 1, natoms
        do jatom = iatom + 1, natoms
          dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
        end do
      end do
    else
      ! For a 2D system, all atoms are on the plane, so the distance is zero
      dz_max = m_zero
    end if
 
    !get a converged value for the cutoff in g
    rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
    if (dz_max > tol) then
      do
        if (rcut * dz_max >= m_max_exp_arg) exit  !Maximum double precision number
        erfc1 = erfc(this%alpha*dz_max + m_half*rcut/this%alpha)
        if (erfc1 * exp(rcut*dz_max) < 1.e-10_real64) exit
        rcut = rcut * 1.414_real64
      end do
    end if
 
    ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
    iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
 
    safe_allocate(force_tmp(1:space%dim, 1:natoms))
    force_tmp = m_zero
 
    ! First the G = 0 term
    efourier = m_zero
    factor = m_pi/latt%rcell_volume
    !$omp parallel do private(jatom, dz_ij, tmp_erf, factor1, factor2) reduction(+:efourier,force_tmp) &
    !$omp& collapse(2)
    do iatom = this%dist%start, this%dist%end
      do jatom = 1, natoms
        ! efourier
        if (space%dim == 3) then
          dz_ij = pos(3, iatom) - pos(3, jatom)
        else
          dz_ij = m_zero
        end if
 
        tmp_erf = erf(this%alpha*dz_ij)
        factor1 = dz_ij*tmp_erf
        factor2 = exp(-(this%alpha*dz_ij)**2)/(this%alpha*sqrt(m_pi))
 
        efourier = efourier - factor &
          * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
 
        ! force
        if (iatom == jatom)cycle
        if (abs(tmp_erf) < m_epsilon) cycle
 
        if (space%dim == 3) then
          force_tmp(3, iatom) = force_tmp(3, iatom) - (- m_two*factor) &
            * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * tmp_erf
        end if
 
      end do
    end do
 
    ! Getting the G-shell cutoff
    gmax_squared = ix_max**2 * sum(latt%klattice(:, 1)**2)
    gmax_squared = min(gmax_squared, iy_max**2 * sum(latt%klattice(:, 2)**2))
 
    !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, gx, dz_ij, erfc1, factor1, erfc2, factor2, coeff) &
    !$omp& collapse(2) reduction(+:efourier, force_tmp)
    do ix = -ix_max, ix_max
      do iy = -iy_max, iy_max
 
        gvec = ix*latt%klattice(:, 1) + iy*latt%klattice(:, 2)
        gg2 = sum(gvec**2)
 
        ! g=0 must be removed from the sum
        if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
        gg_abs = sqrt(gg2)
        factor = m_half*m_pi/(latt%rcell_volume*gg_abs)
 
        do iatom = this%dist%start, this%dist%end
          do jatom = iatom, natoms
            ! efourier
            gx = sum(gvec(1:2) * (pos(1:2, iatom) - pos(1:2, jatom)))
            gx = gvec(1)*(pos(1, iatom) - pos(1, jatom)) + gvec(2)*(pos(2, iatom) - pos(2, jatom))
            if (space%dim == 3) then
              dz_ij = pos(3, iatom) - pos(3, jatom)
            else
              dz_ij = m_zero
            end if
 
            erfc1 = erfc(this%alpha*dz_ij + m_half*gg_abs/this%alpha)
            if (abs(erfc1) > m_epsilon) then
              factor1 = exp(gg_abs*dz_ij)*erfc1
            else
              factor1 = m_zero
            end if
            erfc2 = erfc(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)
            if (abs(erfc2) > m_epsilon) then
              factor2 = exp(-gg_abs*dz_ij)*erfc2
            else
              factor2 = m_zero
            end if
 
            if (iatom == jatom) then
              coeff = m_one
            else
              coeff = m_two
            end if
 
            efourier = efourier &
              + factor * coeff &
              * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
              * cos(gx)* ( factor1 + factor2)
 
            ! force
            if (iatom == jatom) cycle
 
            force_tmp(1:2, iatom) = force_tmp(1:2, iatom) &
              + m_two * factor * gvec(1:2) &
              * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
              *sin(gx)*(factor1 + factor2)
 
            force_tmp(1:2, jatom) = force_tmp(1:2, jatom) &
              - m_two * factor * gvec(1:2) &
              * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
              *sin(gx)*(factor1 + factor2)
 
            factor1 = gg_abs*erfc1 &
              - m_two*this%alpha/sqrt(m_pi)*exp(-(this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
            if (abs(factor1) > m_epsilon) then
              factor1 = factor1*exp(gg_abs*dz_ij)
            else
              factor1 = m_zero
            end if
 
            factor2 = gg_abs*erfc2 &
              - m_two*this%alpha/sqrt(m_pi)*exp(-(-this%alpha*dz_ij + m_half*gg_abs/this%alpha)**2)
            if (abs(factor2) > m_epsilon) then
              factor2 = factor2*exp(-gg_abs*dz_ij)
            else
              factor2 = m_zero
            end if
 
            if (space%dim == 3) then
              force_tmp(3, iatom) = force_tmp(3, iatom) &
                - m_two*factor &
                * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
                * cos(gx)* ( factor1 - factor2)
              force_tmp(3, jatom) = force_tmp(3, jatom) &
                + m_two*factor &
                * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() &
                * cos(gx)* ( factor1 - factor2)
            end if
 
          end do
        end do
 
 
      end do
    end do
 
    call comm_allreduce(this%dist%mpi_grp, efourier)
    call comm_allreduce(this%dist%mpi_grp, force_tmp)
 
    force = force + force_tmp
 
    safe_deallocate_a(force_tmp)
 
    pop_sub(ewald_long_2d)
  end subroutine ewald_long_2d
 
  subroutine ewald_1d(this, space, latt, atom, natoms, pos, ereal, eself, efourier, force, force_real)
    type(ion_interaction_t),   intent(in)    :: this
    class(space_t),            intent(in)    :: space
    type(lattice_vectors_t),   intent(in)    :: latt
    type(atom_t),              intent(in)    :: atom(:)
    integer,                   intent(in)    :: natoms
    real(real64),              intent(in)    :: pos(1:space%dim,1:natoms)
    real(real64),              intent(out)   :: ereal
    real(real64),              intent(out)   :: eself
    real(real64),              intent(out)   :: efourier
    real(real64),              intent(inout) :: force(:, :)
    real(real64),              intent(out)   :: force_real(:, :)
 
    real(real64), parameter :: tol = 1e-10_real64
    ! Internal soft-Coulomb splitting length. The total energy and forces are
    ! independent of its value; it only controls the real/reciprocal convergence rate.
    ! A value of 1 ensures most of the work is performed by the short-range sum (while L is small),
    ! which is cheaper to evaluate than the Bessel functions of the long-range sum.
    real(real64), parameter :: sigma = m_one
    integer :: iatom, jatom, n, n_max, icopy
    real(real64) :: length, gk, rcut
    real(real64) :: uaxis(space%dim), rij(space%dim), rperp(space%dim), xi(space%dim), gvec(space%dim)
    real(real64) :: ri(space%dim), rj(space%dim)
    real(real64) :: xpar, rho, rbar, zi, zj, zz, rnorm, rt2s, phase, gmag
    real(real64) :: dcoeff, k0v, k1v, ctrans, calong
    real(real64) :: fvec(space%dim)
    real(real64), allocatable :: force_recip(:,:)
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(ewald_1d)
 
    assert(space%periodic_dim == 1)
    assert(space%dim >= 1 .and. space%dim <= 3)
 
    length = norm2(latt%rlattice(1:space%dim, 1))
    uaxis(1:space%dim) = latt%rlattice(1:space%dim, 1) / length
    ! Reciprocal vectors are taken along the periodic axis, g_n = (2 pi n / L) uaxis.
    ! This does not rely on klattice(:,1), which is parallel to uaxis only if the
    ! non-periodic lattice vectors are orthogonal to the periodic one.
    gk = m_two * m_pi / length
 
    ! Self-energy correction (remove soft-Coulomb self)
    eself = m_zero
    do iatom = this%dist%start, this%dist%end
      if (atom(iatom)%species%is_user_defined()) cycle
      eself = eself - m_half * atom(iatom)%species%get_zval()**2 / sigma
    end do
    call comm_allreduce(this%dist%mpi_grp, eself)
 
    ! Short-range real-space sum
    ! The summand decays as sigma^2/(2 r^3); choose the cutoff so the 1/t^3 tail < tol.
    rcut = sqrt(sigma**2 / (m_two * length * tol))
    latt_iter = lattice_iterator_t(latt, rcut)
 
    ereal = m_zero
    force_real(1:space%dim, 1:natoms) = m_zero
    !$omp parallel do default(shared) &
    !$omp& private(iatom, zi, ri, jatom, zj, zz, rj, icopy, xi, rnorm, rt2s, dcoeff, fvec) &
    !$omp& reduction(+:ereal, force_real)
    do iatom = this%dist%start, this%dist%end
      if (atom(iatom)%species%is_user_defined()) cycle
      zi = atom(iatom)%species%get_zval()
      ri = pos(:, iatom)
      ! Diagonal self-images (i==i, t/=0). The net force on an atom from its own,
      ! symmetric images vanishes, so no force contribution here.
      do icopy = 1, latt_iter%n_cells
        rnorm = norm2(latt_iter%get(icopy))
        if (rnorm < r_min_atom_dist) cycle
        if (rnorm > rcut) cycle
        ereal = ereal + m_half * zi*zi * (m_one/rnorm - m_one/sqrt(rnorm**2 + sigma**2))
      end do
      ! Upper-triangle pairs (i<j), summed over all images (both orderings folded in).
      do jatom = iatom + 1, natoms
        if (atom(jatom)%species%is_user_defined()) cycle
        zj = atom(jatom)%species%get_zval()
        zz = zi*zj
        rj = pos(:, jatom)
        do icopy = 1, latt_iter%n_cells
          xi = ri + latt_iter%get(icopy)
          rnorm = norm2(xi - rj)
          if (rnorm < r_min_atom_dist) cycle
          if (rnorm > rcut) cycle
          rt2s = sqrt(rnorm**2 + sigma**2)
          ereal = ereal + zz * (m_one/rnorm - m_one/rt2s)
          dcoeff = m_one/rnorm**3 - m_one/rt2s**3
          fvec(1:space%dim) = zz * (xi(1:space%dim) - rj(1:space%dim)) * dcoeff
          force_real(1:space%dim, iatom) = force_real(1:space%dim, iatom) + fvec(1:space%dim)
          force_real(1:space%dim, jatom) = force_real(1:space%dim, jatom) - fvec(1:space%dim)
        end do
      end do
    end do
    !$omp end parallel do
    call comm_allreduce(this%dist%mpi_grp, ereal)
    call comm_allreduce(this%dist%mpi_grp, force_real)
 
    ! Long-range reciprocal sum (full K_0) + k_z=0 term
    ! K_0(rbar*g_n) ~ exp(-rbar*g_n) with rbar >= sigma, so bound the number of shells.
    n_max = max(1, ceiling(-log(tol)/(sigma*gk)) + 1)
 
    safe_allocate(force_recip(1:space%dim, 1:natoms))
    force_recip(1:space%dim, 1:natoms) = m_zero
    efourier = m_zero
    !$omp parallel do default(shared) &
    !$omp& private(iatom, zi, jatom, zj, zz, rij, xpar, rperp, rho, rbar, &
    !$omp&         n, gmag, gvec, phase, k0v, k1v, ctrans, calong, fvec) &
    !$omp& reduction(+:efourier, force_recip)
    do iatom = this%dist%start, this%dist%end
      if (atom(iatom)%species%is_user_defined()) cycle
      zi = atom(iatom)%species%get_zval()
      do jatom = iatom, natoms
        if (atom(jatom)%species%is_user_defined()) cycle
        zj = atom(jatom)%species%get_zval()
        zz = zi*zj
        rij(1:space%dim) = pos(1:space%dim, iatom) - pos(1:space%dim, jatom)
        xpar = dot_product(rij(1:space%dim), uaxis(1:space%dim))
        rperp(1:space%dim) = rij(1:space%dim) - xpar*uaxis(1:space%dim)
        rho = norm2(rperp(1:space%dim))
        rbar = sqrt(rho**2 + sigma**2)   ! always > 0, so K_0 is never singular
 
        if (iatom == jatom) then
          ! Self term (rho=0, z_ij=0): k_z=0 log term; reciprocal self; no force.
          ! Here rbar = sigma, so the log term vanishes for sigma = 1; kept for generality.
          efourier = efourier - (m_one/length) * zz * log(rbar)
          do n = 1, n_max
            gmag = n*gk
            efourier = efourier + (m_two/length) * zz * loct_bessel_k0(rbar*gmag)
          end do
        else
          ! Ordered pair {i,j} counted once here with a factor of 2 (both orderings).
          ! k_z=0 log term and its transverse force.
          efourier = efourier - (m_two/length) * zz * log(rbar)
          fvec(1:space%dim) = (m_two/length) * zz * rperp(1:space%dim) / rbar**2
          force_recip(1:space%dim, iatom) = force_recip(1:space%dim, iatom) + fvec(1:space%dim)
          force_recip(1:space%dim, jatom) = force_recip(1:space%dim, jatom) - fvec(1:space%dim)
          ! k_z /= 0 terms.
          do n = 1, n_max
            gmag = n*gk
            gvec(1:space%dim) = gmag*uaxis(1:space%dim)
            phase = gmag*xpar
            k0v = loct_bessel_k0(rbar*gmag)
            k1v = loct_bessel_k1(rbar*gmag)
            efourier = efourier + (m_four/length) * zz * k0v * cos(phase)
            ! transverse force (through rbar) and along-axis force (through the phase)
            ctrans = (m_four/length) * zz * gmag * k1v * cos(phase) / rbar
            calong = (m_four/length) * zz * k0v * sin(phase)
            fvec(1:space%dim) = ctrans * rperp(1:space%dim) + calong * gvec(1:space%dim)
            force_recip(1:space%dim, iatom) = force_recip(1:space%dim, iatom) + fvec(1:space%dim)
            force_recip(1:space%dim, jatom) = force_recip(1:space%dim, jatom) - fvec(1:space%dim)
          end do
        end if
      end do
    end do
    !$omp end parallel do
    call comm_allreduce(this%dist%mpi_grp, efourier)
    call comm_allreduce(this%dist%mpi_grp, force_recip)
 
    force(1:space%dim, 1:natoms) = force(1:space%dim, 1:natoms) &
      + force_real(1:space%dim, 1:natoms) + force_recip(1:space%dim, 1:natoms)
 
    safe_deallocate_a(force_recip)
 
    pop_sub(ewald_1d)
  end subroutine ewald_1d
 
  !TODO(Alex/Nicolas) Issue #950. Refactor: Move G=0 correction from ion-ion energy to pseudopotential energy
  subroutine pseudopotential_correction_3d(dist, latt, atom, charge, epseudo)
    type(distributed_t), intent(in) :: dist
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    real(real64), intent(in) :: charge
    real(real64), intent(out) :: epseudo
 
    real(real64)             :: zi
    integer                  :: iatom
 
    push_sub(pseudopotential_correction_3d)
 
    epseudo = m_zero
    do iatom = dist%start, dist%end
      select type(spec => atom(iatom)%species)
      class is(pseudopotential_t)
        zi = spec%get_zval()
        epseudo = epseudo + m_pi *zi * &
          (spec%ps%sigma_erf * sqrt(m_two))**2 / latt%rcell_volume * charge
      end select
    end do
    call comm_allreduce(dist%mpi_grp, epseudo)
 
    pop_sub(pseudopotential_correction_3d)
 
  end subroutine pseudopotential_correction_3d
 
  subroutine ion_interaction_stress(this, space, latt, atom, natoms, pos, stress_ii)
    type(ion_interaction_t), intent(inout) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(out) :: stress_ii(space%dim, space%dim)
 
    real(real64) :: stress_short(1:space%dim, 1:space%dim), stress_Ewald(1:space%dim, 1:space%dim)
 
    push_sub(ion_interaction_stress)
 
    stress_ii = m_zero
 
    ! Only implemented in the periodic case
    assert(space%is_periodic())
 
    ! Short range part in real space
    call ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
 
    ! Long range part in Fourier space
    select case(space%periodic_dim)
    case(3)
      call ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
    case(2)
      call ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_ewald)
    case default
      assert(.false.)
    end select
 
    stress_ii = stress_short + stress_ewald
 
    pop_sub(ion_interaction_stress)
  end subroutine ion_interaction_stress
 
  ! ---------------------------------------------------------
 
  subroutine ion_interaction_stress_short(this, space, latt, atom, natoms, pos, stress_short)
    type(ion_interaction_t), intent(inout) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(out) :: stress_short(1:space%dim, 1:space%dim)
 
    real(real64) :: xi(space%dim)
    real(real64) :: r_ij, zi, zj, Hp, factor
    integer :: iatom, jatom, icopy, idir, jdir
    real(real64) :: alpha, rcut
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(ion_interaction_stress_short)
    call profiling_in("ION_ION_STRESS_SHORT")
 
    ! Only implemented in the periodic case
    assert(space%is_periodic())
 
    alpha = this%alpha
 
    ! See the code for the energy above to understand this parameter
    rcut = 6.0_real64/alpha
 
    ! the short-range part is calculated directly
    stress_short = m_zero
    latt_iter = lattice_iterator_t(latt, rcut)
 
    do iatom = this%dist%start, this%dist%end
      select type(spec => atom(iatom)%species)
      class is(jellium_t)
        cycle
      end select
      zi = atom(iatom)%species%get_zval()
 
      do icopy = 1, latt_iter%n_cells
        xi = pos(:, iatom) + latt_iter%get(icopy)
 
        do jatom = 1, natoms
          zj = atom(jatom)%species%get_zval()
          r_ij = norm2(xi - pos(:, jatom))
 
          if (r_ij < r_min_atom_dist) cycle
 
          hp = -m_two/sqrt(m_pi)*exp(-(alpha*r_ij)**2) - erfc(alpha*r_ij)/(alpha*r_ij)
          factor = m_half*zj*zi*alpha*hp
          do idir = 1, space%periodic_dim
            do jdir = 1, space%periodic_dim
              stress_short(idir, jdir) = stress_short(idir, jdir) &
                - factor*(xi(idir) - pos(idir, jatom))*(xi(jdir) - pos(jdir, jatom))/(r_ij**2)
            end do
          end do
 
        end do
      end do
    end do
 
    if (this%dist%parallel) then
      call comm_allreduce(this%dist%mpi_grp, stress_short)
    end if
 
    stress_short = stress_short/latt%rcell_volume
 
    call profiling_out("ION_ION_STRESS_SHORT")
 
    pop_sub(ion_interaction_stress_short)
  end subroutine ion_interaction_stress_short
 
 
 
  ! ---------------------------------------------------------
  subroutine ewald_3d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
    type(ion_interaction_t), intent(inout) :: this
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(out) :: stress_ewald(3, 3)
 
    real(real64)   :: zi, rcut, gmax_squared
    integer :: iatom
    integer :: ix, iy, iz, isph, idim, idir, jdir
    real(real64)   :: gred(3), gvec(3), gg2, gx
    real(real64)   :: factor, charge, charge_sq, off_diagonal_weight
    complex(real64)   :: sumatoms, aa
 
    call profiling_in("STRESS_3D_EWALD")
    push_sub(ewald_3d_stress)
 
    ! Currently this is only implemented for 3D
    assert(space%dim == 3)
    assert(space%periodic_dim == 3) ! Not working for mixed periodicity
    !                                 (klattice along the non-periodic directions is wrong)
    !                                 Anyway gg/gg2 is not correct for mixed periodicity
 
    stress_ewald = m_zero
 
    ! And the long-range part, using an Ewald sum
    charge = m_zero
    charge_sq = m_zero
    do iatom = 1, natoms
      zi = atom(iatom)%species%get_zval()
      charge = charge + zi
      charge_sq = charge_sq + zi**2
    end do
 
    ! get a converged value for the cutoff in g
    rcut = huge(rcut)
    do idim = 1, space%periodic_dim
      rcut = min(rcut, sum(latt%klattice(1:space%periodic_dim, idim)**2))
    end do
 
    rcut = sqrt(rcut)
 
    isph = ceiling(9.5_real64*this%alpha/rcut)
 
    ! Getting the G-shell cutoff
    gmax_squared = isph**2 * minval(sum(latt%klattice**2, dim=1))
 
    do ix = -isph, isph
      do iy = -isph, isph
        do iz = -isph, isph
 
          ! Exploit k <-> -k symmetry
          ! Only process one half of reciprocal space.
          ! g=0 must also be removed from the sum
          if (ix < 0) cycle
          if (ix == 0 .and. iy < 0) cycle
          if (ix == 0 .and. iy == 0 .and. iz <= 0) cycle
 
          gred = [ix, iy, iz]
          call kpoints_to_absolute(latt, gred, gvec)
          gg2 = sum(gvec**2)
 
          ! g=0 must be removed from the sum
          if (gg2 > gmax_squared*1.001_real64) cycle
 
          gx = -0.25_real64*gg2/this%alpha**2
 
          if (gx < -36.0_real64) cycle
 
          ! We have used the k-> -k symmetry, hence the factor 4
          factor = m_four*m_pi*exp(gx)/(latt%rcell_volume*gg2)
 
          if (factor < epsilon(factor)) cycle
 
          sumatoms = m_z0
 
          do iatom = 1, natoms
            gx = sum(gvec*pos(:, iatom))
            aa = atom(iatom)%species%get_zval()*cmplx(cos(gx), sin(gx), real64)
            sumatoms = sumatoms + aa
          end do
 
          factor = factor*abs(sumatoms)**2
          off_diagonal_weight = - m_two*factor/gg2*(0.25_real64*gg2/this%alpha**2+m_one)
 
          do idir = 1, 3
            do jdir = 1, 3
              stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
                + gvec(idir) * gvec(jdir) * off_diagonal_weight
            end do
            stress_ewald(idir, idir) = stress_ewald(idir, idir) + factor
          end do
 
        end do
      end do
    end do
 
 
    ! The G = 0 term of the Ewald summation
    factor = m_half*m_pi*charge**2/(latt%rcell_volume*this%alpha**2)
    do idir = 1,3
      stress_ewald(idir,idir) = stress_ewald(idir,idir) - factor
    end do
 
    stress_ewald = stress_ewald / latt%rcell_volume
 
 
    call profiling_out("STRESS_3D_EWALD")
    pop_sub(ewald_3d_stress)
 
  end subroutine ewald_3d_stress
 
  ! ---------------------------------------------------------
  subroutine ewald_2d_stress(this, space, latt, atom, natoms, pos, stress_Ewald)
    type(ion_interaction_t), intent(inout) :: this
    type(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), intent(out) :: stress_Ewald(3, 3)
 
    real(real64) :: rcut, efourier
    integer :: iatom, jatom, idir, jdir
    integer :: ix, iy, ix_max, iy_max
    real(real64)   :: gvec(3), gred(3), gg2, cos_gx, gg_abs, gmax_squared
    real(real64)   :: factor,factor1,factor2, coeff, e_ewald
    real(real64)   :: dz_max, z_ij, erfc1, erfc2, diff(3)
    real(real64), parameter   :: tol = 1e-10_real64
 
    push_sub(ewald_2d_stress)
 
    assert(space%periodic_dim == 2)
    assert(space%dim == 3)
 
    stress_ewald = m_zero
 
    ! Searching maximum distance
    dz_max = m_zero
    do iatom = 1, natoms
      do jatom = iatom + 1, natoms
        dz_max = max(dz_max, abs(pos(3, iatom) - pos(3, jatom)))
      end do
    end do
 
    !get a converged value for the cutoff in g
    ! Note: to understand these numbers, one needs to look into the energy routine for Ewald 2D
    rcut = m_two*this%alpha*4.6_real64 + m_two*this%alpha**2*dz_max
    if (dz_max > tol) then ! Else the code here does not work properly
      do
        if (rcut * dz_max >= m_max_exp_arg) exit  !Maximum double precision number
        erfc1 = erfc(this%alpha*dz_max + m_half*rcut/this%alpha)
        if (erfc1 * exp(rcut*dz_max) < tol) exit
        rcut = rcut * 1.414_real64
      end do
    end if
 
    ! First the G = 0 term
    efourier = m_zero
    factor = m_pi/latt%rcell_volume
    !$omp parallel do private(jatom, z_ij, factor1, factor2) reduction(+:efourier) collapse(2)
    do iatom = 1, natoms
      do jatom = 1, natoms
        z_ij = pos(3, iatom) - pos(3, jatom)
 
        factor1 = z_ij * erf(this%alpha*z_ij)
        factor2 = exp(-(this%alpha*z_ij)**2)/(this%alpha*sqrt(m_pi))
 
        efourier = efourier - factor &
          * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval() * (factor1 + factor2)
      end do
    end do
 
    ! Adding the G=0 term
    do idir = 1, 2
      stress_ewald(idir, idir) = efourier
    end do
 
    ! Getting the G-shell cutoff
    ix_max = ceiling(rcut/norm2(latt%klattice(:, 1)))
    iy_max = ceiling(rcut/norm2(latt%klattice(:, 2)))
    gmax_squared = ix_max**2 * sum(latt%klattice(:, 1)**2)
    gmax_squared = min(gmax_squared, iy_max**2 * sum(latt%klattice(:, 2)**2))
 
    !$omp parallel do private(iy, gvec, gg2, gg_abs, factor, iatom, jatom, diff, cos_gx, z_ij, idir, jdir, erfc1, factor1) &
    !$omp& private(erfc2, factor2, coeff, e_ewald, gred) &
    !$omp& collapse(2) reduction(+:stress_Ewald)
    do ix = -ix_max, ix_max
      do iy = -iy_max, iy_max
 
        gred = [ix, iy, 0]
        call kpoints_to_absolute(latt, gred, gvec)
        gg2 = dot_product(gvec,gvec)
 
        ! g=0 must be removed from the sum
        if (gg2 < m_epsilon .or. gg2 > gmax_squared*1.001_real64) cycle
 
        gg_abs = sqrt(gg2)
        factor = m_fourth*m_pi/(latt%rcell_volume*this%alpha*gg2)
 
        do iatom = 1, natoms
          do jatom = iatom, natoms
            diff = pos(:, iatom) - pos(:, jatom)
            cos_gx = cos(sum(gvec(1:2) * diff(1:2)))
            z_ij = diff(3)
 
            factor1 = screening_function_2d(this%alpha, z_ij, gg_abs, erfc1)
            factor2 = screening_function_2d(this%alpha,-z_ij, gg_abs, erfc2)
 
            if (iatom == jatom) then
              coeff = m_one
            else
              coeff = m_two
            end if
 
            do idir = 1, 2
              do jdir = 1, 2
                stress_ewald(idir, jdir) = stress_ewald(idir, jdir) &
                  - factor*gvec(idir)*gvec(jdir) * cos_gx * (factor1 + factor2) * coeff&
                  * atom(iatom)%species%get_zval()*atom(jatom)%species%get_zval()
              end do
            end do
 
            if (abs(erfc1) > m_epsilon) then
              factor1 = exp(-gg_abs*z_ij)*erfc1
            else
              factor1 = m_zero
            end if
            if (abs(erfc2) > m_epsilon) then
              factor2 = exp(gg_abs*z_ij)*erfc2
            else
              factor2 = m_zero
            end if
 
            e_ewald = m_half * m_pi/latt%rcell_volume * coeff &
              * atom(iatom)%species%get_zval() * atom(jatom)%species%get_zval() &
              * cos_gx / gg_abs * (factor1 + factor2)
 
            do idir = 1, 2
              stress_ewald(idir, idir) = stress_ewald(idir, idir) + e_ewald
            end do
 
          end do !jatom
        end do !iatom
      end do !iy
    end do !ix
 
    !call comm_allreduce(this%dist%mpi_grp, stress_Ewald)
 
    stress_ewald = stress_ewald / latt%rcell_volume
 
    pop_sub(ewald_2d_stress)
  end subroutine ewald_2d_stress
 
  ! ---------------------------------------------------------
  real(real64) function screening_function_2d(alpha, z_ij, gg_abs, erfc_arg) result(factor)
    real(real64),  intent(in)  :: alpha
    real(real64),  intent(in)  :: z_ij
    real(real64),  intent(in)  :: gg_abs
    real(real64),  intent(out) :: erfc_arg
 
    real(real64) :: arg
 
    arg = -alpha*z_ij + m_half*gg_abs/alpha
    erfc_arg = erfc(arg)
    factor = m_two*alpha*(m_one/gg_abs + z_ij)*erfc_arg - m_two/sqrt(m_pi)*exp(-arg**2)
    factor = factor*exp(-gg_abs*z_ij)
 
  end function screening_function_2d
 
  ! ---------------------------------------------------------
 
  subroutine ion_interaction_test(space, latt, atom, natoms, pos, lsize, &
    namespace, mc)
    class(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(atom_t), contiguous, intent(in) :: atom(:)
    integer, intent(in) :: natoms
    real(real64), intent(in) :: pos(1:space%dim,1:natoms)
    real(real64), contiguous, intent(in) :: lsize(:)
    type(namespace_t), intent(in) :: namespace
    type(multicomm_t), intent(in) :: mc
 
    type(ion_interaction_t) :: ion_interaction
    real(real64) :: energy
    real(real64), allocatable :: force(:, :), force_components(:, :, :)
    real(real64) :: energy_components(1:ION_NUM_COMPONENTS)
    integer :: iatom, idir
 
    push_sub(ion_interaction_test)
 
    call ion_interaction_init(ion_interaction, namespace, space, natoms)
    call ion_interaction_init_parallelization(ion_interaction, natoms, mc)
 
    safe_allocate(force(1:space%dim, 1:natoms))
    safe_allocate(force_components(1:space%dim, 1:natoms, 1:ion_num_components))
 
    call ion_interaction_calculate(ion_interaction, space, latt, atom, natoms, pos, lsize, energy, force, &
      energy_components = energy_components, force_components = force_components)
 
    call messages_write('Ionic energy        =')
    call messages_write(energy, fmt = '(f20.10)')
    call messages_info(namespace=namespace)
 
    call messages_write('Real space energy   =')
    call messages_write(energy_components(ion_component_real), fmt = '(f20.10)')
    call messages_info(namespace=namespace)
 
    call messages_write('Self energy         =')
    call messages_write(energy_components(ion_component_self), fmt = '(f20.10)')
    call messages_info(namespace=namespace)
 
    call messages_write('Fourier energy      =')
    call messages_write(energy_components(ion_component_fourier), fmt = '(f20.10)')
    call messages_info(namespace=namespace)
 
    call messages_info(namespace=namespace)
 
    do iatom = 1, natoms
      call messages_write('Ionic force         atom')
      call messages_write(iatom)
      call messages_write(' =')
      do idir = 1, space%dim
        call messages_write(force(idir, iatom), fmt = '(f20.10)')
      end do
      call messages_info(namespace=namespace)
 
      call messages_write('Real space force    atom')
      call messages_write(iatom)
      call messages_write(' =')
      do idir = 1, space%dim
        call messages_write(force_components(idir, iatom, ion_component_real), fmt = '(f20.10)')
      end do
      call messages_info(namespace=namespace)
 
      call messages_write('Fourier space force atom')
      call messages_write(iatom)
      call messages_write(' =')
      do idir = 1, space%dim
        call messages_write(force_components(idir, iatom, ion_component_fourier), fmt = '(f20.10)')
      end do
      call messages_info(namespace=namespace)
 
      call messages_info(namespace=namespace)
    end do
 
    safe_deallocate_a(force)
    safe_deallocate_a(force_components)
 
    call ion_interaction_end(ion_interaction)
 
    pop_sub(ion_interaction_test)
  end subroutine ion_interaction_test
 
end module ion_interaction_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: