libvdwxc.F90

Go to the documentation of this file.

#include "global.h"
 
! Interface to libvdwxc.
! Naming:
!  * Functions that start with libvdwxc_ are public, to be called from other parts of Octopus.
!  * Interfaces that start with vdwxc_ are actual functions of libvdwxc.
 
 
module libvdwxc_oct_m
  use box_parallelepiped_oct_m
  use cube_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use fft_oct_m
  use global_oct_m
  use lalg_basic_oct_m
  use grid_oct_m
  use io_function_oct_m
  use mesh_oct_m
  use mesh_cube_parallel_map_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pfft_oct_m
  use profiling_oct_m
  use space_oct_m
  use unit_system_oct_m
 
  implicit none
 
  integer, parameter ::       &
    LIBVDWXC_MODE_AUTO = 1,   &
    libvdwxc_mode_serial = 2, &
    libvdwxc_mode_mpi = 3!, &
  !LIBVDWXC_MODE_PFFT = 3
 
  private
 
  public ::                &
    libvdwxc_t,            &
    libvdwxc_init,         &
    libvdwxc_print,        &
    libvdwxc_write_info,   &
    libvdwxc_set_geometry, &
    libvdwxc_calculate,    &
    libvdwxc_end
 
  type libvdwxc_t
    private
    integer, pointer               :: libvdwxc_ptr
    type(mesh_t)                   :: mesh
    type(cube_t)                   :: cube
    type(mesh_cube_parallel_map_t) :: mesh_cube_map
    integer                        :: functional
    logical                        :: debug
    real(real64), public                  :: energy
    real(real64)                   :: vdw_factor
  end type libvdwxc_t
 
#ifdef HAVE_LIBVDWXC
  include "vdwxcfort.f90"
#endif
 
contains
 
  subroutine libvdwxc_init(libvdwxc, namespace, functional)
    type(libvdwxc_t),          intent(out) :: libvdwxc
    type(namespace_t),         intent(in)  :: namespace
    integer,                   intent(in)  :: functional
 
    push_sub(libvdwxc_init)
#ifdef HAVE_LIBVDWXC
    call vdwxc_new(functional, libvdwxc%libvdwxc_ptr)
#else
    message(1) = "Octopus not compiled with libvdwxc"
    call messages_fatal(1, namespace=namespace)
#endif
    assert(associated(libvdwxc%libvdwxc_ptr))
    libvdwxc%functional = functional
 
    !%Variable libvdwxcDebug
    !%Type logical
    !%Section Hamiltonian::XC
    !%Description
    !% Dump libvdwxc inputs and outputs to files.
    !%End
    call parse_variable(namespace, 'libvdwxcDebug', .false., libvdwxc%debug)
    pop_sub(libvdwxc_init)
 
    !%Variable libvdwxcVDWFactor
    !%Type float
    !%Section Hamiltonian::XC
    !%Description
    !% Prefactor of non-local van der Waals functional.
    !% Setting a prefactor other than one is wrong, but useful
    !% for debugging.
    !%End
    call parse_variable(namespace, 'libvdwxcVDWFactor', m_one, libvdwxc%vdw_factor)
  end subroutine libvdwxc_init
 
  subroutine libvdwxc_print(this)
    type(libvdwxc_t), intent(in) :: this
    push_sub(libvdwxc_print)
 
#ifdef HAVE_LIBVDWXC
    call vdwxc_print(this%libvdwxc_ptr)
#endif
 
    pop_sub(libvdwxc_print)
  end subroutine libvdwxc_print
 
  subroutine libvdwxc_write_info(this, iunit, namespace)
    type(libvdwxc_t),            intent(in) :: this
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    push_sub(libvdwxc_write_info)
 
    write(message(1), '(2x,a)') 'Correlation'
    if (this%functional == 1) then
      write(message(2), '(4x,a)') 'vdW-DF from libvdwxc'
    else if (this%functional == 2) then
      write(message(2), '(4x,a)') 'vdW-DF2 from libvdwxc'
    else if (this%functional == 3) then
      write(message(2), '(4x,a)') 'vdW-DF-cx from libvdwxc'
    else
      write(message(2), '(4x,a)') 'unknown libvdwxc functional'
    end if
    call messages_info(2, iunit=iunit, namespace=namespace)
 
    pop_sub(libvdwxc_write_info)
  end subroutine libvdwxc_write_info
 
  subroutine libvdwxc_set_geometry(this, namespace, space, mesh)
    type(libvdwxc_t),  intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
    class(space_t),    intent(in)    :: space
    class(mesh_t),     intent(in)    :: mesh
 
    integer :: blocksize
    integer :: libvdwxc_mode
    type(space_t) :: cube_space
 
    push_sub(libvdwxc_set_geometry)
    this%mesh = mesh
 
    ! libvdwxc can use either FFTW-MPI or PFFT and Octopus should not
    ! care too much, as long as this particular cube is properly
    ! configured.  Unfortunately most of the args to cube_init are
    ! very entangled with FFT libraries.  All we want and need is to
    ! specify our own decomposition and pass that to libvdwxc, but
    ! we can only say "we want to use PFFT" and then it does the
    ! decomposition.
    libvdwxc_mode = libvdwxc_mode_auto
 
    !%Variable libvdwxcMode
    !%Type integer
    !%Section Hamiltonian::XC
    !%Description
    !% Whether libvdwxc should run with serial fftw3, fftw3-mpi, or pfft.
    !% to specify fftw3-mpi in serial for debugging.
    !% pfft is not implemented at the moment.
    !%Option libvdwxc_mode_auto 1
    !% Use serial fftw3 if actually running in serial, else fftw3-mpi.
    !%Option libvdwxc_mode_serial 2
    !% Run with serial fftw3.  Works only when not parallelizing over domains.
    !%Option libvdwxc_mode_mpi 3
    !% Run with fftw3-mpi.  Works only if Octopus is compiled with MPI.
    !%End
    call parse_variable(namespace, 'libvdwxcMode', libvdwxc_mode_auto, libvdwxc_mode)
 
    if (libvdwxc_mode == libvdwxc_mode_auto) then
      if (mesh%mpi_grp%size == 1) then
        libvdwxc_mode = libvdwxc_mode_serial
      else
        libvdwxc_mode = libvdwxc_mode_mpi
      end if
    end if
 
    ! TODO implement.  Should dump quantities to files.
 
    blocksize = mesh%idx%ll(3) / mesh%mpi_grp%size
    if (mod(mesh%idx%ll(3), mesh%mpi_grp%size) /= 0) then
      blocksize = blocksize + 1
    end if
 
    ! Here we need to assume that the space is periodic as libvdwxc assumes a
    ! periodic space. This affects the conversion from lattice vectors to the spacing
    cube_space%dim = space%dim
    cube_space%periodic_dim = space%dim
 
    if (libvdwxc_mode == libvdwxc_mode_serial) then
      call cube_init(this%cube, mesh%idx%ll, namespace, cube_space, mesh%spacing, &
        mesh%coord_system)
      call cube_init_cube_map(this%cube, mesh)
    else
      call cube_init(this%cube, mesh%idx%ll, namespace, cube_space, mesh%spacing, &
        mesh%coord_system, mpi_grp = mesh%mpi_grp, &
        need_partition = .true., blocksize = blocksize)
      call cube_init_cube_map(this%cube, mesh)
      call mesh_cube_parallel_map_init(this%mesh_cube_map, mesh, this%cube)
    end if
 
    ! There is some low-level implementation issue where with PFFT,
    ! the leading dimension is one smaller than normal for some reason.
    ! Therefore we cannot use the PFFT stuff without a frightful mess.
 
#ifdef HAVE_LIBVDWXC
    call vdwxc_set_unit_cell(this%libvdwxc_ptr, &
      this%cube%rs_n_global(3), this%cube%rs_n_global(2), this%cube%rs_n_global(1), &
      this%cube%latt%rlattice(3, 3), this%cube%latt%rlattice(2, 3), this%cube%latt%rlattice(1, 3), &
      this%cube%latt%rlattice(3, 2), this%cube%latt%rlattice(2, 2), this%cube%latt%rlattice(1, 2), &
      this%cube%latt%rlattice(3, 1), this%cube%latt%rlattice(2, 1), this%cube%latt%rlattice(1, 1))
 
    if (libvdwxc_mode == libvdwxc_mode_serial) then
      call vdwxc_init_serial(this%libvdwxc_ptr)
    else
#ifdef HAVE_LIBVDWXC_MPI
      call vdwxc_init_mpi(this%libvdwxc_ptr, mesh%mpi_grp%comm%MPI_VAL)
#else
      message(1) = "libvdwxc was not compiled with MPI"
      message(2) = "Recompile libvdwxc with MPI for vdW with domain decomposition"
      call messages_fatal(2, namespace=namespace)
#endif
    end if
    call vdwxc_print(this%libvdwxc_ptr)
