poisson_psolver.F90

Go to the documentation of this file.

!! Copyright (C) 2013 J. Alberdi-Rodriguez
!!
!! 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_psolver_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use debug_oct_m
  use fourier_space_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_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 profiling_oct_m
  use space_oct_m
  use submesh_oct_m
  ! Support for PSolver from BigDFT
#ifdef HAVE_PSOLVER
  use mesh_function_oct_m
  use poisson_solver
  use dictionaries, dict_set => set
  use yaml_output, only: yaml_map
  use wrapper_mpi, only: mpi_environment, mpi_environment_set
#endif
#ifdef HAVE_PSOLVER_NEW_API
  use at_domain
#endif
 
  implicit none
 
  private
  public ::                         &
    poisson_psolver_t,              &
    poisson_psolver_init,           &
    poisson_psolver_end,            &
    poisson_psolver_reinit,         &
    poisson_psolver_global_solve,   &
    poisson_psolver_parallel_solve, &
    poisson_psolver_get_dims
 
  type poisson_psolver_t
    private
    type(fourier_space_op_t) :: coulb
#ifdef HAVE_PSOLVER
    type(coulomb_operator) :: kernel
    type(dictionary), pointer :: inputs
#endif
 
    character(len = 1) :: geocode  = "F"
    character(len = 1), public :: datacode = "G"
 
    integer :: isf_order
    integer :: localnscatterarr(5)
    real(real64)          :: offset
  end type poisson_psolver_t
 
#ifdef HAVE_PSOLVER
  logical, save :: flib_initialized = .false.
#endif
 
  real(real64), parameter :: TOL_VANISHING_Q = 1e-6_real64
 
contains
 
  !-----------------------------------------------------------------
  subroutine poisson_psolver_init(this, namespace, space, cube, mu, qq, force_isolated)
    type(poisson_psolver_t), intent(out)   :: this
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(cube_t),            intent(inout) :: cube
    real(real64),            intent(in)    :: mu
    real(real64),            intent(in)    :: qq(:)
    logical, optional,       intent(in)    :: force_isolated
 
    logical data_is_parallel
#ifdef HAVE_PSOLVER
    real(real64) :: alpha, beta, gamma
    real(real64) :: modq2
    type(mpi_environment) mpi_env
#ifdef HAVE_PSOLVER_NEW_API
    type(domain) :: dom
#endif
#endif
 
    push_sub(poisson_psolver_init)
 
#ifdef HAVE_PSOLVER
    if (.not. flib_initialized) then
      call f_lib_initialize()
      flib_initialized = .true.
    end if
    call dict_init(this%inputs)
#endif
 
    if (optional_default(force_isolated, .false.)) then
      this%geocode = "F"
    else
      select case (space%periodic_dim)
      case (0)
        ! Free BC
        this%geocode = "F"
      case (1)
        ! Wire BC
        this%geocode = "W"
        call messages_not_implemented("PSolver support for 1D periodic boundary conditions", namespace=namespace)
      case (2)
        ! Surface BC
        this%geocode = "S"
        call messages_not_implemented("PSolver support for 2D periodic boundary conditions", namespace=namespace)
      case (3)
        ! Periodic BC
        this%geocode = "P"
        call messages_experimental("PSolver support for 3D periodic boundary conditions", namespace=namespace)
      end select
    end if
 
