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