#endif
 
    pop_sub(libvdwxc_set_geometry)
  end subroutine libvdwxc_set_geometry
 
  subroutine libvdwxc_calculate(this, namespace, space, rho, gradrho, dedd, dedgd)
    type(libvdwxc_t),                     intent(inout) :: this
    type(namespace_t),                    intent(in)    :: namespace
    class(space_t),                       intent(in)    :: space
    real(real64), dimension(:,:),   contiguous,  intent(inout) :: rho !!! data type
    real(real64), dimension(:,:,:), contiguous,  intent(in)    :: gradrho
    real(real64), dimension(:,:),   contiguous,  intent(inout) :: dedd
    real(real64), dimension(:,:,:), contiguous,  intent(inout) :: dedgd
 
    type(cube_function_t) :: cf!rhocf, sigmacf, dedrhocf, dedsigmacf
 
    real(real64), allocatable :: workbuffer(:)
    real(real64), allocatable :: cube_rho(:,:,:), cube_sigma(:,:,:), cube_dedrho(:,:,:), cube_dedsigma(:,:,:)
    real(real64), dimension(3) :: energy_and_integrals_buffer
    real(real64) :: vsum
    integer :: ii, is, idir
 
    push_sub(libvdwxc_calculate)
 
    assert(size(rho, 2) == 1)
    assert(size(gradrho, 3) == 1)
    assert(size(dedd, 2) == 1)
    assert(size(dedgd, 3) == 1)
 
    ! Well.  I thought we would be using four different cube functions
    ! for rho, sigma, dedrho and dedsigma.
    !
    ! But since the cube has PFFT associated, any attempt to use it
    ! (update: We do not actually use PFFT anymore, so maybe this will
    !  not be that broken now)
    ! results in some kind of behind-the-scenes use of a possibly
    ! global FFT-related buffer.
    ! (Note: actually we now use cubes without FFT lib)
    !
    ! cube functions take a force_alloc variable which appears to make
    ! them allocate their own buffer (like we want), but then Octopus
    ! segfaults on cube function free, which I suspect is a bug in the
    ! code, or at the very least due to something so undocumented that
    ! I cannot reasonably figure it out.
    !
    ! We could just disable PFFT for the cube (since we will never
    ! actually call any FFT from Octopus) and do our own
    ! redistribution, but the redistribution code is loaded with
    ! references to FFT library, so this would probably be a bit
    ! optimistic.
    !
    ! So we create here our own arrays over which we have reasonable control.
 
    safe_allocate(workbuffer(1:this%mesh%np))
    safe_allocate(cube_rho(1:this%cube%rs_n(1), 1:this%cube%rs_n(2), 1:this%cube%rs_n(3)))
    safe_allocate(cube_sigma(1:this%cube%rs_n(1), 1:this%cube%rs_n(2), 1:this%cube%rs_n(3)))
    safe_allocate(cube_dedrho(1:this%cube%rs_n(1), 1:this%cube%rs_n(2), 1:this%cube%rs_n(3)))
    safe_allocate(cube_dedsigma(1:this%cube%rs_n(1), 1:this%cube%rs_n(2), 1:this%cube%rs_n(3)))
    ! This is sigma, the absolute-squared density gradient:
    workbuffer(:) = sum(gradrho(:, :, 1)**2, 2)
 
    if (this%debug) then
      call libvdwxc_write_array(rho(:, 1), 'rho')
      call libvdwxc_write_array(workbuffer, 'gradrho')
    end if
 
    cube_rho = m_zero
    cube_sigma = m_zero
    cube_dedrho = m_zero
    cube_dedsigma = m_zero
 
    call dcube_function_alloc_rs(this%cube, cf, in_device = .false.)
 
    call tocube(rho(:, 1), cube_rho)
    call tocube(workbuffer, cube_sigma)
 
    this%energy = m_zero
#ifdef HAVE_LIBVDWXC
    call vdwxc_calculate(this%libvdwxc_ptr, cube_rho, cube_sigma, cube_dedrho, cube_dedsigma, this%energy)
