poisson.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_oct_m
  use accel_oct_m
  use batch_oct_m
  use box_minimum_oct_m
  use calc_mode_par_oct_m
  use comm_oct_m
  use cube_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use fft_oct_m
  use fourier_space_oct_m
  use global_oct_m
  use index_oct_m
  use io_oct_m
  use io_function_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_cube_parallel_map_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
#ifdef HAVE_OPENMP
  use omp_lib
#endif
  use par_vec_oct_m
  use parser_oct_m
  use partition_oct_m
  use photon_mode_oct_m
  use poisson_cg_oct_m
  use poisson_corrections_oct_m
  use poisson_isf_oct_m
  use poisson_fft_oct_m
  use poisson_fmm_oct_m
  use poisson_psolver_oct_m
  use poisson_multigrid_oct_m
  use poisson_no_oct_m
  use profiling_oct_m
  use space_oct_m
  use stencil_oct_m
  use submesh_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                      &
    poisson_t,                   &
    poisson_fmm_t,               &
    poisson_get_solver,          &
    poisson_init,                &
    poisson_init_sm,             &
    dpoisson_solve,              &
    zpoisson_solve,              &
    dpoisson_solve_sm,           &
    zpoisson_solve_sm,           &
    poisson_solve_batch,         &
    poisson_solver_is_iterative, &
    poisson_solver_has_free_bc,  &
    poisson_end,                 &
    poisson_test,                &
    poisson_is_multigrid,        &
    poisson_slave_work,          &
    poisson_async_init,          &
    poisson_async_end,           &
    dpoisson_solve_start,        &
    dpoisson_solve_finish,       &
    zpoisson_solve_start,        &
    zpoisson_solve_finish,       &
    poisson_build_kernel,        &
    poisson_is_async,            &
    poisson_get_full_range_weight
 
  integer, public, parameter ::         &
    POISSON_DIRECT_SUM    = -1,         &
    poisson_fmm           = -4,         &
    poisson_fft           =  0,         &
    poisson_cg            =  5,         &
    poisson_cg_corrected  =  6,         &
    poisson_multigrid     =  7,         &
    poisson_isf           =  8,         &
    poisson_psolver       = 10,         &
    poisson_no            = -99,        &
    poisson_null          = -999
 
  type poisson_t
    private
    type(derivatives_t), pointer, public :: der
    integer, public           :: method = poisson_null
    integer, public           :: kernel
    type(cube_t), public      :: cube
    type(mesh_cube_parallel_map_t), public :: mesh_cube_map
    type(poisson_mg_solver_t) :: mg
    type(poisson_fft_t), public :: fft_solver
    real(real64), public   :: poisson_soft_coulomb_param
    logical :: all_nodes_default
    type(poisson_corr_t) :: corrector
    type(poisson_isf_t)  :: isf_solver
    type(poisson_psolver_t) :: psolver_solver
    type(poisson_no_t) :: no_solver
    integer :: nslaves
    logical, public :: is_dressed = .false.
    type(photon_mode_t), public :: photons
    type(poisson_fmm_t)  :: params_fmm
#ifdef HAVE_MPI
    type(MPI_Comm)  :: intercomm
    type(mpi_grp_t) :: local_grp
    logical         :: root
#endif
  end type poisson_t
 
  integer, parameter ::             &
    CMD_FINISH = 1,                 &
    cmd_poisson_solve = 2
 
contains
 
  !-----------------------------------------------------------------
  subroutine poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
    type(poisson_t),             intent(inout) :: this
    class(space_t),              intent(in)    :: space
    type(namespace_t),           intent(in)    :: namespace
    type(derivatives_t), target, intent(in)    :: der
    type(multicomm_t),           intent(in)    :: mc
    type(stencil_t),             intent(in)    :: stencil
    real(real64),      optional, intent(in)    :: qtot
    character(len=*),  optional, intent(in)    :: label
    integer,           optional, intent(in)    :: solver
    logical,           optional, intent(in)    :: verbose
    logical,           optional, intent(in)    :: force_serial
    logical,           optional, intent(in)    :: force_cmplx
 
    logical :: need_cube, isf_data_is_parallel
    integer :: default_solver, default_kernel, box(space%dim), fft_type, fft_library
    real(real64) :: fft_alpha
    character(len=60) :: str
 
    ! Make sure we do not try to initialize an already initialized solver
    assert(this%method == poisson_null)
 
    push_sub(poisson_init)
 
    if (optional_default(verbose,.true.)) then
      str = "Hartree"
      if (present(label)) str = trim(label)
      call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
    end if
 
    this%nslaves = 0
    this%der => der
 
    !%Variable DressedOrbitals
    !%Type logical
    !%Default false
    !%Section Hamiltonian::Poisson
    !%Description
    !% Allows for the calculation of coupled elecron-photon problems
    !% by applying the dressed orbital approach. Details can be found in
    !% https://arxiv.org/abs/1812.05562
    !% At the moment, N electrons in d (<=3) spatial dimensions, coupled
    !% to one photon mode can be described. The photon mode is included by
    !% raising the orbital dimension to d+1 and changing the particle interaction
    !% kernel and the local potential, where the former is included automatically,
    !% but the latter needs to by added by hand as a user_defined_potential!
    !% Coordinate 1-d: electron; coordinate d+1: photon.
    !%End
    call parse_variable(namespace, 'DressedOrbitals', .false., this%is_dressed)
    call messages_print_var_value('DressedOrbitals', this%is_dressed, namespace=namespace)
    if (this%is_dressed) then
      assert(present(qtot))
      call messages_experimental('Dressed Orbitals', namespace=namespace)
      assert(qtot > m_zero)
      call photon_mode_init(this%photons, namespace, der%dim-1)
      call photon_mode_set_n_electrons(this%photons, qtot)
      if(.not.allocated(this%photons%pol_dipole)) then
        call photon_mode_compute_dipoles(this%photons, der%mesh)
      end if
      if (this%photons%nmodes > 1) then
        call messages_not_implemented('DressedOrbitals for more than one photon mode', namespace=namespace)
      end if
    end if
 
    this%all_nodes_default = .false.
#ifdef HAVE_MPI
    if (.not. optional_default(force_serial, .false.)) then
      !%Variable ParallelizationPoissonAllNodes
      !%Type logical
      !%Default true
      !%Section Execution::Parallelization
      !%Description
      !% When running in parallel, this variable selects whether the
      !% Poisson solver should divide the work among all nodes or only
      !% among the parallelization-in-domains groups.
      !%End
 
      call parse_variable(namespace, 'ParallelizationPoissonAllNodes', .true., this%all_nodes_default)
    end if
