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