eigen_evolution.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2025 S. Ohlmann
!!
!! 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 eigen_evolution_oct_m
  use debug_oct_m
  use exponential_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use ks_potential_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public :: &
    eigensolver_evolution
 
  integer, public, parameter :: &
    IT_FORWARD_EULER = 1, &
    it_expmid = 2
 
contains
  
  subroutine eigensolver_evolution(namespace, mesh, st, hm, space, ext_partners, te, iter, niter, converged, &
    tau0, tau, time, propagator, variable_timestep, vks_old, normalization_energies, normalization_energies_prev, residual)
    type(namespace_t),                intent(in)    :: namespace
    type(mesh_t),             target, intent(in)    :: mesh
    type(states_elec_t),      target, intent(inout) :: st
    type(hamiltonian_elec_t), target, intent(inout) :: hm
    class(space_t),                   intent(in)    :: space
    type(partner_list_t),             intent(in)    :: ext_partners
    type(exponential_t),              intent(inout) :: te
    integer,                          intent(in)    :: iter
    integer,                          intent(inout) :: niter
    integer,                          intent(inout) :: converged(:)
    real(real64),                     intent(inout) :: tau0
    real(real64),                     intent(inout) :: tau
    real(real64),                     intent(inout) :: time
    integer,                          intent(in)    :: propagator
    logical,                          intent(in)    :: variable_timestep
    type(potential_interpolation_t),  intent(inout) :: vks_old
    real(real64),                     intent(inout) :: normalization_energies(:, :)
    real(real64),                     intent(inout) :: normalization_energies_prev(:, :)
    real(real64),                     intent(inout) :: residual(:, :)
 
    real(real64) :: delta_normalization, delta_expectation, delta_residual
    real(real64) :: weights(1:st%nst, 1:st%nik)
    real(real64), parameter :: gam = 1.0_real64
    real(real64), parameter :: safety = 1.0e-2_real64
    integer :: ist, ik
 
    push_sub(eigensolver_evolution)
 
    ! Warning: it seems that the algorithm is improved if some extra states are added -- states
    ! whose convergence should not be checked.
 
    write(message(1), "(A, I5, A, F10.4, A, F10.4)") "Evolution eigensolver: starting iteration ", iter, &
      " at imaginary time ", time, " with timestep ", tau
    call messages_info(1)
 
    if (variable_timestep .and. iter > 2) then
      do ik = 1, st%nik
        do ist = 1, st%nst
          weights(ist, ik) = st%kweights(ik)*st%occ(ist, ik)
        end do
      end do
      ! The variable timestep method is implemented as described in
      ! Aichinger, M. & Krotscheck, E., Computational Materials Science 34, 188–212 (2005), 10.1016/j.commatsci.2004.11.002
      ! section 2.5
      ! eq. 24
      delta_normalization = norm2(weights*(normalization_energies-normalization_energies_prev))/sqrt(st%qtot)
      ! eq. 25
      delta_expectation = norm2(weights*(normalization_energies-st%eigenval))/sqrt(st%qtot)
      ! also take the residual into account; not in the paper
      delta_residual = norm2(weights*residual)/sqrt(st%qtot)
 
      write(message(1), "(A, ES15.2)") "Normalization energy error: ", delta_normalization
      write(message(2), "(A, ES15.2)") "Expectation energy error:   ", delta_expectation
      write(message(3), "(A, ES15.2)") "Residual:                   ", delta_residual
      call messages_info(3, debug_only=.true.)
 
      ! reduce timestep if normalization energy changes less after one step compared to the
      ! expectation energy; moreover, it is not allowed to be too small compared to the residual
      ! to avoid getting "convergence" due to very small timesteps, but with large residuals
      ! the first condition is eq. 26 in the paper, the second is added by us
      ! Moreover, let the timestep never be less than the initial timestep divided by 20
      if (delta_normalization < gam*delta_expectation .and. delta_normalization > delta_residual * safety &
        .and. tau/m_two >= tau0/20._real64) then
        tau = tau/m_two
        write(message(1), "(A, F10.4)") "Evolution eigensolver: reducing timestep to ", tau
        call messages_info(1)
      end if
    end if
 
    if (variable_timestep) then
      normalization_energies_prev = normalization_energies
    end if
 
    select case(propagator)
    case(it_forward_euler)
      call exponential_imaginary_time(namespace, mesh, st, hm, te, converged, tau, normalization_energies)
    case(it_expmid)
      call hm%ks_pot%interpolation_new(vks_old, time+tau, tau)
      call hm%ks_pot%interpolate_potentials(vks_old, 3, time+tau, tau, time + tau/m_two)
      call hm%update(mesh, namespace, space, ext_partners, time=time+tau/m_two)
      call exponential_imaginary_time(namespace, mesh, st, hm, te, converged, tau, normalization_energies)
    case default
      message(1) = "Unknown propagator for imaginary-time evolution."
      call messages_fatal(1)
    end select
 
 
    ! we cannot estimate the number of matrix-vector products because it depends on the algorithm
    niter = 0
 
    time = time + tau
 
    pop_sub(eigensolver_evolution)
  end subroutine eigensolver_evolution
 
  subroutine exponential_imaginary_time(namespace, mesh, st, hm, te, converged, tau, normalization_energies)
    type(namespace_t),                intent(in)    :: namespace
    type(mesh_t),             target, intent(in)    :: mesh
    type(states_elec_t),      target, intent(inout) :: st
    type(hamiltonian_elec_t), target, intent(inout) :: hm
    type(exponential_t),              intent(inout) :: te
    integer,                          intent(inout) :: converged(:)
    real(real64),                     intent(in)    :: tau
    real(real64),                     intent(inout) :: normalization_energies(:, :)
 
    integer :: ib, convb, ik
    type(wfs_elec_t), pointer :: zpsib
    complex(real64), allocatable :: dot(:, :)
 
    push_sub(exponential_imaginary_time)
 
    safe_allocate(dot(1:st%nst, 1:st%nik))
    dot = m_zero
 
    do ik = st%d%kpt%start, st%d%kpt%end
      convb = get_first_unconverged_batch(st, hm, ik, converged)
 
      do ib = convb + 1, st%group%block_end
        if (st%group%psib(ib, ik)%type() == type_float) then
          ! we need to copy to a complex batch for real states to apply the exponential
          safe_allocate(zpsib)
          call st%group%psib(ib, ik)%copy_to(zpsib, copy_data=.true., dest_type=type_cmplx)
        else
          zpsib => st%group%psib(ib, ik)
        end if
 
        ! simply apply the exponential here
        ! normalization is done in the subspace diagonalization call
        call hm%phase%set_phase_corr(mesh, zpsib)
        call te%apply_batch(namespace, mesh, hm, zpsib, tau, imag_time=.true.)
        call hm%phase%unset_phase_corr(mesh, zpsib)
 
        ! estimate normalization energy for variable timesteps
        call zmesh_batch_dotp_vector(mesh, zpsib, zpsib, &
          dot(states_elec_block_min(st, ib):states_elec_block_max(st, ib), ik), reduce=.false.)
 
        if (st%group%psib(ib, ik)%type() == type_float) then
          ! we need to copy back the data for real states
          call zpsib%copy_data_to(mesh%np, st%group%psib(ib, ik))
          call zpsib%end()
          safe_deallocate_p(zpsib)
        else
          nullify(zpsib)
        end if
      end do
    end do
 
    call st%dom_st_kpt_mpi_grp%allreduce_inplace(dot, st%nst*st%nik, mpi_double_complex, mpi_sum)
    ! estimate energies from change in normalization
    normalization_energies = -log(abs(dot))/(m_two*tau)
 
    safe_deallocate_a(dot)
 
    pop_sub(exponential_imaginary_time)
  end subroutine exponential_imaginary_time
 
  function get_first_unconverged_batch(st, hm, ik, converged) result(convb)
    type(states_elec_t),      intent(in) :: st
    type(hamiltonian_elec_t), intent(in) :: hm
    integer,                  intent(in) :: ik
    integer,                  intent(in) :: converged(:)
    integer :: convb
 
    integer :: conv, ib
 
    push_sub(get_first_unconverged_batch)
 
    convb = st%group%block_start - 1
    ! In case of independent particles, the Hamiltonian does not depends upon the density and we can freeze
    ! the lowest converged states
    if (hm%theory_level == independent_particles) then
      conv = 0
      do ib = st%group%block_start, st%group%block_end
        if (conv + st%group%psib(ib, ik)%nst <= converged(ik)) then
          conv = conv + st%group%psib(ib, ik)%nst
          convb = convb + 1
        else
          exit
        end if
      end do
    end if
 
    pop_sub(get_first_unconverged_batch)
  end function get_first_unconverged_batch
end module eigen_evolution_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: