(git:495eafe)
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-2026 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_region_thermal, do_thermo_al, do_thermo_communication, &
32 do_thermo_csvr, do_thermo_gle, do_thermo_no_communication, do_thermo_nose, &
33 isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
34 nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
35 npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, 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 "A thermostat type DEFINED is requested but no thermostat "// &
654 "regions are defined in THERMOSTAT/DEFINE_REGION.")
655 CASE (do_region_thermal)
656 ! Similar to defined region above, but in THERMAL_REGION%DEFINE_REGION
657 nointer = .false.
658 ! Determine the number of thermostats defined in the input
659 CALL section_vals_get(region_sections, n_repetition=sum_of_thermostats)
660 IF (sum_of_thermostats < 1) &
661 CALL cp_abort(__location__, &
662 "A thermostat type THERMAL is requested but no thermal "// &
663 "regions are defined in THERMAL_REGION/DEFINE_REGION.")
664 END SELECT
665
666 ! Here we decide which parallel algorithm to use.
667 ! if there are only massive and molecule type thermostats we can use
668 ! a local scheme, in cases involving any combination with a
669 ! global thermostat we assume a coupling of degrees of freedom
670 ! from different processors
671 IF (nointer) THEN
672 ! Distributed thermostats, no interaction
673 dis_type = do_thermo_no_communication
674 ! we only count thermostats on this processor
675 number = 0
676 DO ikind = 1, nkind
677 nmol_local = local_molecules%n_el(ikind)
678 molecule_kind => molecule_kind_set(ikind)
679 CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
680 IF (do_shell) THEN
681 natom = nshell
682 IF (nshell == 0) nmol_local = 0
683 END IF
684 IF (region == do_region_molecule) THEN
685 number = number + nmol_local
686 ELSE IF (region == do_region_massive) THEN
687 number = number + 3*nmol_local*natom
688 ELSE
689 cpabort('Invalid region setup')
690 END IF
691 END DO
692 ELSE
693 ! REPlicated thermostats, INTERacting via communication
694 dis_type = do_thermo_communication
695 IF ((region == do_region_global) .OR. (region == do_region_molecule)) THEN
696 number = 1
697 ELSE IF ((region == do_region_defined) .OR. (region == do_region_thermal)) THEN
698 CALL get_defined_region_info(region_sections, number, sum_of_thermostats, &
699 map_loc_thermo_gen=thermostat_info%map_loc_thermo_gen, &
700 local_molecules=local_molecules, molecule_kind_set=molecule_kind_set, &
701 molecules=molecules, particles=particles, qmmm_env=qmmm_env)
702 END IF
703 END IF
704
705 IF (PRESENT(nfree)) THEN
706 IF ((sum_of_thermostats > 1) .OR. (dis_type == do_thermo_no_communication)) THEN
707 ! re-initializing simpar%nfree to zero because of multiple thermostats
708 nfree = 0
709 END IF
710 END IF
711 END IF
712 END SELECT
713
714 ! Saving information about thermostats
715 thermostat_info%sum_of_thermostats = sum_of_thermostats
716 thermostat_info%number_of_thermostats = number
717 thermostat_info%dis_type = dis_type
718 END SUBROUTINE setup_thermostat_info
719
720! **************************************************************************************************
721!> \brief ...
722!> \param region_sections ...
723!> \param number ...
724!> \param sum_of_thermostats ...
725!> \param map_loc_thermo_gen ...
726!> \param local_molecules ...
727!> \param molecule_kind_set ...
728!> \param molecules ...
729!> \param particles ...
730!> \param qmmm_env ...
731!> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
732! **************************************************************************************************
733 SUBROUTINE get_defined_region_info(region_sections, number, sum_of_thermostats, &
734 map_loc_thermo_gen, local_molecules, molecule_kind_set, molecules, particles, &
735 qmmm_env)
736 TYPE(section_vals_type), POINTER :: region_sections
737 INTEGER, INTENT(OUT), OPTIONAL :: number
738 INTEGER, INTENT(INOUT), OPTIONAL :: sum_of_thermostats
739 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
740 TYPE(distribution_1d_type), POINTER :: local_molecules
741 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
742 TYPE(molecule_list_type), POINTER :: molecules
743 TYPE(particle_list_type), POINTER :: particles
744 TYPE(qmmm_env_type), POINTER :: qmmm_env
745
746 CHARACTER(LEN=default_string_length), &
747 DIMENSION(:), POINTER :: tmpstringlist
748 INTEGER :: first_atom, i, ig, ikind, ilist, imol, ipart, jg, last_atom, mregions, n_rep, &
749 natom_local, nmol_per_kind, nregions, output_unit
750 INTEGER, DIMENSION(:), POINTER :: thermolist, tmp, tmplist
751 TYPE(cp_logger_type), POINTER :: logger
752 TYPE(molecule_kind_type), POINTER :: molecule_kind
753 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
754 TYPE(molecule_type), POINTER :: molecule
755
756 NULLIFY (tmplist, tmpstringlist, thermolist, molecule_kind, molecule, molecule_set)
757 NULLIFY (logger)
758 logger => cp_get_default_logger()
759 output_unit = cp_logger_get_default_io_unit(logger)
760 cpassert(.NOT. (ASSOCIATED(map_loc_thermo_gen)))
761 CALL section_vals_get(region_sections, n_repetition=nregions)
762 ALLOCATE (thermolist(particles%n_els))
763 thermolist = huge(0)
764 molecule_set => molecules%els
765 mregions = nregions
766 DO ig = 1, mregions
767 CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, n_rep_val=n_rep)
768 IF (n_rep > 0) THEN
769 DO jg = 1, n_rep
770 CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, i_rep_val=jg, i_vals=tmplist)
771 DO i = 1, SIZE(tmplist)
772 ipart = tmplist(i)
773 cpassert(((ipart > 0) .AND. (ipart <= particles%n_els)))
774 IF (thermolist(ipart) == huge(0) .OR. thermolist(ipart) == ig) THEN
775 thermolist(ipart) = ig
776 ELSE
777 CALL cp_abort(__location__, &
778 "The atom "//cp_to_string(ipart)//" has been "// &
779 "assigned to different thermostat regions "// &
780 cp_to_string(thermolist(ipart))//" and "// &
781 cp_to_string(ig)//" which is not allowed!")
782 END IF
783 END DO
784 END DO
785 END IF
786 CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, n_rep_val=n_rep)
787 IF (n_rep > 0) THEN
788 DO jg = 1, n_rep
789 CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, i_rep_val=jg, c_vals=tmpstringlist)
790 DO ilist = 1, SIZE(tmpstringlist)
791 DO ikind = 1, SIZE(molecule_kind_set)
792 molecule_kind => molecule_kind_set(ikind)
793 IF (molecule_kind%name == tmpstringlist(ilist)) THEN
794 DO imol = 1, SIZE(molecule_kind%molecule_list)
795 molecule => molecule_set(molecule_kind%molecule_list(imol))
796 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
797 DO ipart = first_atom, last_atom
798 IF (thermolist(ipart) == huge(0) .OR. thermolist(ipart) == ig) THEN
799 thermolist(ipart) = ig
800 ELSE
801 CALL cp_abort(__location__, &
802 "The atom "//cp_to_string(ipart)//" has been "// &
803 "assigned to different thermostat regions "// &
804 cp_to_string(thermolist(ipart))//" and "// &
805 cp_to_string(ig)//" which is not allowed!")
806 END IF
807 END DO
808 END DO
809 END IF
810 END DO
811 END DO
812 END DO
813 END IF
814 CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
815 subsys_qm=.false., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
816 CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
817 subsys_qm=.true., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
818 END DO
819
820 ! Dump IO warning for not thermalized particles
821 IF (any(thermolist == huge(0))) THEN
822 nregions = nregions + 1
823 sum_of_thermostats = sum_of_thermostats + 1
824 ALLOCATE (tmp(count(thermolist == huge(0))))
825 ilist = 0
826 DO i = 1, SIZE(thermolist)
827 IF (thermolist(i) == huge(0)) THEN
828 ilist = ilist + 1
829 tmp(ilist) = i
830 thermolist(i) = nregions
831 END IF
832 END DO
833 IF (ilist > 0) THEN
834 IF (output_unit > 0) THEN
835 WRITE (output_unit, '(/,T2,A)') &
836 "THERMOSTAT| Warning: No thermostats defined for the following atoms:"
837 DO i = 1, ilist, 8
838 WRITE (output_unit, '(T2,A,T17,8I8)') "THERMOSTAT|", tmp(i:min(i + 7, ilist))
839 END DO
840 WRITE (output_unit, '(T2,A)') &
841 "THERMOSTAT| They will be included in a further unique thermostat!"
842 END IF
843 END IF
844 DEALLOCATE (tmp)
845 END IF
846 cpassert(all(thermolist /= huge(0)))
847
848 ! Output thermostat region mapping to particles
849 ! The region indices are assumed to be 0-999
850 IF (output_unit > 0) THEN
851 WRITE (output_unit, '(/,T2,A)') &
852 "THERMOSTAT| Mapping of thermostat region indices to particles"
853 DO ipart = 1, particles%n_els, 16
854 WRITE (output_unit, '(T2,A,T17,16(" ",I3))') &
855 "THERMOSTAT|", thermolist(ipart:min(ipart + 15, particles%n_els))
856 END DO
857 END IF
858
859 ! Now identify the local number of thermostats
860 ALLOCATE (tmp(nregions))
861 tmp = 0
862 natom_local = 0
863 DO ikind = 1, SIZE(molecule_kind_set)
864 nmol_per_kind = local_molecules%n_el(ikind)
865 DO imol = 1, nmol_per_kind
866 i = local_molecules%list(ikind)%array(imol)
867 molecule => molecule_set(i)
868 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
869 DO ipart = first_atom, last_atom
870 natom_local = natom_local + 1
871 tmp(thermolist(ipart)) = 1
872 END DO
873 END DO
874 END DO
875 number = sum(tmp)
876 DEALLOCATE (tmp)
877
878 ! Now map the local atoms with the corresponding thermostat
879 ALLOCATE (map_loc_thermo_gen(natom_local))
880 natom_local = 0
881 DO ikind = 1, SIZE(molecule_kind_set)
882 nmol_per_kind = local_molecules%n_el(ikind)
883 DO imol = 1, nmol_per_kind
884 i = local_molecules%list(ikind)%array(imol)
885 molecule => molecule_set(i)
886 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
887 DO ipart = first_atom, last_atom
888 natom_local = natom_local + 1
889 map_loc_thermo_gen(natom_local) = thermolist(ipart)
890 END DO
891 END DO
892 END DO
893
894 DEALLOCATE (thermolist)
895 END SUBROUTINE get_defined_region_info
896
897! **************************************************************************************************
898!> \brief ...
899!> \param region_sections ...
900!> \param qmmm_env ...
901!> \param thermolist ...
902!> \param molecule_set ...
903!> \param subsys_qm ...
904!> \param ig ...
905!> \param sum_of_thermostats ...
906!> \param nregions ...
907!> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
908! **************************************************************************************************
909 SUBROUTINE setup_thermostat_subsys(region_sections, qmmm_env, thermolist, &
910 molecule_set, subsys_qm, ig, sum_of_thermostats, nregions)
911 TYPE(section_vals_type), POINTER :: region_sections
912 TYPE(qmmm_env_type), POINTER :: qmmm_env
913 INTEGER, DIMENSION(:), POINTER :: thermolist
914 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
915 LOGICAL, INTENT(IN) :: subsys_qm
916 INTEGER, INTENT(IN) :: ig
917 INTEGER, INTENT(INOUT) :: sum_of_thermostats, nregions
918
919 CHARACTER(LEN=default_string_length) :: label1, label2
920 INTEGER :: first_atom, i, imolecule, ipart, &
921 last_atom, nrep, thermo1
922 INTEGER, DIMENSION(:), POINTER :: atom_index1
923 LOGICAL :: explicit
924 TYPE(molecule_type), POINTER :: molecule
925
926 label1 = "MM_SUBSYS"
927 label2 = "QM_SUBSYS"
928 IF (subsys_qm) THEN
929 label1 = "QM_SUBSYS"
930 label2 = "MM_SUBSYS"
931 END IF
932 CALL section_vals_val_get(region_sections, trim(label1), i_rep_section=ig, &
933 n_rep_val=nrep, explicit=explicit)
934 IF (nrep == 1 .AND. explicit) THEN
935 IF (ASSOCIATED(qmmm_env)) THEN
936 atom_index1 => qmmm_env%qm%mm_atom_index
937 IF (subsys_qm) THEN
938 atom_index1 => qmmm_env%qm%qm_atom_index
939 END IF
940 CALL section_vals_val_get(region_sections, trim(label1), i_val=thermo1, i_rep_section=ig)
941 SELECT CASE (thermo1)
942 CASE (do_constr_atomic)
943 DO i = 1, SIZE(atom_index1)
944 ipart = atom_index1(i)
945 IF (subsys_qm .AND. qmmm_env%qm%qmmm_link .AND. ASSOCIATED(qmmm_env%qm%mm_link_atoms)) THEN
946 IF (any(ipart == qmmm_env%qm%mm_link_atoms)) cycle
947 END IF
948 IF (thermolist(ipart) == huge(0)) THEN
949 thermolist(ipart) = ig
950 ELSE
951 CALL cp_abort(__location__, &
952 'One atom ('//cp_to_string(ipart)//') of the '// &
953 trim(label1)//' was already assigned to'// &
954 ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
955 '. Please check the input for inconsistencies!')
956 END IF
957 END DO
958 CASE (do_constr_molec)
959 DO imolecule = 1, SIZE(molecule_set)
960 molecule => molecule_set(imolecule)
961 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
962 IF (any(atom_index1 >= first_atom .AND. atom_index1 <= last_atom)) THEN
963 DO ipart = first_atom, last_atom
964 IF (thermolist(ipart) == huge(0)) THEN
965 thermolist(ipart) = ig
966 ELSE
967 CALL cp_abort(__location__, &
968 'One atom ('//cp_to_string(ipart)//') of the '// &
969 trim(label1)//' was already assigned to'// &
970 ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
971 '. Please check the input for inconsistencies!')
972 END IF
973 END DO
974 END IF
975 END DO
976 END SELECT
977 ELSE
978 sum_of_thermostats = sum_of_thermostats - 1
979 nregions = nregions - 1
980 END IF
981 END IF
982 END SUBROUTINE setup_thermostat_subsys
983
984! **************************************************************************************************
985!> \brief ...
986!> \param map_info ...
987!> \param npt ...
988!> \param group ...
989!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
990! **************************************************************************************************
991 SUBROUTINE ke_region_baro(map_info, npt, group)
992 TYPE(map_info_type), POINTER :: map_info
993 TYPE(npt_info_type), DIMENSION(:, :), &
994 INTENT(INOUT) :: npt
995 TYPE(mp_comm_type), INTENT(IN) :: group
996
997 INTEGER :: i, j, ncoef
998
999 map_info%v_scale = 1.0_dp
1000 map_info%s_kin = 0.0_dp
1001 ncoef = 0
1002 DO i = 1, SIZE(npt, 1)
1003 DO j = 1, SIZE(npt, 2)
1004 ncoef = ncoef + 1
1005 map_info%p_kin(1, ncoef)%point = map_info%p_kin(1, ncoef)%point &
1006 + npt(i, j)%mass*npt(i, j)%v**2
1007 END DO
1008 END DO
1009
1010 IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1011
1012 END SUBROUTINE ke_region_baro
1013
1014! **************************************************************************************************
1015!> \brief ...
1016!> \param map_info ...
1017!> \param npt ...
1018!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1019! **************************************************************************************************
1020 SUBROUTINE vel_rescale_baro(map_info, npt)
1021 TYPE(map_info_type), POINTER :: map_info
1022 TYPE(npt_info_type), DIMENSION(:, :), &
1023 INTENT(INOUT) :: npt
1024
1025 INTEGER :: i, j, ncoef
1026
1027 ncoef = 0
1028 DO i = 1, SIZE(npt, 1)
1029 DO j = 1, SIZE(npt, 2)
1030 ncoef = ncoef + 1
1031 npt(i, j)%v = npt(i, j)%v*map_info%p_scale(1, ncoef)%point
1032 END DO
1033 END DO
1034
1035 END SUBROUTINE vel_rescale_baro
1036
1037! **************************************************************************************************
1038!> \brief ...
1039!> \param map_info ...
1040!> \param particle_set ...
1041!> \param molecule_kind_set ...
1042!> \param local_molecules ...
1043!> \param molecule_set ...
1044!> \param group ...
1045!> \param vel ...
1046!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1047! **************************************************************************************************
1048 SUBROUTINE ke_region_particles(map_info, particle_set, molecule_kind_set, &
1049 local_molecules, molecule_set, group, vel)
1050
1051 TYPE(map_info_type), POINTER :: map_info
1052 TYPE(particle_type), POINTER :: particle_set(:)
1053 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1054 TYPE(distribution_1d_type), POINTER :: local_molecules
1055 TYPE(molecule_type), POINTER :: molecule_set(:)
1056 TYPE(mp_comm_type), INTENT(IN) :: group
1057 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
1058
1059 INTEGER :: first_atom, ii, ikind, imol, imol_local, &
1060 ipart, last_atom, nmol_local
1061 LOGICAL :: present_vel
1062 REAL(kind=dp) :: mass
1063 TYPE(atomic_kind_type), POINTER :: atomic_kind
1064 TYPE(molecule_type), POINTER :: molecule
1065
1066 map_info%v_scale = 1.0_dp
1067 map_info%s_kin = 0.0_dp
1068 present_vel = PRESENT(vel)
1069 ii = 0
1070 DO ikind = 1, SIZE(molecule_kind_set)
1071 nmol_local = local_molecules%n_el(ikind)
1072 DO imol_local = 1, nmol_local
1073 imol = local_molecules%list(ikind)%array(imol_local)
1074 molecule => molecule_set(imol)
1075 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1076 DO ipart = first_atom, last_atom
1077 ii = ii + 1
1078 atomic_kind => particle_set(ipart)%atomic_kind
1079 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1080 IF (present_vel) THEN
1081 IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
1082 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*vel(1, ipart)**2
1083 IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
1084 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*vel(2, ipart)**2
1085 IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
1086 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*vel(3, ipart)**2
1087 ELSE
1088 IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
1089 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*particle_set(ipart)%v(1)**2
1090 IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
1091 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*particle_set(ipart)%v(2)**2
1092 IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
1093 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*particle_set(ipart)%v(3)**2
1094 END IF
1095 END DO
1096 END DO
1097 END DO
1098
1099 IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1100
1101 END SUBROUTINE ke_region_particles
1102
1103! **************************************************************************************************
1104!> \brief ...
1105!> \param map_info ...
1106!> \param particle_set ...
1107!> \param molecule_kind_set ...
1108!> \param local_molecules ...
1109!> \param molecule_set ...
1110!> \param group ...
1111!> \param vel ...
1112!> \author 07.2009 MI
1113! **************************************************************************************************
1114 SUBROUTINE momentum_region_particles(map_info, particle_set, molecule_kind_set, &
1115 local_molecules, molecule_set, group, vel)
1116
1117 TYPE(map_info_type), POINTER :: map_info
1118 TYPE(particle_type), POINTER :: particle_set(:)
1119 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1120 TYPE(distribution_1d_type), POINTER :: local_molecules
1121 TYPE(molecule_type), POINTER :: molecule_set(:)
1122 TYPE(mp_comm_type), INTENT(IN) :: group
1123 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
1124
1125 INTEGER :: first_atom, ii, ikind, imol, imol_local, &
1126 ipart, last_atom, nmol_local
1127 LOGICAL :: present_vel
1128 REAL(kind=dp) :: mass
1129 TYPE(atomic_kind_type), POINTER :: atomic_kind
1130 TYPE(molecule_type), POINTER :: molecule
1131
1132 map_info%v_scale = 1.0_dp
1133 map_info%s_kin = 0.0_dp
1134 present_vel = PRESENT(vel)
1135 ii = 0
1136 DO ikind = 1, SIZE(molecule_kind_set)
1137 nmol_local = local_molecules%n_el(ikind)
1138 DO imol_local = 1, nmol_local
1139 imol = local_molecules%list(ikind)%array(imol_local)
1140 molecule => molecule_set(imol)
1141 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1142 DO ipart = first_atom, last_atom
1143 ii = ii + 1
1144 atomic_kind => particle_set(ipart)%atomic_kind
1145 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1146 IF (present_vel) THEN
1147 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + sqrt(mass)*vel(1, ipart)
1148 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + sqrt(mass)*vel(2, ipart)
1149 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + sqrt(mass)*vel(3, ipart)
1150 ELSE
1151 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + sqrt(mass)*particle_set(ipart)%v(1)
1152 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + sqrt(mass)*particle_set(ipart)%v(2)
1153 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + sqrt(mass)*particle_set(ipart)%v(3)
1154 END IF
1155 END DO
1156 END DO
1157 END DO
1158
1159 IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1160
1161 END SUBROUTINE momentum_region_particles
1162
1163! **************************************************************************************************
1164!> \brief ...
1165!> \param map_info ...
1166!> \param molecule_kind_set ...
1167!> \param molecule_set ...
1168!> \param particle_set ...
1169!> \param local_molecules ...
1170!> \param shell_adiabatic ...
1171!> \param shell_particle_set ...
1172!> \param core_particle_set ...
1173!> \param vel ...
1174!> \param shell_vel ...
1175!> \param core_vel ...
1176!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1177! **************************************************************************************************
1178 SUBROUTINE vel_rescale_particles(map_info, molecule_kind_set, molecule_set, &
1179 particle_set, local_molecules, shell_adiabatic, shell_particle_set, &
1180 core_particle_set, vel, shell_vel, core_vel)
1181
1182 TYPE(map_info_type), POINTER :: map_info
1183 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1184 TYPE(molecule_type), POINTER :: molecule_set(:)
1185 TYPE(particle_type), POINTER :: particle_set(:)
1186 TYPE(distribution_1d_type), POINTER :: local_molecules
1187 LOGICAL, INTENT(IN) :: shell_adiabatic
1188 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
1189 core_particle_set(:)
1190 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
1191 core_vel(:, :)
1192
1193 INTEGER :: first_atom, ii, ikind, imol, imol_local, &
1194 ipart, jj, last_atom, nmol_local, &
1195 shell_index
1196 LOGICAL :: present_vel
1197 REAL(kind=dp) :: fac_massc, fac_masss, mass, vc(3), vs(3)
1198 TYPE(atomic_kind_type), POINTER :: atomic_kind
1199 TYPE(molecule_type), POINTER :: molecule
1200 TYPE(shell_kind_type), POINTER :: shell
1201
1202 ii = 0
1203 jj = 0
1204 present_vel = PRESENT(vel)
1205 ! Just few checks for consistency
1206 IF (present_vel) THEN
1207 IF (shell_adiabatic) THEN
1208 cpassert(PRESENT(shell_vel))
1209 cpassert(PRESENT(core_vel))
1210 END IF
1211 ELSE
1212 IF (shell_adiabatic) THEN
1213 cpassert(PRESENT(shell_particle_set))
1214 cpassert(PRESENT(core_particle_set))
1215 END IF
1216 END IF
1217 kind: DO ikind = 1, SIZE(molecule_kind_set)
1218 nmol_local = local_molecules%n_el(ikind)
1219 mol_local: DO imol_local = 1, nmol_local
1220 imol = local_molecules%list(ikind)%array(imol_local)
1221 molecule => molecule_set(imol)
1222 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1223 particle: DO ipart = first_atom, last_atom
1224 ii = ii + 1
1225 IF (present_vel) THEN
1226 vel(1, ipart) = vel(1, ipart)*map_info%p_scale(1, ii)%point
1227 vel(2, ipart) = vel(2, ipart)*map_info%p_scale(2, ii)%point
1228 vel(3, ipart) = vel(3, ipart)*map_info%p_scale(3, ii)%point
1229 ELSE
1230 particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*map_info%p_scale(1, ii)%point
1231 particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*map_info%p_scale(2, ii)%point
1232 particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*map_info%p_scale(3, ii)%point
1233 END IF
1234 ! If Shell Adiabatic then apply the NHC thermostat also to the Shells
1235 IF (shell_adiabatic) THEN
1236 shell_index = particle_set(ipart)%shell_index
1237 IF (shell_index /= 0) THEN
1238 jj = jj + 2
1239 atomic_kind => particle_set(ipart)%atomic_kind
1240 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
1241 fac_masss = shell%mass_shell/mass
1242 fac_massc = shell%mass_core/mass
1243 IF (present_vel) THEN
1244 vs(1:3) = shell_vel(1:3, shell_index)
1245 vc(1:3) = core_vel(1:3, shell_index)
1246 shell_vel(1, shell_index) = vel(1, ipart) + fac_massc*(vs(1) - vc(1))
1247 shell_vel(2, shell_index) = vel(2, ipart) + fac_massc*(vs(2) - vc(2))
1248 shell_vel(3, shell_index) = vel(3, ipart) + fac_massc*(vs(3) - vc(3))
1249 core_vel(1, shell_index) = vel(1, ipart) + fac_masss*(vc(1) - vs(1))
1250 core_vel(2, shell_index) = vel(2, ipart) + fac_masss*(vc(2) - vs(2))
1251 core_vel(3, shell_index) = vel(3, ipart) + fac_masss*(vc(3) - vs(3))
1252 ELSE
1253 vs(1:3) = shell_particle_set(shell_index)%v(1:3)
1254 vc(1:3) = core_particle_set(shell_index)%v(1:3)
1255 shell_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_massc*(vs(1) - vc(1))
1256 shell_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_massc*(vs(2) - vc(2))
1257 shell_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_massc*(vs(3) - vc(3))
1258 core_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_masss*(vc(1) - vs(1))
1259 core_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_masss*(vc(2) - vs(2))
1260 core_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_masss*(vc(3) - vs(3))
1261 END IF
1262 END IF
1263 END IF
1264 END DO particle
1265 END DO mol_local
1266 END DO kind
1267
1268 END SUBROUTINE vel_rescale_particles
1269
1270! **************************************************************************************************
1271!> \brief ...
1272!> \param map_info ...
1273!> \param particle_set ...
1274!> \param atomic_kind_set ...
1275!> \param local_particles ...
1276!> \param group ...
1277!> \param core_particle_set ...
1278!> \param shell_particle_set ...
1279!> \param core_vel ...
1280!> \param shell_vel ...
1281!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1282! **************************************************************************************************
1283 SUBROUTINE ke_region_shells(map_info, particle_set, atomic_kind_set, &
1284 local_particles, group, core_particle_set, shell_particle_set, &
1285 core_vel, shell_vel)
1286
1287 TYPE(map_info_type), POINTER :: map_info
1288 TYPE(particle_type), POINTER :: particle_set(:)
1289 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
1290 TYPE(distribution_1d_type), POINTER :: local_particles
1291 TYPE(mp_comm_type), INTENT(IN) :: group
1292 TYPE(particle_type), OPTIONAL, POINTER :: core_particle_set(:), &
1293 shell_particle_set(:)
1294 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: core_vel(:, :), shell_vel(:, :)
1295
1296 INTEGER :: ii, iparticle, iparticle_kind, &
1297 iparticle_local, nparticle_kind, &
1298 nparticle_local, shell_index
1299 LOGICAL :: is_shell, present_vel
1300 REAL(dp) :: mass, mu_mass, v_sc(3)
1301 TYPE(atomic_kind_type), POINTER :: atomic_kind
1302 TYPE(shell_kind_type), POINTER :: shell
1303
1304 present_vel = PRESENT(shell_vel)
1305 ! Preliminary checks for consistency usage
1306 IF (present_vel) THEN
1307 cpassert(PRESENT(core_vel))
1308 ELSE
1309 cpassert(PRESENT(shell_particle_set))
1310 cpassert(PRESENT(core_particle_set))
1311 END IF
1312 ! get force on first thermostat for all the chains in the system.
1313 map_info%v_scale = 1.0_dp
1314 map_info%s_kin = 0.0_dp
1315 ii = 0
1316
1317 nparticle_kind = SIZE(atomic_kind_set)
1318 DO iparticle_kind = 1, nparticle_kind
1319 atomic_kind => atomic_kind_set(iparticle_kind)
1320 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
1321 IF (is_shell) THEN
1322 mu_mass = shell%mass_shell*shell%mass_core/mass
1323 nparticle_local = local_particles%n_el(iparticle_kind)
1324 DO iparticle_local = 1, nparticle_local
1325 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1326 shell_index = particle_set(iparticle)%shell_index
1327 ii = ii + 1
1328 IF (present_vel) THEN
1329 v_sc(1) = core_vel(1, shell_index) - shell_vel(1, shell_index)
1330 v_sc(2) = core_vel(2, shell_index) - shell_vel(2, shell_index)
1331 v_sc(3) = core_vel(3, shell_index) - shell_vel(3, shell_index)
1332 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
1333 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
1334 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
1335 ELSE
1336 v_sc(1) = core_particle_set(shell_index)%v(1) - shell_particle_set(shell_index)%v(1)
1337 v_sc(2) = core_particle_set(shell_index)%v(2) - shell_particle_set(shell_index)%v(2)
1338 v_sc(3) = core_particle_set(shell_index)%v(3) - shell_particle_set(shell_index)%v(3)
1339 map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
1340 map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
1341 map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
1342 END IF
1343 END DO
1344 END IF
1345 END DO
1346 IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
1347
1348 END SUBROUTINE ke_region_shells
1349
1350! **************************************************************************************************
1351!> \brief ...
1352!> \param map_info ...
1353!> \param atomic_kind_set ...
1354!> \param particle_set ...
1355!> \param local_particles ...
1356!> \param shell_particle_set ...
1357!> \param core_particle_set ...
1358!> \param shell_vel ...
1359!> \param core_vel ...
1360!> \param vel ...
1361!> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
1362! **************************************************************************************************
1363 SUBROUTINE vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
1364 shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
1365
1366 TYPE(map_info_type), POINTER :: map_info
1367 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
1368 TYPE(particle_type), POINTER :: particle_set(:)
1369 TYPE(distribution_1d_type), POINTER :: local_particles
1370 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
1371 core_particle_set(:)
1372 REAL(kind=dp), INTENT(INOUT), OPTIONAL :: shell_vel(:, :), core_vel(:, :), &
1373 vel(:, :)
1374
1375 INTEGER :: ii, iparticle, iparticle_kind, &
1376 iparticle_local, nparticle_kind, &
1377 nparticle_local, shell_index
1378 LOGICAL :: is_shell, present_vel
1379 REAL(dp) :: mass, massc, masss, umass, v(3), vc(3), &
1380 vs(3)
1381 TYPE(atomic_kind_type), POINTER :: atomic_kind
1382 TYPE(shell_kind_type), POINTER :: shell
1383
1384 present_vel = PRESENT(vel)
1385 ! Preliminary checks for consistency usage
1386 IF (present_vel) THEN
1387 cpassert(PRESENT(shell_vel))
1388 cpassert(PRESENT(core_vel))
1389 ELSE
1390 cpassert(PRESENT(shell_particle_set))
1391 cpassert(PRESENT(core_particle_set))
1392 END IF
1393 ii = 0
1394 nparticle_kind = SIZE(atomic_kind_set)
1395 ! now scale the core-shell velocities
1396 kind: DO iparticle_kind = 1, nparticle_kind
1397 atomic_kind => atomic_kind_set(iparticle_kind)
1398 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
1399 IF (is_shell) THEN
1400 umass = 1.0_dp/mass
1401 masss = shell%mass_shell*umass
1402 massc = shell%mass_core*umass
1403
1404 nparticle_local = local_particles%n_el(iparticle_kind)
1405 particles: DO iparticle_local = 1, nparticle_local
1406 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1407 shell_index = particle_set(iparticle)%shell_index
1408 ii = ii + 1
1409 IF (present_vel) THEN
1410 vc(1:3) = core_vel(1:3, shell_index)
1411 vs(1:3) = shell_vel(1:3, shell_index)
1412 v(1:3) = vel(1:3, iparticle)
1413 shell_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
1414 shell_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
1415 shell_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
1416 core_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
1417 core_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
1418 core_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
1419 ELSE
1420 vc(1:3) = core_particle_set(shell_index)%v(1:3)
1421 vs(1:3) = shell_particle_set(shell_index)%v(1:3)
1422 v(1:3) = particle_set(iparticle)%v(1:3)
1423 shell_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
1424 shell_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
1425 shell_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
1426 core_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
1427 core_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
1428 core_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
1429 END IF
1430 END DO particles
1431 END IF
1432 END DO kind
1433
1434 END SUBROUTINE vel_rescale_shells
1435
1436! **************************************************************************************************
1437!> \brief Calculates kinetic energy and potential energy of the nhc variables
1438!> \param nhc ...
1439!> \param nhc_pot ...
1440!> \param nhc_kin ...
1441!> \param para_env ...
1442!> \param array_kin ...
1443!> \param array_pot ...
1444!> \par History
1445!> none
1446!> \author CJM
1447! **************************************************************************************************
1448 SUBROUTINE get_nhc_energies(nhc, nhc_pot, nhc_kin, para_env, array_kin, array_pot)
1449 TYPE(lnhc_parameters_type), POINTER :: nhc
1450 REAL(kind=dp), INTENT(OUT) :: nhc_pot, nhc_kin
1451 TYPE(mp_para_env_type), POINTER :: para_env
1452 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_kin, array_pot
1453
1454 INTEGER :: imap, l, n, number
1455 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin, vpot
1456
1457 number = nhc%glob_num_nhc
1458 ALLOCATE (akin(number))
1459 ALLOCATE (vpot(number))
1460 akin = 0.0_dp
1461 vpot = 0.0_dp
1462 DO n = 1, nhc%loc_num_nhc
1463 imap = nhc%map_info%index(n)
1464 DO l = 1, nhc%nhc_len
1465 akin(imap) = akin(imap) + 0.5_dp*nhc%nvt(l, n)%mass*nhc%nvt(l, n)%v**2
1466 vpot(imap) = vpot(imap) + nhc%nvt(l, n)%nkt*nhc%nvt(l, n)%eta
1467 END DO
1468 END DO
1469
1470 ! Handle the thermostat distribution
1471 IF (nhc%map_info%dis_type == do_thermo_no_communication) THEN
1472 CALL para_env%sum(akin)
1473 CALL para_env%sum(vpot)
1474 ELSE IF (nhc%map_info%dis_type == do_thermo_communication) THEN
1475 CALL communication_thermo_low1(akin, number, para_env)
1476 CALL communication_thermo_low1(vpot, number, para_env)
1477 END IF
1478 nhc_kin = sum(akin)
1479 nhc_pot = sum(vpot)
1480
1481 ! Possibly give back kinetic or potential energy arrays
1482 IF (PRESENT(array_pot)) THEN
1483 IF (ASSOCIATED(array_pot)) THEN
1484 cpassert(SIZE(array_pot) == number)
1485 ELSE
1486 ALLOCATE (array_pot(number))
1487 END IF
1488 array_pot = vpot
1489 END IF
1490 IF (PRESENT(array_kin)) THEN
1491 IF (ASSOCIATED(array_kin)) THEN
1492 cpassert(SIZE(array_kin) == number)
1493 ELSE
1494 ALLOCATE (array_kin(number))
1495 END IF
1496 array_kin = akin
1497 END IF
1498 DEALLOCATE (akin)
1499 DEALLOCATE (vpot)
1500 END SUBROUTINE get_nhc_energies
1501
1502! **************************************************************************************************
1503!> \brief Calculates kinetic energy and potential energy
1504!> of the csvr and gle thermostats
1505!> \param map_info ...
1506!> \param loc_num ...
1507!> \param glob_num ...
1508!> \param thermo_energy ...
1509!> \param thermostat_kin ...
1510!> \param para_env ...
1511!> \param array_pot ...
1512!> \param array_kin ...
1513!> \par History generalized MI [07.2009]
1514!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1515! **************************************************************************************************
1516 SUBROUTINE get_kin_energies(map_info, loc_num, glob_num, thermo_energy, thermostat_kin, &
1517 para_env, array_pot, array_kin)
1518
1519 TYPE(map_info_type), POINTER :: map_info
1520 INTEGER, INTENT(IN) :: loc_num, glob_num
1521 REAL(dp), DIMENSION(:), INTENT(IN) :: thermo_energy
1522 REAL(kind=dp), INTENT(OUT) :: thermostat_kin
1523 TYPE(mp_para_env_type), POINTER :: para_env
1524 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_pot, array_kin
1525
1526 INTEGER :: imap, n, number
1527 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin
1528
1529 number = glob_num
1530 ALLOCATE (akin(number))
1531 akin = 0.0_dp
1532 DO n = 1, loc_num
1533 imap = map_info%index(n)
1534 akin(imap) = thermo_energy(n)
1535 END DO
1536
1537 ! Handle the thermostat distribution
1538 IF (map_info%dis_type == do_thermo_no_communication) THEN
1539 CALL para_env%sum(akin)
1540 ELSE IF (map_info%dis_type == do_thermo_communication) THEN
1541 CALL communication_thermo_low1(akin, number, para_env)
1542 END IF
1543 thermostat_kin = sum(akin)
1544
1545 ! Possibly give back kinetic or potential energy arrays
1546 IF (PRESENT(array_pot)) THEN
1547 IF (ASSOCIATED(array_pot)) THEN
1548 cpassert(SIZE(array_pot) == number)
1549 ELSE
1550 ALLOCATE (array_pot(number))
1551 END IF
1552 array_pot = 0.0_dp
1553 END IF
1554 IF (PRESENT(array_kin)) THEN
1555 IF (ASSOCIATED(array_kin)) THEN
1556 cpassert(SIZE(array_kin) == number)
1557 ELSE
1558 ALLOCATE (array_kin(number))
1559 END IF
1560 array_kin = akin
1561 END IF
1562 DEALLOCATE (akin)
1563 END SUBROUTINE get_kin_energies
1564
1565! **************************************************************************************************
1566!> \brief Calculates the temperatures of the regions when a thermostat is
1567!> applied
1568!> \param map_info ...
1569!> \param loc_num ...
1570!> \param glob_num ...
1571!> \param nkt ...
1572!> \param dof ...
1573!> \param para_env ...
1574!> \param temp_tot ...
1575!> \param array_temp ...
1576!> \par History generalized MI [07.2009]
1577!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1578! **************************************************************************************************
1579 SUBROUTINE get_temperatures(map_info, loc_num, glob_num, nkt, dof, para_env, &
1580 temp_tot, array_temp)
1581 TYPE(map_info_type), POINTER :: map_info
1582 INTEGER, INTENT(IN) :: loc_num, glob_num
1583 REAL(dp), DIMENSION(:), INTENT(IN) :: nkt, dof
1584 TYPE(mp_para_env_type), POINTER :: para_env
1585 REAL(kind=dp), INTENT(OUT) :: temp_tot
1586 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_temp
1587
1588 INTEGER :: i, imap, imap2, n, number
1589 REAL(kind=dp) :: fdeg_of_free
1590 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: akin, deg_of_free
1591
1592 number = glob_num
1593 ALLOCATE (akin(number))
1594 ALLOCATE (deg_of_free(number))
1595 akin = 0.0_dp
1596 deg_of_free = 0.0_dp
1597 DO n = 1, loc_num
1598 imap = map_info%index(n)
1599 imap2 = map_info%map_index(n)
1600 IF (nkt(n) == 0.0_dp) cycle
1601 deg_of_free(imap) = real(dof(n), kind=dp)
1602 akin(imap) = map_info%s_kin(imap2)
1603 END DO
1604
1605 ! Handle the thermostat distribution
1606 IF (map_info%dis_type == do_thermo_no_communication) THEN
1607 CALL para_env%sum(akin)
1608 CALL para_env%sum(deg_of_free)
1609 ELSE IF (map_info%dis_type == do_thermo_communication) THEN
1610 CALL communication_thermo_low1(akin, number, para_env)
1611 CALL communication_thermo_low1(deg_of_free, number, para_env)
1612 END IF
1613 temp_tot = sum(akin)
1614 fdeg_of_free = sum(deg_of_free)
1615
1616 temp_tot = temp_tot/fdeg_of_free
1617 temp_tot = cp_unit_from_cp2k(temp_tot, "K_temp")
1618 ! Possibly give back temperatures of the full set of regions
1619 IF (PRESENT(array_temp)) THEN
1620 IF (ASSOCIATED(array_temp)) THEN
1621 cpassert(SIZE(array_temp) == number)
1622 ELSE
1623 ALLOCATE (array_temp(number))
1624 END IF
1625 DO i = 1, number
1626 array_temp(i) = akin(i)/deg_of_free(i)
1627 array_temp(i) = cp_unit_from_cp2k(array_temp(i), "K_temp")
1628 END DO
1629 END IF
1630 DEALLOCATE (akin)
1631 DEALLOCATE (deg_of_free)
1632 END SUBROUTINE get_temperatures
1633
1634! **************************************************************************************************
1635!> \brief Calculates energy associated with a thermostat
1636!> \param thermostat ...
1637!> \param thermostat_pot ...
1638!> \param thermostat_kin ...
1639!> \param para_env ...
1640!> \param array_pot ...
1641!> \param array_kin ...
1642!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1643! **************************************************************************************************
1644 SUBROUTINE get_thermostat_energies(thermostat, thermostat_pot, thermostat_kin, para_env, &
1645 array_pot, array_kin)
1646 TYPE(thermostat_type), POINTER :: thermostat
1647 REAL(kind=dp), INTENT(OUT) :: thermostat_pot, thermostat_kin
1648 TYPE(mp_para_env_type), POINTER :: para_env
1649 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_pot, array_kin
1650
1651 INTEGER :: i
1652 REAL(dp), ALLOCATABLE, DIMENSION(:) :: thermo_energy
1653
1654 thermostat_pot = 0.0_dp
1655 thermostat_kin = 0.0_dp
1656 IF (ASSOCIATED(thermostat)) THEN
1657 IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
1658 ! Energy associated with the Nose-Hoover thermostat
1659 cpassert(ASSOCIATED(thermostat%nhc))
1660 CALL get_nhc_energies(thermostat%nhc, thermostat_pot, thermostat_kin, para_env, &
1661 array_pot, array_kin)
1662 ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
1663 ! Energy associated with the CSVR thermostat
1664 cpassert(ASSOCIATED(thermostat%csvr))
1665 ALLOCATE (thermo_energy(thermostat%csvr%loc_num_csvr))
1666 DO i = 1, thermostat%csvr%loc_num_csvr
1667 thermo_energy(i) = thermostat%csvr%nvt(i)%thermostat_energy
1668 END DO
1669 CALL get_kin_energies(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
1670 thermostat%csvr%glob_num_csvr, thermo_energy, &
1671 thermostat_kin, para_env, array_pot, array_kin)
1672 DEALLOCATE (thermo_energy)
1673
1674 ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
1675 ! Energy associated with the GLE thermostat
1676 cpassert(ASSOCIATED(thermostat%gle))
1677 ALLOCATE (thermo_energy(thermostat%gle%loc_num_gle))
1678 DO i = 1, thermostat%gle%loc_num_gle
1679 thermo_energy(i) = thermostat%gle%nvt(i)%thermostat_energy
1680 END DO
1681 CALL get_kin_energies(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
1682 thermostat%gle%glob_num_gle, thermo_energy, &
1683 thermostat_kin, para_env, array_pot, array_kin)
1684 DEALLOCATE (thermo_energy)
1685
1686 ![NB] nothing to do for Ad-Langevin?
1687
1688 END IF
1689 END IF
1690
1691 END SUBROUTINE get_thermostat_energies
1692
1693! **************************************************************************************************
1694!> \brief Calculates the temperatures for each region associated to a thermostat
1695!> \param thermostat ...
1696!> \param tot_temperature ...
1697!> \param para_env ...
1698!> \param array_temp ...
1699!> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1700! **************************************************************************************************
1701 SUBROUTINE get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
1702 TYPE(thermostat_type), POINTER :: thermostat
1703 REAL(kind=dp), INTENT(OUT) :: tot_temperature
1704 TYPE(mp_para_env_type), POINTER :: para_env
1705 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: array_temp
1706
1707 INTEGER :: i
1708 REAL(dp), ALLOCATABLE, DIMENSION(:) :: dof, nkt
1709
1710 IF (ASSOCIATED(thermostat)) THEN
1711 IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
1712 ! Energy associated with the Nose-Hoover thermostat
1713 cpassert(ASSOCIATED(thermostat%nhc))
1714 ALLOCATE (nkt(thermostat%nhc%loc_num_nhc))
1715 ALLOCATE (dof(thermostat%nhc%loc_num_nhc))
1716 DO i = 1, thermostat%nhc%loc_num_nhc
1717 nkt(i) = thermostat%nhc%nvt(1, i)%nkt
1718 dof(i) = real(thermostat%nhc%nvt(1, i)%degrees_of_freedom, kind=dp)
1719 END DO
1720 CALL get_temperatures(thermostat%nhc%map_info, thermostat%nhc%loc_num_nhc, &
1721 thermostat%nhc%glob_num_nhc, nkt, dof, para_env, tot_temperature, array_temp)
1722 DEALLOCATE (nkt)
1723 DEALLOCATE (dof)
1724 ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
1725 ! Energy associated with the CSVR thermostat
1726 cpassert(ASSOCIATED(thermostat%csvr))
1727
1728 ALLOCATE (nkt(thermostat%csvr%loc_num_csvr))
1729 ALLOCATE (dof(thermostat%csvr%loc_num_csvr))
1730 DO i = 1, thermostat%csvr%loc_num_csvr
1731 nkt(i) = thermostat%csvr%nvt(i)%nkt
1732 dof(i) = real(thermostat%csvr%nvt(i)%degrees_of_freedom, kind=dp)
1733 END DO
1734 CALL get_temperatures(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
1735 thermostat%csvr%glob_num_csvr, nkt, dof, para_env, tot_temperature, array_temp)
1736 DEALLOCATE (nkt)
1737 DEALLOCATE (dof)
1738 ELSE IF (thermostat%type_of_thermostat == do_thermo_al) THEN
1739 ! Energy associated with the AD_LANGEVIN thermostat
1740 cpassert(ASSOCIATED(thermostat%al))
1741
1742 ALLOCATE (nkt(thermostat%al%loc_num_al))
1743 ALLOCATE (dof(thermostat%al%loc_num_al))
1744 DO i = 1, thermostat%al%loc_num_al
1745 nkt(i) = thermostat%al%nvt(i)%nkt
1746 dof(i) = real(thermostat%al%nvt(i)%degrees_of_freedom, kind=dp)
1747 END DO
1748 CALL get_temperatures(thermostat%al%map_info, thermostat%al%loc_num_al, &
1749 thermostat%al%glob_num_al, nkt, dof, para_env, tot_temperature, array_temp)
1750 DEALLOCATE (nkt)
1751 DEALLOCATE (dof)
1752 ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
1753 ! Energy associated with the GLE thermostat
1754 cpassert(ASSOCIATED(thermostat%gle))
1755
1756 ALLOCATE (nkt(thermostat%gle%loc_num_gle))
1757 ALLOCATE (dof(thermostat%gle%loc_num_gle))
1758 DO i = 1, thermostat%gle%loc_num_gle
1759 nkt(i) = thermostat%gle%nvt(i)%nkt
1760 dof(i) = real(thermostat%gle%nvt(i)%degrees_of_freedom, kind=dp)
1761 END DO
1762 CALL get_temperatures(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
1763 thermostat%gle%glob_num_gle, nkt, dof, para_env, tot_temperature, array_temp)
1764 DEALLOCATE (nkt)
1765 DEALLOCATE (dof)
1766 END IF
1767 END IF
1768
1769 END SUBROUTINE get_region_temperatures
1770
1771! **************************************************************************************************
1772!> \brief Prints status of all thermostats during an MD run
1773!> \param thermostats ...
1774!> \param para_env ...
1775!> \param my_pos ...
1776!> \param my_act ...
1777!> \param itimes ...
1778!> \param time ...
1779!> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1780! **************************************************************************************************
1781 SUBROUTINE print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
1782 TYPE(thermostats_type), POINTER :: thermostats
1783 TYPE(mp_para_env_type), POINTER :: para_env
1784 CHARACTER(LEN=default_string_length) :: my_pos, my_act
1785 INTEGER, INTENT(IN) :: itimes
1786 REAL(kind=dp), INTENT(IN) :: time
1787
1788 IF (ASSOCIATED(thermostats)) THEN
1789 IF (ASSOCIATED(thermostats%thermostat_part)) THEN
1790 CALL print_thermostat_status(thermostats%thermostat_part, para_env, my_pos, my_act, itimes, time)
1791 END IF
1792 IF (ASSOCIATED(thermostats%thermostat_shell)) THEN
1793 CALL print_thermostat_status(thermostats%thermostat_shell, para_env, my_pos, my_act, itimes, time)
1794 END IF
1795 IF (ASSOCIATED(thermostats%thermostat_coef)) THEN
1796 CALL print_thermostat_status(thermostats%thermostat_coef, para_env, my_pos, my_act, itimes, time)
1797 END IF
1798 IF (ASSOCIATED(thermostats%thermostat_baro)) THEN
1799 CALL print_thermostat_status(thermostats%thermostat_baro, para_env, my_pos, my_act, itimes, time)
1800 END IF
1801 END IF
1802 END SUBROUTINE print_thermostats_status
1803
1804! **************************************************************************************************
1805!> \brief Prints status of a specific thermostat
1806!> \param thermostat ...
1807!> \param para_env ...
1808!> \param my_pos ...
1809!> \param my_act ...
1810!> \param itimes ...
1811!> \param time ...
1812!> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
1813! **************************************************************************************************
1814 SUBROUTINE print_thermostat_status(thermostat, para_env, my_pos, my_act, itimes, time)
1815 TYPE(thermostat_type), POINTER :: thermostat
1816 TYPE(mp_para_env_type), POINTER :: para_env
1817 CHARACTER(LEN=default_string_length) :: my_pos, my_act
1818 INTEGER, INTENT(IN) :: itimes
1819 REAL(kind=dp), INTENT(IN) :: time
1820
1821 INTEGER :: i, unit
1822 LOGICAL :: new_file
1823 REAL(kind=dp) :: thermo_kin, thermo_pot, tot_temperature
1824 REAL(kind=dp), DIMENSION(:), POINTER :: array_kin, array_pot, array_temp
1825 TYPE(cp_logger_type), POINTER :: logger
1826 TYPE(section_vals_type), POINTER :: print_key
1827
1828 NULLIFY (logger, print_key, array_pot, array_kin, array_temp)
1829 logger => cp_get_default_logger()
1830
1831 IF (ASSOCIATED(thermostat)) THEN
1832 ! Print Energies
1833 print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%ENERGY")
1834 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
1835 CALL get_thermostat_energies(thermostat, thermo_pot, thermo_kin, para_env, array_pot, array_kin)
1836 unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%ENERGY", &
1837 extension="."//trim(thermostat%label)//".tener", file_position=my_pos, &
1838 file_action=my_act, is_new_file=new_file)
1839 IF (unit > 0) THEN
1840 IF (new_file) THEN
1841 WRITE (unit, '(A)') "# Thermostat Potential and Kinetic Energies - Total and per Region"
1842 WRITE (unit, '("#",3X,A,2X,A,13X,A,10X,A)') "Step Nr.", "Time[fs]", "Kin.[a.u.]", "Pot.[a.u.]"
1843 END IF
1844 WRITE (unit=unit, fmt="(I8, F12.3,6X,2F20.10)") itimes, time*femtoseconds, thermo_kin, thermo_pot
1845 WRITE (unit, '(A,4F20.10)') "# KINETIC ENERGY REGIONS: ", array_kin(1:min(4, SIZE(array_kin)))
1846 DO i = 5, SIZE(array_kin), 4
1847 WRITE (unit=unit, fmt='("#",25X,4F20.10)') array_kin(i:min(i + 3, SIZE(array_kin)))
1848 END DO
1849 WRITE (unit, '(A,4F20.10)') "# POTENT. ENERGY REGIONS: ", array_pot(1:min(4, SIZE(array_pot)))
1850 DO i = 5, SIZE(array_pot), 4
1851 WRITE (unit=unit, fmt='("#",25X,4F20.10)') array_pot(i:min(i + 3, SIZE(array_pot)))
1852 END DO
1853 CALL m_flush(unit)
1854 END IF
1855 DEALLOCATE (array_kin)
1856 DEALLOCATE (array_pot)
1857 CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%ENERGY")
1858 END IF
1859 ! Print Temperatures of the regions
1860 print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%TEMPERATURE")
1861 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
1862 CALL get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
1863 unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%TEMPERATURE", &
1864 extension="."//trim(thermostat%label)//".temp", file_position=my_pos, &
1865 file_action=my_act, is_new_file=new_file)
1866 IF (unit > 0) THEN
1867 IF (new_file) THEN
1868 WRITE (unit, '(A)') "# Temperature Total and per Region"
1869 WRITE (unit, '("#",3X,A,2X,A,10X,A)') "Step Nr.", "Time[fs]", "Temp.[K]"
1870 END IF
1871 WRITE (unit=unit, fmt="(I8, F12.3,3X,F20.10)") itimes, time*femtoseconds, tot_temperature
1872 WRITE (unit, '(A,I10)') "# TEMPERATURE REGIONS: ", SIZE(array_temp)
1873 DO i = 1, SIZE(array_temp), 4
1874 WRITE (unit=unit, fmt='("#",22X,4F20.10)') array_temp(i:min(i + 3, SIZE(array_temp)))
1875 END DO
1876 CALL m_flush(unit)
1877 END IF
1878 DEALLOCATE (array_temp)
1879 CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%TEMPERATURE")
1880 END IF
1881 END IF
1882 END SUBROUTINE print_thermostat_status
1883
1884! **************************************************************************************************
1885!> \brief Handles the communication for thermostats (1D array)
1886!> \param array ...
1887!> \param number ...
1888!> \param para_env ...
1889!> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
1890! **************************************************************************************************
1891 SUBROUTINE communication_thermo_low1(array, number, para_env)
1892 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: array
1893 INTEGER, INTENT(IN) :: number
1894 TYPE(mp_para_env_type), POINTER :: para_env
1895
1896 INTEGER :: i, icheck, ncheck
1897 REAL(kind=dp), DIMENSION(:), POINTER :: work, work2
1898
1899 ALLOCATE (work(para_env%num_pe))
1900 DO i = 1, number
1901 work = 0.0_dp
1902 work(para_env%mepos + 1) = array(i)
1903 CALL para_env%sum(work)
1904 ncheck = count(work /= 0.0_dp)
1905 array(i) = 0.0_dp
1906 IF (ncheck /= 0) THEN
1907 ALLOCATE (work2(ncheck))
1908 ncheck = 0
1909 DO icheck = 1, para_env%num_pe
1910 IF (work(icheck) /= 0.0_dp) THEN
1911 ncheck = ncheck + 1
1912 work2(ncheck) = work(icheck)
1913 END IF
1914 END DO
1915 cpassert(ncheck == SIZE(work2))
1916 cpassert(all(work2 == work2(1)))
1917
1918 array(i) = work2(1)
1919 DEALLOCATE (work2)
1920 END IF
1921 END DO
1922 DEALLOCATE (work)
1923 END SUBROUTINE communication_thermo_low1
1924
1925! **************************************************************************************************
1926!> \brief Handles the communication for thermostats (2D array)
1927!> \param array ...
1928!> \param number1 ...
1929!> \param number2 ...
1930!> \param para_env ...
1931!> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
1932! **************************************************************************************************
1933 SUBROUTINE communication_thermo_low2(array, number1, number2, para_env)
1934 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: array
1935 INTEGER, INTENT(IN) :: number1, number2
1936 TYPE(mp_para_env_type), POINTER :: para_env
1937
1938 INTEGER :: i, icheck, j, ncheck
1939 INTEGER, DIMENSION(:, :), POINTER :: work, work2
1940
1941 ALLOCATE (work(number1, para_env%num_pe))
1942 DO i = 1, number2
1943 work = 0
1944 work(:, para_env%mepos + 1) = array(:, i)
1945 CALL para_env%sum(work)
1946 ncheck = 0
1947 DO j = 1, para_env%num_pe
1948 IF (any(work(:, j) /= 0)) THEN
1949 ncheck = ncheck + 1
1950 END IF
1951 END DO
1952 array(:, i) = 0
1953 IF (ncheck /= 0) THEN
1954 ALLOCATE (work2(number1, ncheck))
1955 ncheck = 0
1956 DO icheck = 1, para_env%num_pe
1957 IF (any(work(:, icheck) /= 0)) THEN
1958 ncheck = ncheck + 1
1959 work2(:, ncheck) = work(:, icheck)
1960 END IF
1961 END DO
1962 cpassert(ncheck == SIZE(work2, 2))
1963 DO j = 1, ncheck
1964 cpassert(all(work2(:, j) == work2(:, 1)))
1965 END DO
1966 array(:, i) = work2(:, 1)
1967 DEALLOCATE (work2)
1968 END IF
1969 END DO
1970 DEALLOCATE (work)
1971 END SUBROUTINE communication_thermo_low2
1972
1973END 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:113
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:302
input_constants::nvt_adiabatic_ensemble
integer, parameter, public nvt_adiabatic_ensemble
Definition input_constants.F:98
input_constants::nph_uniaxial_ensemble
integer, parameter, public nph_uniaxial_ensemble
Definition input_constants.F:98
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:98
input_constants::isokin_ensemble
integer, parameter, public isokin_ensemble
Definition input_constants.F:98
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:98
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:98
input_constants::langevin_ensemble
integer, parameter, public langevin_ensemble
Definition input_constants.F:98
input_constants::do_region_molecule
integer, parameter, public do_region_molecule
Definition input_constants.F:419
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:98
input_constants::do_thermo_al
integer, parameter, public do_thermo_al
Definition input_constants.F:113
input_constants::do_constr_molec
integer, parameter, public do_constr_molec
Definition input_constants.F:80
input_constants::do_region_thermal
integer, parameter, public do_region_thermal
Definition input_constants.F:419
input_constants::do_thermo_csvr
integer, parameter, public do_thermo_csvr
Definition input_constants.F:113
input_constants::do_thermo_gle
integer, parameter, public do_thermo_gle
Definition input_constants.F:113
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:98
input_constants::nve_ensemble
integer, parameter, public nve_ensemble
Definition input_constants.F:98
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:98
input_constants::do_region_massive
integer, parameter, public do_region_massive
Definition input_constants.F:419
input_constants::do_region_global
integer, parameter, public do_region_global
Definition input_constants.F:419
input_constants::do_region_defined
integer, parameter, public do_region_defined
Definition input_constants.F:419
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition input_constants.F:98
input_constants::nvt_ensemble
integer, parameter, public nvt_ensemble
Definition input_constants.F:98
input_constants::do_thermo_communication
integer, parameter, public do_thermo_communication
Definition input_constants.F:302
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:1116
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:1365
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:1934
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:1782
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:1181
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:1286
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:1646
thermostat_utils::ke_region_baro
subroutine, public ke_region_baro(map_info, npt, group)
...
Definition thermostat_utils.F:992
thermostat_utils::vel_rescale_baro
subroutine, public vel_rescale_baro(map_info, npt)
...
Definition thermostat_utils.F:1021
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:1050
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:1449
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:1518
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:155
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
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