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