#endif
 
    !%Variable PoissonSolver
    !%Type integer
    !%Section Hamiltonian::Poisson
    !%Description
    !% Defines which method to use to solve the Poisson equation. Some incompatibilities apply depending on
    !% dimensionality, periodicity, etc.
    !% For a comparison of the accuracy and performance of the methods in Octopus, see P Garcia-Risue&ntilde;o,
    !% J Alberdi-Rodriguez <i>et al.</i>, <i>J. Comp. Chem.</i> <b>35</b>, 427-444 (2014)
    !% or <a href=http://arxiv.org/abs/1211.2092>arXiV</a>.
    !% Defaults:
    !% <br> 1D and 2D: <tt>fft</tt>.
    !% <br> 3D: <tt>cg_corrected</tt> if curvilinear, <tt>isf</tt> if not periodic, <tt>fft</tt> if periodic.
    !% <br> Dressed orbitals: <tt>direct_sum</tt>.
    !%Option NoPoisson -99
    !% Do not use a Poisson solver at all.
    !%Option FMM -4
    !% (Experimental) Fast multipole method. Requires FMM library.
    !%Option direct_sum -1
    !% Direct evaluation of the Hartree potential (only for finite systems).
    !%Option fft 0
    !% The Poisson equation is solved using FFTs. A cutoff technique
    !% for the Poisson kernel is selected so the proper boundary
    !% conditions are imposed according to the periodicity of the
    !% system. This can be overridden by the <tt>PoissonFFTKernel</tt>
    !% variable. To choose the FFT library use <tt>FFTLibrary</tt>
    !%Option cg 5
    !% Conjugate gradients (only for finite systems).
    !%Option cg_corrected 6
    !% Conjugate gradients, corrected for boundary conditions (only for finite systems).
    !%Option multigrid 7
    !% Multigrid method (only for finite systems).
    !%Option isf 8
    !% Interpolating Scaling Functions Poisson solver (only for finite systems).
    !%Option psolver 10
    !% Solver based on Interpolating Scaling Functions as implemented in the PSolver library.
    !% Parallelization in k-points requires <tt>PoissonSolverPSolverParallelData</tt> = no.
    !% Requires the PSolver external library.
    !%End
 
    default_solver = poisson_fft
 
    if (space%dim == 3 .and. .not. space%is_periodic()) default_solver = poisson_isf
 
#ifdef HAVE_CLFFT
    ! this is disabled, since the difference between solvers are big
    ! enough to cause problems with the tests.
    ! if (accel_is_enabled()) default_solver = POISSON_FFT
#endif
#ifdef HAVE_CUDA
    if(accel_is_enabled()) default_solver = poisson_fft
