(git:374b731)
Loading...
Searching...
No Matches
integrator_utils.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
104CONTAINS
105
106! **************************************************************************************************
107!> \brief ...
108!> \param old ...
109!> \param particle_set ...
110!> \param npt ...
111!> \par History
112!> none
113!> \author CJM
114! **************************************************************************************************
115 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 natoms = SIZE(particle_set)
123 idim = SIZE(npt, 1)
124 jdim = SIZE(npt, 2)
125 cpassert(.NOT. ASSOCIATED(old))
126 ALLOCATE (old)
127
128 ALLOCATE (old%v(natoms, 3))
129 old%v = 0.0_dp
130 ALLOCATE (old%r(natoms, 3))
131 old%r = 0.0_dp
132 ALLOCATE (old%eps(idim, jdim))
133 old%eps = 0.0_dp
134 ALLOCATE (old%veps(idim, jdim))
135 old%veps = 0.0_dp
136 ALLOCATE (old%h(3, 3))
137 old%h = 0.0_dp
138
139 END SUBROUTINE allocate_old
140
141! **************************************************************************************************
142!> \brief ...
143!> \param old ...
144!> \par History
145!> none
146!> \author CJM
147! **************************************************************************************************
148 SUBROUTINE deallocate_old(old)
149 TYPE(old_variables_type), POINTER :: old
150
151 IF (ASSOCIATED(old)) THEN
152 IF (ASSOCIATED(old%v)) THEN
153 DEALLOCATE (old%v)
154 END IF
155 IF (ASSOCIATED(old%r)) THEN
156 DEALLOCATE (old%r)
157 END IF
158 IF (ASSOCIATED(old%eps)) THEN
159 DEALLOCATE (old%eps)
160 END IF
161 IF (ASSOCIATED(old%veps)) THEN
162 DEALLOCATE (old%veps)
163 END IF
164 IF (ASSOCIATED(old%h)) THEN
165 DEALLOCATE (old%h)
166 END IF
167 DEALLOCATE (old)
168 END IF
169
170 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 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 cpassert(.NOT. ASSOCIATED(tmp))
191 ALLOCATE (tmp)
192
193 ! Nullify Components
194 NULLIFY (tmp%pos)
195 NULLIFY (tmp%vel)
196 NULLIFY (tmp%shell_pos)
197 NULLIFY (tmp%shell_vel)
198 NULLIFY (tmp%core_pos)
199 NULLIFY (tmp%core_vel)
200 NULLIFY (tmp%itimes)
201
202 ! *** Allocate work storage for positions and velocities ***
203 ALLOCATE (tmp%pos(3, nparticle))
204
205 ALLOCATE (tmp%vel(3, nparticle))
206
207 tmp%pos(:, :) = 0.0_dp
208 tmp%vel(:, :) = 0.0_dp
209
210 IF (shell_adiabatic) THEN
211 ! *** Allocate work storage for positions and velocities ***
212 ALLOCATE (tmp%shell_pos(3, nshell))
213
214 ALLOCATE (tmp%core_pos(3, nshell))
215
216 ALLOCATE (tmp%shell_vel(3, nshell))
217
218 ALLOCATE (tmp%core_vel(3, nshell))
219
220 tmp%shell_pos(:, :) = 0.0_dp
221 tmp%shell_vel(:, :) = 0.0_dp
222 tmp%core_pos(:, :) = 0.0_dp
223 tmp%core_vel(:, :) = 0.0_dp
224 END IF
225
226 tmp%arg_r = 0.0_dp
227 tmp%arg_v = 0.0_dp
228 tmp%u = 0.0_dp
229 tmp%e_val = 0.0_dp
230 tmp%max_vel = 0.0_dp
231 tmp%max_vel_sc = 0.0_dp
232 tmp%max_dvel = 0.0_dp
233 tmp%max_dvel_sc = 0.0_dp
234 tmp%scale_r = 1.0_dp
235 tmp%scale_v = 1.0_dp
236 tmp%poly_r = 1.0_dp
237 tmp%poly_v = 1.0_dp
238
239 CALL get_md_env(md_env=md_env, itimes=tmp%itimes)
240
241 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 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 cpassert(ASSOCIATED(tmp))
272 my_pos = .false.
273 my_vel = .false.
274 my_deall = .true.
275 IF (PRESENT(pos)) my_pos = pos
276 IF (PRESENT(vel)) my_vel = vel
277 IF (PRESENT(should_deall_vel)) my_deall = should_deall_vel
278
279 ! *** Broadcast the new particle positions ***
280 IF (my_pos) THEN
281 CALL update_particle_set(particle_set, para_env, pos=tmp%pos)
282 DEALLOCATE (tmp%pos)
283
284 IF (shell_adiabatic) THEN
285 CALL update_particle_set(shell_particle_set, para_env, pos=tmp%shell_pos)
286 CALL update_particle_set(core_particle_set, para_env, pos=tmp%core_pos)
287 DEALLOCATE (tmp%shell_pos)
288 DEALLOCATE (tmp%core_pos)
289 END IF
290 END IF
291
292 ! *** Broadcast the new particle velocities ***
293 IF (my_vel) THEN
294 CALL update_particle_set(particle_set, para_env, vel=tmp%vel)
295 IF (shell_adiabatic) THEN
296 CALL update_particle_set(shell_particle_set, para_env, vel=tmp%shell_vel)
297 CALL update_particle_set(core_particle_set, para_env, vel=tmp%core_vel)
298 END IF
299 IF (my_deall) THEN
300 DEALLOCATE (tmp%vel)
301 IF (shell_adiabatic) THEN
302 DEALLOCATE (tmp%shell_vel)
303 DEALLOCATE (tmp%core_vel)
304 END IF
305 cpassert(.NOT. ASSOCIATED(tmp%pos))
306 cpassert(.NOT. ASSOCIATED(tmp%shell_pos))
307 cpassert(.NOT. ASSOCIATED(tmp%core_pos))
308 DEALLOCATE (tmp)
309 END IF
310 END IF
311
312 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 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 my_tmpv = .false.
346 IF (PRESENT(tmpv)) my_tmpv = tmpv
347
348 k = 0.0_dp
349 a = 0.0_dp
350 b = 0.0_dp
351 DO iparticle_kind = 1, nparticle_kind
352 atomic_kind => atomic_kind_set(iparticle_kind)
353 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
354 IF (mass /= 0.0_dp) THEN
355 nparticle_local = local_particles%n_el(iparticle_kind)
356 IF (my_tmpv) THEN
357 DO iparticle_local = 1, nparticle_local
358 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
359 k = k + 0.5_dp*mass*dot_product(tmp%vel(:, iparticle), tmp%vel(:, iparticle))
360 a = a + dot_product(tmp%vel(:, iparticle), particle_set(iparticle)%f(:))
361 b = b + (1.0_dp/mass)*dot_product(particle_set(iparticle)%f(:), particle_set(iparticle)%f(:))
362 END DO
363 ELSE
364 DO iparticle_local = 1, nparticle_local
365 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
366 k = k + 0.5_dp*mass*dot_product(particle_set(iparticle)%v(:), particle_set(iparticle)%v(:))
367 a = a + dot_product(particle_set(iparticle)%v(:), particle_set(iparticle)%f(:))
368 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 CALL para_env%sum(k)
374 CALL para_env%sum(a)
375 CALL para_env%sum(b)
376 a = a/(2.0_dp*k)
377 b = b/(2.0_dp*k)
378 rb = sqrt(b)
379 tmp%s = (a/b)*(cosh(dt*rb/2.0_dp) - 1) + sinh(dt*rb/2.0_dp)/rb
380 tmp%ds = (a/b)*(sinh(dt*rb/2.0_dp)*rb) + cosh(dt*rb/2.0_dp)
381
382 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 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 nparticle_kind = SIZE(atomic_kind_set)
411 SELECT CASE (char)
412 CASE ('F') ! forward assigning the old
413 DO iparticle_kind = 1, nparticle_kind
414 nparticle_local = local_particles%n_el(iparticle_kind)
415 DO iparticle_local = 1, nparticle_local
416 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
417 DO idim = 1, 3
418 old%v(iparticle, idim) = particle_set(iparticle)%v(idim)
419 old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
420 END DO
421 END DO
422 END DO
423 old%eps(:, :) = npt(:, :)%eps
424 old%veps(:, :) = npt(:, :)%v
425 old%h(:, :) = cell%hmat(:, :)
426 CASE ('B') ! back assigning the original variables
427 DO iparticle_kind = 1, nparticle_kind
428 nparticle_local = local_particles%n_el(iparticle_kind)
429 DO iparticle_local = 1, nparticle_local
430 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
431 DO idim = 1, 3
432 particle_set(iparticle)%v(idim) = old%v(iparticle, idim)
433 particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
434 END DO
435 END DO
436 END DO
437 npt(:, :)%eps = old%eps(:, :)
438 npt(:, :)%v = old%veps(:, :)
439 cell%hmat(:, :) = old%h(:, :)
440 END SELECT
441
442 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 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 nparticle_kind = SIZE(atomic_kind_set)
473 SELECT CASE (char)
474 CASE ('F') ! forward assigning the old
475 DO iparticle_kind = 1, nparticle_kind
476 nparticle_local = local_particles%n_el(iparticle_kind)
477 DO iparticle_local = 1, nparticle_local
478 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
479 DO idim = 1, 3
480 old%v(iparticle, idim) = vel(idim, iparticle)
481 old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
482 END DO
483 END DO
484 END DO
485 old%eps(:, :) = npt(:, :)%eps
486 old%veps(:, :) = npt(:, :)%v
487 old%h(:, :) = cell%hmat(:, :)
488 CASE ('B') ! back assigning the original variables
489 DO iparticle_kind = 1, nparticle_kind
490 nparticle_local = local_particles%n_el(iparticle_kind)
491 DO iparticle_local = 1, nparticle_local
492 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
493 DO idim = 1, 3
494 vel(idim, iparticle) = old%v(iparticle, idim)
495 particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
496 END DO
497 END DO
498 END DO
499 npt(:, :)%eps = old%eps(:, :)
500 npt(:, :)%v = old%veps(:, :)
501 cell%hmat(:, :) = old%h(:, :)
502 END SELECT
503
504 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 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 kin = 0.0_dp
540 DO ikind = 1, SIZE(molecule_kind_set)
541 ! Compute the total kinetic energy
542 nmol_local = local_molecules%n_el(ikind)
543 DO imol_local = 1, nmol_local
544 imol = local_molecules%list(ikind)%array(imol_local)
545 molecule => molecule_set(imol)
546 CALL get_molecule(molecule, first_atom=first_atom, &
547 last_atom=last_atom)
548 DO ipart = first_atom, last_atom
549 atomic_kind => particle_set(ipart)%atomic_kind
550 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
551 kin = kin + mass*particle_set(ipart)%v(1)* &
552 particle_set(ipart)%v(1)
553 kin = kin + mass*particle_set(ipart)%v(2)* &
554 particle_set(ipart)%v(2)
555 kin = kin + mass*particle_set(ipart)%v(3)* &
556 particle_set(ipart)%v(3)
557 END DO
558 END DO
559 END DO
560 !
561 CALL group%sum(kin)
562 !
563 ikin = 1.0_dp/kin
564 scale = 1.0_dp
565 alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
566 scale = scale*sqrt(1.0_dp + alpha*0.5_dp*dt)
567 ! Scale
568 DO ikind = 1, SIZE(molecule_kind_set)
569 nmol_local = local_molecules%n_el(ikind)
570 DO imol_local = 1, nmol_local
571 imol = local_molecules%list(ikind)%array(imol_local)
572 molecule => molecule_set(imol)
573 CALL get_molecule(molecule, first_atom=first_atom, &
574 last_atom=last_atom)
575 DO ipart = first_atom, last_atom
576 particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*scale
577 particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*scale
578 particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*scale
579 END DO
580 END DO
581 END DO
582
583 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 SUBROUTINE damp_v_velocity(molecule_kind_set, molecule_set, &
601 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 kin = 0.0_dp
621 ! Compute the total kinetic energy
622 DO ikind = 1, SIZE(molecule_kind_set)
623 nmol_local = local_molecules%n_el(ikind)
624 DO imol_local = 1, nmol_local
625 imol = local_molecules%list(ikind)%array(imol_local)
626 molecule => molecule_set(imol)
627 CALL get_molecule(molecule, first_atom=first_atom, &
628 last_atom=last_atom)
629 DO ipart = first_atom, last_atom
630 atomic_kind => particle_set(ipart)%atomic_kind
631 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
632 kin = kin + mass*vel(1, ipart)*vel(1, ipart)
633 kin = kin + mass*vel(2, ipart)*vel(2, ipart)
634 kin = kin + mass*vel(3, ipart)*vel(3, ipart)
635 END DO
636 END DO
637 END DO
638 !
639 CALL group%sum(kin)
640 !
641 ikin = 1.0_dp/kin
642 scale = 1.0_dp
643 alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
644 scale = scale*sqrt(1.0_dp + alpha*0.5_dp*dt)
645 ! Scale
646 DO ikind = 1, SIZE(molecule_kind_set)
647 nmol_local = local_molecules%n_el(ikind)
648 DO imol_local = 1, nmol_local
649 imol = local_molecules%list(ikind)%array(imol_local)
650 molecule => molecule_set(imol)
651 CALL get_molecule(molecule, first_atom=first_atom, &
652 last_atom=last_atom)
653 DO ipart = first_atom, last_atom
654 vel(1, ipart) = vel(1, ipart)*scale
655 vel(2, ipart) = vel(2, ipart)*scale
656 vel(3, ipart) = vel(3, ipart)*scale
657 END DO
658 END DO
659 END DO
660
661 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 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 scale = 1.0_dp
680 scale = scale*exp(-gamma1*0.25_dp*dt)
681 ! Scale
682 npt%v = npt%v*scale
683
684 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 SUBROUTINE rattle_roll_setup(old, gci, atomic_kind_set, particle_set, local_particles, &
714 molecule_kind_set, molecule_set, local_molecules, vel, dt, cell, npt, simpar, virial, &
715 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 IF (first) THEN
744 CALL update_pv(gci, simpar, atomic_kind_set, vel, particle_set, &
745 local_molecules, molecule_set, molecule_kind_set, &
746 local_particles, kin, pv_kin, virial, para_env)
747 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
748
749 END IF
750 first = .false.
751
752 ! assigning local variable
753 SELECT CASE (simpar%ensemble)
754 CASE (npt_i_ensemble, npt_ia_ensemble)
755 npt_loc(:, :)%v = 0.0_dp
756 npt_loc(:, :)%mass = 0.0_dp
757 npt_loc(1, 1)%v = npt(1, 1)%v
758 npt_loc(2, 2)%v = npt(1, 1)%v
759 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 npt_loc(3, 3)%mass = npt(1, 1)%mass
763 CASE (npt_f_ensemble)
764 npt_loc = npt
765 END SELECT
766
767 ! resetting
768 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 para_env, u, cell, local_particles)
772
773 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 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 pv_kin = 0.0_dp
817 kin = 0.0_dp
818 DO iparticle_kind = 1, SIZE(atomic_kind_set)
819 atomic_kind => atomic_kind_set(iparticle_kind)
820 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
821 nparticle_local = local_particles%n_el(iparticle_kind)
822 DO iparticle_local = 1, nparticle_local
823 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
824 DO i = 1, 3
825 DO j = 1, 3
826 pv_kin(i, j) = pv_kin(i, j) + &
827 mass*particle_set(iparticle)%v(i)* &
828 particle_set(iparticle)%v(j)
829 END DO
830 kin = kin + mass*particle_set(iparticle)%v(i)* &
831 particle_set(iparticle)%v(i)
832 END DO
833 END DO
834 END DO
835
836 CALL int_group%sum(pv_kin)
837 CALL int_group%sum(kin)
838
839 ! updating the constraint virial
840 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
841 molecule_set, &
842 molecule_kind_set, &
843 particle_set, virial, &
844 int_group)
845
846 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 SUBROUTINE update_pv_velocity(gci, simpar, atomic_kind_set, vel, particle_set, &
869 local_molecules, molecule_set, molecule_kind_set, &
870 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 pv_kin = 0.0_dp
892 kin = 0._dp
893 DO iparticle_kind = 1, SIZE(atomic_kind_set)
894 atomic_kind => atomic_kind_set(iparticle_kind)
895 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
896 nparticle_local = local_particles%n_el(iparticle_kind)
897 DO iparticle_local = 1, nparticle_local
898 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
899 DO i = 1, 3
900 DO j = 1, 3
901 pv_kin(i, j) = pv_kin(i, j) + &
902 mass*vel(i, iparticle)*vel(j, iparticle)
903 END DO
904 kin = kin + mass*vel(i, iparticle)*vel(i, iparticle)
905 END DO
906 END DO
907 END DO
908
909 CALL int_group%sum(pv_kin)
910 CALL int_group%sum(kin)
911
912 ! updating the constraint virial
913 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
914 molecule_set, &
915 molecule_kind_set, &
916 particle_set, virial, &
917 int_group)
918
919 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 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 cpassert(ASSOCIATED(box))
952
953 ! Initialize unit
954
955 unit = 0.0_dp
956 unit(1, 1) = 1.0_dp
957 unit(2, 2) = 1.0_dp
958 unit(3, 3) = 1.0_dp
959
960 IF (simpar%ensemble == npt_i_ensemble .OR. &
961 simpar%ensemble == npt_ia_ensemble .OR. &
962 simpar%ensemble == npe_i_ensemble) THEN
963 ! get force on barostat
964 fdotr = 0.0_dp
965 DO ii = 1, 3
966 fdotr = fdotr + virial%pv_virial(ii, ii) + &
967 virial%pv_constraint(ii, ii)
968 END DO
969
970 npt(:, :)%f = (1.0_dp + (3.0_dp*infree))*kin + fdotr - &
971 3.0_dp*simpar%p_ext*box%deth
972 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 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 ELSEIF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
984 simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
985 v = box%deth
986 vi = 1._dp/v
987 v0 = simpar%v0
988 v0i = 1._dp/v0
989 ! orthorhombic box ONLY
990 ! Chooses only the compressive solution
991 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 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 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 0.5_dp*simpar%dt*npt(:, :)%f/npt(:, :)%mass
1012
1013 ! Screen the dynamics of the barostat according user request
1014 IF (PRESENT(virial_components)) THEN
1015 SELECT CASE (virial_components)
1016 CASE (do_clv_xyz)
1017 ! Do nothing..
1018 CASE (do_clv_x)
1019 npt(1, 2)%v = 0.0_dp
1020 npt(1, 3)%v = 0.0_dp
1021 npt(1, 2)%f = 0.0_dp
1022 npt(1, 3)%f = 0.0_dp
1023 npt(2, 1)%v = 0.0_dp
1024 npt(2, 2)%v = 0.0_dp
1025 npt(2, 3)%v = 0.0_dp
1026 npt(2, 1)%f = 0.0_dp
1027 npt(2, 2)%f = 0.0_dp
1028 npt(2, 3)%f = 0.0_dp
1029 npt(3, 1)%v = 0.0_dp
1030 npt(3, 2)%v = 0.0_dp
1031 npt(3, 3)%v = 0.0_dp
1032 npt(3, 1)%f = 0.0_dp
1033 npt(3, 2)%f = 0.0_dp
1034 npt(3, 3)%f = 0.0_dp
1035 CASE (do_clv_y)
1036 npt(1, 1)%v = 0.0_dp
1037 npt(1, 2)%v = 0.0_dp
1038 npt(1, 3)%v = 0.0_dp
1039 npt(1, 1)%f = 0.0_dp
1040 npt(1, 2)%f = 0.0_dp
1041 npt(1, 3)%f = 0.0_dp
1042 npt(2, 1)%v = 0.0_dp
1043 npt(2, 3)%v = 0.0_dp
1044 npt(2, 1)%f = 0.0_dp
1045 npt(2, 3)%f = 0.0_dp
1046 npt(3, 1)%v = 0.0_dp
1047 npt(3, 2)%v = 0.0_dp
1048 npt(3, 3)%v = 0.0_dp
1049 npt(3, 1)%f = 0.0_dp
1050 npt(3, 2)%f = 0.0_dp
1051 npt(3, 3)%f = 0.0_dp
1052 CASE (do_clv_z)
1053 npt(1, 1)%v = 0.0_dp
1054 npt(1, 2)%v = 0.0_dp
1055 npt(1, 3)%v = 0.0_dp
1056 npt(1, 1)%f = 0.0_dp
1057 npt(1, 2)%f = 0.0_dp
1058 npt(1, 3)%f = 0.0_dp
1059 npt(2, 1)%v = 0.0_dp
1060 npt(2, 2)%v = 0.0_dp
1061 npt(2, 3)%v = 0.0_dp
1062 npt(2, 1)%f = 0.0_dp
1063 npt(2, 2)%f = 0.0_dp
1064 npt(2, 3)%f = 0.0_dp
1065 npt(3, 1)%v = 0.0_dp
1066 npt(3, 2)%v = 0.0_dp
1067 npt(3, 1)%f = 0.0_dp
1068 npt(3, 2)%f = 0.0_dp
1069 CASE (do_clv_xy)
1070 npt(1, 3)%v = 0.0_dp
1071 npt(1, 3)%f = 0.0_dp
1072 npt(2, 3)%v = 0.0_dp
1073 npt(2, 3)%f = 0.0_dp
1074 npt(3, 1)%v = 0.0_dp
1075 npt(3, 2)%v = 0.0_dp
1076 npt(3, 3)%v = 0.0_dp
1077 npt(3, 1)%f = 0.0_dp
1078 npt(3, 2)%f = 0.0_dp
1079 npt(3, 3)%f = 0.0_dp
1080 CASE (do_clv_xz)
1081 npt(1, 2)%v = 0.0_dp
1082 npt(1, 2)%f = 0.0_dp
1083 npt(2, 1)%v = 0.0_dp
1084 npt(2, 2)%v = 0.0_dp
1085 npt(2, 3)%v = 0.0_dp
1086 npt(2, 1)%f = 0.0_dp
1087 npt(2, 2)%f = 0.0_dp
1088 npt(2, 3)%f = 0.0_dp
1089 npt(3, 2)%v = 0.0_dp
1090 npt(3, 2)%f = 0.0_dp
1091 CASE (do_clv_yz)
1092 npt(1, 1)%v = 0.0_dp
1093 npt(1, 2)%v = 0.0_dp
1094 npt(1, 3)%v = 0.0_dp
1095 npt(1, 1)%f = 0.0_dp
1096 npt(1, 2)%f = 0.0_dp
1097 npt(1, 3)%f = 0.0_dp
1098 npt(2, 1)%v = 0.0_dp
1099 npt(2, 1)%f = 0.0_dp
1100 npt(3, 1)%v = 0.0_dp
1101 npt(3, 1)%f = 0.0_dp
1102 END SELECT
1103 END IF
1104
1105 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 SUBROUTINE vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1126 core_particle_set, shell_particle_set, nparticle_kind, &
1127 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 NULLIFY (atomic_kind, shell)
1151 tmp%max_vel = 0.0_dp
1152 tmp%max_vel_sc = 0.0_dp
1153 tmp%max_dvel = 0.0_dp
1154 tmp%max_dvel_sc = 0.0_dp
1155 tmp%max_dr = 0.0_dp
1156 tmp%max_dsc = 0.0_dp
1157 dsc = 0.0_dp
1158 dvsc = 0.0_dp
1159
1160 ! *** Velocity Verlet (first part) ***
1161 DO iparticle_kind = 1, nparticle_kind
1162 atomic_kind => atomic_kind_set(iparticle_kind)
1163 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell)
1164 IF (mass /= 0.0_dp) THEN
1165 dm = 0.5_dp*dt/mass
1166 IF (is_shell .AND. shell_adiabatic) THEN
1167 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1168 dms = 0.5_dp*dt/shell%mass_shell
1169 dmc = 0.5_dp*dt/shell%mass_core
1170 fac_masss = shell%mass_shell/mass
1171 fac_massc = shell%mass_core/mass
1172 nparticle_local = local_particles%n_el(iparticle_kind)
1173 IF (PRESENT(u)) THEN
1174 DO iparticle_local = 1, nparticle_local
1175 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1176 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 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 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 fac_massc*tmp%core_vel(:, shell_index)
1188 tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1189 fac_massc*tmp%core_pos(:, shell_index)
1190
1191 tmp%max_vel = max(tmp%max_vel, abs(tmp%vel(1, iparticle)), &
1192 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 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 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 abs(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1205
1206 dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1207 shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1208 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 shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1211 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 DO iparticle_local = 1, nparticle_local
1215 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1216 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 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 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 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 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 fac_massc*tmp%core_vel(:, shell_index)
1232 tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
1233 fac_massc*tmp%core_pos(:, shell_index)
1234 tmp%max_vel = max(tmp%max_vel, &
1235 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 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 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 abs(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1248 dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
1249 shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
1250 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 shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
1253 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 nparticle_local = local_particles%n_el(iparticle_kind)
1258 IF (PRESENT(u)) THEN
1259 DO iparticle_local = 1, nparticle_local
1260 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 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 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 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 abs(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
1274 END DO ! iparticle_local
1275 ELSE
1276 DO iparticle_local = 1, nparticle_local
1277 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 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 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 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 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 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 tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
1296
1297 ! overwrite positions of frozen particles
1298 IF (PRESENT(lfixd_list)) THEN
1299 IF (any(lfixd_list(:)%ifixd_index == iparticle)) THEN
1300 tmp%pos(1, iparticle) = particle_set(iparticle)%r(1)
1301 tmp%pos(2, iparticle) = particle_set(iparticle)%r(2)
1302 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 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 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 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 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 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 NULLIFY (atomic_kind, shell)
1363
1364 ! *** Velocity Verlet (second part) ***
1365 DO iparticle_kind = 1, nparticle_kind
1366 atomic_kind => atomic_kind_set(iparticle_kind)
1367 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
1368 shell_active=is_shell)
1369 IF (mass /= 0.0_dp) THEN
1370 dm = 0.5_dp*dt/mass
1371 IF (is_shell .AND. shell_adiabatic) THEN
1372 CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
1373 dms = 0.5_dp*dt/shell%mass_shell
1374 dmc = 0.5_dp*dt/shell%mass_core
1375 fac_masss = shell%mass_shell/mass
1376 fac_massc = shell%mass_core/mass
1377 nparticle_local = local_particles%n_el(iparticle_kind)
1378 IF (PRESENT(u)) THEN
1379 DO iparticle_local = 1, nparticle_local
1380 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1381 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 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 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 fac_massc*tmp%core_vel(1, shell_index)
1393 tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1394 fac_massc*tmp%core_vel(2, shell_index)
1395 tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1396 fac_massc*tmp%core_vel(3, shell_index)
1397 END DO ! iparticle_local
1398 ELSE
1399 DO iparticle_local = 1, nparticle_local
1400 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1401 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 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 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 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 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 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 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 fac_massc*tmp%core_vel(1, shell_index)
1423 tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
1424 fac_massc*tmp%core_vel(2, shell_index)
1425 tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
1426 fac_massc*tmp%core_vel(3, shell_index)
1427 END DO ! iparticle_local
1428 END IF
1429 ELSE
1430 nparticle_local = local_particles%n_el(iparticle_kind)
1431 IF (PRESENT(u)) THEN
1432 DO iparticle_local = 1, nparticle_local
1433 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1434 CALL transform_second(particle_set, tmp%vel, iparticle, &
1435 u, dm, tmp%poly_v, tmp%scale_v)
1436 END DO ! iparticle_local
1437 ELSE
1438 DO iparticle_local = 1, nparticle_local
1439 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 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 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 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 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 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 ur = matmul(transpose(u), particle_set(index)%r(:))
1487 uv = matmul(transpose(u), particle_set(index)%v(:))
1488 uf = matmul(transpose(u), particle_set(index)%f(:))
1489
1490 uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1491 uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1492 uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1493
1494 ur(1) = ur(1)*scale_r(1)*scale_r(1) + uv(1)*scale_r(1)*poly_r(1)*dt
1495 ur(2) = ur(2)*scale_r(2)*scale_r(2) + uv(2)*scale_r(2)*poly_r(2)*dt
1496 ur(3) = ur(3)*scale_r(3)*scale_r(3) + uv(3)*scale_r(3)*poly_r(3)*dt
1497
1498 pos(:, index) = matmul(u, ur)
1499 vel(:, index) = matmul(u, uv)
1500
1501 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 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 uv = matmul(transpose(u), vel(:, index))
1527 uf = matmul(transpose(u), particle_set(index)%f(:))
1528
1529 uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
1530 uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
1531 uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
1532
1533 vel(:, index) = matmul(u, uv)
1534
1535 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 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 dt_fact = 1.0_dp
1579 dt_fact_sc = 1.0_dp
1580 dt_fact_f = 1.0_dp
1581 dt_fact_v = 1.0_dp
1582 dt_old = dt
1583 dt_fact_old = simpar%dt_fact
1584 simpar%dt_fact = 1.0_dp
1585 NULLIFY (thermostats)
1586
1587 itime = itime + 1
1588 CALL para_env%max(tmp%max_dr)
1589 IF (tmp%max_dr > simpar%dr_tol) THEN
1590 CALL para_env%max(tmp%max_dvel)
1591 dt_fact = sqrt(simpar%dr_tol*dt/tmp%max_dvel)/dt
1592 END IF
1593
1594 IF (shell_adiabatic) THEN
1595 CALL para_env%max(tmp%max_dsc)
1596 IF (tmp%max_dsc > simpar%dsc_tol) THEN
1597 CALL para_env%max(tmp%max_dvel_sc)
1598 dt_fact_sc = sqrt(simpar%dsc_tol*dt/tmp%max_dvel_sc)/dt
1599 END IF
1600 END IF
1601
1602 dt_fact_f = min(dt_fact, dt_fact_sc, 1.0_dp)
1603 IF (dt_fact_f < 1.0_dp) THEN
1604 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 dt = dt*dt_fact_f
1607 simpar%dt_fact = dt_fact_f
1608 CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1609 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1610 END IF
1611
1612 dt_fact = 1.0_dp
1613 dt_fact_sc = 1.0_dp
1614 CALL para_env%max(tmp%max_dr)
1615 IF (tmp%max_dr > simpar%dr_tol) THEN
1616 CALL para_env%max(tmp%max_vel)
1617 dt_fact = simpar%dr_tol/tmp%max_vel/dt
1618 END IF
1619
1620 IF (shell_adiabatic) THEN
1621 CALL para_env%max(tmp%max_dsc)
1622 IF (tmp%max_dsc > simpar%dsc_tol) THEN
1623 CALL para_env%max(tmp%max_vel_sc)
1624 dt_fact_sc = simpar%dsc_tol/tmp%max_vel_sc/dt
1625 END IF
1626 END IF
1627 dt_fact_v = min(dt_fact, dt_fact_sc, 1.0_dp)
1628
1629 IF (dt_fact_v < 1.0_dp) THEN
1630 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 dt = dt*dt_fact_v
1633 simpar%dt_fact = dt_fact_f*dt_fact_v
1634 CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
1635 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
1636 END IF
1637
1638 simpar%dt_fact = dt_fact_f*dt_fact_v
1639 IF (simpar%dt_fact < 1.0_dp) THEN
1640 CALL get_md_env(md_env, t=time, thermostats=thermostats)
1641 time = time - dt_old + dt_old*simpar%dt_fact
1642 IF (ASSOCIATED(thermostats)) THEN
1643 CALL set_thermostats(thermostats, dt_fact=simpar%dt_fact)
1644 END IF
1645 END IF
1646
1647 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 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 infree = 1.0_dp/real(simpar%nfree, dp)
1688 arg_r = tmp%arg_r
1689 arg_v = tmp%arg_v
1690 u = tmp%u
1691 e_val = tmp%e_val
1692
1693 SELECT CASE (simpar%ensemble)
1694 CASE (nve_ensemble, nvt_ensemble)
1695 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
1696 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
1697 CASE (isokin_ensemble)
1698 tmp%scale_v(1:3) = sqrt(1.0_dp/tmp%ds)
1699 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 shell_adiabatic, dt)
1703
1704 CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
1705 arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1706 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 arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
1708 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 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 (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 shell_adiabatic, dt)
1715
1716 CASE (npt_f_ensemble, npe_f_ensemble)
1717 trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
1718 arg_r(:) = arg_r(:)*simpar%dt_fact*simpar%dt_fact
1719 tmp%poly_r = 1._dp + e2*arg_r + e4*arg_r*arg_r + e6*arg_r**3 + e8*arg_r**4
1720 tmp%scale_r(:) = exp(0.5_dp*dt*e_val(:))
1721 arg_v(:) = arg_v(:)*simpar%dt_fact*simpar%dt_fact
1722 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 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 shell_adiabatic, dt, u)
1729
1730 CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
1731 arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
1732 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 arg_v(2) = arg_v(2)*simpar%dt_fact*simpar%dt_fact
1734 arg_v(1) = arg_v(1)*simpar%dt_fact*simpar%dt_fact
1735 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 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 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 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 (1._dp + infree))
1741 tmp%scale_v(2) = exp(-0.25_dp*dt*npt(1, 1)%v*infree)
1742 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 shell_adiabatic, dt)
1746
1747 END SELECT
1748
1749 END SUBROUTINE rescaled_vv_first
1750
1751END MODULE integrator_utils
integrator_utils::damp_v
Definition integrator_utils.F:95
integrator_utils::set
Definition integrator_utils.F:87
integrator_utils::update_pv
Definition integrator_utils.F:91
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
barostat_types
Barostat structure: module containing barostat available for MD.
Definition barostat_types.F:12
barostat_types::do_clv_y
integer, parameter, public do_clv_y
Definition barostat_types.F:36
barostat_types::do_clv_xyz
integer, parameter, public do_clv_xyz
Definition barostat_types.F:36
barostat_types::do_clv_yz
integer, parameter, public do_clv_yz
Definition barostat_types.F:36
barostat_types::do_clv_xy
integer, parameter, public do_clv_xy
Definition barostat_types.F:36
barostat_types::do_clv_z
integer, parameter, public do_clv_z
Definition barostat_types.F:36
barostat_types::do_clv_xz
integer, parameter, public do_clv_xz
Definition barostat_types.F:36
barostat_types::do_clv_x
integer, parameter, public do_clv_x
Definition barostat_types.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
constraint_util
Contains routines useful for the application of constraints during MD.
Definition constraint_util.F:13
constraint_util::pv_constraint
subroutine, public pv_constraint(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, virial, group)
...
Definition constraint_util.F:198
constraint::rattle_roll_control
subroutine, public rattle_roll_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, vel, dt, simpar, vector, veps, roll_tol, iroll, para_env, u, cell, local_particles)
...
Definition constraint.F:493
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
extended_system_types
Lumps all possible extended system variables into one type for easy access and passing.
Definition extended_system_types.F:17
extended_system_types::debug_uniaxial_limit
logical, parameter, public debug_uniaxial_limit
Definition extended_system_types.F:39
extended_system_types::debug_isotropic_limit
logical, parameter, public debug_isotropic_limit
Definition extended_system_types.F:38
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::nph_uniaxial_ensemble
integer, parameter, public nph_uniaxial_ensemble
Definition input_constants.F:95
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:95
input_constants::isokin_ensemble
integer, parameter, public isokin_ensemble
Definition input_constants.F:95
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:95
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:95
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:95
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:95
input_constants::nve_ensemble
integer, parameter, public nve_ensemble
Definition input_constants.F:95
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:95
input_constants::nvt_ensemble
integer, parameter, public nvt_ensemble
Definition input_constants.F:95
integrator_utils
Provides integrator utility routines for the integrators.
Definition integrator_utils.F:15
integrator_utils::variable_timestep
subroutine, public variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, local_particles, particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
Compute the timestep rescaling factor.
Definition integrator_utils.F:1558
integrator_utils::rattle_roll_setup
subroutine, public rattle_roll_setup(old, gci, atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, vel, dt, cell, npt, simpar, virial, vector_v, roll_tol, iroll, infree, first, para_env, u)
update veps using multiplier obtained from SHAKE
Definition integrator_utils.F:716
integrator_utils::allocate_tmp
subroutine, public allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
allocate temporary variables to store positions and velocities used by the velocity-verlet integrator
Definition integrator_utils.F:185
integrator_utils::vv_second
subroutine, public vv_second(tmp, atomic_kind_set, local_particles, particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt, u)
Second half of the velocity-verlet algorithm : update velocity by half step using the new forces.
Definition integrator_utils.F:1343
integrator_utils::allocate_old
subroutine, public allocate_old(old, particle_set, npt)
...
Definition integrator_utils.F:116
integrator_utils::update_dealloc_tmp
subroutine, public update_dealloc_tmp(tmp, particle_set, shell_particle_set, core_particle_set, para_env, shell_adiabatic, pos, vel, should_deall_vel)
update positions and deallocate temporary variable
Definition integrator_utils.F:261
integrator_utils::damp_veps
elemental subroutine, public damp_veps(npt, gamma1, dt)
provides damping for barostat via nph_uniaxial_damped dynamics
Definition integrator_utils.F:673
integrator_utils::update_veps
subroutine, public update_veps(box, npt, simpar, pv_kin, kin, virial, infree, virial_components)
Routine to compute veps.
Definition integrator_utils.F:936
integrator_utils::get_s_ds
subroutine, public get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, dt, para_env, tmpv)
...
Definition integrator_utils.F:330
integrator_utils::vv_first
subroutine, public vv_first(tmp, atomic_kind_set, local_particles, particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt, u, lfixd_list)
First half of the velocity-verlet algorithm : update velocity by half step and positions by full step...
Definition integrator_utils.F:1128
integrator_utils::deallocate_old
subroutine, public deallocate_old(old)
...
Definition integrator_utils.F:149
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
md_environment_types
Definition md_environment_types.F:15
md_environment_types::get_md_env
subroutine, public get_md_env(md_env, itimes, constant, used_time, cell, simpar, npt, force_env, para_env, reftraj, t, init, first_time, fe_env, thermostats, barostat, thermostat_coeff, thermostat_part, thermostat_shell, thermostat_baro, thermostat_fast, thermostat_slow, md_ener, averages, thermal_regions, ehrenfest_md)
get components of MD environment type
Definition md_environment_types.F:195
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::update_particle_set
subroutine, public update_particle_set(particle_set, int_group, pos, vel, for, add)
...
Definition particle_types.F:106
shell_potential_types
Definition shell_potential_types.F:11
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
thermostat_types
Thermostat structure: module containing thermostat available for MD.
Definition thermostat_types.F:12
thermostat_types::set_thermostats
subroutine, public set_thermostats(thermostats, dt_fact)
access internal structures of thermostats
Definition thermostat_types.F:282
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
extended_system_types::npt_info_type
Definition extended_system_types.F:42
integrator_utils::old_variables_type
Definition integrator_utils.F:62
integrator_utils::tmp_variables_type
Definition integrator_utils.F:71
md_environment_types::md_environment_type
Definition md_environment_types.F:53
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_kind_types::local_fixd_constraint_type
Definition molecule_kind_types.F:156
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
thermostat_types::thermostats_type
Define thermostat types.
Definition thermostat_types.F:61
virial_types::virial_type
Definition virial_types.F:26