LCOV - code coverage report
Current view: top level - src/motion - integrator.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 875 1031 84.9 %
Date: 2023-03-30 11:55:16 Functions: 10 11 90.9 %

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

Generated by: LCOV version 1.15