#endif
 
    if (space%dim > 3) default_solver = poisson_no ! Kernel for higher dimensions is not implemented.
 
    if (der%mesh%use_curvilinear) then
      select case (space%dim)
      case (1)
        default_solver = poisson_direct_sum
      case (2)
        default_solver = poisson_direct_sum
      case (3)
        default_solver = poisson_cg_corrected
      end select
    end if
 
    if (this%is_dressed) default_solver = poisson_direct_sum
 
    if (.not. present(solver)) then
      call parse_variable(namespace, 'PoissonSolver', default_solver, this%method)
    else
      this%method = solver
    end if
    if (.not. varinfo_valid_option('PoissonSolver', this%method)) call messages_input_error(namespace, 'PoissonSolver')
    if (optional_default(verbose, .true.)) then
      select case (this%method)
      case (poisson_direct_sum)
        str = "direct sum"
      case (poisson_fmm)
        str = "fast multipole method"
      case (poisson_fft)
        str = "fast Fourier transform"
      case (poisson_cg)
        str = "conjugate gradients"
      case (poisson_cg_corrected)
        str = "conjugate gradients, corrected"
      case (poisson_multigrid)
        str = "multigrid"
      case (poisson_isf)
        str = "interpolating scaling functions"
      case (poisson_psolver)
        str = "interpolating scaling functions (from BigDFT)"
      case (poisson_no)
        str = "no Poisson solver - Hartree set to 0"
      end select
      write(message(1),'(a,a,a)') "The chosen Poisson solver is '", trim(str), "'"
      call messages_info(1, namespace=namespace)
    end if
 
    if (space%dim > 3 .and. this%method /= poisson_no) then
      call messages_input_error(namespace, 'PoissonSolver', 'Currently no Poisson solver is available for Dimensions > 3')
    end if
 
    if (this%method /= poisson_fft) then
      this%kernel = poisson_fft_kernel_none
    else
 
      ! Documentation in cube.F90
      call parse_variable(namespace, 'FFTLibrary', fftlib_fftw, fft_library)
 
      !%Variable PoissonFFTKernel
      !%Type integer
      !%Section Hamiltonian::Poisson
      !%Description
      !% Defines which kernel is used to impose the correct boundary
      !% conditions when using FFTs to solve the Poisson equation. The
      !% default is selected depending on the dimensionality and
      !% periodicity of the system:
      !% <br>In 1D, <tt>spherical</tt> if finite, <tt>fft_nocut</tt> if periodic.
      !% <br>In 2D, <tt>spherical</tt> if finite, <tt>cylindrical</tt> if 1D-periodic, <tt>fft_nocut</tt> if 2D-periodic.
      !% <br>In 3D, <tt>spherical</tt> if finite, <tt>cylindrical</tt> if 1D-periodic, <tt>planar</tt> if 2D-periodic,
      !% <tt>fft_nocut</tt> if 3D-periodic.
      !% See C. A. Rozzi et al., <i>Phys. Rev. B</i> <b>73</b>, 205119 (2006) for 3D implementation and
      !% A. Castro et al., <i>Phys. Rev. B</i> <b>80</b>, 033102 (2009) for 2D implementation.
      !%Option spherical 0
      !% FFTs using spherical cutoff (in 2D or 3D).
      !%Option cylindrical 1
      !% FFTs using cylindrical cutoff (in 2D or 3D).
      !%Option planar 2
      !% FFTs using planar cutoff (in 3D).
      !%Option fft_nocut 3
      !% FFTs without using a cutoff (for fully periodic systems).
      !%Option multipole_correction 4
      !% The boundary conditions are imposed by using a multipole expansion. Only appropriate for finite systems.
      !% Further specification occurs with variables <tt>PoissonSolverBoundaries</tt> and <tt>PoissonSolverMaxMultipole</tt>.
      !%End
 
      select case (space%dim)
      case (1)
        if (.not. space%is_periodic()) then
          default_kernel = poisson_fft_kernel_sph
        else
          default_kernel = poisson_fft_kernel_nocut
        end if
      case (2)
        if (space%periodic_dim == 2) then
          default_kernel = poisson_fft_kernel_nocut
        else if (space%is_periodic()) then
          default_kernel = space%periodic_dim
        else
          default_kernel = poisson_fft_kernel_sph
        end if
      case (3)
        default_kernel = space%periodic_dim
      end select
 
      call parse_variable(namespace, 'PoissonFFTKernel', default_kernel, this%kernel)
      if (.not. varinfo_valid_option('PoissonFFTKernel', this%kernel)) call messages_input_error(namespace, 'PoissonFFTKernel')
 
      if (optional_default(verbose,.true.)) then
        call messages_print_var_option("PoissonFFTKernel", this%kernel, namespace=namespace)
      end if
 
      ! the multipole correction kernel does not work on GPUs
      if(this%kernel == poisson_fft_kernel_corrected .and. fft_default_lib == fftlib_accel) then
        fft_default_lib = fftlib_fftw
        message(1) = 'PoissonFFTKernel=multipole_correction is not supported on GPUs'
        message(2) = 'Using FFTW to compute the FFTs on the CPU'
        call messages_info(2, namespace=namespace)
      end if
 
    end if
 
    !We assume the developer knows what he is doing by providing the solver option
    if (.not. present(solver)) then
      if (space%is_periodic() .and. this%method == poisson_direct_sum) then
        message(1) = 'A periodic system may not use the direct_sum Poisson solver.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (space%is_periodic() .and. this%method == poisson_cg_corrected) then
        message(1) = 'A periodic system may not use the cg_corrected Poisson solver.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (space%is_periodic() .and. this%method == poisson_cg) then
        message(1) = 'A periodic system may not use the cg Poisson solver.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (space%is_periodic() .and. this%method == poisson_multigrid) then
        message(1) = 'A periodic system may not use the multigrid Poisson solver.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      select case (space%dim)
      case (1)
 
        select case (space%periodic_dim)
        case (0)
          if ((this%method /= poisson_fft) .and. (this%method /= poisson_direct_sum)) then
            message(1) = 'A finite 1D system may only use fft or direct_sum Poisson solvers.'
            call messages_fatal(1, namespace=namespace)
          end if
        case (1)
          if (this%method /= poisson_fft) then
            message(1) = 'A periodic 1D system may only use the fft Poisson solver.'
            call messages_fatal(1, namespace=namespace)
          end if
        end select
 
        if (der%mesh%use_curvilinear .and. this%method /= poisson_direct_sum) then
          message(1) = 'If curvilinear coordinates are used in 1D, then the only working'
          message(2) = 'Poisson solver is direct_sum.'
          call messages_fatal(2, namespace=namespace)
        end if
 
      case (2)
 
        if ((this%method /= poisson_fft) .and. (this%method /= poisson_direct_sum)) then
          message(1) = 'A 2D system may only use fft or direct_sum solvers.'
          call messages_fatal(1, namespace=namespace)
        end if
 
        if (der%mesh%use_curvilinear .and. (this%method /= poisson_direct_sum)) then
          message(1) = 'If curvilinear coordinates are used in 2D, then the only working'
          message(2) = 'Poisson solver is direct_sum.'
          call messages_fatal(2, namespace=namespace)
        end if
 
      case (3)
 
        if (space%is_periodic() .and. this%method == poisson_fmm) then
          call messages_not_implemented('FMM for periodic systems', namespace=namespace)
        end if
 
        if (space%is_periodic() .and. this%method == poisson_isf) then
          call messages_write('The ISF solver can only be used for finite systems.')
          call messages_fatal()
        end if
 
        if (space%is_periodic() .and. this%method == poisson_fft .and. &
          this%kernel /= space%periodic_dim .and. this%kernel >= 0 .and. this%kernel <= 3) then
          write(message(1), '(a,i1,a)')'The system is periodic in ', space%periodic_dim ,' dimension(s),'
          write(message(2), '(a,i1,a)')'but Poisson solver is set for ', this%kernel, ' dimensions.'
          call messages_warning(2, namespace=namespace)
        end if
 
        if (space%is_periodic() .and. this%method == poisson_fft .and. this%kernel == poisson_fft_kernel_corrected) then
          write(message(1), '(a,i1,a)')'PoissonFFTKernel = multipole_correction cannot be used for periodic systems.'
          call messages_fatal(1, namespace=namespace)
        end if
 
        if (der%mesh%use_curvilinear .and. .not. any(this%method == [poisson_cg_corrected, poisson_multigrid])) then
          message(1) = 'If curvilinear coordinates are used, then the only working'
          message(2) = 'Poisson solver are cg_corrected and multigrid.'
          call messages_fatal(2, namespace=namespace)
        end if
        if (der%mesh%use_curvilinear .and. this%method == poisson_multigrid .and. accel_is_enabled()) then
          call messages_not_implemented('Multigrid Poisson solver with curvilinear coordinates on GPUs')
        end if
 
        select type (box => der%mesh%box)
        type is (box_minimum_t)
          if (this%method == poisson_cg_corrected) then
            message(1) = 'When using the "minimum" box shape and the "cg_corrected"'
            message(2) = 'Poisson solver, we have observed "sometimes" some non-'
            message(3) = 'negligible error. You may want to check that the "fft" or "cg"'
            message(4) = 'solver are providing, in your case, the same results.'
            call messages_warning(4, namespace=namespace)
          end if
        end select
 
        if (this%method == poisson_fmm) then
          call messages_experimental('FMM Poisson solver', namespace=namespace)
          message(1) = "The use of the FMM Poisson solver"
          message(2) = "is deprecated and will be removed in the next major release."
          call messages_warning(2, namespace=namespace)
        end if
      end select
    end if
 
    if (this%method == poisson_psolver) then
#if !(defined HAVE_PSOLVER)
      message(1) = "The PSolver Poisson solver cannot be used since the code was not compiled with the PSolver library."
      call messages_fatal(1, namespace=namespace)
