LCOV - code coverage report
Current view: top level - src/motion - rt_propagation.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:392804b) Lines: 212 218 97.2 %
Date: 2025-02-09 06:50:14 Functions: 6 6 100.0 %

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

Generated by: LCOV version 1.15