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