#ifdef HAVE_PSOLVER
    !Verbosity switch
    call dict_set(this%inputs//'setup'//'verbose', .false.)
    !Order of the Interpolating Scaling Function family
    call dict_set(this%inputs//'kernel'//'isf_order', 16)
    !Mu screening parameter
    call dict_set(this%inputs//'kernel'//'screening', mu)
    !Calculation of the stress tensor
    call dict_set(this%inputs//'kernel'//'stress_tensor', .false.)
#else
    this%isf_order = 16
#endif
 
    !%Variable PoissonSolverPSolverParallelData
    !%Type logical
    !%Section Hamiltonian::Poisson::PSolver
    !%Default yes
    !%Description
    !% Indicates whether data is partitioned within the PSolver library.
    !% If data is distributed among processes, Octopus uses parallel data-structures
    !% and, thus, less memory.
    !% If "yes", data is parallelized. The <i>z</i>-axis of the input vector
    !% is split among the MPI processes.
    !% If "no", entire input and output vector is saved in all the MPI processes.
    !% If k-points parallelization is used, "no" must be selected.
    !%End
    call parse_variable(namespace, 'PoissonSolverPSolverParallelData', .true., data_is_parallel)
    if (.not. cube%parallel_in_domains) then
      data_is_parallel = .false.
    end if
 
    call messages_obsolete_variable(namespace, 'PoissonSolverISFParallelData', 'PoissonSolverPSolverParallelData')
 
#ifdef HAVE_PSOLVER
    if (data_is_parallel) then
      call dict_set(this%inputs//'setup'//'global_data', .false.)
    else
      call dict_set(this%inputs//'setup'//'global_data', .true.)
    end if
#else
    if (data_is_parallel) then
      this%datacode = "D"
    else
      this%datacode = "G"
    end if
#endif
 
#ifdef HAVE_PSOLVER
    call dict_set(this%inputs//'setup'//'verbose', debug%info)
 
    alpha = cube%latt%alpha*m_pi/(180.0_real64)
    beta  = cube%latt%beta*m_pi/(180.0_real64)
    gamma = cube%latt%gamma*m_pi/(180.0_real64)
 
    ! Previously, pkernel_init set the communicator used within PSolver to comm_world.
    ! This can be overwritten by passing an optional argument of type(mpi_environment)
    ! to pkernel_init(). This data type is defined within the wrapper_MPI module of
    ! the Futile library. Futile is a prerequisit for PSolver.
 
    ! TODO: check that cube%mpi_grp corresponds to the correct mpi group, when parallelizing, e.g., over systems !!
 
    call mpi_environment_set(mpi_env, cube%mpi_grp%rank, cube%mpi_grp%size, cube%mpi_grp%comm%MPI_VAL, cube%mpi_grp%size)
 
#ifdef HAVE_PSOLVER_NEW_API
    dom=domain_new(units=atomic_units, bc=geocode_to_bc_enum(this%geocode),&
      alpha_bc=alpha, beta_ac=beta, gamma_ab=gamma, acell=cube%rs_n_global*cube%spacing)
 
 
    this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs, dom, cube%rs_n_global, &
      cube%spacing, alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma, mpi_env = mpi_env)
#else
    this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs, this%geocode, cube%rs_n_global, &
      cube%spacing, alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma, mpi_env = mpi_env)
#endif
 
    call pkernel_set(this%kernel, verbose=debug%info)
 
    !G=0 component
    modq2 = sum(qq(1:space%periodic_dim)**2)
    if (modq2 > m_epsilon) then
      this%offset = m_one/modq2
    else
      this%offset = m_zero
    end if
 
    !Screened coulomb potential (erfc function)
    if (mu > m_epsilon) then
      if (modq2 > m_epsilon) then
        this%offset = this%offset*(m_one - exp(-modq2/((m_two*mu)**2)))
      else
        !Analytical limit of 1/|q|^2*(1-exp(-|q|^2/4mu^2))
        this%offset = m_one/((m_two*mu)**2)
      end if
    end if
    this%offset = this%offset*m_four*m_pi
#endif
 
    pop_sub(poisson_psolver_init)
  end subroutine poisson_psolver_init
 
  !-----------------------------------------------------------------
  subroutine poisson_psolver_end(this)
    type(poisson_psolver_t), intent(inout) :: this
 
    push_sub(poisson_psolver_end)
 
#ifdef HAVE_PSOLVER
    call pkernel_free(this%kernel)
    call f_lib_finalize()
    call dict_free(this%inputs)
#endif
 
    pop_sub(poisson_psolver_end)
  end subroutine poisson_psolver_end
 
  !-----------------------------------------------------------------
  subroutine poisson_psolver_reinit(this, space, cube, coulb, qq_in)
    type(poisson_psolver_t),  intent(inout) :: this
    class(space_t),           intent(in)    :: space
    type(cube_t),             intent(inout) :: cube
    type(fourier_space_op_t), intent(inout) :: coulb
    real(real64),             intent(in)    :: qq_in(:)
 
    real(real64) :: alpha, beta, gamma
    real(real64) :: qq_abs(1:space%dim), qq(1:space%dim)
    real(real64) :: modq2
    integer :: idim
#ifdef HAVE_PSOLVER_NEW_API
    type(domain) :: dom
#endif
 
    push_sub(poisson_psolver_reinit)
 
    ! We only support short range only or long-range only at the moment
    if(coulb%mu > m_epsilon) then
      assert(coulb%alpha < m_epsilon .or. coulb%beta < m_epsilon)
    end if
 
#ifdef HAVE_PSOLVER
    call pkernel_free(this%kernel)
#endif
 
    !We might change the cell angles
    alpha = cube%latt%alpha*m_pi/(180.0_real64)
    beta  = cube%latt%beta*m_pi/(180.0_real64)
    gamma = cube%latt%gamma*m_pi/(180.0_real64)
 
#ifdef HAVE_PSOLVER
    call dict_set(this%inputs//'kernel'//'screening', coulb%mu)
#ifdef HAVE_PSOLVER_NEW_API
    dom=domain_new(units=atomic_units, bc=geocode_to_bc_enum(this%geocode),&
      alpha_bc=alpha, beta_ac=beta, gamma_ab=gamma, acell=cube%rs_n_global*cube%spacing)
 
    this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs,&
      dom,cube%rs_n_global, cube%spacing, &
      alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma)
#else
    this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs,&
      this%geocode, cube%rs_n_global, cube%spacing, &
      alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma)
#endif
    call pkernel_set(this%kernel, verbose=debug%info)
#endif
 
    !G=0 component
    !We remove potential umklapp
    do idim = 1, space%periodic_dim
      qq(idim) = qq_in(idim) - anint(qq_in(idim) + m_half*1e-8_real64)
    end do
    qq(space%periodic_dim + 1:space%dim) = m_zero
    call kpoints_to_absolute(cube%latt, qq, qq_abs)
    modq2 = norm2(qq_abs)
    if (modq2 > tol_vanishing_q) then
      this%offset = m_four*m_pi/modq2
    else
      this%offset = coulb%singularity
    end if
 
    !Screened coulomb potential (erfc function)
    if (coulb%mu > m_epsilon) then
      if (modq2 > tol_vanishing_q) then
        this%offset = this%offset*(m_one - exp(-modq2/((m_two*coulb%mu)**2)))
      else
        !Analytical limit of 1/|q|^2*(1-exp(-|q|^2/4mu^2))
        this%offset = m_four*m_pi/((m_two*coulb%mu)**2)
      end if
    end if
 
    pop_sub(poisson_psolver_reinit)
  end subroutine poisson_psolver_reinit
 
  !-----------------------------------------------------------------
  subroutine poisson_psolver_parallel_solve(this, mesh, cube, pot, rho,  mesh_cube_map)
    type(poisson_psolver_t),        intent(in), target :: this
    type(mesh_t),                   intent(in)         :: mesh
    type(cube_t),                   intent(in)         :: cube
    real(real64),                   intent(out)        :: pot(:)
    real(real64),                   intent(in)         :: rho(:)
    type(mesh_cube_parallel_map_t), intent(in)         :: mesh_cube_map
 
    type(cube_function_t) :: cf
 
#ifdef HAVE_PSOLVER
    type(coulomb_operator), pointer :: kernel_pointer
    real(real64) :: hartree_energy
    real(real64) :: offset
    !                                 To be used only in the periodic case, ignored for other boundary conditions.
#endif
 
    real(real64), allocatable :: pot_ion(:,:,:)
    character(len=3) :: quiet
 
    push_sub(poisson_psolver_parallel_solve)
 
    call dcube_function_alloc_rs(cube, cf)
    call dmesh_to_cube_parallel(mesh, rho, cube, cf, mesh_cube_map)
    safe_allocate(pot_ion(1:cube%rs_n(1),1:cube%rs_n(2),1:cube%rs_n(3)))
 
    quiet = "YES"
    if (debug%info) quiet = "NO"
 
#ifdef HAVE_PSOLVER
    !The offset is the integral over space of the potential
    !this%offset is the G=0 component of the (screened) Coulomb potential
    !The G=0 component of the Hartree therefore needs to be
    ! multiplied by the G=0 component of the density
    if (this%offset > m_epsilon) then
      offset = this%offset*dmf_integrate(mesh,rho)
    end if
#endif
 
    call profiling_in("PSOLVER_LIBRARY")
#ifdef HAVE_PSOLVER
    kernel_pointer => this%kernel
    call h_potential(this%datacode, kernel_pointer, cf%dRS, pot_ion, hartree_energy, offset, .false., quiet = quiet)
#endif
    call profiling_out("PSOLVER_LIBRARY")
    safe_deallocate_a(pot_ion)
 
    call dcube_to_mesh_parallel(cube, cf, mesh, pot, mesh_cube_map)
 
    call dcube_function_free_rs(cube, cf)
 
    pop_sub(poisson_psolver_parallel_solve)
  end subroutine poisson_psolver_parallel_solve
 
  !-----------------------------------------------------------------
  subroutine poisson_psolver_global_solve(this, mesh, cube, pot, rho, sm)
    type(poisson_psolver_t), intent(in), target :: this
    type(mesh_t),            intent(in)         :: mesh
    type(cube_t),            intent(in)         :: cube
    real(real64), contiguous,       intent(out)        :: pot(:)
    real(real64), contiguous,       intent(in)         :: rho(:)
    type(submesh_t),     optional,  intent(in)    :: sm
 
    character(len=3) :: quiet
    type(cube_function_t) :: cf
 
#ifdef HAVE_PSOLVER
    type(coulomb_operator), pointer :: kernel_pointer
    real(real64) :: hartree_energy
    real(real64) :: offset
    !                                 To be used only in the periodic case, ignored for other boundary conditions.
#endif
 
    real(real64), allocatable ::  pot_ion(:,:,:)
 
    push_sub(poisson_psolver_global_solve)
 
    call dcube_function_alloc_rs(cube, cf)
 
    if (present(sm)) then
      call dsubmesh_to_cube(sm, rho, cube, cf)
    else
      call dmesh_to_cube(mesh, rho, cube, cf)
    end if
 
    safe_allocate(pot_ion(1:cube%rs_n(1), 1:cube%rs_n(2), 1:cube%rs_n(3)))
 
    quiet = "YES"
    if (debug%info) quiet = "NO"
 
#ifdef HAVE_PSOLVER
    !The offset is the integral over space of the potential
    !this%offset is the G=0 component of the (screened) Coulomb potential
    !The G=0 component of the Hartree therefore needs to be
    ! multiplied by the G=0 component of the density
    if (this%offset > m_epsilon) then
      offset = this%offset*dmf_integrate(mesh, rho)
    end if
#endif
 
    call profiling_in("PSOLVER_LIBRARY")
 
#ifdef HAVE_PSOLVER
    kernel_pointer => this%kernel
    call h_potential(this%datacode, kernel_pointer, &
      cf%dRS,  pot_ion, hartree_energy, offset, .false., &
      quiet = quiet) !optional argument
#endif
    safe_deallocate_a(pot_ion)
 
    if (present(sm)) then
      call dcube_to_submesh(cube, cf, sm, pot)
    else
      call dcube_to_mesh(cube, cf, mesh, pot)
    end if
 
    call dcube_function_free_rs(cube, cf)
 
    pop_sub(poisson_psolver_global_solve)
  end subroutine poisson_psolver_global_solve
 
  subroutine poisson_psolver_get_dims(this, cube)
    type(poisson_psolver_t), intent(inout) :: this
    type(cube_t),            intent(inout) :: cube
 
#ifdef HAVE_PSOLVER
 
    integer :: n3d
    integer :: n3p
    integer :: n3pi
    integer :: i3xcsh
    integer :: i3s
 
    logical :: use_gradient
    logical :: use_wb_corr
#endif
 
    push_sub(poisson_psolver_get_dims)
 
    ! Get the dimensions of the cube
#ifdef HAVE_PSOLVER
    use_gradient = .false.
    use_wb_corr = .false.
    call ps_dim4allocation(this%geocode, this%datacode, cube%mpi_grp%rank, cube%mpi_grp%size, &
      cube%rs_n_global(1), cube%rs_n_global(2), cube%rs_n_global(3), &
      use_gradient, use_wb_corr, &
      0, n3d, n3p, n3pi, i3xcsh, i3s)
    this%localnscatterarr(:) = (/ n3d, n3p, n3pi, i3xcsh, i3s /)
#endif
    cube%rs_n(1:2)      = cube%rs_n_global(1:2)
    cube%rs_n(3)        = this%localnscatterarr(1)
    cube%rs_istart(1:2) = 1
    cube%rs_istart(3)   = this%localnscatterarr(5)
 
    !! With PSolver we don`t care about the Fourier space and its dimensions
    !! We`ll put as in RS
    cube%fs_n_global(1) = cube%rs_n_global(1)
    cube%fs_n_global(2) = cube%rs_n_global(2)
    cube%fs_n_global(3) = cube%rs_n_global(3)
    cube%fs_n(1:2)      = cube%rs_n_global(1:2)
    cube%fs_n(3)        = this%localnscatterarr(1)
    cube%fs_istart(1:2) = 1
    cube%fs_istart(3)   = this%localnscatterarr(5)
 
    pop_sub(poisson_psolver_get_dims)
  end subroutine poisson_psolver_get_dims
 
end module poisson_psolver_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: