(git:e7e05ae)
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 ! **************************************************************************************************
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  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 
1751 END MODULE integrator_utils
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