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