#endif
    this%energy = this%energy * this%vdw_factor
    cube_dedrho = cube_dedrho * this%vdw_factor
    cube_dedsigma = cube_dedsigma * this%vdw_factor
 
    call fromcube(cube_dedrho, workbuffer)
    ! dedd is 1:mesh%np_part for some reason
    if (this%debug) then
      call libvdwxc_write_array(workbuffer, 'dedrho')
    end if
    call lalg_axpy(this%mesh%np, m_one, workbuffer, dedd(:, 1))
    call fromcube(cube_dedsigma, workbuffer)
    if (this%debug) then
      call libvdwxc_write_array(workbuffer, 'dedsigma')
    end if
    do ii = 1, this%mesh%np
      dedgd(ii, :, 1) = dedgd(ii, :, 1) + m_two * workbuffer(ii) * gradrho(ii, :, 1)
    end do
 
    energy_and_integrals_buffer(1) = this%energy
 
    vsum = m_zero
    !$omp parallel do reduction(+:vsum) collapse(2)
    do is = 1, size(rho, 2)
      do ii = 1, this%mesh%np
        vsum = vsum + rho(ii, is) * dedd(ii, is)
      end do
    end do
    !$omp end parallel do
    energy_and_integrals_buffer(2) = vsum * this%mesh%volume_element
 
    vsum = m_zero
    !$omp parallel do reduction(+:vsum) collapse(3)
    do idir = 1, size(gradrho, 3)
      do is = 1, size(gradrho, 2)
        do ii = 1, this%mesh%np
          vsum = vsum + gradrho(ii, is, idir) * dedgd(ii, is, idir)
        end do
      end do
    end do
    !$omp end parallel do
    energy_and_integrals_buffer(3) = vsum * this%mesh%volume_element
 
    call this%mesh%mpi_grp%allreduce_inplace(energy_and_integrals_buffer(1), 3, mpi_double_precision, mpi_sum)
    this%energy = energy_and_integrals_buffer(1)
    write(message(1), '(a,f18.10,a)') 'libvdwxc non-local correlation energy: ', energy_and_integrals_buffer(1), ' Ha'
    write(message(2), '(a,f18.10)')   '                      n-dedn integral: ', energy_and_integrals_buffer(2)
    write(message(3), '(a,f18.10)')   '              gradn-dedgradn integral: ', energy_and_integrals_buffer(3)
    call messages_info(3, namespace=namespace)
 
    safe_deallocate_a(workbuffer)
    safe_deallocate_a(cube_rho)
    safe_deallocate_a(cube_sigma)
    safe_deallocate_a(cube_dedrho)
    safe_deallocate_a(cube_dedsigma)
    call dcube_function_free_rs(this%cube, cf)
 
    pop_sub(libvdwxc_calculate)
 
  contains
 
    subroutine tocube(array, cubearray)
      real(real64), contiguous,     intent(in)    :: array(:)
      real(real64),          intent(out)   :: cubearray(:,:,:)
 
      push_sub(libvdwxc_calculate.tocube)
 
      if (this%cube%parallel_in_domains) then
        call dmesh_to_cube_parallel(this%mesh, array, this%cube, cf, this%mesh_cube_map)
      else
        call dmesh_to_cube(this%mesh, array, this%cube, cf)
      end if
      cubearray(:,:,:) = cf%dRS
      pop_sub(libvdwxc_calculate.tocube)
    end subroutine tocube
 
    subroutine fromcube(cubearray, array)
      real(real64),          intent(in)  :: cubearray(:,:,:)
      real(real64), contiguous,     intent(out) :: array(:)
 
      push_sub(libvdwxc_calculate.fromcube)
      cf%dRS = cubearray
      if (this%cube%parallel_in_domains) then
        call dcube_to_mesh_parallel(this%cube, cf, this%mesh, array, this%mesh_cube_map)
      else
        call dcube_to_mesh(this%cube, cf, this%mesh, array)
      end if
      pop_sub(libvdwxc_calculate.fromcube)
    end subroutine fromcube
 
    subroutine libvdwxc_write_array(arr, fname)
      real(real64),     intent(in) :: arr(:)
      character(len=*), intent(in) :: fname
      integer :: ierr
 
      call dio_function_output(option__outputformat__dx,  'libvdwxc-debug', &
        fname, namespace, space, this%mesh, arr, unit_one, ierr)
    end subroutine libvdwxc_write_array
 
  end subroutine libvdwxc_calculate
 
  subroutine libvdwxc_end(this)
    type(libvdwxc_t), intent(inout) :: this
    push_sub(libvdwxc_end)
 
#ifdef HAVE_LIBVDWXC
    call vdwxc_finalize(this%libvdwxc_ptr)
#endif
 
    call cube_end(this%cube)
    call mesh_cube_parallel_map_end(this%mesh_cube_map)
 
    pop_sub(libvdwxc_end)
  end subroutine libvdwxc_end
 
end module libvdwxc_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: