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