propagator_factory.F90

Go to the documentation of this file.

!! Copyright (C) 2019 N. Tancogne-Dejean
!! Copyright (C) 2020 M. Oliveira, Heiko Appel
!! Copyright (C) 2021 S. Ohlmann
!! Copyright (C) 2023 M. Oliveira
!!
!! 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 propagator_factory_oct_m
  use, intrinsic :: iso_fortran_env
  use algorithm_factory_oct_m
  use algorithm_oct_m
  use clock_oct_m
  use debug_oct_m
  use global_oct_m
  use interaction_partner_oct_m
  use math_oct_m
  use messages_oct_m
  use multisystem_basic_oct_m
  use namespace_oct_m
  use parser_oct_m
  use propagator_aetrs_oct_m
  use propagator_beeman_oct_m
  use propagator_bomd_oct_m
  use propagator_exp_mid_oct_m
  use propagator_exp_mid_2step_oct_m
  use propagator_exp_gauss1_oct_m
  use propagator_exp_gauss2_oct_m
  use propagator_leapfrog_oct_m
  use propagator_oct_m
  use propagator_rk4_oct_m
  use propagator_static_oct_m
  use propagator_verlet_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use system_oct_m
  implicit none
 
  private
  public :: propagator_factory_t
 
  ! Known multisystem propagators
  integer, public, parameter ::       &
    PROP_STATIC                 = 0,  &
    prop_verlet                 = 1,  &
    prop_beeman                 = 2,  &
    prop_beeman_scf             = 3,  &
    prop_expmid_2step           = 4,  &
    prop_expmid_2step_scf       = 5,  &
    prop_aetrs_ms               = 6,  &
    prop_rk4                    = 7,  &
    prop_expmid                 = 8,  &
    prop_leapfrog               = 9,  &
    prop_bomd                   = 10, &
    prop_expgauss1              = 11, &
    prop_expgauss2              = 12
 
 
  type, extends(algorithm_factory_t) :: propagator_factory_t
    private
    real(real64) :: final_time
  contains
    procedure :: create => propagator_factory_create               
    procedure :: create_static => propagator_factory_create_static 
  end type propagator_factory_t
 
  interface propagator_factory_t
    module procedure propagator_factory_constructor
  end interface propagator_factory_t
 
 
