(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
12MODULE 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
86CONTAINS
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
1903END MODULE thermostat_utils
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
extended_system_types::lnhc_parameters_type
Definition extended_system_types.F:75
extended_system_types::map_info_type
Definition extended_system_types.F:65
extended_system_types::npt_info_type
Definition extended_system_types.F:42
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::molecule_type
Definition molecule_types.F:110
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
qmmm_types::qmmm_env_type
Definition qmmm_types.F:34
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
thermostat_types::thermostat_info_type
Definition thermostat_types.F:90
thermostat_types::thermostat_type
Definition thermostat_types.F:78
thermostat_types::thermostats_type
Define thermostat types.
Definition thermostat_types.F:61