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