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