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_nullify(this%dist, natoms)
 
    if (space%periodic_dim == 1) then
      call messages_write('For systems that are periodic in 1D, the interaction between', new_line = .true.)
      call messages_write('ions is not implemented. This affects the calculation', new_line = .true.)
      call messages_write('of total energy and forces, so both are zeroed.')
      call messages_warning(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),             intent(in)    :: atom(:)
    integer,                  intent(in)    :: natoms
    real(real64),             intent(in)    :: pos(1:space%dim,1:natoms)
    real(real64),             intent(in)    :: lsize(:)
    real(real64),             intent(out)   :: energy
    real(real64),             intent(out)   :: force(:, :)
    real(real64), optional,          intent(out)   :: energy_components(:)
    real(real64), 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),             intent(in)    :: atom(:)
    real(real64),             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),             intent(in)      :: atom(:)
    real(real64),             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),             intent(in)    :: atom(:)
    real(real64),             intent(in)    :: pos(:,:)
    real(real64),             intent(in)    :: lsize(:)
    real(real64),             intent(out)   :: energy
    real(real64),             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),              intent(in)    :: atom(:)
    integer,                   intent(in)    :: natoms
    real(real64),              intent(in)    :: pos(1:space%dim,1:natoms)
    real(real64),              intent(out)   :: energy
    real(real64),              intent(out)   :: force(:, :)
    real(real64), optional,           intent(out)   :: energy_components(:)
    real(real64), optional,           intent(out)   :: force_components(:, :, :)
 
    real(real64) :: ereal, efourier, epseudo, eself
    real(real64) :: charge
 
    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)
      ! Warning added in init routine, such that it is not displayed per SCF step
      efourier = m_zero
      ! Do not confuse the user and set to zero all the other components
      ereal = m_zero
      eself = m_zero
      force = m_zero
      epseudo = m_zero
      if (present(force_components)) then
        force_components(1:space%dim, 1:natoms, ion_component_real) = m_zero
      end if
    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),             intent(in)      :: atom(:)
    real(real64),             intent(in)      :: pos(:, :)
 
    real(real64),             intent(in)      :: alpha
    real(real64),             intent(out)     :: ereal
    real(real64),             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),             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),              intent(in)    :: atom(:)
    integer,                   intent(in)    :: natoms
    real(real64),              intent(in)    :: pos(:,:)
    real(real64),              intent(inout) :: efourier
    real(real64),              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),              intent(in)    :: atom(:)
    integer,                   intent(in)    :: natoms
    real(real64),              intent(in)    :: pos(1:space%dim,1:natoms)
    real(real64),              intent(inout) :: efourier
    real(real64),              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 = sum(ix_max*latt%klattice(:, 1)**2)
    gmax_squared = min(gmax_squared, sum(iy_max*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
 
  !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),             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),              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),              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),              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),              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 = sum(ix_max*latt%klattice(:, 1)**2)
    gmax_squared = min(gmax_squared, sum(iy_max*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) &
    !$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),             intent(in)    :: atom(:)
    integer,                  intent(in)    :: natoms
    real(real64),             intent(in)    :: pos(1:space%dim,1:natoms)
    real(real64),             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: