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