LCOV - code coverage report
Current view: top level - src/motion - rt_propagation.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:904d571) Lines: 198 204 97.1 %
Date: 2024-02-28 06:49:54 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines for the real time propagation.
      10             : !> \author Florian Schiffmann (02.09)
      11             : ! **************************************************************************************************
      12             : 
      13             : MODULE rt_propagation
      14             :    USE bibliography,                    ONLY: Andermatt2016,&
      15             :                                               cite_reference
      16             :    USE cp_control_types,                ONLY: dft_control_type,&
      17             :                                               rtp_control_type
      18             :    USE cp_external_control,             ONLY: external_control
      19             :    USE cp_fm_types,                     ONLY: cp_fm_set_all,&
      20             :                                               cp_fm_to_fm,&
      21             :                                               cp_fm_type
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_type,&
      24             :                                               cp_to_string
      25             :    USE cp_output_handling,              ONLY: cp_add_iter_level,&
      26             :                                               cp_iterate,&
      27             :                                               cp_print_key_unit_nr,&
      28             :                                               cp_rm_iter_level
      29             :    USE dbcsr_api,                       ONLY: dbcsr_copy,&
      30             :                                               dbcsr_p_type
      31             :    USE efield_utils,                    ONLY: calculate_ecore_efield
      32             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      33             :    USE force_env_types,                 ONLY: force_env_get,&
      34             :                                               force_env_type
      35             :    USE global_types,                    ONLY: global_environment_type
      36             :    USE hfx_admm_utils,                  ONLY: hfx_admm_init
      37             :    USE input_constants,                 ONLY: real_time_propagation,&
      38             :                                               use_restart_wfn,&
      39             :                                               use_rt_restart,&
      40             :                                               use_scf_wfn
      41             :    USE input_cp2k_restarts,             ONLY: write_restart
      42             :    USE input_section_types,             ONLY: section_vals_get,&
      43             :                                               section_vals_get_subs_vals,&
      44             :                                               section_vals_type,&
      45             :                                               section_vals_val_get,&
      46             :                                               section_vals_val_set
      47             :    USE kinds,                           ONLY: dp
      48             :    USE machine,                         ONLY: m_walltime
      49             :    USE md_environment_types,            ONLY: md_environment_type
      50             :    USE pw_env_types,                    ONLY: pw_env_type
      51             :    USE qs_core_hamiltonian,             ONLY: qs_matrix_h_allocate_imag_from_real
      52             :    USE qs_energy_init,                  ONLY: qs_energies_init
      53             :    USE qs_energy_types,                 ONLY: qs_energy_type
      54             :    USE qs_environment_types,            ONLY: get_qs_env,&
      55             :                                               qs_environment_type
      56             :    USE qs_external_potential,           ONLY: external_c_potential,&
      57             :                                               external_e_potential
      58             :    USE qs_ks_methods,                   ONLY: qs_ks_allocate_basics,&
      59             :                                               qs_ks_update_qs_env
      60             :    USE qs_ks_types,                     ONLY: qs_ks_did_change,&
      61             :                                               qs_ks_env_type,&
      62             :                                               set_ks_env
      63             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      64             :                                               init_mo_set,&
      65             :                                               mo_set_type
      66             :    USE qs_rho_methods,                  ONLY: allocate_rho_ao_imag_from_real
      67             :    USE qs_rho_types,                    ONLY: qs_rho_set,&
      68             :                                               qs_rho_type
      69             :    USE rt_delta_pulse,                  ONLY: apply_delta_pulse
      70             :    USE rt_hfx_utils,                    ONLY: rtp_hfx_rebuild
      71             :    USE rt_projection_mo_utils,          ONLY: init_mo_projection
      72             :    USE rt_propagation_methods,          ONLY: propagation_step
      73             :    USE rt_propagation_output,           ONLY: rt_prop_output
      74             :    USE rt_propagation_types,            ONLY: get_rtp,&
      75             :                                               rt_prop_create,&
      76             :                                               rt_prop_type,&
      77             :                                               rtp_create_SinvH_imag,&
      78             :                                               rtp_history_create
      79             :    USE rt_propagation_utils,            ONLY: calc_S_derivs,&
      80             :                                               calc_update_rho,&
      81             :                                               calc_update_rho_sparse,&
      82             :                                               get_restart_wfn
      83             :    USE rt_propagation_velocity_gauge,   ONLY: velocity_gauge_ks_matrix
      84             :    USE rt_propagator_init,              ONLY: init_propagators,&
      85             :                                               rt_initialize_rho_from_mos
      86             : #include "../base/base_uses.f90"
      87             : 
      88             :    IMPLICIT NONE
      89             : 
      90             :    PRIVATE
      91             : 
      92             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation'
      93             : 
      94             :    PUBLIC :: rt_prop_setup
      95             : 
      96             : CONTAINS
      97             : 
      98             : ! **************************************************************************************************
      99             : !> \brief creates rtp_type, gets the initial state, either by reading MO's
     100             : !>        from file or calling SCF run
     101             : !> \param force_env ...
     102             : !> \author Florian Schiffmann (02.09)
     103             : ! **************************************************************************************************
     104             : 
     105         594 :    SUBROUTINE rt_prop_setup(force_env)
     106             :       TYPE(force_env_type), POINTER                      :: force_env
     107             : 
     108             :       INTEGER                                            :: aspc_order
     109             :       LOGICAL                                            :: magnetic, vel_reprs
     110             :       TYPE(dft_control_type), POINTER                    :: dft_control
     111             :       TYPE(global_environment_type), POINTER             :: globenv
     112             :       TYPE(qs_energy_type), POINTER                      :: energy
     113             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     114             :       TYPE(rt_prop_type), POINTER                        :: rtp
     115             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     116             :       TYPE(section_vals_type), POINTER                   :: hfx_sections, input, ls_scf_section, &
     117             :                                                             md_section, motion_section, &
     118             :                                                             print_moments_section
     119             : 
     120         198 :       NULLIFY (qs_env, rtp_control, dft_control)
     121             : 
     122         198 :       CALL cite_reference(Andermatt2016)
     123             : 
     124         198 :       CALL force_env_get(force_env=force_env, qs_env=qs_env, globenv=globenv)
     125         198 :       CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
     126         198 :       rtp_control => dft_control%rtp_control
     127             : 
     128             :       ! Takes care that an initial wavefunction/density is available
     129             :       ! Can either be by performing an scf loop or reading a restart
     130         198 :       CALL rt_initial_guess(qs_env, force_env, rtp_control)
     131             : 
     132             :       ! Initializes the extrapolation
     133         198 :       NULLIFY (rtp)
     134         198 :       CALL get_qs_env(qs_env=qs_env, rtp=rtp, input=input)
     135         198 :       aspc_order = rtp_control%aspc_order
     136         198 :       CALL rtp_history_create(rtp, aspc_order)
     137             : 
     138             :       ! Reads the simulation parameters from the input
     139         198 :       motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
     140         198 :       md_section => section_vals_get_subs_vals(motion_section, "MD")
     141         198 :       hfx_sections => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%XC%HF")
     142         198 :       print_moments_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%PRINT%MOMENTS")
     143         198 :       CALL section_vals_val_get(md_section, "TIMESTEP", r_val=qs_env%rtp%dt)
     144         198 :       CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=qs_env%rtp%i_start)
     145         198 :       CALL section_vals_val_get(md_section, "STEPS", i_val=rtp%nsteps)
     146         198 :       CALL section_vals_val_get(md_section, "MAX_STEPS", i_val=rtp%max_steps)
     147             : 
     148         198 :       ls_scf_section => section_vals_get_subs_vals(input, "DFT%LS_SCF")
     149         198 :       CALL section_vals_val_get(ls_scf_section, "EPS_FILTER", r_val=rtp%filter_eps)
     150         198 :       IF (.NOT. qs_env%rtp%linear_scaling) rtp%filter_eps = 0.0_dp
     151         198 :       IF (rtp_control%acc_ref < 1) rtp_control%acc_ref = 1
     152         198 :       rtp%filter_eps_small = rtp%filter_eps/rtp_control%acc_ref
     153         198 :       CALL section_vals_val_get(ls_scf_section, "EPS_LANCZOS", r_val=rtp%lanzcos_threshold)
     154         198 :       CALL section_vals_val_get(ls_scf_section, "MAX_ITER_LANCZOS", i_val=rtp%lanzcos_max_iter)
     155         198 :       CALL section_vals_val_get(ls_scf_section, "SIGN_SQRT_ORDER", i_val=rtp%newton_schulz_order)
     156         198 :       CALL section_vals_get(hfx_sections, explicit=rtp%do_hfx)
     157         198 :       CALL section_vals_val_get(print_moments_section, "MAGNETIC", l_val=magnetic)
     158         198 :       CALL section_vals_val_get(print_moments_section, "VEL_REPRS", l_val=vel_reprs)
     159             : 
     160             :       rtp%track_imag_density = (magnetic) .OR. (vel_reprs) .OR. (rtp_control%velocity_gauge) &
     161         198 :                                .OR. (rtp%do_hfx) .OR. (.NOT. rtp_control%fixed_ions)
     162         198 :       rtp%propagate_complex_ks = rtp%do_hfx .OR. rtp_control%velocity_gauge
     163             : 
     164             :       CALL rt_init_complex_quantities(qs_env, imag_p=rtp%track_imag_density, &
     165         198 :                                       imag_ks=rtp%propagate_complex_ks, imag_h=rtp_control%velocity_gauge)
     166             : 
     167             :       ! Hmm, not really like to initialize with the structure of S but I reckon it is
     168             :       ! done everywhere like this
     169         198 :       IF (rtp%do_hfx) CALL rtp_hfx_rebuild(qs_env)
     170             : 
     171             :       ! Setup the MO projection environment if required
     172         198 :       IF (rtp_control%is_proj_mo) CALL init_mo_projection(qs_env, rtp_control)
     173             : 
     174         198 :       CALL init_propagation_run(qs_env)
     175         198 :       IF (.NOT. rtp_control%fixed_ions) THEN
     176             :          !derivativs of the overlap needed for EMD
     177          72 :          CALL calc_S_derivs(qs_env)
     178             :          ! a bit hidden, but computes SinvH and SinvB (calc_SinvH for CN,EM and ARNOLDI)
     179             :          ! make_etrs_exp in case of ETRS in combination with TAYLOR and PADE
     180             :       END IF
     181         198 :       CALL init_propagators(qs_env)
     182         198 :       IF (rtp_control%fixed_ions) THEN
     183         126 :          CALL run_propagation(qs_env, force_env, globenv)
     184             :       ELSE
     185          72 :          rtp_control%initial_step = .TRUE.
     186          72 :          CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
     187          72 :          rtp_control%initial_step = .FALSE.
     188          72 :          rtp%energy_old = energy%total
     189             :       END IF
     190             : 
     191         198 :    END SUBROUTINE rt_prop_setup
     192             : 
     193             : ! **************************************************************************************************
     194             : !> \brief calculates the matrices needed in the first step of EMD/RTP
     195             : !> \param qs_env ...
     196             : !> \author Florian Schiffmann (02.09)
     197             : ! **************************************************************************************************
     198             : 
     199         198 :    SUBROUTINE init_propagation_run(qs_env)
     200             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     201             : 
     202             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
     203             : 
     204             :       INTEGER                                            :: i, ispin, re
     205             :       INTEGER, DIMENSION(2)                              :: nelectron_spin
     206         198 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new, mos_old
     207         198 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, rho_new, rho_old
     208             :       TYPE(dft_control_type), POINTER                    :: dft_control
     209         198 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     210             :       TYPE(rt_prop_type), POINTER                        :: rtp
     211             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     212             : 
     213         198 :       NULLIFY (dft_control, rtp, rtp_control)
     214             : 
     215         198 :       CALL cite_reference(Andermatt2016)
     216             : 
     217             :       CALL get_qs_env(qs_env, &
     218             :                       rtp=rtp, &
     219         198 :                       dft_control=dft_control)
     220         198 :       rtp_control => dft_control%rtp_control
     221             : 
     222         198 :       IF (rtp_control%initial_wfn == use_scf_wfn) THEN
     223         162 :          IF (rtp_control%apply_delta_pulse .OR. rtp_control%apply_delta_pulse_mag) THEN
     224          54 :             CALL apply_delta_pulse(qs_env, rtp, rtp_control)
     225             :          ELSE
     226         108 :             IF (.NOT. rtp%linear_scaling) THEN
     227          74 :                CALL get_rtp(rtp=rtp, mos_old=mos_old)
     228          74 :                CALL get_qs_env(qs_env, mos=mos)
     229         162 :                DO i = 1, SIZE(mos)
     230          88 :                   CALL cp_fm_to_fm(mos(i)%mo_coeff, mos_old(2*i - 1))
     231         162 :                   CALL cp_fm_set_all(mos_old(2*i), zero, zero)
     232             :                END DO
     233             :             END IF
     234             :          END IF
     235             :       END IF
     236             : 
     237         198 :       IF (.NOT. rtp%linear_scaling) THEN
     238         108 :          CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
     239         384 :          DO i = 1, SIZE(mos_old)
     240         384 :             CALL cp_fm_to_fm(mos_old(i), mos_new(i))
     241             :          END DO
     242         108 :          CALL calc_update_rho(qs_env)
     243             :       ELSE
     244          90 :          IF (rtp_control%initial_wfn == use_scf_wfn) THEN
     245             :             CALL get_qs_env(qs_env, &
     246             :                             matrix_ks=matrix_ks, &
     247             :                             mos=mos, &
     248          74 :                             nelectron_spin=nelectron_spin)
     249          74 :             IF (ASSOCIATED(mos)) THEN
     250             :                !The wavefunction was minimized by an mo based algorith. P is therefore calculated from the mos
     251          64 :                IF (ASSOCIATED(rtp%mos)) THEN
     252          40 :                   IF (ASSOCIATED(rtp%mos%old)) THEN
     253             :                      ! Delta kick was applied and the results is in rtp%mos%old
     254          40 :                      CALL rt_initialize_rho_from_mos(rtp, mos, mos_old=rtp%mos%old)
     255             :                   ELSE
     256           0 :                      CALL rt_initialize_rho_from_mos(rtp, mos)
     257             :                   END IF
     258             :                ELSE
     259          24 :                   CALL rt_initialize_rho_from_mos(rtp, mos)
     260             :                END IF
     261             :             ELSE
     262             :                !The wavefunction was minimized using a linear scaling method. The density matrix is therefore taken from the ls_scf_env.
     263          10 :                CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
     264          24 :                DO ispin = 1, SIZE(rho_old)/2
     265          14 :                   re = 2*ispin - 1
     266          14 :                   CALL dbcsr_copy(rho_old(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
     267          24 :                   CALL dbcsr_copy(rho_new(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
     268             :                END DO
     269             :             END IF
     270          74 :             CALL calc_update_rho_sparse(qs_env)
     271             :          END IF
     272             :       END IF
     273             :       ! Modify KS matrix to include the additional terms in the velocity gauge
     274         198 :       IF (rtp_control%velocity_gauge) THEN
     275             :          ! As matrix_h and matrix_h_im are not updated by qs_ks_update_qs_env()
     276             :          ! the non-gauge transformed non-local part has to be subtracted here
     277           8 :          CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.TRUE.)
     278             :       END IF
     279         198 :       CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
     280             : 
     281         198 :    END SUBROUTINE init_propagation_run
     282             : 
     283             : ! **************************************************************************************************
     284             : !> \brief performs the real RTP run, gets information from MD section
     285             : !>        uses MD as iteration level
     286             : !> \param qs_env ...
     287             : !> \param force_env ...
     288             : !> \param globenv ...
     289             : !> \author Florian Schiffmann (02.09)
     290             : ! **************************************************************************************************
     291             : 
     292         126 :    SUBROUTINE run_propagation(qs_env, force_env, globenv)
     293             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     294             :       TYPE(force_env_type), POINTER                      :: force_env
     295             :       TYPE(global_environment_type), POINTER             :: globenv
     296             : 
     297             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'run_propagation'
     298             : 
     299             :       INTEGER                                            :: aspc_order, handle, i_iter, i_step, &
     300             :                                                             max_iter, max_steps, output_unit
     301             :       LOGICAL                                            :: should_stop
     302             :       REAL(Kind=dp)                                      :: eps_ener, time_iter_start, &
     303             :                                                             time_iter_stop, used_time
     304             :       TYPE(cp_logger_type), POINTER                      :: logger
     305             :       TYPE(dft_control_type), POINTER                    :: dft_control
     306             :       TYPE(pw_env_type), POINTER                         :: pw_env
     307             :       TYPE(qs_energy_type), POINTER                      :: energy
     308             :       TYPE(rt_prop_type), POINTER                        :: rtp
     309             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     310             :       TYPE(section_vals_type), POINTER                   :: input, rtp_section
     311             : 
     312         126 :       should_stop = .FALSE.
     313         126 :       CALL timeset(routineN, handle)
     314             : 
     315         126 :       CALL cite_reference(Andermatt2016)
     316             : 
     317         126 :       NULLIFY (logger, dft_control, energy, rtp, rtp_control, input, rtp_section)
     318         126 :       logger => cp_get_default_logger()
     319             : 
     320         126 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, rtp=rtp, energy=energy, input=input)
     321             : 
     322         126 :       rtp_control => dft_control%rtp_control
     323         126 :       max_steps = MIN(rtp%nsteps, rtp%max_steps)
     324         126 :       max_iter = rtp_control%max_iter
     325         126 :       eps_ener = rtp_control%eps_ener
     326             : 
     327         126 :       aspc_order = rtp_control%aspc_order
     328             : 
     329         126 :       rtp%energy_old = energy%total
     330         126 :       time_iter_start = m_walltime()
     331         126 :       CALL cp_add_iter_level(logger%iter_info, "MD")
     332         126 :       CALL cp_iterate(logger%iter_info, iter_nr=0)
     333         126 :       IF (rtp%i_start >= max_steps) CALL cp_abort(__LOCATION__, &
     334           0 :                                                   "maximum step number smaller than initial step value")
     335             : 
     336         126 :       rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
     337             :       output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
     338         126 :                                          extension=".scfLog")
     339             : 
     340         470 :       DO i_step = rtp%i_start + 1, max_steps
     341         344 :          IF (output_unit > 0) THEN
     342             :             WRITE (output_unit, FMT="(/,(T2,A,T40,I6))") &
     343         172 :                "Real time propagation step:", i_step
     344             :          END IF
     345         344 :          energy%efield_core = 0.0_dp
     346         344 :          qs_env%sim_time = REAL(i_step, dp)*rtp%dt
     347         344 :          CALL get_qs_env(qs_env, pw_env=pw_env)
     348         344 :          pw_env%poisson_env%parameters%dbc_params%time = qs_env%sim_time
     349         344 :          qs_env%sim_step = i_step
     350         344 :          rtp%istep = i_step - rtp%i_start
     351         344 :          CALL calculate_ecore_efield(qs_env, .FALSE.)
     352         344 :          IF (dft_control%apply_external_potential) THEN
     353           0 :             IF (.NOT. dft_control%expot_control%static) THEN
     354           0 :                dft_control%eval_external_potential = .TRUE.
     355             :             END IF
     356             :          END IF
     357         344 :          CALL external_c_potential(qs_env, calculate_forces=.FALSE.)
     358         344 :          CALL external_e_potential(qs_env)
     359         344 :          CALL cp_iterate(logger%iter_info, last=(i_step == max_steps), iter_nr=i_step)
     360         344 :          rtp%converged = .FALSE.
     361        1142 :          DO i_iter = 1, max_iter
     362        1142 :             IF (i_step == rtp%i_start + 1 .AND. i_iter == 2 .AND. rtp_control%hfx_redistribute) &
     363           0 :                CALL qs_ks_did_change(qs_env%ks_env, s_mstruct_changed=.TRUE.)
     364        1142 :             rtp%iter = i_iter
     365        1142 :             CALL propagation_step(qs_env, rtp, rtp_control)
     366        1142 :             CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
     367        1142 :             rtp%energy_new = energy%total
     368        1142 :             IF (rtp%converged) EXIT
     369        1142 :             CALL rt_prop_output(qs_env, real_time_propagation, rtp%delta_iter)
     370             :          END DO
     371         470 :          IF (rtp%converged) THEN
     372         344 :             CALL external_control(should_stop, "MD", globenv=globenv)
     373         344 :             IF (should_stop) CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=i_step)
     374         344 :             time_iter_stop = m_walltime()
     375         344 :             used_time = time_iter_stop - time_iter_start
     376         344 :             time_iter_start = time_iter_stop
     377         344 :             CALL rt_prop_output(qs_env, real_time_propagation, delta_iter=rtp%delta_iter, used_time=used_time)
     378         344 :             CALL rt_write_input_restart(force_env=force_env, qs_env=qs_env)
     379         344 :             IF (should_stop) EXIT
     380             :          ELSE
     381             :             EXIT
     382             :          END IF
     383             :       END DO
     384         126 :       CALL cp_rm_iter_level(logger%iter_info, "MD")
     385             : 
     386         126 :       IF (.NOT. rtp%converged) &
     387             :          CALL cp_abort(__LOCATION__, "propagation did not converge, "// &
     388           0 :                        "either increase MAX_ITER or use a smaller TIMESTEP")
     389             : 
     390         126 :       CALL timestop(handle)
     391             : 
     392         126 :    END SUBROUTINE run_propagation
     393             : 
     394             : ! **************************************************************************************************
     395             : !> \brief overwrites some values in the input file such that the .restart
     396             : !>        file will contain the appropriate information
     397             : !> \param md_env ...
     398             : !> \param qs_env ...
     399             : !> \param force_env ...
     400             : !> \author Florian Schiffmann (02.09)
     401             : ! **************************************************************************************************
     402             : 
     403         344 :    SUBROUTINE rt_write_input_restart(md_env, qs_env, force_env)
     404             :       TYPE(md_environment_type), OPTIONAL, POINTER       :: md_env
     405             :       TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
     406             :       TYPE(force_env_type), POINTER                      :: force_env
     407             : 
     408         344 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_vals
     409             :       TYPE(dft_control_type), POINTER                    :: dft_control
     410             :       TYPE(section_vals_type), POINTER                   :: efield_section, motion_section, &
     411             :                                                             root_section, rt_section
     412             : 
     413         344 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     414         344 :       root_section => force_env%root_section
     415         344 :       motion_section => section_vals_get_subs_vals(root_section, "MOTION")
     416         344 :       rt_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION")
     417         344 :       CALL section_vals_val_set(rt_section, "INITIAL_WFN", i_val=use_rt_restart)
     418             :       ! coming from RTP
     419         344 :       IF (.NOT. PRESENT(md_env)) THEN
     420         344 :          CALL section_vals_val_set(motion_section, "MD%STEP_START_VAL", i_val=force_env%qs_env%sim_step)
     421             :       END IF
     422             : 
     423         344 :       IF (dft_control%apply_vector_potential) THEN
     424          22 :          efield_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%EFIELD")
     425             :          NULLIFY (tmp_vals)
     426          22 :          ALLOCATE (tmp_vals(3))
     427          88 :          tmp_vals = dft_control%efield_fields(1)%efield%vec_pot_initial
     428             :          CALL section_vals_val_set(efield_section, "VEC_POT_INITIAL", &
     429             :                                    r_vals_ptr=tmp_vals, &
     430          22 :                                    i_rep_section=1)
     431             :       END IF
     432             : 
     433         344 :       CALL write_restart(md_env=md_env, root_section=root_section)
     434             : 
     435         344 :    END SUBROUTINE rt_write_input_restart
     436             : 
     437             : ! **************************************************************************************************
     438             : !> \brief Creates the initial electronic states and allocates the necessary
     439             : !>        matrices
     440             : !> \param qs_env ...
     441             : !> \param force_env ...
     442             : !> \param rtp_control ...
     443             : !> \author Florian Schiffmann (02.09)
     444             : ! **************************************************************************************************
     445             : 
     446         198 :    SUBROUTINE rt_initial_guess(qs_env, force_env, rtp_control)
     447             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     448             :       TYPE(force_env_type), POINTER                      :: force_env
     449             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     450             : 
     451             :       INTEGER                                            :: homo, ispin
     452             :       LOGICAL                                            :: energy_consistency
     453             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     454         198 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     455             :       TYPE(dft_control_type), POINTER                    :: dft_control
     456             : 
     457         198 :       NULLIFY (matrix_s, dft_control)
     458         198 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     459             : 
     460         360 :       SELECT CASE (rtp_control%initial_wfn)
     461             :       CASE (use_scf_wfn)
     462         162 :          qs_env%sim_time = 0.0_dp
     463         162 :          qs_env%sim_step = 0
     464         162 :          energy_consistency = .TRUE.
     465             :          !in the linear scaling case we need a correct kohn-sham matrix, which we cannot get with consistent energies
     466         162 :          IF (rtp_control%linear_scaling) energy_consistency = .FALSE.
     467             :          CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., &
     468         162 :                                           consistent_energies=energy_consistency)
     469         162 :          qs_env%run_rtp = .TRUE.
     470         162 :          ALLOCATE (qs_env%rtp)
     471         162 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     472         162 :          IF (dft_control%do_admm) THEN
     473           6 :             CALL hfx_admm_init(qs_env)
     474             :             CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     475           6 :                                 rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
     476             :          ELSE
     477             :             CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     478         156 :                                 rtp_control%linear_scaling)
     479             :          END IF
     480             : 
     481             :       CASE (use_restart_wfn, use_rt_restart)
     482          36 :          CALL qs_energies_init(qs_env, .FALSE.)
     483          36 :          IF (.NOT. rtp_control%linear_scaling .OR. rtp_control%initial_wfn == use_restart_wfn) THEN
     484          86 :             DO ispin = 1, SIZE(qs_env%mos)
     485          52 :                CALL get_mo_set(qs_env%mos(ispin), mo_coeff=mo_coeff, homo=homo)
     486          86 :                IF (.NOT. ASSOCIATED(mo_coeff)) THEN
     487             :                   CALL init_mo_set(qs_env%mos(ispin), &
     488             :                                    qs_env%mpools%ao_mo_fm_pools(ispin)%pool, &
     489          52 :                                    name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
     490             :                END IF
     491             :             END DO
     492          34 :             IF (dft_control%do_admm) CALL hfx_admm_init(qs_env)
     493             :          END IF
     494          36 :          ALLOCATE (qs_env%rtp)
     495          36 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     496             :          CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     497          36 :                              rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
     498          36 :          CALL get_restart_wfn(qs_env)
     499             : 
     500         234 :          qs_env%run_rtp = .TRUE.
     501             :       END SELECT
     502             : 
     503         198 :    END SUBROUTINE rt_initial_guess
     504             : 
     505             : ! **************************************************************************************************
     506             : !> \brief ...
     507             : !> \param qs_env ...
     508             : !> \param imag_p ...
     509             : !> \param imag_ks ...
     510             : !> \param imag_h ...
     511             : ! **************************************************************************************************
     512         198 :    SUBROUTINE rt_init_complex_quantities(qs_env, imag_p, imag_ks, imag_h)
     513             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     514             :       LOGICAL, INTENT(in)                                :: imag_p, imag_ks, imag_h
     515             : 
     516             :       TYPE(dft_control_type), POINTER                    :: dft_control
     517             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     518             :       TYPE(qs_rho_type), POINTER                         :: rho
     519             :       TYPE(rt_prop_type), POINTER                        :: rtp
     520             : 
     521         198 :       NULLIFY (ks_env, rho, dft_control)
     522             : 
     523             :       CALL get_qs_env(qs_env, &
     524             :                       dft_control=dft_control, &
     525             :                       ks_env=ks_env, &
     526             :                       rho=rho, &
     527         198 :                       rtp=rtp)
     528             : 
     529             :       ! rho
     530         198 :       CALL qs_rho_set(rho, complex_rho_ao=imag_p)
     531         198 :       IF (imag_p) CALL allocate_rho_ao_imag_from_real(rho, qs_env)
     532             : 
     533             :       ! ks
     534         198 :       CALL set_ks_env(ks_env, complex_ks=imag_ks)
     535         198 :       IF (imag_ks) THEN
     536          40 :          CALL qs_ks_allocate_basics(qs_env, is_complex=imag_ks)
     537          40 :          IF (.NOT. dft_control%rtp_control%fixed_ions) &
     538          22 :             CALL rtp_create_SinvH_imag(rtp, dft_control%nspins)
     539             :       END IF
     540             : 
     541             :       ! h
     542         198 :       IF (imag_h) CALL qs_matrix_h_allocate_imag_from_real(qs_env)
     543             : 
     544         198 :    END SUBROUTINE rt_init_complex_quantities
     545             : 
     546             : END MODULE rt_propagation

Generated by: LCOV version 1.15