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