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