#endif
    end if
 
    if (optional_default(verbose,.true.)) then
      call messages_print_with_emphasis(namespace=namespace)
    end if
 
    ! Now that we know the method, we check if we need a cube and its dimentions
    need_cube = .false.
    fft_type = fft_real
    if (optional_default(force_cmplx, .false.)) fft_type = fft_complex
 
    if (this%method == poisson_isf .or. this%method == poisson_psolver) then
      fft_type = fft_none
      box(:) = der%mesh%idx%ll(:)
      need_cube = .true.
    end if
 
    if (this%method == poisson_psolver .and. multicomm_have_slaves(mc)) then
      call messages_not_implemented('Task parallelization with PSolver Poisson solver', namespace=namespace)
    end if
 
    if (multicomm_strategy_is_parallel(mc, p_strategy_kpoints)) then
      ! Documentation in poisson_psolver.F90
      call parse_variable(namespace, 'PoissonSolverPSolverParallelData', .true., isf_data_is_parallel)
      if (this%method == poisson_psolver .and. isf_data_is_parallel) then
        call messages_not_implemented("k-point parallelization with PSolver library and", namespace=namespace)
        call messages_not_implemented("PoissonSolverPSolverParallelData = yes", namespace=namespace)
      end if
      if (this%method == poisson_fft .and. fft_library == fftlib_pfft) then
        call messages_not_implemented("k-point parallelization with PFFT library for", namespace=namespace)
        call messages_not_implemented("PFFT library for Poisson solver", namespace=namespace)
      end if
    end if
 
    if (this%method == poisson_fft) then
 
      need_cube = .true.
 
      !%Variable DoubleFFTParameter
      !%Type float
      !%Default 2.0
      !%Section Mesh::FFTs
      !%Description
      !% For solving the Poisson equation in Fourier space, and for applying the local potential
      !% in Fourier space, an auxiliary cubic mesh is built. This mesh will be larger than
      !% the circumscribed cube of the usual mesh by a factor <tt>DoubleFFTParameter</tt>. See
      !% the section that refers to Poisson equation, and to the local potential for details
      !% [the default value of two is typically good].
      !%End
      call parse_variable(namespace, 'DoubleFFTParameter', m_two, fft_alpha)
      if (fft_alpha < m_one .or. fft_alpha > m_three) then
        write(message(1), '(a,f12.5,a)') "Input: '", fft_alpha, &
          "' is not a valid DoubleFFTParameter"
        message(2) = '1.0 <= DoubleFFTParameter <= 3.0'
        call messages_fatal(2, namespace=namespace)
      end if
 
      if (space%dim /= 3 .and. fft_library == fftlib_pfft) then
        call messages_not_implemented('PFFT support for dimensionality other than 3', namespace=namespace)
      end if
 
      select case (space%dim)
 
      case (1)
        select case (this%kernel)
        case (poisson_fft_kernel_sph)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
        case (poisson_fft_kernel_nocut)
          box = der%mesh%idx%ll
        end select
 
      case (2)
        select case (this%kernel)
        case (poisson_fft_kernel_sph)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
          box(1:2) = maxval(box)
        case (poisson_fft_kernel_cyl)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
        case (poisson_fft_kernel_nocut)
          box(:) = der%mesh%idx%ll(:)
        end select
 
      case (3)
        select case (this%kernel)
        case (poisson_fft_kernel_sph)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
          box(:) = maxval(box)
        case (poisson_fft_kernel_cyl)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
          box(2) = maxval(box(2:3)) ! max of finite directions
          box(3) = maxval(box(2:3)) ! max of finite directions
        case (poisson_fft_kernel_corrected)
          box(:) = der%mesh%idx%ll(:)
        case (poisson_fft_kernel_pla, poisson_fft_kernel_nocut)
          call mesh_double_box(space, der%mesh, fft_alpha, box)
        end select
 
      end select
 
    end if
 
    ! Create the cube
    if (need_cube) then
      call cube_init(this%cube, box, namespace, space, der%mesh%spacing, &
        der%mesh%coord_system, fft_type = fft_type, &
        need_partition=.not.der%mesh%parallel_in_domains)
      call cube_init_cube_map(this%cube, der%mesh)
      if (this%cube%parallel_in_domains .and. this%method == poisson_fft) then
        call mesh_cube_parallel_map_init(this%mesh_cube_map, der%mesh, this%cube)
      end if
    end if
 
    if (this%is_dressed .and. .not. this%method == poisson_direct_sum) then
      write(message(1), '(a)')'Dressed Orbital calculation currently only working with direct sum Poisson solver.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    call poisson_kernel_init(this, namespace, space, mc, stencil)
 
    pop_sub(poisson_init)
  end subroutine poisson_init
 
  !-----------------------------------------------------------------
  subroutine poisson_end(this)
    type(poisson_t), intent(inout) :: this
 
    logical :: has_cube
 
    push_sub(poisson_end)
 
    has_cube = .false.
 
    select case (this%method)
    case (poisson_fft)
      call poisson_fft_end(this%fft_solver)
      if (this%kernel == poisson_fft_kernel_corrected) call poisson_corrections_end(this%corrector)
      has_cube = .true.
 
    case (poisson_cg_corrected, poisson_cg)
      call poisson_cg_end()
      call poisson_corrections_end(this%corrector)
 
    case (poisson_multigrid)
      call poisson_multigrid_end(this%mg)
 
    case (poisson_isf)
      call poisson_isf_end(this%isf_solver)
      has_cube = .true.
 
    case (poisson_psolver)
      call poisson_psolver_end(this%psolver_solver)
      has_cube = .true.
 
    case (poisson_fmm)
      call poisson_fmm_end(this%params_fmm)
 
    case (poisson_no)
      call poisson_no_end(this%no_solver)
 
    end select
    this%method = poisson_null
 
    if (has_cube) then
      if (this%cube%parallel_in_domains) then
        call mesh_cube_parallel_map_end(this%mesh_cube_map)
      end if
      call cube_end(this%cube)
    end if
 
    if (this%is_dressed) then
      call photon_mode_end(this%photons)
    end if
    this%is_dressed = .false.
 
    pop_sub(poisson_end)
  end subroutine poisson_end
 
  !-----------------------------------------------------------------
 
  subroutine zpoisson_solve_real_and_imag_separately(this, namespace, pot, rho, all_nodes, kernel)
    type(poisson_t),                    intent(in)    :: this
    type(namespace_t),                  intent(in)    :: namespace
    complex(real64),                    intent(inout) :: pot(:)
    complex(real64),                    intent(in)    :: rho(:)
    logical, optional,                  intent(in)    :: all_nodes
    type(fourier_space_op_t), optional, intent(in)    :: kernel
 
    real(real64), allocatable :: aux1(:), aux2(:)
    type(derivatives_t), pointer :: der
    logical :: all_nodes_value
 
 
    der => this%der
 
    push_sub(zpoisson_solve_real_and_imag_separately)
 
    call profiling_in('POISSON_RE_IM_SOLVE')
 
    if (present(kernel)) then
      assert(.not. any(abs(kernel%qq(:))>1e-8_real64))
    end if
 
    all_nodes_value = optional_default(all_nodes, this%all_nodes_default)
 
    safe_allocate(aux1(1:der%mesh%np))
    safe_allocate(aux2(1:der%mesh%np))
    ! first the real part
    aux1(1:der%mesh%np) = real(rho(1:der%mesh%np))
    aux2(1:der%mesh%np) = real(pot(1:der%mesh%np))
    call dpoisson_solve(this, namespace, aux2, aux1, all_nodes=all_nodes_value, kernel=kernel)
    pot(1:der%mesh%np)  = aux2(1:der%mesh%np)
 
    ! now the imaginary part
    aux1(1:der%mesh%np) = aimag(rho(1:der%mesh%np))
    aux2(1:der%mesh%np) = aimag(pot(1:der%mesh%np))
    call dpoisson_solve(this, namespace, aux2, aux1, all_nodes=all_nodes_value, kernel=kernel)
    pot(1:der%mesh%np) = pot(1:der%mesh%np) + m_zi*aux2(1:der%mesh%np)
 
    safe_deallocate_a(aux1)
    safe_deallocate_a(aux2)
 
    call profiling_out('POISSON_RE_IM_SOLVE')
 
    pop_sub(zpoisson_solve_real_and_imag_separately)
  end subroutine zpoisson_solve_real_and_imag_separately
 
  !-----------------------------------------------------------------
 
  subroutine zpoisson_solve(this, namespace, pot, rho, all_nodes, kernel)
    type(poisson_t),                    intent(in)    :: this
    type(namespace_t),                  intent(in)    :: namespace
    complex(real64), contiguous,                  intent(inout) :: pot(:)
    complex(real64), contiguous,                  intent(in)    :: rho(:)
    logical, optional,                  intent(in)    :: all_nodes
    type(fourier_space_op_t), optional, intent(in)    :: kernel
 
    logical :: all_nodes_value
 
    push_sub(zpoisson_solve)
 
    all_nodes_value = optional_default(all_nodes, this%all_nodes_default)
 
    assert(ubound(pot, dim = 1) == this%der%mesh%np_part .or. ubound(pot, dim = 1) == this%der%mesh%np)
    assert(ubound(rho, dim = 1) == this%der%mesh%np_part .or. ubound(rho, dim = 1) == this%der%mesh%np)
 
    assert(this%method /= poisson_null)
 
    if (this%method == poisson_fft .and. this%kernel /= poisson_fft_kernel_corrected  &
      .and. .not. this%is_dressed) then
      !The default (real) Poisson solver is used for OEP and Sternheimer calls were we do not need
      !a complex-to-xomplex FFT as these parts use the normal Coulomb potential
      if (this%cube%fft%type == fft_complex) then
        !We add the profiling here, as the other path uses dpoisson_solve
        call profiling_in('ZPOISSON_SOLVE')
        call zpoisson_fft_solve(this%fft_solver, this%der%mesh, this%cube, pot, rho, this%mesh_cube_map, kernel=kernel)
        call profiling_out('ZPOISSON_SOLVE')
      else
        call zpoisson_solve_real_and_imag_separately(this, namespace, pot, rho, all_nodes_value, kernel=kernel)
      end if
    else
      call zpoisson_solve_real_and_imag_separately(this, namespace, pot, rho, all_nodes_value, kernel = kernel)
    end if
 
    pop_sub(zpoisson_solve)
  end subroutine zpoisson_solve
 
 
  !-----------------------------------------------------------------
 
  subroutine poisson_solve_batch(this, namespace, potb, rhob, all_nodes, kernel)
    type(poisson_t),                    intent(inout) :: this
    type(namespace_t),                  intent(in)    :: namespace
    type(batch_t),                      intent(inout) :: potb
    type(batch_t),                      intent(inout) :: rhob
    logical, optional,                  intent(in)    :: all_nodes
    type(fourier_space_op_t), optional, intent(in)    :: kernel
 
    integer :: ii
 
    push_sub(poisson_solve_batch)
 
    assert(potb%nst_linear == rhob%nst_linear)
    assert(potb%type() == rhob%type())
 
    if (potb%type() == type_float) then
      do ii = 1, potb%nst_linear
        call dpoisson_solve(this, namespace, potb%dff_linear(:, ii), rhob%dff_linear(:, ii), all_nodes, kernel=kernel)
      end do
    else
      do ii = 1, potb%nst_linear
        call zpoisson_solve(this, namespace, potb%zff_linear(:, ii), rhob%zff_linear(:, ii), all_nodes, kernel=kernel)
      end do
    end if
 
    pop_sub(poisson_solve_batch)
  end subroutine poisson_solve_batch
 
  !-----------------------------------------------------------------
 
  subroutine dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel)
    type(poisson_t),                    intent(in)    :: this
    type(namespace_t),                  intent(in)    :: namespace
    real(real64), contiguous,                  intent(inout) :: pot(:)
    real(real64), contiguous,                  intent(in)    :: rho(:)
    logical, optional,                  intent(in)    :: all_nodes
    type(fourier_space_op_t), optional, intent(in)    :: kernel
 
    type(derivatives_t), pointer :: der
    real(real64), allocatable :: rho_corrected(:), vh_correction(:)
    logical               :: all_nodes_value
 
    call profiling_in('POISSON_SOLVE')
    push_sub(dpoisson_solve)
 
    der => this%der
 
    assert(ubound(pot, dim = 1) == der%mesh%np_part .or. ubound(pot, dim = 1) == der%mesh%np)
    assert(ubound(rho, dim = 1) == der%mesh%np_part .or. ubound(rho, dim = 1) == der%mesh%np)
 
    ! Check optional argument and set to default if necessary.
    all_nodes_value = optional_default(all_nodes, this%all_nodes_default)
 
    assert(this%method /= poisson_null)
 
    if (present(kernel)) then
      assert(this%method == poisson_fft)
    end if
 
    select case (this%method)
    case (poisson_direct_sum)
      if ((this%is_dressed .and. this%der%dim - 1 > 3) .or. this%der%dim > 3) then
        message(1) = "Direct sum Poisson solver only available for 1, 2, or 3 dimensions."
        call messages_fatal(1, namespace=namespace)
      end if
      call poisson_solve_direct(this, namespace, pot, rho)
 
    case (poisson_fmm)
      call poisson_fmm_solve(this%params_fmm, this%der, pot, rho)
 
    case (poisson_cg)
      call poisson_cg1(namespace, der, this%corrector, pot, rho)
 
    case (poisson_cg_corrected)
      safe_allocate(rho_corrected(1:der%mesh%np))
      safe_allocate(vh_correction(1:der%mesh%np_part))
 
      call correct_rho(this%corrector, der, rho, rho_corrected, vh_correction)
 
      pot(1:der%mesh%np) = pot(1:der%mesh%np) - vh_correction(1:der%mesh%np)
      call poisson_cg2(namespace, der, pot, rho_corrected)
      pot(1:der%mesh%np) = pot(1:der%mesh%np) + vh_correction(1:der%mesh%np)
 
      safe_deallocate_a(rho_corrected)
      safe_deallocate_a(vh_correction)
 
    case (poisson_multigrid)
      call poisson_multigrid_solver(this%mg, namespace, der, pot, rho)
 
    case (poisson_fft)
      if (this%kernel /= poisson_fft_kernel_corrected) then
        call dpoisson_fft_solve(this%fft_solver, der%mesh, this%cube, pot, rho, this%mesh_cube_map, kernel=kernel)
      else
        safe_allocate(rho_corrected(1:der%mesh%np))
        safe_allocate(vh_correction(1:der%mesh%np_part))
 
        call correct_rho(this%corrector, der, rho, rho_corrected, vh_correction)
        call dpoisson_fft_solve(this%fft_solver, der%mesh, this%cube, pot, rho_corrected, this%mesh_cube_map, &
          average_to_zero = .true., kernel=kernel)
 
        pot(1:der%mesh%np) = pot(1:der%mesh%np) + vh_correction(1:der%mesh%np)
        safe_deallocate_a(rho_corrected)
        safe_deallocate_a(vh_correction)
      end if
 
    case (poisson_isf)
      call poisson_isf_solve(this%isf_solver, der%mesh, this%cube, pot, rho, all_nodes_value)
 
 
    case (poisson_psolver)
      if (this%psolver_solver%datacode == "G") then
        ! Global version
        call poisson_psolver_global_solve(this%psolver_solver, der%mesh, this%cube, pot, rho)
      else ! "D" Distributed version
        call poisson_psolver_parallel_solve(this%psolver_solver, der%mesh, this%cube, pot, rho, this%mesh_cube_map)
      end if
 
    case (poisson_no)
      call poisson_no_solve(this%no_solver, der%mesh, pot, rho)
    end select
 
 
    ! Add extra terms for dressed interaction
    if (this%is_dressed .and. this%method /= poisson_no) then
      call photon_mode_add_poisson_terms(this%photons, der%mesh, rho, pot)
    end if
 
    pop_sub(dpoisson_solve)
    call profiling_out('POISSON_SOLVE')
  end subroutine dpoisson_solve
 
  !-----------------------------------------------------------------
  subroutine poisson_init_sm(this, namespace, space, main, der, sm, method, force_cmplx)
    type(poisson_t),             intent(inout) :: this
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    type(poisson_t),             intent(in)    :: main
    type(derivatives_t), target, intent(in)    :: der
    type(submesh_t),             intent(inout) :: sm
    integer, optional,           intent(in)    :: method
    logical, optional,           intent(in)    :: force_cmplx
 
    integer :: default_solver, idir, iter, maxl
    integer :: box(space%dim)
    real(real64)   :: qq(der%dim), threshold
 
    if (this%method /= poisson_null) return ! already initialized
 
    push_sub(poisson_init_sm)
 
    this%is_dressed = .false.
    !TODO: To be implemented as an option
    this%all_nodes_default = .false.
 
    this%nslaves = 0
    this%der => der
 
#ifdef HAVE_MPI
    this%all_nodes_default = main%all_nodes_default
#endif
 
    default_solver = poisson_direct_sum
    this%method = default_solver
    if (present(method)) this%method = method
 
    if (der%mesh%use_curvilinear) then
      call messages_not_implemented("Submesh Poisson solver with curvilinear mesh", namespace=namespace)
    end if
 
    this%kernel = poisson_fft_kernel_none
 
    select case (this%method)
    case (poisson_direct_sum)
      !Nothing to be done
 
    case (poisson_isf)
      !TODO: Add support for domain parrallelization
      assert(.not. der%mesh%parallel_in_domains)
      call submesh_get_cube_dim(sm, space, box)
      call submesh_init_cube_map(sm, space)
      call cube_init(this%cube, box, namespace, space, sm%mesh%spacing, sm%mesh%coord_system, &
        fft_type = fft_none, need_partition=.not.der%mesh%parallel_in_domains)
      call cube_init_cube_map(this%cube, sm%mesh)
      call poisson_isf_init(this%isf_solver, namespace, der%mesh, this%cube, mpi_world%comm, init_world = this%all_nodes_default)
 
    case (poisson_psolver)
      !TODO: Add support for domain parrallelization
      assert(.not. der%mesh%parallel_in_domains)
      if (this%all_nodes_default) then
        this%cube%mpi_grp = mpi_world
      else
        this%cube%mpi_grp = this%der%mesh%mpi_grp
      end if
      call submesh_get_cube_dim(sm, space, box)
      call submesh_init_cube_map(sm, space)
      call cube_init(this%cube, box, namespace, space, sm%mesh%spacing, sm%mesh%coord_system, &
        fft_type = fft_none, need_partition=.not.der%mesh%parallel_in_domains)
      call cube_init_cube_map(this%cube, sm%mesh)
      qq = m_zero
      call poisson_psolver_init(this%psolver_solver, namespace, space, this%cube, m_zero, qq, force_isolated=.true.)
      call poisson_psolver_get_dims(this%psolver_solver, this%cube)
    case (poisson_fft)
      !Here we impose zero boundary conditions
      this%kernel = poisson_fft_kernel_sph
      !We need to parse this, in case this routine is called before poisson_init
      call parse_variable(namespace, 'FFTLibrary', fftlib_fftw, fft_default_lib)
 
      call submesh_get_cube_dim(sm, space, box)
      call submesh_init_cube_map(sm, space)
      !We double the size of the cell
      !Maybe the factor of two should be controlled as a variable
      do idir = 1, space%dim
        box(idir) = nint(m_two * (box(idir) - 1)) + 1
      end do
      if (optional_default(force_cmplx, .false.)) then
        call cube_init(this%cube, box, namespace, space, sm%mesh%spacing, sm%mesh%coord_system, &
          fft_type = fft_complex, need_partition=.not.der%mesh%parallel_in_domains)
      else
        call cube_init(this%cube, box, namespace, space, sm%mesh%spacing, sm%mesh%coord_system, &
          fft_type = fft_real, need_partition=.not.der%mesh%parallel_in_domains)
      end if
      call poisson_fft_init(this%fft_solver, namespace, space, this%cube, this%kernel)
    case (poisson_cg)
      call parse_variable(namespace, 'PoissonSolverMaxMultipole', 4, maxl)
      write(message(1),'(a,i2)')'Info: Boundary conditions fixed up to L =',  maxl
      call messages_info(1, namespace=namespace)
      call parse_variable(namespace, 'PoissonSolverMaxIter', 500, iter)
      call parse_variable(namespace, 'PoissonSolverThreshold', 1.0e-6_real64, threshold)
      call poisson_corrections_init(this%corrector, namespace, space, maxl, this%der%mesh)
      call poisson_cg_init(threshold, iter)
    end select
 
    pop_sub(poisson_init_sm)
  end subroutine poisson_init_sm
 
  !-----------------------------------------------------------------
  subroutine poisson_test(this, space, mesh, latt, namespace, repetitions)
    type(poisson_t),         intent(in) :: this
    class(space_t),          intent(in) :: space
    class(mesh_t),           intent(in) :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    type(namespace_t),       intent(in) :: namespace
    integer,                 intent(in) :: repetitions
 
    real(real64), allocatable :: rho(:), vh(:), vh_exact(:), xx(:, :), xx_per(:)
    real(real64) :: alpha, beta, rr, delta, ralpha, hartree_nrg_num, &
      hartree_nrg_analyt, lcl_hartree_nrg
    real(real64) :: total_charge
    integer :: ip, ierr, iunit, nn, n_gaussians, itime, icell
    real(real64) :: threshold
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(poisson_test)
 
    if (mesh%box%dim == 1) then
      call messages_not_implemented('Poisson test for 1D case', namespace=namespace)
    end if
 
    !%Variable PoissonTestPeriodicThreshold
    !%Type float
    !%Default 1e-5
    !%Section Hamiltonian::Poisson
    !%Description
    !% This threshold determines the accuracy of the periodic copies of
    !% the Gaussian charge distribution that are taken into account when
    !% computing the analytical solution for periodic systems.
    !% Be aware that the default leads to good results for systems
    !% that are periodic in 1D - for 3D it is very costly because of the
    !% large number of copies needed.
    !%End
    call parse_variable(namespace, 'PoissonTestPeriodicThreshold', 1e-5_real64, threshold)
 
    ! Use two gaussians with different sign
    n_gaussians = 2
 
    safe_allocate(     rho(1:mesh%np))
    safe_allocate(      vh(1:mesh%np))
    safe_allocate(vh_exact(1:mesh%np))
    safe_allocate(xx(1:space%dim, 1:n_gaussians))
    safe_allocate(xx_per(1:space%dim))
 
    rho = m_zero
    vh = m_zero
    vh_exact = m_zero
 
    alpha = m_four*mesh%spacing(1)
    write(message(1),'(a,f14.6)')  "Info: The alpha value is ", alpha
    write(message(2),'(a)')        "      Higher values of alpha lead to more physical densities and more reliable results."
    call messages_info(2, namespace=namespace)
    beta = m_one / (alpha**space%dim * sqrt(m_pi)**space%dim)
 
    write(message(1), '(a)') 'Building the Gaussian distribution of charge...'
    call messages_info(1, namespace=namespace)
 
    ! Set the centers of the Gaussians by hand
    xx(1, 1) = m_one
    xx(2, 1) = -m_half
    if (space%dim == 3) xx(3, 1) = m_two
    xx(1, 2) = -m_two
    xx(2, 2) = m_zero
    if (space%dim == 3) xx(3, 2) = -m_one
    xx = xx * alpha
 
    ! Density as sum of Gaussians
    rho = m_zero
    do nn = 1, n_gaussians
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, rr, origin = xx(:, nn))
        rho(ip) = rho(ip) + (-1)**nn * beta*exp(-(rr/alpha)**2)
      end do
    end do
 
    total_charge = dmf_integrate(mesh, rho)
 
    write(message(1), '(a,e23.16)') 'Total charge of the Gaussian distribution', total_charge
    call messages_info(1, namespace=namespace)
 
    write(message(1), '(a)') 'Computing exact potential.'
    call messages_info(1, namespace=namespace)
 
    ! This builds analytically its potential
    vh_exact = m_zero
    latt_iter = lattice_iterator_t(latt, m_one/threshold)
    do nn = 1, n_gaussians
      ! sum over all periodic copies for each Gaussian
      write(message(1), '(a,i2,a,i9,a)') 'Computing Gaussian ', nn, ' for ', latt_iter%n_cells, ' periodic copies.'
      call messages_info(1, namespace=namespace)
 
      do icell = 1, latt_iter%n_cells
        xx_per = xx(:, nn) + latt_iter%get(icell)
        !$omp parallel do private(rr, ralpha)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin=xx_per)
          select case (space%dim)
          case (3)
            if (rr > r_small) then
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * loct_erf(rr/alpha)/rr
            else
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * (m_two/sqrt(m_pi))/alpha
            end if
          case (2)
            ralpha = rr**2/(m_two*alpha**2)
            if (ralpha < 100.0_real64) then
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (m_pi)**(m_three*m_half) * alpha * exp(-rr**2/(m_two*alpha**2)) * &
                loct_bessel_in(0, rr**2/(m_two*alpha**2))
            else
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (m_pi)**(m_three*m_half) * alpha * &
                (m_one/sqrt(m_two*m_pi*ralpha))
            end if
          end select
        end do
      end do
    end do
 
    ! This calculates the numerical potential
    do itime = 1, repetitions
      call dpoisson_solve(this, namespace, vh, rho)
    end do
 
    ! Output results
    iunit = io_open("hartree_results", namespace, action='write')
    delta = dmf_nrm2(mesh, vh-vh_exact)
    write(iunit, '(a,e23.16)') 'Hartree test (abs.) = ', delta
    delta = delta/dmf_nrm2(mesh, vh_exact)
    write(iunit, '(a,e23.16)') 'Hartree test (rel.) = ', delta
 
    ! Calculate the numerical Hartree energy (serially)
    lcl_hartree_nrg = m_zero
    do ip = 1, mesh%np
      lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh(ip)
    end do
    lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3)/m_two
