LCOV - code coverage report
Current view: top level - src/motion - integrator.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 881 1037 85.0 %
Date: 2024-04-20 06:29:22 Functions: 10 11 90.9 %

          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 Provides integrator routines (velocity verlet) for all the
      10             : !>      ensemble types
      11             : !> \par History
      12             : !>      JGH (15-Mar-2001) : Pass logical for box change to force routine
      13             : !>      Harald Forbert (Apr-2001): added path integral routine nvt_pimd
      14             : !>      CJM (15-Apr-2001) : added coef integrators and energy routines
      15             : !>      Joost VandeVondele (Juli-2003): simple version of isokinetic ensemble
      16             : !>      Teodoro Laino [tlaino] 10.2007 - University of Zurich: Generalization to
      17             : !>                                       different kind of thermostats
      18             : !>      Teodoro Laino [tlaino] 11.2007 - Metadynamics: now part of the MD modules
      19             : !>      Marcella Iannuzzi      02.2008 - Collecting common code (VV and creation of
      20             : !>                                       a temporary type)
      21             : !>      Teodoro Laino [tlaino] 02.2008 - Splitting integrator module and keeping in
      22             : !>                                       integrator only the INTEGRATORS
      23             : !>      Lianheng Tong [LT]     12.2013 - Added regions to Langevin MD
      24             : !> \author CJM
      25             : ! **************************************************************************************************
      26             : MODULE integrator
      27             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      28             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      29             :                                               get_atomic_kind,&
      30             :                                               get_atomic_kind_set
      31             :    USE barostat_types,                  ONLY: barostat_type
      32             :    USE cell_methods,                    ONLY: init_cell,&
      33             :                                               set_cell_param
      34             :    USE cell_types,                      ONLY: cell_type,&
      35             :                                               parse_cell_line
      36             :    USE constraint,                      ONLY: rattle_control,&
      37             :                                               shake_control,&
      38             :                                               shake_roll_control,&
      39             :                                               shake_update_targets
      40             :    USE constraint_fxd,                  ONLY: create_local_fixd_list,&
      41             :                                               fix_atom_control,&
      42             :                                               release_local_fixd_list
      43             :    USE constraint_util,                 ONLY: getold,&
      44             :                                               pv_constraint
      45             :    USE cp_control_types,                ONLY: dft_control_type
      46             :    USE cp_log_handling,                 ONLY: cp_to_string
      47             :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      48             :                                               parser_read_line
      49             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      50             :                                               cp_subsys_type
      51             :    USE cp_units,                        ONLY: cp_unit_to_cp2k
      52             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      53             :    USE eigenvalueproblems,              ONLY: diagonalise
      54             :    USE extended_system_dynamics,        ONLY: shell_scale_comv
      55             :    USE extended_system_types,           ONLY: npt_info_type
      56             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      57             :    USE force_env_types,                 ONLY: force_env_get,&
      58             :                                               force_env_type
      59             :    USE global_types,                    ONLY: global_environment_type
      60             :    USE input_constants,                 ONLY: ehrenfest,&
      61             :                                               npe_f_ensemble,&
      62             :                                               npe_i_ensemble,&
      63             :                                               npt_ia_ensemble
      64             :    USE integrator_utils,                ONLY: &
      65             :         allocate_old, allocate_tmp, damp_v, damp_veps, deallocate_old, get_s_ds, &
      66             :         old_variables_type, rattle_roll_setup, set, tmp_variables_type, update_dealloc_tmp, &
      67             :         update_pv, update_veps, variable_timestep, vv_first, vv_second
      68             :    USE kinds,                           ONLY: dp,&
      69             :                                               max_line_length
      70             :    USE md_environment_types,            ONLY: get_md_env,&
      71             :                                               md_environment_type,&
      72             :                                               set_md_env
      73             :    USE message_passing,                 ONLY: mp_para_env_type
      74             :    USE metadynamics,                    ONLY: metadyn_integrator,&
      75             :                                               metadyn_velocities_colvar
      76             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      77             :    USE molecule_kind_types,             ONLY: local_fixd_constraint_type,&
      78             :                                               molecule_kind_type
      79             :    USE molecule_list_types,             ONLY: molecule_list_type
      80             :    USE molecule_types,                  ONLY: global_constraint_type,&
      81             :                                               molecule_type
      82             :    USE particle_list_types,             ONLY: particle_list_type
      83             :    USE particle_types,                  ONLY: particle_type,&
      84             :                                               update_particle_set
      85             :    USE physcon,                         ONLY: femtoseconds
      86             :    USE qmmm_util,                       ONLY: apply_qmmm_walls_reflective
      87             :    USE qmmmx_update,                    ONLY: qmmmx_update_force_env
      88             :    USE qs_environment_types,            ONLY: get_qs_env
      89             :    USE reftraj_types,                   ONLY: REFTRAJ_EVAL_ENERGY_FORCES,&
      90             :                                               REFTRAJ_EVAL_NONE,&
      91             :                                               reftraj_type
      92             :    USE reftraj_util,                    ONLY: compute_msd_reftraj
      93             :    USE rt_propagation_methods,          ONLY: propagation_step
      94             :    USE rt_propagation_output,           ONLY: rt_prop_output
      95             :    USE rt_propagation_types,            ONLY: rt_prop_type
      96             :    USE shell_opt,                       ONLY: optimize_shell_core
      97             :    USE simpar_types,                    ONLY: simpar_type
      98             :    USE string_utilities,                ONLY: uppercase
      99             :    USE thermal_region_types,            ONLY: thermal_region_type,&
     100             :                                               thermal_regions_type
     101             :    USE thermostat_methods,              ONLY: apply_thermostat_baro,&
     102             :                                               apply_thermostat_particles,&
     103             :                                               apply_thermostat_shells
     104             :    USE thermostat_types,                ONLY: thermostat_type
     105             :    USE virial_methods,                  ONLY: virial_evaluate
     106             :    USE virial_types,                    ONLY: virial_type
     107             : #include "../base/base_uses.f90"
     108             : 
     109             :    IMPLICIT NONE
     110             : 
     111             :    PRIVATE
     112             : 
     113             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator'
     114             : 
     115             :    PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_f, nve_respa
     116             :    PUBLIC :: nph_uniaxial_damped, nph_uniaxial, nvt_adiabatic, reftraj
     117             : 
     118             : CONTAINS
     119             : 
     120             : ! **************************************************************************************************
     121             : !> \brief Langevin integrator for particle positions & momenta (Brownian dynamics)
     122             : !> \param md_env ...
     123             : !> \par Literature
     124             : !>      - A. Ricci and G. Ciccotti, Mol. Phys. 101, 1927-1931 (2003)
     125             : !>      - For langevin regions:
     126             : !>        - L. Kantorovich, Phys. Rev. B 78, 094304 (2008)
     127             : !>        - L. Kantorovich and N. Rompotis, Phys. Rev. B 78, 094305 (2008)
     128             : !> \par History
     129             : !>   - Created (01.07.2005,MK)
     130             : !>   - Added support for only performing Langevin MD on a region of atoms
     131             : !>     (01.12.2013, LT)
     132             : !> \author Matthias Krack
     133             : ! **************************************************************************************************
     134         222 :    SUBROUTINE langevin(md_env)
     135             : 
     136             :       TYPE(md_environment_type), POINTER                 :: md_env
     137             : 
     138             :       INTEGER :: iparticle, iparticle_kind, iparticle_local, iparticle_reg, ireg, nparticle, &
     139             :          nparticle_kind, nparticle_local, nshell
     140             :       INTEGER, POINTER                                   :: itimes
     141         222 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: do_langevin
     142             :       REAL(KIND=dp)                                      :: c, c1, c2, c3, c4, dm, dt, gam, mass, &
     143             :                                                             noisy_gamma_region, reg_temp, sigma
     144         222 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: var_w
     145         222 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pos, vel, w
     146             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     147         222 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     148             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     149             :       TYPE(cell_type), POINTER                           :: cell
     150             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     151             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     152             :       TYPE(force_env_type), POINTER                      :: force_env
     153             :       TYPE(global_constraint_type), POINTER              :: gci
     154             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     155         222 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     156             :       TYPE(molecule_list_type), POINTER                  :: molecules
     157         222 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     158             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     159             :       TYPE(particle_list_type), POINTER                  :: particles
     160         222 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     161             :       TYPE(simpar_type), POINTER                         :: simpar
     162             :       TYPE(thermal_region_type), POINTER                 :: thermal_region
     163             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     164             :       TYPE(virial_type), POINTER                         :: virial
     165             : 
     166         222 :       NULLIFY (cell, para_env, gci, force_env)
     167         222 :       NULLIFY (atomic_kinds, local_particles, subsys, local_molecules, molecule_kinds, molecules)
     168         222 :       NULLIFY (molecule_kind_set, molecule_set, particles, particle_set, simpar, virial)
     169         222 :       NULLIFY (thermal_region, thermal_regions, itimes)
     170             : 
     171             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
     172             :                       para_env=para_env, thermal_regions=thermal_regions, &
     173         222 :                       itimes=itimes)
     174             : 
     175         222 :       dt = simpar%dt
     176         222 :       gam = simpar%gamma + simpar%shadow_gamma
     177             :       nshell = 0
     178             : 
     179         222 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
     180             : 
     181             :       ! Do some checks on coordinates and box
     182         222 :       CALL apply_qmmm_walls_reflective(force_env)
     183             : 
     184             :       CALL cp_subsys_get(subsys=subsys, &
     185             :                          atomic_kinds=atomic_kinds, &
     186             :                          gci=gci, &
     187             :                          local_particles=local_particles, &
     188             :                          local_molecules=local_molecules, &
     189             :                          molecules=molecules, &
     190             :                          molecule_kinds=molecule_kinds, &
     191             :                          nshell=nshell, &
     192             :                          particles=particles, &
     193         222 :                          virial=virial)
     194         222 :       IF (nshell /= 0) &
     195           0 :          CPABORT("Langevin dynamics is not yet implemented for core-shell models")
     196             : 
     197         222 :       nparticle_kind = atomic_kinds%n_els
     198         222 :       atomic_kind_set => atomic_kinds%els
     199         222 :       molecule_kind_set => molecule_kinds%els
     200             : 
     201         222 :       nparticle = particles%n_els
     202         222 :       particle_set => particles%els
     203         222 :       molecule_set => molecules%els
     204             : 
     205             :       ! Setup the langevin regions information
     206         666 :       ALLOCATE (do_langevin(nparticle))
     207         222 :       IF (simpar%do_thermal_region) THEN
     208         392 :          DO iparticle = 1, nparticle
     209         392 :             do_langevin(iparticle) = thermal_regions%do_langevin(iparticle)
     210             :          END DO
     211             :       ELSE
     212       15604 :          do_langevin(1:nparticle) = .TRUE.
     213             :       END IF
     214             : 
     215             :       ! Allocate the temperature dependent variance (var_w) of the
     216             :       ! random variable for each atom. It may be different for different
     217             :       ! atoms because of the possibility of Langevin regions, and var_w
     218             :       ! for each region should depend on the temperature defined in the
     219             :       ! region
     220             :       ! RZK explains: sigma is the variance of the Wiener process associated
     221             :       ! with the stochastic term, sigma = m*var_w = m*(2*k_B*T*gamma*dt),
     222             :       ! noisy_gamma adds excessive noise that is not balanced by the damping term
     223         666 :       ALLOCATE (var_w(nparticle))
     224       15996 :       var_w(1:nparticle) = simpar%var_w
     225         222 :       IF (simpar%do_thermal_region) THEN
     226         136 :          DO ireg = 1, thermal_regions%nregions
     227          80 :             thermal_region => thermal_regions%thermal_region(ireg)
     228          80 :             noisy_gamma_region = thermal_region%noisy_gamma_region
     229         384 :             DO iparticle_reg = 1, thermal_region%npart
     230         248 :                iparticle = thermal_region%part_index(iparticle_reg)
     231         248 :                reg_temp = thermal_region%temp_expected
     232         328 :                var_w(iparticle) = 2.0_dp*reg_temp*simpar%dt*(simpar%gamma + noisy_gamma_region)
     233             :             END DO
     234             :          END DO
     235             :       END IF
     236             : 
     237             :       ! Allocate work storage
     238         666 :       ALLOCATE (pos(3, nparticle))
     239       63318 :       pos(:, :) = 0.0_dp
     240             : 
     241         666 :       ALLOCATE (vel(3, nparticle))
     242       63318 :       vel(:, :) = 0.0_dp
     243             : 
     244         666 :       ALLOCATE (w(3, nparticle))
     245       63318 :       w(:, :) = 0.0_dp
     246             : 
     247         222 :       IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, &
     248           4 :                                          molecule_kind_set, particle_set, cell)
     249             : 
     250             :       ! Generate random variables
     251         666 :       DO iparticle_kind = 1, nparticle_kind
     252         444 :          atomic_kind => atomic_kind_set(iparticle_kind)
     253         444 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     254         444 :          nparticle_local = local_particles%n_el(iparticle_kind)
     255        8553 :          DO iparticle_local = 1, nparticle_local
     256        7887 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     257        8331 :             IF (do_langevin(iparticle)) THEN
     258        7863 :                sigma = var_w(iparticle)*mass
     259             :                ASSOCIATE (rng_stream => local_particles%local_particle_set(iparticle_kind)% &
     260             :                           rng(iparticle_local))
     261       15726 :                   w(1, iparticle) = rng_stream%stream%next(variance=sigma)
     262        7863 :                   w(2, iparticle) = rng_stream%stream%next(variance=sigma)
     263       15726 :                   w(3, iparticle) = rng_stream%stream%next(variance=sigma)
     264             :                END ASSOCIATE
     265             :             END IF
     266             :          END DO
     267             :       END DO
     268             : 
     269         222 :       DEALLOCATE (var_w)
     270             : 
     271             :       ! Apply fix atom constraint
     272         222 :       CALL fix_atom_control(force_env, w)
     273             : 
     274             :       ! Velocity Verlet (first part)
     275         222 :       c = EXP(-0.25_dp*dt*gam)
     276         222 :       c2 = c*c
     277         222 :       c4 = c2*c2
     278         222 :       c1 = dt*c2
     279             : 
     280         666 :       DO iparticle_kind = 1, nparticle_kind
     281         444 :          atomic_kind => atomic_kind_set(iparticle_kind)
     282         444 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     283         444 :          nparticle_local = local_particles%n_el(iparticle_kind)
     284         444 :          dm = 0.5_dp*dt/mass
     285         444 :          c3 = dm/c2
     286        8553 :          DO iparticle_local = 1, nparticle_local
     287        7887 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     288        8331 :             IF (do_langevin(iparticle)) THEN
     289             :                vel(:, iparticle) = particle_set(iparticle)%v(:) + &
     290       31452 :                                    c3*particle_set(iparticle)%f(:)
     291             :                pos(:, iparticle) = particle_set(iparticle)%r(:) + &
     292             :                                    c1*particle_set(iparticle)%v(:) + &
     293             :                                    c*dm*(dt*particle_set(iparticle)%f(:) + &
     294       31452 :                                          w(:, iparticle))
     295             :             ELSE
     296             :                vel(:, iparticle) = particle_set(iparticle)%v(:) + &
     297          96 :                                    dm*particle_set(iparticle)%f(:)
     298             :                pos(:, iparticle) = particle_set(iparticle)%r(:) + &
     299             :                                    dt*particle_set(iparticle)%v(:) + &
     300          96 :                                    dm*dt*particle_set(iparticle)%f(:)
     301             :             END IF
     302             :          END DO
     303             :       END DO
     304             : 
     305         222 :       IF (simpar%constraint) THEN
     306             :          ! Possibly update the target values
     307             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
     308           4 :                                    molecule_kind_set, dt, force_env%root_section)
     309             : 
     310             :          CALL shake_control(gci, local_molecules, molecule_set, molecule_kind_set, &
     311             :                             particle_set, pos, vel, dt, simpar%shake_tol, &
     312             :                             simpar%info_constraint, simpar%lagrange_multipliers, &
     313           4 :                             simpar%dump_lm, cell, para_env, local_particles)
     314             :       END IF
     315             : 
     316             :       ! Broadcast the new particle positions
     317         222 :       CALL update_particle_set(particle_set, para_env, pos=pos)
     318             : 
     319         222 :       DEALLOCATE (pos)
     320             : 
     321             :       ! Update forces
     322         222 :       CALL force_env_calc_energy_force(force_env)
     323             : 
     324             :       ! Metadynamics
     325         222 :       CALL metadyn_integrator(force_env, itimes, vel)
     326             : 
     327             :       ! Update Verlet (second part)
     328         666 :       DO iparticle_kind = 1, nparticle_kind
     329         444 :          atomic_kind => atomic_kind_set(iparticle_kind)
     330         444 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     331         444 :          dm = 0.5_dp*dt/mass
     332         444 :          c3 = dm/c2
     333         444 :          nparticle_local = local_particles%n_el(iparticle_kind)
     334        8553 :          DO iparticle_local = 1, nparticle_local
     335        7887 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     336        8331 :             IF (do_langevin(iparticle)) THEN
     337        7863 :                vel(1, iparticle) = vel(1, iparticle) + c3*particle_set(iparticle)%f(1)
     338        7863 :                vel(2, iparticle) = vel(2, iparticle) + c3*particle_set(iparticle)%f(2)
     339        7863 :                vel(3, iparticle) = vel(3, iparticle) + c3*particle_set(iparticle)%f(3)
     340        7863 :                vel(1, iparticle) = c4*vel(1, iparticle) + c2*w(1, iparticle)/mass
     341        7863 :                vel(2, iparticle) = c4*vel(2, iparticle) + c2*w(2, iparticle)/mass
     342        7863 :                vel(3, iparticle) = c4*vel(3, iparticle) + c2*w(3, iparticle)/mass
     343             :             ELSE
     344          24 :                vel(1, iparticle) = vel(1, iparticle) + dm*particle_set(iparticle)%f(1)
     345          24 :                vel(2, iparticle) = vel(2, iparticle) + dm*particle_set(iparticle)%f(2)
     346          24 :                vel(3, iparticle) = vel(3, iparticle) + dm*particle_set(iparticle)%f(3)
     347             :             END IF
     348             :          END DO
     349             :       END DO
     350             : 
     351         222 :       IF (simpar%temperature_annealing) THEN
     352          40 :          simpar%temp_ext = simpar%temp_ext*simpar%f_temperature_annealing
     353          40 :          simpar%var_w = simpar%var_w*simpar%f_temperature_annealing
     354             :       END IF
     355             : 
     356         222 :       IF (simpar%constraint) THEN
     357             :          CALL rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, &
     358             :                              particle_set, vel, dt, simpar%shake_tol, &
     359             :                              simpar%info_constraint, simpar%lagrange_multipliers, &
     360           4 :                              simpar%dump_lm, cell, para_env, local_particles)
     361             :       END IF
     362             : 
     363             :       ! Broadcast the new particle velocities
     364         222 :       CALL update_particle_set(particle_set, para_env, vel=vel)
     365             : 
     366         222 :       DEALLOCATE (vel)
     367             : 
     368         222 :       DEALLOCATE (w)
     369             : 
     370         222 :       DEALLOCATE (do_langevin)
     371             : 
     372             :       ! Update virial
     373         222 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, molecule_set, &
     374           4 :                                                 molecule_kind_set, particle_set, virial, para_env)
     375             : 
     376             :       CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
     377         222 :                            virial, para_env)
     378             : 
     379         444 :    END SUBROUTINE langevin
     380             : 
     381             : ! **************************************************************************************************
     382             : !> \brief nve integrator for particle positions & momenta
     383             : !> \param md_env ...
     384             : !> \param globenv ...
     385             : !> \par History
     386             : !>   - the local particle lists are used instead of pnode (Sep. 2003,MK)
     387             : !>   - usage of fragments retrieved from the force environment (Oct. 2003,MK)
     388             : !> \author CJM
     389             : ! **************************************************************************************************
     390       32737 :    SUBROUTINE nve(md_env, globenv)
     391             : 
     392             :       TYPE(md_environment_type), POINTER                 :: md_env
     393             :       TYPE(global_environment_type), POINTER             :: globenv
     394             : 
     395             :       INTEGER                                            :: i_iter, n_iter, nparticle, &
     396             :                                                             nparticle_kind, nshell
     397             :       INTEGER, POINTER                                   :: itimes
     398             :       LOGICAL                                            :: deallocate_vel, ehrenfest_md, &
     399             :                                                             shell_adiabatic, shell_check_distance, &
     400             :                                                             shell_present
     401             :       REAL(KIND=dp)                                      :: dt
     402       32737 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: v_old
     403             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     404       32737 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     405             :       TYPE(cell_type), POINTER                           :: cell
     406             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     407             :       TYPE(dft_control_type), POINTER                    :: dft_control
     408             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     409             :       TYPE(force_env_type), POINTER                      :: force_env
     410             :       TYPE(global_constraint_type), POINTER              :: gci
     411             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     412       32737 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     413             :       TYPE(molecule_list_type), POINTER                  :: molecules
     414       32737 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     415             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     416             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     417             :                                                             shell_particles
     418       32737 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     419       32737 :                                                             shell_particle_set
     420             :       TYPE(rt_prop_type), POINTER                        :: rtp
     421             :       TYPE(simpar_type), POINTER                         :: simpar
     422             :       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_shell
     423             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     424             :       TYPE(virial_type), POINTER                         :: virial
     425             : 
     426       32737 :       NULLIFY (thermostat_coeff, tmp)
     427       32737 :       NULLIFY (subsys, simpar, para_env, cell, gci, force_env, virial)
     428       32737 :       NULLIFY (atomic_kinds, local_particles, molecules, molecule_kind_set, molecule_set, particle_set)
     429       32737 :       NULLIFY (shell_particles, shell_particle_set, core_particles, &
     430       32737 :                core_particle_set, thermostat_shell, dft_control, itimes)
     431             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
     432             :                       thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell, &
     433       32737 :                       para_env=para_env, ehrenfest_md=ehrenfest_md, itimes=itimes)
     434       32737 :       dt = simpar%dt
     435       32737 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
     436             : 
     437             :       ! Do some checks on coordinates and box
     438       32737 :       CALL apply_qmmm_walls_reflective(force_env)
     439             : 
     440             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     441             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
     442       32737 :                          molecule_kinds=molecule_kinds, gci=gci, virial=virial)
     443             : 
     444       32737 :       nparticle_kind = atomic_kinds%n_els
     445       32737 :       atomic_kind_set => atomic_kinds%els
     446       32737 :       molecule_kind_set => molecule_kinds%els
     447             : 
     448       32737 :       nparticle = particles%n_els
     449       32737 :       particle_set => particles%els
     450       32737 :       molecule_set => molecules%els
     451             : 
     452             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     453             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
     454       32737 :                                shell_check_distance=shell_check_distance)
     455             : 
     456       32737 :       IF (shell_present) THEN
     457             :          CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, &
     458         600 :                             core_particles=core_particles)
     459         600 :          shell_particle_set => shell_particles%els
     460         600 :          nshell = SIZE(shell_particles%els)
     461             : 
     462         600 :          IF (shell_adiabatic) THEN
     463         600 :             core_particle_set => core_particles%els
     464             :          END IF
     465             :       END IF
     466             : 
     467       32737 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
     468             : 
     469             :       ! Apply thermostat over the full set of shells if required
     470       32737 :       IF (shell_adiabatic) THEN
     471             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
     472             :                                       local_particles, para_env, shell_particle_set=shell_particle_set, &
     473         600 :                                       core_particle_set=core_particle_set)
     474             :       END IF
     475             : 
     476       32737 :       IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, &
     477       13228 :                                          molecule_kind_set, particle_set, cell)
     478             : 
     479             :       ! Velocity Verlet (first part)
     480             :       CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
     481       32737 :                     core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
     482             : 
     483       32737 :       IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, &
     484             :                                                      local_particles, particle_set, core_particle_set, shell_particle_set, &
     485         280 :                                                      nparticle_kind, shell_adiabatic)
     486             : 
     487       32737 :       IF (simpar%constraint) THEN
     488             :          ! Possibly update the target values
     489             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
     490       13228 :                                    molecule_kind_set, dt, force_env%root_section)
     491             : 
     492             :          CALL shake_control(gci, local_molecules, molecule_set, &
     493             :                             molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, &
     494             :                             simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
     495       13228 :                             cell, para_env, local_particles)
     496             :       END IF
     497             : 
     498             :       ! Broadcast the new particle positions and deallocate pos part of temporary
     499             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     500       32737 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
     501             : 
     502       32737 :       IF (shell_adiabatic .AND. shell_check_distance) THEN
     503             :          CALL optimize_shell_core(force_env, particle_set, &
     504         180 :                                   shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
     505             :       END IF
     506             : 
     507             :       ! Update forces
     508             :       ! In case of ehrenfest dynamics, velocities need to be iterated
     509       32737 :       IF (ehrenfest_md) THEN
     510         810 :          ALLOCATE (v_old(3, SIZE(tmp%vel, 2)))
     511        3414 :          v_old(:, :) = tmp%vel
     512             :          CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
     513         270 :                         core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
     514             :          CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     515             :                                  core_particle_set, para_env, shell_adiabatic, vel=.TRUE., &
     516         270 :                                  should_deall_vel=.FALSE.)
     517        3414 :          tmp%vel = v_old
     518         270 :          CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
     519         270 :          n_iter = dft_control%rtp_control%max_iter
     520             :       ELSE
     521             :          n_iter = 1
     522             :       END IF
     523             : 
     524       66084 :       DO i_iter = 1, n_iter
     525             : 
     526       33617 :          IF (ehrenfest_md) THEN
     527        1150 :             CALL get_qs_env(qs_env=force_env%qs_env, rtp=rtp)
     528        1150 :             rtp%iter = i_iter
     529       14742 :             tmp%vel = v_old
     530        1150 :             CALL propagation_step(force_env%qs_env, rtp, dft_control%rtp_control)
     531             :          END IF
     532             : 
     533             :          ![NB] let nve work with force mixing which does not have consistent energies and forces
     534       33617 :          CALL force_env_calc_energy_force(force_env, require_consistent_energy_force=.FALSE.)
     535             : 
     536       33617 :          IF (ehrenfest_md) THEN
     537        1150 :             CALL rt_prop_output(force_env%qs_env, ehrenfest, delta_iter=force_env%qs_env%rtp%delta_iter)
     538             :          END IF
     539             : 
     540             :          ! Metadynamics
     541       33617 :          CALL metadyn_integrator(force_env, itimes, tmp%vel)
     542             : 
     543             :          ! Velocity Verlet (second part)
     544             :          CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
     545       33617 :                         core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
     546             : 
     547       33617 :          IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, &
     548             :                                                     molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, &
     549             :                                                     simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
     550       13228 :                                                     cell, para_env, local_particles)
     551             : 
     552             :          ! Apply thermostat over the full set of shell if required
     553       33617 :          IF (shell_adiabatic) THEN
     554             :             CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
     555             :                                          local_particles, para_env, vel=tmp%vel, &
     556         600 :                                          shell_vel=tmp%shell_vel, core_vel=tmp%core_vel)
     557             :          END IF
     558             : 
     559       33617 :          IF (simpar%annealing) THEN
     560           0 :             tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
     561           0 :             IF (shell_adiabatic) THEN
     562             :                CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
     563           0 :                                      tmp%vel, tmp%shell_vel, tmp%core_vel)
     564             :             END IF
     565             :          END IF
     566             : 
     567       33617 :          IF (ehrenfest_md) deallocate_vel = force_env%qs_env%rtp%converged
     568       33617 :          IF (i_iter .EQ. n_iter) deallocate_vel = .TRUE.
     569             :          ! Broadcast the new particle velocities and deallocate the full temporary
     570             :          CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     571             :                                  core_particle_set, para_env, shell_adiabatic, vel=.TRUE., &
     572       33617 :                                  should_deall_vel=deallocate_vel)
     573       66354 :          IF (ehrenfest_md) THEN
     574        1150 :             IF (force_env%qs_env%rtp%converged) EXIT
     575             :          END IF
     576             : 
     577             :       END DO
     578             : 
     579             :       ! Update virial
     580       32737 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
     581       13228 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
     582             : 
     583             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
     584       32737 :                            local_particles, virial, para_env)
     585             : 
     586       65474 :    END SUBROUTINE nve
     587             : 
     588             : ! **************************************************************************************************
     589             : !> \brief simplest version of the isokinetic gaussian thermostat
     590             : !> \param md_env ...
     591             : !> \par History
     592             : !>   - Created [2004-07]
     593             : !> \author Joost VandeVondele
     594             : !> \note
     595             : !>      - time reversible, and conserves the kinetic energy to machine precision
     596             : !>      - is not yet supposed to work for e.g. constraints, our the extended version
     597             : !>        of this thermostat
     598             : !>        see:
     599             : !>         - Zhang F. , JCP 106, 6102 (1997)
     600             : !>         - Minary P. et al, JCP 118, 2510 (2003)
     601             : ! **************************************************************************************************
     602          12 :    SUBROUTINE isokin(md_env)
     603             : 
     604             :       TYPE(md_environment_type), POINTER                 :: md_env
     605             : 
     606             :       INTEGER                                            :: nparticle, nparticle_kind, nshell
     607             :       INTEGER, POINTER                                   :: itimes
     608             :       LOGICAL                                            :: shell_adiabatic, shell_present
     609             :       REAL(KIND=dp)                                      :: dt
     610             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     611           6 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     612             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     613             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     614             :       TYPE(force_env_type), POINTER                      :: force_env
     615             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     616             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     617             :                                                             shell_particles
     618           6 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     619           6 :                                                             shell_particle_set
     620             :       TYPE(simpar_type), POINTER                         :: simpar
     621             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     622             : 
     623           6 :       NULLIFY (force_env, tmp, simpar, itimes)
     624           6 :       NULLIFY (atomic_kinds, para_env, subsys, local_particles)
     625           6 :       NULLIFY (core_particles, particles, shell_particles)
     626           6 :       NULLIFY (core_particle_set, particle_set, shell_particle_set)
     627             : 
     628             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
     629           6 :                       para_env=para_env, itimes=itimes)
     630             : 
     631           6 :       dt = simpar%dt
     632             : 
     633           6 :       CALL force_env_get(force_env=force_env, subsys=subsys)
     634             : 
     635             :       ! Do some checks on coordinates and box
     636           6 :       CALL apply_qmmm_walls_reflective(force_env)
     637             : 
     638           6 :       IF (simpar%constraint) THEN
     639           0 :          CPABORT("Constraints not yet implemented")
     640             :       END IF
     641             : 
     642             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, &
     643             :                          local_particles=local_particles, &
     644           6 :                          particles=particles)
     645             : 
     646           6 :       nparticle_kind = atomic_kinds%n_els
     647           6 :       atomic_kind_set => atomic_kinds%els
     648           6 :       nparticle = particles%n_els
     649           6 :       particle_set => particles%els
     650             : 
     651             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     652           6 :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic)
     653             : 
     654           6 :       IF (shell_present) THEN
     655             :          CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, &
     656           0 :                             core_particles=core_particles)
     657           0 :          shell_particle_set => shell_particles%els
     658           0 :          nshell = SIZE(shell_particles%els)
     659             : 
     660           0 :          IF (shell_adiabatic) THEN
     661           0 :             core_particle_set => core_particles%els
     662             :          END IF
     663             :       END IF
     664             : 
     665           6 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
     666             : 
     667             :       ! compute s,ds
     668             :       CALL get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, &
     669           6 :                     dt, para_env)
     670             : 
     671             :       ! Velocity Verlet (first part)
     672          24 :       tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds)
     673          24 :       tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt
     674             :       CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
     675             :                     core_particle_set, shell_particle_set, nparticle_kind, &
     676           6 :                     shell_adiabatic, dt)
     677             : 
     678           6 :       IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, &
     679             :                                                      local_particles, particle_set, core_particle_set, shell_particle_set, &
     680           0 :                                                      nparticle_kind, shell_adiabatic)
     681             : 
     682             :       ! Broadcast the new particle positions and deallocate the pos components of temporary
     683             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     684           6 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
     685             : 
     686           6 :       CALL force_env_calc_energy_force(force_env)
     687             : 
     688             :       ! Metadynamics
     689           6 :       CALL metadyn_integrator(force_env, itimes, tmp%vel)
     690             : 
     691             :       ! compute s,ds
     692             :       CALL get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, &
     693           6 :                     dt, para_env, tmpv=.TRUE.)
     694             : 
     695             :       ! Velocity Verlet (second part)
     696          24 :       tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds)
     697          24 :       tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt
     698             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
     699             :                      core_particle_set, shell_particle_set, nparticle_kind, &
     700           6 :                      shell_adiabatic, dt)
     701             : 
     702           6 :       IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
     703             : 
     704             :       !  Broadcast the new particle velocities and deallocate the temporary
     705             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     706           6 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
     707             : 
     708           6 :    END SUBROUTINE isokin
     709             : ! **************************************************************************************************
     710             : !> \brief nvt adiabatic integrator for particle positions & momenta
     711             : !> \param md_env ...
     712             : !> \param globenv ...
     713             : !> \par History
     714             : !>   - the local particle lists are used instead of pnode (Sep. 2003,MK)
     715             : !>   - usage of fragments retrieved from the force environment (Oct. 2003,MK)
     716             : !> \author CJM
     717             : ! **************************************************************************************************
     718           0 :    SUBROUTINE nvt_adiabatic(md_env, globenv)
     719             : 
     720             :       TYPE(md_environment_type), POINTER                 :: md_env
     721             :       TYPE(global_environment_type), POINTER             :: globenv
     722             : 
     723             :       INTEGER                                            :: ivar, nparticle, nparticle_kind, nshell
     724             :       INTEGER, POINTER                                   :: itimes
     725             :       LOGICAL                                            :: shell_adiabatic, shell_check_distance, &
     726             :                                                             shell_present
     727             :       REAL(KIND=dp)                                      :: dt
     728           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rand
     729             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     730           0 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     731             :       TYPE(cell_type), POINTER                           :: cell
     732             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     733             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     734             :       TYPE(force_env_type), POINTER                      :: force_env
     735             :       TYPE(global_constraint_type), POINTER              :: gci
     736             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     737           0 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     738             :       TYPE(molecule_list_type), POINTER                  :: molecules
     739           0 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     740             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     741             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     742             :                                                             shell_particles
     743           0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     744           0 :                                                             shell_particle_set
     745             :       TYPE(simpar_type), POINTER                         :: simpar
     746             :       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_fast, &
     747             :                                                             thermostat_shell, thermostat_slow
     748             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     749             :       TYPE(virial_type), POINTER                         :: virial
     750             : 
     751           0 :       NULLIFY (gci, force_env, thermostat_coeff, tmp, &
     752           0 :                thermostat_fast, thermostat_slow, thermostat_shell, cell, shell_particles, &
     753           0 :                shell_particle_set, core_particles, core_particle_set, rand)
     754           0 :       NULLIFY (para_env, subsys, local_molecules, local_particles, molecule_kinds, &
     755           0 :                molecules, molecule_kind_set, molecule_set, atomic_kinds, particles)
     756           0 :       NULLIFY (simpar, itimes)
     757             : 
     758             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
     759             :                       thermostat_fast=thermostat_fast, thermostat_slow=thermostat_slow, &
     760             :                       thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell, &
     761           0 :                       para_env=para_env, itimes=itimes)
     762           0 :       dt = simpar%dt
     763             : 
     764           0 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
     765             : 
     766             :       ! Do some checks on coordinates and box
     767           0 :       CALL apply_qmmm_walls_reflective(force_env)
     768             : 
     769             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     770             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
     771           0 :                          molecule_kinds=molecule_kinds, gci=gci, virial=virial)
     772             : 
     773           0 :       nparticle_kind = atomic_kinds%n_els
     774           0 :       atomic_kind_set => atomic_kinds%els
     775           0 :       molecule_kind_set => molecule_kinds%els
     776             : 
     777           0 :       nparticle = particles%n_els
     778           0 :       particle_set => particles%els
     779           0 :       molecule_set => molecules%els
     780             : 
     781             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     782             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
     783           0 :                                shell_check_distance=shell_check_distance)
     784             : 
     785           0 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     786             :          ! Allocate random number for Langevin Thermostat acting on COLVARS
     787           0 :          IF (force_env%meta_env%langevin) THEN
     788           0 :             ALLOCATE (rand(force_env%meta_env%n_colvar))
     789           0 :             rand(:) = 0.0_dp
     790             :          END IF
     791             :       END IF
     792             : 
     793             :       !  Allocate work storage for positions and velocities
     794           0 :       IF (shell_present) THEN
     795             :          CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, &
     796           0 :                             core_particles=core_particles)
     797           0 :          shell_particle_set => shell_particles%els
     798           0 :          nshell = SIZE(shell_particles%els)
     799             : 
     800           0 :          IF (shell_adiabatic) THEN
     801           0 :             core_particle_set => core_particles%els
     802             :          END IF
     803             :       END IF
     804             : 
     805           0 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
     806             : 
     807             :       ! Apply Thermostat over the full set of particles
     808           0 :       IF (shell_adiabatic) THEN
     809             : !       CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,&
     810             : !            particle_set, local_molecules, para_env, shell_adiabatic=shell_adiabatic,&
     811             : !            shell_particle_set=shell_particle_set, core_particle_set=core_particle_set)
     812             : 
     813             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
     814             :                                       local_particles, para_env, shell_particle_set=shell_particle_set, &
     815           0 :                                       core_particle_set=core_particle_set)
     816             :       ELSE
     817             :          CALL apply_thermostat_particles(thermostat_fast, force_env, molecule_kind_set, molecule_set, &
     818           0 :                                          particle_set, local_molecules, local_particles, para_env)
     819             : 
     820             :          CALL apply_thermostat_particles(thermostat_slow, force_env, molecule_kind_set, molecule_set, &
     821           0 :                                          particle_set, local_molecules, local_particles, para_env)
     822             :       END IF
     823             : 
     824           0 :       IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, &
     825           0 :                                          molecule_kind_set, particle_set, cell)
     826             : 
     827             :       !    *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
     828           0 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     829           0 :          IF (force_env%meta_env%langevin) THEN
     830           0 :             DO ivar = 1, force_env%meta_env%n_colvar
     831           0 :                rand(ivar) = force_env%meta_env%rng(ivar)%next()
     832             :             END DO
     833           0 :             CALL metadyn_velocities_colvar(force_env, rand)
     834             :          END IF
     835             :       END IF
     836             : 
     837             :       ! Velocity Verlet (first part)
     838             :       CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
     839           0 :                     core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
     840             : 
     841           0 :       IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, &
     842             :                                                      local_particles, particle_set, core_particle_set, shell_particle_set, &
     843           0 :                                                      nparticle_kind, shell_adiabatic)
     844             : 
     845           0 :       IF (simpar%constraint) THEN
     846             :          ! Possibly update the target values
     847             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
     848           0 :                                    molecule_kind_set, dt, force_env%root_section)
     849             : 
     850             :          CALL shake_control(gci, local_molecules, molecule_set, &
     851             :                             molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, &
     852             :                             simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
     853           0 :                             cell, para_env, local_particles)
     854             :       END IF
     855             : 
     856             :       ! Broadcast the new particle positions and deallocate pos components of temporary
     857             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     858           0 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
     859             : 
     860           0 :       IF (shell_adiabatic .AND. shell_check_distance) THEN
     861             :          CALL optimize_shell_core(force_env, particle_set, &
     862           0 :                                   shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
     863             :       END IF
     864             : 
     865             :       ! Update forces
     866           0 :       CALL force_env_calc_energy_force(force_env)
     867             : 
     868             :       ! Metadynamics
     869           0 :       CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand)
     870             : 
     871             :       ! Velocity Verlet (second part)
     872             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
     873           0 :                      core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
     874             : 
     875           0 :       IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, &
     876             :                                                  molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, &
     877             :                                                  simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
     878           0 :                                                  cell, para_env, local_particles)
     879             : 
     880             :       ! Apply Thermostat over the full set of particles
     881           0 :       IF (shell_adiabatic) THEN
     882             :          !  CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, &
     883             :          !       particle_set, local_molecules, para_env, shell_adiabatic=shell_adiabatic,&
     884             :          !       vel= tmp%vel, shell_vel= tmp%shell_vel, core_vel= tmp%core_vel)
     885             : 
     886             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
     887             :                                       local_particles, para_env, vel=tmp%vel, shell_vel=tmp%shell_vel, &
     888           0 :                                       core_vel=tmp%core_vel)
     889             :       ELSE
     890             :          CALL apply_thermostat_particles(thermostat_slow, force_env, molecule_kind_set, molecule_set, &
     891           0 :                                          particle_set, local_molecules, local_particles, para_env, vel=tmp%vel)
     892             : 
     893             :          CALL apply_thermostat_particles(thermostat_fast, force_env, molecule_kind_set, molecule_set, &
     894           0 :                                          particle_set, local_molecules, local_particles, para_env, vel=tmp%vel)
     895             :       END IF
     896             : 
     897             :       ! Broadcast the new particle velocities and deallocate temporary
     898             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     899           0 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
     900             : 
     901           0 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     902           0 :          IF (force_env%meta_env%langevin) THEN
     903           0 :             DEALLOCATE (rand)
     904             :          END IF
     905             :       END IF
     906             : 
     907             :       ! Update constraint virial
     908           0 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
     909           0 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
     910             : 
     911             :       !     **  Evaluate Virial
     912             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
     913           0 :                            local_particles, virial, para_env)
     914             : 
     915           0 :    END SUBROUTINE nvt_adiabatic
     916             : 
     917             : ! **************************************************************************************************
     918             : !> \brief nvt integrator for particle positions & momenta
     919             : !> \param md_env ...
     920             : !> \param globenv ...
     921             : !> \par History
     922             : !>   - the local particle lists are used instead of pnode (Sep. 2003,MK)
     923             : !>   - usage of fragments retrieved from the force environment (Oct. 2003,MK)
     924             : !> \author CJM
     925             : ! **************************************************************************************************
     926       22044 :    SUBROUTINE nvt(md_env, globenv)
     927             : 
     928             :       TYPE(md_environment_type), POINTER                 :: md_env
     929             :       TYPE(global_environment_type), POINTER             :: globenv
     930             : 
     931             :       INTEGER                                            :: ivar, nparticle, nparticle_kind, nshell
     932             :       INTEGER, POINTER                                   :: itimes
     933             :       LOGICAL                                            :: shell_adiabatic, shell_check_distance, &
     934             :                                                             shell_present
     935             :       REAL(KIND=dp)                                      :: dt
     936        7348 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rand
     937             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     938        7348 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     939             :       TYPE(cell_type), POINTER                           :: cell
     940             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     941             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     942             :       TYPE(force_env_type), POINTER                      :: force_env
     943             :       TYPE(global_constraint_type), POINTER              :: gci
     944             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     945        7348 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     946             :       TYPE(molecule_list_type), POINTER                  :: molecules
     947        7348 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     948             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     949             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     950             :                                                             shell_particles
     951        7348 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     952        7348 :                                                             shell_particle_set
     953             :       TYPE(simpar_type), POINTER                         :: simpar
     954             :       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_part, &
     955             :                                                             thermostat_shell
     956             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     957             :       TYPE(virial_type), POINTER                         :: virial
     958             : 
     959        7348 :       NULLIFY (gci, force_env, thermostat_coeff, tmp, &
     960        7348 :                thermostat_part, thermostat_shell, cell, shell_particles, &
     961        7348 :                shell_particle_set, core_particles, core_particle_set, rand)
     962        7348 :       NULLIFY (para_env, subsys, local_molecules, local_particles, molecule_kinds, &
     963        7348 :                molecules, molecule_kind_set, molecule_set, atomic_kinds, particles)
     964        7348 :       NULLIFY (simpar, thermostat_coeff, thermostat_part, thermostat_shell, itimes)
     965             : 
     966             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
     967             :                       thermostat_part=thermostat_part, thermostat_coeff=thermostat_coeff, &
     968             :                       thermostat_shell=thermostat_shell, para_env=para_env, &
     969        7348 :                       itimes=itimes)
     970        7348 :       dt = simpar%dt
     971             : 
     972        7348 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
     973             : 
     974             :       ! Do some checks on coordinates and box
     975        7348 :       CALL apply_qmmm_walls_reflective(force_env)
     976             : 
     977             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     978             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
     979        7348 :                          molecule_kinds=molecule_kinds, gci=gci, virial=virial)
     980             : 
     981        7348 :       nparticle_kind = atomic_kinds%n_els
     982        7348 :       atomic_kind_set => atomic_kinds%els
     983        7348 :       molecule_kind_set => molecule_kinds%els
     984             : 
     985        7348 :       nparticle = particles%n_els
     986        7348 :       particle_set => particles%els
     987        7348 :       molecule_set => molecules%els
     988             : 
     989             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     990             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
     991        7348 :                                shell_check_distance=shell_check_distance)
     992             : 
     993        7348 :       IF (ASSOCIATED(force_env%meta_env)) THEN
     994             :          ! Allocate random number for Langevin Thermostat acting on COLVARS
     995        1014 :          IF (force_env%meta_env%langevin) THEN
     996         720 :             ALLOCATE (rand(force_env%meta_env%n_colvar))
     997         720 :             rand(:) = 0.0_dp
     998             :          END IF
     999             :       END IF
    1000             : 
    1001             :       !  Allocate work storage for positions and velocities
    1002        7348 :       IF (shell_present) THEN
    1003             :          CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, &
    1004         920 :                             core_particles=core_particles)
    1005         920 :          shell_particle_set => shell_particles%els
    1006         920 :          nshell = SIZE(shell_particles%els)
    1007             : 
    1008         920 :          IF (shell_adiabatic) THEN
    1009         920 :             core_particle_set => core_particles%els
    1010             :          END IF
    1011             :       END IF
    1012             : 
    1013        7348 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
    1014             : 
    1015             :       ! Apply Thermostat over the full set of particles
    1016        7348 :       IF (shell_adiabatic) THEN
    1017             :          CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1018             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    1019         920 :                                          shell_particle_set=shell_particle_set, core_particle_set=core_particle_set)
    1020             : 
    1021             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    1022             :                                       local_particles, para_env, shell_particle_set=shell_particle_set, &
    1023         920 :                                       core_particle_set=core_particle_set)
    1024             :       ELSE
    1025             :          CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1026        6428 :                                          particle_set, local_molecules, local_particles, para_env)
    1027             :       END IF
    1028             : 
    1029        7348 :       IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, &
    1030        2970 :                                          molecule_kind_set, particle_set, cell)
    1031             : 
    1032             :       !    *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
    1033        7348 :       IF (ASSOCIATED(force_env%meta_env)) THEN
    1034        1014 :          IF (force_env%meta_env%langevin) THEN
    1035         720 :             DO ivar = 1, force_env%meta_env%n_colvar
    1036         720 :                rand(ivar) = force_env%meta_env%rng(ivar)%next()
    1037             :             END DO
    1038         240 :             CALL metadyn_velocities_colvar(force_env, rand)
    1039             :          END IF
    1040             :       END IF
    1041             : 
    1042             :       ! Velocity Verlet (first part)
    1043             :       CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1044        7348 :                     core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
    1045             : 
    1046        7348 :       IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, &
    1047             :                                                      local_particles, particle_set, core_particle_set, shell_particle_set, &
    1048           0 :                                                      nparticle_kind, shell_adiabatic)
    1049             : 
    1050        7348 :       IF (simpar%constraint) THEN
    1051             :          ! Possibly update the target values
    1052             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
    1053        2970 :                                    molecule_kind_set, dt, force_env%root_section)
    1054             : 
    1055             :          CALL shake_control(gci, local_molecules, molecule_set, &
    1056             :                             molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, &
    1057             :                             simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
    1058        2970 :                             cell, para_env, local_particles)
    1059             :       END IF
    1060             : 
    1061             :       ! Broadcast the new particle positions and deallocate pos components of temporary
    1062             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1063        7348 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
    1064             : 
    1065        7348 :       IF (shell_adiabatic .AND. shell_check_distance) THEN
    1066             :          CALL optimize_shell_core(force_env, particle_set, &
    1067         280 :                                   shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
    1068             :       END IF
    1069             : 
    1070             :       ![ADAPT] update input structure with new coordinates, make new labels
    1071        7348 :       CALL qmmmx_update_force_env(force_env, force_env%root_section)
    1072             : 
    1073             :       ![NB] recreate pointers changed by creation of new subsys in qmmm_update_force_mixing_env
    1074             :       ![NB] ugly hack, which is why adaptivity isn't implemented in most other ensembles
    1075             :       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1076        7348 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
    1077             : 
    1078             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    1079             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
    1080        7348 :                          molecule_kinds=molecule_kinds, gci=gci, virial=virial)
    1081             : 
    1082        7348 :       nparticle_kind = atomic_kinds%n_els
    1083        7348 :       atomic_kind_set => atomic_kinds%els
    1084        7348 :       molecule_kind_set => molecule_kinds%els
    1085             : 
    1086        7348 :       nparticle = particles%n_els
    1087        7348 :       particle_set => particles%els
    1088        7348 :       molecule_set => molecules%els
    1089             : 
    1090             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    1091             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
    1092        7348 :                                shell_check_distance=shell_check_distance)
    1093             : 
    1094             :       !  Allocate work storage for positions and velocities
    1095        7348 :       IF (shell_present) THEN
    1096             :          CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, &
    1097         920 :                             core_particles=core_particles)
    1098         920 :          shell_particle_set => shell_particles%els
    1099         920 :          nshell = SIZE(shell_particles%els)
    1100             : 
    1101         920 :          IF (shell_adiabatic) THEN
    1102         920 :             core_particle_set => core_particles%els
    1103             :          END IF
    1104             :       END IF
    1105             :       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1106             : 
    1107             :       ! Update forces
    1108             :       ![NB] let nvt work with force mixing which does not have consistent energies and forces
    1109        7348 :       CALL force_env_calc_energy_force(force_env, require_consistent_energy_force=.FALSE.)
    1110             : 
    1111             :       ! Metadynamics
    1112        7348 :       CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand)
    1113             : 
    1114             :       ! Velocity Verlet (second part)
    1115             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    1116        7348 :                      core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
    1117             : 
    1118        7348 :       IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, &
    1119             :                                                  molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, &
    1120             :                                                  simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, &
    1121        2970 :                                                  cell, para_env, local_particles)
    1122             : 
    1123             :       ! Apply Thermostat over the full set of particles
    1124        7348 :       IF (shell_adiabatic) THEN
    1125             :          CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1126             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    1127         920 :                                          vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel)
    1128             : 
    1129             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    1130             :                                       local_particles, para_env, vel=tmp%vel, shell_vel=tmp%shell_vel, &
    1131         920 :                                       core_vel=tmp%core_vel)
    1132             :       ELSE
    1133             :          CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1134        6428 :                                          particle_set, local_molecules, local_particles, para_env, vel=tmp%vel)
    1135             :       END IF
    1136             : 
    1137             :       ! Broadcast the new particle velocities and deallocate temporary
    1138             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1139        7348 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
    1140             : 
    1141        7348 :       IF (ASSOCIATED(force_env%meta_env)) THEN
    1142        1014 :          IF (force_env%meta_env%langevin) THEN
    1143         240 :             DEALLOCATE (rand)
    1144             :          END IF
    1145             :       END IF
    1146             : 
    1147             :       ! Update constraint virial
    1148        7348 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
    1149        2970 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
    1150             : 
    1151             :       !     **  Evaluate Virial
    1152             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
    1153        7348 :                            local_particles, virial, para_env)
    1154             : 
    1155        7348 :    END SUBROUTINE nvt
    1156             : 
    1157             : ! **************************************************************************************************
    1158             : !> \brief npt_i integrator for particle positions & momenta
    1159             : !>      isotropic box changes
    1160             : !> \param md_env ...
    1161             : !> \param globenv ...
    1162             : !> \par History
    1163             : !>      none
    1164             : !> \author CJM
    1165             : ! **************************************************************************************************
    1166        3128 :    SUBROUTINE npt_i(md_env, globenv)
    1167             : 
    1168             :       TYPE(md_environment_type), POINTER                 :: md_env
    1169             :       TYPE(global_environment_type), POINTER             :: globenv
    1170             : 
    1171             :       REAL(KIND=dp), PARAMETER                           :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
    1172             :                                                             e6 = e4/42.0_dp, e8 = e6/72.0_dp
    1173             : 
    1174             :       INTEGER                                            :: iroll, ivar, nkind, nparticle, &
    1175             :                                                             nparticle_kind, nshell
    1176             :       INTEGER, POINTER                                   :: itimes
    1177             :       LOGICAL                                            :: first, first_time, shell_adiabatic, &
    1178             :                                                             shell_check_distance, shell_present
    1179             :       REAL(KIND=dp)                                      :: dt, infree, kin, roll_tol, roll_tol_thrs
    1180             :       REAL(KIND=dp), DIMENSION(3)                        :: vector_r, vector_v
    1181             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_kin
    1182        1564 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rand
    1183             :       REAL(KIND=dp), SAVE                                :: eps_0
    1184             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1185        1564 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1186             :       TYPE(cell_type), POINTER                           :: cell
    1187             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1188             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    1189             :       TYPE(force_env_type), POINTER                      :: force_env
    1190             :       TYPE(global_constraint_type), POINTER              :: gci
    1191             :       TYPE(local_fixd_constraint_type), DIMENSION(:), &
    1192        1564 :          POINTER                                         :: lfixd_list
    1193             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1194        1564 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1195             :       TYPE(molecule_list_type), POINTER                  :: molecules
    1196        1564 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1197             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1198        1564 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
    1199             :       TYPE(old_variables_type), POINTER                  :: old
    1200             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1201             :                                                             shell_particles
    1202        1564 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
    1203        1564 :                                                             shell_particle_set
    1204             :       TYPE(simpar_type), POINTER                         :: simpar
    1205             :       TYPE(thermostat_type), POINTER                     :: thermostat_baro, thermostat_part, &
    1206             :                                                             thermostat_shell
    1207             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1208             :       TYPE(virial_type), POINTER                         :: virial
    1209             : 
    1210        1564 :       NULLIFY (gci, thermostat_baro, thermostat_part, thermostat_shell, force_env)
    1211        1564 :       NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles)
    1212        1564 :       NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt)
    1213        1564 :       NULLIFY (core_particles, particles, shell_particles, tmp, old)
    1214        1564 :       NULLIFY (core_particle_set, particle_set, shell_particle_set)
    1215        1564 :       NULLIFY (simpar, virial, rand, itimes, lfixd_list)
    1216             : 
    1217             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
    1218             :                       thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
    1219             :                       thermostat_shell=thermostat_shell, npt=npt, first_time=first_time, &
    1220        1564 :                       para_env=para_env, itimes=itimes)
    1221        1564 :       dt = simpar%dt
    1222        1564 :       infree = 1.0_dp/REAL(simpar%nfree, KIND=dp)
    1223             : 
    1224        1564 :       CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
    1225             : 
    1226             :       ! Do some checks on coordinates and box
    1227        1564 :       CALL apply_qmmm_walls_reflective(force_env)
    1228             : 
    1229             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    1230             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
    1231        1564 :                          gci=gci, molecule_kinds=molecule_kinds, virial=virial)
    1232             : 
    1233        1564 :       nparticle_kind = atomic_kinds%n_els
    1234        1564 :       nkind = molecule_kinds%n_els
    1235        1564 :       atomic_kind_set => atomic_kinds%els
    1236        1564 :       molecule_kind_set => molecule_kinds%els
    1237             : 
    1238        1564 :       nparticle = particles%n_els
    1239        1564 :       particle_set => particles%els
    1240        1564 :       molecule_set => molecules%els
    1241             : 
    1242             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    1243             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
    1244        1564 :                                shell_check_distance=shell_check_distance)
    1245             : 
    1246        1564 :       IF (first_time) THEN
    1247             :          CALL virial_evaluate(atomic_kind_set, particle_set, &
    1248         108 :                               local_particles, virial, para_env)
    1249             :       END IF
    1250             : 
    1251             :       ! Allocate work storage for positions and velocities
    1252        1564 :       CALL allocate_old(old, particle_set, npt)
    1253             : 
    1254        1564 :       IF (ASSOCIATED(force_env%meta_env)) THEN
    1255             :          ! Allocate random number for Langevin Thermostat acting on COLVARS
    1256           0 :          IF (force_env%meta_env%langevin) THEN
    1257           0 :             ALLOCATE (rand(force_env%meta_env%n_colvar))
    1258           0 :             rand(:) = 0.0_dp
    1259             :          END IF
    1260             :       END IF
    1261             : 
    1262        1564 :       IF (shell_present) THEN
    1263             :          CALL cp_subsys_get(subsys=subsys, &
    1264         120 :                             shell_particles=shell_particles, core_particles=core_particles)
    1265         120 :          shell_particle_set => shell_particles%els
    1266         120 :          nshell = SIZE(shell_particles%els)
    1267         120 :          IF (shell_adiabatic) THEN
    1268         120 :             core_particle_set => core_particles%els
    1269             :          END IF
    1270             :       END IF
    1271             : 
    1272        1564 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
    1273             : 
    1274             :       ! Initialize eps_0 the first time through
    1275        1564 :       IF (first_time) eps_0 = npt(1, 1)%eps
    1276             : 
    1277             :       ! Apply thermostat to  barostat
    1278        1564 :       CALL apply_thermostat_baro(thermostat_baro, npt, para_env)
    1279             : 
    1280             :       ! Apply Thermostat over the full set of particles
    1281        1564 :       IF (simpar%ensemble /= npe_i_ensemble) THEN
    1282        1524 :          IF (shell_adiabatic) THEN
    1283             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1284             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    1285          80 :                                             shell_particle_set=shell_particle_set, core_particle_set=core_particle_set)
    1286             : 
    1287             :          ELSE
    1288             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1289        1444 :                                             particle_set, local_molecules, local_particles, para_env)
    1290             :          END IF
    1291             :       END IF
    1292             : 
    1293             :       ! Apply Thermostat over the core-shell motion
    1294             :       CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    1295             :                                    local_particles, para_env, shell_particle_set=shell_particle_set, &
    1296        1564 :                                    core_particle_set=core_particle_set)
    1297             : 
    1298        1564 :       IF (simpar%constraint) THEN
    1299             :          ! Possibly update the target values
    1300             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
    1301         668 :                                    molecule_kind_set, dt, force_env%root_section)
    1302             :       END IF
    1303             : 
    1304             :       ! setting up for ROLL: saving old variables
    1305        1564 :       IF (simpar%constraint) THEN
    1306         668 :          roll_tol_thrs = simpar%roll_tol
    1307         668 :          iroll = 1
    1308         668 :          CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F')
    1309             :          CALL getold(gci, local_molecules, molecule_set, &
    1310         668 :                      molecule_kind_set, particle_set, cell)
    1311             :       ELSE
    1312             :          roll_tol_thrs = EPSILON(0.0_dp)
    1313             :       END IF
    1314        1564 :       roll_tol = -roll_tol_thrs
    1315             : 
    1316             :       !    *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
    1317        1564 :       IF (ASSOCIATED(force_env%meta_env)) THEN
    1318           0 :          IF (force_env%meta_env%langevin) THEN
    1319           0 :             DO ivar = 1, force_env%meta_env%n_colvar
    1320           0 :                rand(ivar) = force_env%meta_env%rng(ivar)%next()
    1321             :             END DO
    1322           0 :             CALL metadyn_velocities_colvar(force_env, rand)
    1323             :          END IF
    1324             :       END IF
    1325             : 
    1326        4266 :       SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP
    1327             : 
    1328        2702 :          IF (simpar%constraint) THEN
    1329        1806 :             CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B')
    1330             :          END IF
    1331             : 
    1332             :          CALL update_pv(gci, simpar, atomic_kind_set, particle_set, &
    1333             :                         local_molecules, molecule_set, molecule_kind_set, &
    1334        2702 :                         local_particles, kin, pv_kin, virial, para_env)
    1335        2702 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    1336             : 
    1337             :          tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* &
    1338        2702 :                         (0.5_dp*npt(1, 1)%v*dt)
    1339             :          tmp%poly_r(1:3) = 1.0_dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + &
    1340       10808 :                            e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4
    1341             : 
    1342             :          tmp%arg_v(1) = (0.25_dp*npt(1, 1)%v*dt* &
    1343             :                          (1.0_dp + 3.0_dp*infree))*(0.25_dp*npt(1, 1)%v* &
    1344        2702 :                                                     dt*(1.0_dp + 3.0_dp*infree))
    1345             :          tmp%poly_v(1:3) = 1.0_dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + &
    1346       10808 :                            e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4
    1347             : 
    1348       10808 :          tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
    1349             :          tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
    1350       10808 :                                 (1.0_dp + 3.0_dp*infree))
    1351             : 
    1352             :          ! first half of velocity verlet
    1353        2702 :          IF (simpar%ensemble == npt_ia_ensemble) THEN
    1354          20 :             CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
    1355             :             CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1356             :                           core_particle_set, shell_particle_set, nparticle_kind, &
    1357          20 :                           shell_adiabatic, dt, lfixd_list=lfixd_list)
    1358          20 :             CALL release_local_fixd_list(lfixd_list)
    1359             :          ELSE
    1360             :             CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1361             :                           core_particle_set, shell_particle_set, nparticle_kind, &
    1362        2682 :                           shell_adiabatic, dt)
    1363             :          END IF
    1364             : 
    1365        2702 :          IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
    1366             :                                                         atomic_kind_set, local_particles, particle_set, core_particle_set, &
    1367           0 :                                                         shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt)
    1368             : 
    1369        2702 :          roll_tol = 0.0_dp
    1370       10808 :          vector_r(:) = tmp%scale_r(:)*tmp%poly_r(:)
    1371       10808 :          vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:)
    1372             : 
    1373        2702 :          IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, &
    1374             :                                                       molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, &
    1375             :                                                         roll_tol, iroll, vector_r, vector_v, para_env, cell=cell, &
    1376        3370 :                                                         local_particles=local_particles)
    1377             :       END DO SR
    1378             : 
    1379             :       ! Update eps:
    1380        4692 :       npt(:, :)%eps = npt(:, :)%eps + dt*npt(:, :)%v
    1381             : 
    1382             :       ! Update h_mat
    1383       20332 :       cell%hmat(:, :) = cell%hmat(:, :)*EXP(npt(1, 1)%eps - eps_0)
    1384             : 
    1385        1564 :       eps_0 = npt(1, 1)%eps
    1386             : 
    1387             :       ! Update the inverse
    1388        1564 :       CALL init_cell(cell)
    1389             : 
    1390             :       ! Broadcast the new particle positions and deallocate the pos components of temporary
    1391             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1392        1564 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
    1393             : 
    1394        1564 :       IF (shell_adiabatic .AND. shell_check_distance) THEN
    1395             :          CALL optimize_shell_core(force_env, particle_set, &
    1396           0 :                                   shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
    1397             :       END IF
    1398             : 
    1399             :       ! Update forces
    1400        1564 :       CALL force_env_calc_energy_force(force_env)
    1401             : 
    1402             :       ! Metadynamics
    1403        1564 :       CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand)
    1404             : 
    1405             :       ! Velocity Verlet (second part)
    1406             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    1407             :                      core_particle_set, shell_particle_set, nparticle_kind, &
    1408        1564 :                      shell_adiabatic, dt)
    1409             : 
    1410        1564 :       IF (simpar%constraint) THEN
    1411         668 :          roll_tol_thrs = simpar%roll_tol
    1412         668 :          first = .TRUE.
    1413         668 :          iroll = 1
    1414         668 :          CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F')
    1415             :       ELSE
    1416             :          roll_tol_thrs = EPSILON(0.0_dp)
    1417             :       END IF
    1418        1564 :       roll_tol = -roll_tol_thrs
    1419             : 
    1420        4234 :       RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP
    1421        2670 :          roll_tol = 0.0_dp
    1422        2670 :          IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, &
    1423             :                                                        particle_set, local_particles, molecule_kind_set, molecule_set, &
    1424             :                                                        local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, &
    1425        1774 :                                                        roll_tol, iroll, infree, first, para_env)
    1426             : 
    1427             :          CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
    1428             :                         local_molecules, molecule_set, molecule_kind_set, &
    1429        2670 :                         local_particles, kin, pv_kin, virial, para_env)
    1430        4234 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    1431             :       END DO RR
    1432             : 
    1433             :       ! Apply Thermostat over the full set of particles
    1434        1564 :       IF (simpar%ensemble /= npe_i_ensemble) THEN
    1435        1524 :          IF (shell_adiabatic) THEN
    1436             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1437             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    1438          80 :                                             vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel)
    1439             :          ELSE
    1440             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    1441        1444 :                                             particle_set, local_molecules, local_particles, para_env, vel=tmp%vel)
    1442             :          END IF
    1443             :       END IF
    1444             : 
    1445             :       ! Apply Thermostat over the core-shell motion
    1446        1564 :       IF (ASSOCIATED(thermostat_shell)) THEN
    1447             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    1448             :                                       local_particles, para_env, vel=tmp%vel, shell_vel=tmp%shell_vel, &
    1449          40 :                                       core_vel=tmp%core_vel)
    1450             :       END IF
    1451             : 
    1452             :       ! Apply Thermostat to Barostat
    1453        1564 :       CALL apply_thermostat_baro(thermostat_baro, npt, para_env)
    1454             : 
    1455             :       ! Annealing of particle velocities is only possible when no thermostat is active
    1456        1564 :       IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing) THEN
    1457           0 :          tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
    1458           0 :          IF (shell_adiabatic) THEN
    1459             :             CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
    1460           0 :                                   tmp%vel, tmp%shell_vel, tmp%core_vel)
    1461             :          END IF
    1462             :       END IF
    1463             :       ! Annealing of CELL velocities is only possible when no thermostat is active
    1464        1564 :       IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing_cell) THEN
    1465           0 :          npt(1, 1)%v = npt(1, 1)%v*simpar%f_annealing_cell
    1466             :       END IF
    1467             : 
    1468             :       ! Broadcast the new particle velocities and deallocate temporary
    1469             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1470        1564 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
    1471             : 
    1472             :       ! Update constraint virial
    1473        1564 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
    1474         668 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
    1475             : 
    1476             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
    1477        1564 :                            local_particles, virial, para_env)
    1478             : 
    1479             :       ! Deallocate old variables
    1480        1564 :       CALL deallocate_old(old)
    1481             : 
    1482        1564 :       IF (ASSOCIATED(force_env%meta_env)) THEN
    1483           0 :          IF (force_env%meta_env%langevin) THEN
    1484           0 :             DEALLOCATE (rand)
    1485             :          END IF
    1486             :       END IF
    1487             : 
    1488        1564 :       IF (first_time) THEN
    1489         108 :          first_time = .FALSE.
    1490         108 :          CALL set_md_env(md_env, first_time=first_time)
    1491             :       END IF
    1492             : 
    1493        1564 :    END SUBROUTINE npt_i
    1494             : 
    1495             : ! **************************************************************************************************
    1496             : !> \brief uses coordinates in a file and generates frame after frame of these
    1497             : !> \param md_env ...
    1498             : !> \par History
    1499             : !>   - 04.2005 created [Joost VandeVondele]
    1500             : !>   - modified to make it more general [MI]
    1501             : !> \note
    1502             : !>     it can be used to compute some properties on already available trajectories
    1503             : ! **************************************************************************************************
    1504         288 :    SUBROUTINE reftraj(md_env)
    1505             :       TYPE(md_environment_type), POINTER                 :: md_env
    1506             : 
    1507             :       CHARACTER(LEN=2)                                   :: element_kind_ref0, element_symbol, &
    1508             :                                                             element_symbol_ref0
    1509             :       CHARACTER(LEN=max_line_length)                     :: errmsg
    1510             :       INTEGER                                            :: cell_itimes, i, nparticle, Nread, &
    1511             :                                                             trj_itimes
    1512             :       INTEGER, POINTER                                   :: itimes
    1513             :       LOGICAL                                            :: init, my_end, test_ok, traj_has_cell_info
    1514             :       REAL(KIND=dp)                                      :: cell_time, h(3, 3), trj_epot, trj_time, &
    1515             :                                                             vol
    1516             :       REAL(KIND=dp), DIMENSION(3)                        :: abc, albega
    1517             :       REAL(KIND=dp), POINTER                             :: time
    1518             :       TYPE(cell_type), POINTER                           :: cell
    1519             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1520             :       TYPE(force_env_type), POINTER                      :: force_env
    1521             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1522             :       TYPE(particle_list_type), POINTER                  :: particles
    1523         288 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1524             :       TYPE(reftraj_type), POINTER                        :: reftraj_env
    1525             :       TYPE(simpar_type), POINTER                         :: simpar
    1526             : 
    1527         288 :       NULLIFY (reftraj_env, particle_set, particles, force_env, subsys, simpar, para_env, cell, itimes, time)
    1528             :       CALL get_md_env(md_env=md_env, init=init, reftraj=reftraj_env, force_env=force_env, &
    1529         288 :                       para_env=para_env, simpar=simpar)
    1530             : 
    1531         288 :       CALL force_env_get(force_env=force_env, cell=cell, subsys=subsys)
    1532         288 :       reftraj_env%isnap = reftraj_env%isnap + reftraj_env%info%stride
    1533             : 
    1534             :       ! Do some checks on coordinates and box
    1535         288 :       CALL apply_qmmm_walls_reflective(force_env)
    1536         288 :       CALL cp_subsys_get(subsys=subsys, particles=particles)
    1537         288 :       nparticle = particles%n_els
    1538         288 :       particle_set => particles%els
    1539         288 :       abc(:) = 0.0_dp
    1540         288 :       albega(:) = 0.0_dp
    1541             : 
    1542             :       ! SnapShots read from an external file (parsers calls are buffered! please
    1543             :       ! don't put any additional MPI call!) [tlaino]
    1544         288 :       CALL parser_read_line(reftraj_env%info%traj_parser, 1)
    1545         288 :       READ (reftraj_env%info%traj_parser%input_line, FMT="(I8)") nread
    1546         288 :       CALL parser_read_line(reftraj_env%info%traj_parser, 1)
    1547         288 :       test_ok = .FALSE.
    1548         288 :       IF (INDEX(reftraj_env%info%traj_parser%input_line, ", a = ") > 60) THEN
    1549           0 :          traj_has_cell_info = .TRUE.
    1550             :          READ (reftraj_env%info%traj_parser%input_line, &
    1551             :                FMT="(T6,I8,T23,F12.3,T41,F20.10,T67,F14.6,T87,F14.6,T107,F14.6,T131,F8.3,T149,F8.3,T167,F8.3)", &
    1552           0 :                ERR=999) trj_itimes, trj_time, trj_epot, abc(1:3), albega(1:3)
    1553             :          ! Convert cell parameters from angstrom and degree to the internal CP2K units
    1554           0 :          DO i = 1, 3
    1555           0 :             abc(i) = cp_unit_to_cp2k(abc(i), "angstrom")
    1556           0 :             albega(i) = cp_unit_to_cp2k(albega(i), "deg")
    1557             :          END DO
    1558             :       ELSE
    1559         288 :          traj_has_cell_info = .FALSE.
    1560             :          READ (reftraj_env%info%traj_parser%input_line, FMT="(T6,I8,T23,F12.3,T41,F20.10)", ERR=999) &
    1561         288 :             trj_itimes, trj_time, trj_epot
    1562             :       END IF
    1563             :       test_ok = .TRUE.
    1564         288 : 999   IF (.NOT. test_ok) THEN
    1565             :          ! Handling properly the error when reading the title of an XYZ
    1566          50 :          CALL get_md_env(md_env, itimes=itimes)
    1567          50 :          trj_itimes = itimes
    1568          50 :          trj_time = 0.0_dp
    1569          50 :          trj_epot = 0.0_dp
    1570             :       END IF
    1571             :       ! Delayed print of error message until the step number is known
    1572         288 :       IF (nread /= nparticle) THEN
    1573             :          errmsg = "Number of atoms for step "//TRIM(ADJUSTL(cp_to_string(trj_itimes)))// &
    1574             :                   " in the trajectory file does not match the reference configuration: "// &
    1575           0 :                   TRIM(ADJUSTL(cp_to_string(nread)))//" != "//TRIM(ADJUSTL(cp_to_string(nparticle)))
    1576           0 :          CPABORT(errmsg)
    1577             :       END IF
    1578       10392 :       DO i = 1, nread - 1
    1579       10104 :          CALL parser_read_line(reftraj_env%info%traj_parser, 1)
    1580             :          READ (UNIT=reftraj_env%info%traj_parser%input_line(1:LEN_TRIM(reftraj_env%info%traj_parser%input_line)), FMT=*) &
    1581       40416 :             element_symbol, particle_set(i)%r
    1582       10104 :          CALL uppercase(element_symbol)
    1583       10104 :          element_symbol_ref0 = particle_set(i)%atomic_kind%element_symbol
    1584       10104 :          element_kind_ref0 = particle_set(i)%atomic_kind%name
    1585       10104 :          CALL uppercase(element_symbol_ref0)
    1586       10104 :          CALL uppercase(element_kind_ref0)
    1587       10104 :          IF (element_symbol /= element_symbol_ref0) THEN
    1588             :             ! Make sure the kind also does not match a potential alias name
    1589          14 :             IF (element_symbol /= element_kind_ref0) THEN
    1590             :                errmsg = "Atomic configuration from trajectory file does not match the reference configuration: Check atom "// &
    1591           0 :                         TRIM(ADJUSTL(cp_to_string(i)))//" of step "//TRIM(ADJUSTL(cp_to_string(trj_itimes)))
    1592           0 :                CPABORT(errmsg)
    1593             :             END IF
    1594             :          END IF
    1595       10104 :          particle_set(i)%r(1) = cp_unit_to_cp2k(particle_set(i)%r(1), "angstrom")
    1596       10104 :          particle_set(i)%r(2) = cp_unit_to_cp2k(particle_set(i)%r(2), "angstrom")
    1597       10392 :          particle_set(i)%r(3) = cp_unit_to_cp2k(particle_set(i)%r(3), "angstrom")
    1598             :       END DO
    1599             :       ! End of file is properly addressed in the previous call..
    1600             :       ! Let's check directly (providing some info) also for the last
    1601             :       ! line of this frame..
    1602         288 :       CALL parser_read_line(reftraj_env%info%traj_parser, 1, at_end=my_end)
    1603        1152 :       READ (UNIT=reftraj_env%info%traj_parser%input_line, FMT=*) element_symbol, particle_set(i)%r
    1604         288 :       CALL uppercase(element_symbol)
    1605         288 :       element_symbol_ref0 = particle_set(i)%atomic_kind%element_symbol
    1606         288 :       element_kind_ref0 = particle_set(i)%atomic_kind%name
    1607         288 :       CALL uppercase(element_symbol_ref0)
    1608         288 :       CALL uppercase(element_kind_ref0)
    1609         288 :       IF (element_symbol /= element_symbol_ref0) THEN
    1610             :          ! Make sure the kind also does not match a potential alias name
    1611           2 :          IF (element_symbol /= element_kind_ref0) THEN
    1612             :             errmsg = "Atomic configuration from trajectory file does not match the reference configuration: Check atom "// &
    1613           0 :                      TRIM(ADJUSTL(cp_to_string(i)))//" of step "//TRIM(ADJUSTL(cp_to_string(trj_itimes)))
    1614           0 :             CPABORT(errmsg)
    1615             :          END IF
    1616             :       END IF
    1617         288 :       particle_set(i)%r(1) = cp_unit_to_cp2k(particle_set(i)%r(1), "angstrom")
    1618         288 :       particle_set(i)%r(2) = cp_unit_to_cp2k(particle_set(i)%r(2), "angstrom")
    1619         288 :       particle_set(i)%r(3) = cp_unit_to_cp2k(particle_set(i)%r(3), "angstrom")
    1620             : 
    1621             :       ! Check if we reached the end of the file and provide some info..
    1622         288 :       IF (my_end) THEN
    1623           0 :          IF (reftraj_env%isnap /= (simpar%nsteps - 1)) &
    1624             :             CALL cp_abort(__LOCATION__, &
    1625             :                           "Reached the end of the Trajectory  frames in the TRAJECTORY file. Number of "// &
    1626           0 :                           "missing frames ("//cp_to_string((simpar%nsteps - 1) - reftraj_env%isnap)//").")
    1627             :       END IF
    1628             : 
    1629             :       ! Read cell parameters from cell file if requested and if not yet available
    1630         288 :       IF (reftraj_env%info%variable_volume .AND. (.NOT. traj_has_cell_info)) THEN
    1631          62 :          CALL parser_get_next_line(reftraj_env%info%cell_parser, 1, at_end=my_end)
    1632          62 :          CALL parse_cell_line(reftraj_env%info%cell_parser%input_line, cell_itimes, cell_time, h, vol)
    1633          62 :          CPASSERT(trj_itimes == cell_itimes)
    1634             :          ! Check if we reached the end of the file and provide some info..
    1635          62 :          IF (my_end) THEN
    1636           0 :             IF (reftraj_env%isnap /= (simpar%nsteps - 1)) &
    1637             :                CALL cp_abort(__LOCATION__, &
    1638             :                              "Reached the end of the cell info frames in the CELL file. Number of "// &
    1639           0 :                              "missing frames ("//cp_to_string((simpar%nsteps - 1) - reftraj_env%isnap)//").")
    1640             :          END IF
    1641             :       END IF
    1642             : 
    1643         288 :       IF (init) THEN
    1644          38 :          reftraj_env%time0 = trj_time
    1645          38 :          reftraj_env%epot0 = trj_epot
    1646          38 :          reftraj_env%itimes0 = trj_itimes
    1647             :       END IF
    1648             : 
    1649         288 :       IF (trj_itimes /= 0.0_dp .AND. trj_time /= 0.0_dp) simpar%dt = (trj_time/femtoseconds)/REAL(trj_itimes, KIND=dp)
    1650             : 
    1651         288 :       reftraj_env%epot = trj_epot
    1652         288 :       reftraj_env%itimes = trj_itimes
    1653         288 :       reftraj_env%time = trj_time/femtoseconds
    1654         288 :       CALL get_md_env(md_env, t=time)
    1655         288 :       time = reftraj_env%time
    1656             : 
    1657         288 :       IF (traj_has_cell_info) THEN
    1658           0 :          CALL set_cell_param(cell, cell_length=abc, cell_angle=albega, do_init_cell=.TRUE.)
    1659         288 :       ELSE IF (reftraj_env%info%variable_volume) THEN
    1660         806 :          cell%hmat = h
    1661          62 :          CALL init_cell(cell)
    1662             :       END IF
    1663             : 
    1664             :       ![ADAPT] update input structure with new coordinates, make new labels
    1665         288 :       CALL qmmmx_update_force_env(force_env, force_env%root_section)
    1666             :       ! no pointers into force_env%subsys to update
    1667             : 
    1668             :       ! Task to perform on the reference trajectory
    1669             :       ! Compute energy and forces
    1670             :       ![NB] let reftraj work with force mixing which does not have consistent energies and forces
    1671             :       CALL force_env_calc_energy_force(force_env, &
    1672             :                                        calc_force=(reftraj_env%info%eval == REFTRAJ_EVAL_ENERGY_FORCES), &
    1673             :                                        eval_energy_forces=(reftraj_env%info%eval /= REFTRAJ_EVAL_NONE), &
    1674         288 :                                        require_consistent_energy_force=.FALSE.)
    1675             : 
    1676             :       ! Metadynamics
    1677         288 :       CALL metadyn_integrator(force_env, trj_itimes)
    1678             : 
    1679             :       ! Compute MSD with respect to a reference configuration
    1680         288 :       IF (reftraj_env%info%msd) THEN
    1681          14 :          CALL compute_msd_reftraj(reftraj_env, md_env, particle_set)
    1682             :       END IF
    1683             : 
    1684             :       ! Skip according the stride both Trajectory and Cell (if possible)
    1685         288 :       CALL parser_get_next_line(reftraj_env%info%traj_parser, (reftraj_env%info%stride - 1)*(nparticle + 2))
    1686         288 :       IF (reftraj_env%info%variable_volume) THEN
    1687          62 :          CALL parser_get_next_line(reftraj_env%info%cell_parser, (reftraj_env%info%stride - 1))
    1688             :       END IF
    1689         288 :    END SUBROUTINE reftraj
    1690             : 
    1691             : ! **************************************************************************************************
    1692             : !> \brief nph_uniaxial integrator (non-Hamiltonian version)
    1693             : !>      for particle positions & momenta undergoing
    1694             : !>      uniaxial stress ( in x-direction of orthorhombic cell)
    1695             : !>      due to a shock compression:
    1696             : !>      Reed et. al. Physical Review Letters 90, 235503 (2003).
    1697             : !> \param md_env ...
    1698             : !> \par History
    1699             : !>      none
    1700             : !> \author CJM
    1701             : ! **************************************************************************************************
    1702          80 :    SUBROUTINE nph_uniaxial(md_env)
    1703             : 
    1704             :       TYPE(md_environment_type), POINTER                 :: md_env
    1705             : 
    1706             :       REAL(dp), PARAMETER                                :: e2 = 1._dp/6._dp, e4 = e2/20._dp, &
    1707             :                                                             e6 = e4/42._dp, e8 = e6/72._dp
    1708             : 
    1709             :       INTEGER                                            :: iroll, nparticle, nparticle_kind, nshell
    1710             :       INTEGER, POINTER                                   :: itimes
    1711             :       LOGICAL                                            :: first, first_time, shell_adiabatic, &
    1712             :                                                             shell_present
    1713             :       REAL(KIND=dp)                                      :: dt, infree, kin, roll_tol, roll_tol_thrs
    1714             :       REAL(KIND=dp), DIMENSION(3)                        :: vector_r, vector_v
    1715             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_kin
    1716             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1717          40 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1718             :       TYPE(cell_type), POINTER                           :: cell
    1719             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1720             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    1721             :       TYPE(force_env_type), POINTER                      :: force_env
    1722             :       TYPE(global_constraint_type), POINTER              :: gci
    1723             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1724          40 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1725             :       TYPE(molecule_list_type), POINTER                  :: molecules
    1726          40 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1727             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1728          40 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
    1729             :       TYPE(old_variables_type), POINTER                  :: old
    1730             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1731             :                                                             shell_particles
    1732          40 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
    1733          40 :                                                             shell_particle_set
    1734             :       TYPE(simpar_type), POINTER                         :: simpar
    1735             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1736             :       TYPE(virial_type), POINTER                         :: virial
    1737             : 
    1738          40 :       NULLIFY (gci, force_env)
    1739          40 :       NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles)
    1740          40 :       NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt)
    1741          40 :       NULLIFY (core_particles, particles, shell_particles, tmp, old)
    1742          40 :       NULLIFY (core_particle_set, particle_set, shell_particle_set)
    1743          40 :       NULLIFY (simpar, virial, itimes)
    1744             : 
    1745             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, npt=npt, &
    1746          40 :                       first_time=first_time, para_env=para_env, itimes=itimes)
    1747          40 :       dt = simpar%dt
    1748          40 :       infree = 1.0_dp/REAL(simpar%nfree, dp)
    1749             : 
    1750          40 :       CALL force_env_get(force_env, subsys=subsys, cell=cell)
    1751             : 
    1752             :       ! Do some checks on coordinates and box
    1753          40 :       CALL apply_qmmm_walls_reflective(force_env)
    1754             : 
    1755             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    1756             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, gci=gci, &
    1757          40 :                          molecule_kinds=molecule_kinds, virial=virial)
    1758             : 
    1759          40 :       nparticle_kind = atomic_kinds%n_els
    1760          40 :       atomic_kind_set => atomic_kinds%els
    1761          40 :       molecule_kind_set => molecule_kinds%els
    1762             : 
    1763          40 :       nparticle = particles%n_els
    1764          40 :       particle_set => particles%els
    1765          40 :       molecule_set => molecules%els
    1766             : 
    1767          40 :       IF (first_time) THEN
    1768             :          CALL virial_evaluate(atomic_kind_set, particle_set, &
    1769           4 :                               local_particles, virial, para_env)
    1770             :       END IF
    1771             : 
    1772             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    1773          40 :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic)
    1774             : 
    1775             :       ! Allocate work storage for positions and velocities
    1776          40 :       CALL allocate_old(old, particle_set, npt)
    1777             : 
    1778          40 :       IF (shell_present) THEN
    1779             :          CALL cp_subsys_get(subsys=subsys, &
    1780           0 :                             shell_particles=shell_particles, core_particles=core_particles)
    1781           0 :          shell_particle_set => shell_particles%els
    1782           0 :          nshell = SIZE(shell_particles%els)
    1783           0 :          IF (shell_adiabatic) THEN
    1784           0 :             core_particle_set => core_particles%els
    1785             :          END IF
    1786             :       END IF
    1787             : 
    1788          40 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
    1789             : 
    1790          40 :       IF (simpar%constraint) THEN
    1791             :          ! Possibly update the target values
    1792             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
    1793           0 :                                    molecule_kind_set, dt, force_env%root_section)
    1794             :       END IF
    1795             : 
    1796             :       ! setting up for ROLL: saving old variables
    1797          40 :       IF (simpar%constraint) THEN
    1798           0 :          roll_tol_thrs = simpar%roll_tol
    1799           0 :          iroll = 1
    1800           0 :          CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F')
    1801             :          CALL getold(gci, local_molecules, molecule_set, &
    1802           0 :                      molecule_kind_set, particle_set, cell)
    1803             :       ELSE
    1804             :          roll_tol_thrs = EPSILON(0.0_dp)
    1805             :       END IF
    1806          40 :       roll_tol = -roll_tol_thrs
    1807             : 
    1808          80 :       SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP
    1809             : 
    1810          40 :          IF (simpar%constraint) THEN
    1811           0 :             CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B')
    1812             :          END IF
    1813             :          CALL update_pv(gci, simpar, atomic_kind_set, particle_set, &
    1814             :                         local_molecules, molecule_set, molecule_kind_set, &
    1815          40 :                         local_particles, kin, pv_kin, virial, para_env)
    1816          40 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    1817             : 
    1818             :          tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* &
    1819          40 :                         (0.5_dp*npt(1, 1)%v*dt)
    1820             :          tmp%poly_r(1) = 1._dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + &
    1821          40 :                          e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4
    1822          40 :          tmp%poly_r(2) = 1.0_dp
    1823          40 :          tmp%poly_r(3) = 1.0_dp
    1824             : 
    1825             :          tmp%arg_v(1) = (0.25_dp*npt(1, 1)%v*dt* &
    1826             :                          (1._dp + infree))*(0.25_dp*npt(1, 1)%v* &
    1827          40 :                                             dt*(1._dp + infree))
    1828             :          tmp%arg_v(2) = (0.25_dp*npt(1, 1)%v*dt*infree)* &
    1829          40 :                         (0.25_dp*npt(1, 1)%v*dt*infree)
    1830             :          tmp%poly_v(1) = 1._dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + &
    1831          40 :                          e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4
    1832             :          tmp%poly_v(2) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
    1833          40 :                          e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
    1834             :          tmp%poly_v(3) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
    1835          40 :                          e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
    1836             : 
    1837          40 :          tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v)
    1838          40 :          tmp%scale_r(2) = 1.0_dp
    1839          40 :          tmp%scale_r(3) = 1.0_dp
    1840             : 
    1841             :          tmp%scale_v(1) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
    1842          40 :                               (1._dp + infree))
    1843          40 :          tmp%scale_v(2) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
    1844          40 :          tmp%scale_v(3) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
    1845             : 
    1846             :          ! first half of velocity verlet
    1847             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1848             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1849          40 :                        shell_adiabatic, dt)
    1850             : 
    1851          40 :          IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
    1852             :                                                         atomic_kind_set, local_particles, particle_set, core_particle_set, &
    1853           0 :                                                         shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt)
    1854             : 
    1855          40 :          roll_tol = 0._dp
    1856          40 :          vector_r(:) = 0._dp
    1857         160 :          vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:)
    1858          40 :          vector_r(1) = tmp%scale_r(1)*tmp%poly_r(1)
    1859             : 
    1860          40 :          IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, &
    1861             :                                                       molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, &
    1862             :                                                         roll_tol, iroll, vector_r, vector_v, para_env, cell=cell, &
    1863          40 :                                                         local_particles=local_particles)
    1864             :       END DO SR
    1865             : 
    1866             :       ! Update h_mat
    1867          40 :       cell%hmat(1, 1) = cell%hmat(1, 1)*tmp%scale_r(1)*tmp%scale_r(1)
    1868             : 
    1869             :       ! Update the cell
    1870          40 :       CALL init_cell(cell)
    1871             : 
    1872             :       ! Broadcast the new particle positions and deallocate the pos component of temporary
    1873             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1874          40 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
    1875             : 
    1876             :       ! Update forces (and stress)
    1877          40 :       CALL force_env_calc_energy_force(force_env)
    1878             : 
    1879             :       ! Metadynamics
    1880          40 :       CALL metadyn_integrator(force_env, itimes, tmp%vel)
    1881             : 
    1882             :       ! Velocity Verlet (second part)
    1883             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    1884             :                      core_particle_set, shell_particle_set, nparticle_kind, &
    1885          40 :                      shell_adiabatic, dt)
    1886             : 
    1887          40 :       IF (simpar%constraint) THEN
    1888           0 :          roll_tol_thrs = simpar%roll_tol
    1889           0 :          first = .TRUE.
    1890           0 :          iroll = 1
    1891           0 :          CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F')
    1892             :       ELSE
    1893             :          roll_tol_thrs = EPSILON(0.0_dp)
    1894             :       END IF
    1895          40 :       roll_tol = -roll_tol_thrs
    1896             : 
    1897          80 :       RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP
    1898          40 :          roll_tol = 0._dp
    1899          40 :          IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, &
    1900             :                                                        particle_set, local_particles, molecule_kind_set, molecule_set, &
    1901             :                                                        local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, &
    1902           0 :                                                        roll_tol, iroll, infree, first, para_env)
    1903             : 
    1904             :          CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
    1905             :                         local_molecules, molecule_set, molecule_kind_set, &
    1906          40 :                         local_particles, kin, pv_kin, virial, para_env)
    1907          80 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    1908             :       END DO RR
    1909             : 
    1910          40 :       IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
    1911             : 
    1912             :       ! Broadcast the new particle velocities and deallocate the temporary
    1913             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    1914          40 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
    1915             : 
    1916             :       ! Update constraint virial
    1917          40 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
    1918           0 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
    1919             : 
    1920             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
    1921          40 :                            local_particles, virial, para_env)
    1922             : 
    1923             :       ! Deallocate old variables
    1924          40 :       CALL deallocate_old(old)
    1925             : 
    1926          40 :       IF (first_time) THEN
    1927           4 :          first_time = .FALSE.
    1928           4 :          CALL set_md_env(md_env, first_time=first_time)
    1929             :       END IF
    1930             : 
    1931          40 :    END SUBROUTINE nph_uniaxial
    1932             : 
    1933             : ! **************************************************************************************************
    1934             : !> \brief nph_uniaxial integrator (non-Hamiltonian version)
    1935             : !>      for particle positions & momenta undergoing
    1936             : !>      uniaxial stress ( in x-direction of orthorhombic cell)
    1937             : !>      due to a shock compression:
    1938             : !>      Reed et. al. Physical Review Letters 90, 235503 (2003).
    1939             : !>      Added damping (e.g. thermostat to barostat)
    1940             : !> \param md_env ...
    1941             : !> \par History
    1942             : !>      none
    1943             : !> \author CJM
    1944             : ! **************************************************************************************************
    1945          40 :    SUBROUTINE nph_uniaxial_damped(md_env)
    1946             : 
    1947             :       TYPE(md_environment_type), POINTER                 :: md_env
    1948             : 
    1949             :       REAL(dp), PARAMETER                                :: e2 = 1._dp/6._dp, e4 = e2/20._dp, &
    1950             :                                                             e6 = e4/42._dp, e8 = e6/72._dp
    1951             : 
    1952             :       INTEGER                                            :: iroll, nparticle, nparticle_kind, nshell
    1953             :       INTEGER, POINTER                                   :: itimes
    1954             :       LOGICAL                                            :: first, first_time, shell_adiabatic, &
    1955             :                                                             shell_present
    1956             :       REAL(KIND=dp)                                      :: aa, aax, dt, gamma1, infree, kin, &
    1957             :                                                             roll_tol, roll_tol_thrs
    1958             :       REAL(KIND=dp), DIMENSION(3)                        :: vector_r, vector_v
    1959             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_kin
    1960             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1961          20 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1962             :       TYPE(cell_type), POINTER                           :: cell
    1963             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1964             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    1965             :       TYPE(force_env_type), POINTER                      :: force_env
    1966             :       TYPE(global_constraint_type), POINTER              :: gci
    1967             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1968          20 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1969             :       TYPE(molecule_list_type), POINTER                  :: molecules
    1970          20 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1971             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1972          20 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
    1973             :       TYPE(old_variables_type), POINTER                  :: old
    1974             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1975             :                                                             shell_particles
    1976          20 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
    1977          20 :                                                             shell_particle_set
    1978             :       TYPE(simpar_type), POINTER                         :: simpar
    1979             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1980             :       TYPE(virial_type), POINTER                         :: virial
    1981             : 
    1982          20 :       NULLIFY (gci, force_env)
    1983          20 :       NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles)
    1984          20 :       NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt)
    1985          20 :       NULLIFY (core_particles, particles, shell_particles, tmp, old)
    1986          20 :       NULLIFY (core_particle_set, particle_set, shell_particle_set)
    1987          20 :       NULLIFY (simpar, virial, itimes)
    1988             : 
    1989             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, npt=npt, &
    1990          20 :                       first_time=first_time, para_env=para_env, itimes=itimes)
    1991          20 :       dt = simpar%dt
    1992          20 :       infree = 1.0_dp/REAL(simpar%nfree, dp)
    1993          20 :       gamma1 = simpar%gamma_nph
    1994             : 
    1995          20 :       CALL force_env_get(force_env, subsys=subsys, cell=cell)
    1996             : 
    1997             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    1998             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, gci=gci, &
    1999          20 :                          molecule_kinds=molecule_kinds, virial=virial)
    2000             : 
    2001          20 :       nparticle_kind = atomic_kinds%n_els
    2002          20 :       atomic_kind_set => atomic_kinds%els
    2003          20 :       molecule_kind_set => molecule_kinds%els
    2004             : 
    2005          20 :       nparticle = particles%n_els
    2006          20 :       particle_set => particles%els
    2007          20 :       molecule_set => molecules%els
    2008             : 
    2009          20 :       IF (first_time) THEN
    2010             :          CALL virial_evaluate(atomic_kind_set, particle_set, &
    2011           2 :                               local_particles, virial, para_env)
    2012             :       END IF
    2013             : 
    2014             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    2015          20 :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic)
    2016             : 
    2017             :       ! Allocate work storage for positions and velocities
    2018          20 :       CALL allocate_old(old, particle_set, npt)
    2019             : 
    2020          20 :       IF (shell_present) THEN
    2021             :          CALL cp_subsys_get(subsys=subsys, &
    2022           0 :                             shell_particles=shell_particles, core_particles=core_particles)
    2023           0 :          shell_particle_set => shell_particles%els
    2024           0 :          nshell = SIZE(shell_particles%els)
    2025           0 :          IF (shell_adiabatic) THEN
    2026           0 :             core_particle_set => core_particles%els
    2027             :          END IF
    2028             :       END IF
    2029             : 
    2030          20 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
    2031             : 
    2032             :       ! perform damping on velocities
    2033             :       CALL damp_v(molecule_kind_set, molecule_set, particle_set, local_molecules, &
    2034          20 :                   gamma1, npt(1, 1), dt, para_env)
    2035             : 
    2036          20 :       IF (simpar%constraint) THEN
    2037             :          ! Possibly update the target values
    2038             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2039           0 :                                    molecule_kind_set, dt, force_env%root_section)
    2040             :       END IF
    2041             : 
    2042             :       ! setting up for ROLL: saving old variables
    2043          20 :       IF (simpar%constraint) THEN
    2044           0 :          roll_tol_thrs = simpar%roll_tol
    2045           0 :          iroll = 1
    2046           0 :          CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F')
    2047             :          CALL getold(gci, local_molecules, molecule_set, &
    2048           0 :                      molecule_kind_set, particle_set, cell)
    2049             :       ELSE
    2050             :          roll_tol_thrs = EPSILON(0.0_dp)
    2051             :       END IF
    2052          20 :       roll_tol = -roll_tol_thrs
    2053             : 
    2054          40 :       SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP
    2055             : 
    2056             :          ! perform damping on the barostat momentum
    2057          20 :          CALL damp_veps(npt(1, 1), gamma1, dt)
    2058             : 
    2059          20 :          IF (simpar%constraint) THEN
    2060           0 :             CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B')
    2061             :          END IF
    2062             :          CALL update_pv(gci, simpar, atomic_kind_set, particle_set, &
    2063             :                         local_molecules, molecule_set, molecule_kind_set, &
    2064          20 :                         local_particles, kin, pv_kin, virial, para_env)
    2065          20 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    2066             : 
    2067             :          ! perform damping on the barostat momentum
    2068          20 :          CALL damp_veps(npt(1, 1), gamma1, dt)
    2069             : 
    2070             :          tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* &
    2071          20 :                         (0.5_dp*npt(1, 1)%v*dt)
    2072             :          tmp%poly_r(1) = 1._dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + &
    2073          20 :                          e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4
    2074             : 
    2075          20 :          aax = npt(1, 1)%v*(1.0_dp + infree)
    2076          20 :          tmp%arg_v(1) = (0.25_dp*dt*aax)*(0.25_dp*dt*aax)
    2077             :          tmp%poly_v(1) = 1._dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + &
    2078          20 :                          e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4
    2079             : 
    2080          20 :          aa = npt(1, 1)%v*infree
    2081          20 :          tmp%arg_v(2) = (0.25_dp*dt*aa)*(0.25_dp*dt*aa)
    2082             :          tmp%poly_v(2) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
    2083          20 :                          e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
    2084             :          tmp%poly_v(3) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + &
    2085          20 :                          e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4
    2086             : 
    2087          20 :          tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v)
    2088          20 :          tmp%scale_v(1) = EXP(-0.25_dp*dt*aax)
    2089          20 :          tmp%scale_v(2) = EXP(-0.25_dp*dt*aa)
    2090          20 :          tmp%scale_v(3) = EXP(-0.25_dp*dt*aa)
    2091             : 
    2092             :          ! first half of velocity verlet
    2093             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    2094             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    2095          20 :                        shell_adiabatic, dt)
    2096             : 
    2097          20 :          IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
    2098             :                                                         atomic_kind_set, local_particles, particle_set, core_particle_set, &
    2099           0 :                                                         shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt)
    2100             : 
    2101          20 :          roll_tol = 0._dp
    2102          20 :          vector_r(:) = 0._dp
    2103          80 :          vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:)
    2104          20 :          vector_r(1) = tmp%scale_r(1)*tmp%poly_r(1)
    2105             : 
    2106          20 :          IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, &
    2107             :                                                       molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, &
    2108             :                                                         roll_tol, iroll, vector_r, vector_v, para_env, cell=cell, &
    2109          20 :                                                         local_particles=local_particles)
    2110             :       END DO SR
    2111             : 
    2112             :       ! Update h_mat
    2113          20 :       cell%hmat(1, 1) = cell%hmat(1, 1)*tmp%scale_r(1)*tmp%scale_r(1)
    2114             : 
    2115             :       ! Update the inverse
    2116          20 :       CALL init_cell(cell)
    2117             : 
    2118             :       ! Broadcast the new particle positions and deallocate the pos components of temporary
    2119             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    2120          20 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
    2121             : 
    2122             :       ! Update forces
    2123          20 :       CALL force_env_calc_energy_force(force_env)
    2124             : 
    2125             :       ! Metadynamics
    2126          20 :       CALL metadyn_integrator(force_env, itimes, tmp%vel)
    2127             : 
    2128             :       ! Velocity Verlet (second part)
    2129             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    2130             :                      core_particle_set, shell_particle_set, nparticle_kind, &
    2131          20 :                      shell_adiabatic, dt)
    2132             : 
    2133          20 :       IF (simpar%constraint) THEN
    2134           0 :          roll_tol_thrs = simpar%roll_tol
    2135           0 :          first = .TRUE.
    2136           0 :          iroll = 1
    2137           0 :          CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F')
    2138             :       ELSE
    2139             :          roll_tol_thrs = EPSILON(0.0_dp)
    2140             :       END IF
    2141          20 :       roll_tol = -roll_tol_thrs
    2142             : 
    2143          40 :       RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP
    2144          20 :          roll_tol = 0._dp
    2145          20 :          IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, &
    2146             :                                                   particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, &
    2147             :                                                  tmp%vel, dt, cell, npt, simpar, virial, vector_v, roll_tol, iroll, infree, first, &
    2148           0 :                                                        para_env)
    2149             :          ! perform damping on the barostat momentum
    2150          20 :          CALL damp_veps(npt(1, 1), gamma1, dt)
    2151             : 
    2152             :          CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
    2153             :                         local_molecules, molecule_set, molecule_kind_set, &
    2154          20 :                         local_particles, kin, pv_kin, virial, para_env)
    2155          20 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
    2156             : 
    2157             :          ! perform damping on the barostat momentum
    2158          20 :          CALL damp_veps(npt(1, 1), gamma1, dt)
    2159             : 
    2160             :       END DO RR
    2161             : 
    2162             :       ! perform damping on velocities
    2163             :       CALL damp_v(molecule_kind_set, molecule_set, particle_set, local_molecules, &
    2164          20 :                   tmp%vel, gamma1, npt(1, 1), dt, para_env)
    2165             : 
    2166          20 :       IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
    2167             : 
    2168             :       ! Broadcast the new particle velocities and deallocate temporary
    2169             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    2170          20 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
    2171             : 
    2172             :       ! Update constraint virial
    2173          20 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
    2174           0 :                                                 molecule_set, molecule_kind_set, particle_set, virial, para_env)
    2175             : 
    2176             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
    2177          20 :                            local_particles, virial, para_env)
    2178             : 
    2179             :       ! Deallocate old variables
    2180          20 :       CALL deallocate_old(old)
    2181             : 
    2182          20 :       IF (first_time) THEN
    2183           2 :          first_time = .FALSE.
    2184           2 :          CALL set_md_env(md_env, first_time=first_time)
    2185             :       END IF
    2186             : 
    2187          20 :    END SUBROUTINE nph_uniaxial_damped
    2188             : 
    2189             : ! **************************************************************************************************
    2190             : !> \brief Velocity Verlet integrator for the NPT ensemble with fully flexible cell
    2191             : !> \param md_env ...
    2192             : !> \param globenv ...
    2193             : !> \par History
    2194             : !>      none
    2195             : !> \author CJM
    2196             : ! **************************************************************************************************
    2197         896 :    SUBROUTINE npt_f(md_env, globenv)
    2198             : 
    2199             :       TYPE(md_environment_type), POINTER                 :: md_env
    2200             :       TYPE(global_environment_type), POINTER             :: globenv
    2201             : 
    2202             :       REAL(KIND=dp), PARAMETER                           :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
    2203             :                                                             e6 = e4/42.0_dp, e8 = e6/72.0_dp
    2204             : 
    2205             :       INTEGER                                            :: i, iroll, j, nparticle, nparticle_kind, &
    2206             :                                                             nshell
    2207             :       INTEGER, POINTER                                   :: itimes
    2208             :       LOGICAL                                            :: first, first_time, shell_adiabatic, &
    2209             :                                                             shell_check_distance, shell_present
    2210             :       REAL(KIND=dp)                                      :: dt, infree, kin, roll_tol, &
    2211             :                                                             roll_tol_thrs, trvg
    2212             :       REAL(KIND=dp), DIMENSION(3)                        :: vector_r, vector_v
    2213             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_kin, uh
    2214             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    2215         896 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2216             :       TYPE(barostat_type), POINTER                       :: barostat
    2217             :       TYPE(cell_type), POINTER                           :: cell
    2218             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    2219             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    2220             :       TYPE(force_env_type), POINTER                      :: force_env
    2221             :       TYPE(global_constraint_type), POINTER              :: gci
    2222             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    2223         896 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    2224             :       TYPE(molecule_list_type), POINTER                  :: molecules
    2225         896 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    2226             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2227         896 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
    2228             :       TYPE(old_variables_type), POINTER                  :: old
    2229             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    2230             :                                                             shell_particles
    2231         896 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
    2232         896 :                                                             shell_particle_set
    2233             :       TYPE(simpar_type), POINTER                         :: simpar
    2234             :       TYPE(thermostat_type), POINTER                     :: thermostat_baro, thermostat_part, &
    2235             :                                                             thermostat_shell
    2236             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    2237             :       TYPE(virial_type), POINTER                         :: virial
    2238             : 
    2239         896 :       NULLIFY (gci, thermostat_baro, thermostat_part, thermostat_shell, force_env)
    2240         896 :       NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles)
    2241         896 :       NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt, barostat)
    2242         896 :       NULLIFY (core_particles, particles, shell_particles, tmp, old)
    2243         896 :       NULLIFY (core_particle_set, particle_set, shell_particle_set)
    2244         896 :       NULLIFY (simpar, virial, itimes)
    2245             : 
    2246             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
    2247             :                       thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
    2248             :                       thermostat_shell=thermostat_shell, npt=npt, first_time=first_time, &
    2249         896 :                       para_env=para_env, barostat=barostat, itimes=itimes)
    2250         896 :       dt = simpar%dt
    2251         896 :       infree = 1.0_dp/REAL(simpar%nfree, KIND=dp)
    2252             : 
    2253         896 :       CALL force_env_get(force_env, subsys=subsys, cell=cell)
    2254             : 
    2255             :       ! Do some checks on coordinates and box
    2256         896 :       CALL apply_qmmm_walls_reflective(force_env)
    2257             : 
    2258             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    2259             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
    2260         896 :                          gci=gci, molecule_kinds=molecule_kinds, virial=virial)
    2261             : 
    2262         896 :       nparticle_kind = atomic_kinds%n_els
    2263         896 :       atomic_kind_set => atomic_kinds%els
    2264         896 :       molecule_kind_set => molecule_kinds%els
    2265             : 
    2266         896 :       nparticle = particles%n_els
    2267         896 :       particle_set => particles%els
    2268         896 :       molecule_set => molecules%els
    2269             : 
    2270             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
    2271             :                                shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
    2272         896 :                                shell_check_distance=shell_check_distance)
    2273             : 
    2274         896 :       IF (first_time) THEN
    2275             :          CALL virial_evaluate(atomic_kind_set, particle_set, &
    2276          58 :                               local_particles, virial, para_env)
    2277             :       END IF
    2278             : 
    2279             :       ! Allocate work storage for positions and velocities
    2280         896 :       CALL allocate_old(old, particle_set, npt)
    2281             : 
    2282         896 :       IF (shell_present) THEN
    2283             :          CALL cp_subsys_get(subsys=subsys, &
    2284         650 :                             shell_particles=shell_particles, core_particles=core_particles)
    2285         650 :          shell_particle_set => shell_particles%els
    2286         650 :          nshell = SIZE(shell_particles%els)
    2287         650 :          IF (shell_adiabatic) THEN
    2288         650 :             core_particle_set => core_particles%els
    2289             :          END IF
    2290             :       END IF
    2291             : 
    2292         896 :       CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
    2293             : 
    2294             :       ! Apply Thermostat to Barostat
    2295         896 :       CALL apply_thermostat_baro(thermostat_baro, npt, para_env)
    2296             : 
    2297             :       ! Apply Thermostat over the full set of particles
    2298         896 :       IF (simpar%ensemble /= npe_f_ensemble) THEN
    2299         656 :          IF (shell_adiabatic) THEN
    2300             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    2301             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    2302         410 :                                             shell_particle_set=shell_particle_set, core_particle_set=core_particle_set)
    2303             :          ELSE
    2304             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    2305         246 :                                             particle_set, local_molecules, local_particles, para_env)
    2306             :          END IF
    2307             :       END IF
    2308             : 
    2309             :       ! Apply Thermostat over the core-shell motion
    2310             :       CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    2311             :                                    local_particles, para_env, shell_particle_set=shell_particle_set, &
    2312         896 :                                    core_particle_set=core_particle_set)
    2313             : 
    2314         896 :       IF (simpar%constraint) THEN
    2315             :          ! Possibly update the target values
    2316             :          CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2317          10 :                                    molecule_kind_set, dt, force_env%root_section)
    2318             :       END IF
    2319             : 
    2320             :       ! setting up for ROLL: saving old variables
    2321         896 :       IF (simpar%constraint) THEN
    2322          10 :          roll_tol_thrs = simpar%roll_tol
    2323          10 :          iroll = 1
    2324          10 :          CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F')
    2325             :          CALL getold(gci, local_molecules, molecule_set, &
    2326          10 :                      molecule_kind_set, particle_set, cell)
    2327             :       ELSE
    2328             :          roll_tol_thrs = EPSILON(0.0_dp)
    2329             :       END IF
    2330         896 :       roll_tol = -roll_tol_thrs
    2331             : 
    2332        1802 :       SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP
    2333             : 
    2334         906 :          IF (simpar%constraint) THEN
    2335          20 :             CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B')
    2336             :          END IF
    2337             :          CALL update_pv(gci, simpar, atomic_kind_set, particle_set, &
    2338             :                         local_molecules, molecule_set, molecule_kind_set, &
    2339         906 :                         local_particles, kin, pv_kin, virial, para_env)
    2340             :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree, &
    2341         906 :                           virial_components=barostat%virial_components)
    2342             : 
    2343         906 :          trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
    2344             :          !
    2345             :          ! find eigenvalues and eigenvectors of npt ( :, : )%v
    2346             :          !
    2347             : 
    2348             :          CALL diagonalise(matrix=npt(:, :)%v, mysize=3, &
    2349       11778 :                           storageform="UPPER", eigenvalues=tmp%e_val, eigenvectors=tmp%u)
    2350             : 
    2351             :          tmp%arg_r(:) = 0.5_dp*tmp%e_val(:)*dt* &
    2352        3624 :                         0.5_dp*tmp%e_val(:)*dt
    2353             :          tmp%poly_r = 1.0_dp + e2*tmp%arg_r + e4*tmp%arg_r*tmp%arg_r + &
    2354        3624 :                       e6*tmp%arg_r**3 + e8*tmp%arg_r**4
    2355        3624 :          tmp%scale_r(:) = EXP(0.5_dp*dt*tmp%e_val(:))
    2356             : 
    2357             :          tmp%arg_v(:) = 0.25_dp*dt*(tmp%e_val(:) + trvg*infree)* &
    2358        3624 :                         0.25_dp*dt*(tmp%e_val(:) + trvg*infree)
    2359             :          tmp%poly_v = 1.0_dp + e2*tmp%arg_v + e4*tmp%arg_v*tmp%arg_v + &
    2360        3624 :                       e6*tmp%arg_v**3 + e8*tmp%arg_v**4
    2361        3624 :          tmp%scale_v(:) = EXP(-0.25_dp*dt*(tmp%e_val(:) + trvg*infree))
    2362             : 
    2363             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    2364             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    2365         906 :                        shell_adiabatic, dt, u=tmp%u)
    2366             : 
    2367         906 :          IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
    2368             :                                                         atomic_kind_set, local_particles, particle_set, core_particle_set, &
    2369         200 :                                                         shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt)
    2370             : 
    2371         906 :          roll_tol = 0.0_dp
    2372        3624 :          vector_r = tmp%scale_r*tmp%poly_r
    2373        3624 :          vector_v = tmp%scale_v*tmp%poly_v
    2374             : 
    2375         906 :          IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, &
    2376             :                                                         molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, &
    2377             :                                                         simpar, roll_tol, iroll, vector_r, vector_v, &
    2378             :                                                         para_env, u=tmp%u, cell=cell, &
    2379         916 :                                                         local_particles=local_particles)
    2380             :       END DO SR
    2381             : 
    2382             :       ! Update h_mat
    2383       35840 :       uh = MATMUL(TRANSPOSE(tmp%u), cell%hmat)
    2384             : 
    2385        3584 :       DO i = 1, 3
    2386       11648 :          DO j = 1, 3
    2387       10752 :             uh(i, j) = uh(i, j)*tmp%scale_r(i)*tmp%scale_r(i)
    2388             :          END DO
    2389             :       END DO
    2390             : 
    2391       57344 :       cell%hmat = MATMUL(tmp%u, uh)
    2392             :       ! Update the inverse
    2393         896 :       CALL init_cell(cell)
    2394             : 
    2395             :       ! Broadcast the new particle positions and deallocate the pos components of temporary
    2396             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    2397         896 :                               core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
    2398             : 
    2399         896 :       IF (shell_adiabatic .AND. shell_check_distance) THEN
    2400             :          CALL optimize_shell_core(force_env, particle_set, &
    2401         170 :                                   shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
    2402             :       END IF
    2403             : 
    2404             :       ! Update forces
    2405         896 :       CALL force_env_calc_energy_force(force_env)
    2406             : 
    2407             :       ! Metadynamics
    2408         896 :       CALL metadyn_integrator(force_env, itimes, tmp%vel)
    2409             : 
    2410             :       ! Velocity Verlet (second part)
    2411             :       CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    2412             :                      core_particle_set, shell_particle_set, nparticle_kind, &
    2413         896 :                      shell_adiabatic, dt, tmp%u)
    2414             : 
    2415         896 :       IF (simpar%constraint) THEN
    2416          10 :          roll_tol_thrs = simpar%roll_tol
    2417          10 :          first = .TRUE.
    2418          10 :          iroll = 1
    2419          10 :          CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F')
    2420             :       ELSE
    2421             :          roll_tol_thrs = EPSILON(0.0_dp)
    2422             :       END IF
    2423         896 :       roll_tol = -roll_tol_thrs
    2424             : 
    2425        1802 :       RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP
    2426         906 :          roll_tol = 0.0_dp
    2427         906 :          IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, &
    2428             :                                                        particle_set, local_particles, molecule_kind_set, molecule_set, &
    2429             :                                                        local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, &
    2430          20 :                                                        roll_tol, iroll, infree, first, para_env, u=tmp%u)
    2431             : 
    2432             :          CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
    2433             :                         local_molecules, molecule_set, molecule_kind_set, &
    2434         906 :                         local_particles, kin, pv_kin, virial, para_env)
    2435             :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree, &
    2436        1802 :                           virial_components=barostat%virial_components)
    2437             :       END DO RR
    2438             : 
    2439             :       ! Apply Thermostat over the full set of particles
    2440         896 :       IF (simpar%ensemble /= npe_f_ensemble) THEN
    2441         656 :          IF (shell_adiabatic) THEN
    2442             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    2443             :                                         particle_set, local_molecules, local_particles, para_env, shell_adiabatic=shell_adiabatic, &
    2444         410 :                                             vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel)
    2445             : 
    2446             :          ELSE
    2447             :             CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
    2448         246 :                                             particle_set, local_molecules, local_particles, para_env, vel=tmp%vel)
    2449             :          END IF
    2450             :       END IF
    2451             : 
    2452             :       ! Apply Thermostat over the core-shell motion
    2453         896 :       IF (ASSOCIATED(thermostat_shell)) THEN
    2454             :          CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
    2455             :                                       local_particles, para_env, vel=tmp%vel, shell_vel=tmp%shell_vel, &
    2456         320 :                                       core_vel=tmp%core_vel)
    2457             :       END IF
    2458             : 
    2459             :       ! Apply Thermostat to Barostat
    2460         896 :       CALL apply_thermostat_baro(thermostat_baro, npt, para_env)
    2461             : 
    2462             :       ! Annealing of particle velocities is only possible when no thermostat is active
    2463         896 :       IF (simpar%ensemble == npe_f_ensemble .AND. simpar%annealing) THEN
    2464       30800 :          tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
    2465          80 :          IF (shell_adiabatic) THEN
    2466             :             CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
    2467          80 :                                   tmp%vel, tmp%shell_vel, tmp%core_vel)
    2468             :          END IF
    2469             :       END IF
    2470             :       ! Annealing of CELL velocities is only possible when no thermostat is active
    2471         896 :       IF (simpar%ensemble == npe_f_ensemble .AND. simpar%annealing_cell) THEN
    2472         520 :          npt(:, :)%v = npt(:, :)%v*simpar%f_annealing_cell
    2473             :       END IF
    2474             : 
    2475             :       ! Broadcast the new particle velocities and deallocate temporary
    2476             :       CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
    2477         896 :                               core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
    2478             : 
    2479             :       ! Update constraint virial
    2480         896 :       IF (simpar%constraint) &
    2481             :          CALL pv_constraint(gci, local_molecules, molecule_set, &
    2482          10 :                             molecule_kind_set, particle_set, virial, para_env)
    2483             : 
    2484             :       CALL virial_evaluate(atomic_kind_set, particle_set, &
    2485         896 :                            local_particles, virial, para_env)
    2486             : 
    2487             :       ! Deallocate old variables
    2488         896 :       CALL deallocate_old(old)
    2489             : 
    2490         896 :       IF (first_time) THEN
    2491          58 :          first_time = .FALSE.
    2492          58 :          CALL set_md_env(md_env, first_time=first_time)
    2493             :       END IF
    2494             : 
    2495        1792 :    END SUBROUTINE npt_f
    2496             : 
    2497             : ! **************************************************************************************************
    2498             : !> \brief RESPA integrator for nve ensemble for particle positions & momenta
    2499             : !> \param md_env ...
    2500             : !> \author FS
    2501             : ! **************************************************************************************************
    2502          14 :    SUBROUTINE nve_respa(md_env)
    2503             : 
    2504             :       TYPE(md_environment_type), POINTER                 :: md_env
    2505             : 
    2506             :       INTEGER                                            :: i_step, iparticle, iparticle_kind, &
    2507             :                                                             iparticle_local, n_time_steps, &
    2508             :                                                             nparticle, nparticle_kind, &
    2509             :                                                             nparticle_local
    2510             :       INTEGER, POINTER                                   :: itimes
    2511             :       REAL(KIND=dp)                                      :: dm, dt, mass
    2512             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pos, vel
    2513             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    2514          14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2515             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    2516             :       TYPE(cell_type), POINTER                           :: cell
    2517             :       TYPE(cp_subsys_type), POINTER                      :: subsys, subsys_respa
    2518             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
    2519             :       TYPE(force_env_type), POINTER                      :: force_env
    2520             :       TYPE(global_constraint_type), POINTER              :: gci
    2521             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    2522          14 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    2523             :       TYPE(molecule_list_type), POINTER                  :: molecules
    2524          14 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    2525             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2526             :       TYPE(particle_list_type), POINTER                  :: particles, particles_respa
    2527          14 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, particle_set_respa
    2528             :       TYPE(simpar_type), POINTER                         :: simpar
    2529             : 
    2530          14 :       NULLIFY (para_env, cell, subsys_respa, particles_respa, particle_set_respa, gci, force_env, atomic_kinds)
    2531          14 :       NULLIFY (atomic_kind_set, simpar, subsys, particles, particle_set)
    2532          14 :       NULLIFY (local_molecules, molecule_kinds, molecules, molecule_kind_set, local_particles, itimes)
    2533             :       CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
    2534          14 :                       para_env=para_env, itimes=itimes)
    2535          14 :       dt = simpar%dt
    2536             : 
    2537          14 :       n_time_steps = simpar%n_time_steps
    2538             : 
    2539          14 :       CALL force_env_get(force_env, subsys=subsys, cell=cell)
    2540          14 :       CALL force_env_get(force_env%sub_force_env(1)%force_env, subsys=subsys_respa)
    2541             : 
    2542             :       ! Do some checks on coordinates and box
    2543          14 :       CALL apply_qmmm_walls_reflective(force_env)
    2544             : 
    2545             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
    2546             :                          particles=particles, local_molecules=local_molecules, molecules=molecules, &
    2547          14 :                          gci=gci, molecule_kinds=molecule_kinds)
    2548             : 
    2549          14 :       CALL cp_subsys_get(subsys=subsys_respa, particles=particles_respa)
    2550          14 :       particle_set_respa => particles_respa%els
    2551             : 
    2552          14 :       nparticle_kind = atomic_kinds%n_els
    2553          14 :       atomic_kind_set => atomic_kinds%els
    2554          14 :       molecule_kind_set => molecule_kinds%els
    2555             : 
    2556          14 :       nparticle = particles%n_els
    2557          14 :       particle_set => particles%els
    2558          14 :       molecule_set => molecules%els
    2559             : 
    2560             :       ! Allocate work storage for positions and velocities
    2561          42 :       ALLOCATE (pos(3, nparticle))
    2562          42 :       ALLOCATE (vel(3, nparticle))
    2563       21590 :       vel(:, :) = 0.0_dp
    2564             : 
    2565          14 :       IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, &
    2566           0 :                                          molecule_kind_set, particle_set, cell)
    2567             : 
    2568             :       ! Multiple time step (first part)
    2569          58 :       DO iparticle_kind = 1, nparticle_kind
    2570          44 :          atomic_kind => atomic_kind_set(iparticle_kind)
    2571          44 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    2572          44 :          dm = 0.5_dp*dt/mass
    2573          44 :          nparticle_local = local_particles%n_el(iparticle_kind)
    2574        2755 :          DO iparticle_local = 1, nparticle_local
    2575        2697 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    2576             :             vel(:, iparticle) = particle_set(iparticle)%v(:) + &
    2577             :                                 dm*(particle_set(iparticle)%f(:) - &
    2578       10832 :                                     particle_set_respa(iparticle)%f(:))
    2579             :          END DO
    2580             :       END DO
    2581             : 
    2582             :       ! Velocity Verlet (first part)
    2583          84 :       DO i_step = 1, n_time_steps
    2584      107950 :          pos(:, :) = 0.0_dp
    2585         290 :          DO iparticle_kind = 1, nparticle_kind
    2586         220 :             atomic_kind => atomic_kind_set(iparticle_kind)
    2587         220 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    2588         220 :             dm = 0.5_dp*dt/(n_time_steps*mass)
    2589         220 :             nparticle_local = local_particles%n_el(iparticle_kind)
    2590       13775 :             DO iparticle_local = 1, nparticle_local
    2591       13485 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    2592             :                vel(:, iparticle) = vel(:, iparticle) + &
    2593       53940 :                                    dm*particle_set_respa(iparticle)%f(:)
    2594             :                pos(:, iparticle) = particle_set(iparticle)%r(:) + &
    2595       54160 :                                    (dt/n_time_steps)*vel(:, iparticle)
    2596             :             END DO
    2597             :          END DO
    2598             : 
    2599          70 :          IF (simpar%constraint) THEN
    2600             :             ! Possibly update the target values
    2601             :             CALL shake_update_targets(gci, local_molecules, molecule_set, &
    2602           0 :                                       molecule_kind_set, dt, force_env%root_section)
    2603             : 
    2604             :             CALL shake_control(gci, local_molecules, molecule_set, &
    2605             :                                molecule_kind_set, particle_set, pos, vel, dt, simpar%shake_tol, &
    2606             :                                simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, cell, &
    2607           0 :                                para_env, local_particles)
    2608             :          END IF
    2609             : 
    2610             :          ! Broadcast the new particle positions
    2611          70 :          CALL update_particle_set(particle_set, para_env, pos=pos)
    2612       27040 :          DO iparticle = 1, SIZE(particle_set)
    2613      215830 :             particle_set_respa(iparticle)%r = particle_set(iparticle)%r
    2614             :          END DO
    2615             : 
    2616             :          ! Update forces
    2617          70 :          CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env)
    2618             : 
    2619             :          ! Metadynamics
    2620          70 :          CALL metadyn_integrator(force_env, itimes, vel)
    2621             : 
    2622             :          ! Velocity Verlet (second part)
    2623         290 :          DO iparticle_kind = 1, nparticle_kind
    2624         220 :             atomic_kind => atomic_kind_set(iparticle_kind)
    2625         220 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    2626         220 :             dm = 0.5_dp*dt/(n_time_steps*mass)
    2627         220 :             nparticle_local = local_particles%n_el(iparticle_kind)
    2628       13775 :             DO iparticle_local = 1, nparticle_local
    2629       13485 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    2630       13485 :                vel(1, iparticle) = vel(1, iparticle) + dm*particle_set_respa(iparticle)%f(1)
    2631       13485 :                vel(2, iparticle) = vel(2, iparticle) + dm*particle_set_respa(iparticle)%f(2)
    2632       13705 :                vel(3, iparticle) = vel(3, iparticle) + dm*particle_set_respa(iparticle)%f(3)
    2633             :             END DO
    2634             :          END DO
    2635             : 
    2636          70 :          IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, &
    2637             :                                                     molecule_kind_set, particle_set, vel, dt, simpar%shake_tol, &
    2638             :                                                     simpar%info_constraint, simpar%lagrange_multipliers, &
    2639           0 :                                                     simpar%dump_lm, cell, para_env, local_particles)
    2640             : 
    2641          84 :          IF (simpar%annealing) vel(:, :) = vel(:, :)*simpar%f_annealing
    2642             :       END DO
    2643          14 :       DEALLOCATE (pos)
    2644             : 
    2645             :       ! Multiple time step (second part)
    2646             :       ! Compute forces for respa force_env
    2647          14 :       CALL force_env_calc_energy_force(force_env)
    2648             : 
    2649             :       ! Metadynamics
    2650          14 :       CALL metadyn_integrator(force_env, itimes, vel)
    2651             : 
    2652          58 :       DO iparticle_kind = 1, nparticle_kind
    2653          44 :          atomic_kind => atomic_kind_set(iparticle_kind)
    2654          44 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    2655          44 :          dm = 0.5_dp*dt/mass
    2656          44 :          nparticle_local = local_particles%n_el(iparticle_kind)
    2657        2755 :          DO iparticle_local = 1, nparticle_local
    2658        2697 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    2659        2697 :             vel(1, iparticle) = vel(1, iparticle) + dm*(particle_set(iparticle)%f(1) - particle_set_respa(iparticle)%f(1))
    2660        2697 :             vel(2, iparticle) = vel(2, iparticle) + dm*(particle_set(iparticle)%f(2) - particle_set_respa(iparticle)%f(2))
    2661        2741 :             vel(3, iparticle) = vel(3, iparticle) + dm*(particle_set(iparticle)%f(3) - particle_set_respa(iparticle)%f(3))
    2662             :          END DO
    2663             :       END DO
    2664             : 
    2665             :       ! Broadcast the new particle velocities
    2666          14 :       CALL update_particle_set(particle_set, para_env, vel=vel)
    2667             : 
    2668          14 :       DEALLOCATE (vel)
    2669             : 
    2670          14 :    END SUBROUTINE nve_respa
    2671             : 
    2672             : END MODULE integrator

Generated by: LCOV version 1.15