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