#ifdef HAVE_MPI
    call mpi_reduce(lcl_hartree_nrg, hartree_nrg_num, 1, &
      mpi_double_precision, mpi_sum, 0, mpi_comm_world, mpi_err)
    if (mpi_err /= 0) then
      write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90"
      call messages_warning(1, namespace=namespace)
    end if
#else
    hartree_nrg_num = lcl_hartree_nrg
#endif
 
    ! Calculate the analytical Hartree energy (serially, discrete - not exactly exact)
    lcl_hartree_nrg = m_zero
    do ip = 1, mesh%np
      lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh_exact(ip)
    end do
    lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3) / m_two
#ifdef HAVE_MPI
    call mpi_reduce(lcl_hartree_nrg, hartree_nrg_analyt, 1, &
      mpi_double_precision, mpi_sum, 0, mpi_comm_world, mpi_err)
    if (mpi_err /= 0) then
      write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90"
      call messages_warning(1, namespace=namespace)
    end if
#else
    hartree_nrg_analyt = lcl_hartree_nrg
#endif
 
    write(iunit, '(a)')
 
    if (mpi_world%rank == 0) then
      write(iunit,'(a,e23.16)') 'Hartree Energy (numerical)  =', hartree_nrg_num, 'Hartree Energy (analytical) =',hartree_nrg_analyt
    end if
 
    call io_close(iunit)
 
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_rho", namespace, space, &
      mesh, rho, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_exact", namespace, space, &
      mesh, vh_exact, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_numerical", namespace, space, &
      mesh, vh, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_rho", namespace, space, &
      mesh, rho, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_exact", namespace, space, &
      mesh, vh_exact, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_numerical", namespace, space, &
      mesh, vh, unit_one, ierr)
    ! not dimensionless, but no need for unit conversion for a test routine
 
    safe_deallocate_a(rho)
    safe_deallocate_a(vh)
    safe_deallocate_a(vh_exact)
    safe_deallocate_a(xx)
 
    pop_sub(poisson_test)
  end subroutine poisson_test
 
  ! -----------------------------------------------------------------
 
  logical pure function poisson_solver_is_iterative(this) result(iterative)
    type(poisson_t), intent(in) :: this
 
    iterative = this%method == poisson_cg .or. this%method == poisson_cg_corrected .or. this%method == poisson_multigrid
  end function poisson_solver_is_iterative
 
  ! -----------------------------------------------------------------
 
  logical pure function poisson_is_multigrid(this) result(is_multigrid)
    type(poisson_t), intent(in) :: this
 
    is_multigrid = (this%method == poisson_multigrid)
 
  end function poisson_is_multigrid
 
  ! -----------------------------------------------------------------
 
  logical pure function poisson_solver_has_free_bc(this) result(free_bc)
    type(poisson_t), intent(in) :: this
 
    free_bc = .true.
 
    if (this%method == poisson_fft .and. &
      this%kernel /= poisson_fft_kernel_sph .and. this%kernel /= poisson_fft_kernel_corrected) then
      free_bc = .false.
    end if
 
  end function poisson_solver_has_free_bc
 
  !-----------------------------------------------------------------
 
  integer pure function poisson_get_solver(this) result (solver)
    type(poisson_t), intent(in) :: this
 
    solver = this%method
  end function poisson_get_solver
 
  !-----------------------------------------------------------------
 
  subroutine poisson_async_init(this, mc)
    type(poisson_t), intent(inout) :: this
    type(multicomm_t), intent(in)  :: mc
 
    push_sub(poisson_async_init)
 
#ifdef HAVE_MPI
    if (multicomm_have_slaves(mc)) then
 
      call mpi_grp_init(this%local_grp, mc%group_comm(p_strategy_states))
 
      this%root = (this%local_grp%rank == 0)
 
      this%intercomm = mc%slave_intercomm
      call mpi_comm_remote_size(this%intercomm, this%nslaves, mpi_err)
 
    end if
#endif
 
    pop_sub(poisson_async_init)
 
  end subroutine poisson_async_init
 
  !-----------------------------------------------------------------
 
  subroutine poisson_async_end(this, mc)
    type(poisson_t), intent(inout) :: this
    type(multicomm_t), intent(in)  :: mc
 
#ifdef HAVE_MPI
    integer :: islave
#endif
 
    push_sub(poisson_async_end)
 
#ifdef HAVE_MPI
    if (multicomm_have_slaves(mc)) then
 
      ! send the finish signal
      do islave = this%local_grp%rank, this%nslaves - 1, this%local_grp%size
        call mpi_send(m_one, 1, mpi_double_precision, islave, cmd_finish, this%intercomm, mpi_err)
      end do
 
    end if
#endif
 
    pop_sub(poisson_async_end)
 
  end subroutine poisson_async_end
 
  !-----------------------------------------------------------------
 
  subroutine poisson_slave_work(this, namespace)
    type(poisson_t),   intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
 
