electrons_ground_state.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 electrons_ground_state_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_function_oct_m
  use ions_oct_m
  use ks_potential_oct_m
  use lcao_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use output_low_oct_m
  use pcm_oct_m
  use rdmft_oct_m
  use restart_oct_m
  use scf_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use v_ks_oct_m
  use mpi_oct_m
 
  implicit none
 
  private
  public ::                                 &
    electrons_ground_state_run,             &
    electrons_gs_allocate_wavefunctions,    &
    electrons_gs_load_from_restart,         &
    electrons_gs_initialize,                &
    electrons_gs_cleanup
 
contains
 
  subroutine electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, fromScratch)
    type(namespace_t),        intent(in)    :: namespace
    type(multicomm_t),        intent(in)    :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(output_t),           intent(in)    :: outp
    type(electron_space_t),   intent(in)    :: space
    logical,                  intent(inout) :: fromScratch
 
    type(rdm_t)     :: rdm
    type(scf_t)     :: scf
 
    push_sub(electrons_ground_state_run)
 
    call electrons_gs_allocate_wavefunctions(namespace, gr, st, hm, scf, ks, ions, mc, space)
 
    if (.not. fromscratch) then
      call electrons_gs_load_from_restart(namespace, scf, gr, mc, st, hm, ks, space, ions, ext_partners, fromscratch)
    end if
 
    call electrons_gs_initialize(namespace, scf, rdm, gr, mc, st, hm, ions, ks, space, ext_partners, fromscratch)
 
    call electrons_gs_run(namespace, scf, rdm, mc, gr, ions, ext_partners, st, ks, hm, outp, space)
 
    call electrons_gs_cleanup(ks, scf, rdm, st, hm)
 
    pop_sub(electrons_ground_state_run)
 
  end subroutine electrons_ground_state_run
 
 
  subroutine electrons_gs_cleanup(ks, scf, rdm, st, hm)
    type(v_ks_t),             intent(inout) :: ks
    type(scf_t),              intent(inout) :: scf
    type(rdm_t),              intent(inout) :: rdm
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    push_sub(electrons_gs_cleanup)
 
    if (ks%theory_level == rdmft) then
      call rdmft_end(rdm)
    else
      call scf_end(scf)
    end if
 
    if (scf%restart_dump%get_data_type() /= restart_undefined) then
      call scf%restart_dump%end()
    end if
 
    if (st%pack_states .and. hm%apply_packed()) then
      call st%unpack()
    end if
 
    ! clean up
    call states_elec_deallocate_wfns(st)
 
    pop_sub(electrons_gs_cleanup)
  end subroutine electrons_gs_cleanup
 
  ! ---------------------------------------------------------
  subroutine electrons_gs_allocate_wavefunctions(namespace, gr, st, hm, scf, ks, ions, mc, space)
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(inout) :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(scf_t),              intent(inout) :: scf
    type(v_ks_t),             intent(inout) :: ks
    type(ions_t),             intent(inout) :: ions
    type(multicomm_t),        intent(in)    :: mc
    type(electron_space_t),   intent(in)    :: space
 
    push_sub(electrons_gs_allocate_wavefunctions)
 
    call messages_write('Info: Allocating ground state wave-functions')
    call messages_info(namespace=namespace)
 
    if (st%parallel_in_states) then
      call messages_experimental('State parallelization for ground state calculations', namespace=namespace)
    end if
 
    if (hm%pcm%run_pcm) then
      if (.not. is_close(hm%pcm%epsilon_infty, hm%pcm%epsilon_0) .and. hm%pcm%tdlevel /= pcm_td_eq) then
        message(1) = 'Non-equilbrium PCM is not active in a time-independent run.'
        message(2) = 'You set epsilon_infty /= epsilon_0, but epsilon_infty is not relevant for CalculationMode = gs.'
        message(3) = 'By definition, the ground state is in equilibrium with the solvent.'
        message(4) = 'Therefore, the only relevant dielectric constant is the static one.'
        message(5) = 'Nevertheless, the dynamical PCM response matrix is evaluated for benchamarking purposes.'
        call messages_warning(5, namespace=namespace)
      end if
    end if
 
    call states_elec_allocate_wfns(st, gr, packed=.true.)
 
    ! sometimes a deadlock can occur here (if some nodes can allocate and other cannot)
    if (st%dom_st_kpt_mpi_grp%comm /= mpi_comm_undefined &
      .and. st%dom_st_kpt_mpi_grp%comm /= mpi_comm_null) then
      call st%dom_st_kpt_mpi_grp%barrier()
    end if
    call messages_write('Info: Ground-state allocation done.')
    call messages_info(namespace=namespace)
 
    pop_sub(electrons_gs_allocate_wavefunctions)
  end subroutine electrons_gs_allocate_wavefunctions
 
  ! ---------------------------------------------------------
  subroutine electrons_gs_load_from_restart(namespace, scf, gr, mc, st, hm, ks, space, ions, ext_partners, fromScratch)
    type(namespace_t),        intent(in)    :: namespace
    type(scf_t),              intent(inout) :: scf
    type(grid_t),             intent(inout) :: gr
    type(multicomm_t),        intent(in)    :: mc
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(v_ks_t),             intent(inout) :: ks
    type(electron_space_t),   intent(in)    :: space
    type(ions_t),             intent(in)    :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    logical,                  intent(inout) :: fromscratch
 
    integer :: ierr
 
    push_sub(electrons_gs_load_from_restart)
 
    ! load wavefunctions
    ! in RDMFT we need the full ground state
    call scf%restart_load%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, &
      exact = (ks%theory_level == rdmft))
    if (ierr == 0) then
      call states_elec_load(scf%restart_load, namespace, space, st, gr, hm%kpoints, ierr)
    end if
 
    if (ierr /= 0) then
      call messages_write("Unable to read wavefunctions.")
      call messages_new_line()
      call messages_write("Starting from scratch!")
      call messages_warning(namespace=namespace)
      fromscratch = .true.
    end if
 
    pop_sub(electrons_gs_load_from_restart)
  end subroutine electrons_gs_load_from_restart
 
  ! ---------------------------------------------------------
  subroutine electrons_gs_initialize(namespace, scf, rdm, gr, mc, st, hm, ions, ks, space, ext_partners, fromScratch)
    type(namespace_t),        intent(in)    :: namespace
    type(scf_t),              intent(inout) :: scf
    type(rdm_t),              intent(inout) :: rdm
    type(grid_t),             intent(inout) :: gr
    type(multicomm_t),        intent(in)    :: mc
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(ions_t),             intent(in)    :: ions
    type(v_ks_t),             intent(inout) :: ks
    type(electron_space_t),   intent(in)    :: space
    type(partner_list_t),     intent(in)    :: ext_partners
    logical,                  intent(in)    :: fromscratch
 
    integer :: ierr
    logical :: restart_init_data_type_dump
 
    push_sub(electrons_gs_initialize)
 
    call write_canonicalized_xyz_file("exec", "initial_coordinates", space, ions%latt, ions%pos, ions%atom, &
      gr%box, namespace)
 
    if (ks%theory_level /= rdmft) then
      call scf_init(scf, namespace, gr, ions, st, mc, hm, space)
      ! only initialize dumping restart files for more than one iteration
      restart_init_data_type_dump = scf%max_iter > 0
    else
      restart_init_data_type_dump = .true.
    end if
 
    if (fromscratch .and. ks%theory_level /= rdmft) then
      call lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, lmm_r = scf%lmm_r)
    else
      ! setup Hamiltonian
      call messages_write('Info: Setting up Hamiltonian.')
      call messages_info(namespace=namespace)
 
      ! v_ks_h_setup calls states_fermi, that uses the eigenvalues, which are uninitialized at this stage
      if (fromscratch .and. ks%theory_level == rdmft) st%eigenval = default_eigenval
 
      call v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, &
        calc_eigenval = .false., calc_current = .false.)
    end if
 
    if (restart_init_data_type_dump) then
      call scf%restart_dump%init(namespace, restart_gs, restart_type_dump, mc, ierr, mesh=gr)
    end if
 
    ! run self-consistency
    call scf_state_info(namespace, st)
 
    if (st%pack_states .and. hm%apply_packed()) then
      call st%pack()
    end if
 
    if (ks%theory_level == rdmft) then
      call rdmft_init(rdm, namespace, gr, st, hm, mc, space, fromscratch)
    end if
 
    if (.not. fromscratch) then
      call scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, scf%restart_load)
      call scf%restart_load%end()
    end if
    pop_sub(electrons_gs_initialize)
  end subroutine electrons_gs_initialize
 
  ! ---------------------------------------------------------
  subroutine electrons_gs_run(namespace, scf, rdm, mc, gr, ions, ext_partners, st, ks, hm, outp, space)
    type(namespace_t),        intent(in)    :: namespace
    type(scf_t),              intent(inout) :: scf
    type(rdm_t),              intent(inout) :: rdm
    type(multicomm_t),        intent(in)    :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(output_t),           intent(in)    :: outp
    type(electron_space_t),   intent(in)    :: space
 
    push_sub(electrons_gs_run)
 
    ! self-consistency for occupation numbers and natural orbitals in RDMFT
    if (ks%theory_level == rdmft) then
      call scf_rdmft(rdm, namespace, space, gr, ions, ext_partners, st, ks, hm, outp, scf%restart_dump)
    else
      if (scf%restart_dump%get_data_type() /= restart_undefined) then
        call scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp=outp, &
          restart_dump=scf%restart_dump)
      else
        call scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp=outp)
      end if
    end if
 
    pop_sub(electrons_gs_run)
  end subroutine electrons_gs_run
 
 
end module electrons_ground_state_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: