LCOV - code coverage report
Current view: top level - src/motion - md_run.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 179 191 93.7 %
Date: 2024-03-28 07:31:50 Functions: 2 2 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 Perform a molecular dynamics (MD) run using QUICKSTEP
      10             : !> \par History
      11             : !>   - Added support for Langevin regions (2014/02/05, LT)
      12             : !> \author Matthias Krack (07.11.2002)
      13             : ! **************************************************************************************************
      14             : MODULE md_run
      15             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      16             :    USE averages_types,                  ONLY: average_quantities_type
      17             :    USE barostat_types,                  ONLY: barostat_type,&
      18             :                                               create_barostat_type
      19             :    USE cell_types,                      ONLY: cell_type
      20             :    USE cp_external_control,             ONLY: external_control
      21             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      22             :                                               cp_logger_type
      23             :    USE cp_output_handling,              ONLY: cp_add_iter_level,&
      24             :                                               cp_iterate,&
      25             :                                               cp_p_file,&
      26             :                                               cp_print_key_finished_output,&
      27             :                                               cp_print_key_should_output,&
      28             :                                               cp_print_key_unit_nr,&
      29             :                                               cp_rm_iter_level
      30             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      31             :                                               cp_subsys_type
      32             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      33             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      34             :    USE force_env_types,                 ONLY: force_env_get,&
      35             :                                               force_env_type
      36             :    USE free_energy_methods,             ONLY: free_energy_evaluate
      37             :    USE free_energy_types,               ONLY: fe_env_create,&
      38             :                                               free_energy_type
      39             :    USE global_types,                    ONLY: global_environment_type
      40             :    USE input_constants,                 ONLY: &
      41             :         ehrenfest, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
      42             :         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
      43             :         npt_ia_ensemble, reftraj_ensemble
      44             :    USE input_cp2k_check,                ONLY: remove_restart_info
      45             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      46             :                                               section_vals_remove_values,&
      47             :                                               section_vals_type,&
      48             :                                               section_vals_val_get
      49             :    USE kinds,                           ONLY: default_string_length,&
      50             :                                               dp
      51             :    USE machine,                         ONLY: m_walltime
      52             :    USE md_ener_types,                   ONLY: create_md_ener,&
      53             :                                               md_ener_type
      54             :    USE md_energies,                     ONLY: initialize_md_ener,&
      55             :                                               md_ener_reftraj,&
      56             :                                               md_energy,&
      57             :                                               md_write_output
      58             :    USE md_environment_types,            ONLY: get_md_env,&
      59             :                                               md_env_create,&
      60             :                                               md_env_release,&
      61             :                                               md_environment_type,&
      62             :                                               need_per_atom_wiener_process,&
      63             :                                               set_md_env
      64             :    USE md_util,                         ONLY: md_output
      65             :    USE md_vel_utils,                    ONLY: angvel_control,&
      66             :                                               comvel_control,&
      67             :                                               setup_velocities,&
      68             :                                               temperature_control
      69             :    USE mdctrl_methods,                  ONLY: mdctrl_callback
      70             :    USE mdctrl_types,                    ONLY: mdctrl_type
      71             :    USE message_passing,                 ONLY: mp_para_env_type
      72             :    USE metadynamics,                    ONLY: metadyn_finalise_plumed,&
      73             :                                               metadyn_forces,&
      74             :                                               metadyn_initialise_plumed,&
      75             :                                               metadyn_write_colvar
      76             :    USE metadynamics_types,              ONLY: set_meta_env
      77             :    USE particle_list_types,             ONLY: particle_list_type
      78             :    USE qs_environment_methods,          ONLY: qs_env_time_update
      79             :    USE reftraj_types,                   ONLY: create_reftraj,&
      80             :                                               reftraj_type
      81             :    USE reftraj_util,                    ONLY: initialize_reftraj,&
      82             :                                               write_output_reftraj
      83             :    USE rt_propagation,                  ONLY: rt_prop_setup
      84             :    USE simpar_methods,                  ONLY: read_md_section
      85             :    USE simpar_types,                    ONLY: create_simpar_type,&
      86             :                                               release_simpar_type,&
      87             :                                               simpar_type
      88             :    USE thermal_region_types,            ONLY: thermal_regions_type
      89             :    USE thermal_region_utils,            ONLY: create_thermal_regions,&
      90             :                                               print_thermal_regions_langevin
      91             :    USE thermostat_methods,              ONLY: create_thermostats
      92             :    USE thermostat_types,                ONLY: thermostats_type
      93             :    USE velocity_verlet_control,         ONLY: velocity_verlet
      94             :    USE virial_methods,                  ONLY: virial_evaluate
      95             :    USE virial_types,                    ONLY: virial_type
      96             :    USE wiener_process,                  ONLY: create_wiener_process,&
      97             :                                               create_wiener_process_cv
      98             : #include "../base/base_uses.f90"
      99             : 
     100             :    IMPLICIT NONE
     101             : 
     102             :    PRIVATE
     103             : 
     104             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_run'
     105             : 
     106             :    PUBLIC :: qs_mol_dyn
     107             : 
     108             : CONTAINS
     109             : 
     110             : ! **************************************************************************************************
     111             : !> \brief Main driver module for Molecular Dynamics
     112             : !> \param force_env ...
     113             : !> \param globenv ...
     114             : !> \param averages ...
     115             : !> \param rm_restart_info ...
     116             : !> \param hmc_e_initial ...
     117             : !> \param hmc_e_final ...
     118             : !> \param mdctrl ...
     119             : ! **************************************************************************************************
     120        1790 :    SUBROUTINE qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
     121             : 
     122             :       TYPE(force_env_type), POINTER                      :: force_env
     123             :       TYPE(global_environment_type), POINTER             :: globenv
     124             :       TYPE(average_quantities_type), OPTIONAL, POINTER   :: averages
     125             :       LOGICAL, INTENT(IN), OPTIONAL                      :: rm_restart_info
     126             :       REAL(KIND=dp), OPTIONAL                            :: hmc_e_initial, hmc_e_final
     127             :       TYPE(mdctrl_type), OPTIONAL, POINTER               :: mdctrl
     128             : 
     129             :       LOGICAL                                            :: my_rm_restart_info
     130             :       TYPE(md_environment_type), POINTER                 :: md_env
     131             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     132             :       TYPE(section_vals_type), POINTER                   :: md_section, motion_section
     133             : 
     134        1790 :       my_rm_restart_info = .TRUE.
     135        1790 :       IF (PRESENT(rm_restart_info)) my_rm_restart_info = rm_restart_info
     136        1790 :       NULLIFY (md_env, para_env)
     137        1790 :       para_env => force_env%para_env
     138        1790 :       motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
     139        1790 :       md_section => section_vals_get_subs_vals(motion_section, "MD")
     140             : 
     141             :       ! Real call to MD driver - Low Level
     142        1790 :       ALLOCATE (md_env)
     143        1790 :       CALL md_env_create(md_env, md_section, para_env, force_env=force_env)
     144        1790 :       CALL set_md_env(md_env, averages=averages)
     145        1790 :       IF (PRESENT(hmc_e_initial) .AND. PRESENT(hmc_e_final)) THEN
     146             :          CALL qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, &
     147          28 :                              hmc_e_initial=hmc_e_initial, hmc_e_final=hmc_e_final)
     148             :       ELSE
     149        1762 :          CALL qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, mdctrl=mdctrl)
     150             :       END IF
     151        1790 :       CALL md_env_release(md_env)
     152        1790 :       DEALLOCATE (md_env)
     153             : 
     154             :       ! Clean restartable sections..
     155        1790 :       IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
     156        1790 :    END SUBROUTINE qs_mol_dyn
     157             : 
     158             : ! **************************************************************************************************
     159             : !> \brief Purpose: Driver routine for MD run using QUICKSTEP.
     160             : !> \param md_env ...
     161             : !> \param md_section ...
     162             : !> \param motion_section ...
     163             : !> \param force_env ...
     164             : !> \param globenv ...
     165             : !> \param hmc_e_initial ...
     166             : !> \param hmc_e_final ...
     167             : !> \param mdctrl ...
     168             : !> \par History
     169             : !>   - Cleaning (09.2007) Teodoro Laino [tlaino] - University of Zurich
     170             : !>   - Added lines to print out langevin regions (2014/02/04, LT)
     171             : !> \author Creation (07.11.2002,MK)
     172             : ! **************************************************************************************************
     173        7160 :    SUBROUTINE qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, hmc_e_initial, hmc_e_final, mdctrl)
     174             : 
     175             :       TYPE(md_environment_type), POINTER                 :: md_env
     176             :       TYPE(section_vals_type), POINTER                   :: md_section, motion_section
     177             :       TYPE(force_env_type), POINTER                      :: force_env
     178             :       TYPE(global_environment_type), POINTER             :: globenv
     179             :       REAL(KIND=dp), OPTIONAL                            :: hmc_e_initial, hmc_e_final
     180             :       TYPE(mdctrl_type), OPTIONAL, POINTER               :: mdctrl
     181             : 
     182             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_mol_dyn_low'
     183             : 
     184             :       CHARACTER(LEN=default_string_length)               :: my_act, my_pos
     185             :       INTEGER                                            :: handle, i, istep, md_stride, run_type_id
     186             :       INTEGER, POINTER                                   :: itimes
     187             :       LOGICAL                                            :: check, ehrenfest_md, save_mem, &
     188             :                                                             should_stop, write_binary_restart_file
     189             :       REAL(KIND=dp)                                      :: dummy, time_iter_start, time_iter_stop
     190             :       REAL(KIND=dp), POINTER                             :: constant, time, used_time
     191             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     192             :       TYPE(barostat_type), POINTER                       :: barostat
     193             :       TYPE(cell_type), POINTER                           :: cell
     194             :       TYPE(cp_logger_type), POINTER                      :: logger
     195             :       TYPE(cp_subsys_type), POINTER                      :: subsys, subsys_i
     196             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     197             :       TYPE(free_energy_type), POINTER                    :: fe_env
     198             :       TYPE(md_ener_type), POINTER                        :: md_ener
     199             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     200             :       TYPE(particle_list_type), POINTER                  :: particles
     201             :       TYPE(reftraj_type), POINTER                        :: reftraj
     202             :       TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, &
     203             :          free_energy_section, global_section, reftraj_section, subsys_section, work_section
     204             :       TYPE(simpar_type), POINTER                         :: simpar
     205             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     206             :       TYPE(thermostats_type), POINTER                    :: thermostats
     207             :       TYPE(virial_type), POINTER                         :: virial
     208             : 
     209        1790 :       CALL timeset(routineN, handle)
     210        1790 :       CPASSERT(ASSOCIATED(globenv))
     211        1790 :       CPASSERT(ASSOCIATED(force_env))
     212             : 
     213        1790 :       NULLIFY (particles, cell, simpar, itimes, used_time, subsys, &
     214        1790 :                md_ener, thermostats, barostat, reftraj, force_env_section, &
     215        1790 :                reftraj_section, work_section, atomic_kinds, &
     216        1790 :                local_particles, time, fe_env, free_energy_section, &
     217        1790 :                constraint_section, thermal_regions, virial, subsys_i)
     218        1790 :       logger => cp_get_default_logger()
     219        1790 :       para_env => force_env%para_env
     220             : 
     221        1790 :       global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL")
     222        1790 :       free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY")
     223        1790 :       constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT")
     224        1790 :       CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
     225             : 
     226        1790 :       CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id)
     227        1790 :       IF (run_type_id == ehrenfest) CALL set_md_env(md_env, ehrenfest_md=.TRUE.)
     228             : 
     229        1790 :       CALL create_simpar_type(simpar)
     230        1790 :       force_env_section => force_env%force_env_section
     231        1790 :       subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
     232        1790 :       CALL cp_add_iter_level(logger%iter_info, "MD")
     233        1790 :       CALL cp_iterate(logger%iter_info, iter_nr=0)
     234             :       ! Read MD section
     235        1790 :       CALL read_md_section(simpar, motion_section, md_section)
     236             :       ! Setup print_keys
     237             :       simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, &
     238        1790 :                                                     "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.FALSE.)
     239             :       simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, &
     240        1790 :                                                          "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
     241             :       simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
     242        1790 :                                                         "LAGRANGE_MULTIPLIERS"), cp_p_file)
     243             : 
     244             :       ! Create the structure for the md energies
     245        7160 :       ALLOCATE (md_ener)
     246        1790 :       CALL create_md_ener(md_ener)
     247        1790 :       CALL set_md_env(md_env, md_ener=md_ener)
     248        1790 :       NULLIFY (md_ener)
     249             : 
     250             :       ! If requested setup Thermostats
     251             :       CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, &
     252        1790 :                               globenv, global_section)
     253             : 
     254             :       ! If requested setup Barostat
     255        1790 :       CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv)
     256             : 
     257             :       ! If requested setup different thermal regions
     258        1790 :       CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env)
     259             : 
     260             :       ! If doing langevin_ensemble, then print out langevin_regions information upon request
     261        1790 :       IF (simpar%ensemble == langevin_ensemble) THEN
     262          42 :          my_pos = "REWIND"
     263          42 :          my_act = "WRITE"
     264             :          CALL print_thermal_regions_langevin(thermal_regions, simpar, &
     265          42 :                                              pos=my_pos, act=my_act)
     266             :       END IF
     267             : 
     268        1790 :       CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions)
     269             : 
     270        1790 :       CALL get_md_env(md_env, ehrenfest_md=ehrenfest_md)
     271             : 
     272             :       !If requested set up the REFTRAJ run
     273        1790 :       IF (simpar%ensemble == reftraj_ensemble .AND. ehrenfest_md) &
     274           0 :          CPABORT("Ehrenfest MD does not support reftraj ensemble ")
     275        1790 :       IF (simpar%ensemble == reftraj_ensemble) THEN
     276          36 :          reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ")
     277          36 :          ALLOCATE (reftraj)
     278          36 :          CALL create_reftraj(reftraj, reftraj_section, para_env)
     279          36 :          CALL set_md_env(md_env, reftraj=reftraj)
     280             :       END IF
     281             : 
     282             :       CALL force_env_get(force_env, subsys=subsys, cell=cell, &
     283        1790 :                          force_env_section=force_env_section)
     284        1790 :       CALL cp_subsys_get(subsys, virial=virial)
     285             : 
     286             :       ! Set V0 if needed
     287        1790 :       IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
     288           6 :          IF (simpar%v0 == 0._dp) simpar%v0 = cell%deth
     289             :       END IF
     290             : 
     291             :       ! Initialize velocities possibly applying constraints at the zeroth MD step
     292             :       CALL section_vals_val_get(motion_section, "PRINT%RESTART%SPLIT_RESTART_FILE", &
     293        1790 :                                 l_val=write_binary_restart_file)
     294             :       CALL setup_velocities(force_env, simpar, globenv, md_env, md_section, constraint_section, &
     295        1790 :                             write_binary_restart_file)
     296             : 
     297             :       ! Setup Free Energy Calculation (if required)
     298        1790 :       CALL fe_env_create(fe_env, free_energy_section)
     299             : 
     300             :       CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, &
     301        1790 :                       force_env=force_env)
     302             : 
     303             :       ! Possibly initialize Wiener processes
     304             :       ![NB] Tested again within create_wiener_process.  Why??
     305        1790 :       IF (need_per_atom_wiener_process(md_env)) CALL create_wiener_process(md_env)
     306             : 
     307        1790 :       time_iter_start = m_walltime()
     308             : 
     309             :       CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
     310        1790 :                       md_ener=md_ener, t=time, used_time=used_time)
     311             : 
     312             :       ! Attach the time counter of the meta_env to the one of the MD
     313        1790 :       CALL set_meta_env(force_env%meta_env, time=time)
     314             : 
     315             :       ! Initialize the md_ener structure
     316        1790 :       CALL initialize_md_ener(md_ener, force_env, simpar)
     317             : 
     318             :       ! Check for ensembles requiring the stress tensor - takes into account the possibility for
     319             :       ! multiple force_evals
     320             :       IF ((simpar%ensemble == npt_i_ensemble) .OR. &
     321             :           (simpar%ensemble == npt_ia_ensemble) .OR. &
     322             :           (simpar%ensemble == npt_f_ensemble) .OR. &
     323             :           (simpar%ensemble == npe_f_ensemble) .OR. &
     324             :           (simpar%ensemble == npe_i_ensemble) .OR. &
     325        1790 :           (simpar%ensemble == nph_uniaxial_ensemble) .OR. &
     326             :           (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
     327         172 :          check = virial%pv_availability
     328         172 :          IF (.NOT. check) &
     329             :             CALL cp_abort(__LOCATION__, &
     330             :                           "Virial evaluation not requested for this run in the input file!"// &
     331             :                           " You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR."// &
     332           0 :                           " Be sure the method you are using can compute the virial!")
     333         172 :          IF (ASSOCIATED(force_env%sub_force_env)) THEN
     334          26 :             DO i = 1, SIZE(force_env%sub_force_env)
     335          26 :                IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN
     336          10 :                   CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i)
     337          10 :                   CALL cp_subsys_get(subsys_i, virial=virial)
     338          10 :                   check = check .AND. virial%pv_availability
     339             :                END IF
     340             :             END DO
     341             :          END IF
     342         172 :          IF (.NOT. check) &
     343             :             CALL cp_abort(__LOCATION__, &
     344             :                           "Virial evaluation not requested for all the force_eval sections present in"// &
     345             :                           " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR"// &
     346           0 :                           " in each force_eval section. Be sure the method you are using can compute the virial!")
     347             :       END IF
     348             : 
     349             :       ! Computing Forces at zero MD step
     350        1790 :       IF (simpar%ensemble /= reftraj_ensemble) THEN
     351        1754 :          CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
     352        1754 :          CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
     353        1754 :          CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
     354        1754 :          CALL cp_iterate(logger%iter_info, iter_nr=itimes)
     355        1754 :          IF (save_mem) THEN
     356           2 :             work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
     357           2 :             CALL section_vals_remove_values(work_section)
     358           2 :             work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
     359           2 :             CALL section_vals_remove_values(work_section)
     360           2 :             work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
     361           2 :             CALL section_vals_remove_values(work_section)
     362             :          END IF
     363             : 
     364        1754 :          IF (ehrenfest_md) THEN
     365          72 :             CALL rt_prop_setup(force_env)
     366          72 :             force_env%qs_env%rtp%dt = simpar%dt
     367             :          ELSE
     368             :             ![NB] Lets let all methods, even ones without consistent energies, succeed here.
     369             :             !     They'll fail in actual integrator if needed
     370             :             ! consistent_energies=.FALSE. by default
     371        1682 :             CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
     372             :          END IF
     373             : 
     374        1754 :          IF (ASSOCIATED(force_env%qs_env)) THEN
     375         592 :             CALL qs_env_time_update(force_env%qs_env, time, itimes)
     376             : !deb        force_env%qs_env%sim_time = time
     377             : !deb        force_env%qs_env%sim_step = itimes
     378             :          END IF
     379             :          ! Warm-up engines for metadynamics
     380        1754 :          IF (ASSOCIATED(force_env%meta_env)) THEN
     381             :             ! Setup stuff for plumed if needed
     382         144 :             IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
     383           2 :                CALL metadyn_initialise_plumed(force_env, simpar, itimes)
     384             :             ELSE
     385         142 :                IF (force_env%meta_env%langevin) THEN
     386           4 :                   CALL create_wiener_process_cv(force_env%meta_env)
     387             :                END IF
     388         142 :                IF (force_env%meta_env%well_tempered) THEN
     389           2 :                   force_env%meta_env%wttemperature = simpar%temp_ext
     390           2 :                   IF (force_env%meta_env%wtgamma > EPSILON(1._dp)) THEN
     391           0 :                      dummy = force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma - 1._dp)
     392           0 :                      IF (force_env%meta_env%delta_t > EPSILON(1._dp)) THEN
     393           0 :                         check = ABS(force_env%meta_env%delta_t - dummy) < 1.E+3_dp*EPSILON(1._dp)
     394           0 :                         IF (.NOT. check) &
     395             :                            CALL cp_abort(__LOCATION__, &
     396             :                                          "Inconsistency between DELTA_T and WTGAMMA (both specified):"// &
     397           0 :                                          " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE")
     398             :                      ELSE
     399           0 :                         force_env%meta_env%delta_t = dummy
     400             :                      END IF
     401             :                   ELSE
     402             :                      force_env%meta_env%wtgamma = 1._dp &
     403           2 :                                                   + force_env%meta_env%delta_t/force_env%meta_env%wttemperature
     404             :                   END IF
     405           2 :                   force_env%meta_env%invdt = 1._dp/force_env%meta_env%delta_t
     406             :                END IF
     407         142 :                CALL metadyn_forces(force_env)
     408         142 :                CALL metadyn_write_colvar(force_env)
     409             :             END IF
     410             :          END IF
     411             : 
     412        1754 :          IF (simpar%do_respa) THEN
     413             :             CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
     414           6 :                                              calc_force=.TRUE.)
     415             :          END IF
     416             : 
     417        1754 :          CALL force_env_get(force_env, subsys=subsys)
     418             : 
     419             :          CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     420        1754 :                             particles=particles, virial=virial)
     421             : 
     422             :          CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
     423        1754 :                               virial, force_env%para_env)
     424             : 
     425        1754 :          CALL md_energy(md_env, md_ener)
     426        1754 :          CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories
     427        1754 :          md_stride = 1
     428             :       ELSE
     429          36 :          CALL get_md_env(md_env, reftraj=reftraj)
     430          36 :          CALL initialize_reftraj(reftraj, reftraj_section, md_env)
     431          36 :          itimes = reftraj%info%first_snapshot - 1
     432          36 :          md_stride = reftraj%info%stride
     433          36 :          IF (ASSOCIATED(force_env%meta_env)) THEN
     434           4 :             IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
     435           0 :                CALL metadyn_initialise_plumed(force_env, simpar, itimes)
     436             :             END IF
     437             :          END IF
     438             :       END IF
     439             : 
     440             :       CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
     441        1790 :                                         constraint_section, "CONSTRAINT_INFO")
     442             :       CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
     443        1790 :                                         constraint_section, "LAGRANGE_MULTIPLIERS")
     444             : 
     445             : ! if we need the initial kinetic energy for Hybrid Monte Carlo
     446        1790 :       IF (PRESENT(hmc_e_initial)) hmc_e_initial = md_ener%ekin
     447             : 
     448        1790 :       IF (itimes >= simpar%max_steps) CALL cp_abort(__LOCATION__, &
     449           0 :                                                     "maximum step number smaller than initial step value")
     450             : 
     451             :       ! Real MD Loop
     452       44855 :       DO istep = 1, simpar%nsteps, md_stride
     453             :          ! Increase counters
     454       43131 :          itimes = itimes + 1
     455       43131 :          time = time + simpar%dt
     456             :          !needed when electric field fields are applied
     457       43131 :          IF (ASSOCIATED(force_env%qs_env)) THEN
     458        3066 :             CALL qs_env_time_update(force_env%qs_env, time, itimes)
     459             : !deb        force_env%qs_env%sim_time = time
     460             : !deb        force_env%qs_env%sim_step = itimes
     461             :          END IF
     462       43131 :          IF (ehrenfest_md) force_env%qs_env%rtp%istep = istep
     463             : 
     464       43131 :          IF (.NOT. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) THEN
     465       43131 :             CALL cp_iterate(logger%iter_info, last=(istep == simpar%nsteps), iter_nr=itimes)
     466             :          ELSE
     467           0 :             CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
     468             :          END IF
     469             : 
     470             :          ! Open possible Shake output units
     471             :          simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", &
     472       43131 :                                                        extension=".shakeLog", log_filename=.FALSE.)
     473             :          simpar%lagrange_multipliers = cp_print_key_unit_nr( &
     474             :                                        logger, constraint_section, &
     475       43131 :                                        "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
     476             :          simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
     477       43131 :                                                            "LAGRANGE_MULTIPLIERS"), cp_p_file)
     478             : 
     479             :          ! Velocity Verlet Integrator
     480       43131 :          CALL velocity_verlet(md_env, globenv)
     481             : 
     482             :          ! Close Shake output if requested...
     483             :          CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
     484       43131 :                                            constraint_section, "CONSTRAINT_INFO")
     485             :          CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
     486       43131 :                                            constraint_section, "LAGRANGE_MULTIPLIERS")
     487             : 
     488             :          ! Free Energy calculation
     489       43131 :          CALL free_energy_evaluate(md_env, should_stop, free_energy_section)
     490             : 
     491       43131 :          IF (should_stop) EXIT
     492             : 
     493             :          ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit
     494             :          ! Default:
     495             :          ! IF so we don't overwrite the restart or append to the trajectory
     496             :          ! because the execution could in principle stop inside the SCF where energy
     497             :          ! and forces are not converged.
     498             :          ! But:
     499             :          ! You can force to print the last step (for example if the method used
     500             :          ! to compute energy and forces is not SCF based) activating the print_key
     501             :          ! MOTION%MD%PRINT%FORCE_LAST.
     502       43131 :          CALL external_control(should_stop, "MD", globenv=globenv)
     503             : 
     504             :          !check if upper bound of total steps has been reached
     505       43131 :          IF (.NOT. (istep == simpar%nsteps) .AND. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) should_stop = .TRUE.
     506       43131 :          IF (itimes >= simpar%max_steps) should_stop = .TRUE.
     507             : 
     508             :          ! call external hook e.g. from global optimization
     509       43131 :          IF (PRESENT(mdctrl)) &
     510        5625 :             CALL mdctrl_callback(mdctrl, md_env, should_stop)
     511             : 
     512       43131 :          IF (should_stop) THEN
     513          66 :             CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
     514             :             !In Ehrenfest molecular dynamics the external control is only checked after a converged propagation
     515             :             !The restart needs to be written in order to be consistent with the mos/density matrix for the restart
     516          66 :             IF (run_type_id == ehrenfest) THEN
     517           2 :                CALL md_output(md_env, md_section, force_env%root_section, .FALSE.)
     518             :             ELSE
     519          64 :                CALL md_output(md_env, md_section, force_env%root_section, should_stop)
     520             :             END IF
     521             :             EXIT
     522             :          END IF
     523             : 
     524       43065 :          IF (simpar%ensemble /= reftraj_ensemble) THEN
     525       42779 :             CALL md_energy(md_env, md_ener)
     526       42779 :             CALL temperature_control(simpar, md_env, md_ener, force_env, logger)
     527       42779 :             CALL comvel_control(md_ener, force_env, md_section, logger)
     528       42779 :             CALL angvel_control(md_ener, force_env, md_section, logger)
     529             :          ELSE
     530         286 :             CALL md_ener_reftraj(md_env, md_ener)
     531             :          END IF
     532             : 
     533       43065 :          time_iter_stop = m_walltime()
     534       43065 :          used_time = time_iter_stop - time_iter_start
     535       43065 :          time_iter_start = time_iter_stop
     536             : 
     537       43065 :          CALL md_output(md_env, md_section, force_env%root_section, should_stop)
     538      131051 :          IF (simpar%ensemble == reftraj_ensemble) THEN
     539         286 :             CALL write_output_reftraj(md_env)
     540             :          END IF
     541             :       END DO
     542             : 
     543             : ! if we need the final kinetic energy for Hybrid Monte Carlo
     544        1790 :       IF (PRESENT(hmc_e_final)) hmc_e_final = md_ener%ekin
     545             : 
     546             :       ! Remove the iteration level
     547        1790 :       CALL cp_rm_iter_level(logger%iter_info, "MD")
     548             : 
     549             :       ! Clean up PLUMED
     550        1790 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     551         148 :          IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
     552           2 :             CALL metadyn_finalise_plumed()
     553             :          END IF
     554             :       END IF
     555             : 
     556             :       ! Deallocate Thermostats and Barostats
     557        1790 :       CALL release_simpar_type(simpar)
     558        1790 :       CALL timestop(handle)
     559             : 
     560        1790 :    END SUBROUTINE qs_mol_dyn_low
     561             : 
     562             : END MODULE md_run

Generated by: LCOV version 1.15