#ifdef HAVE_MPI
    real(real64), allocatable :: rho(:), pot(:)
    logical :: done
    type(mpi_status) :: status
    integer :: bcast_root
 
    push_sub(poisson_slave_work)
    call profiling_in("SLAVE_WORK")
 
    safe_allocate(rho(1:this%der%mesh%np))
    safe_allocate(pot(1:this%der%mesh%np))
    done = .false.
 
    do while(.not. done)
 
      call profiling_in("SLAVE_WAIT")
      call mpi_recv(rho(1), this%der%mesh%np, mpi_double_precision, mpi_any_source, mpi_any_tag, this%intercomm, status, mpi_err)
      call profiling_out("SLAVE_WAIT")
 
      ! The tag of the message tells us what we have to do.
      select case (status%MPI_TAG)
 
      case (cmd_finish)
        done = .true.
 
      case (cmd_poisson_solve)
        call dpoisson_solve(this, namespace, pot, rho)
 
        call profiling_in("SLAVE_BROADCAST")
        bcast_root = mpi_proc_null
        if (this%root) bcast_root = mpi_root
        call mpi_bcast(pot(1), this%der%mesh%np, mpi_double_precision, bcast_root, this%intercomm, mpi_err)
        call profiling_out("SLAVE_BROADCAST")
 
      end select
 
    end do
 
    safe_deallocate_a(pot)
    safe_deallocate_a(rho)
 
    call profiling_out("SLAVE_WORK")
    pop_sub(poisson_slave_work)
#endif
  end subroutine poisson_slave_work
 
  !----------------------------------------------------------------
 
  logical pure function poisson_is_async(this) result(async)
    type(poisson_t),  intent(in) :: this
 
    async = (this%nslaves > 0)
 
  end function poisson_is_async
 
  !----------------------------------------------------------------
 
  subroutine poisson_build_kernel(this, namespace, space, coulb, qq, alpha, beta, mu, singul)
    type(poisson_t),          intent(in)    :: this
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(fourier_space_op_t), intent(inout) :: coulb
    real(real64),             intent(in)    :: qq(:)
    real(real64),             intent(in)    :: alpha, beta, mu
    real(real64),   optional, intent(in)    :: singul
 
    real(real64), parameter :: vanishing_q = 1.0e-5_real64
    logical :: reinit
 
    push_sub(poisson_build_kernel)
 
    if (space%is_periodic()) then
      assert(ubound(qq, 1) >= space%periodic_dim)
      assert(this%method == poisson_fft)
    end if
 
    if (mu > m_epsilon) then
      if (this%method /= poisson_fft) then
        write(message(1),'(a)') "Poisson solver with range separation is only implemented with FFT."
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    !We only reinitialize the poisson solver if needed
    reinit = .false.
    if (allocated(coulb%qq)) then
      reinit = any(abs(coulb%qq(1:space%periodic_dim) - qq(1:space%periodic_dim)) > m_epsilon)
    end if
    reinit = reinit .or. (abs(coulb%mu-mu) > m_epsilon .and. mu > m_epsilon)
    reinit = reinit .or. (abs(coulb%alpha-alpha) > m_epsilon .and. alpha > m_epsilon)
    reinit = reinit .or. (abs(coulb%beta-beta) > m_epsilon .and. beta > m_epsilon)
 
 
    if (reinit) then
      !TODO: this should be a select case supporting other kernels.
      ! This means that we need an abstract object for kernels.
      select case (this%method)
      case (poisson_fft)
        ! Check that we are consistent: the Poisson solver supports must return 1 here
        assert(is_close(poisson_get_full_range_weight(this, alpha, beta, mu), m_one))
 
        call fourier_space_op_end(coulb)
 
        safe_allocate(coulb%qq(1:space%dim))
        coulb%qq(1:space%periodic_dim) = qq(1:space%periodic_dim)
        coulb%qq(space%periodic_dim+1:space%dim) = vanishing_q
        !We must define the singularity if we specify a q vector and we do not use the short-range Coulomb potential
        coulb%singularity = optional_default(singul, m_zero)
        coulb%mu = mu
        coulb%alpha = alpha
        coulb%beta = beta
        call poisson_fft_get_kernel(namespace, space, this%cube, coulb, this%kernel, &
          this%poisson_soft_coulomb_param)
      case default
        call messages_not_implemented("poisson_build_kernel with other methods than FFT", namespace=namespace)
      end select
    end if
 
    pop_sub(poisson_build_kernel)
  end subroutine poisson_build_kernel
 
  !----------------------------------------------------------------
  real(real64) function poisson_get_full_range_weight(this, alpha, beta, mu) result(weight)
    type(poisson_t),          intent(in) :: this
    real(real64),             intent(in) :: alpha, beta, mu
 
    select case (this%method)
    case (poisson_fft)
      weight = m_one
    case default
      if(mu < m_epsilon) then
        weight = alpha
      else if(alpha > m_epsilon .and. beta < m_epsilon) then
        weight = alpha
      else if(alpha < m_epsilon .and. beta > m_epsilon) then
        weight = beta
      else
        assert(.false.)
      end if
    end select
  end function poisson_get_full_range_weight
 
#include "poisson_init_inc.F90"
#include "poisson_direct_inc.F90"
#include "poisson_direct_sm_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "poisson_inc.F90"
#include "undef.F90"
#include "complex.F90"
#include "poisson_inc.F90"
 
end module poisson_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: