(git:34ef472)
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 ! **************************************************************************************************
11 MODULE 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 
54 CONTAINS
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 
1212 END MODULE thermostat_mapping
1213 
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