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