LCOV - code coverage report
Current view: top level - src/motion - md_run.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 93.7 % 191 179
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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 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         3534 :    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         1767 :       my_rm_restart_info = .TRUE.
     135         1767 :       IF (PRESENT(rm_restart_info)) my_rm_restart_info = rm_restart_info
     136         1767 :       NULLIFY (md_env, para_env)
     137         1767 :       para_env => force_env%para_env
     138         1767 :       motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
     139         1767 :       md_section => section_vals_get_subs_vals(motion_section, "MD")
     140              : 
     141              :       ! Real call to MD driver - Low Level
     142         1767 :       ALLOCATE (md_env)
     143         1767 :       CALL md_env_create(md_env, md_section, para_env, force_env=force_env)
     144         1767 :       CALL set_md_env(md_env, averages=averages)
     145         1767 :       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         1739 :          CALL qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, mdctrl=mdctrl)
     150              :       END IF
     151         1767 :       CALL md_env_release(md_env)
     152         1767 :       DEALLOCATE (md_env)
     153              : 
     154              :       ! Clean restartable sections..
     155         1767 :       IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
     156         1767 :    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         8835 :    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         1767 :       CALL timeset(routineN, handle)
     210         1767 :       CPASSERT(ASSOCIATED(globenv))
     211         1767 :       CPASSERT(ASSOCIATED(force_env))
     212              : 
     213         1767 :       NULLIFY (particles, cell, simpar, itimes, used_time, subsys, &
     214         1767 :                md_ener, thermostats, barostat, reftraj, force_env_section, &
     215         1767 :                reftraj_section, work_section, atomic_kinds, &
     216         1767 :                local_particles, time, fe_env, free_energy_section, &
     217         1767 :                constraint_section, thermal_regions, virial, subsys_i)
     218         1767 :       logger => cp_get_default_logger()
     219         1767 :       para_env => force_env%para_env
     220              : 
     221         1767 :       global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL")
     222         1767 :       free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY")
     223         1767 :       constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT")
     224         1767 :       CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
     225              : 
     226         1767 :       CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id)
     227         1767 :       IF (run_type_id == ehrenfest) CALL set_md_env(md_env, ehrenfest_md=.TRUE.)
     228              : 
     229         1767 :       CALL create_simpar_type(simpar)
     230         1767 :       force_env_section => force_env%force_env_section
     231         1767 :       subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
     232         1767 :       CALL cp_add_iter_level(logger%iter_info, "MD")
     233         1767 :       CALL cp_iterate(logger%iter_info, iter_nr=0)
     234              :       ! Read MD section
     235         1767 :       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         1767 :                                                     "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.FALSE.)
     239              :       simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, &
     240         1767 :                                                          "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
     241              :       simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
     242         1767 :                                                         "LAGRANGE_MULTIPLIERS"), cp_p_file)
     243              : 
     244              :       ! Create the structure for the md energies
     245         7068 :       ALLOCATE (md_ener)
     246         1767 :       CALL create_md_ener(md_ener)
     247         1767 :       CALL set_md_env(md_env, md_ener=md_ener)
     248         1767 :       NULLIFY (md_ener)
     249              : 
     250              :       ! If requested setup Thermostats
     251              :       CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, &
     252         1767 :                               globenv, global_section)
     253              : 
     254              :       ! If requested setup Barostat
     255         1767 :       CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv)
     256              : 
     257              :       ! If requested setup different thermal regions
     258         1767 :       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         1767 :       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         1767 :       CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions)
     269              : 
     270         1767 :       CALL get_md_env(md_env, ehrenfest_md=ehrenfest_md)
     271              : 
     272              :       !If requested set up the REFTRAJ run
     273         1767 :       IF (simpar%ensemble == reftraj_ensemble .AND. ehrenfest_md) &
     274            0 :          CPABORT("Ehrenfest MD does not support reftraj ensemble ")
     275         1767 :       IF (simpar%ensemble == reftraj_ensemble) THEN
     276           38 :          reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ")
     277           38 :          ALLOCATE (reftraj)
     278           38 :          CALL create_reftraj(reftraj, reftraj_section, para_env)
     279           38 :          CALL set_md_env(md_env, reftraj=reftraj)
     280              :       END IF
     281              : 
     282              :       CALL force_env_get(force_env, subsys=subsys, cell=cell, &
     283         1767 :                          force_env_section=force_env_section)
     284         1767 :       CALL cp_subsys_get(subsys, virial=virial)
     285              : 
     286              :       ! Set V0 if needed
     287         1767 :       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         1767 :                                 l_val=write_binary_restart_file)
     294              :       CALL setup_velocities(force_env, simpar, globenv, md_env, md_section, constraint_section, &
     295         1767 :                             write_binary_restart_file)
     296              : 
     297              :       ! Setup Free Energy Calculation (if required)
     298         1767 :       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         1767 :                       force_env=force_env)
     302              : 
     303              :       ! Possibly initialize Wiener processes
     304              :       ![NB] Tested again within create_wiener_process.  Why??
     305         1767 :       IF (need_per_atom_wiener_process(md_env)) CALL create_wiener_process(md_env)
     306              : 
     307         1767 :       time_iter_start = m_walltime()
     308              : 
     309              :       CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
     310         1767 :                       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         1767 :       CALL set_meta_env(force_env%meta_env, time=time)
     314              : 
     315              :       ! Initialize the md_ener structure
     316         1767 :       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         1767 :           (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         1767 :       IF (simpar%ensemble /= reftraj_ensemble) THEN
     351         1729 :          CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
     352         1729 :          CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
     353         1729 :          CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
     354         1729 :          CALL cp_iterate(logger%iter_info, iter_nr=itimes)
     355         1729 :          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         1729 :          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         1657 :             CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
     372              :          END IF
     373              : 
     374         1729 :          IF (ASSOCIATED(force_env%qs_env)) THEN
     375          596 :             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         1729 :          IF (ASSOCIATED(force_env%meta_env)) THEN
     381              :             ! Setup stuff for plumed if needed
     382          146 :             IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
     383            2 :                CALL metadyn_initialise_plumed(force_env, simpar, itimes)
     384              :             ELSE
     385          144 :                IF (force_env%meta_env%langevin) THEN
     386            4 :                   CALL create_wiener_process_cv(force_env%meta_env)
     387              :                END IF
     388          144 :                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          144 :                CALL metadyn_forces(force_env)
     408          144 :                CALL metadyn_write_colvar(force_env)
     409              :             END IF
     410              :          END IF
     411              : 
     412         1729 :          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         1729 :          CALL force_env_get(force_env, subsys=subsys)
     418              : 
     419              :          CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     420         1729 :                             particles=particles, virial=virial)
     421              : 
     422              :          CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
     423         1729 :                               virial, force_env%para_env)
     424              : 
     425         1729 :          CALL md_energy(md_env, md_ener)
     426         1729 :          CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories
     427         1729 :          md_stride = 1
     428              :       ELSE
     429           38 :          CALL get_md_env(md_env, reftraj=reftraj)
     430           38 :          CALL initialize_reftraj(reftraj, reftraj_section, md_env)
     431           38 :          itimes = reftraj%info%first_snapshot - 1
     432           38 :          md_stride = reftraj%info%stride
     433           38 :          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         1767 :                                         constraint_section, "CONSTRAINT_INFO")
     442              :       CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
     443         1767 :                                         constraint_section, "LAGRANGE_MULTIPLIERS")
     444              : 
     445              : ! if we need the initial kinetic energy for Hybrid Monte Carlo
     446         1767 :       IF (PRESENT(hmc_e_initial)) hmc_e_initial = md_ener%ekin
     447              : 
     448         1767 :       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        42367 :       DO istep = 1, simpar%nsteps, md_stride
     453              :          ! Increase counters
     454        40631 :          itimes = itimes + 1
     455        40631 :          time = time + simpar%dt
     456              :          !needed when electric field fields are applied
     457        40631 :          IF (ASSOCIATED(force_env%qs_env)) THEN
     458         3198 :             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        40631 :          IF (ehrenfest_md) force_env%qs_env%rtp%istep = istep
     463              : 
     464        40631 :          IF (.NOT. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) THEN
     465        40631 :             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        40631 :                                                        extension=".shakeLog", log_filename=.FALSE.)
     473              :          simpar%lagrange_multipliers = cp_print_key_unit_nr( &
     474              :                                        logger, constraint_section, &
     475        40631 :                                        "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
     476              :          simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
     477        40631 :                                                            "LAGRANGE_MULTIPLIERS"), cp_p_file)
     478              : 
     479              :          ! Velocity Verlet Integrator
     480        40631 :          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        40631 :                                            constraint_section, "CONSTRAINT_INFO")
     485              :          CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
     486        40631 :                                            constraint_section, "LAGRANGE_MULTIPLIERS")
     487              : 
     488              :          ! Free Energy calculation
     489        40631 :          CALL free_energy_evaluate(md_env, should_stop, free_energy_section)
     490              : 
     491        40631 :          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        40631 :          CALL external_control(should_stop, "MD", globenv=globenv)
     503              : 
     504              :          !check if upper bound of total steps has been reached
     505        40631 :          IF (.NOT. (istep == simpar%nsteps) .AND. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) should_stop = .TRUE.
     506        40631 :          IF (itimes >= simpar%max_steps) should_stop = .TRUE.
     507              : 
     508              :          ! call external hook e.g. from global optimization
     509        40631 :          IF (PRESENT(mdctrl)) &
     510         2689 :             CALL mdctrl_callback(mdctrl, md_env, should_stop)
     511              : 
     512        40631 :          IF (should_stop) THEN
     513           31 :             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           31 :             IF (run_type_id == ehrenfest) THEN
     517            2 :                CALL md_output(md_env, md_section, force_env%root_section, .FALSE.)
     518              :             ELSE
     519           29 :                CALL md_output(md_env, md_section, force_env%root_section, should_stop)
     520              :             END IF
     521              :             EXIT
     522              :          END IF
     523              : 
     524        40600 :          IF (simpar%ensemble /= reftraj_ensemble) THEN
     525        40312 :             CALL md_energy(md_env, md_ener)
     526        40312 :             CALL temperature_control(simpar, md_env, md_ener, force_env, logger)
     527        40312 :             CALL comvel_control(md_ener, force_env, md_section, logger)
     528        40312 :             CALL angvel_control(md_ener, force_env, md_section, logger)
     529              :          ELSE
     530          288 :             CALL md_ener_reftraj(md_env, md_ener)
     531              :          END IF
     532              : 
     533        40600 :          time_iter_stop = m_walltime()
     534        40600 :          used_time = time_iter_stop - time_iter_start
     535        40600 :          time_iter_start = time_iter_stop
     536              : 
     537        40600 :          CALL md_output(md_env, md_section, force_env%root_section, should_stop)
     538       123598 :          IF (simpar%ensemble == reftraj_ensemble) THEN
     539          288 :             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         1767 :       IF (PRESENT(hmc_e_final)) hmc_e_final = md_ener%ekin
     545              : 
     546              :       ! Remove the iteration level
     547         1767 :       CALL cp_rm_iter_level(logger%iter_info, "MD")
     548              : 
     549              :       ! Clean up PLUMED
     550         1767 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     551          150 :          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         1767 :       CALL release_simpar_type(simpar)
     558         1767 :       CALL timestop(handle)
     559              : 
     560         1767 :    END SUBROUTINE qs_mol_dyn_low
     561              : 
     562              : END MODULE md_run
        

Generated by: LCOV version 2.0-1