(git:3a22a26)
Loading...
Searching...
No Matches
thermostat_mapping.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
10! **************************************************************************************************
11MODULE thermostat_mapping
12
13 USE cell_types, ONLY: use_perd_x,&
14 use_perd_xy,&
15 use_perd_xyz,&
16 use_perd_xz,&
17 use_perd_y,&
18 use_perd_yz,&
19 use_perd_z
20 USE distribution_1d_types, ONLY: distribution_1d_type
21 USE extended_system_types, ONLY: map_info_type
22 USE input_constants, ONLY: do_region_defined,&
23 do_region_global,&
24 do_region_massive,&
25 do_region_molecule,&
26 do_region_thermal,&
27 do_thermo_communication,&
28 do_thermo_no_communication
29 USE kinds, ONLY: dp
30 USE memory_utilities, ONLY: reallocate
31 USE message_passing, ONLY: mp_para_env_type
32 USE molecule_kind_types, ONLY: colvar_constraint_type,&
33 fixd_constraint_type,&
34 g3x3_constraint_type,&
35 g4x6_constraint_type,&
36 get_molecule_kind,&
37 molecule_kind_type
38 USE molecule_types, ONLY: get_molecule,&
39 global_constraint_type,&
40 molecule_type
41 USE simpar_types, ONLY: simpar_type
42 USE util, ONLY: locate,&
43 sort
44#include "../../base/base_uses.f90"
45
46 IMPLICIT NONE
47
48 PUBLIC :: thermostat_mapping_region, &
49 adiabatic_mapping_region, &
50 init_baro_map_info
51
52 PRIVATE
53 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermostat_mapping'
54
55CONTAINS
56! **************************************************************************************************
57!> \brief Main general setup for adiabatic thermostat regions (Nose only)
58!> \param map_info ...
59!> \param deg_of_freedom ...
60!> \param massive_atom_list ...
61!> \param molecule_kind_set ...
62!> \param local_molecules ...
63!> \param molecule_set ...
64!> \param para_env ...
65!> \param natoms_local ...
66!> \param simpar ...
67!> \param number ...
68!> \param region ...
69!> \param gci ...
70!> \param shell ...
71!> \param map_loc_thermo_gen ...
72!> \param sum_of_thermostats ...
73!> \author CJM - PNNL
74! **************************************************************************************************
75 SUBROUTINE adiabatic_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
76 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, &
77 number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats)
78
79 TYPE(map_info_type), POINTER :: map_info
80 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
81 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
82 TYPE(distribution_1d_type), POINTER :: local_molecules
83 TYPE(molecule_type), POINTER :: molecule_set(:)
84 TYPE(mp_para_env_type), POINTER :: para_env
85 INTEGER, INTENT(OUT) :: natoms_local
86 TYPE(simpar_type), POINTER :: simpar
87 INTEGER, INTENT(INOUT) :: number
88 INTEGER, INTENT(IN) :: region
89 TYPE(global_constraint_type), POINTER :: gci
90 LOGICAL, INTENT(IN) :: shell
91 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
92 INTEGER, INTENT(INOUT) :: sum_of_thermostats
93
94 CHARACTER(LEN=*), PARAMETER :: routinen = 'adiabatic_mapping_region'
95
96 INTEGER :: handle, nkind, nmol_local, nsize, &
97 number_of_thermostats
98 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
99 INTEGER, DIMENSION(:, :), POINTER :: point
100
101 CALL timeset(routinen, handle)
102
103 NULLIFY (const_mol, tot_const, point)
104 cpassert(.NOT. ASSOCIATED(deg_of_freedom))
105 cpassert(.NOT. ASSOCIATED(massive_atom_list))
106
107 nkind = SIZE(molecule_kind_set)
108 CALL adiabatic_region_evaluate(map_info%dis_type, natoms_local, nmol_local, &
109 const_mol, tot_const, point, local_molecules, molecule_kind_set, molecule_set, &
110 simpar, shell)
111
112 ! Now we can allocate the target array s_kin and p_kin..
113 SELECT CASE (region)
114 CASE (do_region_global, do_region_molecule, do_region_massive)
115 nsize = number
116 CASE DEFAULT
117 ! STOP PROGRAM
118 END SELECT
119 ALLOCATE (map_info%s_kin(nsize))
120 ALLOCATE (map_info%v_scale(nsize))
121 ALLOCATE (map_info%p_kin(3, natoms_local))
122 ALLOCATE (map_info%p_scale(3, natoms_local))
123! nullify thermostat pointers
124 ! Allocate index array to 1
125 ALLOCATE (map_info%index(1))
126 ALLOCATE (map_info%map_index(1))
127 ALLOCATE (deg_of_freedom(1))
128
129 CALL massive_list_generate(molecule_set, molecule_kind_set, &
130 local_molecules, para_env, massive_atom_list, region, shell)
131
132 CALL adiabatic_mapping_region_low(region, map_info, nkind, point, &
133 deg_of_freedom, local_molecules, const_mol, massive_atom_list, &
134 tot_const, molecule_set, number_of_thermostats, shell, gci, &
135 map_loc_thermo_gen)
136
137 number = number_of_thermostats
138 sum_of_thermostats = number
139 CALL para_env%sum(sum_of_thermostats)
140
141! check = (number==number_of_thermostats)
142! CPPrecondition(check,cp_fatal_level,routineP,failure)
143 DEALLOCATE (const_mol)
144 DEALLOCATE (tot_const)
145 DEALLOCATE (point)
146
147 CALL timestop(handle)
148
149 END SUBROUTINE adiabatic_mapping_region
150
151! **************************************************************************************************
152!> \brief Performs the real mapping for the thermostat region
153!> \param region ...
154!> \param map_info ...
155!> \param nkind ...
156!> \param point ...
157!> \param deg_of_freedom ...
158!> \param local_molecules ...
159!> \param const_mol ...
160!> \param massive_atom_list ...
161!> \param tot_const ...
162!> \param molecule_set ...
163!> \param ntherm ...
164!> \param shell ...
165!> \param gci ...
166!> \param map_loc_thermo_gen ...
167!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
168! **************************************************************************************************
169 SUBROUTINE adiabatic_mapping_region_low(region, map_info, nkind, point, &
170 deg_of_freedom, local_molecules, const_mol, massive_atom_list, tot_const, &
171 molecule_set, ntherm, shell, gci, map_loc_thermo_gen)
172
173 INTEGER, INTENT(IN) :: region
174 TYPE(map_info_type), POINTER :: map_info
175 INTEGER :: nkind
176 INTEGER, DIMENSION(:, :), POINTER :: point
177 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom
178 TYPE(distribution_1d_type), POINTER :: local_molecules
179 INTEGER, DIMENSION(:), POINTER :: const_mol, massive_atom_list, tot_const
180 TYPE(molecule_type), POINTER :: molecule_set(:)
181 INTEGER, INTENT(OUT) :: ntherm
182 LOGICAL, INTENT(IN) :: shell
183 TYPE(global_constraint_type), POINTER :: gci
184 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
185
186 CHARACTER(LEN=*), PARAMETER :: routinen = 'adiabatic_mapping_region_low'
187
188 INTEGER :: first_atom, first_shell, glob_therm_num, handle, icount, ielement, ii, ikind, &
189 imol, imol_local, ipart, jj, k, kk, last_atom, last_shell, nglob_cns, nmol_local, number
190 LOGICAL :: check, global_constraints, &
191 have_thermostat
192 REAL(dp), SAVE, TARGET :: unity
193 TYPE(molecule_type), POINTER :: molecule
194
195 CALL timeset(routinen, handle)
196 unity = 1.0_dp
197 global_constraints = ASSOCIATED(gci)
198 deg_of_freedom = 0
199 icount = 0
200 number = 0
201 ntherm = 0
202 nglob_cns = 0
203 IF (global_constraints) nglob_cns = gci%ntot - gci%nrestraint
204 IF (region == do_region_global) THEN
205 ! Global Region
206 check = (map_info%dis_type == do_thermo_communication)
207 cpassert(check)
208 DO ikind = 1, nkind
209 DO jj = point(1, ikind), point(2, ikind)
210 IF (map_loc_thermo_gen(jj) /= huge(0)) THEN
211 DO ii = 1, 3
212 map_info%p_kin(ii, jj)%point => map_info%s_kin(1)
213 map_info%p_scale(ii, jj)%point => map_info%v_scale(1)
214 END DO
215 ELSE
216 DO ii = 1, 3
217 NULLIFY (map_info%p_kin(ii, jj)%point)
218 map_info%p_scale(ii, jj)%point => unity
219 END DO
220 END IF
221 END DO
222 deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind)
223 map_info%index(1) = 1
224 map_info%map_index(1) = 1
225 number = 1
226 END DO
227 deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns
228 ELSE IF (region == do_region_molecule) THEN
229 ! Molecular Region
230 IF (map_info%dis_type == do_thermo_no_communication) THEN
231 ! This is the standard case..
232 DO ikind = 1, nkind
233 nmol_local = local_molecules%n_el(ikind)
234 DO imol_local = 1, nmol_local
235 imol = local_molecules%list(ikind)%array(imol_local)
236 number = number + 1
237 have_thermostat = .true.
238! determine if the local molecule belongs to a thermostat
239 DO kk = point(1, number), point(2, number)
240 ! WRITE ( *, * ) 'kk', size(map_loc_thermo_gen), kk, map_loc_thermo_gen ( kk )
241 IF (map_loc_thermo_gen(kk) == huge(0)) THEN
242 have_thermostat = .false.
243 EXIT
244 END IF
245 END DO
246! If molecule has a thermostat, then map
247 IF (have_thermostat) THEN
248! glob_therm_num is the global thermostat number attached to the local molecule
249! We can test to make sure all atoms in the local molecule belong to the same
250! global thermostat as a way to detect errors.
251 glob_therm_num = map_loc_thermo_gen(point(1, number))
252 ntherm = ntherm + 1
253 CALL reallocate(map_info%index, 1, ntherm)
254 CALL reallocate(map_info%map_index, 1, ntherm)
255 CALL reallocate(deg_of_freedom, 1, ntherm)
256 map_info%index(ntherm) = glob_therm_num
257 map_info%map_index(ntherm) = ntherm
258 deg_of_freedom(ntherm) = const_mol(number)
259 DO kk = point(1, number), point(2, number)
260 DO jj = 1, 3
261 map_info%p_kin(jj, kk)%point => map_info%s_kin(ntherm)
262 map_info%p_scale(jj, kk)%point => map_info%v_scale(ntherm)
263 END DO
264 END DO
265! If molecule has no thermostat, then nullify pointers to that atom
266 ELSE
267 DO kk = point(1, number), point(2, number)
268 DO jj = 1, 3
269 NULLIFY (map_info%p_kin(jj, kk)%point)
270 map_info%p_scale(jj, kk)%point => unity
271 END DO
272 END DO
273 END IF
274 END DO
275 END DO
276 ELSE IF (map_info%dis_type == do_thermo_communication) THEN
277 ! This case is quite rare and happens only when we have one molecular
278 ! kind and one molecule..
279 cpassert(nkind == 1)
280 number = number + 1
281 ntherm = ntherm + 1
282 map_info%index(ntherm) = ntherm
283 map_info%map_index(ntherm) = ntherm
284 deg_of_freedom(ntherm) = deg_of_freedom(ntherm) + tot_const(nkind)
285 DO kk = point(1, nkind), point(2, nkind)
286 IF (map_loc_thermo_gen(kk) /= huge(0)) THEN
287 DO jj = 1, 3
288 map_info%p_kin(jj, kk)%point => map_info%s_kin(ntherm)
289 map_info%p_scale(jj, kk)%point => map_info%v_scale(ntherm)
290 END DO
291 ELSE
292 END IF
293 DO jj = 1, 3
294 NULLIFY (map_info%p_kin(jj, kk)%point)
295 map_info%p_scale(jj, kk)%point => unity
296 END DO
297 END DO
298 ELSE
299 cpabort("")
300 END IF
301 IF (nglob_cns /= 0) THEN
302 cpabort("Molecular thermostats with global constraints are impossible!")
303 END IF
304 ELSE IF (region == do_region_massive) THEN
305 ! Massive Region
306 check = (map_info%dis_type == do_thermo_no_communication)
307 cpassert(check)
308 DO ikind = 1, nkind
309 nmol_local = local_molecules%n_el(ikind)
310 DO imol_local = 1, nmol_local
311 icount = icount + 1
312 imol = local_molecules%list(ikind)%array(imol_local)
313 molecule => molecule_set(imol)
314 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom, &
315 first_shell=first_shell, last_shell=last_shell)
316 IF (shell) THEN
317 first_atom = first_shell
318 last_atom = last_shell
319 ELSE
320 IF ((tot_const(icount) > 0) .OR. (nglob_cns /= 0)) THEN
321 cpabort("Massive thermostats with constraints are impossible!")
322 END IF
323 END IF
324 k = 0
325
326 have_thermostat = .true.
327 DO ii = point(1, icount), point(2, icount)
328 IF (map_loc_thermo_gen(ii) /= 1) THEN
329 have_thermostat = .false.
330 EXIT
331 END IF
332 END DO
333
334 IF (have_thermostat) THEN
335 DO ii = point(1, icount), point(2, icount)
336 ipart = first_atom + k
337 ielement = locate(massive_atom_list, ipart)
338 k = k + 1
339 DO jj = 1, 3
340 ntherm = ntherm + 1
341 CALL reallocate(map_info%index, 1, ntherm)
342 CALL reallocate(map_info%map_index, 1, ntherm)
343 map_info%index(ntherm) = (ielement - 1)*3 + jj
344 map_info%map_index(ntherm) = ntherm
345 map_info%p_kin(jj, ii)%point => map_info%s_kin(ntherm)
346 map_info%p_scale(jj, ii)%point => map_info%v_scale(ntherm)
347 END DO
348 END DO
349 ELSE
350 DO ii = point(1, icount), point(2, icount)
351 DO jj = 1, 3
352 NULLIFY (map_info%p_kin(jj, ii)%point)
353 map_info%p_scale(jj, ii)%point => unity
354 END DO
355 END DO
356 END IF
357 IF (first_atom + k - 1 /= last_atom) THEN
358 cpabort("Inconsistent mapping of particles")
359 END IF
360 END DO
361 END DO
362 ELSE
363 cpabort("Invalid region!")
364 END IF
365
366 CALL timestop(handle)
367
368 END SUBROUTINE adiabatic_mapping_region_low
369
370! **************************************************************************************************
371!> \brief creates the mapping between the system and the thermostats
372!> \param dis_type ...
373!> \param natoms_local ...
374!> \param nmol_local ...
375!> \param const_mol ...
376!> \param tot_const ...
377!> \param point ...
378!> \param local_molecules ...
379!> \param molecule_kind_set ...
380!> \param molecule_set ...
381!> \param simpar ...
382!> \param shell ...
383!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
384! **************************************************************************************************
385 SUBROUTINE adiabatic_region_evaluate(dis_type, natoms_local, nmol_local, const_mol, &
386 tot_const, point, local_molecules, molecule_kind_set, molecule_set, simpar, shell)
387 INTEGER, INTENT(IN) :: dis_type
388 INTEGER, INTENT(OUT) :: natoms_local, nmol_local
389 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
390 INTEGER, DIMENSION(:, :), POINTER :: point
391 TYPE(distribution_1d_type), POINTER :: local_molecules
392 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
393 TYPE(molecule_type), POINTER :: molecule_set(:)
394 TYPE(simpar_type), POINTER :: simpar
395 LOGICAL, INTENT(IN) :: shell
396
397 CHARACTER(LEN=*), PARAMETER :: routinen = 'adiabatic_region_evaluate'
398
399 INTEGER :: atm_offset, first_atom, handle, icount, ikind, ilist, imol, imol_local, katom, &
400 last_atom, natom, nc, nfixd, nkind, nmol_per_kind, nmolecule, nshell
401 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
402 TYPE(molecule_kind_type), POINTER :: molecule_kind
403 TYPE(molecule_type), POINTER :: molecule
404
405 CALL timeset(routinen, handle)
406
407 natoms_local = 0
408 nmol_local = 0
409 nkind = SIZE(molecule_kind_set)
410 NULLIFY (fixd_list, molecule_kind, molecule)
411 ! Compute the TOTAL number of molecules and atoms on THIS PROC and
412 ! TOTAL number of molecules of IKIND on THIS PROC
413 DO ikind = 1, nkind
414 molecule_kind => molecule_kind_set(ikind)
415 CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
416 IF (shell) THEN
417 IF (nshell /= 0) THEN
418 natoms_local = natoms_local + nshell*local_molecules%n_el(ikind)
419 nmol_local = nmol_local + local_molecules%n_el(ikind)
420 END IF
421 ELSE
422 natoms_local = natoms_local + natom*local_molecules%n_el(ikind)
423 nmol_local = nmol_local + local_molecules%n_el(ikind)
424 END IF
425 END DO
426
427 cpassert(.NOT. ASSOCIATED(const_mol))
428 cpassert(.NOT. ASSOCIATED(tot_const))
429 cpassert(.NOT. ASSOCIATED(point))
430 IF (dis_type == do_thermo_no_communication) THEN
431 ALLOCATE (const_mol(nmol_local))
432 ALLOCATE (tot_const(nmol_local))
433 ALLOCATE (point(2, nmol_local))
434
435 point(:, :) = 0
436 atm_offset = 0
437 icount = 0
438 DO ikind = 1, nkind
439 nmol_per_kind = local_molecules%n_el(ikind)
440 molecule_kind => molecule_kind_set(ikind)
441 CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
442 fixd_list=fixd_list, nshell=nshell)
443 IF (shell) natom = nshell
444 DO imol_local = 1, nmol_per_kind
445 icount = icount + 1
446 point(1, icount) = atm_offset + 1
447 point(2, icount) = atm_offset + natom
448 IF (.NOT. shell) THEN
449 ! nc keeps track of all constraints but not fixed ones..
450 ! Let's identify fixed atoms for this molecule
451 nfixd = 0
452 imol = local_molecules%list(ikind)%array(imol_local)
453 molecule => molecule_set(imol)
454 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
455 IF (ASSOCIATED(fixd_list)) THEN
456 DO katom = first_atom, last_atom
457 DO ilist = 1, SIZE(fixd_list)
458 IF ((katom == fixd_list(ilist)%fixd) .AND. &
459 (.NOT. fixd_list(ilist)%restraint%active)) THEN
460 SELECT CASE (fixd_list(ilist)%itype)
461 CASE (use_perd_x, use_perd_y, use_perd_z)
462 nfixd = nfixd + 1
463 CASE (use_perd_xy, use_perd_xz, use_perd_yz)
464 nfixd = nfixd + 2
465 CASE (use_perd_xyz)
466 nfixd = nfixd + 3
467 END SELECT
468 END IF
469 END DO
470 END DO
471 END IF
472 const_mol(icount) = nc + nfixd
473 tot_const(icount) = const_mol(icount)
474 END IF
475 atm_offset = point(2, icount)
476 END DO
477 END DO
478 ELSE IF (dis_type == do_thermo_communication) THEN
479 ALLOCATE (const_mol(nkind))
480 ALLOCATE (tot_const(nkind))
481 ALLOCATE (point(2, nkind))
482 point(:, :) = 0
483 atm_offset = 0
484 ! nc keeps track of all constraints but not fixed ones..
485 DO ikind = 1, nkind
486 nmol_per_kind = local_molecules%n_el(ikind)
487 molecule_kind => molecule_kind_set(ikind)
488 CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
489 nmolecule=nmolecule, nconstraint_fixd=nfixd, nshell=nshell)
490 IF (shell) natom = nshell
491 IF (.NOT. shell) THEN
492 const_mol(ikind) = nc
493 ! Let's consider the fixed atoms only for the total number of constraints
494 ! in case we are in REPLICATED/INTERACTING thermostats
495 tot_const(ikind) = const_mol(ikind)*nmolecule + nfixd
496 END IF
497 point(1, ikind) = atm_offset + 1
498 point(2, ikind) = atm_offset + natom*nmol_per_kind
499 atm_offset = point(2, ikind)
500 END DO
501 END IF
502 IF ((.NOT. simpar%constraint) .OR. shell) THEN
503 const_mol = 0
504 tot_const = 0
505 END IF
506
507 CALL timestop(handle)
508
509 END SUBROUTINE adiabatic_region_evaluate
510
511! **************************************************************************************************
512!> \brief Main general setup thermostat regions (thermostat independent)
513!> \param map_info ...
514!> \param deg_of_freedom ...
515!> \param massive_atom_list ...
516!> \param molecule_kind_set ...
517!> \param local_molecules ...
518!> \param molecule_set ...
519!> \param para_env ...
520!> \param natoms_local ...
521!> \param simpar ...
522!> \param number ...
523!> \param region ...
524!> \param gci ...
525!> \param shell ...
526!> \param map_loc_thermo_gen ...
527!> \param sum_of_thermostats ...
528!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
529! **************************************************************************************************
530 SUBROUTINE thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
531 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, &
532 number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats)
533
534 TYPE(map_info_type), POINTER :: map_info
535 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
536 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
537 TYPE(distribution_1d_type), POINTER :: local_molecules
538 TYPE(molecule_type), POINTER :: molecule_set(:)
539 TYPE(mp_para_env_type), POINTER :: para_env
540 INTEGER, INTENT(OUT) :: natoms_local
541 TYPE(simpar_type), POINTER :: simpar
542 INTEGER, INTENT(IN) :: number, region
543 TYPE(global_constraint_type), POINTER :: gci
544 LOGICAL, INTENT(IN) :: shell
545 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
546 INTEGER, INTENT(IN) :: sum_of_thermostats
547
548 CHARACTER(LEN=*), PARAMETER :: routinen = 'thermostat_mapping_region'
549
550 INTEGER :: handle, nkind, nmol_local, nsize, &
551 number_of_thermostats
552 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
553 INTEGER, DIMENSION(:, :), POINTER :: point
554 LOGICAL :: check
555
556 CALL timeset(routinen, handle)
557
558 NULLIFY (const_mol, tot_const, point)
559 cpassert(.NOT. ASSOCIATED(deg_of_freedom))
560 cpassert(.NOT. ASSOCIATED(massive_atom_list))
561
562 nkind = SIZE(molecule_kind_set)
563 CALL mapping_region_evaluate(map_info%dis_type, natoms_local, nmol_local, &
564 const_mol, tot_const, point, local_molecules, molecule_kind_set, molecule_set, &
565 region, simpar, shell, map_loc_thermo_gen, sum_of_thermostats, para_env)
566
567 ! Now we can allocate the target array s_kin and p_kin..
568 SELECT CASE (region)
569 CASE (do_region_global, do_region_molecule, do_region_massive)
570 nsize = number
571 CASE (do_region_defined, do_region_thermal)
572 nsize = sum_of_thermostats
573 END SELECT
574 ALLOCATE (map_info%s_kin(nsize))
575 ALLOCATE (map_info%v_scale(nsize))
576 ALLOCATE (map_info%p_kin(3, natoms_local))
577 ALLOCATE (map_info%p_scale(3, natoms_local))
578 ! Allocate index array
579 ALLOCATE (map_info%index(number))
580 ALLOCATE (map_info%map_index(number))
581 ALLOCATE (deg_of_freedom(number))
582
583 CALL massive_list_generate(molecule_set, molecule_kind_set, &
584 local_molecules, para_env, massive_atom_list, region, shell)
585
586 CALL thermostat_mapping_region_low(region, map_info, nkind, point, &
587 deg_of_freedom, local_molecules, const_mol, massive_atom_list, &
588 tot_const, molecule_set, number_of_thermostats, shell, gci, &
589 map_loc_thermo_gen)
590
591 check = (number == number_of_thermostats)
592 cpassert(check)
593 DEALLOCATE (const_mol)
594 DEALLOCATE (tot_const)
595 DEALLOCATE (point)
596
597 CALL timestop(handle)
598
599 END SUBROUTINE thermostat_mapping_region
600
601! **************************************************************************************************
602!> \brief Performs the real mapping for the thermostat region
603!> \param region ...
604!> \param map_info ...
605!> \param nkind ...
606!> \param point ...
607!> \param deg_of_freedom ...
608!> \param local_molecules ...
609!> \param const_mol ...
610!> \param massive_atom_list ...
611!> \param tot_const ...
612!> \param molecule_set ...
613!> \param number ...
614!> \param shell ...
615!> \param gci ...
616!> \param map_loc_thermo_gen ...
617!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
618! **************************************************************************************************
619 SUBROUTINE thermostat_mapping_region_low(region, map_info, nkind, point, &
620 deg_of_freedom, local_molecules, const_mol, massive_atom_list, tot_const, &
621 molecule_set, number, shell, gci, map_loc_thermo_gen)
622
623 INTEGER, INTENT(IN) :: region
624 TYPE(map_info_type), POINTER :: map_info
625 INTEGER :: nkind
626 INTEGER, DIMENSION(:, :), POINTER :: point
627 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom
628 TYPE(distribution_1d_type), POINTER :: local_molecules
629 INTEGER, DIMENSION(:), POINTER :: const_mol, massive_atom_list, tot_const
630 TYPE(molecule_type), POINTER :: molecule_set(:)
631 INTEGER, INTENT(OUT) :: number
632 LOGICAL, INTENT(IN) :: shell
633 TYPE(global_constraint_type), POINTER :: gci
634 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
635
636 CHARACTER(LEN=*), PARAMETER :: routinen = 'thermostat_mapping_region_low'
637
638 INTEGER :: first_atom, first_shell, handle, i, icount, ielement, ii, ikind, imap, imol, &
639 imol_local, ipart, itmp, jj, k, kk, last_atom, last_shell, nglob_cns, nmol_local
640 INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp, wrk
641 LOGICAL :: check, global_constraints
642 TYPE(molecule_type), POINTER :: molecule
643
644 CALL timeset(routinen, handle)
645
646 global_constraints = ASSOCIATED(gci)
647 deg_of_freedom = 0
648 icount = 0
649 number = 0
650 nglob_cns = 0
651 IF (global_constraints) nglob_cns = gci%ntot - gci%nrestraint
652 IF (region == do_region_global) THEN
653 ! Global Region
654 check = (map_info%dis_type == do_thermo_communication)
655 cpassert(check)
656 DO ikind = 1, nkind
657 DO jj = point(1, ikind), point(2, ikind)
658 DO ii = 1, 3
659 map_info%p_kin(ii, jj)%point => map_info%s_kin(1)
660 map_info%p_scale(ii, jj)%point => map_info%v_scale(1)
661 END DO
662 END DO
663 deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind)
664 map_info%index(1) = 1
665 map_info%map_index(1) = 1
666 number = 1
667 END DO
668 deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns
669 ELSE IF ((region == do_region_defined) .OR. (region == do_region_thermal)) THEN
670 ! User defined Region to thermostat
671 check = (map_info%dis_type == do_thermo_communication)
672 cpassert(check)
673 ! Lets' identify the matching of the local thermostat w.r.t. the global one
674 itmp = SIZE(map_loc_thermo_gen)
675 ALLOCATE (tmp(itmp))
676 ALLOCATE (wrk(itmp))
677 tmp(:) = map_loc_thermo_gen
678 CALL sort(tmp, itmp, wrk)
679 number = 1
680 map_info%index(number) = tmp(1)
681 map_info%map_index(number) = tmp(1)
682 deg_of_freedom(number) = tot_const(tmp(1))
683 DO i = 2, itmp
684 IF (tmp(i) /= tmp(i - 1)) THEN
685 number = number + 1
686 map_info%index(number) = tmp(i)
687 map_info%map_index(number) = tmp(i)
688 deg_of_freedom(number) = tot_const(tmp(i))
689 END IF
690 END DO
691 DEALLOCATE (tmp)
692 DEALLOCATE (wrk)
693 DO jj = 1, SIZE(map_loc_thermo_gen)
694 DO ii = 1, 3
695 imap = map_loc_thermo_gen(jj)
696 map_info%p_kin(ii, jj)%point => map_info%s_kin(imap)
697 map_info%p_scale(ii, jj)%point => map_info%v_scale(imap)
698 END DO
699 END DO
700 IF (nglob_cns /= 0) THEN
701 CALL cp_abort(__location__, &
702 "User Defined thermostats with global constraints not implemented!")
703 END IF
704 ELSE IF (region == do_region_molecule) THEN
705 ! Molecular Region
706 IF (map_info%dis_type == do_thermo_no_communication) THEN
707 ! This is the standard case..
708 DO ikind = 1, nkind
709 nmol_local = local_molecules%n_el(ikind)
710 DO imol_local = 1, nmol_local
711 imol = local_molecules%list(ikind)%array(imol_local)
712 number = number + 1
713 map_info%index(number) = imol
714 map_info%map_index(number) = number
715 deg_of_freedom(number) = const_mol(number)
716 DO kk = point(1, number), point(2, number)
717 DO jj = 1, 3
718 map_info%p_kin(jj, kk)%point => map_info%s_kin(number)
719 map_info%p_scale(jj, kk)%point => map_info%v_scale(number)
720 END DO
721 END DO
722 END DO
723 END DO
724 ELSE IF (map_info%dis_type == do_thermo_communication) THEN
725 ! This case is quite rare and happens only when we have one molecular
726 ! kind and one molecule..
727 cpassert(nkind == 1)
728 number = number + 1
729 map_info%index(number) = number
730 map_info%map_index(number) = number
731 deg_of_freedom(number) = deg_of_freedom(number) + tot_const(nkind)
732 DO kk = point(1, nkind), point(2, nkind)
733 DO jj = 1, 3
734 map_info%p_kin(jj, kk)%point => map_info%s_kin(number)
735 map_info%p_scale(jj, kk)%point => map_info%v_scale(number)
736 END DO
737 END DO
738 ELSE
739 cpabort("")
740 END IF
741 IF (nglob_cns /= 0) THEN
742 cpabort("Molecular thermostats with global constraints are impossible!")
743 END IF
744 ELSE IF (region == do_region_massive) THEN
745 ! Massive Region
746 check = (map_info%dis_type == do_thermo_no_communication)
747 cpassert(check)
748 DO ikind = 1, nkind
749 nmol_local = local_molecules%n_el(ikind)
750 DO imol_local = 1, nmol_local
751 icount = icount + 1
752 imol = local_molecules%list(ikind)%array(imol_local)
753 molecule => molecule_set(imol)
754 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom, &
755 first_shell=first_shell, last_shell=last_shell)
756 IF (shell) THEN
757 first_atom = first_shell
758 last_atom = last_shell
759 ELSE
760 IF ((tot_const(icount) > 0) .OR. (nglob_cns /= 0)) THEN
761 cpabort("Massive thermostats with constraints are impossible!")
762 END IF
763 END IF
764 k = 0
765 DO ii = point(1, icount), point(2, icount)
766 ipart = first_atom + k
767 ielement = locate(massive_atom_list, ipart)
768 k = k + 1
769 DO jj = 1, 3
770 number = number + 1
771 map_info%index(number) = (ielement - 1)*3 + jj
772 map_info%map_index(number) = number
773 map_info%p_kin(jj, ii)%point => map_info%s_kin(number)
774 map_info%p_scale(jj, ii)%point => map_info%v_scale(number)
775 END DO
776 END DO
777 IF (first_atom + k - 1 /= last_atom) THEN
778 cpabort("Inconsistent mapping of particles")
779 END IF
780 END DO
781 END DO
782 ELSE
783 cpabort("Invalid region!")
784 END IF
785
786 CALL timestop(handle)
787
788 END SUBROUTINE thermostat_mapping_region_low
789
790! **************************************************************************************************
791!> \brief creates the mapping between the system and the thermostats
792!> \param dis_type ...
793!> \param natoms_local ...
794!> \param nmol_local ...
795!> \param const_mol ...
796!> \param tot_const ...
797!> \param point ...
798!> \param local_molecules ...
799!> \param molecule_kind_set ...
800!> \param molecule_set ...
801!> \param region ...
802!> \param simpar ...
803!> \param shell ...
804!> \param map_loc_thermo_gen ...
805!> \param sum_of_thermostats ...
806!> \param para_env ...
807!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
808! **************************************************************************************************
809 SUBROUTINE mapping_region_evaluate(dis_type, natoms_local, nmol_local, const_mol, &
810 tot_const, point, local_molecules, molecule_kind_set, molecule_set, region, &
811 simpar, shell, map_loc_thermo_gen, sum_of_thermostats, para_env)
812 INTEGER, INTENT(IN) :: dis_type
813 INTEGER, INTENT(OUT) :: natoms_local, nmol_local
814 INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
815 INTEGER, DIMENSION(:, :), POINTER :: point
816 TYPE(distribution_1d_type), POINTER :: local_molecules
817 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
818 TYPE(molecule_type), POINTER :: molecule_set(:)
819 INTEGER, INTENT(IN) :: region
820 TYPE(simpar_type), POINTER :: simpar
821 LOGICAL, INTENT(IN) :: shell
822 INTEGER, DIMENSION(:), POINTER :: map_loc_thermo_gen
823 INTEGER, INTENT(IN) :: sum_of_thermostats
824 TYPE(mp_para_env_type), POINTER :: para_env
825
826 CHARACTER(LEN=*), PARAMETER :: routinen = 'mapping_region_evaluate'
827
828 INTEGER :: atm_offset, first_atom, handle, i, iatm, icount, id_region, ikind, ilist, imol, &
829 imol_local, j, jatm, katom, last_atom, natom, nc, nfixd, nkind, nmol_per_kind, nmolecule, &
830 nshell
831 TYPE(colvar_constraint_type), DIMENSION(:), &
832 POINTER :: colv_list
833 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
834 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
835 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
836 TYPE(molecule_kind_type), POINTER :: molecule_kind
837 TYPE(molecule_type), POINTER :: molecule
838
839 CALL timeset(routinen, handle)
840
841 natoms_local = 0
842 nmol_local = 0
843 nkind = SIZE(molecule_kind_set)
844 NULLIFY (fixd_list, molecule_kind, molecule, colv_list, g3x3_list, g4x6_list)
845 ! Compute the TOTAL number of molecules and atoms on THIS PROC and
846 ! TOTAL number of molecules of IKIND on THIS PROC
847 DO ikind = 1, nkind
848 molecule_kind => molecule_kind_set(ikind)
849 CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
850 IF (shell) THEN
851 IF (nshell /= 0) THEN
852 natoms_local = natoms_local + nshell*local_molecules%n_el(ikind)
853 nmol_local = nmol_local + local_molecules%n_el(ikind)
854 END IF
855 ELSE
856 natoms_local = natoms_local + natom*local_molecules%n_el(ikind)
857 nmol_local = nmol_local + local_molecules%n_el(ikind)
858 END IF
859 END DO
860
861 cpassert(.NOT. ASSOCIATED(const_mol))
862 cpassert(.NOT. ASSOCIATED(tot_const))
863 cpassert(.NOT. ASSOCIATED(point))
864 IF (dis_type == do_thermo_no_communication) THEN
865 ALLOCATE (const_mol(nmol_local))
866 ALLOCATE (tot_const(nmol_local))
867 ALLOCATE (point(2, nmol_local))
868
869 point(:, :) = 0
870 atm_offset = 0
871 icount = 0
872 DO ikind = 1, nkind
873 nmol_per_kind = local_molecules%n_el(ikind)
874 molecule_kind => molecule_kind_set(ikind)
875 CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
876 fixd_list=fixd_list, nshell=nshell)
877 IF (shell) natom = nshell
878 DO imol_local = 1, nmol_per_kind
879 icount = icount + 1
880 point(1, icount) = atm_offset + 1
881 point(2, icount) = atm_offset + natom
882 IF (.NOT. shell) THEN
883 ! nc keeps track of all constraints but not fixed ones..
884 ! Let's identify fixed atoms for this molecule
885 nfixd = 0
886 imol = local_molecules%list(ikind)%array(imol_local)
887 molecule => molecule_set(imol)
888 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
889 IF (ASSOCIATED(fixd_list)) THEN
890 DO katom = first_atom, last_atom
891 DO ilist = 1, SIZE(fixd_list)
892 IF ((katom == fixd_list(ilist)%fixd) .AND. &
893 (.NOT. fixd_list(ilist)%restraint%active)) THEN
894 SELECT CASE (fixd_list(ilist)%itype)
895 CASE (use_perd_x, use_perd_y, use_perd_z)
896 nfixd = nfixd + 1
897 CASE (use_perd_xy, use_perd_xz, use_perd_yz)
898 nfixd = nfixd + 2
899 CASE (use_perd_xyz)
900 nfixd = nfixd + 3
901 END SELECT
902 END IF
903 END DO
904 END DO
905 END IF
906 const_mol(icount) = nc + nfixd
907 tot_const(icount) = const_mol(icount)
908 END IF
909 atm_offset = point(2, icount)
910 END DO
911 END DO
912 ELSE IF (dis_type == do_thermo_communication) THEN
913 IF ((region == do_region_defined) .OR. (region == do_region_thermal)) THEN
914 ! Setup of the arbitrary region
915 ALLOCATE (tot_const(sum_of_thermostats))
916 ALLOCATE (point(2, 0))
917 ALLOCATE (const_mol(0))
918 atm_offset = 0
919 tot_const = 0
920 const_mol = 0
921 point = 0
922 DO ikind = 1, nkind
923 nmol_per_kind = local_molecules%n_el(ikind)
924 molecule_kind => molecule_kind_set(ikind)
925 CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
926 fixd_list=fixd_list, colv_list=colv_list, g3x3_list=g3x3_list, &
927 g4x6_list=g4x6_list, nshell=nshell)
928 IF (shell) natom = nshell
929 DO imol_local = 1, nmol_per_kind
930 IF (.NOT. shell) THEN
931 ! First if nc is not zero let's check if all atoms of a molecule
932 ! are in the same thermostatting region..
933 imol = local_molecules%list(ikind)%array(imol_local)
934 molecule => molecule_set(imol)
935 id_region = map_loc_thermo_gen(atm_offset + 1)
936 IF (all(map_loc_thermo_gen(atm_offset + 1:atm_offset + natom) == id_region)) THEN
937 ! All the atoms of a molecule are within the same thermostatting
938 ! region.. this is the easy case..
939 tot_const(id_region) = tot_const(id_region) + nc
940 ELSE
941 ! If not let's check the single constraints defined for this molecule
942 ! and continue only when atoms involved in the constraint belong to
943 ! the same thermostatting region
944 IF (ASSOCIATED(colv_list)) THEN
945 DO i = 1, SIZE(colv_list)
946 IF (.NOT. colv_list(i)%restraint%active) THEN
947 iatm = atm_offset + colv_list(i)%i_atoms(1)
948 DO j = 2, SIZE(colv_list(i)%i_atoms)
949 jatm = atm_offset + colv_list(i)%i_atoms(j)
950 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
951 CALL cp_abort(__location__, &
952 "User Defined Region: "// &
953 "A constraint (COLV) was defined between two thermostatting regions! "// &
954 "This is not allowed!")
955 END IF
956 END DO
957 id_region = map_loc_thermo_gen(iatm)
958 tot_const(id_region) = tot_const(id_region) + 1
959 END IF
960 END DO
961 END IF
962 IF (ASSOCIATED(g3x3_list)) THEN
963 DO i = 1, SIZE(g3x3_list)
964 IF (.NOT. g3x3_list(i)%restraint%active) THEN
965 iatm = atm_offset + g3x3_list(i)%a
966 jatm = atm_offset + g3x3_list(i)%b
967 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
968 CALL cp_abort(__location__, &
969 "User Defined Region: "// &
970 "A constraint (G3X3) was defined between two thermostatting regions! "// &
971 "This is not allowed!")
972 END IF
973 jatm = atm_offset + g3x3_list(i)%c
974 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
975 CALL cp_abort(__location__, &
976 "User Defined Region: "// &
977 "A constraint (G3X3) was defined between two thermostatting regions! "// &
978 "This is not allowed!")
979 END IF
980 END IF
981 id_region = map_loc_thermo_gen(iatm)
982 tot_const(id_region) = tot_const(id_region) + 3
983 END DO
984 END IF
985 IF (ASSOCIATED(g4x6_list)) THEN
986 DO i = 1, SIZE(g4x6_list)
987 IF (.NOT. g4x6_list(i)%restraint%active) THEN
988 iatm = atm_offset + g4x6_list(i)%a
989 jatm = atm_offset + g4x6_list(i)%b
990 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
991 CALL cp_abort(__location__, &
992 " User Defined Region: "// &
993 "A constraint (G4X6) was defined between two thermostatting regions! "// &
994 "This is not allowed!")
995 END IF
996 jatm = atm_offset + g4x6_list(i)%c
997 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
998 CALL cp_abort(__location__, &
999 " User Defined Region: "// &
1000 "A constraint (G4X6) was defined between two thermostatting regions! "// &
1001 "This is not allowed!")
1002 END IF
1003 jatm = atm_offset + g4x6_list(i)%d
1004 IF (map_loc_thermo_gen(iatm) /= map_loc_thermo_gen(jatm)) THEN
1005 CALL cp_abort(__location__, &
1006 " User Defined Region: "// &
1007 "A constraint (G4X6) was defined between two thermostatting regions! "// &
1008 "This is not allowed!")
1009 END IF
1010 END IF
1011 id_region = map_loc_thermo_gen(iatm)
1012 tot_const(id_region) = tot_const(id_region) + 6
1013 END DO
1014 END IF
1015 END IF
1016 ! Here we handle possibly fixed atoms
1017 IF (ASSOCIATED(fixd_list)) THEN
1018 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1019 iatm = 0
1020 DO katom = first_atom, last_atom
1021 iatm = iatm + 1
1022 DO ilist = 1, SIZE(fixd_list)
1023 IF ((katom == fixd_list(ilist)%fixd) .AND. &
1024 (.NOT. fixd_list(ilist)%restraint%active)) THEN
1025 id_region = map_loc_thermo_gen(atm_offset + iatm)
1026 SELECT CASE (fixd_list(ilist)%itype)
1027 CASE (use_perd_x, use_perd_y, use_perd_z)
1028 tot_const(id_region) = tot_const(id_region) + 1
1029 CASE (use_perd_xy, use_perd_xz, use_perd_yz)
1030 tot_const(id_region) = tot_const(id_region) + 2
1031 CASE (use_perd_xyz)
1032 tot_const(id_region) = tot_const(id_region) + 3
1033 END SELECT
1034 END IF
1035 END DO
1036 END DO
1037 END IF
1038 END IF
1039 atm_offset = atm_offset + natom
1040 END DO
1041 END DO
1042 CALL para_env%sum(tot_const)
1043 ELSE
1044 ALLOCATE (const_mol(nkind))
1045 ALLOCATE (tot_const(nkind))
1046 ALLOCATE (point(2, nkind))
1047 point(:, :) = 0
1048 atm_offset = 0
1049 ! nc keeps track of all constraints but not fixed ones..
1050 DO ikind = 1, nkind
1051 nmol_per_kind = local_molecules%n_el(ikind)
1052 molecule_kind => molecule_kind_set(ikind)
1053 CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
1054 nmolecule=nmolecule, nconstraint_fixd=nfixd, nshell=nshell)
1055 IF (shell) natom = nshell
1056 IF (.NOT. shell) THEN
1057 const_mol(ikind) = nc
1058 ! Let's consider the fixed atoms only for the total number of constraints
1059 ! in case we are in REPLICATED/INTERACTING thermostats
1060 tot_const(ikind) = const_mol(ikind)*nmolecule + nfixd
1061 END IF
1062 point(1, ikind) = atm_offset + 1
1063 point(2, ikind) = atm_offset + natom*nmol_per_kind
1064 atm_offset = point(2, ikind)
1065 END DO
1066 END IF
1067 END IF
1068 IF ((.NOT. simpar%constraint) .OR. shell) THEN
1069 const_mol = 0
1070 tot_const = 0
1071 END IF
1072
1073 CALL timestop(handle)
1074
1075 END SUBROUTINE mapping_region_evaluate
1076
1077! **************************************************************************************************
1078!> \brief ...
1079!> \param molecule_set ...
1080!> \param molecule_kind_set ...
1081!> \param local_molecules ...
1082!> \param para_env ...
1083!> \param massive_atom_list ...
1084!> \param region ...
1085!> \param shell ...
1086! **************************************************************************************************
1087 SUBROUTINE massive_list_generate(molecule_set, molecule_kind_set, &
1088 local_molecules, para_env, massive_atom_list, region, shell)
1089
1090 TYPE(molecule_type), POINTER :: molecule_set(:)
1091 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
1092 TYPE(distribution_1d_type), POINTER :: local_molecules
1093 TYPE(mp_para_env_type), POINTER :: para_env
1094 INTEGER, POINTER :: massive_atom_list(:)
1095 INTEGER, INTENT(IN) :: region
1096 LOGICAL, INTENT(IN) :: shell
1097
1098 CHARACTER(LEN=*), PARAMETER :: routinen = 'massive_list_generate'
1099
1100 INTEGER :: first_atom, first_shell, handle, i, ikind, imol, iproc, j, natom, ncount, nkind, &
1101 nmol_per_kind, nshell, num_massive_atm, num_massive_atm_local, offset
1102 INTEGER, DIMENSION(:), POINTER :: array_num_massive_atm, local_atm_list, &
1103 work
1104 TYPE(molecule_kind_type), POINTER :: molecule_kind
1105 TYPE(molecule_type), POINTER :: molecule
1106
1107 CALL timeset(routinen, handle)
1108
1109 num_massive_atm_local = 0
1110 NULLIFY (local_atm_list)
1111 CALL reallocate(local_atm_list, 1, num_massive_atm_local)
1112
1113 nkind = SIZE(molecule_kind_set)
1114 DO ikind = 1, nkind
1115 nmol_per_kind = local_molecules%n_el(ikind)
1116 DO imol = 1, nmol_per_kind
1117 i = local_molecules%list(ikind)%array(imol)
1118 molecule => molecule_set(i)
1119 molecule_kind => molecule%molecule_kind
1120 CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
1121 IF (region == do_region_massive) THEN
1122 IF (shell) THEN
1123 natom = nshell
1124 END IF
1125 num_massive_atm_local = num_massive_atm_local + natom
1126 CALL reallocate(local_atm_list, 1, num_massive_atm_local)
1127 CALL get_molecule(molecule, first_atom=first_atom, first_shell=first_shell)
1128 IF (shell) THEN
1129 first_atom = first_shell
1130 END IF
1131 DO j = 1, natom
1132 local_atm_list(num_massive_atm_local - natom + j) = first_atom - 1 + j
1133 END DO
1134 END IF
1135 END DO
1136 END DO
1137
1138 ALLOCATE (array_num_massive_atm(para_env%num_pe))
1139 CALL para_env%allgather(num_massive_atm_local, array_num_massive_atm)
1140
1141 num_massive_atm = sum(array_num_massive_atm)
1142 ALLOCATE (massive_atom_list(num_massive_atm))
1143
1144 offset = 0
1145 DO iproc = 1, para_env%num_pe
1146 ncount = array_num_massive_atm(iproc)
1147 ALLOCATE (work(ncount))
1148 IF (para_env%mepos == (iproc - 1)) THEN
1149 DO i = 1, ncount
1150 work(i) = local_atm_list(i)
1151 END DO
1152 ELSE
1153 work(:) = 0
1154 END IF
1155 CALL para_env%bcast(work, iproc - 1)
1156 DO i = 1, ncount
1157 massive_atom_list(offset + i) = work(i)
1158 END DO
1159 DEALLOCATE (work)
1160 offset = offset + array_num_massive_atm(iproc)
1161 END DO
1162
1163 ! Sort atom list
1164 ALLOCATE (work(num_massive_atm))
1165 CALL sort(massive_atom_list, num_massive_atm, work)
1166 DEALLOCATE (work)
1167
1168 DEALLOCATE (local_atm_list)
1169 DEALLOCATE (array_num_massive_atm)
1170
1171 CALL timestop(handle)
1172
1173 END SUBROUTINE massive_list_generate
1174
1175! **************************************************************************************************
1176!> \brief Initialize the map_info for barostat thermostat
1177!> \param map_info ...
1178!> \param ndeg ...
1179!> \param num_thermo ...
1180!> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
1181! **************************************************************************************************
1182 SUBROUTINE init_baro_map_info(map_info, ndeg, num_thermo)
1183
1184 TYPE(map_info_type), POINTER :: map_info
1185 INTEGER, INTENT(IN) :: ndeg, num_thermo
1186
1187 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_baro_map_info'
1188
1189 INTEGER :: handle, i
1190
1191 CALL timeset(routinen, handle)
1192
1193 ALLOCATE (map_info%s_kin(num_thermo))
1194 ALLOCATE (map_info%v_scale(num_thermo))
1195 ALLOCATE (map_info%p_kin(1, ndeg))
1196 ALLOCATE (map_info%p_scale(1, ndeg))
1197 ! Allocate the index array
1198 ALLOCATE (map_info%index(1))
1199 ALLOCATE (map_info%map_index(1))
1200
1201 ! Begin the mapping loop
1202 DO i = 1, ndeg
1203 map_info%p_kin(1, i)%point => map_info%s_kin(1)
1204 map_info%p_scale(1, i)%point => map_info%v_scale(1)
1205 END DO
1206 map_info%index(1) = 1
1207 map_info%map_index(1) = 1
1208
1209 CALL timestop(handle)
1210
1211 END SUBROUTINE init_baro_map_info
1212
1213END MODULE thermostat_mapping
1214
memory_utilities::reallocate
Definition memory_utilities.F:27
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition cell_types.F:42
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition cell_types.F:42
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition cell_types.F:42
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition cell_types.F:42
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition cell_types.F:42
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
extended_system_types
Lumps all possible extended system variables into one type for easy access and passing.
Definition extended_system_types.F:17
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_thermo_no_communication
integer, parameter, public do_thermo_no_communication
Definition input_constants.F:302
input_constants::do_region_molecule
integer, parameter, public do_region_molecule
Definition input_constants.F:419
input_constants::do_region_thermal
integer, parameter, public do_region_thermal
Definition input_constants.F:419
input_constants::do_region_massive
integer, parameter, public do_region_massive
Definition input_constants.F:419
input_constants::do_region_global
integer, parameter, public do_region_global
Definition input_constants.F:419
input_constants::do_region_defined
integer, parameter, public do_region_defined
Definition input_constants.F:419
input_constants::do_thermo_communication
integer, parameter, public do_thermo_communication
Definition input_constants.F:302
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:512
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:326
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
thermostat_mapping
Definition thermostat_mapping.F:11
thermostat_mapping::adiabatic_mapping_region
subroutine, public adiabatic_mapping_region(map_info, deg_of_freedom, massive_atom_list, molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats)
Main general setup for adiabatic thermostat regions (Nose only)
Definition thermostat_mapping.F:78
thermostat_mapping::thermostat_mapping_region
subroutine, public thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats)
Main general setup thermostat regions (thermostat independent)
Definition thermostat_mapping.F:533
thermostat_mapping::init_baro_map_info
subroutine, public init_baro_map_info(map_info, ndeg, num_thermo)
Initialize the map_info for barostat thermostat.
Definition thermostat_mapping.F:1183
util
All kind of helpful little routines.
Definition util.F:14
util::locate
pure integer function, public locate(array, x)
Purpose: Given an array array(1:n), and given a value x, a value x_index is returned which is the ind...
Definition util.F:61
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
extended_system_types::map_info_type
Definition extended_system_types.F:65
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
molecule_kind_types::colvar_constraint_type
Definition molecule_kind_types.F:111
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:140
molecule_kind_types::g3x3_constraint_type
Definition molecule_kind_types.F:121
molecule_kind_types::g4x6_constraint_type
Definition molecule_kind_types.F:127
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:151
molecule_types::global_constraint_type
Definition molecule_types.F:134
molecule_types::molecule_type
Definition molecule_types.F:159
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30