contains
  ! ---------------------------------------------------------------------------------------
  !
  function propagator_factory_constructor(namespace) result(factory)
    type(namespace_t), intent(in) :: namespace
    type(propagator_factory_t) :: factory
 
    push_sub(propagator_factory_constructor)
 
    ! Get final propagation time from input
    ! This variable is also defined (and properly documented) in td/td.F90.
    ! This is temporary, until all the propagators are moved to the new framework.
    call parse_variable(namespace, 'TDPropagationTime', -1.0_real64, factory%final_time, unit = units_inp%time)
    if (factory%final_time <= m_zero) then
      call messages_input_error(namespace, 'TDPropagationTime', 'must be greater than zero')
    end if
    call messages_print_var_value('TDPropagationTime', factory%final_time)
 
    pop_sub(propagator_factory_constructor)
  end function propagator_factory_constructor
 
  ! ---------------------------------------------------------------------------------------
  !
  function propagator_factory_create(this, system) result(algorithm)
    class(propagator_factory_t),  intent(in) :: this
    class(interaction_partner_t), intent(in), target :: system
    class(algorithm_t), pointer :: algorithm
 
    integer :: prop_type
    real(real64) :: dt
    class(propagator_t), pointer :: propagator
 
    push_sub(propagator_factory_create)
 
    !%Variable TDSystemPropagator
    !%Type integer
    !%Default static
    !%Section Time-Dependent::Propagation
    !%Description
    !% A variable to set the propagator in the multisystem framework.
    !% This is a temporary solution, and should be replaced by the
    !% TDPropagator variable.
    !%Option static 0
    !% (Experimental) Do not propagate the system in time.
    !%Option verlet 1
    !% (Experimental) Verlet propagator.
    !%Option beeman 2
    !% (Experimental) Beeman propagator without predictor-corrector.
    !%Option beeman_scf 3
    !% (Experimental) Beeman propagator with predictor-corrector scheme.
    !%Option exp_mid_2step 4
    !% (Experimental) Exponential midpoint propagator without predictor-corrector.
    !%Option exp_mid_2step_scf 5
    !% (Experimental) Exponential midpoint propagator with predictor-corrector scheme.
    !%Option prop_aetrs 6
    !% (Experimental) Approximate ETRS propagator
    !%Option prop_rk4 7
    !% (Experimental) RK4 propagator
    !%Option prop_expmid 8
    !% (Experimental) Exponential midpoint propagator with extrapolation.
    !%Option prop_leapfrog 9
    !% (Experimental) Leap frog algorithm
    !%Option prop_bomd 10
    !% (Experimental) Born-Oppenheimer MD propagator for the matter system.
    !%Option prop_exp_gauss1 11
    !% (Experimental) Exponential RK method with a Gauss quadrature rule with one point (second order)
    !% see also Hochbruck, M. & Ostermann, A.: Exponential Runge–Kutta methods for parabolic
    !% problems. Applied Numerical Mathematics 53, 323–339 (2005).
    !%Option prop_exp_gauss2 12
    !% (Experimental) Exponential RK method with a Gauss quadrature rule with two points (fourth order)
    !% see also Hochbruck, M. & Ostermann, A.: Exponential Runge–Kutta methods for parabolic
    !% problems. Applied Numerical Mathematics 53, 323–339 (2005).
    !%End
    call parse_variable(system%namespace, 'TDSystemPropagator', prop_static, prop_type)
    if (.not. varinfo_valid_option('TDSystemPropagator', prop_type)) then
      call messages_input_error(system%namespace, 'TDSystemPropagator')
    end if
    call messages_print_with_emphasis(msg='System propagator', namespace=system%namespace)
    call messages_print_var_option('TDSystemPropagator', prop_type, namespace=system%namespace)
 
    dt = propagator_factory_read_dt(this, system%namespace)
 
    select case (prop_type)
    case (prop_static)
      propagator => propagator_static_t(dt, 1)
    case (prop_verlet)
      propagator => propagator_verlet_t(dt)
    case (prop_beeman)
      propagator => propagator_beeman_t(dt, predictor_corrector=.false.)
    case (prop_beeman_scf)
      propagator => propagator_beeman_t(dt, predictor_corrector=.true.)
    case (prop_expmid_2step)
      propagator => propagator_exp_mid_2step_t(dt, predictor_corrector=.false.)
    case (prop_expmid_2step_scf)
      propagator => propagator_exp_mid_2step_t(dt, predictor_corrector=.true.)
    case (prop_aetrs_ms)
      propagator => propagator_aetrs_t(dt)
    case (prop_rk4)
      propagator => propagator_rk4_t(dt)
    case (prop_expmid)
      propagator => propagator_exp_mid_t(dt)
    case (prop_leapfrog)
      propagator => propagator_leapfrog_t(dt)
    case (prop_bomd)
      propagator => propagator_bomd_t(dt)
    case (prop_expgauss1)
      propagator => propagator_exp_gauss1_t(dt)
    case (prop_expgauss2)
      propagator => propagator_exp_gauss2_t(dt)
    case default
      call messages_input_error(system%namespace, 'TDSystemPropagator')
    end select
 
    call messages_print_with_emphasis(namespace=system%namespace)
 
    propagator%final_time = this%final_time
 
    select type (system)
    class is (system_t)
      propagator%system => system
    class default
      assert(.false.)
    end select
 
    algorithm => propagator
 
    pop_sub(propagator_factory_create)
  end function propagator_factory_create
 
  ! ---------------------------------------------------------------------------------------
  !
  function propagator_factory_create_static(this, system) result(algorithm)
    class(propagator_factory_t),  intent(in)         :: this
    class(interaction_partner_t), intent(in), target :: system
    class(algorithm_t), pointer :: algorithm
 
    real(real64) :: largest_dt, smallest_dt
    type(system_iterator_t) :: iter
    class(system_t), pointer :: subsystem
    class(propagator_t), pointer :: propagator
 
    push_sub(propagator_factory_create_static)
 
    ! Note: We need to pass the system as an interaction_partner_t, but we can only create a propagator
    ! for systems.
    select type (system)
    type is (multisystem_basic_t)
      ! The multisystem container is an exception, as its time-step and the
      ! number of algorithmic steps is determined by the subsystems.
      largest_dt = m_zero
      smallest_dt = m_huge
      call iter%start(system%list)
      do while (iter%has_next())
        subsystem => iter%get_next()
        select type (subalgorithm => subsystem%algo)
        class is (propagator_t)
          largest_dt = max(largest_dt, subalgorithm%dt)
          smallest_dt = min(smallest_dt, subalgorithm%dt/subalgorithm%algo_steps)
        end select
      end do
      propagator => propagator_static_t(largest_dt, nsteps=int(largest_dt/smallest_dt))
 
    class default
      propagator => propagator_static_t(propagator_factory_read_dt(this, system%namespace), 1)
    end select
 
    propagator%final_time = this%final_time
 
    select type (system)
    class is (system_t)
      propagator%system => system
    class default
      assert(.false.)
    end select
 
    algorithm => propagator
 
    pop_sub(propagator_factory_create_static)
  end function propagator_factory_create_static
 
  ! ---------------------------------------------------------------------------------------
  !
  real(real64) function propagator_factory_read_dt(this, namespace) result(dt)
    class(propagator_factory_t),  intent(in) :: this
    type(namespace_t),            intent(in) :: namespace
 
    real(real64) :: dt_adjusted
 
    ! This variable is also defined (and properly documented) in td/td.F90.
    ! This is temporary, until all the propagators are moved to the new framework.
    call parse_variable(namespace, 'TDTimeStep', 10.0_real64, dt)
    if (dt <= m_zero) then
      call messages_input_error(namespace, 'TDTimeStep', "must be greater than zero")
    end if
    call messages_print_var_value('TDTimeStep', dt, namespace=namespace)
 
    ! Time-step needs to be a multiple of the minimum clock time-step
    dt_adjusted = int(dt/clock_minimum_dt, int64) * clock_minimum_dt
    message(1) = "The time-step was adjusted to make it commensurable with the minimum time-step:"
    write(message(2), '(a,es21.14)') "  requested time-step = ", dt
    write(message(3), '(a,es21.14)') "  adjusted  time-step = ", dt_adjusted
    call messages_info(3, namespace=namespace)
    dt = dt_adjusted
 
  end function propagator_factory_read_dt
 
end module propagator_factory_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: