LCOV - code coverage report
Current view: top level - src/motion - rt_propagation.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.3 % 220 214
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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         1140 :          DO i_iter = 1, max_iter
     382         1140 :             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         1140 :             rtp%iter = i_iter
     385         1140 :             CALL propagation_step(qs_env, rtp, rtp_control)
     386         1140 :             CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
     387         1140 :             rtp%energy_new = energy%total
     388         1140 :             IF (rtp%converged) EXIT
     389         1140 :             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          198 :       CPASSERT(ASSOCIATED(qs_env))
     500              : 
     501          360 :       SELECT CASE (rtp_control%initial_wfn)
     502              :       CASE (use_scf_wfn)
     503          162 :          qs_env%sim_time = 0.0_dp
     504          162 :          qs_env%sim_step = 0
     505          162 :          energy_consistency = .TRUE.
     506              :          !in the linear scaling case we need a correct kohn-sham matrix, which we cannot get with consistent energies
     507          162 :          IF (rtp_control%linear_scaling) energy_consistency = .FALSE.
     508              :          CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., &
     509          162 :                                           consistent_energies=energy_consistency)
     510          162 :          qs_env%run_rtp = .TRUE.
     511          162 :          ALLOCATE (qs_env%rtp)
     512          162 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     513          162 :          IF (dft_control%do_admm) THEN
     514            6 :             CALL hfx_admm_init(qs_env)
     515              :             CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     516            6 :                                 rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
     517              :          ELSE
     518              :             CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     519          156 :                                 rtp_control%linear_scaling)
     520              :          END IF
     521              : 
     522              :       CASE (use_restart_wfn, use_rt_restart)
     523           36 :          CALL qs_energies_init(qs_env, .FALSE.)
     524           36 :          IF (.NOT. rtp_control%linear_scaling .OR. rtp_control%initial_wfn == use_restart_wfn) THEN
     525           86 :             DO ispin = 1, SIZE(qs_env%mos)
     526           52 :                CALL get_mo_set(qs_env%mos(ispin), mo_coeff=mo_coeff, homo=homo)
     527           86 :                IF (.NOT. ASSOCIATED(mo_coeff)) THEN
     528              :                   CALL init_mo_set(qs_env%mos(ispin), &
     529              :                                    qs_env%mpools%ao_mo_fm_pools(ispin)%pool, &
     530           52 :                                    name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
     531              :                END IF
     532              :             END DO
     533           34 :             IF (dft_control%do_admm) CALL hfx_admm_init(qs_env)
     534              :          END IF
     535           36 :          ALLOCATE (qs_env%rtp)
     536           36 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     537              :          CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
     538           36 :                              rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
     539           36 :          CALL get_restart_wfn(qs_env)
     540           36 :          CPASSERT(ASSOCIATED(qs_env))
     541              : 
     542          234 :          qs_env%run_rtp = .TRUE.
     543              :       END SELECT
     544              : 
     545          198 :    END SUBROUTINE rt_initial_guess
     546              : 
     547              : ! **************************************************************************************************
     548              : !> \brief ...
     549              : !> \param qs_env ...
     550              : !> \param imag_p ...
     551              : !> \param imag_ks ...
     552              : !> \param imag_h ...
     553              : ! **************************************************************************************************
     554          198 :    SUBROUTINE rt_init_complex_quantities(qs_env, imag_p, imag_ks, imag_h)
     555              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     556              :       LOGICAL, INTENT(in)                                :: imag_p, imag_ks, imag_h
     557              : 
     558              :       TYPE(dft_control_type), POINTER                    :: dft_control
     559              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     560              :       TYPE(qs_rho_type), POINTER                         :: rho
     561              :       TYPE(rt_prop_type), POINTER                        :: rtp
     562              : 
     563          198 :       NULLIFY (ks_env, rho, dft_control)
     564              : 
     565              :       CALL get_qs_env(qs_env, &
     566              :                       dft_control=dft_control, &
     567              :                       ks_env=ks_env, &
     568              :                       rho=rho, &
     569          198 :                       rtp=rtp)
     570              : 
     571              :       ! rho
     572          198 :       CALL qs_rho_set(rho, complex_rho_ao=imag_p)
     573          198 :       IF (imag_p) CALL allocate_rho_ao_imag_from_real(rho, qs_env)
     574              : 
     575              :       ! ks
     576          198 :       CALL set_ks_env(ks_env, complex_ks=imag_ks)
     577          198 :       IF (imag_ks) THEN
     578           40 :          CALL qs_ks_allocate_basics(qs_env, is_complex=imag_ks)
     579           40 :          IF (.NOT. dft_control%rtp_control%fixed_ions) &
     580           22 :             CALL rtp_create_SinvH_imag(rtp, dft_control%nspins)
     581              :       END IF
     582              : 
     583              :       ! h
     584          198 :       IF (imag_h) CALL qs_matrix_h_allocate_imag_from_real(qs_env)
     585              : 
     586          198 :    END SUBROUTINE rt_init_complex_quantities
     587              : 
     588              : END MODULE rt_propagation
        

Generated by: LCOV version 2.0-1