vdw.F90

Go to the documentation of this file.

!! Copyright (C) 2006 Hyllios
!!
!! 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 vdw_oct_m
  use debug_oct_m
  use em_resp_calc_oct_m
  use gauss_legendre_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use linear_response_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_basic_oct_m
  use parser_oct_m
  use perturbation_electric_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use sternheimer_oct_m
  use electrons_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
 
  implicit none
 
  private
  public :: &
    vdw_run
 
contains
 
  ! ---------------------------------------------------------
  subroutine vdw_run(system, from_scratch)
    class(*),        intent(inout) :: system
    logical,         intent(in)    :: from_scratch
 
    push_sub(vdw_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = vdw not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call vdw_run_legacy(system, from_scratch)
    end select
 
    pop_sub(vdw_run)
  end subroutine vdw_run
 
  ! ---------------------------------------------------------
  subroutine vdw_run_legacy(sys, fromScratch)
    type(electrons_t),    intent(inout) :: sys
    logical,              intent(in)    :: fromScratch
 
    type(lr_t) :: lr(sys%space%dim, 1)
    type(sternheimer_t)     :: sh
 
    integer :: dir, i, iunit, gauss_start, ndir
    complex(real64) :: omega
    real(real64) :: domega, pol, c3, c6, cat
 
    integer :: gaus_leg_n
    real(real64), allocatable :: gaus_leg_points(:), gaus_leg_weights(:)
    real(real64), parameter :: omega0 = 0.3_real64
 
    type(restart_t) :: restart_dump
 
    push_sub(vdw_run_legacy)
 
    if (sys%hm%pcm%run_pcm) then
      call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
    end if
 
    if (sys%space%is_periodic()) then
      call messages_not_implemented('Van der Waals calculation for periodic system', namespace=sys%namespace)
    end if
 
    if (sys%kpoints%use_symmetries) call messages_experimental("KPoints symmetries with CalculationMode = vdw", &
      namespace=sys%namespace)
 
    call input()
    call init_()
    call sternheimer_init(sh, sys%namespace, sys%space, sys%gr, sys%st, sys%hm, sys%ks%xc, sys%mc, wfs_are_cplx = .true.)
 
    if (gauss_start == 1 .and. mpi_grp_is_root(mpi_world)) then
      iunit = io_open(vdw_dir//'vdw_c6', sys%namespace, action='write')
      write(iunit, '(a,i3)') '# npoints = ', gaus_leg_n
      write(iunit, '(a1,a12,2a20)') '#', 'omega', 'domega', 'pol'
      call io_close(iunit)
    end if
 
    do i = gauss_start, gaus_leg_n
      omega  = m_zi*omega0*(m_one - gaus_leg_points(i))/(m_one + gaus_leg_points(i))
      domega = gaus_leg_weights(i) * omega0 * (m_two)/(m_one + gaus_leg_points(i))**2
 
      pol = get_pol(omega)
      if (mpi_grp_is_root(mpi_world)) then
        iunit = io_open(vdw_dir//'vdw_c6', sys%namespace, action='write', position='append')
        write(iunit, '(3es20.12)') aimag(omega), domega, pol
        call io_close(iunit)
      end if
 
      c3  = c3  + m_three/m_pi * domega * pol
      c6  = c6  + m_three/m_pi * domega * pol**2
      cat = cat + m_three/m_pi * domega * pol**3
    end do
 
    if ((gauss_start  <=  gaus_leg_n).and.mpi_grp_is_root(mpi_world)) then
      iunit = io_open(vdw_dir//'vdw_c6', sys%namespace, action='write', position='append')
      write(iunit, '(1x)')
      write(iunit, '(a,es20.12)') "C_3  [a.u. ] = ", c3
      write(iunit, '(a,es20.12)') "C_6  [a.u. ] = ", c6
      write(iunit, '(a,es20.12)') "C_AT [a.u. ] = ", cat
      write(iunit, '(1x)')
 
      write(iunit, '(3a,es20.12)') "C_3  [", &
        trim(units_abbrev(units_out%energy * units_out%length**sys%space%dim)), "] = ", &
        units_from_atomic(units_out%energy * units_out%length**sys%space%dim, c3)
      write(iunit, '(3a,es20.12)') "C_6  [", &
        trim(units_abbrev(units_out%energy * units_out%length**(2*sys%space%dim))), "] = ", &
        units_from_atomic(units_out%energy * units_out%length**(2*sys%space%dim), c6)
      write(iunit, '(3a,es20.12)') "C_AT [", &
        trim(units_abbrev(units_out%energy * units_out%length**(3*sys%space%dim))), "] = ", &
        units_from_atomic(units_out%energy * units_out%length**(3*sys%space%dim), cat)
 
      call io_close(iunit)
    end if
 
    call sternheimer_end(sh)
    call end_()
 
    pop_sub(vdw_run_legacy)
  contains
 
    ! --------------------------------------------------------------------
    subroutine input()
      integer :: equiv_axes
 
      push_sub(vdw_run_legacy.input)
 
      !%Variable vdWNPoints
      !%Type integer
      !%Default 6
      !%Section Linear Response::Polarizabilities
      !%Description
      !% How many points to use in the Gauss-Legendre integration to obtain the
      !% van der Waals coefficients.
      !%End
      call messages_obsolete_variable(sys%namespace, 'vdW_npoints', 'vdWNPoints')
      call parse_variable(sys%namespace, 'vdWNPoints', 6, gaus_leg_n)
 
      ! \todo symmetry stuff should be general
      call parse_variable(sys%namespace, 'TDPolarizationEquivAxes', 0, equiv_axes)
 
      select case (equiv_axes)
      case (3)
        ndir = 1
      case (2)
        ndir = min(2, sys%space%dim)
      case default
        ndir = min(3, sys%space%dim)
      end select
 
      pop_sub(vdw_run_legacy.input)
    end subroutine input
 
 
    ! --------------------------------------------------------------------
    subroutine init_()
      integer :: ierr, iunit, ii
      logical :: file_exists
      character(len=80) :: dirname
      real(real64) :: iomega, domega, pol
 
      type(restart_t) :: restart_load, gs_restart
 
      push_sub(vdw_run_legacy.init_)
 
      ! make some space for static polarizability
      gaus_leg_n = gaus_leg_n + 1
 
      ! get Gauss-Legendre points
      safe_allocate(gaus_leg_points(1:gaus_leg_n))
      safe_allocate(gaus_leg_weights(1:gaus_leg_n))
 
      call gauss_legendre_points(gaus_leg_n-1, gaus_leg_points, gaus_leg_weights)
      c3  = m_zero
      c6  = m_zero
      cat = m_zero
      gauss_start = 1
      gaus_leg_points(gaus_leg_n) = 0.99999_real64
      gaus_leg_weights(gaus_leg_n) = m_zero
 
      ! FIXME: this should be part of the restart framework
      ! check if we can restart
      inquire(file=vdw_dir//'vdw_c6', exist=file_exists)
      if (.not. fromscratch .and. file_exists) then
        iunit = io_open(vdw_dir//'vdw_c6', sys%namespace, action='read')
        read(iunit, '(a12,i3)', iostat=ierr) dirname, ii
        if (ii /= gaus_leg_n) then
          message(1) = "Invalid restart of van der Waals calculation."
          message(2) = "The number of points in the Gauss-Legendre integration changed."
          write(message(3), '(i3,a,i3,a)') gaus_leg_n, " (input) != ", ii, "(restart)"
          call messages_fatal(3, namespace=sys%namespace)
        end if
        read(iunit,*) ! skip comment line
        do
          read(iunit, *, iostat=ierr) iomega, domega, pol
          if (ierr /= 0) exit
          gauss_start = gauss_start + 1
          c3  = c3  + m_three/m_pi * domega * pol
          c6  = c6  + m_three/m_pi * domega * pol**2
          cat = cat + m_three/m_pi * domega * pol**3
        end do
        call io_close(iunit)
      end if
 
      ! we always need complex response
      call restart_init(gs_restart, sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr, exact=.true.)
      if (ierr == 0) then
        call states_elec_look_and_load(gs_restart, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
          is_complex = .true.)
        call restart_end(gs_restart)
      else
        message(1) = "Previous gs calculation required."
        call messages_fatal(1, namespace=sys%namespace)
      end if
 
      ! setup Hamiltonian
      message(1) = 'Info: Setting up Hamiltonian for linear response.'
      call messages_info(1, namespace=sys%namespace)
      call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
 
      do dir = 1, ndir
        call lr_init(lr(dir,1))
        call lr_allocate(lr(dir,1), sys%st, sys%gr)
      end do
 
      ! load wavefunctions
      if (.not. fromscratch) then
        call restart_init(restart_load, sys%namespace, restart_vdw, restart_type_load, sys%mc, ierr, mesh=sys%gr)
 
        do dir = 1, ndir
          write(dirname,'(a,i1,a)') "wfs_", dir, "_1_1"
          call restart_open_dir(restart_load, dirname, ierr)
          if (ierr == 0) then
            call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, ierr, lr=lr(dir,1))
          end if
          if (ierr /= 0) then
            message(1) = "Unable to read response wavefunctions from '"//trim(dirname)//"'."
            call messages_warning(1, namespace=sys%namespace)
          end if
          call restart_close_dir(restart_load)
        end do
 
        call restart_end(restart_load)
      end if
 
      if (mpi_grp_is_root(mpi_world)) then
        call io_mkdir(vdw_dir, sys%namespace)               ! output data
      end if
 
      call restart_init(restart_dump, sys%namespace, restart_vdw, restart_type_dump, sys%mc, ierr, mesh=sys%gr)
 
      pop_sub(vdw_run_legacy.init_)
    end subroutine init_
 
    ! --------------------------------------------------------------------
    subroutine end_()
      integer :: dir
 
      push_sub(vdw_run_legacy.end_)
 
      safe_deallocate_a(gaus_leg_points)
      safe_deallocate_a(gaus_leg_weights)
 
      do dir = 1, ndir
        call lr_dealloc(lr(dir, 1))
      end do
 
      call restart_end(restart_dump)
 
      pop_sub(vdw_run_legacy.end_)
    end subroutine end_
 
 
    ! --------------------------------------------------------------------
    real(real64) function get_pol(omega)
      complex(real64), intent(in) :: omega
 
      complex(real64)   :: alpha(1:sys%space%dim, 1:sys%space%dim)
      type(perturbation_electric_t), pointer :: pert
 
      push_sub(vdw_run_legacy.get_pol)
 
      pert => perturbation_electric_t(sys%namespace)
      do dir = 1, ndir
        write(message(1), '(3a,f7.3)') 'Info: Calculating response for the ', index2axis(dir), &
          '-direction and imaginary frequency ', units_from_atomic(units_out%energy, aimag(omega))
        call messages_info(1, namespace=sys%namespace)
 
        call pert%setup_dir(dir)
        call zsternheimer_solve(sh, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, sys%mc, &
          lr(dir, :), 1, omega, pert, restart_dump, em_rho_tag(real(omega, real64) ,dir), em_wfs_tag(dir,1))
      end do
 
      call zcalc_polarizability_finite(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, lr(:,:), 1, pert, &
        alpha(:,:), ndir = ndir)
 
      get_pol = m_zero
      do dir = 1, ndir
        get_pol = get_pol + real(alpha(dir, dir), real64)
      end do
      do dir = ndir+1, sys%space%dim
        get_pol = get_pol + real(alpha(ndir, ndir), real64)
      end do
 
      get_pol = get_pol / real(sys%space%dim, real64)
 
      safe_deallocate_p(pert)
      pop_sub(vdw_run_legacy.get_pol)
    end function get_pol
 
  end subroutine vdw_run_legacy
 
end module vdw_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: