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