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