Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Provides integrator utility routines for the integrators
10 : !> \par History
11 : !> Teodoro Laino [tlaino] 02.2008 - Splitting from integrator and creation
12 : !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
13 : !> (patch by Marcel Baer)
14 : ! **************************************************************************************************
15 : MODULE integrator_utils
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE barostat_types, ONLY: do_clv_x,&
19 : do_clv_xy,&
20 : do_clv_xyz,&
21 : do_clv_xz,&
22 : do_clv_y,&
23 : do_clv_yz,&
24 : do_clv_z
25 : USE cell_types, ONLY: cell_type
26 : USE constraint, ONLY: rattle_roll_control
27 : USE constraint_util, ONLY: pv_constraint
28 : USE distribution_1d_types, ONLY: distribution_1d_type
29 : USE extended_system_types, ONLY: debug_isotropic_limit,&
30 : debug_uniaxial_limit,&
31 : npt_info_type
32 : USE input_constants, ONLY: &
33 : isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
34 : nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, nve_ensemble, &
35 : nvt_ensemble
36 : USE kinds, ONLY: dp
37 : USE md_environment_types, ONLY: get_md_env,&
38 : md_environment_type
39 : USE message_passing, ONLY: mp_comm_type,&
40 : mp_para_env_type
41 : USE molecule_kind_types, ONLY: local_fixd_constraint_type,&
42 : molecule_kind_type
43 : USE molecule_types, ONLY: get_molecule,&
44 : global_constraint_type,&
45 : molecule_type
46 : USE particle_types, ONLY: particle_type,&
47 : update_particle_set
48 : USE shell_potential_types, ONLY: shell_kind_type
49 : USE simpar_types, ONLY: simpar_type
50 : USE thermostat_types, ONLY: set_thermostats,&
51 : thermostats_type
52 : USE virial_types, ONLY: virial_type
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator_utils'
60 :
61 : ! **************************************************************************************************
62 : TYPE old_variables_type
63 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: v
64 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: r
65 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: eps
66 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: veps
67 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: h
68 : END TYPE old_variables_type
69 :
70 : ! **************************************************************************************************
71 : TYPE tmp_variables_type
72 : INTEGER, POINTER :: itimes
73 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: pos
74 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: vel
75 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: shell_pos
76 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: shell_vel
77 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: core_pos
78 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: core_vel
79 : REAL(KIND=dp) :: max_vel, max_vel_sc
80 : REAL(KIND=dp) :: max_dvel, max_dvel_sc
81 : REAL(KIND=dp) :: max_dr, max_dsc
82 : REAL(KIND=dp) :: arg_r(3), arg_v(3), u(3, 3), e_val(3), s, ds
83 : REAL(KIND=dp), DIMENSION(3) :: poly_r, &
84 : poly_v, scale_r, scale_v
85 : END TYPE tmp_variables_type
86 :
87 : INTERFACE set
88 : MODULE PROCEDURE set_particle_set, set_vel
89 : END INTERFACE
90 :
91 : INTERFACE update_pv
92 : MODULE PROCEDURE update_pv_particle_set, update_pv_velocity
93 : END INTERFACE
94 :
95 : INTERFACE damp_v
96 : MODULE PROCEDURE damp_v_particle_set, damp_v_velocity
97 : END INTERFACE
98 :
99 : PUBLIC :: old_variables_type, tmp_variables_type, allocate_old, deallocate_old, &
100 : allocate_tmp, update_dealloc_tmp, damp_v, set, update_pv, get_s_ds, &
101 : rattle_roll_setup, damp_veps, update_veps, vv_first, &
102 : vv_second, variable_timestep
103 :
104 : CONTAINS
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param old ...
109 : !> \param particle_set ...
110 : !> \param npt ...
111 : !> \par History
112 : !> none
113 : !> \author CJM
114 : ! **************************************************************************************************
115 2520 : SUBROUTINE allocate_old(old, particle_set, npt)
116 : TYPE(old_variables_type), POINTER :: old
117 : TYPE(particle_type), POINTER :: particle_set(:)
118 : TYPE(npt_info_type), POINTER :: npt(:, :)
119 :
120 : INTEGER :: idim, jdim, natoms
121 :
122 2520 : natoms = SIZE(particle_set)
123 2520 : idim = SIZE(npt, 1)
124 2520 : jdim = SIZE(npt, 2)
125 2520 : CPASSERT(.NOT. ASSOCIATED(old))
126 2520 : ALLOCATE (old)
127 :
128 7560 : ALLOCATE (old%v(natoms, 3))
129 2163108 : old%v = 0.0_dp
130 7560 : ALLOCATE (old%r(natoms, 3))
131 2163108 : old%r = 0.0_dp
132 10080 : ALLOCATE (old%eps(idim, jdim))
133 16520 : old%eps = 0.0_dp
134 10080 : ALLOCATE (old%veps(idim, jdim))
135 16520 : old%veps = 0.0_dp
136 2520 : ALLOCATE (old%h(3, 3))
137 32760 : old%h = 0.0_dp
138 :
139 2520 : END SUBROUTINE allocate_old
140 :
141 : ! **************************************************************************************************
142 : !> \brief ...
143 : !> \param old ...
144 : !> \par History
145 : !> none
146 : !> \author CJM
147 : ! **************************************************************************************************
148 2520 : SUBROUTINE deallocate_old(old)
149 : TYPE(old_variables_type), POINTER :: old
150 :
151 2520 : IF (ASSOCIATED(old)) THEN
152 2520 : IF (ASSOCIATED(old%v)) THEN
153 2520 : DEALLOCATE (old%v)
154 : END IF
155 2520 : IF (ASSOCIATED(old%r)) THEN
156 2520 : DEALLOCATE (old%r)
157 : END IF
158 2520 : IF (ASSOCIATED(old%eps)) THEN
159 2520 : DEALLOCATE (old%eps)
160 : END IF
161 2520 : IF (ASSOCIATED(old%veps)) THEN
162 2520 : DEALLOCATE (old%veps)
163 : END IF
164 2520 : IF (ASSOCIATED(old%h)) THEN
165 2520 : DEALLOCATE (old%h)
166 : END IF
167 2520 : DEALLOCATE (old)
168 : END IF
169 :
170 2520 : END SUBROUTINE deallocate_old
171 :
172 : ! **************************************************************************************************
173 : !> \brief allocate temporary variables to store positions and velocities
174 : !> used by the velocity-verlet integrator
175 : !> \param md_env ...
176 : !> \param tmp ...
177 : !> \param nparticle ...
178 : !> \param nshell ...
179 : !> \param shell_adiabatic ...
180 : !> \par History
181 : !> none
182 : !> \author MI (February 2008)
183 : ! **************************************************************************************************
184 42611 : SUBROUTINE allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
185 : TYPE(md_environment_type), POINTER :: md_env
186 : TYPE(tmp_variables_type), POINTER :: tmp
187 : INTEGER, INTENT(IN) :: nparticle, nshell
188 : LOGICAL, INTENT(IN) :: shell_adiabatic
189 :
190 42611 : CPASSERT(.NOT. ASSOCIATED(tmp))
191 42611 : ALLOCATE (tmp)
192 :
193 : ! Nullify Components
194 42611 : NULLIFY (tmp%pos)
195 42611 : NULLIFY (tmp%vel)
196 42611 : NULLIFY (tmp%shell_pos)
197 42611 : NULLIFY (tmp%shell_vel)
198 42611 : NULLIFY (tmp%core_pos)
199 42611 : NULLIFY (tmp%core_vel)
200 42611 : NULLIFY (tmp%itimes)
201 :
202 : ! *** Allocate work storage for positions and velocities ***
203 127833 : ALLOCATE (tmp%pos(3, nparticle))
204 :
205 127833 : ALLOCATE (tmp%vel(3, nparticle))
206 :
207 23778595 : tmp%pos(:, :) = 0.0_dp
208 23778595 : tmp%vel(:, :) = 0.0_dp
209 :
210 42611 : IF (shell_adiabatic) THEN
211 : ! *** Allocate work storage for positions and velocities ***
212 6870 : ALLOCATE (tmp%shell_pos(3, nshell))
213 :
214 6870 : ALLOCATE (tmp%core_pos(3, nshell))
215 :
216 6870 : ALLOCATE (tmp%shell_vel(3, nshell))
217 :
218 6870 : ALLOCATE (tmp%core_vel(3, nshell))
219 :
220 906050 : tmp%shell_pos(:, :) = 0.0_dp
221 906050 : tmp%shell_vel(:, :) = 0.0_dp
222 906050 : tmp%core_pos(:, :) = 0.0_dp
223 906050 : tmp%core_vel(:, :) = 0.0_dp
224 : END IF
225 :
226 170444 : tmp%arg_r = 0.0_dp
227 170444 : tmp%arg_v = 0.0_dp
228 553943 : tmp%u = 0.0_dp
229 170444 : tmp%e_val = 0.0_dp
230 42611 : tmp%max_vel = 0.0_dp
231 42611 : tmp%max_vel_sc = 0.0_dp
232 42611 : tmp%max_dvel = 0.0_dp
233 42611 : tmp%max_dvel_sc = 0.0_dp
234 170444 : tmp%scale_r = 1.0_dp
235 170444 : tmp%scale_v = 1.0_dp
236 170444 : tmp%poly_r = 1.0_dp
237 170444 : tmp%poly_v = 1.0_dp
238 :
239 42611 : CALL get_md_env(md_env=md_env, itimes=tmp%itimes)
240 :
241 42611 : END SUBROUTINE allocate_tmp
242 :
243 : ! **************************************************************************************************
244 : !> \brief update positions and deallocate temporary variable
245 : !> \param tmp ...
246 : !> \param particle_set ...
247 : !> \param shell_particle_set ...
248 : !> \param core_particle_set ...
249 : !> \param para_env ...
250 : !> \param shell_adiabatic ...
251 : !> \param pos ...
252 : !> \param vel ...
253 : !> \param should_deall_vel ...
254 : !> \par History
255 : !> none
256 : !> \author MI (February 2008)
257 : ! **************************************************************************************************
258 86372 : SUBROUTINE update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
259 : core_particle_set, para_env, shell_adiabatic, pos, vel, &
260 : should_deall_vel)
261 :
262 : TYPE(tmp_variables_type), POINTER :: tmp
263 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, shell_particle_set, &
264 : core_particle_set
265 : TYPE(mp_para_env_type), POINTER :: para_env
266 : LOGICAL, INTENT(IN) :: shell_adiabatic
267 : LOGICAL, INTENT(IN), OPTIONAL :: pos, vel, should_deall_vel
268 :
269 : LOGICAL :: my_deall, my_pos, my_vel
270 :
271 86372 : CPASSERT(ASSOCIATED(tmp))
272 86372 : my_pos = .FALSE.
273 86372 : my_vel = .FALSE.
274 86372 : my_deall = .TRUE.
275 86372 : IF (PRESENT(pos)) my_pos = pos
276 86372 : IF (PRESENT(vel)) my_vel = vel
277 86372 : IF (PRESENT(should_deall_vel)) my_deall = should_deall_vel
278 :
279 : ! *** Broadcast the new particle positions ***
280 86372 : IF (my_pos) THEN
281 42611 : CALL update_particle_set(particle_set, para_env, pos=tmp%pos)
282 42611 : DEALLOCATE (tmp%pos)
283 :
284 42611 : IF (shell_adiabatic) THEN
285 2290 : CALL update_particle_set(shell_particle_set, para_env, pos=tmp%shell_pos)
286 2290 : CALL update_particle_set(core_particle_set, para_env, pos=tmp%core_pos)
287 2290 : DEALLOCATE (tmp%shell_pos)
288 2290 : DEALLOCATE (tmp%core_pos)
289 : END IF
290 : END IF
291 :
292 : ! *** Broadcast the new particle velocities ***
293 86372 : IF (my_vel) THEN
294 43761 : CALL update_particle_set(particle_set, para_env, vel=tmp%vel)
295 43761 : IF (shell_adiabatic) THEN
296 2290 : CALL update_particle_set(shell_particle_set, para_env, vel=tmp%shell_vel)
297 2290 : CALL update_particle_set(core_particle_set, para_env, vel=tmp%core_vel)
298 : END IF
299 43761 : IF (my_deall) THEN
300 42611 : DEALLOCATE (tmp%vel)
301 42611 : IF (shell_adiabatic) THEN
302 2290 : DEALLOCATE (tmp%shell_vel)
303 2290 : DEALLOCATE (tmp%core_vel)
304 : END IF
305 42611 : CPASSERT(.NOT. ASSOCIATED(tmp%pos))
306 42611 : CPASSERT(.NOT. ASSOCIATED(tmp%shell_pos))
307 42611 : CPASSERT(.NOT. ASSOCIATED(tmp%core_pos))
308 42611 : DEALLOCATE (tmp)
309 : END IF
310 : END IF
311 :
312 86372 : END SUBROUTINE update_dealloc_tmp
313 :
314 : ! **************************************************************************************************
315 : !> \brief ...
316 : !> \param tmp ...
317 : !> \param nparticle_kind ...
318 : !> \param atomic_kind_set ...
319 : !> \param local_particles ...
320 : !> \param particle_set ...
321 : !> \param dt ...
322 : !> \param para_env ...
323 : !> \param tmpv ...
324 : !> \par History
325 : !> - Created [2004-07]
326 : !> \author Joost VandeVondele
327 : ! **************************************************************************************************
328 12 : SUBROUTINE get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, &
329 : dt, para_env, tmpv)
330 : TYPE(tmp_variables_type), POINTER :: tmp
331 : INTEGER :: nparticle_kind
332 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
333 : TYPE(distribution_1d_type), POINTER :: local_particles
334 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
335 : REAL(KIND=dp) :: dt
336 : TYPE(mp_para_env_type), POINTER :: para_env
337 : LOGICAL, INTENT(IN), OPTIONAL :: tmpv
338 :
339 : INTEGER :: iparticle, iparticle_kind, &
340 : iparticle_local, nparticle_local
341 : LOGICAL :: my_tmpv
342 : REAL(KIND=dp) :: a, b, K, mass, rb
343 : TYPE(atomic_kind_type), POINTER :: atomic_kind
344 :
345 12 : my_tmpv = .FALSE.
346 12 : IF (PRESENT(tmpv)) my_tmpv = tmpv
347 :
348 12 : K = 0.0_dp
349 12 : a = 0.0_dp
350 12 : b = 0.0_dp
351 36 : DO iparticle_kind = 1, nparticle_kind
352 24 : atomic_kind => atomic_kind_set(iparticle_kind)
353 24 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
354 36 : IF (mass /= 0.0_dp) THEN
355 24 : nparticle_local = local_particles%n_el(iparticle_kind)
356 24 : IF (my_tmpv) THEN
357 21 : DO iparticle_local = 1, nparticle_local
358 9 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
359 36 : K = K + 0.5_dp*mass*DOT_PRODUCT(tmp%vel(:, iparticle), tmp%vel(:, iparticle))
360 36 : a = a + DOT_PRODUCT(tmp%vel(:, iparticle), particle_set(iparticle)%f(:))
361 48 : b = b + (1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:), particle_set(iparticle)%f(:))
362 : END DO
363 : ELSE
364 21 : DO iparticle_local = 1, nparticle_local
365 9 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
366 36 : K = K + 0.5_dp*mass*DOT_PRODUCT(particle_set(iparticle)%v(:), particle_set(iparticle)%v(:))
367 36 : a = a + DOT_PRODUCT(particle_set(iparticle)%v(:), particle_set(iparticle)%f(:))
368 48 : b = b + (1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:), particle_set(iparticle)%f(:))
369 : END DO
370 : END IF
371 : END IF
372 : END DO
373 12 : CALL para_env%sum(K)
374 12 : CALL para_env%sum(a)
375 12 : CALL para_env%sum(b)
376 12 : a = a/(2.0_dp*K)
377 12 : b = b/(2.0_dp*K)
378 12 : rb = SQRT(b)
379 12 : tmp%s = (a/b)*(COSH(dt*rb/2.0_dp) - 1) + SINH(dt*rb/2.0_dp)/rb
380 12 : tmp%ds = (a/b)*(SINH(dt*rb/2.0_dp)*rb) + COSH(dt*rb/2.0_dp)
381 :
382 12 : END SUBROUTINE get_s_ds
383 :
384 : ! **************************************************************************************************
385 : !> \brief ...
386 : !> \param old ...
387 : !> \param atomic_kind_set ...
388 : !> \param particle_set ...
389 : !> \param local_particles ...
390 : !> \param cell ...
391 : !> \param npt ...
392 : !> \param char ...
393 : !> \par History
394 : !> none
395 : !> \author CJM
396 : ! **************************************************************************************************
397 2504 : SUBROUTINE set_particle_set(old, atomic_kind_set, particle_set, local_particles, cell, npt, char)
398 : TYPE(old_variables_type), POINTER :: old
399 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
400 : TYPE(particle_type), POINTER :: particle_set(:)
401 : TYPE(distribution_1d_type), POINTER :: local_particles
402 : TYPE(cell_type), POINTER :: cell
403 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt
404 : CHARACTER(LEN=*), INTENT(IN) :: char
405 :
406 : INTEGER :: idim, iparticle, iparticle_kind, &
407 : iparticle_local, nparticle_kind, &
408 : nparticle_local
409 :
410 2504 : nparticle_kind = SIZE(atomic_kind_set)
411 : SELECT CASE (char)
412 : CASE ('F') ! forward assigning the old
413 2040 : DO iparticle_kind = 1, nparticle_kind
414 1362 : nparticle_local = local_particles%n_el(iparticle_kind)
415 124267 : DO iparticle_local = 1, nparticle_local
416 122227 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
417 490270 : DO idim = 1, 3
418 366681 : old%v(iparticle, idim) = particle_set(iparticle)%v(idim)
419 488908 : old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
420 : END DO
421 : END DO
422 : END DO
423 2134 : old%eps(:, :) = npt(:, :)%eps
424 2134 : old%veps(:, :) = npt(:, :)%v
425 16950 : old%h(:, :) = cell%hmat(:, :)
426 : CASE ('B') ! back assigning the original variables
427 5498 : DO iparticle_kind = 1, nparticle_kind
428 3672 : nparticle_local = local_particles%n_el(iparticle_kind)
429 361572 : DO iparticle_local = 1, nparticle_local
430 356074 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
431 1427968 : DO idim = 1, 3
432 1068222 : particle_set(iparticle)%v(idim) = old%v(iparticle, idim)
433 1424296 : particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
434 : END DO
435 : END DO
436 : END DO
437 5678 : npt(:, :)%eps = old%eps(:, :)
438 5678 : npt(:, :)%v = old%veps(:, :)
439 48154 : cell%hmat(:, :) = old%h(:, :)
440 : END SELECT
441 :
442 2504 : END SUBROUTINE set_particle_set
443 :
444 : ! **************************************************************************************************
445 : !> \brief ...
446 : !> \param old ...
447 : !> \param atomic_kind_set ...
448 : !> \param particle_set ...
449 : !> \param vel ...
450 : !> \param local_particles ...
451 : !> \param cell ...
452 : !> \param npt ...
453 : !> \param char ...
454 : !> \par History
455 : !> none
456 : !> \author CJM
457 : ! **************************************************************************************************
458 2472 : SUBROUTINE set_vel(old, atomic_kind_set, particle_set, vel, local_particles, cell, npt, char)
459 : TYPE(old_variables_type), POINTER :: old
460 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
461 : TYPE(particle_type), POINTER :: particle_set(:)
462 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
463 : TYPE(distribution_1d_type), POINTER :: local_particles
464 : TYPE(cell_type), POINTER :: cell
465 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt
466 : CHARACTER(LEN=*), INTENT(IN) :: char
467 :
468 : INTEGER :: idim, iparticle, iparticle_kind, &
469 : iparticle_local, nparticle_kind, &
470 : nparticle_local
471 :
472 2472 : nparticle_kind = SIZE(atomic_kind_set)
473 : SELECT CASE (char)
474 : CASE ('F') ! forward assigning the old
475 2040 : DO iparticle_kind = 1, nparticle_kind
476 1362 : nparticle_local = local_particles%n_el(iparticle_kind)
477 124267 : DO iparticle_local = 1, nparticle_local
478 122227 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
479 490270 : DO idim = 1, 3
480 366681 : old%v(iparticle, idim) = vel(idim, iparticle)
481 488908 : old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
482 : END DO
483 : END DO
484 : END DO
485 2134 : old%eps(:, :) = npt(:, :)%eps
486 2134 : old%veps(:, :) = npt(:, :)%v
487 16950 : old%h(:, :) = cell%hmat(:, :)
488 : CASE ('B') ! back assigning the original variables
489 5400 : DO iparticle_kind = 1, nparticle_kind
490 3606 : nparticle_local = local_particles%n_el(iparticle_kind)
491 317027 : DO iparticle_local = 1, nparticle_local
492 311627 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
493 1250114 : DO idim = 1, 3
494 934881 : vel(idim, iparticle) = old%v(iparticle, idim)
495 1246508 : particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
496 : END DO
497 : END DO
498 : END DO
499 5582 : npt(:, :)%eps = old%eps(:, :)
500 5582 : npt(:, :)%v = old%veps(:, :)
501 47322 : cell%hmat(:, :) = old%h(:, :)
502 : END SELECT
503 :
504 2472 : END SUBROUTINE set_vel
505 :
506 : ! **************************************************************************************************
507 : !> \brief overloaded routine provides damping for particles via nph_uniaxial_damped dynamics
508 : !> \param molecule_kind_set ...
509 : !> \param molecule_set ...
510 : !> \param particle_set ...
511 : !> \param local_molecules ...
512 : !> \param gamma1 ...
513 : !> \param npt ...
514 : !> \param dt ...
515 : !> \param group ...
516 : !> \par History
517 : !> none
518 : !> \author CJM
519 : ! **************************************************************************************************
520 20 : SUBROUTINE damp_v_particle_set(molecule_kind_set, molecule_set, &
521 : particle_set, local_molecules, gamma1, npt, dt, group)
522 :
523 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
524 : TYPE(molecule_type), POINTER :: molecule_set(:)
525 : TYPE(particle_type), POINTER :: particle_set(:)
526 : TYPE(distribution_1d_type), POINTER :: local_molecules
527 : REAL(KIND=dp), INTENT(IN) :: gamma1
528 : TYPE(npt_info_type), INTENT(IN) :: npt
529 : REAL(KIND=dp), INTENT(IN) :: dt
530 :
531 : CLASS(mp_comm_type), INTENT(IN) :: group
532 :
533 : INTEGER :: first_atom, ikind, imol, imol_local, &
534 : ipart, last_atom, nmol_local
535 : REAL(KIND=dp) :: alpha, ikin, kin, mass, scale
536 : TYPE(atomic_kind_type), POINTER :: atomic_kind
537 : TYPE(molecule_type), POINTER :: molecule
538 :
539 20 : kin = 0.0_dp
540 1420 : DO ikind = 1, SIZE(molecule_kind_set)
541 : ! Compute the total kinetic energy
542 1400 : nmol_local = local_molecules%n_el(ikind)
543 2120 : DO imol_local = 1, nmol_local
544 700 : imol = local_molecules%list(ikind)%array(imol_local)
545 700 : molecule => molecule_set(imol)
546 : CALL get_molecule(molecule, first_atom=first_atom, &
547 700 : last_atom=last_atom)
548 2800 : DO ipart = first_atom, last_atom
549 700 : atomic_kind => particle_set(ipart)%atomic_kind
550 700 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
551 : kin = kin + mass*particle_set(ipart)%v(1)* &
552 700 : particle_set(ipart)%v(1)
553 : kin = kin + mass*particle_set(ipart)%v(2)* &
554 700 : particle_set(ipart)%v(2)
555 : kin = kin + mass*particle_set(ipart)%v(3)* &
556 1400 : particle_set(ipart)%v(3)
557 : END DO
558 : END DO
559 : END DO
560 : !
561 20 : CALL group%sum(kin)
562 : !
563 20 : ikin = 1.0_dp/kin
564 20 : scale = 1.0_dp
565 20 : alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
566 20 : scale = scale*SQRT(1.0_dp + alpha*0.5_dp*dt)
567 : ! Scale
568 1420 : DO ikind = 1, SIZE(molecule_kind_set)
569 1400 : nmol_local = local_molecules%n_el(ikind)
570 2120 : DO imol_local = 1, nmol_local
571 700 : imol = local_molecules%list(ikind)%array(imol_local)
572 700 : molecule => molecule_set(imol)
573 : CALL get_molecule(molecule, first_atom=first_atom, &
574 700 : last_atom=last_atom)
575 2800 : DO ipart = first_atom, last_atom
576 700 : particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*scale
577 700 : particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*scale
578 1400 : particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*scale
579 : END DO
580 : END DO
581 : END DO
582 :
583 20 : END SUBROUTINE damp_v_particle_set
584 :
585 : ! **************************************************************************************************
586 : !> \brief overloaded provides damping for particles via nph_uniaxial_damped dynamics
587 : !> \param molecule_kind_set ...
588 : !> \param molecule_set ...
589 : !> \param particle_set ...
590 : !> \param local_molecules ...
591 : !> \param vel ...
592 : !> \param gamma1 ...
593 : !> \param npt ...
594 : !> \param dt ...
595 : !> \param group ...
596 : !> \par History
597 : !> none
598 : !> \author CJM
599 : ! **************************************************************************************************
600 20 : SUBROUTINE damp_v_velocity(molecule_kind_set, molecule_set, &
601 20 : particle_set, local_molecules, vel, gamma1, npt, dt, group)
602 :
603 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
604 : TYPE(molecule_type), POINTER :: molecule_set(:)
605 : TYPE(particle_type), POINTER :: particle_set(:)
606 : TYPE(distribution_1d_type), POINTER :: local_molecules
607 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
608 : REAL(KIND=dp), INTENT(IN) :: gamma1
609 : TYPE(npt_info_type), INTENT(IN) :: npt
610 : REAL(KIND=dp), INTENT(IN) :: dt
611 :
612 : CLASS(mp_comm_type), INTENT(IN) :: group
613 :
614 : INTEGER :: first_atom, ikind, imol, imol_local, &
615 : ipart, last_atom, nmol_local
616 : REAL(KIND=dp) :: alpha, ikin, kin, mass, scale
617 : TYPE(atomic_kind_type), POINTER :: atomic_kind
618 : TYPE(molecule_type), POINTER :: molecule
619 :
620 20 : kin = 0.0_dp
621 : ! Compute the total kinetic energy
622 1420 : DO ikind = 1, SIZE(molecule_kind_set)
623 1400 : nmol_local = local_molecules%n_el(ikind)
624 2120 : DO imol_local = 1, nmol_local
625 700 : imol = local_molecules%list(ikind)%array(imol_local)
626 700 : molecule => molecule_set(imol)
627 : CALL get_molecule(molecule, first_atom=first_atom, &
628 700 : last_atom=last_atom)
629 2800 : DO ipart = first_atom, last_atom
630 700 : atomic_kind => particle_set(ipart)%atomic_kind
631 700 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
632 700 : kin = kin + mass*vel(1, ipart)*vel(1, ipart)
633 700 : kin = kin + mass*vel(2, ipart)*vel(2, ipart)
634 1400 : kin = kin + mass*vel(3, ipart)*vel(3, ipart)
635 : END DO
636 : END DO
637 : END DO
638 : !
639 20 : CALL group%sum(kin)
640 : !
641 20 : ikin = 1.0_dp/kin
642 20 : scale = 1.0_dp
643 20 : alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
644 20 : scale = scale*SQRT(1.0_dp + alpha*0.5_dp*dt)
645 : ! Scale
646 1420 : DO ikind = 1, SIZE(molecule_kind_set)
647 1400 : nmol_local = local_molecules%n_el(ikind)
648 2120 : DO imol_local = 1, nmol_local
649 700 : imol = local_molecules%list(ikind)%array(imol_local)
650 700 : molecule => molecule_set(imol)
651 : CALL get_molecule(molecule, first_atom=first_atom, &
652 700 : last_atom=last_atom)
653 2800 : DO ipart = first_atom, last_atom
654 700 : vel(1, ipart) = vel(1, ipart)*scale
655 700 : vel(2, ipart) = vel(2, ipart)*scale
656 1400 : vel(3, ipart) = vel(3, ipart)*scale
657 : END DO
658 : END DO
659 : END DO
660 :
661 20 : END SUBROUTINE damp_v_velocity
662 :
663 : ! **************************************************************************************************
664 : !> \brief provides damping for barostat via nph_uniaxial_damped dynamics
665 : !> \param npt ...
666 : !> \param gamma1 ...
667 : !> \param dt ...
668 : !> \par History
669 : !> none
670 : !> \author CJM
671 : ! **************************************************************************************************
672 80 : ELEMENTAL SUBROUTINE damp_veps(npt, gamma1, dt)
673 :
674 : TYPE(npt_info_type), INTENT(INOUT) :: npt
675 : REAL(KIND=dp), INTENT(IN) :: gamma1, dt
676 :
677 : REAL(KIND=dp) :: scale
678 :
679 80 : scale = 1.0_dp
680 80 : scale = scale*EXP(-gamma1*0.25_dp*dt)
681 : ! Scale
682 80 : npt%v = npt%v*scale
683 :
684 80 : END SUBROUTINE damp_veps
685 :
686 : ! **************************************************************************************************
687 : !> \brief update veps using multiplier obtained from SHAKE
688 : !> \param old ...
689 : !> \param gci ...
690 : !> \param atomic_kind_set ...
691 : !> \param particle_set ...
692 : !> \param local_particles ...
693 : !> \param molecule_kind_set ...
694 : !> \param molecule_set ...
695 : !> \param local_molecules ...
696 : !> \param vel ...
697 : !> \param dt ...
698 : !> \param cell ...
699 : !> \param npt ...
700 : !> \param simpar ...
701 : !> \param virial ...
702 : !> \param vector_v ...
703 : !> \param roll_tol ...
704 : !> \param iroll ...
705 : !> \param infree ...
706 : !> \param first ...
707 : !> \param para_env ...
708 : !> \param u ...
709 : !> \par History
710 : !> none
711 : !> \author CJM
712 : ! **************************************************************************************************
713 1794 : SUBROUTINE rattle_roll_setup(old, gci, atomic_kind_set, particle_set, local_particles, &
714 1794 : molecule_kind_set, molecule_set, local_molecules, vel, dt, cell, npt, simpar, virial, &
715 1794 : vector_v, roll_tol, iroll, infree, first, para_env, u)
716 :
717 : TYPE(old_variables_type), POINTER :: old
718 : TYPE(global_constraint_type), POINTER :: gci
719 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
720 : TYPE(particle_type), POINTER :: particle_set(:)
721 : TYPE(distribution_1d_type), POINTER :: local_particles
722 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
723 : TYPE(molecule_type), POINTER :: molecule_set(:)
724 : TYPE(distribution_1d_type), POINTER :: local_molecules
725 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
726 : REAL(KIND=dp), INTENT(IN) :: dt
727 : TYPE(cell_type), POINTER :: cell
728 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt
729 : TYPE(simpar_type), INTENT(IN) :: simpar
730 : TYPE(virial_type), POINTER :: virial
731 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector_v
732 : REAL(KIND=dp), INTENT(OUT) :: roll_tol
733 : INTEGER, INTENT(INOUT) :: iroll
734 : REAL(KIND=dp), INTENT(IN) :: infree
735 : LOGICAL, INTENT(INOUT) :: first
736 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
737 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: u(:, :)
738 :
739 : REAL(KIND=dp) :: kin
740 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin
741 : TYPE(npt_info_type), DIMENSION(3, 3) :: npt_loc
742 :
743 1794 : IF (first) THEN
744 : CALL update_pv(gci, simpar, atomic_kind_set, vel, particle_set, &
745 : local_molecules, molecule_set, molecule_kind_set, &
746 678 : local_particles, kin, pv_kin, virial, para_env)
747 678 : CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
748 :
749 : END IF
750 1794 : first = .FALSE.
751 :
752 : ! assigning local variable
753 1794 : SELECT CASE (simpar%ensemble)
754 : CASE (npt_i_ensemble, npt_ia_ensemble)
755 23062 : npt_loc(:, :)%v = 0.0_dp
756 23062 : npt_loc(:, :)%mass = 0.0_dp
757 1774 : npt_loc(1, 1)%v = npt(1, 1)%v
758 1774 : npt_loc(2, 2)%v = npt(1, 1)%v
759 1774 : npt_loc(3, 3)%v = npt(1, 1)%v
760 : npt_loc(1, 1)%mass = npt(1, 1)%mass
761 : npt_loc(2, 2)%mass = npt(1, 1)%mass
762 1774 : npt_loc(3, 3)%mass = npt(1, 1)%mass
763 : CASE (npt_f_ensemble)
764 2034 : npt_loc = npt
765 : END SELECT
766 :
767 : ! resetting
768 1794 : CALL set(old, atomic_kind_set, particle_set, vel, local_particles, cell, npt, 'B')
769 : CALL rattle_roll_control(gci, local_molecules, molecule_set, molecule_kind_set, &
770 : particle_set, vel, dt, simpar, vector_v, npt_loc%v, roll_tol, iroll, &
771 46624 : para_env, u, cell, local_particles)
772 :
773 1794 : END SUBROUTINE rattle_roll_setup
774 :
775 : ! **************************************************************************************************
776 : !> \brief Overloaded routine to compute veps given the particles
777 : !> structure or a local copy of the velocity array
778 : !> \param gci ...
779 : !> \param simpar ...
780 : !> \param atomic_kind_set ...
781 : !> \param particle_set ...
782 : !> \param local_molecules ...
783 : !> \param molecule_set ...
784 : !> \param molecule_kind_set ...
785 : !> \param local_particles ...
786 : !> \param kin ...
787 : !> \param pv_kin ...
788 : !> \param virial ...
789 : !> \param int_group ...
790 : !> \par History
791 : !> none
792 : !> \author CJM
793 : ! **************************************************************************************************
794 3668 : SUBROUTINE update_pv_particle_set(gci, simpar, atomic_kind_set, particle_set, &
795 : local_molecules, molecule_set, molecule_kind_set, &
796 : local_particles, kin, pv_kin, virial, int_group)
797 : TYPE(global_constraint_type), POINTER :: gci
798 : TYPE(simpar_type), INTENT(IN) :: simpar
799 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
800 : TYPE(particle_type), POINTER :: particle_set(:)
801 : TYPE(distribution_1d_type), POINTER :: local_molecules
802 : TYPE(molecule_type), POINTER :: molecule_set(:)
803 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
804 : TYPE(distribution_1d_type), POINTER :: local_particles
805 : REAL(KIND=dp), INTENT(OUT) :: kin
806 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: pv_kin
807 : TYPE(virial_type), INTENT(INOUT) :: virial
808 :
809 : CLASS(mp_comm_type), INTENT(IN) :: int_group
810 :
811 : INTEGER :: i, iparticle, iparticle_kind, &
812 : iparticle_local, j, nparticle_local
813 : REAL(KIND=dp) :: mass
814 : TYPE(atomic_kind_type), POINTER :: atomic_kind
815 :
816 3668 : pv_kin = 0.0_dp
817 3668 : kin = 0.0_dp
818 11584 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
819 7916 : atomic_kind => atomic_kind_set(iparticle_kind)
820 7916 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
821 7916 : nparticle_local = local_particles%n_el(iparticle_kind)
822 604269 : DO iparticle_local = 1, nparticle_local
823 592685 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
824 2378656 : DO i = 1, 3
825 7112220 : DO j = 1, 3
826 : pv_kin(i, j) = pv_kin(i, j) + &
827 : mass*particle_set(iparticle)%v(i)* &
828 7112220 : particle_set(iparticle)%v(j)
829 : END DO
830 : kin = kin + mass*particle_set(iparticle)%v(i)* &
831 2370740 : particle_set(iparticle)%v(i)
832 : END DO
833 : END DO
834 : END DO
835 :
836 3668 : CALL int_group%sum(pv_kin)
837 3668 : CALL int_group%sum(kin)
838 :
839 : ! updating the constraint virial
840 3668 : IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
841 : molecule_set, &
842 : molecule_kind_set, &
843 : particle_set, virial, &
844 1826 : int_group)
845 :
846 3668 : END SUBROUTINE update_pv_particle_set
847 :
848 : ! **************************************************************************************************
849 : !> \brief Overloaded routine to compute kinetic virials given the particles
850 : !> structure or a local copy of the velocity array
851 : !> \param gci ...
852 : !> \param simpar ...
853 : !> \param atomic_kind_set ...
854 : !> \param vel ...
855 : !> \param particle_set ...
856 : !> \param local_molecules ...
857 : !> \param molecule_set ...
858 : !> \param molecule_kind_set ...
859 : !> \param local_particles ...
860 : !> \param kin ...
861 : !> \param pv_kin ...
862 : !> \param virial ...
863 : !> \param int_group ...
864 : !> \par History
865 : !> none
866 : !> \author CJM
867 : ! **************************************************************************************************
868 4314 : SUBROUTINE update_pv_velocity(gci, simpar, atomic_kind_set, vel, particle_set, &
869 : local_molecules, molecule_set, molecule_kind_set, &
870 4314 : local_particles, kin, pv_kin, virial, int_group)
871 : TYPE(global_constraint_type), POINTER :: gci
872 : TYPE(simpar_type), INTENT(IN) :: simpar
873 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
874 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
875 : TYPE(particle_type), POINTER :: particle_set(:)
876 : TYPE(distribution_1d_type), POINTER :: local_molecules
877 : TYPE(molecule_type), POINTER :: molecule_set(:)
878 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
879 : TYPE(distribution_1d_type), POINTER :: local_particles
880 : REAL(KIND=dp), INTENT(OUT) :: kin
881 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pv_kin
882 : TYPE(virial_type), INTENT(INOUT) :: virial
883 :
884 : CLASS(mp_comm_type), INTENT(IN) :: int_group
885 :
886 : INTEGER :: i, iparticle, iparticle_kind, &
887 : iparticle_local, j, nparticle_local
888 : REAL(KIND=dp) :: mass
889 : TYPE(atomic_kind_type), POINTER :: atomic_kind
890 :
891 56082 : pv_kin = 0.0_dp
892 4314 : kin = 0._dp
893 13526 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
894 9212 : atomic_kind => atomic_kind_set(iparticle_kind)
895 9212 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
896 9212 : nparticle_local = local_particles%n_el(iparticle_kind)
897 683991 : DO iparticle_local = 1, nparticle_local
898 670465 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
899 2691072 : DO i = 1, 3
900 8045580 : DO j = 1, 3
901 : pv_kin(i, j) = pv_kin(i, j) + &
902 8045580 : mass*vel(i, iparticle)*vel(j, iparticle)
903 : END DO
904 2681860 : kin = kin + mass*vel(i, iparticle)*vel(i, iparticle)
905 : END DO
906 : END DO
907 : END DO
908 :
909 107850 : CALL int_group%sum(pv_kin)
910 4314 : CALL int_group%sum(kin)
911 :
912 : ! updating the constraint virial
913 4314 : IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
914 : molecule_set, &
915 : molecule_kind_set, &
916 : particle_set, virial, &
917 2472 : int_group)
918 :
919 4314 : END SUBROUTINE update_pv_velocity
920 :
921 : ! **************************************************************************************************
922 : !> \brief Routine to compute veps
923 : !> \param box ...
924 : !> \param npt ...
925 : !> \param simpar ...
926 : !> \param pv_kin ...
927 : !> \param kin ...
928 : !> \param virial ...
929 : !> \param infree ...
930 : !> \param virial_components ...
931 : !> \par History
932 : !> none
933 : !> \author CJM
934 : ! **************************************************************************************************
935 7982 : SUBROUTINE update_veps(box, npt, simpar, pv_kin, kin, virial, infree, virial_components)
936 :
937 : TYPE(cell_type), POINTER :: box
938 : TYPE(npt_info_type), DIMENSION(:, :), &
939 : INTENT(INOUT) :: npt
940 : TYPE(simpar_type), INTENT(IN) :: simpar
941 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pv_kin
942 : REAL(KIND=dp), INTENT(IN) :: kin
943 : TYPE(virial_type), INTENT(INOUT) :: virial
944 : REAL(KIND=dp), INTENT(IN) :: infree
945 : INTEGER, INTENT(IN), OPTIONAL :: virial_components
946 :
947 : INTEGER :: ii
948 : REAL(KIND=dp) :: fdotr, trace, v, v0, v0i, vi
949 : REAL(KIND=dp), DIMENSION(3, 3) :: unit
950 :
951 7982 : CPASSERT(ASSOCIATED(box))
952 :
953 : ! Initialize unit
954 :
955 7982 : unit = 0.0_dp
956 7982 : unit(1, 1) = 1.0_dp
957 7982 : unit(2, 2) = 1.0_dp
958 7982 : unit(3, 3) = 1.0_dp
959 :
960 : IF (simpar%ensemble == npt_i_ensemble .OR. &
961 7982 : simpar%ensemble == npt_ia_ensemble .OR. &
962 : simpar%ensemble == npe_i_ensemble) THEN
963 : ! get force on barostat
964 : fdotr = 0.0_dp
965 24160 : DO ii = 1, 3
966 : fdotr = fdotr + virial%pv_virial(ii, ii) + &
967 24160 : virial%pv_constraint(ii, ii)
968 : END DO
969 :
970 : npt(:, :)%f = (1.0_dp + (3.0_dp*infree))*kin + fdotr - &
971 18120 : 3.0_dp*simpar%p_ext*box%deth
972 1942 : ELSEIF (simpar%ensemble == npt_f_ensemble .OR. &
973 : simpar%ensemble == npe_f_ensemble) THEN
974 : npt(:, :)%f = virial%pv_virial(:, :) + &
975 : pv_kin(:, :) + virial%pv_constraint(:, :) - &
976 : unit(:, :)*simpar%p_ext*box%deth + &
977 23686 : infree*kin*unit(:, :)
978 : IF (debug_isotropic_limit) THEN
979 : trace = npt(1, 1)%f + npt(2, 2)%f + npt(3, 3)%f
980 : trace = trace/3.0_dp
981 : npt(:, :)%f = trace*unit(:, :)
982 : END IF
983 120 : ELSEIF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
984 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
985 120 : v = box%deth
986 120 : vi = 1._dp/v
987 120 : v0 = simpar%v0
988 120 : v0i = 1._dp/v0
989 : ! orthorhombic box ONLY
990 : ! Chooses only the compressive solution
991 120 : IF (v < v0) THEN
992 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
993 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
994 : simpar%p0*v - simpar%v_shock*simpar%v_shock* &
995 0 : v*v0i*(1._dp - v*v0i) + infree*kin
996 : ELSE
997 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
998 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
999 120 : simpar%p0*v + infree*kin
1000 : END IF
1001 : IF (debug_uniaxial_limit) THEN
1002 : ! orthorhombic box ONLY
1003 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
1004 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
1005 : simpar%p0*box%deth + infree*kin
1006 : END IF
1007 : END IF
1008 :
1009 : ! update barostat velocities
1010 : npt(:, :)%v = npt(:, :)%v + &
1011 42166 : 0.5_dp*simpar%dt*npt(:, :)%f/npt(:, :)%mass
1012 :
1013 : ! Screen the dynamics of the barostat according user request
1014 7982 : IF (PRESENT(virial_components)) THEN
1015 1812 : SELECT CASE (virial_components)
1016 : CASE (do_clv_xyz)
1017 : ! Do nothing..
1018 : CASE (do_clv_x)
1019 0 : npt(1, 2)%v = 0.0_dp
1020 0 : npt(1, 3)%v = 0.0_dp
1021 0 : npt(1, 2)%f = 0.0_dp
1022 0 : npt(1, 3)%f = 0.0_dp
1023 0 : npt(2, 1)%v = 0.0_dp
1024 0 : npt(2, 2)%v = 0.0_dp
1025 0 : npt(2, 3)%v = 0.0_dp
1026 0 : npt(2, 1)%f = 0.0_dp
1027 0 : npt(2, 2)%f = 0.0_dp
1028 0 : npt(2, 3)%f = 0.0_dp
1029 0 : npt(3, 1)%v = 0.0_dp
1030 0 : npt(3, 2)%v = 0.0_dp
1031 0 : npt(3, 3)%v = 0.0_dp
1032 0 : npt(3, 1)%f = 0.0_dp
1033 0 : npt(3, 2)%f = 0.0_dp
1034 0 : npt(3, 3)%f = 0.0_dp
1035 : CASE (do_clv_y)
1036 0 : npt(1, 1)%v = 0.0_dp
1037 0 : npt(1, 2)%v = 0.0_dp
1038 0 : npt(1, 3)%v = 0.0_dp
1039 0 : npt(1, 1)%f = 0.0_dp
1040 0 : npt(1, 2)%f = 0.0_dp
1041 0 : npt(1, 3)%f = 0.0_dp
1042 0 : npt(2, 1)%v = 0.0_dp
1043 0 : npt(2, 3)%v = 0.0_dp
1044 0 : npt(2, 1)%f = 0.0_dp
1045 0 : npt(2, 3)%f = 0.0_dp
1046 0 : npt(3, 1)%v = 0.0_dp
1047 0 : npt(3, 2)%v = 0.0_dp
1048 0 : npt(3, 3)%v = 0.0_dp
1049 0 : npt(3, 1)%f = 0.0_dp
1050 0 : npt(3, 2)%f = 0.0_dp
1051 0 : npt(3, 3)%f = 0.0_dp
1052 : CASE (do_clv_z)
1053 80 : npt(1, 1)%v = 0.0_dp
1054 80 : npt(1, 2)%v = 0.0_dp
1055 80 : npt(1, 3)%v = 0.0_dp
1056 80 : npt(1, 1)%f = 0.0_dp
1057 80 : npt(1, 2)%f = 0.0_dp
1058 80 : npt(1, 3)%f = 0.0_dp
1059 80 : npt(2, 1)%v = 0.0_dp
1060 80 : npt(2, 2)%v = 0.0_dp
1061 80 : npt(2, 3)%v = 0.0_dp
1062 80 : npt(2, 1)%f = 0.0_dp
1063 80 : npt(2, 2)%f = 0.0_dp
1064 80 : npt(2, 3)%f = 0.0_dp
1065 80 : npt(3, 1)%v = 0.0_dp
1066 80 : npt(3, 2)%v = 0.0_dp
1067 80 : npt(3, 1)%f = 0.0_dp
1068 80 : npt(3, 2)%f = 0.0_dp
1069 : CASE (do_clv_xy)
1070 0 : npt(1, 3)%v = 0.0_dp
1071 0 : npt(1, 3)%f = 0.0_dp
1072 0 : npt(2, 3)%v = 0.0_dp
1073 0 : npt(2, 3)%f = 0.0_dp
1074 0 : npt(3, 1)%v = 0.0_dp
1075 0 : npt(3, 2)%v = 0.0_dp
1076 0 : npt(3, 3)%v = 0.0_dp
1077 0 : npt(3, 1)%f = 0.0_dp
1078 0 : npt(3, 2)%f = 0.0_dp
1079 0 : npt(3, 3)%f = 0.0_dp
1080 : CASE (do_clv_xz)
1081 0 : npt(1, 2)%v = 0.0_dp
1082 0 : npt(1, 2)%f = 0.0_dp
1083 0 : npt(2, 1)%v = 0.0_dp
1084 0 : npt(2, 2)%v = 0.0_dp
1085 0 : npt(2, 3)%v = 0.0_dp
1086 0 : npt(2, 1)%f = 0.0_dp
1087 0 : npt(2, 2)%f = 0.0_dp
1088 0 : npt(2, 3)%f = 0.0_dp
1089 0 : npt(3, 2)%v = 0.0_dp
1090 0 : npt(3, 2)%f = 0.0_dp
1091 : CASE (do_clv_yz)
1092 0 : npt(1, 1)%v = 0.0_dp
1093 0 : npt(1, 2)%v = 0.0_dp
1094 0 : npt(1, 3)%v = 0.0_dp
1095 0 : npt(1, 1)%f = 0.0_dp
1096 0 : npt(1, 2)%f = 0.0_dp
1097 0 : npt(1, 3)%f = 0.0_dp
1098 0 : npt(2, 1)%v = 0.0_dp
1099 0 : npt(2, 1)%f = 0.0_dp
1100 0 : npt(3, 1)%v = 0.0_dp
1101 1812 : npt(3, 1)%f = 0.0_dp
1102 : END SELECT
1103 : END IF
1104 :
1105 7982 : END SUBROUTINE update_veps
1106 :
1107 : ! **************************************************************************************************
1108 : !> \brief First half of the velocity-verlet algorithm : update velocity by half
1109 : !> step and positions by full step
1110 : !> \param tmp ...
1111 : !> \param atomic_kind_set ...
1112 : !> \param local_particles ...
1113 : !> \param particle_set ...
1114 : !> \param core_particle_set ...
1115 : !> \param shell_particle_set ...
1116 : !> \param nparticle_kind ...
1117 : !> \param shell_adiabatic ...
1118 : !> \param dt ...
1119 : !> \param u ...
1120 : !> \param lfixd_list ...
1121 : !> \par History
1122 : !> none
1123 : !> \author MI (February 2008)
1124 : ! **************************************************************************************************
1125 44235 : SUBROUTINE vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1126 : core_particle_set, shell_particle_set, nparticle_kind, &
1127 44235 : shell_adiabatic, dt, u, lfixd_list)
1128 :
1129 : TYPE(tmp_variables_type), POINTER :: tmp
1130 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1131 : TYPE(distribution_1d_type), POINTER :: local_particles
1132 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1133 : shell_particle_set
1134 : INTEGER, INTENT(IN) :: nparticle_kind
1135 : LOGICAL, INTENT(IN) :: shell_adiabatic
1136 : REAL(KIND=dp) :: dt
1137 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
1138 : TYPE(local_fixd_constraint_type), DIMENSION(:), &
1139 : OPTIONAL :: lfixd_list
1140 :
1141 : INTEGER :: iparticle, iparticle_kind, &
1142 : iparticle_local, nparticle_local, &
1143 : shell_index
1144 : LOGICAL :: is_shell
1145 : REAL(KIND=dp) :: dm, dmc, dms, dsc(3), dvsc(3), &
1146 : fac_massc, fac_masss, mass
1147 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1148 : TYPE(shell_kind_type), POINTER :: shell
1149 :
1150 44235 : NULLIFY (atomic_kind, shell)
1151 44235 : tmp%max_vel = 0.0_dp
1152 44235 : tmp%max_vel_sc = 0.0_dp
1153 44235 : tmp%max_dvel = 0.0_dp
1154 44235 : tmp%max_dvel_sc = 0.0_dp
1155 44235 : tmp%max_dr = 0.0_dp
1156 44235 : tmp%max_dsc = 0.0_dp
1157 44235 : dsc = 0.0_dp
1158 44235 : dvsc = 0.0_dp
1159 :
1160 : ! *** Velocity Verlet (first part) ***
1161 154510 : DO iparticle_kind = 1, nparticle_kind
1162 110275 : atomic_kind => atomic_kind_set(iparticle_kind)
1163 110275 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell)
1164 154510 : IF (mass /= 0.0_dp) THEN
1165 109869 : dm = 0.5_dp*dt/mass
1166 109869 : IF (is_shell .AND. shell_adiabatic) THEN
1167 5412 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1168 5412 : dms = 0.5_dp*dt/shell%mass_shell
1169 5412 : dmc = 0.5_dp*dt/shell%mass_core
1170 5412 : fac_masss = shell%mass_shell/mass
1171 5412 : fac_massc = shell%mass_core/mass
1172 5412 : nparticle_local = local_particles%n_el(iparticle_kind)
1173 5412 : IF (PRESENT(u)) THEN
1174 54620 : DO iparticle_local = 1, nparticle_local
1175 52704 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1176 52704 : shell_index = particle_set(iparticle)%shell_index
1177 : ! Transform positions and velocities and forces of the shell
1178 : CALL transform_first(shell_particle_set, tmp%shell_pos, tmp%shell_vel, &
1179 52704 : shell_index, u, dms, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1180 :
1181 : ! Transform positions and velocities and forces of the core
1182 : CALL transform_first(core_particle_set, tmp%core_pos, tmp%core_vel, &
1183 52704 : shell_index, u, dmc, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1184 :
1185 : ! Derive velocities and positions of the COM
1186 : tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
1187 210816 : fac_massc*tmp%core_vel(:, shell_index)
1188 : tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1189 210816 : fac_massc*tmp%core_pos(:, shell_index)
1190 :
1191 : tmp%max_vel = MAX(tmp%max_vel, ABS(tmp%vel(1, iparticle)), &
1192 52704 : ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1193 : tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
1194 : ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
1195 : ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
1196 52704 : ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
1197 : tmp%max_dr = MAX(tmp%max_dr, &
1198 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1199 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1200 52704 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1201 : tmp%max_dvel = MAX(tmp%max_dvel, &
1202 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1203 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1204 52704 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1205 :
1206 : dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1207 210816 : shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1208 52704 : tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
1209 : dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
1210 210816 : shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1211 54620 : tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
1212 : END DO ! iparticle_local
1213 : ELSE
1214 88530 : DO iparticle_local = 1, nparticle_local
1215 85034 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1216 85034 : shell_index = particle_set(iparticle)%shell_index
1217 : tmp%shell_vel(:, shell_index) = &
1218 : shell_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
1219 680272 : tmp%scale_v(:)*tmp%poly_v(:)*dms*shell_particle_set(shell_index)%f(:)
1220 : tmp%shell_pos(:, shell_index) = &
1221 : shell_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
1222 680272 : tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%shell_vel(:, shell_index)
1223 : tmp%core_vel(:, shell_index) = &
1224 : core_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
1225 680272 : tmp%scale_v(:)*tmp%poly_v(:)*dmc*core_particle_set(shell_index)%f(:)
1226 : tmp%core_pos(:, shell_index) = &
1227 : core_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
1228 680272 : tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%core_vel(:, shell_index)
1229 :
1230 : tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
1231 340136 : fac_massc*tmp%core_vel(:, shell_index)
1232 : tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1233 340136 : fac_massc*tmp%core_pos(:, shell_index)
1234 : tmp%max_vel = MAX(tmp%max_vel, &
1235 85034 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1236 : tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
1237 : ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
1238 : ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
1239 85034 : ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
1240 : tmp%max_dr = MAX(tmp%max_dr, &
1241 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1242 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1243 85034 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1244 : tmp%max_dvel = MAX(tmp%max_dvel, &
1245 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1246 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1247 85034 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1248 : dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1249 340136 : shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1250 85034 : tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
1251 : dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
1252 340136 : shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1253 88530 : tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
1254 : END DO ! iparticle_local
1255 : END IF
1256 : ELSE
1257 104457 : nparticle_local = local_particles%n_el(iparticle_kind)
1258 104457 : IF (PRESENT(u)) THEN
1259 44038 : DO iparticle_local = 1, nparticle_local
1260 43412 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1261 : ! Transform positions and velocities and forces
1262 : CALL transform_first(particle_set, tmp%pos, tmp%vel, &
1263 43412 : iparticle, u, dm, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1264 :
1265 : tmp%max_vel = MAX(tmp%max_vel, &
1266 43412 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1267 : tmp%max_dr = MAX(tmp%max_dr, ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1268 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1269 43412 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1270 : tmp%max_dvel = MAX(tmp%max_dvel, &
1271 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1272 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1273 44038 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1274 : END DO ! iparticle_local
1275 : ELSE
1276 3174720 : DO iparticle_local = 1, nparticle_local
1277 3070889 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1278 : tmp%vel(1, iparticle) = &
1279 : particle_set(iparticle)%v(1)*tmp%scale_v(1)*tmp%scale_v(1) + &
1280 3070889 : tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
1281 : tmp%vel(2, iparticle) = &
1282 : particle_set(iparticle)%v(2)*tmp%scale_v(2)*tmp%scale_v(2) + &
1283 3070889 : tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
1284 : tmp%vel(3, iparticle) = &
1285 : particle_set(iparticle)%v(3)*tmp%scale_v(3)*tmp%scale_v(3) + &
1286 3070889 : tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
1287 : tmp%pos(1, iparticle) = &
1288 : particle_set(iparticle)%r(1)*tmp%scale_r(1)*tmp%scale_r(1) + &
1289 3070889 : tmp%scale_r(1)*tmp%poly_r(1)*dt*tmp%vel(1, iparticle)
1290 : tmp%pos(2, iparticle) = &
1291 : particle_set(iparticle)%r(2)*tmp%scale_r(2)*tmp%scale_r(2) + &
1292 3070889 : tmp%scale_r(2)*tmp%poly_r(2)*dt*tmp%vel(2, iparticle)
1293 : tmp%pos(3, iparticle) = &
1294 : particle_set(iparticle)%r(3)*tmp%scale_r(3)*tmp%scale_r(3) + &
1295 3070889 : tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
1296 :
1297 : ! overwrite positions of frozen particles
1298 3070889 : IF (PRESENT(lfixd_list)) THEN
1299 251570 : IF (ANY(lfixd_list(:)%ifixd_index == iparticle)) THEN
1300 200 : tmp%pos(1, iparticle) = particle_set(iparticle)%r(1)
1301 200 : tmp%pos(2, iparticle) = particle_set(iparticle)%r(2)
1302 200 : tmp%pos(3, iparticle) = particle_set(iparticle)%r(3)
1303 : END IF
1304 : END IF
1305 :
1306 : tmp%max_vel = MAX(tmp%max_vel, &
1307 3070889 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1308 : tmp%max_dr = MAX(tmp%max_dr, &
1309 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1310 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1311 3070889 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1312 : tmp%max_dvel = MAX(tmp%max_dvel, &
1313 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1314 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1315 3174720 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1316 : END DO
1317 : END IF
1318 : END IF
1319 : END IF
1320 : END DO
1321 44235 : END SUBROUTINE vv_first
1322 :
1323 : ! **************************************************************************************************
1324 : !> \brief Second half of the velocity-verlet algorithm : update velocity by half
1325 : !> step using the new forces
1326 : !> \param tmp ...
1327 : !> \param atomic_kind_set ...
1328 : !> \param local_particles ...
1329 : !> \param particle_set ...
1330 : !> \param core_particle_set ...
1331 : !> \param shell_particle_set ...
1332 : !> \param nparticle_kind ...
1333 : !> \param shell_adiabatic ...
1334 : !> \param dt ...
1335 : !> \param u ...
1336 : !> \par History
1337 : !> none
1338 : !> \author MI (February 2008)
1339 : ! **************************************************************************************************
1340 43761 : SUBROUTINE vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
1341 : core_particle_set, shell_particle_set, nparticle_kind, &
1342 : shell_adiabatic, dt, u)
1343 :
1344 : TYPE(tmp_variables_type), POINTER :: tmp
1345 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1346 : TYPE(distribution_1d_type), POINTER :: local_particles
1347 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1348 : shell_particle_set
1349 : INTEGER, INTENT(IN) :: nparticle_kind
1350 : LOGICAL, INTENT(IN) :: shell_adiabatic
1351 : REAL(KIND=dp) :: dt
1352 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
1353 :
1354 : INTEGER :: iparticle, iparticle_kind, &
1355 : iparticle_local, nparticle_local, &
1356 : shell_index
1357 : LOGICAL :: is_shell
1358 : REAL(KIND=dp) :: dm, dmc, dms, fac_massc, fac_masss, mass
1359 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1360 : TYPE(shell_kind_type), POINTER :: shell
1361 :
1362 43761 : NULLIFY (atomic_kind, shell)
1363 :
1364 : ! *** Velocity Verlet (second part) ***
1365 152828 : DO iparticle_kind = 1, nparticle_kind
1366 109067 : atomic_kind => atomic_kind_set(iparticle_kind)
1367 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
1368 109067 : shell_active=is_shell)
1369 152828 : IF (mass /= 0.0_dp) THEN
1370 108661 : dm = 0.5_dp*dt/mass
1371 108661 : IF (is_shell .AND. shell_adiabatic) THEN
1372 4520 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1373 4520 : dms = 0.5_dp*dt/shell%mass_shell
1374 4520 : dmc = 0.5_dp*dt/shell%mass_core
1375 4520 : fac_masss = shell%mass_shell/mass
1376 4520 : fac_massc = shell%mass_core/mass
1377 4520 : nparticle_local = local_particles%n_el(iparticle_kind)
1378 4520 : IF (PRESENT(u)) THEN
1379 35860 : DO iparticle_local = 1, nparticle_local
1380 34560 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1381 34560 : shell_index = particle_set(iparticle)%shell_index
1382 : ! Transform velocities and forces shell
1383 : CALL transform_second(shell_particle_set, tmp%shell_vel, shell_index, &
1384 34560 : u, dms, tmp%poly_v, tmp%scale_v)
1385 :
1386 : ! Transform velocities and forces core
1387 : CALL transform_second(core_particle_set, tmp%core_vel, shell_index, &
1388 34560 : u, dmc, tmp%poly_v, tmp%scale_v)
1389 :
1390 : ! Derive velocties of the COM
1391 : tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
1392 34560 : fac_massc*tmp%core_vel(1, shell_index)
1393 : tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1394 34560 : fac_massc*tmp%core_vel(2, shell_index)
1395 : tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1396 35860 : fac_massc*tmp%core_vel(3, shell_index)
1397 : END DO ! iparticle_local
1398 : ELSE
1399 81630 : DO iparticle_local = 1, nparticle_local
1400 78410 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1401 78410 : shell_index = particle_set(iparticle)%shell_index
1402 : tmp%shell_vel(1, shell_index) = &
1403 : tmp%shell_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
1404 78410 : tmp%scale_v(1)*tmp%poly_v(1)*dms*shell_particle_set(shell_index)%f(1)
1405 : tmp%shell_vel(2, shell_index) = &
1406 : tmp%shell_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
1407 78410 : tmp%scale_v(2)*tmp%poly_v(2)*dms*shell_particle_set(shell_index)%f(2)
1408 : tmp%shell_vel(3, shell_index) = &
1409 : tmp%shell_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
1410 78410 : tmp%scale_v(3)*tmp%poly_v(3)*dms*shell_particle_set(shell_index)%f(3)
1411 : tmp%core_vel(1, shell_index) = &
1412 : tmp%core_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
1413 78410 : tmp%scale_v(1)*tmp%poly_v(1)*dmc*core_particle_set(shell_index)%f(1)
1414 : tmp%core_vel(2, shell_index) = &
1415 : tmp%core_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
1416 78410 : tmp%scale_v(2)*tmp%poly_v(2)*dmc*core_particle_set(shell_index)%f(2)
1417 : tmp%core_vel(3, shell_index) = &
1418 : tmp%core_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
1419 78410 : tmp%scale_v(3)*tmp%poly_v(3)*dmc*core_particle_set(shell_index)%f(3)
1420 :
1421 : tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
1422 78410 : fac_massc*tmp%core_vel(1, shell_index)
1423 : tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1424 78410 : fac_massc*tmp%core_vel(2, shell_index)
1425 : tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1426 81630 : fac_massc*tmp%core_vel(3, shell_index)
1427 : END DO ! iparticle_local
1428 : END IF
1429 : ELSE
1430 104141 : nparticle_local = local_particles%n_el(iparticle_kind)
1431 104141 : IF (PRESENT(u)) THEN
1432 44038 : DO iparticle_local = 1, nparticle_local
1433 43412 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1434 : CALL transform_second(particle_set, tmp%vel, iparticle, &
1435 44038 : u, dm, tmp%poly_v, tmp%scale_v)
1436 : END DO ! iparticle_local
1437 : ELSE
1438 2944176 : DO iparticle_local = 1, nparticle_local
1439 2840661 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1440 : tmp%vel(1, iparticle) = &
1441 : tmp%vel(1, iparticle)*tmp%scale_v(1)*tmp%scale_v(1) + &
1442 2840661 : tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
1443 : tmp%vel(2, iparticle) = &
1444 : tmp%vel(2, iparticle)*tmp%scale_v(2)*tmp%scale_v(2) + &
1445 2840661 : tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
1446 : tmp%vel(3, iparticle) = &
1447 : tmp%vel(3, iparticle)*tmp%scale_v(3)*tmp%scale_v(3) + &
1448 2944176 : tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
1449 : END DO
1450 : END IF
1451 : END IF
1452 : END IF
1453 : END DO
1454 :
1455 43761 : END SUBROUTINE vv_second
1456 :
1457 : ! **************************************************************************************************
1458 : !> \brief Transform position and velocities for the first half of the
1459 : !> Velocity-Verlet integration in the npt_f ensemble
1460 : !> \param particle_set ...
1461 : !> \param pos ...
1462 : !> \param vel ...
1463 : !> \param index ...
1464 : !> \param u ...
1465 : !> \param dm ...
1466 : !> \param dt ...
1467 : !> \param poly_v ...
1468 : !> \param poly_r ...
1469 : !> \param scale_v ...
1470 : !> \param scale_r ...
1471 : !> \par History
1472 : !> none
1473 : !> \author MI (February 2008)
1474 : ! **************************************************************************************************
1475 148820 : SUBROUTINE transform_first(particle_set, pos, vel, index, u, dm, dt, poly_v, &
1476 : poly_r, scale_v, scale_r)
1477 :
1478 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1479 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pos, vel
1480 : INTEGER, INTENT(IN) :: index
1481 : REAL(KIND=dp), INTENT(IN) :: u(3, 3), dm, dt, poly_v(3), poly_r(3), &
1482 : scale_v(3), scale_r(3)
1483 :
1484 : REAL(KIND=dp), DIMENSION(3) :: uf, ur, uv
1485 :
1486 148820 : ur = MATMUL(TRANSPOSE(u), particle_set(index)%r(:))
1487 148820 : uv = MATMUL(TRANSPOSE(u), particle_set(index)%v(:))
1488 148820 : uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
1489 :
1490 148820 : uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1491 148820 : uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1492 148820 : uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1493 :
1494 148820 : ur(1) = ur(1)*scale_r(1)*scale_r(1) + uv(1)*scale_r(1)*poly_r(1)*dt
1495 148820 : ur(2) = ur(2)*scale_r(2)*scale_r(2) + uv(2)*scale_r(2)*poly_r(2)*dt
1496 148820 : ur(3) = ur(3)*scale_r(3)*scale_r(3) + uv(3)*scale_r(3)*poly_r(3)*dt
1497 :
1498 2529940 : pos(:, index) = MATMUL(u, ur)
1499 2529940 : vel(:, index) = MATMUL(u, uv)
1500 :
1501 148820 : END SUBROUTINE transform_first
1502 :
1503 : ! **************************************************************************************************
1504 : !> \brief Transform position and velocities for the second half of the
1505 : !> Velocity-Verlet integration in the npt_f ensemble
1506 : !> \param particle_set ...
1507 : !> \param vel ...
1508 : !> \param index ...
1509 : !> \param u ...
1510 : !> \param dm ...
1511 : !> \param poly_v ...
1512 : !> \param scale_v ...
1513 : !> \par History
1514 : !> none
1515 : !> \author MI (February 2008)
1516 : ! **************************************************************************************************
1517 112532 : SUBROUTINE transform_second(particle_set, vel, index, u, dm, poly_v, scale_v)
1518 :
1519 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1520 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vel
1521 : INTEGER, INTENT(IN) :: index
1522 : REAL(KIND=dp), INTENT(IN) :: u(3, 3), dm, poly_v(3), scale_v(3)
1523 :
1524 : REAL(KIND=dp), DIMENSION(3) :: uf, uv
1525 :
1526 112532 : uv = MATMUL(TRANSPOSE(u), vel(:, index))
1527 112532 : uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
1528 :
1529 112532 : uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1530 112532 : uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1531 112532 : uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1532 :
1533 1913044 : vel(:, index) = MATMUL(u, uv)
1534 :
1535 112532 : END SUBROUTINE transform_second
1536 :
1537 : ! **************************************************************************************************
1538 : !> \brief Compute the timestep rescaling factor
1539 : !> \param md_env ...
1540 : !> \param tmp ...
1541 : !> \param dt ...
1542 : !> \param simpar ...
1543 : !> \param para_env ...
1544 : !> \param atomic_kind_set ...
1545 : !> \param local_particles ...
1546 : !> \param particle_set ...
1547 : !> \param core_particle_set ...
1548 : !> \param shell_particle_set ...
1549 : !> \param nparticle_kind ...
1550 : !> \param shell_adiabatic ...
1551 : !> \param npt ...
1552 : !> \par History
1553 : !> none
1554 : !> \author MI (October 2008)
1555 : ! **************************************************************************************************
1556 480 : SUBROUTINE variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, local_particles, &
1557 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1558 :
1559 : TYPE(md_environment_type), POINTER :: md_env
1560 : TYPE(tmp_variables_type), POINTER :: tmp
1561 : REAL(KIND=dp), INTENT(INOUT) :: dt
1562 : TYPE(simpar_type), POINTER :: simpar
1563 : TYPE(mp_para_env_type), POINTER :: para_env
1564 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1565 : TYPE(distribution_1d_type), POINTER :: local_particles
1566 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1567 : shell_particle_set
1568 : INTEGER, INTENT(IN) :: nparticle_kind
1569 : LOGICAL, INTENT(IN) :: shell_adiabatic
1570 : TYPE(npt_info_type), OPTIONAL, POINTER :: npt(:, :)
1571 :
1572 : INTEGER, SAVE :: itime = 0
1573 : REAL(KIND=dp) :: dt_fact, dt_fact_f, dt_fact_old, &
1574 : dt_fact_sc, dt_fact_v, dt_old
1575 : REAL(KIND=dp), POINTER :: time
1576 : TYPE(thermostats_type), POINTER :: thermostats
1577 :
1578 480 : dt_fact = 1.0_dp
1579 480 : dt_fact_sc = 1.0_dp
1580 480 : dt_fact_f = 1.0_dp
1581 480 : dt_fact_v = 1.0_dp
1582 480 : dt_old = dt
1583 480 : dt_fact_old = simpar%dt_fact
1584 480 : simpar%dt_fact = 1.0_dp
1585 480 : NULLIFY (thermostats)
1586 :
1587 480 : itime = itime + 1
1588 480 : CALL para_env%max(tmp%max_dr)
1589 480 : IF (tmp%max_dr > simpar%dr_tol) THEN
1590 278 : CALL para_env%max(tmp%max_dvel)
1591 278 : dt_fact = SQRT(simpar%dr_tol*dt/tmp%max_dvel)/dt
1592 : END IF
1593 :
1594 480 : IF (shell_adiabatic) THEN
1595 440 : CALL para_env%max(tmp%max_dsc)
1596 440 : IF (tmp%max_dsc > simpar%dsc_tol) THEN
1597 118 : CALL para_env%max(tmp%max_dvel_sc)
1598 118 : dt_fact_sc = SQRT(simpar%dsc_tol*dt/tmp%max_dvel_sc)/dt
1599 : END IF
1600 : END IF
1601 :
1602 480 : dt_fact_f = MIN(dt_fact, dt_fact_sc, 1.0_dp)
1603 480 : IF (dt_fact_f < 1.0_dp) THEN
1604 100 : IF (dt_fact_f < 0.1_dp) dt_fact_f = 0.1_dp
1605 : ! repeat the first vv half-step with the rescaled time step
1606 100 : dt = dt*dt_fact_f
1607 100 : simpar%dt_fact = dt_fact_f
1608 : CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1609 100 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1610 : END IF
1611 :
1612 480 : dt_fact = 1.0_dp
1613 480 : dt_fact_sc = 1.0_dp
1614 480 : CALL para_env%max(tmp%max_dr)
1615 480 : IF (tmp%max_dr > simpar%dr_tol) THEN
1616 278 : CALL para_env%max(tmp%max_vel)
1617 278 : dt_fact = simpar%dr_tol/tmp%max_vel/dt
1618 : END IF
1619 :
1620 480 : IF (shell_adiabatic) THEN
1621 440 : CALL para_env%max(tmp%max_dsc)
1622 440 : IF (tmp%max_dsc > simpar%dsc_tol) THEN
1623 116 : CALL para_env%max(tmp%max_vel_sc)
1624 116 : dt_fact_sc = simpar%dsc_tol/tmp%max_vel_sc/dt
1625 : END IF
1626 : END IF
1627 480 : dt_fact_v = MIN(dt_fact, dt_fact_sc, 1.0_dp)
1628 :
1629 480 : IF (dt_fact_v < 1.0_dp) THEN
1630 376 : IF (dt_fact_v < 0.1_dp) dt_fact_v = 0.1_dp
1631 : ! repeat the first vv half-step with the rescaled time step
1632 376 : dt = dt*dt_fact_v
1633 376 : simpar%dt_fact = dt_fact_f*dt_fact_v
1634 : CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1635 554 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1636 : END IF
1637 :
1638 480 : simpar%dt_fact = dt_fact_f*dt_fact_v
1639 480 : IF (simpar%dt_fact < 1.0_dp) THEN
1640 378 : CALL get_md_env(md_env, t=time, thermostats=thermostats)
1641 378 : time = time - dt_old + dt_old*simpar%dt_fact
1642 378 : IF (ASSOCIATED(thermostats)) THEN
1643 240 : CALL set_thermostats(thermostats, dt_fact=simpar%dt_fact)
1644 : END IF
1645 : END IF
1646 :
1647 480 : END SUBROUTINE variable_timestep
1648 :
1649 : ! **************************************************************************************************
1650 : !> \brief Repeat the first step of velocity verlet with the rescaled time step
1651 : !> \param tmp ...
1652 : !> \param dt ...
1653 : !> \param simpar ...
1654 : !> \param atomic_kind_set ...
1655 : !> \param local_particles ...
1656 : !> \param particle_set ...
1657 : !> \param core_particle_set ...
1658 : !> \param shell_particle_set ...
1659 : !> \param nparticle_kind ...
1660 : !> \param shell_adiabatic ...
1661 : !> \param npt ...
1662 : !> \par History
1663 : !> none
1664 : !> \author MI (October 2008)
1665 : !> I soliti ignoti
1666 : ! **************************************************************************************************
1667 476 : SUBROUTINE rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1668 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1669 :
1670 : TYPE(tmp_variables_type), POINTER :: tmp
1671 : REAL(KIND=dp), INTENT(IN) :: dt
1672 : TYPE(simpar_type), POINTER :: simpar
1673 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1674 : TYPE(distribution_1d_type), POINTER :: local_particles
1675 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1676 : shell_particle_set
1677 : INTEGER, INTENT(IN) :: nparticle_kind
1678 : LOGICAL, INTENT(IN) :: shell_adiabatic
1679 : TYPE(npt_info_type), OPTIONAL, POINTER :: npt(:, :)
1680 :
1681 : REAL(KIND=dp), PARAMETER :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
1682 : e6 = e4/42.0_dp, e8 = e6/72.0_dp
1683 :
1684 : REAL(KIND=dp) :: arg_r(3), arg_v(3), e_val(3), infree, &
1685 : trvg, u(3, 3)
1686 :
1687 476 : infree = 1.0_dp/REAL(simpar%nfree, dp)
1688 1904 : arg_r = tmp%arg_r
1689 1904 : arg_v = tmp%arg_v
1690 6188 : u = tmp%u
1691 1904 : e_val = tmp%e_val
1692 :
1693 654 : SELECT CASE (simpar%ensemble)
1694 : CASE (nve_ensemble, nvt_ensemble)
1695 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1696 178 : core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
1697 : CASE (isokin_ensemble)
1698 0 : tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds)
1699 0 : tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt
1700 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1701 : core_particle_set, shell_particle_set, nparticle_kind, &
1702 0 : shell_adiabatic, dt)
1703 :
1704 : CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
1705 0 : arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1706 0 : tmp%poly_r(1:3) = 1.0_dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
1707 0 : arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
1708 0 : tmp%poly_v(1:3) = 1.0_dp + e2*arg_v(1) + e4*arg_v(1)*arg_v(1) + e6*arg_v(1)**3 + e8*arg_v(1)**4
1709 0 : tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
1710 : tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
1711 0 : (1.0_dp + 3.0_dp*infree))
1712 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1713 : core_particle_set, shell_particle_set, nparticle_kind, &
1714 0 : shell_adiabatic, dt)
1715 :
1716 : CASE (npt_f_ensemble, npe_f_ensemble)
1717 298 : trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
1718 1192 : arg_r(:) = arg_r(:)*simpar%dt_fact*simpar%dt_fact
1719 1192 : tmp%poly_r = 1._dp + e2*arg_r + e4*arg_r*arg_r + e6*arg_r**3 + e8*arg_r**4
1720 1192 : tmp%scale_r(:) = EXP(0.5_dp*dt*e_val(:))
1721 1192 : arg_v(:) = arg_v(:)*simpar%dt_fact*simpar%dt_fact
1722 1192 : tmp%poly_v = 1.0_dp + e2*arg_v + e4*arg_v*arg_v + e6*arg_v**3 + e8*arg_v**4
1723 : tmp%scale_v(:) = EXP(-0.25_dp*dt*( &
1724 1192 : e_val(:) + trvg*infree))
1725 :
1726 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1727 : core_particle_set, shell_particle_set, nparticle_kind, &
1728 298 : shell_adiabatic, dt, u)
1729 :
1730 : CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
1731 0 : arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1732 0 : tmp%poly_r(1) = 1._dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
1733 0 : arg_v(2) = arg_v(2)*simpar%dt_fact*simpar%dt_fact
1734 0 : arg_v(1) = arg_v(1)*simpar%dt_fact*simpar%dt_fact
1735 0 : tmp%poly_v(1) = 1._dp + e2*arg_v(1) + e4*arg_v(1)*arg_v(1) + e6*arg_v(1)**3 + e8*arg_v(1)**4
1736 0 : tmp%poly_v(2) = 1._dp + e2*arg_v(2) + e4*arg_v(2)*arg_v(2) + e6*arg_v(2)**3 + e8*arg_v(2)**4
1737 0 : tmp%poly_v(3) = 1._dp + e2*arg_v(2) + e4*arg_v(2)*arg_v(2) + e6*arg_v(2)**3 + e8*arg_v(2)**4
1738 0 : tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v)
1739 : tmp%scale_v(1) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
1740 0 : (1._dp + infree))
1741 0 : tmp%scale_v(2) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
1742 0 : tmp%scale_v(3) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
1743 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1744 : core_particle_set, shell_particle_set, nparticle_kind, &
1745 476 : shell_adiabatic, dt)
1746 :
1747 : END SELECT
1748 :
1749 476 : END SUBROUTINE rescaled_vv_first
1750 :
1751 0 : END MODULE integrator_utils
|
