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