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