exchange_operator.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2018 M. Marques, A. Castro, A. Rubio, G. Bertsch, N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module exchange_operator_oct_m
  use accel_oct_m
  use accel_blas_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use blas_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use fft_oct_m,           only: fftlib_accel
  use fourier_space_oct_m
  use global_oct_m
  use grid_oct_m
  use isdf_options_oct_m
  use kpoints_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use mesh_batch_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use par_vec_oct_m
  use parser_oct_m
  use phase_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use ring_pattern_oct_m
  use singularity_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_parallel_oct_m
  use stencil_oct_m
  use symmetries_oct_m
  use symmetrizer_oct_m,   only: dsymmetrizer_apply_batch, zsymmetrizer_apply_batch, &
    ssymmetrizer_apply_batch, csymmetrizer_apply_batch
  use states_elec_all_to_all_communications_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
  use xc_cam_oct_m
 
  implicit none
 
  private
  public ::                          &
    exchange_operator_t,             &
    exchange_operator_init,          &
    exchange_operator_reinit,        &
    exchange_operator_end,           &
    dexchange_operator_single,       &
    zexchange_operator_single,       &
    dexchange_operator_apply,        &
    zexchange_operator_apply,        &
    dexchange_operator_hartree_apply,&
    zexchange_operator_hartree_apply,&
    exchange_operator_rdmft_occ_apply,&
    dexchange_operator_compute_potentials, &
    zexchange_operator_compute_potentials, &
    dexchange_operator_commute_r, &
    zexchange_operator_commute_r, &
    dexchange_operator_ace, &
    zexchange_operator_ace, &
    dexchange_operator_compute_ex, &
    zexchange_operator_compute_ex
 
  type ace_t
    integer :: nst
    real(real64),    allocatable :: dchi(:,:,:,:)
    complex(real64), allocatable :: zchi(:,:,:,:)
    type(wfs_elec_t), allocatable :: chib(:)
    !                                                    one batch per k-point (kpt%start:kpt%end), used to
    !                                                    apply the operator with GEMMs while keeping the
    !                                                    wavefunctions resident on the GPU. dchi/zchi above
    !                                                    remain the reference host copy.
  contains
    procedure :: init => ace_init
    procedure :: end  => ace_end
    procedure :: write_info => ace_write_info
  end type ace_t
 
 
  type exchange_operator_t
    type(states_elec_t), public, pointer :: st => null()
    type(xc_cam_t)       :: cam
    type(poisson_t)      :: psolver
    type(singularity_t)  :: singul
    logical              :: useACE
    logical              :: with_isdf
    logical              :: compute_single_precision = .false. 
    type(ACE_t)          :: ace
    type(isdf_options_t) :: isdf
  contains
    procedure            :: write_info => exchange_operator_write_info
  end type exchange_operator_t
 
 
  type(fourier_space_op_t) :: coulb
  !                                    Saved as we avoid then to recompute it,
  !                                    for instance in the case of CAM functionals in isolated systems
 
  type(states_elec_t) :: st_single_shadow, xst_single_shadow
  logical :: single_shadows_init = .false.
  integer :: single_shadows_np = 0
 
  real(real64), parameter, private :: TOL_EXX_WEIGHT = 1.0e-3_real64 
 
contains
 
  subroutine ace_init(this, namespace, st)
    class(ACE_t),        intent(inout) :: this
    type(namespace_t),   intent(in)    :: namespace
    type(states_elec_t), intent(in   ) :: st
 
    push_sub(ace_init)
 
    !%Variable ACESize
    !%Type integer
    !%Default All states
    !%Section Hamiltonian
    !%Description
    !% (Experimental) The number of ACE projection vectors (i.e. the size of the test-orbital set) and hence the dimension of the
    !% subspace on which the low-rank ACE operator is (approximately) exact. By default, Octopus sets this to the
    !% number of states requested for a calculation (all states), which is essential when the band gap is of interest.
    !% For development purposes, if only occupied states are required, the user can set this value manually.
    !% For more information, see Lin, J. Chem. Theory Comput. 2016, 12, 2242.
    !%End
    call parse_variable(namespace, 'ACESize', st%nst, this%nst)
 
    if (this%nst > st%nst) then
      call messages_input_error(namespace, 'ACESize', "Exceeds the total number of states available")
    endif
 
    if (this%nst < st%nst) then
      call messages_experimental("Adaptively-compressed exchange defined with a subset of states", namespace=namespace)
    endif
 
    if (this%nst * st%smear%el_per_state < st%qtot) then
      write(message(1),'(a)') "ACESize should at least equal the number of occupied states."
      call messages_warning(1, namespace=namespace)
    endif
 
    ! dchi/zchi are allocated at their point of use because at this point there is
    ! nothing to indicate whether one needs dchi or zchi
 
    pop_sub(ace_init)
 
  end subroutine ace_init
 
  subroutine ace_end(this)
    class(ACE_t), intent(inout) :: this
 
    integer :: ik
 
    push_sub(ace_end)
 
    this%nst = 0
    safe_deallocate_a(this%dchi)
    safe_deallocate_a(this%zchi)
 
    if (allocated(this%chib)) then
      do ik = lbound(this%chib, 1), ubound(this%chib, 1)
        call this%chib(ik)%end()
      end do
      safe_deallocate_a(this%chib)
    end if
 
    pop_sub(ace_end)
 
  end subroutine ace_end
 
  subroutine ace_write_info(this, namespace)
    class(ACE_t),       intent(in) :: this
    type(namespace_t),  intent(in) :: namespace
 
    push_sub(ace_write_info)
 
    call messages_print_var_value("Dimension of the subspace in which the low-rank ACE operator is"// &
    & " approximately exact (ACESize)", this%nst, namespace=namespace)
 
    pop_sub(ace_write_info)
 
  end subroutine ace_write_info
 
  subroutine exchange_operator_init(this, namespace, space, st, der, mc, stencil, kpoints, cam)
    type(exchange_operator_t), intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    type(states_elec_t),       intent(in)    :: st
    type(derivatives_t),       intent(in)    :: der
    type(multicomm_t),         intent(in)    :: mc
    type(stencil_t),           intent(in)    :: stencil
    type(kpoints_t),           intent(in)    :: kpoints
    type(xc_cam_t),            intent(in)    :: cam
 
    push_sub(exchange_operator_init)
 
    !%Variable AdaptivelyCompressedExchange
    !%Type logical
    !%Default true
    !%Section Hamiltonian
    !%Description
    !% If set to yes, Octopus will use the adaptively compressed exchange
    !% operator (ACE) for HF and hybrid calculations, as defined in
    !%  Lin, J. Chem. Theory Comput. 2016, 12, 2242.
    !%
    !% This is currently ignored for TheoryLevel = rdmft
    !%End
    call parse_variable(namespace, 'AdaptivelyCompressedExchange', .true., this%useACE)
    call messages_print_var_value('AdaptivelyCompressedExchange', this%useACE)
 
    !%Variable ExchangeSinglePrecision
    !%Type logical
    !%Default no
    !%Section Hamiltonian
    !%Description
    !% If set to yes, the computation of the exchange potentials in
    !% <tt>exchange_operator_compute_potentials</tt> is performed in single
    !% precision: the local states are copied to a single-precision shadow and
    !% the codensities, the Poisson solves and the accumulation are all done in
    !% single precision, before the result is converted back to double precision.
    !% This trades some accuracy for speed and reduced memory and communication
    !% traffic. State-, k-point- and domain-parallel runs are all supported.
    !%
    !% Since the whole Poisson solver of the exchange operator is then set up in single
    !% precision, this requires <tt>AdaptivelyCompressedExchange = yes</tt>, <tt>ACEWithISDF = no</tt>
    !% and a batch-capable FFT Poisson solver (the CPU FFTW solver or cuFFT/rocFFT for GPUs; the
    !% multipole-corrected kernel and the non-FFT solvers have no single-precision variant).
    !% It is incompatible with RDMFT, with the BerkeleyGW output and with the two-body
    !% matrix-element output, which need the exchange operator in double precision.
    !%
    !% Both cube types are supported. The FFT cube is sized once, when the Hamiltonian is
    !% built, so a run whose ground state is real (a finite or Gamma-only system) keeps a
    !% real cube even after a time propagation has made the states complex. The complex
    !% codensities are then transformed as their real and imaginary parts separately, in
    !% single precision just as in double.
    !%End
    call parse_variable(namespace, 'ExchangeSinglePrecision', .false., this%compute_single_precision)
    call messages_print_var_value('ExchangeSinglePrecision', this%compute_single_precision)
 
    !%Variable ACEWithISDF
    !%Type logical
    !%Default no
    !%Section ISDF
    !%Description
    !% If set to yes, Octopus will use interpolative separable density fitting (ISDF)
    !% to accelerate the calculation of adaptively compressed exchange in hybrid
    !% functionals. For more details, please refer to J.Chem.TheoryComput.2017, 13, 5420-5431.
    !% ISDF is currently only implemented for spin-unpolarized, molecular systems.
    !%End
    call parse_variable(namespace, 'ACEWithISDF', .false., this%with_isdf)
 
    if (this%with_isdf .and. .not. this%useACE) then
      call messages_input_error(namespace, 'ACEWithISDF', &
        '"ACEWithISDF = yes" must be used with "AdaptivelyCompressedExchange = yes"')
    endif
 
    ! Objs initialised by exchange constructor
    if (this%useACE) then
      call this%ace%init(namespace, st)
      if (this%with_isdf) call this%isdf%init(namespace, space, der%mesh, this%ace%nst)
    endif
    call singularity_init(this%singul, namespace, space, st, kpoints)
    ! exchange_operator_compute_potentials solves the codensities of a whole state block in one
    ! batched call, so size the cube to hold one block. The active count per call (psib%nst) may be
    ! smaller than the capacity (the last/partial block), and solvers that cannot batch fall back to
    ! a per-function loop; both are handled inside X(poisson_solve_batch).
    call poisson_init(this%psolver, namespace, space, der, mc, stencil, st%qtot, &
      force_serial = .true., verbose = .false., force_cmplx = .not. states_are_real(st), &
      fft_batch_size = st%block_size, single_precision = this%compute_single_precision)
 
    ! Single precision requires a batch-capable solver: there is no single-precision
    ! variant of the non-FFT solvers nor of the multipole-corrected FFT kernel, so the
    ! per-function fallback of X(poisson_solve_batch) has no single-precision solver to
    ! call. poisson_is_batch_capable covers both the host and the accelerated FFT path
    ! (poisson_is_device_batch_capable implies it).
    ! ACE is required as well: the solver above is now single precision, and without ACE
    ! the Hamiltonian application (exchange_operator_apply_standard) would use that very
    ! same solver in double precision. ISDF also solves with it in double precision.
    if (this%compute_single_precision .and. &
      .not. (this%useACE .and. .not. this%with_isdf .and. poisson_is_batch_capable(this%psolver))) then
      call messages_not_implemented("ExchangeSinglePrecision without a batched FFT exchange path " // &
        "(needs AdaptivelyCompressedExchange = yes, ACEWithISDF = no and an FFT Poisson solver " // &
        "with a batch-capable kernel)", namespace=namespace)
    end if
 
    ! Objs initialised by the caller
    this%cam = cam
 
    pop_sub(exchange_operator_init)
  end subroutine exchange_operator_init
 
  subroutine exchange_operator_reinit(this, cam, st)
    type(exchange_operator_t),             intent(inout) :: this
    type(xc_cam_t),                        intent(in   ) :: cam
    type(states_elec_t), target, optional, intent(in   ) :: st
 
    push_sub(exchange_operator_reinit)
 
    if (present(st)) then
      this%st => st
    end if
 
    this%cam = cam
 
    pop_sub(exchange_operator_reinit)
  end subroutine exchange_operator_reinit
 
  subroutine exchange_operator_end(this)
    type(exchange_operator_t), intent(inout) :: this
 
    push_sub(exchange_operator_end)
 
    if (associated(this%st) .and. .not. this%useACE) then
      call states_elec_parallel_remote_access_stop(this%st)
      call states_elec_end(this%st)
      safe_deallocate_p(this%st)
    end if
    nullify(this%st)
 
    if (single_shadows_init) then
      call states_elec_end(st_single_shadow)
      call states_elec_end(xst_single_shadow)
      single_shadows_init = .false.
      single_shadows_np = 0
    end if
 
    call this%ace%end()
    call singularity_end(this%singul)
    call coulb%end()
    call poisson_end(this%psolver)
    call this%isdf%end()
 
    pop_sub(exchange_operator_end)
  end subroutine exchange_operator_end
 
  subroutine exchange_operator_rdmft_occ_apply(this, mesh, hpsib)
    type(exchange_operator_t), intent(in) :: this
    type(mesh_t),              intent(in) :: mesh
    class(wfs_elec_t),      intent(inout) :: hpsib
 
    push_sub(exchange_operator_rdmft_occ_apply)
 
    ! multiply linear terms in hamiltonian with occupation number
    ! nonlinear occupation number dependency occurs only in the exchange, which is treated there
    call batch_scal(mesh%np, this%st%occ(:, hpsib%ik), hpsib)
 
    pop_sub(exchange_operator_rdmft_occ_apply)
  end subroutine exchange_operator_rdmft_occ_apply
 
  subroutine exchange_operator_write_info(this, namespace)
    class(exchange_operator_t), intent(in) :: this
    type(namespace_t),          intent(in) :: namespace
 
    push_sub(exchange_operator_write_info)
 
    call messages_print_with_emphasis(msg='Exact Exchange', namespace=namespace)
    call messages_print_var_value("Adaptively compressed exchange", this%useACE, namespace=namespace)
    if (this%useACE) then
      call this%ace%write_info(namespace)
      if (this%with_isdf) then
        call messages_print_var_value("Density Fitting in ACE with ISDF", this%with_isdf, namespace=namespace)
        call this%isdf%write_info(namespace)
      endif
    endif
 
    pop_sub(exchange_operator_write_info)
 
  end subroutine exchange_operator_write_info
 
 
#include "undef.F90"
#include "real.F90"
#include "exchange_operator_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "exchange_operator_inc.F90"
 
! The exchange-potential routines are shared by all precisions (one source).
#include "undef.F90"
#include "real.F90"
#include "exchange_operator_potentials_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "exchange_operator_potentials_inc.F90"
 
#include "undef.F90"
#include "real32.F90"
#include "exchange_operator_potentials_inc.F90"
 
#include "undef.F90"
#include "complex32.F90"
#include "exchange_operator_potentials_inc.F90"
 
end module exchange_operator_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: