(git:1f285aa)
thermostat_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 Utilities for thermostats
10 !> \author teo [tlaino] - University of Zurich - 10.2007
11 ! **************************************************************************************************
12 MODULE thermostat_utils
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  get_atomic_kind
15  USE cell_types, ONLY: cell_type
16  USE cp_log_handling, ONLY: cp_get_default_logger,&
17  cp_logger_get_default_io_unit,&
18  cp_logger_type,&
19  cp_to_string
20  USE cp_output_handling, ONLY: cp_p_file,&
21  cp_print_key_finished_output,&
22  cp_print_key_should_output,&
23  cp_print_key_unit_nr
24  USE cp_units, ONLY: cp_unit_from_cp2k
25  USE distribution_1d_types, ONLY: distribution_1d_type
26  USE extended_system_types, ONLY: lnhc_parameters_type,&
27  map_info_type,&
28  npt_info_type
29  USE input_constants, ONLY: &
30  do_constr_atomic, do_constr_molec, do_region_defined, do_region_global, do_region_massive, &
31  do_region_molecule, do_thermo_al, do_thermo_communication, do_thermo_csvr, do_thermo_gle, &
32  do_thermo_no_communication, do_thermo_nose, isokin_ensemble, langevin_ensemble, &
33  npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
34  npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, &
35  nvt_ensemble, reftraj_ensemble
36  USE input_section_types, ONLY: section_vals_get,&
37  section_vals_get_subs_vals,&
38  section_vals_type,&
39  section_vals_val_get
40  USE kinds, ONLY: default_string_length,&
41  dp
42  USE machine, ONLY: m_flush
43  USE message_passing, ONLY: mp_comm_type,&
44  mp_para_env_type
45  USE molecule_kind_types, ONLY: get_molecule_kind,&
46  get_molecule_kind_set,&
47  molecule_kind_type
48  USE molecule_list_types, ONLY: molecule_list_type
49  USE molecule_types, ONLY: get_molecule,&
50  global_constraint_type,&
51  molecule_type
52  USE motion_utils, ONLY: rot_ana
53  USE particle_list_types, ONLY: particle_list_type
54  USE particle_types, ONLY: particle_type
55  USE physcon, ONLY: femtoseconds
56  USE qmmm_types, ONLY: qmmm_env_type
57  USE shell_potential_types, ONLY: shell_kind_type
58  USE simpar_types, ONLY: simpar_type
59  USE thermostat_types, ONLY: thermostat_info_type,&
60  thermostat_type,&
61  thermostats_type
62 #include "../../base/base_uses.f90"
63 
64  IMPLICIT NONE
65 
66  PRIVATE
67  PUBLIC :: compute_degrees_of_freedom, &
68  compute_nfree, &
69  setup_thermostat_info, &
70  setup_adiabatic_thermostat_info, &
71  ke_region_baro, &
72  ke_region_particles, &
73  ke_region_shells, &
74  vel_rescale_baro, &
75  vel_rescale_particles, &
76  vel_rescale_shells, &
77  get_thermostat_energies, &
78  get_nhc_energies, &
79  get_kin_energies, &
80  communication_thermo_low2, &
81  print_thermostats_status, &
82  momentum_region_particles
83 
84  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermostat_utils'
85 
86 CONTAINS
87 
88 ! **************************************************************************************************
89 !> \brief ...
90 !> \param cell ...
91 !> \param simpar ...
92 !> \param molecule_kind_set ...
93 !> \param print_section ...
94 !> \param particles ...
95 !> \param gci ...
96 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
97 ! **************************************************************************************************
98  SUBROUTINE compute_nfree(cell, simpar, molecule_kind_set, &
99  print_section, particles, gci)
100 
101  TYPE(cell_type), POINTER :: cell
102  TYPE(simpar_type), POINTER :: simpar
103  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
104  TYPE(section_vals_type), POINTER :: print_section
105  TYPE(particle_list_type), POINTER :: particles
106  TYPE(global_constraint_type), POINTER :: gci
107 
108  INTEGER :: natom, nconstraint_ext, nconstraint_int, &
109  nrestraints_int, rot_dof, &
110  roto_trasl_dof
111 
112 ! Retrieve information on number of atoms, constraints (external and internal)
113 
114  CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
115  natom=natom, nconstraint=nconstraint_int, nrestraints=nrestraints_int)
116 
117  ! Compute degrees of freedom
118  CALL rot_ana(particles%els, dof=roto_trasl_dof, rot_dof=rot_dof, &
119  print_section=print_section, keep_rotations=.false., &
120  mass_weighted=.true., natoms=natom)
121 
122  roto_trasl_dof = roto_trasl_dof - min(sum(cell%perd(1:3)), rot_dof)
123 
124  ! Saving this value of simpar preliminar to the real count of constraints..
125  simpar%nfree_rot_transl = roto_trasl_dof
126 
127  ! compute the total number of degrees of freedom for temperature
128  nconstraint_ext = gci%ntot - gci%nrestraint
129  simpar%nfree = 3*natom - nconstraint_int - nconstraint_ext - roto_trasl_dof
130 
131  END SUBROUTINE compute_nfree
132 
133 ! **************************************************************************************************
134 !> \brief ...
135 !> \param thermostats ...
136 !> \param cell ...
137 !> \param simpar ...
138 !> \param molecule_kind_set ...
139 !> \param local_molecules ...
140 !> \param molecules ...
141 !> \param particles ...
142 !> \param print_section ...
143 !> \param region_sections ...
144 !> \param gci ...
145 !> \param region ...
146 !> \param qmmm_env ...
147 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
148 ! **************************************************************************************************
149  SUBROUTINE compute_degrees_of_freedom(thermostats, cell, simpar, molecule_kind_set, &
150  local_molecules, molecules, particles, print_section, region_sections, gci, &
151  region, qmmm_env)
152 
153  TYPE(thermostats_type), POINTER :: thermostats
154  TYPE(cell_type), POINTER :: cell
155  TYPE(simpar_type), POINTER :: simpar
156  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
157  TYPE(distribution_1d_type), POINTER :: local_molecules
158  TYPE(molecule_list_type), POINTER :: molecules
159  TYPE(particle_list_type), POINTER :: particles
160  TYPE(section_vals_type), POINTER :: print_section, region_sections
161  TYPE(global_constraint_type), POINTER :: gci
162  INTEGER, INTENT(IN) :: region
163  TYPE(qmmm_env_type), POINTER :: qmmm_env
164 
165  INTEGER :: iw, natom, nconstraint_ext, &
166  nconstraint_int, nrestraints_int, &
167  rot_dof, roto_trasl_dof
168  TYPE(cp_logger_type), POINTER :: logger
169 
170 ! Retrieve information on number of atoms, constraints (external and internal)
171 
172  CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
173  natom=natom, nconstraint=nconstraint_int, nrestraints=nrestraints_int)
174 
175  ! Compute degrees of freedom
176  CALL rot_ana(particles%els, dof=roto_trasl_dof, rot_dof=rot_dof, &
177  print_section=print_section, keep_rotations=.false., &
178  mass_weighted=.true., natoms=natom)
179 
180  roto_trasl_dof = roto_trasl_dof - min(sum(cell%perd(1:3)), rot_dof)
181 
182  ! Collect info about thermostats
183  CALL setup_thermostat_info(thermostats%thermostat_info_part, molecule_kind_set, &
184  local_molecules, molecules, particles, region, simpar%ensemble, roto_trasl_dof, &
185  region_sections=region_sections, qmmm_env=qmmm_env)
186 
187  ! Saving this value of simpar preliminar to the real count of constraints..
188  simpar%nfree_rot_transl = roto_trasl_dof
189 
190  ! compute the total number of degrees of freedom for temperature
191  nconstraint_ext = gci%ntot - gci%nrestraint
192  simpar%nfree = 3*natom - nconstraint_int - nconstraint_ext - roto_trasl_dof
193 
194  logger => cp_get_default_logger()
195  iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
196  extension=".log")
197  IF (iw > 0) THEN
198  WRITE (iw, '(/,T2,A)') &
199  'DOF| Calculation of degrees of freedom'
200  WRITE (iw, '(T2,A,T71,I10)') &
201  'DOF| Number of atoms', natom, &
202  'DOF| Number of intramolecular constraints', nconstraint_int, &
203  'DOF| Number of intermolecular constraints', nconstraint_ext, &
204  'DOF| Invariants (translations + rotations)', roto_trasl_dof, &
205  'DOF| Degrees of freedom', simpar%nfree
206  WRITE (iw, '(/,T2,A)') &
207  'DOF| Restraints information'
208  WRITE (iw, '(T2,A,T71,I10)') &
209  'DOF| Number of intramolecular restraints', nrestraints_int, &
210  'DOF| Number of intermolecular restraints', gci%nrestraint
211  END IF
212  CALL cp_print_key_finished_output(iw, logger, print_section, &
213  "PROGRAM_RUN_INFO")
214 
215  END SUBROUTINE compute_degrees_of_freedom
216 
217 ! **************************************************************************************************
218 !> \brief ...
219 !> \param thermostat_info ...
220 !> \param molecule_kind_set ...
221 !> \param local_molecules ...
222 !> \param molecules ...
223 !> \param particles ...
224 !> \param region ...
225 !> \param ensemble ...
226 !> \param nfree ...
227 !> \param shell ...
228 !> \param region_sections ...
229 !> \param qmmm_env ...
230 !> \author 10.2011 CJM - PNNL
231 ! **************************************************************************************************
232  SUBROUTINE setup_adiabatic_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, &
233  molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
234  TYPE(thermostat_info_type), POINTER :: thermostat_info
235  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
236  TYPE(distribution_1d_type), POINTER :: local_molecules
237  TYPE(molecule_list_type), POINTER :: molecules
238  TYPE(particle_list_type), POINTER :: particles
239  INTEGER, INTENT(IN) :: region, ensemble
240  INTEGER, INTENT(INOUT), OPTIONAL :: nfree
241  LOGICAL, INTENT(IN), OPTIONAL :: shell
242  TYPE(section_vals_type), POINTER :: region_sections
243  TYPE(qmmm_env_type), POINTER :: qmmm_env
244 
245  INTEGER :: dis_type, first_atom, i, ikind, imol, imol_global, ipart, itherm, katom, &
246  last_atom, natom, natom_local, nkind, nmol_local, nmol_per_kind, nmolecule, nshell, &
247  number, stat, sum_of_thermostats
248  INTEGER, POINTER :: molecule_list(:), thermolist(:)
249  LOGICAL :: check, do_shell, nointer, on_therm
250  TYPE(molecule_kind_type), POINTER :: molecule_kind
251  TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
252 
253  NULLIFY (molecule_kind, molecule, thermostat_info%map_loc_thermo_gen, thermolist)
254  nkind = SIZE(molecule_kind_set)
255  do_shell = .false.
256  IF (PRESENT(shell)) do_shell = shell
257  ! Counting the global number of thermostats
258  sum_of_thermostats = 0
259  ! Variable to denote independent thermostats (no communication necessary)
260  nointer = .true.
261  check = .true.
262  number = 0
263  dis_type = do_thermo_no_communication
264 
265  CALL get_adiabatic_region_info(region_sections, sum_of_thermostats, &
266  thermolist=thermolist, &
267  molecule_kind_set=molecule_kind_set, &
268  molecules=molecules, particles=particles, qmmm_env=qmmm_env)
269 
270 ! map_loc_thermo_gen=>thermostat_info%map_loc_thermo_gen
271  molecule_set => molecules%els
272  SELECT CASE (ensemble)
273  CASE DEFAULT
274  cpabort('Unknown ensemble')
275  CASE (nvt_adiabatic_ensemble)
276  SELECT CASE (region)
277  CASE (do_region_global)
278  ! Global Thermostat
279  nointer = .false.
280  sum_of_thermostats = 1
281  CASE (do_region_molecule)
282  ! Molecular Thermostat
283  itherm = 0
284  DO ikind = 1, nkind
285  molecule_kind => molecule_kind_set(ikind)
286  nmol_per_kind = local_molecules%n_el(ikind)
287  CALL get_molecule_kind(molecule_kind, natom=natom, &
288  molecule_list=molecule_list)
289 ! use thermolist ( ipart ) to get global indexing correct
290  DO imol_global = 1, SIZE(molecule_list)
291  molecule => molecule_set(molecule_list(imol_global))
292  CALL get_molecule(molecule, first_atom=first_atom, &
293  last_atom=last_atom)
294  on_therm = .true.
295  DO katom = first_atom, last_atom
296  IF (thermolist(katom) == huge(0)) THEN
297  on_therm = .false.
298  EXIT
299  END IF
300  END DO
301  IF (on_therm) THEN
302  itherm = itherm + 1
303  DO katom = first_atom, last_atom
304  thermolist(katom) = itherm
305  END DO
306  END IF
307  END DO
308  END DO
309  DO i = 1, nkind
310  molecule_kind => molecule_kind_set(i)
311  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, nshell=nshell)
312  IF ((do_shell) .AND. (nshell == 0)) nmolecule = 0
313  sum_of_thermostats = sum_of_thermostats + nmolecule
314  END DO
315  ! If we have ONE kind and ONE molecule, then effectively we have a GLOBAL thermostat
316  ! and the degrees of freedom will be computed correctly for this special case
317  IF ((nmolecule == 1) .AND. (nkind == 1)) nointer = .false.
318  CASE (do_region_massive)
319  ! Massive Thermostat
320  DO i = 1, nkind
321  molecule_kind => molecule_kind_set(i)
322  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
323  natom=natom, nshell=nshell)
324  IF (do_shell) natom = nshell
325  sum_of_thermostats = sum_of_thermostats + 3*natom*nmolecule
326  END DO
327  END SELECT
328 
329  natom_local = 0
330  DO ikind = 1, SIZE(molecule_kind_set)
331  nmol_per_kind = local_molecules%n_el(ikind)
332  DO imol = 1, nmol_per_kind
333  i = local_molecules%list(ikind)%array(imol)
334  molecule => molecule_set(i)
335  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
336  DO ipart = first_atom, last_atom
337  natom_local = natom_local + 1
338  END DO
339  END DO
340  END DO
341 
342  ! Now map the local atoms with the corresponding thermostat
343  ALLOCATE (thermostat_info%map_loc_thermo_gen(natom_local), stat=stat)
344  thermostat_info%map_loc_thermo_gen = huge(0)
345  cpassert(stat == 0)
346  natom_local = 0
347  DO ikind = 1, SIZE(molecule_kind_set)
348  nmol_per_kind = local_molecules%n_el(ikind)
349  DO imol = 1, nmol_per_kind
350  i = local_molecules%list(ikind)%array(imol)
351  molecule => molecule_set(i)
352  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
353  DO ipart = first_atom, last_atom
354  natom_local = natom_local + 1
355 ! only map the correct region to the thermostat
356  IF (thermolist(ipart) /= huge(0)) &
357  thermostat_info%map_loc_thermo_gen(natom_local) = thermolist(ipart)
358  END DO
359  END DO
360  END DO
361  ! Here we decide which parallel algorithm to use.
362  ! if there are only massive and molecule type thermostats we can use
363  ! a local scheme, in cases involving any combination with a
364  ! global thermostat we assume a coupling of degrees of freedom
365  ! from different processors
366  IF (nointer) THEN
367  ! Distributed thermostats, no interaction
368  dis_type = do_thermo_no_communication
369  ! we only count thermostats on this processor
370  number = 0
371  DO ikind = 1, nkind
372  nmol_local = local_molecules%n_el(ikind)
373  molecule_kind => molecule_kind_set(ikind)
374  CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
375  IF (do_shell) THEN
376  natom = nshell
377  IF (nshell == 0) nmol_local = 0
378  END IF
379  IF (region == do_region_molecule) THEN
380  number = number + nmol_local
381  ELSE IF (region == do_region_massive) THEN
382  number = number + 3*nmol_local*natom
383  ELSE
384  cpabort('Invalid region setup')
385  END IF
386  END DO
387  ELSE
388  ! REPlicated thermostats, INTERacting via communication
389  dis_type = do_thermo_communication
390  IF ((region == do_region_global) .OR. (region == do_region_molecule)) number = 1
391  END IF
392 
393  IF (PRESENT(nfree)) THEN
394  ! re-initializing simpar%nfree to zero because of multiple thermostats in the adiabatic sampling
395  nfree = 0
396  END IF
397  END SELECT
398 
399  ! Saving information about thermostats
400  thermostat_info%sum_of_thermostats = sum_of_thermostats
401  thermostat_info%number_of_thermostats = number
402  thermostat_info%dis_type = dis_type
403 
404  DEALLOCATE (thermolist)
405 
406  END SUBROUTINE setup_adiabatic_thermostat_info
407 
408 ! **************************************************************************************************
409 !> \brief ...
410 !> \param region_sections ...
411 !> \param sum_of_thermostats ...
412 !> \param thermolist ...
413 !> \param molecule_kind_set ...
414 !> \param molecules ...
415 !> \param particles ...
416 !> \param qmmm_env ...
417 !> \author 10.2011 CJM -PNNL
418 ! **************************************************************************************************
419  SUBROUTINE get_adiabatic_region_info(region_sections, sum_of_thermostats, &
420  thermolist, molecule_kind_set, molecules, particles, &
421  qmmm_env)
422  TYPE(section_vals_type), POINTER :: region_sections
423  INTEGER, INTENT(INOUT), OPTIONAL :: sum_of_thermostats
424  INTEGER, DIMENSION(:), POINTER :: thermolist(:)
425  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
426  TYPE(molecule_list_type), POINTER :: molecules
427  TYPE(particle_list_type), POINTER :: particles
428  TYPE(qmmm_env_type), POINTER :: qmmm_env
429 
430  CHARACTER(LEN=default_string_length), &
431  DIMENSION(:), POINTER :: tmpstringlist
432  INTEGER :: first_atom, i, ig, ikind, ilist, imol, &
433  ipart, itherm, jg, last_atom, &
434  mregions, n_rep, nregions, output_unit
435  INTEGER, DIMENSION(:), POINTER :: tmplist
436  TYPE(cp_logger_type), POINTER :: logger
437  TYPE(molecule_kind_type), POINTER :: molecule_kind
438  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
439  TYPE(molecule_type), POINTER :: molecule
440 
441  NULLIFY (tmplist, tmpstringlist, thermolist, molecule_kind, molecule, molecule_set)
442  NULLIFY (logger)
443  logger => cp_get_default_logger()
444  output_unit = cp_logger_get_default_io_unit(logger)
445  ! CPASSERT(.NOT.(ASSOCIATED(map_loc_thermo_gen)))
446  CALL section_vals_get(region_sections, n_repetition=nregions)
447  ALLOCATE (thermolist(particles%n_els))
448  thermolist = huge(0)
449  molecule_set => molecules%els
450  mregions = nregions
451  itherm = 0
452  DO ig = 1, mregions
453  CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, n_rep_val=n_rep)
454  DO jg = 1, n_rep
455  CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, i_rep_val=jg, i_vals=tmplist)
456  DO i = 1, SIZE(tmplist)
457  ipart = tmplist(i)
458  cpassert(((ipart > 0) .AND. (ipart <= particles%n_els)))
459  IF (thermolist(ipart) == huge(0)) THEN
460  itherm = itherm + 1
461  thermolist(ipart) = itherm
462  ELSE
463  cpabort("")
464  END IF
465  END DO
466  END DO
467  CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, n_rep_val=n_rep)
468  DO jg = 1, n_rep
469  CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, i_rep_val=jg, c_vals=tmpstringlist)
470  DO ilist = 1, SIZE(tmpstringlist)
471  DO ikind = 1, SIZE(molecule_kind_set)
472  molecule_kind => molecule_kind_set(ikind)
473  IF (molecule_kind%name == tmpstringlist(ilist)) THEN
474  DO imol = 1, SIZE(molecule_kind%molecule_list)
475  molecule => molecule_set(molecule_kind%molecule_list(imol))
476  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
477  DO ipart = first_atom, last_atom
478  IF (thermolist(ipart) == huge(0)) THEN
479  itherm = itherm + 1
480  thermolist(ipart) = itherm
481  ELSE
482  cpabort("")
483  END IF
484  END DO
485  END DO
486  END IF
487  END DO
488  END DO
489  END DO
490  CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
491  subsys_qm=.false., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
492  CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
493  subsys_qm=.true., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
494  END DO
495 
496  cpassert(.NOT. all(thermolist == huge(0)))
497 
498 ! natom_local = 0
499 ! DO ikind = 1, SIZE(molecule_kind_set)
500 ! nmol_per_kind = local_molecules%n_el(ikind)
501 ! DO imol = 1, nmol_per_kind
502 ! i = local_molecules%list(ikind)%array(imol)
503 ! molecule => molecule_set(i)
504 ! CALL get_molecule ( molecule, first_atom = first_atom, last_atom = last_atom )
505 ! DO ipart = first_atom, last_atom
506 ! natom_local = natom_local + 1
507 ! END DO
508 ! END DO
509 ! END DO
510 
511  ! Now map the local atoms with the corresponding thermostat
512 ! ALLOCATE(map_loc_thermo_gen(natom_local),stat=stat)
513 ! map_loc_thermo_gen = HUGE ( 0 )
514 ! CPPostcondition(stat==0,cp_failure_level,routineP,failure)
515 ! natom_local = 0
516 ! DO ikind = 1, SIZE(molecule_kind_set)
517 ! nmol_per_kind = local_molecules%n_el(ikind)
518 ! DO imol = 1, nmol_per_kind
519 ! i = local_molecules%list(ikind)%array(imol)
520 ! molecule => molecule_set(i)
521 ! CALL get_molecule ( molecule, first_atom = first_atom, last_atom = last_atom )
522 ! DO ipart = first_atom, last_atom
523 ! natom_local = natom_local + 1
524 ! only map the correct region to the thermostat
525 ! IF ( thermolist (ipart ) /= HUGE ( 0 ) ) &
526 ! map_loc_thermo_gen(natom_local) = thermolist(ipart)
527 ! END DO
528 ! END DO
529 ! END DO
530 
531 ! DEALLOCATE(thermolist, stat=stat)
532 ! CPPostcondition(stat==0,cp_failure_level,routineP,failure)
533  END SUBROUTINE get_adiabatic_region_info
534 ! **************************************************************************************************
535 !> \brief ...
536 !> \param thermostat_info ...
537 !> \param molecule_kind_set ...
538 !> \param local_molecules ...
539 !> \param molecules ...
540 !> \param particles ...
541 !> \param region ...
542 !> \param ensemble ...
543 !> \param nfree ...
544 !> \param shell ...
545 !> \param region_sections ...
546 !> \param qmmm_env ...
547 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
548 ! **************************************************************************************************
549  SUBROUTINE setup_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, &
550  molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
551  TYPE(thermostat_info_type), POINTER :: thermostat_info
552  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
553  TYPE(distribution_1d_type), POINTER :: local_molecules
554  TYPE(molecule_list_type), POINTER :: molecules
555  TYPE(particle_list_type), POINTER :: particles
556  INTEGER, INTENT(IN) :: region, ensemble
557  INTEGER, INTENT(INOUT), OPTIONAL :: nfree
558  LOGICAL, INTENT(IN), OPTIONAL :: shell
559  TYPE(section_vals_type), POINTER :: region_sections
560  TYPE(qmmm_env_type), POINTER :: qmmm_env
561 
562  INTEGER :: dis_type, i, ikind, natom, nkind, &
563  nmol_local, nmolecule, nshell, number, &
564  sum_of_thermostats
565  LOGICAL :: check, do_shell, nointer
566  TYPE(molecule_kind_type), POINTER :: molecule_kind
567 
568  NULLIFY (molecule_kind)
569  nkind = SIZE(molecule_kind_set)
570  do_shell = .false.
571  IF (PRESENT(shell)) do_shell = shell
572  ! Counting the global number of thermostats
573  sum_of_thermostats = 0
574  ! Variable to denote independent thermostats (no communication necessary)
575  nointer = .true.
576  check = .true.
577  number = 0
578  dis_type = do_thermo_no_communication
579 
580  SELECT CASE (ensemble)
581  CASE DEFAULT
582  cpabort('Unknown ensemble')
583  CASE (isokin_ensemble, nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble, &
584  reftraj_ensemble, langevin_ensemble)
585  ! Do Nothing
586  CASE (nve_ensemble, nvt_ensemble, nvt_adiabatic_ensemble, npt_i_ensemble, &
587  npt_f_ensemble, npe_i_ensemble, npe_f_ensemble, npt_ia_ensemble)
588  IF (ensemble == nve_ensemble) check = do_shell
589  IF (check) THEN
590  SELECT CASE (region)
591  CASE (do_region_global)
592  ! Global Thermostat
593  nointer = .false.
594  sum_of_thermostats = 1
595  CASE (do_region_molecule)
596  ! Molecular Thermostat
597  DO i = 1, nkind
598  molecule_kind => molecule_kind_set(i)
599  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, nshell=nshell)
600  IF ((do_shell) .AND. (nshell == 0)) nmolecule = 0
601  sum_of_thermostats = sum_of_thermostats + nmolecule
602  END DO
603  ! If we have ONE kind and ONE molecule, then effectively we have a GLOBAL thermostat
604  ! and the degrees of freedom will be computed correctly for this special case
605  IF ((nmolecule == 1) .AND. (nkind == 1)) nointer = .false.
606  CASE (do_region_massive)
607  ! Massive Thermostat
608  DO i = 1, nkind
609  molecule_kind => molecule_kind_set(i)
610  CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
611  natom=natom, nshell=nshell)
612  IF (do_shell) natom = nshell
613  sum_of_thermostats = sum_of_thermostats + 3*natom*nmolecule
614  END DO
615  CASE (do_region_defined)
616  ! User defined region to thermostat..
617  nointer = .false.
618  ! Determine the number of thermostats defined in the input
619  CALL section_vals_get(region_sections, n_repetition=sum_of_thermostats)
620  IF (sum_of_thermostats < 1) &
621  CALL cp_abort(__location__, &
622  "Provide at least 1 region (&DEFINE_REGION) when using the thermostat type DEFINED")
623  END SELECT
624 
625  ! Here we decide which parallel algorithm to use.
626  ! if there are only massive and molecule type thermostats we can use
627  ! a local scheme, in cases involving any combination with a
628  ! global thermostat we assume a coupling of degrees of freedom
629  ! from different processors
630  IF (nointer) THEN
631  ! Distributed thermostats, no interaction
632  dis_type = do_thermo_no_communication
633  ! we only count thermostats on this processor
634  number = 0
635  DO ikind = 1, nkind
636  nmol_local = local_molecules%n_el(ikind)
637  molecule_kind => molecule_kind_set(ikind)
638  CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
639  IF (do_shell) THEN
640  natom = nshell
641  IF (nshell == 0) nmol_local = 0
642  END IF
643  IF (region == do_region_molecule) THEN
644  number = number + nmol_local
645  ELSE IF (region == do_region_massive) THEN
646  number = number + 3*nmol_local*natom
647  ELSE
648  cpabort('Invalid region setup')
649  END IF
650  END DO
651  ELSE
652  ! REPlicated thermostats, INTERacting via communication
653  dis_type = do_thermo_communication
654  IF ((region == do_region_global) .OR. (region == do_region_molecule)) THEN
655  number = 1
656  ELSE IF (region == do_region_defined) THEN
657  CALL get_defined_region_info(region_sections, number, sum_of_thermostats, &
658  map_loc_thermo_gen=thermostat_info%map_loc_thermo_gen, &
659  local_molecules=local_molecules, molecule_kind_set=molecule_kind_set, &
660  molecules=molecules, particles=particles, qmmm_env=qmmm_env)
661  END IF
662  END IF
663 
664  IF (PRESENT(nfree)) THEN
665  IF ((sum_of_thermostats > 1) .OR. (dis_type == do_thermo_no_communication)) THEN
666  ! re-initializing simpar%nfree to zero because of multiple thermostats
667  nfree = 0
668  END IF
669  END IF
670  END IF
671  END SELECT
672 
673  ! Saving information about thermostats
674  thermostat_info%sum_of_thermostats = sum_of_thermostats
675  thermostat_info%number_of_thermostats = number
676  thermostat_info%dis_type = dis_type
677  END SUBROUTINE setup_thermostat_info
678 
679 ! **************************************************************************************************
680 !> \brief ...
681 !> \param region_sections ...
682 !> \param number ...
683 !> \param sum_of_thermostats ...
684 !> \param map_loc_thermo_gen ...
685 !> \param local_molecules ...
686 !> \param molecule_kind_set ...
687 !> \param molecules ...
688 !> \param particles ...
689 !> \param qmmm_env ...
690 !> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
691 ! **************************************************************************************************
692  SUBROUTINE get_defined_region_info(region_sections, number, sum_of_thermostats, &
693  map_loc_thermo_gen, local_molecules, molecule_kind_set, molecules, particles, &
694  qmmm_env)
695  TYPE(section_vals_type), POINTER :: region_sections
696  INTEGER, INTENT(OUT), OPTIONAL :: number
697  INTEGER, INTENT(INOUT), OPTIONAL :: sum_of_thermostats
698  INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
699  TYPE(distribution_1d_type), POINTER :: local_molecules
700  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
701  TYPE(molecule_list_type), POINTER :: molecules
702  TYPE(particle_list_type), POINTER :: particles
703  TYPE(qmmm_env_type), POINTER :: qmmm_env
704 
705  CHARACTER(LEN=default_string_length), &
706  DIMENSION(:), POINTER :: tmpstringlist
707  INTEGER :: first_atom, i, ig, ikind, ilist, imol, ipart, jg, last_atom, mregions, n_rep, &
708  natom_local, nmol_per_kind, nregions, output_unit
709  INTEGER, DIMENSION(:), POINTER :: thermolist, tmp, tmplist
710  TYPE(cp_logger_type), POINTER :: logger
711  TYPE(molecule_kind_type), POINTER :: molecule_kind
712  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
713  TYPE(molecule_type), POINTER :: molecule
714 
715  NULLIFY (tmplist, tmpstringlist, thermolist, molecule_kind, molecule, molecule_set)
716  NULLIFY (logger)
717  logger => cp_get_default_logger()
718  output_unit = cp_logger_get_default_io_unit(logger)
719  cpassert(.NOT. (ASSOCIATED(map_loc_thermo_gen)))
720  CALL section_vals_get(region_sections, n_repetition=nregions)
721  ALLOCATE (thermolist(particles%n_els))
722  thermolist = huge(0)
723  molecule_set => molecules%els
724  mregions = nregions
725  DO ig = 1, mregions
726  CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, n_rep_val=n_rep)
727  DO jg = 1, n_rep
728  CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, i_rep_val=jg, i_vals=tmplist)
729  DO i = 1, SIZE(tmplist)
730  ipart = tmplist(i)
731  cpassert(((ipart > 0) .AND. (ipart <= particles%n_els)))
732  IF (thermolist(ipart) == huge(0)) THEN
733  thermolist(ipart) = ig
734  ELSE
735  cpabort("")
736  END IF
737  END DO
738  END DO
739  CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, n_rep_val=n_rep)
740  DO jg = 1, n_rep
741  CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, i_rep_val=jg, c_vals=tmpstringlist)
742  DO ilist = 1, SIZE(tmpstringlist)
743  DO ikind = 1, SIZE(molecule_kind_set)
744  molecule_kind => molecule_kind_set(ikind)
745  IF (molecule_kind%name == tmpstringlist(ilist)) THEN
746  DO imol = 1, SIZE(molecule_kind%molecule_list)
747  molecule => molecule_set(molecule_kind%molecule_list(imol))
748  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
749  DO ipart = first_atom, last_atom
750  IF (thermolist(ipart) == huge(0)) THEN
751  thermolist(ipart) = ig
752  ELSE
753  cpabort("")
754  END IF
755  END DO
756  END DO
757  END IF
758  END DO
759  END DO
760  END DO
761  CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
762  subsys_qm=.false., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
763  CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
764  subsys_qm=.true., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
765  END DO
766 
767  ! Dump IO warning for not thermalized particles
768  IF (any(thermolist == huge(0))) THEN
769  nregions = nregions + 1
770  sum_of_thermostats = sum_of_thermostats + 1
771  ALLOCATE (tmp(count(thermolist == huge(0))))
772  ilist = 0
773  DO i = 1, SIZE(thermolist)
774  IF (thermolist(i) == huge(0)) THEN
775  ilist = ilist + 1
776  tmp(ilist) = i
777  thermolist(i) = nregions
778  END IF
779  END DO
780  IF (output_unit > 0) THEN
781  WRITE (output_unit, '(A)') " WARNING| No thermostats defined for the following atoms:"
782  IF (ilist > 0) WRITE (output_unit, '(" WARNING|",9I8)') tmp
783  WRITE (output_unit, '(A)') " WARNING| They will be included in a further unique thermostat!"
784  END IF
785  DEALLOCATE (tmp)
786  END IF
787  cpassert(all(thermolist /= huge(0)))
788 
789  ! Now identify the local number of thermostats
790  ALLOCATE (tmp(nregions))
791  tmp = 0
792  natom_local = 0
793  DO ikind = 1, SIZE(molecule_kind_set)
794  nmol_per_kind = local_molecules%n_el(ikind)
795  DO imol = 1, nmol_per_kind
796  i = local_molecules%list(ikind)%array(imol)
797  molecule => molecule_set(i)
798  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
799  DO ipart = first_atom, last_atom
800  natom_local = natom_local + 1
801  tmp(thermolist(ipart)) = 1
802  END DO
803  END DO
804  END DO
805  number = sum(tmp)
806  DEALLOCATE (tmp)
807 
808  ! Now map the local atoms with the corresponding thermostat
809  ALLOCATE (map_loc_thermo_gen(natom_local))
810  natom_local = 0
811  DO ikind = 1, SIZE(molecule_kind_set)
812  nmol_per_kind = local_molecules%n_el(ikind)
813  DO imol = 1, nmol_per_kind
814  i = local_molecules%list(ikind)%array(imol)
815  molecule => molecule_set(i)
816  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
817  DO ipart = first_atom, last_atom
818  natom_local = natom_local + 1
819  map_loc_thermo_gen(natom_local) = thermolist(ipart)
820  END DO
821  END DO
822  END DO
823 
824  DEALLOCATE (thermolist)
825  END SUBROUTINE get_defined_region_info
826 
827 ! **************************************************************************************************
828 !> \brief ...
829 !> \param region_sections ...
830 !> \param qmmm_env ...
831 !> \param thermolist ...
832 !> \param molecule_set ...
833 !> \param subsys_qm ...
834 !> \param ig ...
835 !> \param sum_of_thermostats ...
836 !> \param nregions ...
837 !> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
838 ! **************************************************************************************************
839  SUBROUTINE setup_thermostat_subsys(region_sections, qmmm_env, thermolist, &
840  molecule_set, subsys_qm, ig, sum_of_thermostats, nregions)
841  TYPE(section_vals_type), POINTER :: region_sections
842  TYPE(qmmm_env_type), POINTER :: qmmm_env
843  INTEGER, DIMENSION(:), POINTER :: thermolist
844  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
845  LOGICAL, INTENT(IN) :: subsys_qm
846  INTEGER, INTENT(IN) :: ig
847  INTEGER, INTENT(INOUT) :: sum_of_thermostats, nregions
848 
849  CHARACTER(LEN=default_string_length) :: label1, label2
850  INTEGER :: first_atom, i, imolecule, ipart, &
851  last_atom, nrep, thermo1
852  INTEGER, DIMENSION(:), POINTER :: atom_index1
853  LOGICAL :: explicit
854  TYPE(molecule_type), POINTER :: molecule
855 
856  label1 = "MM_SUBSYS"
857  label2 = "QM_SUBSYS"
858  IF (subsys_qm) THEN
859  label1 = "QM_SUBSYS"
860  label2 = "MM_SUBSYS"
861  END IF
862  CALL section_vals_val_get(region_sections, trim(label1), i_rep_section=ig, &
863  n_rep_val=nrep, explicit=explicit)
864  IF (nrep == 1 .AND. explicit) THEN
865  IF (ASSOCIATED(qmmm_env)) THEN
866  atom_index1 => qmmm_env%qm%mm_atom_index
867  IF (subsys_qm) THEN
868  atom_index1 => qmmm_env%qm%qm_atom_index
869  END IF
870  CALL section_vals_val_get(region_sections, trim(label1), i_val=thermo1, i_rep_section=ig)
871  SELECT CASE (thermo1)
872  CASE (do_constr_atomic)
873  DO i = 1, SIZE(atom_index1)
874  ipart = atom_index1(i)
875  IF (subsys_qm .AND. qmmm_env%qm%qmmm_link .AND. ASSOCIATED(qmmm_env%qm%mm_link_atoms)) THEN
876  IF (any(ipart == qmmm_env%qm%mm_link_atoms)) cycle
877  END IF
878  IF (thermolist(ipart) == huge(0)) THEN
879  thermolist(ipart) = ig
880  ELSE
881  CALL cp_abort(__location__, &
882  'One atom ('//cp_to_string(ipart)//') of the '// &
883  trim(label1)//' was already assigned to'// &
884  ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
885  '. Please check the input for inconsistencies!')
886  END IF
887  END DO
888  CASE (do_constr_molec)
889  DO imolecule = 1, SIZE(molecule_set)
890  molecule => molecule_set(imolecule)
891  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
892  IF (any(atom_index1 >= first_atom .AND. atom_index1 <= last_atom)) THEN
893  DO ipart = first_atom, last_atom
894  IF (thermolist(ipart) == huge(0)) THEN
895  thermolist(ipart) = ig
896  ELSE
897  CALL cp_abort(__location__, &
898  'One atom ('//cp_to_string(ipart)//') of the '// &
899  trim(label1)//' was already assigned to'// &
900  ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
901  '. Please check the input for inconsistencies!')
902  END IF
903  END DO
904  END IF
905  END DO
906  END SELECT
907  ELSE
908  sum_of_thermostats = sum_of_thermostats - 1
909  nregions = nregions - 1
910  END IF
911  END IF
912  END SUBROUTINE setup_thermostat_subsys
913 
914 ! **************************************************************************************************
915 !> \brief ...
916 !> \param map_info ...
917 !> \param npt ...
918 !> \param group ...
919 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
920 ! **************************************************************************************************
921  SUBROUTINE ke_region_baro(map_info, npt, group)
922  TYPE(map_info_type), POINTER :: map_info
923  TYPE(npt_info_type), DIMENSION(:, :), &
924  INTENT(INOUT) :: npt
925  TYPE(mp_comm_type), INTENT(IN) :: group
926 
927  INTEGER :: i, j, ncoef
928 
929  map_info%v_scale = 1.0_dp
930  map_info%s_kin = 0.0_dp
931  ncoef = 0
932  DO i = 1, SIZE(npt, 1)
933  DO j = 1, SIZE(npt, 2)
934  ncoef = ncoef + 1
935  map_info%p_kin(1, ncoef)%point = map_info%p_kin(1, ncoef)%point &
936  + npt(i, j)%mass*npt(i, j)%v**2
937  END DO
938  END DO
939 
940  IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
941 
942  END SUBROUTINE ke_region_baro
943 
944 ! **************************************************************************************************
945 !> \brief ...
946 !> \param map_info ...
947 !> \param npt ...
948 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
949 ! **************************************************************************************************
950  SUBROUTINE vel_rescale_baro(map_info, npt)
951  TYPE(map_info_type), POINTER :: map_info
952  TYPE(npt_info_type), DIMENSION(:, :), &
953  INTENT(INOUT) :: npt
954 
955  INTEGER :: i, j, ncoef
956 
957  ncoef = 0
958  DO i = 1, SIZE(npt, 1)
959  DO j = 1, SIZE(npt, 2)
960  ncoef = ncoef + 1
961  npt(i, j)%v = npt(i, j)%v*map_info%p_scale(1, ncoef)%point
962  END DO
963  END DO
964 
965  END SUBROUTINE vel_rescale_baro
966 
967 ! **************************************************************************************************
968 !> \brief ...
969 !> \param map_info ...
970 !> \param particle_set ...
971 !> \param molecule_kind_set ...
972 !> \param local_molecules ...
973 !> \param molecule_set ...
974 !> \param group ...
975 !> \param vel ...
976 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
977 ! **************************************************************************************************
978  SUBROUTINE ke_region_particles(map_info, particle_set, molecule_kind_set, &
979  local_molecules, molecule_set, group, vel)
980 
981  TYPE(map_info_type), POINTER :: map_info
982  TYPE(particle_type), POINTER :: particle_set(:)
983  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
984  TYPE(distribution_1d_type), POINTER :: local_molecules
985  TYPE(molecule_type), POINTER :: molecule_set(:)
986  TYPE(mp_comm_type), INTENT(IN) :: group
987  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
988 
989  INTEGER :: first_atom, ii, ikind, imol, imol_local, &
990  ipart, last_atom, nmol_local
991  LOGICAL :: present_vel
992  REAL(kind=dp) :: mass
993  TYPE(atomic_kind_type), POINTER :: atomic_kind
994  TYPE(molecule_type), POINTER :: molecule
995 
996  map_info%v_scale = 1.0_dp
997  map_info%s_kin = 0.0_dp
998  present_vel = PRESENT(vel)
999  ii = 0
1000  DO ikind = 1, SIZE(molecule_kind_set)
1001  nmol_local = local_molecules%n_el(ikind)
1002  DO imol_local = 1, nmol_local
1003  imol = local_molecules%list(ikind)%array(imol_local)
1004  molecule => molecule_set(imol)
1005  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1006  DO ipart = first_atom, last_atom
1007  ii = ii + 1
1008  atomic_kind => particle_set(ipart)%atomic_kind
1009  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1010  IF (present_vel) THEN
1011  IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
1012  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*vel(1, ipart)**2
1013  IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
1014  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*vel(2, ipart)**2
1015  IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
1016  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*vel(3, ipart)**2
1017  ELSE
1018  IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
1019  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*particle_set(ipart)%v(1)**2
1020  IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
1021  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*particle_set(ipart)%v(2)**2
1022  IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
1023  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*particle_set(ipart)%v(3)**2
1024  END IF
1025  END DO
1026  END DO
1027  END DO
1028 
1029  IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1030 
1031  END SUBROUTINE ke_region_particles
1032 
1033 ! **************************************************************************************************
1034 !> \brief ...
1035 !> \param map_info ...
1036 !> \param particle_set ...
1037 !> \param molecule_kind_set ...
1038 !> \param local_molecules ...
1039 !> \param molecule_set ...
1040 !> \param group ...
1041 !> \param vel ...
1042 !> \author 07.2009 MI
1043 ! **************************************************************************************************
1044  SUBROUTINE momentum_region_particles(map_info, particle_set, molecule_kind_set, &
1045  local_molecules, molecule_set, group, vel)
1046 
1047  TYPE(map_info_type), POINTER :: map_info
1048  TYPE(particle_type), POINTER :: particle_set(:)
1049  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1050  TYPE(distribution_1d_type), POINTER :: local_molecules
1051  TYPE(molecule_type), POINTER :: molecule_set(:)
1052  TYPE(mp_comm_type), INTENT(IN) :: group
1053  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
1054 
1055  INTEGER :: first_atom, ii, ikind, imol, imol_local, &
1056  ipart, last_atom, nmol_local
1057  LOGICAL :: present_vel
1058  REAL(kind=dp) :: mass
1059  TYPE(atomic_kind_type), POINTER :: atomic_kind
1060  TYPE(molecule_type), POINTER :: molecule
1061 
1062  map_info%v_scale = 1.0_dp
1063  map_info%s_kin = 0.0_dp
1064  present_vel = PRESENT(vel)
1065  ii = 0
1066  DO ikind = 1, SIZE(molecule_kind_set)
1067  nmol_local = local_molecules%n_el(ikind)
1068  DO imol_local = 1, nmol_local
1069  imol = local_molecules%list(ikind)%array(imol_local)
1070  molecule => molecule_set(imol)
1071  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1072  DO ipart = first_atom, last_atom
1073  ii = ii + 1
1074  atomic_kind => particle_set(ipart)%atomic_kind
1075  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1076  IF (present_vel) THEN
1077  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + sqrt(mass)*vel(1, ipart)
1078  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + sqrt(mass)*vel(2, ipart)
1079  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + sqrt(mass)*vel(3, ipart)
1080  ELSE
1081  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + sqrt(mass)*particle_set(ipart)%v(1)
1082  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + sqrt(mass)*particle_set(ipart)%v(2)
1083  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + sqrt(mass)*particle_set(ipart)%v(3)
1084  END IF
1085  END DO
1086  END DO
1087  END DO
1088 
1089  IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1090 
1091  END SUBROUTINE momentum_region_particles
1092 
1093 ! **************************************************************************************************
1094 !> \brief ...
1095 !> \param map_info ...
1096 !> \param molecule_kind_set ...
1097 !> \param molecule_set ...
1098 !> \param particle_set ...
1099 !> \param local_molecules ...
1100 !> \param shell_adiabatic ...
1101 !> \param shell_particle_set ...
1102 !> \param core_particle_set ...
1103 !> \param vel ...
1104 !> \param shell_vel ...
1105 !> \param core_vel ...
1106 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1107 ! **************************************************************************************************
1108  SUBROUTINE vel_rescale_particles(map_info, molecule_kind_set, molecule_set, &
1109  particle_set, local_molecules, shell_adiabatic, shell_particle_set, &
1110  core_particle_set, vel, shell_vel, core_vel)
1111 
1112  TYPE(map_info_type), POINTER :: map_info
1113  TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1114  TYPE(molecule_type), POINTER :: molecule_set(:)
1115  TYPE(particle_type), POINTER :: particle_set(:)
1116  TYPE(distribution_1d_type), POINTER :: local_molecules
1117  LOGICAL, INTENT(IN) :: shell_adiabatic
1118  TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
1119  core_particle_set(:)
1120  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
1121  core_vel(:, :)
1122 
1123  INTEGER :: first_atom, ii, ikind, imol, imol_local, &
1124  ipart, jj, last_atom, nmol_local, &
1125  shell_index
1126  LOGICAL :: present_vel
1127  REAL(kind=dp) :: fac_massc, fac_masss, mass, vc(3), vs(3)
1128  TYPE(atomic_kind_type), POINTER :: atomic_kind
1129  TYPE(molecule_type), POINTER :: molecule
1130  TYPE(shell_kind_type), POINTER :: shell
1131 
1132  ii = 0
1133  jj = 0
1134  present_vel = PRESENT(vel)
1135  ! Just few checks for consistency
1136  IF (present_vel) THEN
1137  IF (shell_adiabatic) THEN
1138  cpassert(PRESENT(shell_vel))
1139  cpassert(PRESENT(core_vel))
1140  END IF
1141  ELSE
1142  IF (shell_adiabatic) THEN
1143  cpassert(PRESENT(shell_particle_set))
1144  cpassert(PRESENT(core_particle_set))
1145  END IF
1146  END IF
1147  kind: DO ikind = 1, SIZE(molecule_kind_set)
1148  nmol_local = local_molecules%n_el(ikind)
1149  mol_local: DO imol_local = 1, nmol_local
1150  imol = local_molecules%list(ikind)%array(imol_local)
1151  molecule => molecule_set(imol)
1152  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1153  particle: DO ipart = first_atom, last_atom
1154  ii = ii + 1
1155  IF (present_vel) THEN
1156  vel(1, ipart) = vel(1, ipart)*map_info%p_scale(1, ii)%point
1157  vel(2, ipart) = vel(2, ipart)*map_info%p_scale(2, ii)%point
1158  vel(3, ipart) = vel(3, ipart)*map_info%p_scale(3, ii)%point
1159  ELSE
1160  particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*map_info%p_scale(1, ii)%point
1161  particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*map_info%p_scale(2, ii)%point
1162  particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*map_info%p_scale(3, ii)%point
1163  END IF
1164  ! If Shell Adiabatic then apply the NHC thermostat also to the Shells
1165  IF (shell_adiabatic) THEN
1166  shell_index = particle_set(ipart)%shell_index
1167  IF (shell_index /= 0) THEN
1168  jj = jj + 2
1169  atomic_kind => particle_set(ipart)%atomic_kind
1170  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
1171  fac_masss = shell%mass_shell/mass
1172  fac_massc = shell%mass_core/mass
1173  IF (present_vel) THEN
1174  vs(1:3) = shell_vel(1:3, shell_index)
1175  vc(1:3) = core_vel(1:3, shell_index)
1176  shell_vel(1, shell_index) = vel(1, ipart) + fac_massc*(vs(1) - vc(1))
1177  shell_vel(2, shell_index) = vel(2, ipart) + fac_massc*(vs(2) - vc(2))
1178  shell_vel(3, shell_index) = vel(3, ipart) + fac_massc*(vs(3) - vc(3))
1179  core_vel(1, shell_index) = vel(1, ipart) + fac_masss*(vc(1) - vs(1))
1180  core_vel(2, shell_index) = vel(2, ipart) + fac_masss*(vc(2) - vs(2))
1181  core_vel(3, shell_index) = vel(3, ipart) + fac_masss*(vc(3) - vs(3))
1182  ELSE
1183  vs(1:3) = shell_particle_set(shell_index)%v(1:3)
1184  vc(1:3) = core_particle_set(shell_index)%v(1:3)
1185  shell_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_massc*(vs(1) - vc(1))
1186  shell_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_massc*(vs(2) - vc(2))
1187  shell_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_massc*(vs(3) - vc(3))
1188  core_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_masss*(vc(1) - vs(1))
1189  core_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_masss*(vc(2) - vs(2))
1190  core_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_masss*(vc(3) - vs(3))
1191  END IF
1192  END IF
1193  END IF
1194  END DO particle
1195  END DO mol_local
1196  END DO kind
1197 
1198  END SUBROUTINE vel_rescale_particles
1199 
1200 ! **************************************************************************************************
1201 !> \brief ...
1202 !> \param map_info ...
1203 !> \param particle_set ...
1204 !> \param atomic_kind_set ...
1205 !> \param local_particles ...
1206 !> \param group ...
1207 !> \param core_particle_set ...
1208 !> \param shell_particle_set ...
1209 !> \param core_vel ...
1210 !> \param shell_vel ...
1211 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1212 ! **************************************************************************************************
1213  SUBROUTINE ke_region_shells(map_info, particle_set, atomic_kind_set, &
1214  local_particles, group, core_particle_set, shell_particle_set, &
1215  core_vel, shell_vel)
1216 
1217  TYPE(map_info_type), POINTER :: map_info
1218  TYPE(particle_type), POINTER :: particle_set(:)
1219  TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
1220  TYPE(distribution_1d_type), POINTER :: local_particles
1221  TYPE(mp_comm_type), INTENT(IN) :: group
1222  TYPE(particle_type), OPTIONAL, POINTER :: core_particle_set(:), &
1223  shell_particle_set(:)
1224  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: core_vel(:, :), shell_vel(:, :)
1225 
1226  INTEGER :: ii, iparticle, iparticle_kind, &
1227  iparticle_local, nparticle_kind, &
1228  nparticle_local, shell_index
1229  LOGICAL :: is_shell, present_vel
1230  REAL(dp) :: mass, mu_mass, v_sc(3)
1231  TYPE(atomic_kind_type), POINTER :: atomic_kind
1232  TYPE(shell_kind_type), POINTER :: shell
1233 
1234  present_vel = PRESENT(shell_vel)
1235  ! Preliminary checks for consistency usage
1236  IF (present_vel) THEN
1237  cpassert(PRESENT(core_vel))
1238  ELSE
1239  cpassert(PRESENT(shell_particle_set))
1240  cpassert(PRESENT(core_particle_set))
1241  END IF
1242  ! get force on first thermostat for all the chains in the system.
1243  map_info%v_scale = 1.0_dp
1244  map_info%s_kin = 0.0_dp
1245  ii = 0
1246 
1247  nparticle_kind = SIZE(atomic_kind_set)
1248  DO iparticle_kind = 1, nparticle_kind
1249  atomic_kind => atomic_kind_set(iparticle_kind)
1250  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
1251  IF (is_shell) THEN
1252  mu_mass = shell%mass_shell*shell%mass_core/mass
1253  nparticle_local = local_particles%n_el(iparticle_kind)
1254  DO iparticle_local = 1, nparticle_local
1255  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1256  shell_index = particle_set(iparticle)%shell_index
1257  ii = ii + 1
1258  IF (present_vel) THEN
1259  v_sc(1) = core_vel(1, shell_index) - shell_vel(1, shell_index)
1260  v_sc(2) = core_vel(2, shell_index) - shell_vel(2, shell_index)
1261  v_sc(3) = core_vel(3, shell_index) - shell_vel(3, shell_index)
1262  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
1263  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
1264  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
1265  ELSE
1266  v_sc(1) = core_particle_set(shell_index)%v(1) - shell_particle_set(shell_index)%v(1)
1267  v_sc(2) = core_particle_set(shell_index)%v(2) - shell_particle_set(shell_index)%v(2)
1268  v_sc(3) = core_particle_set(shell_index)%v(3) - shell_particle_set(shell_index)%v(3)
1269  map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
1270  map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
1271  map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
1272  END IF
1273  END DO
1274  END IF
1275  END DO
1276  IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1277 
1278  END SUBROUTINE ke_region_shells
1279 
1280 ! **************************************************************************************************
1281 !> \brief ...
1282 !> \param map_info ...
1283 !> \param atomic_kind_set ...
1284 !> \param particle_set ...
1285 !> \param local_particles ...
1286 !> \param shell_particle_set ...
1287 !> \param core_particle_set ...
1288 !> \param shell_vel ...
1289 !> \param core_vel ...
1290 !> \param vel ...
1291 !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1292 ! **************************************************************************************************
1293  SUBROUTINE vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
1294  shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
1295 
1296  TYPE(map_info_type), POINTER :: map_info
1297  TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
1298  TYPE(particle_type), POINTER :: particle_set(:)
1299  TYPE(distribution_1d_type), POINTER :: local_particles
1300  TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
1301  core_particle_set(:)
1302  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: shell_vel(:, :), core_vel(:, :), &
1303  vel(:, :)
1304 
1305  INTEGER :: ii, iparticle, iparticle_kind, &
1306  iparticle_local, nparticle_kind, &
1307  nparticle_local, shell_index
1308  LOGICAL :: is_shell, present_vel
1309  REAL(dp) :: mass, massc, masss, umass, v(3), vc(3), &
1310  vs(3)
1311  TYPE(atomic_kind_type), POINTER :: atomic_kind
1312  TYPE(shell_kind_type), POINTER :: shell
1313 
1314  present_vel = PRESENT(vel)
1315  ! Preliminary checks for consistency usage
1316  IF (present_vel) THEN
1317  cpassert(PRESENT(shell_vel))
1318  cpassert(PRESENT(core_vel))
1319  ELSE
1320  cpassert(PRESENT(shell_particle_set))
1321  cpassert(PRESENT(core_particle_set))
1322  END IF
1323  ii = 0
1324  nparticle_kind = SIZE(atomic_kind_set)
1325  ! now scale the core-shell velocities
1326  kind: DO iparticle_kind = 1, nparticle_kind
1327  atomic_kind => atomic_kind_set(iparticle_kind)
1328  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
1329  IF (is_shell) THEN
1330  umass = 1.0_dp/mass
1331  masss = shell%mass_shell*umass
1332  massc = shell%mass_core*umass
1333 
1334  nparticle_local = local_particles%n_el(iparticle_kind)
1335  particles: DO iparticle_local = 1, nparticle_local
1336  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1337  shell_index = particle_set(iparticle)%shell_index
1338  ii = ii + 1
1339  IF (present_vel) THEN
1340  vc(1:3) = core_vel(1:3, shell_index)
1341  vs(1:3) = shell_vel(1:3, shell_index)
1342  v(1:3) = vel(1:3, iparticle)
1343  shell_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
1344  shell_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
1345  shell_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
1346  core_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
1347  core_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
1348  core_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
1349  ELSE
1350  vc(1:3) = core_particle_set(shell_index)%v(1:3)
1351  vs(1:3) = shell_particle_set(shell_index)%v(1:3)
1352  v(1:3) = particle_set(iparticle)%v(1:3)
1353  shell_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
1354  shell_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
1355  shell_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
1356  core_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
1357  core_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
1358  core_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
1359  END IF
1360  END DO particles
1361  END IF
1362  END DO kind
1363 
1364  END SUBROUTINE vel_rescale_shells
1365 
1366 ! **************************************************************************************************
1367 !> \brief Calculates kinetic energy and potential energy of the nhc variables
1368 !> \param nhc ...
1369 !> \param nhc_pot ...
1370 !> \param nhc_kin ...
1371 !> \param para_env ...
1372 !> \param array_kin ...
1373 !> \param array_pot ...
1374 !> \par History
1375 !> none
1376 !> \author CJM
1377 ! **************************************************************************************************
1378  SUBROUTINE get_nhc_energies(nhc, nhc_pot, nhc_kin, para_env, array_kin, array_pot)
1379  TYPE(lnhc_parameters_type), POINTER :: nhc
1380  REAL(kind=dp), INTENT(OUT) :: nhc_pot, nhc_kin
1381  TYPE(mp_para_env_type), POINTER :: para_env
1382  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_kin, array_pot
1383 
1384  INTEGER :: imap, l, n, number
1385  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin, vpot
1386 
1387  number = nhc%glob_num_nhc
1388  ALLOCATE (akin(number))
1389  ALLOCATE (vpot(number))
1390  akin = 0.0_dp
1391  vpot = 0.0_dp
1392  DO n = 1, nhc%loc_num_nhc
1393  imap = nhc%map_info%index(n)
1394  DO l = 1, nhc%nhc_len
1395  akin(imap) = akin(imap) + 0.5_dp*nhc%nvt(l, n)%mass*nhc%nvt(l, n)%v**2
1396  vpot(imap) = vpot(imap) + nhc%nvt(l, n)%nkt*nhc%nvt(l, n)%eta
1397  END DO
1398  END DO
1399 
1400  ! Handle the thermostat distribution
1401  IF (nhc%map_info%dis_type == do_thermo_no_communication) THEN
1402  CALL para_env%sum(akin)
1403  CALL para_env%sum(vpot)
1404  ELSE IF (nhc%map_info%dis_type == do_thermo_communication) THEN
1405  CALL communication_thermo_low1(akin, number, para_env)
1406  CALL communication_thermo_low1(vpot, number, para_env)
1407  END IF
1408  nhc_kin = sum(akin)
1409  nhc_pot = sum(vpot)
1410 
1411  ! Possibly give back kinetic or potential energy arrays
1412  IF (PRESENT(array_pot)) THEN
1413  IF (ASSOCIATED(array_pot)) THEN
1414  cpassert(SIZE(array_pot) == number)
1415  ELSE
1416  ALLOCATE (array_pot(number))
1417  END IF
1418  array_pot = vpot
1419  END IF
1420  IF (PRESENT(array_kin)) THEN
1421  IF (ASSOCIATED(array_kin)) THEN
1422  cpassert(SIZE(array_kin) == number)
1423  ELSE
1424  ALLOCATE (array_kin(number))
1425  END IF
1426  array_kin = akin
1427  END IF
1428  DEALLOCATE (akin)
1429  DEALLOCATE (vpot)
1430  END SUBROUTINE get_nhc_energies
1431 
1432 ! **************************************************************************************************
1433 !> \brief Calculates kinetic energy and potential energy
1434 !> of the csvr and gle thermostats
1435 !> \param map_info ...
1436 !> \param loc_num ...
1437 !> \param glob_num ...
1438 !> \param thermo_energy ...
1439 !> \param thermostat_kin ...
1440 !> \param para_env ...
1441 !> \param array_pot ...
1442 !> \param array_kin ...
1443 !> \par History generalized MI [07.2009]
1444 !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1445 ! **************************************************************************************************
1446  SUBROUTINE get_kin_energies(map_info, loc_num, glob_num, thermo_energy, thermostat_kin, &
1447  para_env, array_pot, array_kin)
1448 
1449  TYPE(map_info_type), POINTER :: map_info
1450  INTEGER, INTENT(IN) :: loc_num, glob_num
1451  REAL(dp), DIMENSION(:), INTENT(IN) :: thermo_energy
1452  REAL(kind=dp), INTENT(OUT) :: thermostat_kin
1453  TYPE(mp_para_env_type), POINTER :: para_env
1454  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_pot, array_kin
1455 
1456  INTEGER :: imap, n, number
1457  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin
1458 
1459  number = glob_num
1460  ALLOCATE (akin(number))
1461  akin = 0.0_dp
1462  DO n = 1, loc_num
1463  imap = map_info%index(n)
1464  akin(imap) = thermo_energy(n)
1465  END DO
1466 
1467  ! Handle the thermostat distribution
1468  IF (map_info%dis_type == do_thermo_no_communication) THEN
1469  CALL para_env%sum(akin)
1470  ELSE IF (map_info%dis_type == do_thermo_communication) THEN
1471  CALL communication_thermo_low1(akin, number, para_env)
1472  END IF
1473  thermostat_kin = sum(akin)
1474 
1475  ! Possibly give back kinetic or potential energy arrays
1476  IF (PRESENT(array_pot)) THEN
1477  IF (ASSOCIATED(array_pot)) THEN
1478  cpassert(SIZE(array_pot) == number)
1479  ELSE
1480  ALLOCATE (array_pot(number))
1481  END IF
1482  array_pot = 0.0_dp
1483  END IF
1484  IF (PRESENT(array_kin)) THEN
1485  IF (ASSOCIATED(array_kin)) THEN
1486  cpassert(SIZE(array_kin) == number)
1487  ELSE
1488  ALLOCATE (array_kin(number))
1489  END IF
1490  array_kin = akin
1491  END IF
1492  DEALLOCATE (akin)
1493  END SUBROUTINE get_kin_energies
1494 
1495 ! **************************************************************************************************
1496 !> \brief Calculates the temperatures of the regions when a thermostat is
1497 !> applied
1498 !> \param map_info ...
1499 !> \param loc_num ...
1500 !> \param glob_num ...
1501 !> \param nkt ...
1502 !> \param dof ...
1503 !> \param para_env ...
1504 !> \param temp_tot ...
1505 !> \param array_temp ...
1506 !> \par History generalized MI [07.2009]
1507 !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1508 ! **************************************************************************************************
1509  SUBROUTINE get_temperatures(map_info, loc_num, glob_num, nkt, dof, para_env, &
1510  temp_tot, array_temp)
1511  TYPE(map_info_type), POINTER :: map_info
1512  INTEGER, INTENT(IN) :: loc_num, glob_num
1513  REAL(dp), DIMENSION(:), INTENT(IN) :: nkt, dof
1514  TYPE(mp_para_env_type), POINTER :: para_env
1515  REAL(kind=dp), INTENT(OUT) :: temp_tot
1516  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_temp
1517 
1518  INTEGER :: i, imap, imap2, n, number
1519  REAL(kind=dp) :: fdeg_of_free
1520  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin, deg_of_free
1521 
1522  number = glob_num
1523  ALLOCATE (akin(number))
1524  ALLOCATE (deg_of_free(number))
1525  akin = 0.0_dp
1526  deg_of_free = 0.0_dp
1527  DO n = 1, loc_num
1528  imap = map_info%index(n)
1529  imap2 = map_info%map_index(n)
1530  IF (nkt(n) == 0.0_dp) cycle
1531  deg_of_free(imap) = real(dof(n), kind=dp)
1532  akin(imap) = map_info%s_kin(imap2)
1533  END DO
1534 
1535  ! Handle the thermostat distribution
1536  IF (map_info%dis_type == do_thermo_no_communication) THEN
1537  CALL para_env%sum(akin)
1538  CALL para_env%sum(deg_of_free)
1539  ELSE IF (map_info%dis_type == do_thermo_communication) THEN
1540  CALL communication_thermo_low1(akin, number, para_env)
1541  CALL communication_thermo_low1(deg_of_free, number, para_env)
1542  END IF
1543  temp_tot = sum(akin)
1544  fdeg_of_free = sum(deg_of_free)
1545 
1546  temp_tot = temp_tot/fdeg_of_free
1547  temp_tot = cp_unit_from_cp2k(temp_tot, "K_temp")
1548  ! Possibly give back temperatures of the full set of regions
1549  IF (PRESENT(array_temp)) THEN
1550  IF (ASSOCIATED(array_temp)) THEN
1551  cpassert(SIZE(array_temp) == number)
1552  ELSE
1553  ALLOCATE (array_temp(number))
1554  END IF
1555  DO i = 1, number
1556  array_temp(i) = akin(i)/deg_of_free(i)
1557  array_temp(i) = cp_unit_from_cp2k(array_temp(i), "K_temp")
1558  END DO
1559  END IF
1560  DEALLOCATE (akin)
1561  DEALLOCATE (deg_of_free)
1562  END SUBROUTINE get_temperatures
1563 
1564 ! **************************************************************************************************
1565 !> \brief Calculates energy associated with a thermostat
1566 !> \param thermostat ...
1567 !> \param thermostat_pot ...
1568 !> \param thermostat_kin ...
1569 !> \param para_env ...
1570 !> \param array_pot ...
1571 !> \param array_kin ...
1572 !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1573 ! **************************************************************************************************
1574  SUBROUTINE get_thermostat_energies(thermostat, thermostat_pot, thermostat_kin, para_env, &
1575  array_pot, array_kin)
1576  TYPE(thermostat_type), POINTER :: thermostat
1577  REAL(kind=dp), INTENT(OUT) :: thermostat_pot, thermostat_kin
1578  TYPE(mp_para_env_type), POINTER :: para_env
1579  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_pot, array_kin
1580 
1581  INTEGER :: i
1582  REAL(dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1583 
1584  thermostat_pot = 0.0_dp
1585  thermostat_kin = 0.0_dp
1586  IF (ASSOCIATED(thermostat)) THEN
1587  IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
1588  ! Energy associated with the Nose-Hoover thermostat
1589  cpassert(ASSOCIATED(thermostat%nhc))
1590  CALL get_nhc_energies(thermostat%nhc, thermostat_pot, thermostat_kin, para_env, &
1591  array_pot, array_kin)
1592  ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
1593  ! Energy associated with the CSVR thermostat
1594  cpassert(ASSOCIATED(thermostat%csvr))
1595  ALLOCATE (thermo_energy(thermostat%csvr%loc_num_csvr))
1596  DO i = 1, thermostat%csvr%loc_num_csvr
1597  thermo_energy(i) = thermostat%csvr%nvt(i)%thermostat_energy
1598  END DO
1599  CALL get_kin_energies(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
1600  thermostat%csvr%glob_num_csvr, thermo_energy, &
1601  thermostat_kin, para_env, array_pot, array_kin)
1602  DEALLOCATE (thermo_energy)
1603 
1604  ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
1605  ! Energy associated with the GLE thermostat
1606  cpassert(ASSOCIATED(thermostat%gle))
1607  ALLOCATE (thermo_energy(thermostat%gle%loc_num_gle))
1608  DO i = 1, thermostat%gle%loc_num_gle
1609  thermo_energy(i) = thermostat%gle%nvt(i)%thermostat_energy
1610  END DO
1611  CALL get_kin_energies(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
1612  thermostat%gle%glob_num_gle, thermo_energy, &
1613  thermostat_kin, para_env, array_pot, array_kin)
1614  DEALLOCATE (thermo_energy)
1615 
1616  ![NB] nothing to do for Ad-Langevin?
1617 
1618  END IF
1619  END IF
1620 
1621  END SUBROUTINE get_thermostat_energies
1622 
1623 ! **************************************************************************************************
1624 !> \brief Calculates the temperatures for each region associated to a thermostat
1625 !> \param thermostat ...
1626 !> \param tot_temperature ...
1627 !> \param para_env ...
1628 !> \param array_temp ...
1629 !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1630 ! **************************************************************************************************
1631  SUBROUTINE get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
1632  TYPE(thermostat_type), POINTER :: thermostat
1633  REAL(kind=dp), INTENT(OUT) :: tot_temperature
1634  TYPE(mp_para_env_type), POINTER :: para_env
1635  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_temp
1636 
1637  INTEGER :: i
1638  REAL(dp), ALLOCATABLE, DIMENSION(:) :: dof, nkt
1639 
1640  IF (ASSOCIATED(thermostat)) THEN
1641  IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
1642  ! Energy associated with the Nose-Hoover thermostat
1643  cpassert(ASSOCIATED(thermostat%nhc))
1644  ALLOCATE (nkt(thermostat%nhc%loc_num_nhc))
1645  ALLOCATE (dof(thermostat%nhc%loc_num_nhc))
1646  DO i = 1, thermostat%nhc%loc_num_nhc
1647  nkt(i) = thermostat%nhc%nvt(1, i)%nkt
1648  dof(i) = real(thermostat%nhc%nvt(1, i)%degrees_of_freedom, kind=dp)
1649  END DO
1650  CALL get_temperatures(thermostat%nhc%map_info, thermostat%nhc%loc_num_nhc, &
1651  thermostat%nhc%glob_num_nhc, nkt, dof, para_env, tot_temperature, array_temp)
1652  DEALLOCATE (nkt)
1653  DEALLOCATE (dof)
1654  ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
1655  ! Energy associated with the CSVR thermostat
1656  cpassert(ASSOCIATED(thermostat%csvr))
1657 
1658  ALLOCATE (nkt(thermostat%csvr%loc_num_csvr))
1659  ALLOCATE (dof(thermostat%csvr%loc_num_csvr))
1660  DO i = 1, thermostat%csvr%loc_num_csvr
1661  nkt(i) = thermostat%csvr%nvt(i)%nkt
1662  dof(i) = real(thermostat%csvr%nvt(i)%degrees_of_freedom, kind=dp)
1663  END DO
1664  CALL get_temperatures(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
1665  thermostat%csvr%glob_num_csvr, nkt, dof, para_env, tot_temperature, array_temp)
1666  DEALLOCATE (nkt)
1667  DEALLOCATE (dof)
1668  ELSE IF (thermostat%type_of_thermostat == do_thermo_al) THEN
1669  ! Energy associated with the AD_LANGEVIN thermostat
1670  cpassert(ASSOCIATED(thermostat%al))
1671 
1672  ALLOCATE (nkt(thermostat%al%loc_num_al))
1673  ALLOCATE (dof(thermostat%al%loc_num_al))
1674  DO i = 1, thermostat%al%loc_num_al
1675  nkt(i) = thermostat%al%nvt(i)%nkt
1676  dof(i) = real(thermostat%al%nvt(i)%degrees_of_freedom, kind=dp)
1677  END DO
1678  CALL get_temperatures(thermostat%al%map_info, thermostat%al%loc_num_al, &
1679  thermostat%al%glob_num_al, nkt, dof, para_env, tot_temperature, array_temp)
1680  DEALLOCATE (nkt)
1681  DEALLOCATE (dof)
1682  ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
1683  ! Energy associated with the GLE thermostat
1684  cpassert(ASSOCIATED(thermostat%gle))
1685 
1686  ALLOCATE (nkt(thermostat%gle%loc_num_gle))
1687  ALLOCATE (dof(thermostat%gle%loc_num_gle))
1688  DO i = 1, thermostat%gle%loc_num_gle
1689  nkt(i) = thermostat%gle%nvt(i)%nkt
1690  dof(i) = real(thermostat%gle%nvt(i)%degrees_of_freedom, kind=dp)
1691  END DO
1692  CALL get_temperatures(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
1693  thermostat%gle%glob_num_gle, nkt, dof, para_env, tot_temperature, array_temp)
1694  DEALLOCATE (nkt)
1695  DEALLOCATE (dof)
1696  END IF
1697  END IF
1698 
1699  END SUBROUTINE get_region_temperatures
1700 
1701 ! **************************************************************************************************
1702 !> \brief Prints status of all thermostats during an MD run
1703 !> \param thermostats ...
1704 !> \param para_env ...
1705 !> \param my_pos ...
1706 !> \param my_act ...
1707 !> \param itimes ...
1708 !> \param time ...
1709 !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1710 ! **************************************************************************************************
1711  SUBROUTINE print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
1712  TYPE(thermostats_type), POINTER :: thermostats
1713  TYPE(mp_para_env_type), POINTER :: para_env
1714  CHARACTER(LEN=default_string_length) :: my_pos, my_act
1715  INTEGER, INTENT(IN) :: itimes
1716  REAL(kind=dp), INTENT(IN) :: time
1717 
1718  IF (ASSOCIATED(thermostats)) THEN
1719  IF (ASSOCIATED(thermostats%thermostat_part)) THEN
1720  CALL print_thermostat_status(thermostats%thermostat_part, para_env, my_pos, my_act, itimes, time)
1721  END IF
1722  IF (ASSOCIATED(thermostats%thermostat_shell)) THEN
1723  CALL print_thermostat_status(thermostats%thermostat_shell, para_env, my_pos, my_act, itimes, time)
1724  END IF
1725  IF (ASSOCIATED(thermostats%thermostat_coef)) THEN
1726  CALL print_thermostat_status(thermostats%thermostat_coef, para_env, my_pos, my_act, itimes, time)
1727  END IF
1728  IF (ASSOCIATED(thermostats%thermostat_baro)) THEN
1729  CALL print_thermostat_status(thermostats%thermostat_baro, para_env, my_pos, my_act, itimes, time)
1730  END IF
1731  END IF
1732  END SUBROUTINE print_thermostats_status
1733 
1734 ! **************************************************************************************************
1735 !> \brief Prints status of a specific thermostat
1736 !> \param thermostat ...
1737 !> \param para_env ...
1738 !> \param my_pos ...
1739 !> \param my_act ...
1740 !> \param itimes ...
1741 !> \param time ...
1742 !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1743 ! **************************************************************************************************
1744  SUBROUTINE print_thermostat_status(thermostat, para_env, my_pos, my_act, itimes, time)
1745  TYPE(thermostat_type), POINTER :: thermostat
1746  TYPE(mp_para_env_type), POINTER :: para_env
1747  CHARACTER(LEN=default_string_length) :: my_pos, my_act
1748  INTEGER, INTENT(IN) :: itimes
1749  REAL(kind=dp), INTENT(IN) :: time
1750 
1751  INTEGER :: i, unit
1752  LOGICAL :: new_file
1753  REAL(kind=dp) :: thermo_kin, thermo_pot, tot_temperature
1754  REAL(kind=dp), DIMENSION(:), POINTER :: array_kin, array_pot, array_temp
1755  TYPE(cp_logger_type), POINTER :: logger
1756  TYPE(section_vals_type), POINTER :: print_key
1757 
1758  NULLIFY (logger, print_key, array_pot, array_kin, array_temp)
1759  logger => cp_get_default_logger()
1760 
1761  IF (ASSOCIATED(thermostat)) THEN
1762  ! Print Energies
1763  print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%ENERGY")
1764  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
1765  CALL get_thermostat_energies(thermostat, thermo_pot, thermo_kin, para_env, array_pot, array_kin)
1766  unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%ENERGY", &
1767  extension="."//trim(thermostat%label)//".tener", file_position=my_pos, &
1768  file_action=my_act, is_new_file=new_file)
1769  IF (unit > 0) THEN
1770  IF (new_file) THEN
1771  WRITE (unit, '(A)') "# Thermostat Potential and Kinetic Energies - Total and per Region"
1772  WRITE (unit, '("#",3X,A,2X,A,13X,A,10X,A)') "Step Nr.", "Time[fs]", "Kin.[a.u.]", "Pot.[a.u.]"
1773  END IF
1774  WRITE (unit=unit, fmt="(I8, F12.3,6X,2F20.10)") itimes, time*femtoseconds, thermo_kin, thermo_pot
1775  WRITE (unit, '(A,4F20.10)') "# KINETIC ENERGY REGIONS: ", array_kin(1:min(4, SIZE(array_kin)))
1776  DO i = 5, SIZE(array_kin), 4
1777  WRITE (unit=unit, fmt='("#",25X,4F20.10)') array_kin(i:min(i + 3, SIZE(array_kin)))
1778  END DO
1779  WRITE (unit, '(A,4F20.10)') "# POTENT. ENERGY REGIONS: ", array_pot(1:min(4, SIZE(array_pot)))
1780  DO i = 5, SIZE(array_pot), 4
1781  WRITE (unit=unit, fmt='("#",25X,4F20.10)') array_pot(i:min(i + 3, SIZE(array_pot)))
1782  END DO
1783  CALL m_flush(unit)
1784  END IF
1785  DEALLOCATE (array_kin)
1786  DEALLOCATE (array_pot)
1787  CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%ENERGY")
1788  END IF
1789  ! Print Temperatures of the regions
1790  print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%TEMPERATURE")
1791  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
1792  CALL get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
1793  unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%TEMPERATURE", &
1794  extension="."//trim(thermostat%label)//".temp", file_position=my_pos, &
1795  file_action=my_act, is_new_file=new_file)
1796  IF (unit > 0) THEN
1797  IF (new_file) THEN
1798  WRITE (unit, '(A)') "# Temperature Total and per Region"
1799  WRITE (unit, '("#",3X,A,2X,A,10X,A)') "Step Nr.", "Time[fs]", "Temp.[K]"
1800  END IF
1801  WRITE (unit=unit, fmt="(I8, F12.3,3X,F20.10)") itimes, time*femtoseconds, tot_temperature
1802  WRITE (unit, '(A,I10)') "# TEMPERATURE REGIONS: ", SIZE(array_temp)
1803  DO i = 1, SIZE(array_temp), 4
1804  WRITE (unit=unit, fmt='("#",22X,4F20.10)') array_temp(i:min(i + 3, SIZE(array_temp)))
1805  END DO
1806  CALL m_flush(unit)
1807  END IF
1808  DEALLOCATE (array_temp)
1809  CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%TEMPERATURE")
1810  END IF
1811  END IF
1812  END SUBROUTINE print_thermostat_status
1813 
1814 ! **************************************************************************************************
1815 !> \brief Handles the communication for thermostats (1D array)
1816 !> \param array ...
1817 !> \param number ...
1818 !> \param para_env ...
1819 !> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
1820 ! **************************************************************************************************
1821  SUBROUTINE communication_thermo_low1(array, number, para_env)
1822  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: array
1823  INTEGER, INTENT(IN) :: number
1824  TYPE(mp_para_env_type), POINTER :: para_env
1825 
1826  INTEGER :: i, icheck, ncheck
1827  REAL(kind=dp), DIMENSION(:), POINTER :: work, work2
1828 
1829  ALLOCATE (work(para_env%num_pe))
1830  DO i = 1, number
1831  work = 0.0_dp
1832  work(para_env%mepos + 1) = array(i)
1833  CALL para_env%sum(work)
1834  ncheck = count(work /= 0.0_dp)
1835  array(i) = 0.0_dp
1836  IF (ncheck /= 0) THEN
1837  ALLOCATE (work2(ncheck))
1838  ncheck = 0
1839  DO icheck = 1, para_env%num_pe
1840  IF (work(icheck) /= 0.0_dp) THEN
1841  ncheck = ncheck + 1
1842  work2(ncheck) = work(icheck)
1843  END IF
1844  END DO
1845  cpassert(ncheck == SIZE(work2))
1846  cpassert(all(work2 == work2(1)))
1847 
1848  array(i) = work2(1)
1849  DEALLOCATE (work2)
1850  END IF
1851  END DO
1852  DEALLOCATE (work)
1853  END SUBROUTINE communication_thermo_low1
1854 
1855 ! **************************************************************************************************
1856 !> \brief Handles the communication for thermostats (2D array)
1857 !> \param array ...
1858 !> \param number1 ...
1859 !> \param number2 ...
1860 !> \param para_env ...
1861 !> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
1862 ! **************************************************************************************************
1863  SUBROUTINE communication_thermo_low2(array, number1, number2, para_env)
1864  INTEGER, DIMENSION(:, :), INTENT(INOUT) :: array
1865  INTEGER, INTENT(IN) :: number1, number2
1866  TYPE(mp_para_env_type), POINTER :: para_env
1867 
1868  INTEGER :: i, icheck, j, ncheck
1869  INTEGER, DIMENSION(:, :), POINTER :: work, work2
1870 
1871  ALLOCATE (work(number1, para_env%num_pe))
1872  DO i = 1, number2
1873  work = 0
1874  work(:, para_env%mepos + 1) = array(:, i)
1875  CALL para_env%sum(work)
1876  ncheck = 0
1877  DO j = 1, para_env%num_pe
1878  IF (any(work(:, j) /= 0)) THEN
1879  ncheck = ncheck + 1
1880  END IF
1881  END DO
1882  array(:, i) = 0
1883  IF (ncheck /= 0) THEN
1884  ALLOCATE (work2(number1, ncheck))
1885  ncheck = 0
1886  DO icheck = 1, para_env%num_pe
1887  IF (any(work(:, icheck) /= 0)) THEN
1888  ncheck = ncheck + 1
1889  work2(:, ncheck) = work(:, icheck)
1890  END IF
1891  END DO
1892  cpassert(ncheck == SIZE(work2, 2))
1893  DO j = 1, ncheck
1894  cpassert(all(work2(:, j) == work2(:, 1)))
1895  END DO
1896  array(:, i) = work2(:, 1)
1897  DEALLOCATE (work2)
1898  END IF
1899  END DO
1900  DEALLOCATE (work)
1901  END SUBROUTINE communication_thermo_low2
1902 
1903 END MODULE thermostat_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
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_thermo_nose
integer, parameter, public do_thermo_nose
Definition: input_constants.F:110
input_constants::do_constr_atomic
integer, parameter, public do_constr_atomic
Definition: input_constants.F:79
input_constants::do_thermo_no_communication
integer, parameter, public do_thermo_no_communication
Definition: input_constants.F:296
input_constants::nvt_adiabatic_ensemble
integer, parameter, public nvt_adiabatic_ensemble
Definition: input_constants.F:95
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::langevin_ensemble
integer, parameter, public langevin_ensemble
Definition: input_constants.F:95
input_constants::do_region_molecule
integer, parameter, public do_region_molecule
Definition: input_constants.F:412
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition: input_constants.F:95
input_constants::do_thermo_al
integer, parameter, public do_thermo_al
Definition: input_constants.F:110
input_constants::do_constr_molec
integer, parameter, public do_constr_molec
Definition: input_constants.F:79
input_constants::do_thermo_csvr
integer, parameter, public do_thermo_csvr
Definition: input_constants.F:110
input_constants::do_thermo_gle
integer, parameter, public do_thermo_gle
Definition: input_constants.F:110
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::do_region_massive
integer, parameter, public do_region_massive
Definition: input_constants.F:412
input_constants::do_region_global
integer, parameter, public do_region_global
Definition: input_constants.F:412
input_constants::do_region_defined
integer, parameter, public do_region_defined
Definition: input_constants.F:412
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition: input_constants.F:95
input_constants::nvt_ensemble
integer, parameter, public nvt_ensemble
Definition: input_constants.F:95
input_constants::do_thermo_communication
integer, parameter, public do_thermo_communication
Definition: input_constants.F:296
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
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_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition: molecule_kind_types.F:510
molecule_kind_types::get_molecule_kind_set
subroutine, public get_molecule_kind_set(molecule_kind_set, maxatom, natom, nbond, nbend, nub, ntorsion, nimpr, nopbend, nconstraint, nconstraint_fixd, nmolecule, nrestraints)
Get informations about a molecule kind set.
Definition: molecule_kind_types.F:658
molecule_list_types
represent a simple array based list of the given type
Definition: molecule_list_types.F:15
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
motion_utils
Output Utilities for MOTION_SECTION.
Definition: motion_utils.F:13
motion_utils::rot_ana
subroutine, public rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weighted, natoms, rot_dof, inertia)
Performs an analysis of the principal inertia axis Getting back the generators of the translating and...
Definition: motion_utils.F:81
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::femtoseconds
real(kind=dp), parameter, public femtoseconds
Definition: physcon.F:153
qmmm_types
Basic container type for QM/MM.
Definition: qmmm_types.F:12
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_utils
Utilities for thermostats.
Definition: thermostat_utils.F:12
thermostat_utils::setup_adiabatic_thermostat_info
subroutine, public setup_adiabatic_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
...
Definition: thermostat_utils.F:234
thermostat_utils::compute_nfree
subroutine, public compute_nfree(cell, simpar, molecule_kind_set, print_section, particles, gci)
...
Definition: thermostat_utils.F:100
thermostat_utils::momentum_region_particles
subroutine, public momentum_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...
Definition: thermostat_utils.F:1046
thermostat_utils::vel_rescale_shells
subroutine, public vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
...
Definition: thermostat_utils.F:1295
thermostat_utils::communication_thermo_low2
subroutine, public communication_thermo_low2(array, number1, number2, para_env)
Handles the communication for thermostats (2D array)
Definition: thermostat_utils.F:1864
thermostat_utils::print_thermostats_status
subroutine, public print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
Prints status of all thermostats during an MD run.
Definition: thermostat_utils.F:1712
thermostat_utils::vel_rescale_particles
subroutine, public vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, local_molecules, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
Definition: thermostat_utils.F:1111
thermostat_utils::ke_region_shells
subroutine, public ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, group, core_particle_set, shell_particle_set, core_vel, shell_vel)
...
Definition: thermostat_utils.F:1216
thermostat_utils::get_thermostat_energies
subroutine, public get_thermostat_energies(thermostat, thermostat_pot, thermostat_kin, para_env, array_pot, array_kin)
Calculates energy associated with a thermostat.
Definition: thermostat_utils.F:1576
thermostat_utils::ke_region_baro
subroutine, public ke_region_baro(map_info, npt, group)
...
Definition: thermostat_utils.F:922
thermostat_utils::vel_rescale_baro
subroutine, public vel_rescale_baro(map_info, npt)
...
Definition: thermostat_utils.F:951
thermostat_utils::ke_region_particles
subroutine, public ke_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...
Definition: thermostat_utils.F:980
thermostat_utils::get_nhc_energies
subroutine, public get_nhc_energies(nhc, nhc_pot, nhc_kin, para_env, array_kin, array_pot)
Calculates kinetic energy and potential energy of the nhc variables.
Definition: thermostat_utils.F:1379
thermostat_utils::setup_thermostat_info
subroutine, public setup_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
...
Definition: thermostat_utils.F:551
thermostat_utils::compute_degrees_of_freedom
subroutine, public compute_degrees_of_freedom(thermostats, cell, simpar, molecule_kind_set, local_molecules, molecules, particles, print_section, region_sections, gci, region, qmmm_env)
...
Definition: thermostat_utils.F:152
thermostat_utils::get_kin_energies
subroutine, public get_kin_energies(map_info, loc_num, glob_num, thermo_energy, thermostat_kin, para_env, array_pot, array_kin)
Calculates kinetic energy and potential energy of the csvr and gle thermostats.
Definition: thermostat_utils.F:1448