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