Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Provides integrator 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 2108 : 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 2108 : natoms = SIZE(particle_set)
123 2108 : idim = SIZE(npt, 1)
124 2108 : jdim = SIZE(npt, 2)
125 2108 : CPASSERT(.NOT. ASSOCIATED(old))
126 2108 : ALLOCATE (old)
127 :
128 6324 : ALLOCATE (old%v(natoms, 3))
129 1383752 : old%v = 0.0_dp
130 6324 : ALLOCATE (old%r(natoms, 3))
131 1383752 : old%r = 0.0_dp
132 8432 : ALLOCATE (old%eps(idim, jdim))
133 15284 : old%eps = 0.0_dp
134 8432 : ALLOCATE (old%veps(idim, jdim))
135 15284 : old%veps = 0.0_dp
136 2108 : ALLOCATE (old%h(3, 3))
137 27404 : old%h = 0.0_dp
138 :
139 2108 : END SUBROUTINE allocate_old
140 :
141 : ! **************************************************************************************************
142 : !> \brief ...
143 : !> \param old ...
144 : !> \par History
145 : !> none
146 : !> \author CJM
147 : ! **************************************************************************************************
148 2108 : SUBROUTINE deallocate_old(old)
149 : TYPE(old_variables_type), POINTER :: old
150 :
151 2108 : IF (ASSOCIATED(old)) THEN
152 2108 : IF (ASSOCIATED(old%v)) THEN
153 2108 : DEALLOCATE (old%v)
154 : END IF
155 2108 : IF (ASSOCIATED(old%r)) THEN
156 2108 : DEALLOCATE (old%r)
157 : END IF
158 2108 : IF (ASSOCIATED(old%eps)) THEN
159 2108 : DEALLOCATE (old%eps)
160 : END IF
161 2108 : IF (ASSOCIATED(old%veps)) THEN
162 2108 : DEALLOCATE (old%veps)
163 : END IF
164 2108 : IF (ASSOCIATED(old%h)) THEN
165 2108 : DEALLOCATE (old%h)
166 : END IF
167 2108 : DEALLOCATE (old)
168 : END IF
169 :
170 2108 : 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 39185 : 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 39185 : CPASSERT(.NOT. ASSOCIATED(tmp))
191 39185 : ALLOCATE (tmp)
192 :
193 : ! Nullify Components
194 39185 : NULLIFY (tmp%pos)
195 39185 : NULLIFY (tmp%vel)
196 39185 : NULLIFY (tmp%shell_pos)
197 39185 : NULLIFY (tmp%shell_vel)
198 39185 : NULLIFY (tmp%core_pos)
199 39185 : NULLIFY (tmp%core_vel)
200 39185 : NULLIFY (tmp%itimes)
201 :
202 : ! *** Allocate work storage for positions and velocities ***
203 117555 : ALLOCATE (tmp%pos(3, nparticle))
204 :
205 117555 : ALLOCATE (tmp%vel(3, nparticle))
206 :
207 22618145 : tmp%pos(:, :) = 0.0_dp
208 22618145 : tmp%vel(:, :) = 0.0_dp
209 :
210 39185 : 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 156740 : tmp%arg_r = 0.0_dp
227 156740 : tmp%arg_v = 0.0_dp
228 509405 : tmp%u = 0.0_dp
229 156740 : tmp%e_val = 0.0_dp
230 39185 : tmp%max_vel = 0.0_dp
231 39185 : tmp%max_vel_sc = 0.0_dp
232 39185 : tmp%max_dvel = 0.0_dp
233 39185 : tmp%max_dvel_sc = 0.0_dp
234 156740 : tmp%scale_r = 1.0_dp
235 156740 : tmp%scale_v = 1.0_dp
236 156740 : tmp%poly_r = 1.0_dp
237 156740 : tmp%poly_v = 1.0_dp
238 :
239 39185 : CALL get_md_env(md_env=md_env, itimes=tmp%itimes)
240 :
241 39185 : 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 79618 : 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 79618 : CPASSERT(ASSOCIATED(tmp))
272 79618 : my_pos = .FALSE.
273 79618 : my_vel = .FALSE.
274 79618 : my_deall = .TRUE.
275 79618 : IF (PRESENT(pos)) my_pos = pos
276 79618 : IF (PRESENT(vel)) my_vel = vel
277 79618 : IF (PRESENT(should_deall_vel)) my_deall = should_deall_vel
278 :
279 : ! *** Broadcast the new particle positions ***
280 79618 : IF (my_pos) THEN
281 39185 : CALL update_particle_set(particle_set, para_env, pos=tmp%pos)
282 39185 : DEALLOCATE (tmp%pos)
283 :
284 39185 : 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 79618 : IF (my_vel) THEN
294 40433 : CALL update_particle_set(particle_set, para_env, vel=tmp%vel)
295 40433 : 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 40433 : IF (my_deall) THEN
300 39185 : DEALLOCATE (tmp%vel)
301 39185 : IF (shell_adiabatic) THEN
302 2290 : DEALLOCATE (tmp%shell_vel)
303 2290 : DEALLOCATE (tmp%core_vel)
304 : END IF
305 39185 : CPASSERT(.NOT. ASSOCIATED(tmp%pos))
306 39185 : CPASSERT(.NOT. ASSOCIATED(tmp%shell_pos))
307 39185 : CPASSERT(.NOT. ASSOCIATED(tmp%core_pos))
308 39185 : DEALLOCATE (tmp)
309 : END IF
310 : END IF
311 :
312 79618 : 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 1348 : 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 1348 : nparticle_kind = SIZE(atomic_kind_set)
411 : SELECT CASE (char)
412 : CASE ('F') ! forward assigning the old
413 1140 : DO iparticle_kind = 1, nparticle_kind
414 762 : nparticle_local = local_particles%n_el(iparticle_kind)
415 26167 : DO iparticle_local = 1, nparticle_local
416 25027 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
417 100870 : DO idim = 1, 3
418 75081 : old%v(iparticle, idim) = particle_set(iparticle)%v(idim)
419 100108 : old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
420 : END DO
421 : END DO
422 : END DO
423 1234 : old%eps(:, :) = npt(:, :)%eps
424 1234 : old%veps(:, :) = npt(:, :)%v
425 9450 : old%h(:, :) = cell%hmat(:, :)
426 : CASE ('B') ! back assigning the original variables
427 2930 : DO iparticle_kind = 1, nparticle_kind
428 1960 : nparticle_local = local_particles%n_el(iparticle_kind)
429 81660 : DO iparticle_local = 1, nparticle_local
430 78730 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
431 316880 : DO idim = 1, 3
432 236190 : particle_set(iparticle)%v(idim) = old%v(iparticle, idim)
433 314920 : particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
434 : END DO
435 : END DO
436 : END DO
437 3110 : npt(:, :)%eps = old%eps(:, :)
438 3110 : npt(:, :)%v = old%veps(:, :)
439 25598 : cell%hmat(:, :) = old%h(:, :)
440 : END SELECT
441 :
442 1348 : 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 1436 : 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 1436 : nparticle_kind = SIZE(atomic_kind_set)
473 : SELECT CASE (char)
474 : CASE ('F') ! forward assigning the old
475 1140 : DO iparticle_kind = 1, nparticle_kind
476 762 : nparticle_local = local_particles%n_el(iparticle_kind)
477 26167 : DO iparticle_local = 1, nparticle_local
478 25027 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
479 100870 : DO idim = 1, 3
480 75081 : old%v(iparticle, idim) = vel(idim, iparticle)
481 100108 : old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
482 : END DO
483 : END DO
484 : END DO
485 1234 : old%eps(:, :) = npt(:, :)%eps
486 1234 : old%veps(:, :) = npt(:, :)%v
487 9450 : old%h(:, :) = cell%hmat(:, :)
488 : CASE ('B') ! back assigning the original variables
489 3192 : DO iparticle_kind = 1, nparticle_kind
490 2134 : nparticle_local = local_particles%n_el(iparticle_kind)
491 76355 : DO iparticle_local = 1, nparticle_local
492 73163 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
493 294786 : DO idim = 1, 3
494 219489 : vel(idim, iparticle) = old%v(iparticle, idim)
495 292652 : particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
496 : END DO
497 : END DO
498 : END DO
499 3374 : npt(:, :)%eps = old%eps(:, :)
500 3374 : npt(:, :)%v = old%veps(:, :)
501 27886 : cell%hmat(:, :) = old%h(:, :)
502 : END SELECT
503 :
504 1436 : 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 1058 : SUBROUTINE rattle_roll_setup(old, gci, atomic_kind_set, particle_set, local_particles, &
714 1058 : molecule_kind_set, molecule_set, local_molecules, vel, dt, cell, npt, simpar, virial, &
715 1058 : 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 1058 : IF (first) THEN
744 : CALL update_pv(gci, simpar, atomic_kind_set, vel, particle_set, &
745 : local_molecules, molecule_set, molecule_kind_set, &
746 378 : local_particles, kin, pv_kin, virial, para_env)
747 378 : CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
748 :
749 : END IF
750 1058 : first = .FALSE.
751 :
752 : ! assigning local variable
753 1058 : SELECT CASE (simpar%ensemble)
754 : CASE (npt_i_ensemble, npt_ia_ensemble)
755 13494 : npt_loc(:, :)%v = 0.0_dp
756 13494 : npt_loc(:, :)%mass = 0.0_dp
757 1038 : npt_loc(1, 1)%v = npt(1, 1)%v
758 1038 : npt_loc(2, 2)%v = npt(1, 1)%v
759 1038 : 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 1038 : npt_loc(3, 3)%mass = npt(1, 1)%mass
763 : CASE (npt_f_ensemble)
764 1298 : npt_loc = npt
765 : END SELECT
766 :
767 : ! resetting
768 1058 : 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 27488 : para_env, u, cell, local_particles)
772 :
773 1058 : 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 2700 : 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 2700 : pv_kin = 0.0_dp
817 2700 : kin = 0.0_dp
818 8680 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
819 5980 : atomic_kind => atomic_kind_set(iparticle_kind)
820 5980 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
821 5980 : nparticle_local = local_particles%n_el(iparticle_kind)
822 291603 : DO iparticle_local = 1, nparticle_local
823 282923 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
824 1137672 : DO i = 1, 3
825 3395076 : DO j = 1, 3
826 : pv_kin(i, j) = pv_kin(i, j) + &
827 : mass*particle_set(iparticle)%v(i)* &
828 3395076 : particle_set(iparticle)%v(j)
829 : END DO
830 : kin = kin + mass*particle_set(iparticle)%v(i)* &
831 1131692 : particle_set(iparticle)%v(i)
832 : END DO
833 : END DO
834 : END DO
835 :
836 2700 : CALL int_group%sum(pv_kin)
837 2700 : CALL int_group%sum(kin)
838 :
839 : ! updating the constraint virial
840 2700 : IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
841 : molecule_set, &
842 : molecule_kind_set, &
843 : particle_set, virial, &
844 970 : int_group)
845 :
846 2700 : 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 3166 : SUBROUTINE update_pv_velocity(gci, simpar, atomic_kind_set, vel, particle_set, &
869 : local_molecules, molecule_set, molecule_kind_set, &
870 3166 : 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 41158 : pv_kin = 0.0_dp
892 3166 : kin = 0._dp
893 10082 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
894 6916 : atomic_kind => atomic_kind_set(iparticle_kind)
895 6916 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
896 6916 : nparticle_local = local_particles%n_el(iparticle_kind)
897 312465 : DO iparticle_local = 1, nparticle_local
898 302383 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
899 1216448 : DO i = 1, 3
900 3628596 : DO j = 1, 3
901 : pv_kin(i, j) = pv_kin(i, j) + &
902 3628596 : mass*vel(i, iparticle)*vel(j, iparticle)
903 : END DO
904 1209532 : kin = kin + mass*vel(i, iparticle)*vel(i, iparticle)
905 : END DO
906 : END DO
907 : END DO
908 :
909 79150 : CALL int_group%sum(pv_kin)
910 3166 : CALL int_group%sum(kin)
911 :
912 : ! updating the constraint virial
913 3166 : IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
914 : molecule_set, &
915 : molecule_kind_set, &
916 : particle_set, virial, &
917 1436 : int_group)
918 :
919 3166 : 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 5866 : SUBROUTINE update_veps(box, npt, simpar, pv_kin, kin, virial, infree, virial_components)
936 :
937 : TYPE(cell_type), INTENT(IN) :: 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 : ! initialize unit
952 :
953 5866 : unit = 0.0_dp
954 5866 : unit(1, 1) = 1.0_dp
955 5866 : unit(2, 2) = 1.0_dp
956 5866 : unit(3, 3) = 1.0_dp
957 :
958 : IF (simpar%ensemble == npt_i_ensemble .OR. &
959 5866 : simpar%ensemble == npt_ia_ensemble .OR. &
960 : simpar%ensemble == npe_i_ensemble) THEN
961 : ! get force on barostat
962 : fdotr = 0.0_dp
963 15696 : DO ii = 1, 3
964 : fdotr = fdotr + virial%pv_virial(ii, ii) + &
965 15696 : virial%pv_constraint(ii, ii)
966 : END DO
967 :
968 : npt(:, :)%f = (1.0_dp + (3.0_dp*infree))*kin + fdotr - &
969 11772 : 3.0_dp*simpar%p_ext*box%deth
970 1942 : ELSEIF (simpar%ensemble == npt_f_ensemble .OR. &
971 : simpar%ensemble == npe_f_ensemble) THEN
972 : npt(:, :)%f = virial%pv_virial(:, :) + &
973 : pv_kin(:, :) + virial%pv_constraint(:, :) - &
974 : unit(:, :)*simpar%p_ext*box%deth + &
975 23686 : infree*kin*unit(:, :)
976 : IF (debug_isotropic_limit) THEN
977 : trace = npt(1, 1)%f + npt(2, 2)%f + npt(3, 3)%f
978 : trace = trace/3.0_dp
979 : npt(:, :)%f = trace*unit(:, :)
980 : END IF
981 120 : ELSEIF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
982 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
983 120 : v = box%deth
984 120 : vi = 1._dp/v
985 120 : v0 = simpar%v0
986 120 : v0i = 1._dp/v0
987 : ! orthorhombic box ONLY
988 : ! Chooses only the compressive solution
989 120 : IF (v < v0) THEN
990 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
991 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
992 : simpar%p0*v - simpar%v_shock*simpar%v_shock* &
993 0 : v*v0i*(1._dp - v*v0i) + infree*kin
994 : ELSE
995 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
996 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
997 120 : simpar%p0*v + infree*kin
998 : END IF
999 : IF (debug_uniaxial_limit) THEN
1000 : ! orthorhombic box ONLY
1001 : npt(1, 1)%f = virial%pv_virial(1, 1) + &
1002 : pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
1003 : simpar%p0*box%deth + infree*kin
1004 : END IF
1005 : END IF
1006 :
1007 : ! update barostat velocities
1008 : npt(:, :)%v = npt(:, :)%v + &
1009 35818 : 0.5_dp*simpar%dt*npt(:, :)%f/npt(:, :)%mass
1010 :
1011 : ! Screen the dynamics of the barostat according user request
1012 5866 : IF (PRESENT(virial_components)) THEN
1013 1812 : SELECT CASE (virial_components)
1014 : CASE (do_clv_xyz)
1015 : ! Do nothing..
1016 : CASE (do_clv_x)
1017 0 : npt(1, 2)%v = 0.0_dp
1018 0 : npt(1, 3)%v = 0.0_dp
1019 0 : npt(1, 2)%f = 0.0_dp
1020 0 : npt(1, 3)%f = 0.0_dp
1021 0 : npt(2, 1)%v = 0.0_dp
1022 0 : npt(2, 2)%v = 0.0_dp
1023 0 : npt(2, 3)%v = 0.0_dp
1024 0 : npt(2, 1)%f = 0.0_dp
1025 0 : npt(2, 2)%f = 0.0_dp
1026 0 : npt(2, 3)%f = 0.0_dp
1027 0 : npt(3, 1)%v = 0.0_dp
1028 0 : npt(3, 2)%v = 0.0_dp
1029 0 : npt(3, 3)%v = 0.0_dp
1030 0 : npt(3, 1)%f = 0.0_dp
1031 0 : npt(3, 2)%f = 0.0_dp
1032 0 : npt(3, 3)%f = 0.0_dp
1033 : CASE (do_clv_y)
1034 0 : npt(1, 1)%v = 0.0_dp
1035 0 : npt(1, 2)%v = 0.0_dp
1036 0 : npt(1, 3)%v = 0.0_dp
1037 0 : npt(1, 1)%f = 0.0_dp
1038 0 : npt(1, 2)%f = 0.0_dp
1039 0 : npt(1, 3)%f = 0.0_dp
1040 0 : npt(2, 1)%v = 0.0_dp
1041 0 : npt(2, 3)%v = 0.0_dp
1042 0 : npt(2, 1)%f = 0.0_dp
1043 0 : npt(2, 3)%f = 0.0_dp
1044 0 : npt(3, 1)%v = 0.0_dp
1045 0 : npt(3, 2)%v = 0.0_dp
1046 0 : npt(3, 3)%v = 0.0_dp
1047 0 : npt(3, 1)%f = 0.0_dp
1048 0 : npt(3, 2)%f = 0.0_dp
1049 0 : npt(3, 3)%f = 0.0_dp
1050 : CASE (do_clv_z)
1051 80 : npt(1, 1)%v = 0.0_dp
1052 80 : npt(1, 2)%v = 0.0_dp
1053 80 : npt(1, 3)%v = 0.0_dp
1054 80 : npt(1, 1)%f = 0.0_dp
1055 80 : npt(1, 2)%f = 0.0_dp
1056 80 : npt(1, 3)%f = 0.0_dp
1057 80 : npt(2, 1)%v = 0.0_dp
1058 80 : npt(2, 2)%v = 0.0_dp
1059 80 : npt(2, 3)%v = 0.0_dp
1060 80 : npt(2, 1)%f = 0.0_dp
1061 80 : npt(2, 2)%f = 0.0_dp
1062 80 : npt(2, 3)%f = 0.0_dp
1063 80 : npt(3, 1)%v = 0.0_dp
1064 80 : npt(3, 2)%v = 0.0_dp
1065 80 : npt(3, 1)%f = 0.0_dp
1066 80 : npt(3, 2)%f = 0.0_dp
1067 : CASE (do_clv_xy)
1068 0 : npt(1, 3)%v = 0.0_dp
1069 0 : npt(1, 3)%f = 0.0_dp
1070 0 : npt(2, 3)%v = 0.0_dp
1071 0 : npt(2, 3)%f = 0.0_dp
1072 0 : npt(3, 1)%v = 0.0_dp
1073 0 : npt(3, 2)%v = 0.0_dp
1074 0 : npt(3, 3)%v = 0.0_dp
1075 0 : npt(3, 1)%f = 0.0_dp
1076 0 : npt(3, 2)%f = 0.0_dp
1077 0 : npt(3, 3)%f = 0.0_dp
1078 : CASE (do_clv_xz)
1079 0 : npt(1, 2)%v = 0.0_dp
1080 0 : npt(1, 2)%f = 0.0_dp
1081 0 : npt(2, 1)%v = 0.0_dp
1082 0 : npt(2, 2)%v = 0.0_dp
1083 0 : npt(2, 3)%v = 0.0_dp
1084 0 : npt(2, 1)%f = 0.0_dp
1085 0 : npt(2, 2)%f = 0.0_dp
1086 0 : npt(2, 3)%f = 0.0_dp
1087 0 : npt(3, 2)%v = 0.0_dp
1088 0 : npt(3, 2)%f = 0.0_dp
1089 : CASE (do_clv_yz)
1090 0 : npt(1, 1)%v = 0.0_dp
1091 0 : npt(1, 2)%v = 0.0_dp
1092 0 : npt(1, 3)%v = 0.0_dp
1093 0 : npt(1, 1)%f = 0.0_dp
1094 0 : npt(1, 2)%f = 0.0_dp
1095 0 : npt(1, 3)%f = 0.0_dp
1096 0 : npt(2, 1)%v = 0.0_dp
1097 0 : npt(2, 1)%f = 0.0_dp
1098 0 : npt(3, 1)%v = 0.0_dp
1099 1812 : npt(3, 1)%f = 0.0_dp
1100 : END SELECT
1101 : END IF
1102 :
1103 5866 : END SUBROUTINE update_veps
1104 :
1105 : ! **************************************************************************************************
1106 : !> \brief First half of the velocity-verlet algorithm : update velocity by half
1107 : !> step and positions by full step
1108 : !> \param tmp ...
1109 : !> \param atomic_kind_set ...
1110 : !> \param local_particles ...
1111 : !> \param particle_set ...
1112 : !> \param core_particle_set ...
1113 : !> \param shell_particle_set ...
1114 : !> \param nparticle_kind ...
1115 : !> \param shell_adiabatic ...
1116 : !> \param dt ...
1117 : !> \param u ...
1118 : !> \param lfixd_list ...
1119 : !> \par History
1120 : !> none
1121 : !> \author MI (February 2008)
1122 : ! **************************************************************************************************
1123 40253 : SUBROUTINE vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1124 : core_particle_set, shell_particle_set, nparticle_kind, &
1125 40253 : shell_adiabatic, dt, u, lfixd_list)
1126 :
1127 : TYPE(tmp_variables_type), POINTER :: tmp
1128 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1129 : TYPE(distribution_1d_type), POINTER :: local_particles
1130 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1131 : shell_particle_set
1132 : INTEGER, INTENT(IN) :: nparticle_kind
1133 : LOGICAL, INTENT(IN) :: shell_adiabatic
1134 : REAL(KIND=dp) :: dt
1135 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
1136 : TYPE(local_fixd_constraint_type), DIMENSION(:), &
1137 : OPTIONAL :: lfixd_list
1138 :
1139 : INTEGER :: iparticle, iparticle_kind, &
1140 : iparticle_local, nparticle_local, &
1141 : shell_index
1142 : LOGICAL :: is_shell
1143 : REAL(KIND=dp) :: dm, dmc, dms, dsc(3), dvsc(3), &
1144 : fac_massc, fac_masss, mass
1145 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1146 : TYPE(shell_kind_type), POINTER :: shell
1147 :
1148 40253 : NULLIFY (atomic_kind, shell)
1149 40253 : tmp%max_vel = 0.0_dp
1150 40253 : tmp%max_vel_sc = 0.0_dp
1151 40253 : tmp%max_dvel = 0.0_dp
1152 40253 : tmp%max_dvel_sc = 0.0_dp
1153 40253 : tmp%max_dr = 0.0_dp
1154 40253 : tmp%max_dsc = 0.0_dp
1155 40253 : dsc = 0.0_dp
1156 40253 : dvsc = 0.0_dp
1157 :
1158 : ! *** Velocity Verlet (first part) ***
1159 145546 : DO iparticle_kind = 1, nparticle_kind
1160 105293 : atomic_kind => atomic_kind_set(iparticle_kind)
1161 105293 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell)
1162 145546 : IF (mass /= 0.0_dp) THEN
1163 104889 : dm = 0.5_dp*dt/mass
1164 104889 : IF (is_shell .AND. shell_adiabatic) THEN
1165 5412 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1166 5412 : dms = 0.5_dp*dt/shell%mass_shell
1167 5412 : dmc = 0.5_dp*dt/shell%mass_core
1168 5412 : fac_masss = shell%mass_shell/mass
1169 5412 : fac_massc = shell%mass_core/mass
1170 5412 : nparticle_local = local_particles%n_el(iparticle_kind)
1171 5412 : IF (PRESENT(u)) THEN
1172 54620 : DO iparticle_local = 1, nparticle_local
1173 52704 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1174 52704 : shell_index = particle_set(iparticle)%shell_index
1175 : ! Transform positions and velocities and forces of the shell
1176 : CALL transform_first(shell_particle_set, tmp%shell_pos, tmp%shell_vel, &
1177 52704 : shell_index, u, dms, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1178 :
1179 : ! Transform positions and velocities and forces of the core
1180 : CALL transform_first(core_particle_set, tmp%core_pos, tmp%core_vel, &
1181 52704 : shell_index, u, dmc, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1182 :
1183 : ! Derive velocities and positions of the COM
1184 : tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
1185 210816 : fac_massc*tmp%core_vel(:, shell_index)
1186 : tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1187 210816 : fac_massc*tmp%core_pos(:, shell_index)
1188 :
1189 : tmp%max_vel = MAX(tmp%max_vel, ABS(tmp%vel(1, iparticle)), &
1190 52704 : ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1191 : tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
1192 : ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
1193 : ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
1194 52704 : ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
1195 : tmp%max_dr = MAX(tmp%max_dr, &
1196 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1197 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1198 52704 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1199 : tmp%max_dvel = MAX(tmp%max_dvel, &
1200 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1201 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1202 52704 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1203 :
1204 : dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1205 210816 : shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1206 52704 : tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
1207 : dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
1208 210816 : shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1209 54620 : tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
1210 : END DO ! iparticle_local
1211 : ELSE
1212 88530 : DO iparticle_local = 1, nparticle_local
1213 85034 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1214 85034 : shell_index = particle_set(iparticle)%shell_index
1215 : tmp%shell_vel(:, shell_index) = &
1216 : shell_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
1217 680272 : tmp%scale_v(:)*tmp%poly_v(:)*dms*shell_particle_set(shell_index)%f(:)
1218 : tmp%shell_pos(:, shell_index) = &
1219 : shell_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
1220 680272 : tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%shell_vel(:, shell_index)
1221 : tmp%core_vel(:, shell_index) = &
1222 : core_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
1223 680272 : tmp%scale_v(:)*tmp%poly_v(:)*dmc*core_particle_set(shell_index)%f(:)
1224 : tmp%core_pos(:, shell_index) = &
1225 : core_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
1226 680272 : tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%core_vel(:, shell_index)
1227 :
1228 : tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
1229 340136 : fac_massc*tmp%core_vel(:, shell_index)
1230 : tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1231 340136 : fac_massc*tmp%core_pos(:, shell_index)
1232 : tmp%max_vel = MAX(tmp%max_vel, &
1233 85034 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1234 : tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
1235 : ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
1236 : ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
1237 85034 : ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
1238 : tmp%max_dr = MAX(tmp%max_dr, &
1239 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1240 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1241 85034 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1242 : tmp%max_dvel = MAX(tmp%max_dvel, &
1243 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1244 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1245 85034 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1246 : dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1247 340136 : shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1248 85034 : tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
1249 : dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
1250 340136 : shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1251 88530 : tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
1252 : END DO ! iparticle_local
1253 : END IF
1254 : ELSE
1255 99477 : nparticle_local = local_particles%n_el(iparticle_kind)
1256 99477 : IF (PRESENT(u)) THEN
1257 44038 : DO iparticle_local = 1, nparticle_local
1258 43412 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1259 : ! Transform positions and velocities and forces
1260 : CALL transform_first(particle_set, tmp%pos, tmp%vel, &
1261 43412 : iparticle, u, dm, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
1262 :
1263 : tmp%max_vel = MAX(tmp%max_vel, &
1264 43412 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1265 : tmp%max_dr = MAX(tmp%max_dr, ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1266 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1267 43412 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1268 : tmp%max_dvel = MAX(tmp%max_dvel, &
1269 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1270 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1271 44038 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1272 : END DO ! iparticle_local
1273 : ELSE
1274 2830010 : DO iparticle_local = 1, nparticle_local
1275 2731159 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1276 : tmp%vel(1, iparticle) = &
1277 : particle_set(iparticle)%v(1)*tmp%scale_v(1)*tmp%scale_v(1) + &
1278 2731159 : tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
1279 : tmp%vel(2, iparticle) = &
1280 : particle_set(iparticle)%v(2)*tmp%scale_v(2)*tmp%scale_v(2) + &
1281 2731159 : tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
1282 : tmp%vel(3, iparticle) = &
1283 : particle_set(iparticle)%v(3)*tmp%scale_v(3)*tmp%scale_v(3) + &
1284 2731159 : tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
1285 : tmp%pos(1, iparticle) = &
1286 : particle_set(iparticle)%r(1)*tmp%scale_r(1)*tmp%scale_r(1) + &
1287 2731159 : tmp%scale_r(1)*tmp%poly_r(1)*dt*tmp%vel(1, iparticle)
1288 : tmp%pos(2, iparticle) = &
1289 : particle_set(iparticle)%r(2)*tmp%scale_r(2)*tmp%scale_r(2) + &
1290 2731159 : tmp%scale_r(2)*tmp%poly_r(2)*dt*tmp%vel(2, iparticle)
1291 : tmp%pos(3, iparticle) = &
1292 : particle_set(iparticle)%r(3)*tmp%scale_r(3)*tmp%scale_r(3) + &
1293 2731159 : tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
1294 :
1295 : ! overwrite positions of frozen particles
1296 2731159 : IF (PRESENT(lfixd_list)) THEN
1297 251570 : IF (ANY(lfixd_list(:)%ifixd_index == iparticle)) THEN
1298 200 : tmp%pos(1, iparticle) = particle_set(iparticle)%r(1)
1299 200 : tmp%pos(2, iparticle) = particle_set(iparticle)%r(2)
1300 200 : tmp%pos(3, iparticle) = particle_set(iparticle)%r(3)
1301 : END IF
1302 : END IF
1303 :
1304 : tmp%max_vel = MAX(tmp%max_vel, &
1305 2731159 : ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
1306 : tmp%max_dr = MAX(tmp%max_dr, &
1307 : ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
1308 : ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
1309 2731159 : ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
1310 : tmp%max_dvel = MAX(tmp%max_dvel, &
1311 : ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
1312 : ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
1313 2830010 : ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1314 : END DO
1315 : END IF
1316 : END IF
1317 : END IF
1318 : END DO
1319 40253 : END SUBROUTINE vv_first
1320 :
1321 : ! **************************************************************************************************
1322 : !> \brief Second half of the velocity-verlet algorithm : update velocity by half
1323 : !> step using the new forces
1324 : !> \param tmp ...
1325 : !> \param atomic_kind_set ...
1326 : !> \param local_particles ...
1327 : !> \param particle_set ...
1328 : !> \param core_particle_set ...
1329 : !> \param shell_particle_set ...
1330 : !> \param nparticle_kind ...
1331 : !> \param shell_adiabatic ...
1332 : !> \param dt ...
1333 : !> \param u ...
1334 : !> \par History
1335 : !> none
1336 : !> \author MI (February 2008)
1337 : ! **************************************************************************************************
1338 40433 : SUBROUTINE vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
1339 : core_particle_set, shell_particle_set, nparticle_kind, &
1340 : shell_adiabatic, dt, u)
1341 :
1342 : TYPE(tmp_variables_type), POINTER :: tmp
1343 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1344 : TYPE(distribution_1d_type), POINTER :: local_particles
1345 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1346 : shell_particle_set
1347 : INTEGER, INTENT(IN) :: nparticle_kind
1348 : LOGICAL, INTENT(IN) :: shell_adiabatic
1349 : REAL(KIND=dp) :: dt
1350 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
1351 :
1352 : INTEGER :: iparticle, iparticle_kind, &
1353 : iparticle_local, nparticle_local, &
1354 : shell_index
1355 : LOGICAL :: is_shell
1356 : REAL(KIND=dp) :: dm, dmc, dms, fac_massc, fac_masss, mass
1357 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1358 : TYPE(shell_kind_type), POINTER :: shell
1359 :
1360 40433 : NULLIFY (atomic_kind, shell)
1361 :
1362 : ! *** Velocity Verlet (second part) ***
1363 145690 : DO iparticle_kind = 1, nparticle_kind
1364 105257 : atomic_kind => atomic_kind_set(iparticle_kind)
1365 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
1366 105257 : shell_active=is_shell)
1367 145690 : IF (mass /= 0.0_dp) THEN
1368 104853 : dm = 0.5_dp*dt/mass
1369 104853 : IF (is_shell .AND. shell_adiabatic) THEN
1370 4520 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1371 4520 : dms = 0.5_dp*dt/shell%mass_shell
1372 4520 : dmc = 0.5_dp*dt/shell%mass_core
1373 4520 : fac_masss = shell%mass_shell/mass
1374 4520 : fac_massc = shell%mass_core/mass
1375 4520 : nparticle_local = local_particles%n_el(iparticle_kind)
1376 4520 : IF (PRESENT(u)) THEN
1377 35860 : DO iparticle_local = 1, nparticle_local
1378 34560 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1379 34560 : shell_index = particle_set(iparticle)%shell_index
1380 : ! Transform velocities and forces shell
1381 : CALL transform_second(shell_particle_set, tmp%shell_vel, shell_index, &
1382 34560 : u, dms, tmp%poly_v, tmp%scale_v)
1383 :
1384 : ! Transform velocities and forces core
1385 : CALL transform_second(core_particle_set, tmp%core_vel, shell_index, &
1386 34560 : u, dmc, tmp%poly_v, tmp%scale_v)
1387 :
1388 : ! Derive velocties of the COM
1389 : tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
1390 34560 : fac_massc*tmp%core_vel(1, shell_index)
1391 : tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1392 34560 : fac_massc*tmp%core_vel(2, shell_index)
1393 : tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1394 35860 : fac_massc*tmp%core_vel(3, shell_index)
1395 : END DO ! iparticle_local
1396 : ELSE
1397 81630 : DO iparticle_local = 1, nparticle_local
1398 78410 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1399 78410 : shell_index = particle_set(iparticle)%shell_index
1400 : tmp%shell_vel(1, shell_index) = &
1401 : tmp%shell_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
1402 78410 : tmp%scale_v(1)*tmp%poly_v(1)*dms*shell_particle_set(shell_index)%f(1)
1403 : tmp%shell_vel(2, shell_index) = &
1404 : tmp%shell_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
1405 78410 : tmp%scale_v(2)*tmp%poly_v(2)*dms*shell_particle_set(shell_index)%f(2)
1406 : tmp%shell_vel(3, shell_index) = &
1407 : tmp%shell_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
1408 78410 : tmp%scale_v(3)*tmp%poly_v(3)*dms*shell_particle_set(shell_index)%f(3)
1409 : tmp%core_vel(1, shell_index) = &
1410 : tmp%core_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
1411 78410 : tmp%scale_v(1)*tmp%poly_v(1)*dmc*core_particle_set(shell_index)%f(1)
1412 : tmp%core_vel(2, shell_index) = &
1413 : tmp%core_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
1414 78410 : tmp%scale_v(2)*tmp%poly_v(2)*dmc*core_particle_set(shell_index)%f(2)
1415 : tmp%core_vel(3, shell_index) = &
1416 : tmp%core_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
1417 78410 : tmp%scale_v(3)*tmp%poly_v(3)*dmc*core_particle_set(shell_index)%f(3)
1418 :
1419 : tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
1420 78410 : fac_massc*tmp%core_vel(1, shell_index)
1421 : tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1422 78410 : fac_massc*tmp%core_vel(2, shell_index)
1423 : tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1424 81630 : fac_massc*tmp%core_vel(3, shell_index)
1425 : END DO ! iparticle_local
1426 : END IF
1427 : ELSE
1428 100333 : nparticle_local = local_particles%n_el(iparticle_kind)
1429 100333 : IF (PRESENT(u)) THEN
1430 44038 : DO iparticle_local = 1, nparticle_local
1431 43412 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1432 : CALL transform_second(particle_set, tmp%vel, iparticle, &
1433 44038 : u, dm, tmp%poly_v, tmp%scale_v)
1434 : END DO ! iparticle_local
1435 : ELSE
1436 2780861 : DO iparticle_local = 1, nparticle_local
1437 2681154 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1438 : tmp%vel(1, iparticle) = &
1439 : tmp%vel(1, iparticle)*tmp%scale_v(1)*tmp%scale_v(1) + &
1440 2681154 : tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
1441 : tmp%vel(2, iparticle) = &
1442 : tmp%vel(2, iparticle)*tmp%scale_v(2)*tmp%scale_v(2) + &
1443 2681154 : tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
1444 : tmp%vel(3, iparticle) = &
1445 : tmp%vel(3, iparticle)*tmp%scale_v(3)*tmp%scale_v(3) + &
1446 2780861 : tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
1447 : END DO
1448 : END IF
1449 : END IF
1450 : END IF
1451 : END DO
1452 :
1453 40433 : END SUBROUTINE vv_second
1454 :
1455 : ! **************************************************************************************************
1456 : !> \brief Transform position and velocities for the first half of the
1457 : !> Velocity-Verlet integration in the npt_f ensemble
1458 : !> \param particle_set ...
1459 : !> \param pos ...
1460 : !> \param vel ...
1461 : !> \param index ...
1462 : !> \param u ...
1463 : !> \param dm ...
1464 : !> \param dt ...
1465 : !> \param poly_v ...
1466 : !> \param poly_r ...
1467 : !> \param scale_v ...
1468 : !> \param scale_r ...
1469 : !> \par History
1470 : !> none
1471 : !> \author MI (February 2008)
1472 : ! **************************************************************************************************
1473 148820 : SUBROUTINE transform_first(particle_set, pos, vel, index, u, dm, dt, poly_v, &
1474 : poly_r, scale_v, scale_r)
1475 :
1476 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1477 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pos, vel
1478 : INTEGER, INTENT(IN) :: index
1479 : REAL(KIND=dp), INTENT(IN) :: u(3, 3), dm, dt, poly_v(3), poly_r(3), &
1480 : scale_v(3), scale_r(3)
1481 :
1482 : REAL(KIND=dp), DIMENSION(3) :: uf, ur, uv
1483 :
1484 148820 : ur = MATMUL(TRANSPOSE(u), particle_set(index)%r(:))
1485 148820 : uv = MATMUL(TRANSPOSE(u), particle_set(index)%v(:))
1486 148820 : uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
1487 :
1488 148820 : uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1489 148820 : uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1490 148820 : uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1491 :
1492 148820 : ur(1) = ur(1)*scale_r(1)*scale_r(1) + uv(1)*scale_r(1)*poly_r(1)*dt
1493 148820 : ur(2) = ur(2)*scale_r(2)*scale_r(2) + uv(2)*scale_r(2)*poly_r(2)*dt
1494 148820 : ur(3) = ur(3)*scale_r(3)*scale_r(3) + uv(3)*scale_r(3)*poly_r(3)*dt
1495 :
1496 2529940 : pos(:, index) = MATMUL(u, ur)
1497 2529940 : vel(:, index) = MATMUL(u, uv)
1498 :
1499 148820 : END SUBROUTINE transform_first
1500 :
1501 : ! **************************************************************************************************
1502 : !> \brief Transform position and velocities for the second half of the
1503 : !> Velocity-Verlet integration in the npt_f ensemble
1504 : !> \param particle_set ...
1505 : !> \param vel ...
1506 : !> \param index ...
1507 : !> \param u ...
1508 : !> \param dm ...
1509 : !> \param poly_v ...
1510 : !> \param scale_v ...
1511 : !> \par History
1512 : !> none
1513 : !> \author MI (February 2008)
1514 : ! **************************************************************************************************
1515 112532 : SUBROUTINE transform_second(particle_set, vel, index, u, dm, poly_v, scale_v)
1516 :
1517 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1518 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vel
1519 : INTEGER, INTENT(IN) :: index
1520 : REAL(KIND=dp), INTENT(IN) :: u(3, 3), dm, poly_v(3), scale_v(3)
1521 :
1522 : REAL(KIND=dp), DIMENSION(3) :: uf, uv
1523 :
1524 112532 : uv = MATMUL(TRANSPOSE(u), vel(:, index))
1525 112532 : uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
1526 :
1527 112532 : uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1528 112532 : uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1529 112532 : uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1530 :
1531 1913044 : vel(:, index) = MATMUL(u, uv)
1532 :
1533 112532 : END SUBROUTINE transform_second
1534 :
1535 : ! **************************************************************************************************
1536 : !> \brief Compute the timestep rescaling factor
1537 : !> \param md_env ...
1538 : !> \param tmp ...
1539 : !> \param dt ...
1540 : !> \param simpar ...
1541 : !> \param para_env ...
1542 : !> \param atomic_kind_set ...
1543 : !> \param local_particles ...
1544 : !> \param particle_set ...
1545 : !> \param core_particle_set ...
1546 : !> \param shell_particle_set ...
1547 : !> \param nparticle_kind ...
1548 : !> \param shell_adiabatic ...
1549 : !> \param npt ...
1550 : !> \par History
1551 : !> none
1552 : !> \author MI (October 2008)
1553 : ! **************************************************************************************************
1554 480 : SUBROUTINE variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, local_particles, &
1555 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1556 :
1557 : TYPE(md_environment_type), POINTER :: md_env
1558 : TYPE(tmp_variables_type), POINTER :: tmp
1559 : REAL(KIND=dp), INTENT(INOUT) :: dt
1560 : TYPE(simpar_type), POINTER :: simpar
1561 : TYPE(mp_para_env_type), POINTER :: para_env
1562 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1563 : TYPE(distribution_1d_type), POINTER :: local_particles
1564 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1565 : shell_particle_set
1566 : INTEGER, INTENT(IN) :: nparticle_kind
1567 : LOGICAL, INTENT(IN) :: shell_adiabatic
1568 : TYPE(npt_info_type), OPTIONAL, POINTER :: npt(:, :)
1569 :
1570 : INTEGER, SAVE :: itime = 0
1571 : REAL(KIND=dp) :: dt_fact, dt_fact_f, dt_fact_old, &
1572 : dt_fact_sc, dt_fact_v, dt_old
1573 : REAL(KIND=dp), POINTER :: time
1574 : TYPE(thermostats_type), POINTER :: thermostats
1575 :
1576 480 : dt_fact = 1.0_dp
1577 480 : dt_fact_sc = 1.0_dp
1578 480 : dt_fact_f = 1.0_dp
1579 480 : dt_fact_v = 1.0_dp
1580 480 : dt_old = dt
1581 480 : dt_fact_old = simpar%dt_fact
1582 480 : simpar%dt_fact = 1.0_dp
1583 480 : NULLIFY (thermostats)
1584 :
1585 480 : itime = itime + 1
1586 480 : CALL para_env%max(tmp%max_dr)
1587 480 : IF (tmp%max_dr > simpar%dr_tol) THEN
1588 278 : CALL para_env%max(tmp%max_dvel)
1589 278 : dt_fact = SQRT(simpar%dr_tol*dt/tmp%max_dvel)/dt
1590 : END IF
1591 :
1592 480 : IF (shell_adiabatic) THEN
1593 440 : CALL para_env%max(tmp%max_dsc)
1594 440 : IF (tmp%max_dsc > simpar%dsc_tol) THEN
1595 118 : CALL para_env%max(tmp%max_dvel_sc)
1596 118 : dt_fact_sc = SQRT(simpar%dsc_tol*dt/tmp%max_dvel_sc)/dt
1597 : END IF
1598 : END IF
1599 :
1600 480 : dt_fact_f = MIN(dt_fact, dt_fact_sc, 1.0_dp)
1601 480 : IF (dt_fact_f < 1.0_dp) THEN
1602 100 : IF (dt_fact_f < 0.1_dp) dt_fact_f = 0.1_dp
1603 : ! repeat the first vv half-step with the rescaled time step
1604 100 : dt = dt*dt_fact_f
1605 100 : simpar%dt_fact = dt_fact_f
1606 : CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1607 100 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1608 : END IF
1609 :
1610 480 : dt_fact = 1.0_dp
1611 480 : dt_fact_sc = 1.0_dp
1612 480 : CALL para_env%max(tmp%max_dr)
1613 480 : IF (tmp%max_dr > simpar%dr_tol) THEN
1614 278 : CALL para_env%max(tmp%max_vel)
1615 278 : dt_fact = simpar%dr_tol/tmp%max_vel/dt
1616 : END IF
1617 :
1618 480 : IF (shell_adiabatic) THEN
1619 440 : CALL para_env%max(tmp%max_dsc)
1620 440 : IF (tmp%max_dsc > simpar%dsc_tol) THEN
1621 116 : CALL para_env%max(tmp%max_vel_sc)
1622 116 : dt_fact_sc = simpar%dsc_tol/tmp%max_vel_sc/dt
1623 : END IF
1624 : END IF
1625 480 : dt_fact_v = MIN(dt_fact, dt_fact_sc, 1.0_dp)
1626 :
1627 480 : IF (dt_fact_v < 1.0_dp) THEN
1628 376 : IF (dt_fact_v < 0.1_dp) dt_fact_v = 0.1_dp
1629 : ! repeat the first vv half-step with the rescaled time step
1630 376 : dt = dt*dt_fact_v
1631 376 : simpar%dt_fact = dt_fact_f*dt_fact_v
1632 : CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1633 554 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1634 : END IF
1635 :
1636 480 : simpar%dt_fact = dt_fact_f*dt_fact_v
1637 480 : IF (simpar%dt_fact < 1.0_dp) THEN
1638 378 : CALL get_md_env(md_env, t=time, thermostats=thermostats)
1639 378 : time = time - dt_old + dt_old*simpar%dt_fact
1640 378 : IF (ASSOCIATED(thermostats)) THEN
1641 240 : CALL set_thermostats(thermostats, dt_fact=simpar%dt_fact)
1642 : END IF
1643 : END IF
1644 :
1645 480 : END SUBROUTINE variable_timestep
1646 :
1647 : ! **************************************************************************************************
1648 : !> \brief Repeat the first step of velocity verlet with the rescaled time step
1649 : !> \param tmp ...
1650 : !> \param dt ...
1651 : !> \param simpar ...
1652 : !> \param atomic_kind_set ...
1653 : !> \param local_particles ...
1654 : !> \param particle_set ...
1655 : !> \param core_particle_set ...
1656 : !> \param shell_particle_set ...
1657 : !> \param nparticle_kind ...
1658 : !> \param shell_adiabatic ...
1659 : !> \param npt ...
1660 : !> \par History
1661 : !> none
1662 : !> \author MI (October 2008)
1663 : !> I soliti ignoti
1664 : ! **************************************************************************************************
1665 476 : SUBROUTINE rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1666 : particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1667 :
1668 : TYPE(tmp_variables_type), POINTER :: tmp
1669 : REAL(KIND=dp), INTENT(IN) :: dt
1670 : TYPE(simpar_type), POINTER :: simpar
1671 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1672 : TYPE(distribution_1d_type), POINTER :: local_particles
1673 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
1674 : shell_particle_set
1675 : INTEGER, INTENT(IN) :: nparticle_kind
1676 : LOGICAL, INTENT(IN) :: shell_adiabatic
1677 : TYPE(npt_info_type), OPTIONAL, POINTER :: npt(:, :)
1678 :
1679 : REAL(KIND=dp), PARAMETER :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
1680 : e6 = e4/42.0_dp, e8 = e6/72.0_dp
1681 :
1682 : REAL(KIND=dp) :: arg_r(3), arg_v(3), e_val(3), infree, &
1683 : trvg, u(3, 3)
1684 :
1685 476 : infree = 1.0_dp/REAL(simpar%nfree, dp)
1686 1904 : arg_r = tmp%arg_r
1687 1904 : arg_v = tmp%arg_v
1688 6188 : u = tmp%u
1689 1904 : e_val = tmp%e_val
1690 :
1691 654 : SELECT CASE (simpar%ensemble)
1692 : CASE (nve_ensemble, nvt_ensemble)
1693 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1694 178 : core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
1695 : CASE (isokin_ensemble)
1696 0 : tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds)
1697 0 : tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt
1698 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1699 : core_particle_set, shell_particle_set, nparticle_kind, &
1700 0 : shell_adiabatic, dt)
1701 :
1702 : CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
1703 0 : arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1704 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
1705 0 : arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
1706 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
1707 0 : tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
1708 : tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
1709 0 : (1.0_dp + 3.0_dp*infree))
1710 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1711 : core_particle_set, shell_particle_set, nparticle_kind, &
1712 0 : shell_adiabatic, dt)
1713 :
1714 : CASE (npt_f_ensemble, npe_f_ensemble)
1715 298 : trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
1716 1192 : arg_r(:) = arg_r(:)*simpar%dt_fact*simpar%dt_fact
1717 1192 : tmp%poly_r = 1._dp + e2*arg_r + e4*arg_r*arg_r + e6*arg_r**3 + e8*arg_r**4
1718 1192 : tmp%scale_r(:) = EXP(0.5_dp*dt*e_val(:))
1719 1192 : arg_v(:) = arg_v(:)*simpar%dt_fact*simpar%dt_fact
1720 1192 : tmp%poly_v = 1.0_dp + e2*arg_v + e4*arg_v*arg_v + e6*arg_v**3 + e8*arg_v**4
1721 : tmp%scale_v(:) = EXP(-0.25_dp*dt*( &
1722 1192 : e_val(:) + trvg*infree))
1723 :
1724 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1725 : core_particle_set, shell_particle_set, nparticle_kind, &
1726 298 : shell_adiabatic, dt, u)
1727 :
1728 : CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
1729 0 : arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1730 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
1731 0 : arg_v(2) = arg_v(2)*simpar%dt_fact*simpar%dt_fact
1732 0 : arg_v(1) = arg_v(1)*simpar%dt_fact*simpar%dt_fact
1733 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
1734 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
1735 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
1736 0 : tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v)
1737 : tmp%scale_v(1) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
1738 0 : (1._dp + infree))
1739 0 : tmp%scale_v(2) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
1740 0 : tmp%scale_v(3) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
1741 : CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1742 : core_particle_set, shell_particle_set, nparticle_kind, &
1743 476 : shell_adiabatic, dt)
1744 :
1745 : END SELECT
1746 :
1747 476 : END SUBROUTINE rescaled_vv_first
1748 :
1749 0 : END MODULE integrator_utils
|
