invert_ks.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 invert_ks_oct_m
  use comm_oct_m
  use debug_oct_m
  use density_oct_m
  use eigensolver_oct_m
  use electron_space_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use lalg_basic_oct_m
  use output_oct_m
  use output_modelmb_oct_m
  use io_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use multisystem_basic_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pcm_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use restart_oct_m
  use states_elec_restart_oct_m
  use electrons_oct_m
  use xc_ks_inversion_oct_m
 
  implicit none
 
  private
  public :: invert_ks_run
 
contains
 
  ! ---------------------------------------------------------
  subroutine invert_ks_run(system)
    class(*), intent(inout) :: system
 
    push_sub(invert_ks_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = invert_ks not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call invert_ks_run_legacy(system)
    end select
 
    pop_sub(invert_ks_run)
  end subroutine invert_ks_run
 
  ! ---------------------------------------------------------
  subroutine invert_ks_run_legacy(sys)
    type(electrons_t), intent(inout) :: sys
 
    integer :: ii, jj, np, nspin
    integer :: err
    real(real64)   :: diffdensity
    real(real64), allocatable :: target_rho(:,:), rho(:)
    type(restart_t) :: restart
 
    push_sub(invert_ks_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%kpoints%use_symmetries) then
      call messages_experimental("KPoints symmetries with CalculationMode = invert_ks", namespace=sys%namespace)
    end if
 
    if (sys%st%d%ispin == spinors) then
      call messages_not_implemented("Kohn-Sham inversion with spinors", namespace=sys%namespace)
    end if
 
    ! initialize KS inversion module
    call xc_ks_inversion_write_info(sys%ks%ks_inversion, namespace=sys%namespace)
 
    !abbreviations
    np      = sys%gr%np
    nspin   = sys%st%d%nspin
 
    ! read target density
    safe_allocate(target_rho(1:np, 1:nspin))
    safe_allocate(rho(1:np))
 
    call read_target_rho()
 
    sys%hm%energy%intnvxc = m_zero
    sys%hm%energy%hartree = m_zero
    sys%hm%energy%exchange = m_zero
    sys%hm%energy%correlation = m_zero
    sys%hm%vxc = m_zero
 
    ! calculate total density
    call lalg_copy(np, target_rho(:, 1), rho)
    do ii = 2, sys%st%d%spin_channels
      call lalg_axpy(np, m_one, target_rho(:, ii), rho)
    end do
 
    ! calculate the Hartree potential
    call dpoisson_solve(sys%hm%psolver, sys%namespace, sys%hm%vhartree, rho)
 
    ! Copy it to vhxc
    do ii = 1, sys%st%d%spin_channels
      call lalg_copy(np, sys%hm%vhartree, sys%hm%vhxc(:, ii))
    end do
 
 
    write(message(1),'(a)') "Calculating KS potential"
    call messages_info(1, namespace=sys%namespace)
 
    ! Update the Hamiltonian
    call invertks_update_hamiltonian(sys%namespace, sys%gr, sys%space, sys%ext_partners, &
      sys%ks%ks_inversion%eigensolver, sys%ks%ks_inversion%aux_st, sys%hm, sys%ks%ks_inversion%max_iter)
 
    ! Perform the KS Inversion
    if (sys%ks%ks_inversion%method == xc_inv_method_two_particle) then ! 2-particle exact inversion
 
      call invertks_2part(sys%ks%ks_inversion, target_rho, nspin, sys%hm, sys%gr, &
        sys%ks%ks_inversion%aux_st, sys%ks%ks_inversion%eigensolver, sys%namespace, sys%space, &
        sys%ext_partners)
 
    else ! iterative case
      if (sys%ks%ks_inversion%method >= xc_inv_method_vs_iter .and. &
        sys%ks%ks_inversion%method <= xc_inv_method_iter_godby) then ! iterative procedure for v_s
        call invertks_iter(sys%ks%ks_inversion, target_rho, sys%namespace, sys%space, sys%ext_partners, nspin, sys%hm, sys%gr, &
          sys%ks%ks_inversion%aux_st, sys%ks%ks_inversion%eigensolver)
      end if
    end if
 
    ! Update the Hamiltonian again
    call invertks_update_hamiltonian(sys%namespace, sys%gr, sys%space, sys%ext_partners, &
      sys%ks%ks_inversion%eigensolver, sys%ks%ks_inversion%aux_st, sys%hm, sys%ks%ks_inversion%max_iter)
 
    ! output quality of KS inversion
    diffdensity = m_zero
    do jj = 1, nspin
      do ii = 1, np
        if (abs(sys%ks%ks_inversion%aux_st%rho(ii,jj)-target_rho(ii,jj)) > diffdensity) then
          diffdensity = abs(sys%ks%ks_inversion%aux_st%rho(ii,jj)-target_rho(ii,jj))
        end if
      end do
    end do
    write (message(1),'(a,F16.6)') 'Achieved difference in densities wrt target:', &
      diffdensity
    call messages_info(1, namespace=sys%namespace)
 
    ! output for all cases
    call output_all(sys%outp, sys%namespace, sys%space, static_dir, sys%gr, sys%ions, -1, &
      sys%ks%ks_inversion%aux_st, sys%hm, sys%ks)
    call output_modelmb(sys%outp, sys%namespace, sys%space, static_dir, sys%gr, sys%ions, -1, &
      sys%ks%ks_inversion%aux_st)
 
    sys%ks%ks_inversion%aux_st%dom_st_kpt_mpi_grp = sys%st%dom_st_kpt_mpi_grp
    ! save files in restart format
    call restart_init(restart, sys%namespace, restart_gs, restart_type_dump, sys%mc, err, mesh = sys%gr)
    call states_elec_dump(restart, sys%space, sys%ks%ks_inversion%aux_st, sys%gr, sys%kpoints, err, 0)
    if (err /= 0) then
      message(1) = "Unable to write states wavefunctions."
      call messages_warning(1, namespace=sys%namespace)
    end if
    call restart_end(restart)
 
    safe_deallocate_a(target_rho)
    safe_deallocate_a(rho)
 
    pop_sub(invert_ks_run_legacy)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine read_target_rho()
      character(len=256) :: filename
      integer :: pass, iunit, ierr, ii, npoints
      real(real64)   :: l_xx(sys%space%dim), l_ff(nspin), rr
      real(real64), allocatable :: xx(:,:), ff(:,:)
      character(len=1)   :: char
 
      push_sub(invert_ks_run_legacy.read_target_rho)
 
      ! FIXME: just use restart/gs/density*.obf file rather than needing to set this and use a different format.
 
      !%Variable InvertKSTargetDensity
      !%Type string
      !%Default <tt>target_density.dat</tt>
      !%Section Calculation Modes::Invert KS
      !%Description
      !% Name of the file that contains the density used as the target in the
      !% inversion of the KS equations.
      !%End
      call parse_variable(sys%namespace, 'InvertKSTargetDensity', "target_density.dat", filename)
 
      iunit = io_open(filename, sys%namespace, action='read', status='old')
 
      npoints = 0
      do pass = 1, 2
        ii = -1
        rewind(iunit)
        do
          if (ii == -1) then
            read(iunit, '(a)', advance='no') char
            ii = 0
            if (char == '#') read(iunit, '(a)') char
          end if
          read(iunit, fmt=*, iostat=ierr) l_xx(:), l_ff(:)
          if (ierr /= 0) exit
          ii = ii + 1
          if (pass == 1) npoints = npoints + 1
          if (pass == 2) then
            xx(ii, :) = l_xx(:)
            ff(ii, :) = l_ff(:)
          end if
        end do
        if (pass == 1) then
          safe_allocate(xx(1:npoints, 1:sys%space%dim))
          safe_allocate(ff(1:npoints, 1:nspin))
        end if
      end do
 
      do ii = 1, nspin
        call dmf_interpolate_points(sys%space%dim, npoints, xx, ff(:,ii), &
          np, sys%gr%x, target_rho(:, ii))
      end do
 
      ! we now renormalize the density (necessary if we have a charged system)
      rr = m_zero
      do ii = 1, sys%st%d%spin_channels
        rr = rr + dmf_integrate(sys%gr, target_rho(:, ii), reduce = .false.)
      end do
      if (sys%gr%parallel_in_domains) then
        call sys%gr%allreduce(rr)
      end if
      rr = sys%st%qtot/rr
      target_rho(:,:) = rr*target_rho(:,:)
 
      safe_deallocate_a(xx)
      safe_deallocate_a(ff)
 
      pop_sub(invert_ks_run_legacy.read_target_rho)
    end subroutine read_target_rho
 
  end subroutine invert_ks_run_legacy
 
 
end module invert_ks_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: