poisson_fft.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!!
!! 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 poisson_fft_oct_m
  use accel_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use debug_oct_m
  use fft_oct_m
  use fourier_space_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use math_oct_m
  use mesh_cube_parallel_map_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_cutoff_oct_m
  use profiling_oct_m
  use space_oct_m
  use splines_oct_m
  use submesh_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use xc_cam_oct_m, only: cam_null
 
  implicit none
  private
  public ::                  &
    poisson_fft_t,           &
    poisson_fft_init,        &
    poisson_fft_get_kernel,  &
    poisson_fft_end,         &
    dpoisson_fft_solve,      &
    zpoisson_fft_solve,      &
    spoisson_fft_solve,      &
    cpoisson_fft_solve,      &
    dpoisson_fft_solve_batch,&
    zpoisson_fft_solve_batch,&
    spoisson_fft_solve_batch,&
    cpoisson_fft_solve_batch
 
  integer, public, parameter ::             &
    POISSON_FFT_KERNEL_NONE      = -1,      &
    poisson_fft_kernel_sph       =  0,      &
    poisson_fft_kernel_cyl       =  1,      &
    poisson_fft_kernel_pla       =  2,      &
    poisson_fft_kernel_nocut     =  3,      &
    poisson_fft_kernel_corrected =  4,      &
    poisson_fft_kernel_hockney   =  5
 
  type poisson_fft_t
    ! Components are public by default
    type(fourier_space_op_t) :: coulb
    integer                  :: kernel
    real(real64)             :: soft_coulb_param
  end type poisson_fft_t
 
  real(real64), parameter :: TOL_VANISHING_Q = 1e-6_real64
contains
 
  subroutine poisson_fft_init(this, namespace, space, cube, kernel, soft_coulb_param, fullcube)
    type(poisson_fft_t), intent(out)   :: this
    type(namespace_t),   intent(in)    :: namespace
    class(space_t),      intent(in)    :: space
    type(cube_t),        intent(inout) :: cube
    integer,             intent(in)    :: kernel
    real(real64), optional, intent(in) :: soft_coulb_param
    type(cube_t), optional, intent(in) :: fullcube
 
    real(real64) :: qvector(1:space%dim)
 
    push_sub(poisson_fft_init)
 
    this%kernel = kernel
    this%soft_coulb_param = optional_default(soft_coulb_param, m_zero)
    qvector = m_zero
    call this%coulb%init(space, qvector, cam_null)
 
    call poisson_fft_get_kernel(namespace, space, cube, this%coulb, kernel, soft_coulb_param, fullcube)
 
    pop_sub(poisson_fft_init)
  end subroutine poisson_fft_init
 
  subroutine poisson_fft_get_kernel(namespace, space, cube, coulb, kernel, soft_coulb_param, fullcube)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
    integer,                  intent(in)    :: kernel
    real(real64), optional,   intent(in)    :: soft_coulb_param
    type(cube_t), optional,   intent(in)    :: fullcube
 
    push_sub(poisson_fft_get_kernel)
 
    if (coulb%mu > m_epsilon) then
      if (space%dim /= 3 .or. kernel /= poisson_fft_kernel_nocut) then
        message(1) = "The screened Coulomb potential is only implemented in 3D for PoissonFFTKernel=fft_nocut."
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
 
    if (kernel == poisson_fft_kernel_hockney) then
      if (.not. present(fullcube)) then
        message(1) = "Hockney's FFT-kernel needs cube of full unit cell "
        call messages_fatal(1, namespace=namespace)
      else
        if (.not. allocated(fullcube%fft)) then
          message(1) = "Hockney's FFT-kernel needs PoissonSolver=fft"
          call messages_fatal(1, namespace=namespace)
        end if
      end if
    end if
 
 
    select case (space%dim)
    case (1)
      assert(present(soft_coulb_param))
      select case (kernel)
      case (poisson_fft_kernel_sph)
        call poisson_fft_build_1d_0d(namespace, cube, coulb, soft_coulb_param)
      case (poisson_fft_kernel_nocut)
        call poisson_fft_build_1d_1d(cube, coulb, soft_coulb_param)
      case default
        message(1) = "Invalid Poisson FFT kernel for 1D."
        call messages_fatal(1, namespace=namespace)
      end select
 
    case (2)
      select case (kernel)
      case (poisson_fft_kernel_sph)
        call poisson_fft_build_2d_0d(namespace, cube, coulb)
      case (poisson_fft_kernel_cyl)
        call poisson_fft_build_2d_1d(namespace, cube, coulb)
      case (poisson_fft_kernel_nocut)
        call poisson_fft_build_2d_2d(cube, coulb)
      case default
        message(1) = "Invalid Poisson FFT kernel for 2D."
        call messages_fatal(1, namespace=namespace)
      end select
 
    case (3)
      select case (kernel)
      case (poisson_fft_kernel_sph, poisson_fft_kernel_corrected)
        call poisson_fft_build_3d_0d(namespace,  cube, kernel, coulb, space%is_periodic())
 
      case (poisson_fft_kernel_cyl)
        call poisson_fft_build_3d_1d(namespace, space, cube, coulb)
 
      case (poisson_fft_kernel_pla)
        call poisson_fft_build_3d_2d(namespace, cube, coulb)
 
      case (poisson_fft_kernel_nocut)
        call poisson_fft_build_3d_3d(cube, coulb)
 
      case (poisson_fft_kernel_hockney)
        call poisson_fft_build_3d_3d_hockney(cube, coulb, fullcube)
 
      case default
        message(1) = "Invalid Poisson FFT kernel for 3D."
        call messages_fatal(1, namespace=namespace)
      end select
    end select
 
    pop_sub(poisson_fft_get_kernel)
  end subroutine poisson_fft_get_kernel
 
  !-----------------------------------------------------------------
 
  subroutine get_cutoff(namespace, default_r_c, r_c)
    type(namespace_t),   intent(in)  :: namespace
    real(real64),        intent(in)  :: default_r_c
    real(real64),        intent(out) :: r_c
 
    push_sub(get_cutoff)
 
    call parse_variable(namespace, 'PoissonCutoffRadius', default_r_c, r_c, units_inp%length)
 
    call messages_write('Info: Poisson Cutoff Radius     =')
    call messages_write(r_c, units = units_out%length, fmt = '(f6.1)')
    call messages_info()
 
    if (r_c > default_r_c + m_epsilon) then
      call messages_write('Poisson cutoff radius is larger than cell size.', new_line = .true.)
      call messages_write('You can see electrons in neighboring cell(s).')
      call messages_warning()
    end if
 
    pop_sub(get_cutoff)
  end subroutine get_cutoff
 
 
  function poisson_fft_singularity_3d(coulb) result(singularity_term)
    type(fourier_space_op_t), intent(in) :: coulb
 
    real(real64) :: singularity_term
 
    real(real64) :: inv_four_mu2
 
    push_sub(poisson_fft_singularity_3d)
 
    if (coulb%mu > m_epsilon) then
      inv_four_mu2 = m_one/((m_two*coulb%mu)**2)
    else
      inv_four_mu2 = m_zero
    end if
 
    ! Screened short-range coulomb potential (erfc function)
    if (coulb%mu > m_epsilon .and. abs(coulb%alpha) < m_epsilon) then
      ! Analytical limit of 4pi/q^2 * (1-beta*exp(-|q|^2/4mu^2))
      singularity_term = m_four * m_pi * inv_four_mu2 * coulb%beta
    else ! We use the user-defined value of the singularity
      ! Long-range screened singularity
      if (abs(coulb%alpha) > m_epsilon) then
        singularity_term = coulb%singularity*coulb%alpha + m_four * m_pi * inv_four_mu2 * coulb%beta
      else
        singularity_term = coulb%singularity
      end if
      ! 4pi/q^2 * alpha
    end if
 
    pop_sub(poisson_fft_singularity_3d)
 
  end function poisson_fft_singularity_3d
 
  subroutine poisson_fft_build_3d_3d(cube, coulb)
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    integer      :: n1, n2, n3, lx, ly, lz, i
    real(real64) :: modg2, modgyz, beta, modg2_cutoff, inv_four_mu2, ecut
    real(real64) :: temp(3), diag_temp(3, 3), metric(3, 3), a(3, 3)
    real(real64) :: q1, q2, q3, ux, uy, uz, a11, two_a12, two_a13, a22, two_a23, a33, four_a11
    real(real64) :: singularity_term
    real(real64), allocatable :: fft_coulb_fs(:,:,:)
 
    push_sub(poisson_fft_build_3d_3d)
 
    if (coulb%mu > m_epsilon) then
      inv_four_mu2 = m_one/((m_two*coulb%mu)**2)
    else
      inv_four_mu2 = m_zero
    end if
 
    n1 = max(1, cube%fs_n(1))
    n2 = max(1, cube%fs_n(2))
    n3 = max(1, cube%fs_n(3))
 
    ! Define q+G = 0 term
    singularity_term = poisson_fft_singularity_3d(coulb)
 
    ! store the Fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:n1, 1:n2, 1:n3))
 
    ! G vector cutoff
    temp = m_two * m_pi / cube%length()
    ecut = fft_get_ecut_from_box(cube%rs_n_global, cube%fs_istart, cube%latt, temp, 3, coulb%qq)
 
    ! metric = B^T B
    metric = matmul(transpose(cube%latt%klattice_primitive), cube%latt%klattice_primitive)
 
    ! A = D (B^T B) D
    diag_temp = m_zero
    do i = 1, 3
      diag_temp(i, i) = temp(i)
    enddo
    a = matmul(diag_temp, matmul(metric, diag_temp))
    a11 = a(1, 1)
    a22 = a(2, 2)
    a33 = a(3, 3)
    two_a12 = m_two * a(1, 2)
    two_a13 = m_two * a(1, 3)
    two_a23 = m_two * a(2, 3)
    four_a11 = m_four * a11
 
    q1 = coulb%qq(1)
    q2 = coulb%qq(2)
    q3 = coulb%qq(3)
 
    modg2_cutoff = m_two * ecut * 1.001_real64
 
    do lz = 1, n3
      ! u=G+q (integer-mode coords)
      uz = real(cube%fs_ifz(lz), real64) + q3
 
      do ly = 1, n2
        uy = real(cube%fs_ify(ly), real64) + q2
        modgyz = a22*uy*uy + two_a23*uy*uz + a33*uz*uz
        beta = two_a12*uy + two_a13*uz
        if (modgyz - beta * beta / four_a11 > modg2_cutoff) then
          fft_coulb_fs(1:n1, ly, lz) = m_zero
          cycle
        end if
 
        do lx = 1, n1
          ux = real(cube%fs_ifx(lx), real64) + q1
          ! Cartesian |G + q|^2
          modg2 = (a11*ux + beta)*ux + modgyz
 
          if (modg2 > modg2_cutoff) then
            fft_coulb_fs(lx, ly, lz) = m_zero
          else if (modg2 > tol_vanishing_q) then
            !Screened coulomb potential (erfc function)
            if (coulb%mu > m_epsilon) then
              if(abs(coulb%alpha) > m_epsilon) then ! CAM
                fft_coulb_fs(lx, ly, lz) = m_four * m_pi / modg2 * (coulb%alpha - coulb%beta * expm1(-modg2*inv_four_mu2))
              else ! purely short-range screened
                fft_coulb_fs(lx, ly, lz) = m_four * m_pi / modg2 * coulb%beta * (-expm1(-modg2 * inv_four_mu2))
              end if
            else
              if (abs(coulb%alpha) > m_epsilon) then ! global screened hybrids
                fft_coulb_fs(lx, ly, lz) = m_four * m_pi / modg2 * coulb%alpha
              else ! Bare interaction
                fft_coulb_fs(lx, ly, lz) = m_four * m_pi / modg2
              end if
            end if
          else ! This is the term q+G = 0
            fft_coulb_fs(lx, ly, lz) = singularity_term
          end if
 
        end do
      end do
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
 
    pop_sub(poisson_fft_build_3d_3d)
 
  end subroutine poisson_fft_build_3d_3d
 
 
  subroutine poisson_fft_build_3d_3d_hockney(cube, coulb, fullcube)
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
    type(cube_t),             intent(in)    :: fullcube
 
    integer :: ix, iy, iz, ixx(3), db(3), nfs(3), nrs(3), nfs_s(3), nrs_s(3), dnrs(3)
    real(real64) :: temp(3), modg2, weight
    real(real64) :: gg(3)
    real(real64), allocatable :: fft_Coulb_small_RS(:,:,:,:)
    real(real64), allocatable :: fft_Coulb_RS(:,:,:,:)
    complex(real64), allocatable :: fft_Coulb_small_FS(:,:,:,:)
    complex(real64), allocatable :: fft_Coulb_FS(:,:,:,:)
    integer, parameter :: howmany = 1
 
    push_sub(poisson_fft_build_3d_3d_hockney)
 
    assert(abs(coulb%mu) < m_epsilon)
    assert(cube%batch_capacity == 1)
    assert(fullcube%batch_capacity == 1)
 
    ! dimensions of large boxes
    nfs(1:3) = fullcube%fs_n_global(1:3)
    nrs(1:3) = fullcube%rs_n_global(1:3)
 
    safe_allocate(fft_coulb_fs(1:nfs(1),1:nfs(2),1:nfs(3), 1:howmany))
    safe_allocate(fft_coulb_rs(1:nrs(1),1:nrs(2),1:nrs(3), 1:howmany))
 
    ! dimensions of small boxes x_s
    nfs_s(1:3) = cube%fs_n_global(1:3)
    nrs_s(1:3) = cube%rs_n_global(1:3)
 
    safe_allocate(fft_coulb_small_fs(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3), 1:howmany))
    safe_allocate(fft_coulb_small_rs(1:nrs_s(1),1:nrs_s(2),1:nrs_s(3), 1:howmany))
 
    ! build full periodic Coulomb potenital in Fourier space
    fft_coulb_fs = m_zero
 
    db(1:3) = fullcube%rs_n_global(1:3)
    temp(1:3) = m_two*m_pi/(db(1:3)*cube%spacing(1:3))
 
    do iz = 1, nfs(3)
      ixx(3) = pad_feq(iz, db(3), .true.)
      do iy = 1, nfs(2)
        ixx(2) = pad_feq(iy, db(2), .true.)
        do ix = 1, nfs(1)
          ixx(1) = pad_feq(ix, db(1), .true.)
 
          call fft_gg_transform(ixx, temp, 3, cube%latt, coulb%qq, gg, modg2)
 
          if (abs(modg2) > tol_vanishing_q) then
            fft_coulb_fs(ix, iy, iz, 1) = m_one/modg2
          else
            fft_coulb_fs(ix, iy, iz, 1) = m_zero
          end if
        end do
      end do
    end do
 
    ! Full range hybrids weight
    weight = m_four*m_pi
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
    do iz = 1, nfs(3)
      do iy = 1, nfs(2)
        do ix = 1, nfs(1)
          fft_coulb_fs(ix, iy, iz, 1) = weight*fft_coulb_fs(ix, iy, iz, 1)
        end do
      end do
    end do
 
    ! get periodic Coulomb potential in real space
    call dfft_backward(fullcube%fft, fft_coulb_fs, fft_coulb_rs)
 
    !  copy to small box by respecting this pattern
    !  full periodic coulomb: |abc--------------------------xyz|
    !                Hockney: |abcxyz|
    dnrs = nrs - nrs_s
 
    do iz = 1, nrs_s(3)
      ixx(3) = iz
      if (iz > nrs_s(3)/2+1) ixx(3) = ixx(3) + dnrs(3)
      do iy = 1, nrs_s(2)
        ixx(2) = iy
        if (iy > nrs_s(2)/2+1) ixx(2) = ixx(2) + dnrs(2)
        do ix = 1, nrs_s(1)
          ixx(1) = ix
          if (ix > nrs_s(1)/2+1) ixx(1) = ixx(1) + dnrs(1)
          fft_coulb_small_rs(ix, iy, iz, 1) = fft_coulb_rs(ixx(1),ixx(2),ixx(3), 1)
        end do
      end do
    end do
    ! make Hockney kernel in Fourier space
    call dfft_forward(cube%fft, fft_coulb_small_rs, fft_coulb_small_fs)
    !dummy copy for type conversion
    fft_coulb_small_rs(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3),1:howmany) = &
      real( fft_Coulb_small_FS(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3),1:howmany), real64)
 
 
    ! Restrict array to local part to support pfft
    ! For FFTW this reduces simply to the full array
    call coulb%dset_op(cube, &
      fft_coulb_small_rs(cube%fs_istart(1):cube%fs_istart(1)+cube%fs_n(1), &
      cube%fs_istart(2):cube%fs_istart(2)+cube%fs_n(2), &
      cube%fs_istart(3):cube%fs_istart(3)+cube%fs_n(3), 1))
 
    safe_deallocate_a(fft_coulb_fs)
    safe_deallocate_a(fft_coulb_rs)
    safe_deallocate_a(fft_coulb_small_fs)
    safe_deallocate_a(fft_coulb_small_rs)
 
    pop_sub(poisson_fft_build_3d_3d_hockney)
 
  end subroutine poisson_fft_build_3d_3d_hockney
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_3d_2d(namespace, cube, coulb)
    type(namespace_t),        intent(in)    :: namespace
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    integer :: ix, iy, iz, ixx(3), db(3)
    integer :: lx, ly, lz, n1, n2, n3
    real(real64) :: temp(3), modg2, ecut, weight
    real(real64) :: gpar, gz, r_c, gg(3), default_r_c
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
 
    push_sub(poisson_fft_build_3d_2d)
 
    db(1:3) = cube%rs_n_global(1:3)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_four*m_pi
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
 
    !%Variable PoissonCutoffRadius
    !%Type float
    !%Section Hamiltonian::Poisson
    !%Description
    !% When <tt>PoissonSolver = fft</tt> and <tt>PoissonFFTKernel</tt> is neither <tt>multipole_corrections</tt>
    !% nor <tt>fft_nocut</tt>,
    !% this variable controls the distance after which the electron-electron interaction goes to zero.
    !% A warning will be written if the value is too large and will cause spurious interactions between images.
    !% The default is half of the FFT box max dimension in a finite direction.
    !%End
 
    default_r_c = db(3)*cube%spacing(3)/m_two
    call get_cutoff(namespace, default_r_c, r_c)
 
    n1 = max(1, cube%fs_n(1))
    n2 = max(1, cube%fs_n(2))
    n3 = max(1, cube%fs_n(3))
    ! store the Fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:n1, 1:n2, 1:n3))
    fft_coulb_fs = m_zero
 
    temp(1:3) = m_two*m_pi/(db(1:3)*cube%spacing(1:3))
 
    ecut = fft_get_ecut_from_box(cube%rs_n_global, cube%fs_istart, cube%latt, temp, 2, coulb%qq)
 
    do lz = 1, n3
      iz = cube%fs_istart(3) + lz - 1
      ixx(3) = pad_feq(iz, db(3), .true.)
      do ly = 1, n2
        iy = cube%fs_istart(2) + ly - 1
        ixx(2) = pad_feq(iy, db(2), .true.)
        do lx = 1, n1
          ix = cube%fs_istart(1) + lx - 1
          ixx(1) = pad_feq(ix, db(1), .true.)
 
          call fft_gg_transform(ixx, temp, 2, cube%latt, coulb%qq, gg, modg2)
 
          if(sum(gg(1:2)**2) > m_two*ecut*1.001_real64) cycle
 
          if (abs(modg2) > tol_vanishing_q) then
            gz = abs(gg(3))
            gpar = hypot(gg(1), gg(2))
            ! note: if gpar = 0, then modg2 = gz**2
            fft_coulb_fs(lx, ly, lz) = poisson_cutoff_3d_2d(gpar,gz,r_c)/modg2
          else
            fft_coulb_fs(lx, ly, lz) = -m_half*r_c**2
          end if
          fft_coulb_fs(lx, ly, lz) = weight*fft_coulb_fs(lx, ly, lz)
        end do
      end do
 
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
    pop_sub(poisson_fft_build_3d_2d)
  end subroutine poisson_fft_build_3d_2d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_3d_1d(namespace, space, cube, coulb)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    type(spline_t)     :: cylinder_cutoff_f
    real(real64), allocatable :: x(:), y(:)
    integer :: ix, iy, iz, ixx(3), db(3), k, ngp
    integer :: lx, ly, lz, n1, n2, n3, lxx(3)
    real(real64) :: temp(3), modg2, xmax, weight
    real(real64) :: gperp, gx, gy, gz, r_c, gg(3), default_r_c
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
 
    push_sub(poisson_fft_build_3d_1d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_four*m_pi
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
 
    db(1:3) = cube%rs_n_global(1:3)
 
    default_r_c = maxval(db(2:3)*cube%spacing(2:3)/m_two)
    call get_cutoff(namespace, default_r_c, r_c)
 
    n1 = max(1, cube%fs_n(1))
    n2 = max(1, cube%fs_n(2))
    n3 = max(1, cube%fs_n(3))
    ! store the Fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:n1, 1:n2, 1:n3))
    fft_coulb_fs = m_zero
 
    temp(1:3) = m_two*m_pi/(db(1:3)*cube%spacing(1:3))
 
    if (.not. space%is_periodic()) then
      ngp = 8*db(2)
      safe_allocate(x(1:ngp))
      safe_allocate(y(1:ngp))
    end if
 
    ! Note(Alex) This loop ordering results in bad memory access for fft_Coulb_FS
    ! It should be refactored
    do lx = 1, n1
      ix = cube%fs_istart(1) + lx - 1
      ixx(1) = pad_feq(ix, db(1), .true.)
      lxx(1) = ixx(1) - cube%fs_istart(1) + 1
      gx = temp(1)*ixx(1)
 
      if (.not. space%is_periodic()) then
        call spline_init(cylinder_cutoff_f)
        xmax = norm2(temp(2:3)*db(2:3))/2
        do k = 1, ngp
          x(k) = (k-1)*(xmax/(ngp-1))
          y(k) = poisson_cutoff_3d_1d_finite(gx, x(k), norm2(cube%latt%rlattice_primitive(:, 1)), &
            maxval(norm2(cube%latt%rlattice_primitive(:, 2:3), dim=1)))
        end do
        call spline_fit(ngp, x, y, cylinder_cutoff_f, m_zero)
      end if
 
      do ly = 1, n2
        iy = cube%fs_istart(2) + ly - 1
        ixx(2) = pad_feq(iy, db(2), .true.)
        lxx(2) = ixx(2) - cube%fs_istart(2) + 1
        do lz = 1, n3
          iz = cube%fs_istart(3) + lz - 1
          ixx(3) = pad_feq(iz, db(3), .true.)
          lxx(3) = ixx(3) - cube%fs_istart(3) + 1
 
          call fft_gg_transform(ixx, temp, 1, cube%latt, coulb%qq, gg, modg2)
 
          if (abs(modg2) > tol_vanishing_q) then
            gperp = hypot(gg(2), gg(3))
            if (space%periodic_dim == 1) then
              if (gperp > r_c) then
                fft_coulb_fs(lx, ly, lz) = m_zero
              else
                fft_coulb_fs(lx, ly, lz) = poisson_cutoff_3d_1d(abs(gx), gperp, r_c)/modg2
              end if
            else if (.not. space%is_periodic()) then
              gy = gg(2)
              gz = gg(3)
              if ((gz >= m_zero) .and. (gy >= m_zero)) then
                fft_coulb_fs(lx, ly, lz) = spline_eval(cylinder_cutoff_f, gperp)
              end if
              if ((gz >= m_zero) .and. (gy < m_zero)) then
                fft_coulb_fs(lx, ly, lz) = fft_coulb_fs(lx, -lxx(2) + 1, lz)
              end if
              if ((gz < m_zero) .and. (gy >= m_zero)) then
                fft_coulb_fs(lx, ly, lz) = fft_coulb_fs(lx, ly, -lxx(3) + 1)
              end if
              if ((gz < m_zero) .and. (gy < m_zero)) then
                fft_coulb_fs(lx, ly, lz) = fft_coulb_fs(lx, -lxx(2) + 1, -lxx(3) + 1)
              end if
            end if
 
          else
            if (space%periodic_dim == 1) then
              fft_coulb_fs(lx, ly, lz) = -(m_half*log(r_c) - m_fourth)*r_c**2
            else if (.not. space%is_periodic()) then
              fft_coulb_fs(lx, ly, lz) = poisson_cutoff_3d_1d_finite(m_zero, m_zero, &
                norm2(cube%latt%rlattice_primitive(:, 1)), maxval(norm2(cube%latt%rlattice_primitive(:, 2:3), dim=1)))
            end if
 
          end if
          fft_coulb_fs(lx, ly, lz) = weight*fft_coulb_fs(lx, ly, lz)
        end do
      end do
 
      if (.not. space%is_periodic()) then
        call spline_end(cylinder_cutoff_f)
      end if
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
    safe_deallocate_a(x)
    safe_deallocate_a(y)
    pop_sub(poisson_fft_build_3d_1d)
  end subroutine poisson_fft_build_3d_1d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_3d_0d(namespace, cube, kernel, coulb, is_periodic)
    type(namespace_t),        intent(in)    :: namespace
    type(cube_t),             intent(in)    :: cube
    integer,                  intent(in)    :: kernel
    type(fourier_space_op_t), intent(inout) :: coulb
    logical,                  intent(in)    :: is_periodic
 
    integer :: ix, iy, iz, ixx(3), db(3), lx, ly, lz, n1, n2, n3
    real(real64) :: temp(3), modg2, ecut, weight
    real(real64) :: r_c, gg(3), default_r_c
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
    real(real64) :: axis(3,3)
 
    push_sub(poisson_fft_build_3d_0d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_four*m_pi
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
 
    db(1:3) = cube%rs_n_global(1:3)
 
    if (kernel /= poisson_fft_kernel_corrected) then
 
      ! This is the real-space cutoff
      do ix = 1, 3
        axis(:,ix) = cube%latt%rlattice_primitive(:, ix) * cube%spacing(ix) * db(ix) / m_two
      end do
 
      default_r_c = m_huge
      do ix = 1, 3
        iy = mod(ix, 3)+1
        iz = mod(ix+1, 3)+1
 
        ! For orthogonal cells, this is determined by the size of the cube
        if (.not. cube%latt%nonorthogonal) then
          temp(1:3) = axis(:, ix)
          default_r_c = min(default_r_c, norm2(temp(1:3)))
        else
          ! At the moment, this codepath is only called for DFT+U submesh Poisson solver
          !
          ! For non-orthogonal cells, the relevant length is the distance between two planes of
          ! the parallelepiped. This ensures that we draw a sphere that touches the borders of the box,
          ! thus avoiding contribution from periodic replicas
          ! This distance is given by the usual formula of the distance from a point to a plan
          temp = dcross_product(axis(:, iy), axis(:, iz))
          temp = temp / norm2(temp)
          default_r_c = min(default_r_c, dot_product(temp, axis(:, ix)-axis(:, iy)))
        end if
      end do
      call get_cutoff(namespace, default_r_c, r_c)
    end if
 
    n1 = max(1, cube%fs_n(1))
    n2 = max(1, cube%fs_n(2))
    n3 = max(1, cube%fs_n(3))
 
    ! store the fourier transform of the Coulomb interaction
    ! store only the relevant part if PFFT is used
    safe_allocate(fft_coulb_fs(1:n1,1:n2,1:n3))
    fft_coulb_fs = m_zero
 
    temp(1:3) = m_two*m_pi/(db(1:3)*cube%spacing(1:3))
 
    ecut = fft_get_ecut_from_box(cube%rs_n_global, cube%fs_istart, cube%latt, temp, 3, coulb%qq)
 
    do lz = 1, n3
      iz = cube%fs_istart(3) + lz - 1
      ixx(3) = pad_feq(iz, db(3), .true.)
      do ly = 1, n2
        iy = cube%fs_istart(2) + ly - 1
        ixx(2) = pad_feq(iy, db(2), .true.)
        do lx = 1, n1
          ix = cube%fs_istart(1) + lx - 1
          ixx(1) = pad_feq(ix, db(1), .true.)
 
          call fft_gg_transform(ixx, temp, 0, cube%latt, coulb%qq, gg, modg2)
 
          ! At the moment this is only done for periodic space, so for DFT+U Coulomb integrals
          if(modg2 > m_two*ecut*1.001_real64 .and. is_periodic) cycle
 
          if (abs(modg2) > tol_vanishing_q) then
            select case (kernel)
            case (poisson_fft_kernel_sph)
              fft_coulb_fs(lx, ly, lz) = weight*poisson_cutoff_3d_0d(sqrt(modg2),r_c)/modg2
            case (poisson_fft_kernel_corrected)
              fft_coulb_fs(lx, ly, lz) = weight/modg2
            end select
          else
            select case (kernel)
            case (poisson_fft_kernel_sph)
              fft_coulb_fs(lx, ly, lz) = weight*r_c**2/m_two
            case (poisson_fft_kernel_corrected)
              fft_coulb_fs(lx, ly, lz) = m_zero
            end select
          end if
        end do
      end do
    end do
 
    call coulb%dset_op(cube, fft_coulb_fs, in_device = (kernel /= poisson_fft_kernel_corrected))
 
    safe_deallocate_a(fft_coulb_fs)
    pop_sub(poisson_fft_build_3d_0d)
  end subroutine poisson_fft_build_3d_0d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_2d_0d(namespace, cube, coulb)
    type(namespace_t),        intent(in)    :: namespace
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    type(spline_t) :: besselintf
    integer :: i, ix, iy, ixx(2), db(2), npoints
    real(real64) :: temp(2), vec, r_c, maxf, dk, default_r_c, weight
    real(real64), allocatable :: x(:), y(:)
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
 
    push_sub(poisson_fft_build_2d_0d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_one
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
    db(1:2) = cube%rs_n_global(1:2)
 
    default_r_c = maxval(db(1:2)*cube%spacing(1:2)/m_two)
    call get_cutoff(namespace, default_r_c, r_c)
 
    call spline_init(besselintf)
 
    ! store the fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3)))
    fft_coulb_fs = m_zero
    temp(1:2) = m_two*m_pi/(db(1:2)*cube%spacing(1:2))
 
    maxf = r_c * norm2(temp(1:2)*db(1:2))/2
    dk = 0.25_real64 ! This seems to be reasonable.
    npoints = nint(maxf/dk)
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
    x(1) = m_zero
    y(1) = m_zero
    do i = 2, npoints
      x(i) = (i-1) * maxf / (npoints-1)
      y(i) = y(i-1) + poisson_cutoff_2d_0d(x(i-1), x(i))
    end do
    call spline_fit(npoints, x, y, besselintf, m_zero)
 
    do iy = 1, cube%fs_n_global(2)
      ixx(2) = pad_feq(iy, db(2), .true.)
      do ix = 1, cube%fs_n_global(1)
        ixx(1) = pad_feq(ix, db(1), .true.)
        vec = norm2(temp(1:2)*ixx(1:2))
        ! extra check to avoid extrapolation which leads to an error in gsl
        if (vec*r_c >= x(npoints)) then
          fft_coulb_fs(ix, iy, 1) = weight * y(npoints)
        else if (vec > m_zero) then
          fft_coulb_fs(ix, iy, 1) = weight * (m_two * m_pi / vec) * spline_eval(besselintf, vec*r_c)
        else
          fft_coulb_fs(ix, iy, 1) = weight * m_two * m_pi * r_c
        end if
      end do
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
    safe_deallocate_a(x)
    safe_deallocate_a(y)
    call spline_end(besselintf)
    pop_sub(poisson_fft_build_2d_0d)
  end subroutine poisson_fft_build_2d_0d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_2d_1d(namespace, cube, coulb)
    type(namespace_t),        intent(in)    :: namespace
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    integer :: ix, iy, ixx(2), db(2)
    real(real64) :: temp(2), r_c, gx, gy, default_r_c, weight
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
 
    push_sub(poisson_fft_build_2d_1d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_one
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
 
    db(1:2) = cube%rs_n_global(1:2)
 
    default_r_c = db(2)*cube%spacing(2)/m_two
    call get_cutoff(namespace, default_r_c, r_c)
 
    ! store the fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3)))
    fft_coulb_fs = m_zero
    temp(1:2) = m_two*m_pi/(db(1:2)*cube%spacing(1:2))
 
    ! First, the term ix = 0 => gx = 0.
    fft_coulb_fs(1, 1, 1) = -m_four * r_c * (log(r_c)-m_one)
    do iy = 2, cube%fs_n_global(2)
      ixx(2) = pad_feq(iy, db(2), .true.)
      gy = temp(2)*ixx(2)
      fft_coulb_fs(1, iy, 1) = -m_four * poisson_cutoff_intcoslog(r_c, gy, m_one) * weight
    end do
 
    do ix = 2, cube%fs_n_global(1)
      ixx(1) = pad_feq(ix, db(1), .true.)
      gx = temp(1)*ixx(1)
      do iy = 1, cube%fs_n_global(2)
        ixx(2) = pad_feq(iy, db(2), .true.)
        gy = temp(2)*ixx(2)
        fft_coulb_fs(ix, iy, 1) = poisson_cutoff_2d_1d(gy, gx, r_c) * weight
      end do
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
 
    pop_sub(poisson_fft_build_2d_1d)
  end subroutine poisson_fft_build_2d_1d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_2d_2d(cube, coulb)
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
 
    integer :: ix, iy, ixx(2), db(2)
    real(real64) :: temp(2), vec, weight
    real(real64), allocatable :: fft_coulb_FS(:,:,:)
 
    push_sub(poisson_fft_build_2d_2d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_one
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
    db(1:2) = cube%rs_n_global(1:2)
 
    ! store the fourier transform of the Coulomb interaction
    safe_allocate(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3)))
    fft_coulb_fs = m_zero
    temp(1:2) = m_two*m_pi/(db(1:2)*cube%spacing(1:2))
 
    do iy = 1, cube%fs_n_global(2)
      ixx(2) = pad_feq(iy, db(2), .true.)
      do ix = 1, cube%fs_n_global(1)
        ixx(1) = pad_feq(ix, db(1), .true.)
        vec = sqrt((temp(1) * ixx(1))**2 + (temp(2) * ixx(2))**2)
        if (vec > m_zero) fft_coulb_fs(ix, iy, 1) = m_two * m_pi / vec * weight
      end do
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
 
    safe_deallocate_a(fft_coulb_fs)
    pop_sub(poisson_fft_build_2d_2d)
  end subroutine poisson_fft_build_2d_2d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_1d_1d(cube, coulb, poisson_soft_coulomb_param)
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
    real(real64),             intent(in)    :: poisson_soft_coulomb_param
 
    integer            :: ix, ixx
    real(real64)       :: g, weight
    real(real64), allocatable :: fft_coulb_fs(:, :, :)
 
    push_sub(poisson_fft_build_1d_1d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_one
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
    safe_allocate(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3)))
    fft_coulb_fs = m_zero
 
    ! Fourier transform of Soft Coulomb interaction.
    do ix = 1, cube%fs_n_global(1)
      ixx = pad_feq(ix, cube%rs_n_global(1), .true.)
      g = (ixx + coulb%qq(1))*m_two*m_pi/abs(cube%latt%rlattice(1,1))
      if (abs(g) > tol_vanishing_q) then
        fft_coulb_fs(ix, 1, 1) = m_two * loct_bessel_k0(poisson_soft_coulomb_param*abs(g)) * weight
      else
        fft_coulb_fs(ix, 1, 1) = coulb%singularity * m_two * weight
      end if
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
    safe_deallocate_a(fft_coulb_fs)
 
    pop_sub(poisson_fft_build_1d_1d)
  end subroutine poisson_fft_build_1d_1d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_build_1d_0d(namespace, cube, coulb, poisson_soft_coulomb_param)
    type(namespace_t),        intent(in)    :: namespace
    type(cube_t),             intent(in)    :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
    real(real64),             intent(in)    :: poisson_soft_coulomb_param
 
    integer            :: box(1), ixx(1), ix
    real(real64)       :: temp(1), g, r_c, default_r_c, weight
    real(real64), allocatable :: fft_coulb_fs(:, :, :)
 
    push_sub(poisson_fft_build_1d_0d)
 
    assert(abs(coulb%mu) < m_epsilon)
    ! Full range hybrids weight
    weight = m_one
    if(coulb%alpha > m_epsilon) weight = weight * coulb%alpha
 
    box(1:1) = cube%rs_n_global(1:1)
 
    default_r_c = box(1)*cube%spacing(1)/m_two
    call get_cutoff(namespace, default_r_c, r_c)
 
    safe_allocate(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3)))
    fft_coulb_fs = m_zero
    temp(1:1) = m_two*m_pi/(box(1:1)*cube%spacing(1:1))
 
    ! Fourier transform of Soft Coulomb interaction.
    do ix = 1, cube%fs_n_global(1)
      ixx(1) = pad_feq(ix, box(1), .true.)
      g = temp(1)*ixx(1)
      fft_coulb_fs(ix, 1, 1) = poisson_cutoff_1d_0d(g, poisson_soft_coulomb_param, r_c) * weight
    end do
 
    call coulb%dset_op(cube, op_move = fft_coulb_fs)
    safe_deallocate_a(fft_coulb_fs)
 
    pop_sub(poisson_fft_build_1d_0d)
  end subroutine poisson_fft_build_1d_0d
  !-----------------------------------------------------------------
 
 
  !-----------------------------------------------------------------
  subroutine poisson_fft_end(this)
    type(poisson_fft_t), intent(inout) :: this
 
    push_sub(poisson_fft_end)
 
    call this%coulb%end()
 
    pop_sub(poisson_fft_end)
  end subroutine poisson_fft_end
 
#include "undef.F90"
#include "real.F90"
#include "poisson_fft_inc.F90"
#include "undef.F90"
#include "complex.F90"
#include "poisson_fft_inc.F90"
 
#include "undef.F90"
#include "real32.F90"
#include "poisson_fft_inc.F90"
#include "undef.F90"
#include "complex32.F90"
#include "poisson_fft_inc.F90"
 
 
end module poisson_fft_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: