Line data Source code
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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
98 0 : INTEGER, DIMENSION(:, :), POINTER :: point
99 :
100 0 : CALL timeset(routineN, handle)
101 :
102 0 : NULLIFY (const_mol, tot_const, point)
103 0 : CPASSERT(.NOT. ASSOCIATED(deg_of_freedom))
104 0 : CPASSERT(.NOT. ASSOCIATED(massive_atom_list))
105 :
106 0 : 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 0 : simpar, shell)
110 :
111 : ! Now we can allocate the target array s_kin and p_kin..
112 0 : SELECT CASE (region)
113 : CASE (do_region_global, do_region_molecule, do_region_massive)
114 0 : nsize = number
115 : CASE DEFAULT
116 : ! STOP PROGRAM
117 : END SELECT
118 0 : ALLOCATE (map_info%s_kin(nsize))
119 0 : ALLOCATE (map_info%v_scale(nsize))
120 0 : ALLOCATE (map_info%p_kin(3, natoms_local))
121 0 : ALLOCATE (map_info%p_scale(3, natoms_local))
122 : ! nullify thermostat pointers
123 : ! Allocate index array to 1
124 0 : ALLOCATE (map_info%index(1))
125 0 : ALLOCATE (map_info%map_index(1))
126 0 : ALLOCATE (deg_of_freedom(1))
127 :
128 : CALL massive_list_generate(molecule_set, molecule_kind_set, &
129 0 : 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 0 : map_loc_thermo_gen)
135 :
136 0 : number = number_of_thermostats
137 0 : sum_of_thermostats = number
138 0 : CALL para_env%sum(sum_of_thermostats)
139 :
140 : ! check = (number==number_of_thermostats)
141 : ! CPPrecondition(check,cp_fatal_level,routineP,failure)
142 0 : DEALLOCATE (const_mol)
143 0 : DEALLOCATE (tot_const)
144 0 : DEALLOCATE (point)
145 :
146 0 : CALL timestop(handle)
147 :
148 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
195 0 : unity = 1.0_dp
196 0 : global_constraints = ASSOCIATED(gci)
197 0 : deg_of_freedom = 0
198 0 : icount = 0
199 0 : number = 0
200 0 : ntherm = 0
201 0 : nglob_cns = 0
202 0 : IF (global_constraints) nglob_cns = gci%ntot - gci%nrestraint
203 0 : IF (region == do_region_global) THEN
204 : ! Global Region
205 0 : check = (map_info%dis_type == do_thermo_communication)
206 0 : CPASSERT(check)
207 0 : DO ikind = 1, nkind
208 0 : DO jj = point(1, ikind), point(2, ikind)
209 0 : IF (map_loc_thermo_gen(jj) /= HUGE(0)) THEN
210 0 : DO ii = 1, 3
211 0 : map_info%p_kin(ii, jj)%point => map_info%s_kin(1)
212 0 : map_info%p_scale(ii, jj)%point => map_info%v_scale(1)
213 : END DO
214 : ELSE
215 0 : DO ii = 1, 3
216 0 : NULLIFY (map_info%p_kin(ii, jj)%point)
217 0 : map_info%p_scale(ii, jj)%point => unity
218 : END DO
219 : END IF
220 : END DO
221 0 : deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind)
222 0 : map_info%index(1) = 1
223 0 : map_info%map_index(1) = 1
224 0 : number = 1
225 : END DO
226 0 : deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns
227 0 : ELSE IF (region == do_region_molecule) THEN
228 : ! Molecular Region
229 0 : IF (map_info%dis_type == do_thermo_no_communication) THEN
230 : ! This is the standard case..
231 0 : DO ikind = 1, nkind
232 0 : nmol_local = local_molecules%n_el(ikind)
233 0 : DO imol_local = 1, nmol_local
234 0 : imol = local_molecules%list(ikind)%array(imol_local)
235 0 : number = number + 1
236 0 : have_thermostat = .TRUE.
237 : ! determine if the local molecule belongs to a thermostat
238 0 : DO kk = point(1, number), point(2, number)
239 : ! WRITE ( *, * ) 'kk', size(map_loc_thermo_gen), kk, map_loc_thermo_gen ( kk )
240 0 : 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 0 : 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 0 : glob_therm_num = map_loc_thermo_gen(point(1, number))
251 0 : ntherm = ntherm + 1
252 0 : CALL reallocate(map_info%index, 1, ntherm)
253 0 : CALL reallocate(map_info%map_index, 1, ntherm)
254 0 : CALL reallocate(deg_of_freedom, 1, ntherm)
255 0 : map_info%index(ntherm) = glob_therm_num
256 0 : map_info%map_index(ntherm) = ntherm
257 0 : deg_of_freedom(ntherm) = const_mol(number)
258 0 : DO kk = point(1, number), point(2, number)
259 0 : DO jj = 1, 3
260 0 : map_info%p_kin(jj, kk)%point => map_info%s_kin(ntherm)
261 0 : 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 0 : DO kk = point(1, number), point(2, number)
267 0 : DO jj = 1, 3
268 0 : NULLIFY (map_info%p_kin(jj, kk)%point)
269 0 : map_info%p_scale(jj, kk)%point => unity
270 : END DO
271 : END DO
272 : END IF
273 : END DO
274 : END DO
275 0 : 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 0 : CPASSERT(nkind == 1)
279 0 : number = number + 1
280 0 : ntherm = ntherm + 1
281 0 : map_info%index(ntherm) = ntherm
282 0 : map_info%map_index(ntherm) = ntherm
283 0 : deg_of_freedom(ntherm) = deg_of_freedom(ntherm) + tot_const(nkind)
284 0 : DO kk = point(1, nkind), point(2, nkind)
285 0 : IF (map_loc_thermo_gen(kk) /= HUGE(0)) THEN
286 0 : DO jj = 1, 3
287 0 : map_info%p_kin(jj, kk)%point => map_info%s_kin(ntherm)
288 0 : map_info%p_scale(jj, kk)%point => map_info%v_scale(ntherm)
289 : END DO
290 : ELSE
291 : END IF
292 0 : DO jj = 1, 3
293 0 : NULLIFY (map_info%p_kin(jj, kk)%point)
294 0 : map_info%p_scale(jj, kk)%point => unity
295 : END DO
296 : END DO
297 : ELSE
298 0 : CPABORT("")
299 : END IF
300 0 : IF (nglob_cns /= 0) THEN
301 0 : CPABORT("Molecular thermostats with global constraints are impossible!")
302 : END IF
303 0 : ELSE IF (region == do_region_massive) THEN
304 : ! Massive Region
305 0 : check = (map_info%dis_type == do_thermo_no_communication)
306 0 : CPASSERT(check)
307 0 : DO ikind = 1, nkind
308 0 : nmol_local = local_molecules%n_el(ikind)
309 0 : DO imol_local = 1, nmol_local
310 0 : icount = icount + 1
311 0 : imol = local_molecules%list(ikind)%array(imol_local)
312 0 : molecule => molecule_set(imol)
313 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom, &
314 0 : first_shell=first_shell, last_shell=last_shell)
315 0 : IF (shell) THEN
316 0 : first_atom = first_shell
317 0 : last_atom = last_shell
318 : ELSE
319 0 : IF ((tot_const(icount) > 0) .OR. (nglob_cns /= 0)) THEN
320 0 : CPABORT("Massive thermostats with constraints are impossible!")
321 : END IF
322 : END IF
323 0 : k = 0
324 :
325 0 : have_thermostat = .TRUE.
326 0 : DO ii = point(1, icount), point(2, icount)
327 0 : IF (map_loc_thermo_gen(ii) /= 1) THEN
328 : have_thermostat = .FALSE.
329 : EXIT
330 : END IF
331 : END DO
332 :
333 0 : IF (have_thermostat) THEN
334 0 : DO ii = point(1, icount), point(2, icount)
335 0 : ipart = first_atom + k
336 0 : ielement = locate(massive_atom_list, ipart)
337 0 : k = k + 1
338 0 : DO jj = 1, 3
339 0 : ntherm = ntherm + 1
340 0 : CALL reallocate(map_info%index, 1, ntherm)
341 0 : CALL reallocate(map_info%map_index, 1, ntherm)
342 0 : map_info%index(ntherm) = (ielement - 1)*3 + jj
343 0 : map_info%map_index(ntherm) = ntherm
344 0 : map_info%p_kin(jj, ii)%point => map_info%s_kin(ntherm)
345 0 : map_info%p_scale(jj, ii)%point => map_info%v_scale(ntherm)
346 : END DO
347 : END DO
348 : ELSE
349 0 : DO ii = point(1, icount), point(2, icount)
350 0 : DO jj = 1, 3
351 0 : NULLIFY (map_info%p_kin(jj, ii)%point)
352 0 : map_info%p_scale(jj, ii)%point => unity
353 : END DO
354 : END DO
355 : END IF
356 0 : IF (first_atom + k - 1 /= last_atom) THEN
357 0 : CPABORT("Inconsistent mapping of particles")
358 : END IF
359 : END DO
360 : END DO
361 : ELSE
362 0 : CPABORT("Invalid region!")
363 : END IF
364 :
365 0 : CALL timestop(handle)
366 :
367 0 : 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 0 : 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 0 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
401 : TYPE(molecule_kind_type), POINTER :: molecule_kind
402 : TYPE(molecule_type), POINTER :: molecule
403 :
404 0 : CALL timeset(routineN, handle)
405 :
406 0 : natoms_local = 0
407 0 : nmol_local = 0
408 0 : nkind = SIZE(molecule_kind_set)
409 0 : 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 0 : DO ikind = 1, nkind
413 0 : molecule_kind => molecule_kind_set(ikind)
414 0 : CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
415 0 : IF (shell) THEN
416 0 : IF (nshell /= 0) THEN
417 0 : natoms_local = natoms_local + nshell*local_molecules%n_el(ikind)
418 0 : nmol_local = nmol_local + local_molecules%n_el(ikind)
419 : END IF
420 : ELSE
421 0 : natoms_local = natoms_local + natom*local_molecules%n_el(ikind)
422 0 : nmol_local = nmol_local + local_molecules%n_el(ikind)
423 : END IF
424 : END DO
425 :
426 0 : CPASSERT(.NOT. ASSOCIATED(const_mol))
427 0 : CPASSERT(.NOT. ASSOCIATED(tot_const))
428 0 : CPASSERT(.NOT. ASSOCIATED(point))
429 0 : IF (dis_type == do_thermo_no_communication) THEN
430 0 : ALLOCATE (const_mol(nmol_local))
431 0 : ALLOCATE (tot_const(nmol_local))
432 0 : ALLOCATE (point(2, nmol_local))
433 :
434 0 : point(:, :) = 0
435 : atm_offset = 0
436 : icount = 0
437 0 : DO ikind = 1, nkind
438 0 : nmol_per_kind = local_molecules%n_el(ikind)
439 0 : molecule_kind => molecule_kind_set(ikind)
440 : CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
441 0 : fixd_list=fixd_list, nshell=nshell)
442 0 : IF (shell) natom = nshell
443 0 : DO imol_local = 1, nmol_per_kind
444 0 : icount = icount + 1
445 0 : point(1, icount) = atm_offset + 1
446 0 : point(2, icount) = atm_offset + natom
447 0 : 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 0 : nfixd = 0
451 0 : imol = local_molecules%list(ikind)%array(imol_local)
452 0 : molecule => molecule_set(imol)
453 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
454 0 : IF (ASSOCIATED(fixd_list)) THEN
455 0 : DO katom = first_atom, last_atom
456 0 : DO ilist = 1, SIZE(fixd_list)
457 0 : IF ((katom == fixd_list(ilist)%fixd) .AND. &
458 0 : (.NOT. fixd_list(ilist)%restraint%active)) THEN
459 0 : SELECT CASE (fixd_list(ilist)%itype)
460 : CASE (use_perd_x, use_perd_y, use_perd_z)
461 0 : nfixd = nfixd + 1
462 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
463 0 : nfixd = nfixd + 2
464 : CASE (use_perd_xyz)
465 0 : nfixd = nfixd + 3
466 : END SELECT
467 : END IF
468 : END DO
469 : END DO
470 : END IF
471 0 : const_mol(icount) = nc + nfixd
472 0 : tot_const(icount) = const_mol(icount)
473 : END IF
474 0 : atm_offset = point(2, icount)
475 : END DO
476 : END DO
477 0 : ELSE IF (dis_type == do_thermo_communication) THEN
478 0 : ALLOCATE (const_mol(nkind))
479 0 : ALLOCATE (tot_const(nkind))
480 0 : ALLOCATE (point(2, nkind))
481 0 : point(:, :) = 0
482 : atm_offset = 0
483 : ! nc keeps track of all constraints but not fixed ones..
484 0 : DO ikind = 1, nkind
485 0 : nmol_per_kind = local_molecules%n_el(ikind)
486 0 : molecule_kind => molecule_kind_set(ikind)
487 : CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
488 0 : nmolecule=nmolecule, nconstraint_fixd=nfixd, nshell=nshell)
489 0 : IF (shell) natom = nshell
490 0 : IF (.NOT. shell) THEN
491 0 : 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 0 : tot_const(ikind) = const_mol(ikind)*nmolecule + nfixd
495 : END IF
496 0 : point(1, ikind) = atm_offset + 1
497 0 : point(2, ikind) = atm_offset + natom*nmol_per_kind
498 0 : atm_offset = point(2, ikind)
499 : END DO
500 : END IF
501 0 : IF ((.NOT. simpar%constraint) .OR. shell) THEN
502 0 : const_mol = 0
503 0 : tot_const = 0
504 : END IF
505 :
506 0 : CALL timestop(handle)
507 :
508 0 : 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 570 : 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 570 : INTEGER, DIMENSION(:), POINTER :: const_mol, tot_const
552 570 : INTEGER, DIMENSION(:, :), POINTER :: point
553 : LOGICAL :: check
554 :
555 570 : CALL timeset(routineN, handle)
556 :
557 570 : NULLIFY (const_mol, tot_const, point)
558 570 : CPASSERT(.NOT. ASSOCIATED(deg_of_freedom))
559 570 : CPASSERT(.NOT. ASSOCIATED(massive_atom_list))
560 :
561 570 : 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 570 : 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 1108 : SELECT CASE (region)
568 : CASE (do_region_global, do_region_molecule, do_region_massive)
569 538 : nsize = number
570 : CASE (do_region_defined)
571 570 : nsize = sum_of_thermostats
572 : END SELECT
573 1710 : ALLOCATE (map_info%s_kin(nsize))
574 1710 : ALLOCATE (map_info%v_scale(nsize))
575 1700 : ALLOCATE (map_info%p_kin(3, natoms_local))
576 1700 : ALLOCATE (map_info%p_scale(3, natoms_local))
577 : ! Allocate index array
578 1710 : ALLOCATE (map_info%index(number))
579 1710 : ALLOCATE (map_info%map_index(number))
580 1710 : ALLOCATE (deg_of_freedom(number))
581 :
582 : CALL massive_list_generate(molecule_set, molecule_kind_set, &
583 570 : 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 570 : map_loc_thermo_gen)
589 :
590 570 : check = (number == number_of_thermostats)
591 570 : CPASSERT(check)
592 570 : DEALLOCATE (const_mol)
593 570 : DEALLOCATE (tot_const)
594 570 : DEALLOCATE (point)
595 :
596 570 : CALL timestop(handle)
597 :
598 1710 : 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 570 : 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 570 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp, wrk
640 : LOGICAL :: check, global_constraints
641 : TYPE(molecule_type), POINTER :: molecule
642 :
643 570 : CALL timeset(routineN, handle)
644 :
645 570 : global_constraints = ASSOCIATED(gci)
646 77936 : deg_of_freedom = 0
647 570 : icount = 0
648 570 : number = 0
649 570 : nglob_cns = 0
650 570 : IF (global_constraints) nglob_cns = gci%ntot - gci%nrestraint
651 570 : IF (region == do_region_global) THEN
652 : ! Global Region
653 270 : check = (map_info%dis_type == do_thermo_communication)
654 270 : CPASSERT(check)
655 15372 : DO ikind = 1, nkind
656 39670 : DO jj = point(1, ikind), point(2, ikind)
657 113374 : DO ii = 1, 3
658 73704 : map_info%p_kin(ii, jj)%point => map_info%s_kin(1)
659 98272 : map_info%p_scale(ii, jj)%point => map_info%v_scale(1)
660 : END DO
661 : END DO
662 15102 : deg_of_freedom(1) = deg_of_freedom(1) + tot_const(ikind)
663 15102 : map_info%index(1) = 1
664 15102 : map_info%map_index(1) = 1
665 15372 : number = 1
666 : END DO
667 270 : deg_of_freedom(1) = deg_of_freedom(1) + nglob_cns
668 300 : ELSE IF (region == do_region_defined) THEN
669 : ! User defined Region to thermostat
670 32 : check = (map_info%dis_type == do_thermo_communication)
671 32 : CPASSERT(check)
672 : ! Lets' identify the matching of the local thermostat w.r.t. the global one
673 32 : itmp = SIZE(map_loc_thermo_gen)
674 96 : ALLOCATE (tmp(itmp))
675 96 : ALLOCATE (wrk(itmp))
676 22001 : tmp(:) = map_loc_thermo_gen
677 32 : CALL sort(tmp, itmp, wrk)
678 32 : number = 1
679 32 : map_info%index(number) = tmp(1)
680 32 : map_info%map_index(number) = tmp(1)
681 32 : deg_of_freedom(number) = tot_const(tmp(1))
682 21969 : DO i = 2, itmp
683 21969 : IF (tmp(i) /= tmp(i - 1)) THEN
684 40 : number = number + 1
685 40 : map_info%index(number) = tmp(i)
686 40 : map_info%map_index(number) = tmp(i)
687 40 : deg_of_freedom(number) = tot_const(tmp(i))
688 : END IF
689 : END DO
690 32 : DEALLOCATE (tmp)
691 32 : DEALLOCATE (wrk)
692 22001 : DO jj = 1, SIZE(map_loc_thermo_gen)
693 87908 : DO ii = 1, 3
694 65907 : imap = map_loc_thermo_gen(jj)
695 65907 : map_info%p_kin(ii, jj)%point => map_info%s_kin(imap)
696 87876 : map_info%p_scale(ii, jj)%point => map_info%v_scale(imap)
697 : END DO
698 : END DO
699 32 : IF (nglob_cns /= 0) THEN
700 : CALL cp_abort(__LOCATION__, &
701 0 : "User Defined thermostats with global constraints not implemented!")
702 : END IF
703 268 : ELSE IF (region == do_region_molecule) THEN
704 : ! Molecular Region
705 136 : IF (map_info%dis_type == do_thermo_no_communication) THEN
706 : ! This is the standard case..
707 6044 : DO ikind = 1, nkind
708 5912 : nmol_local = local_molecules%n_el(ikind)
709 13008 : DO imol_local = 1, nmol_local
710 6964 : imol = local_molecules%list(ikind)%array(imol_local)
711 6964 : number = number + 1
712 6964 : map_info%index(number) = imol
713 6964 : map_info%map_index(number) = number
714 6964 : deg_of_freedom(number) = const_mol(number)
715 28554 : DO kk = point(1, number), point(2, number)
716 69676 : DO jj = 1, 3
717 47034 : map_info%p_kin(jj, kk)%point => map_info%s_kin(number)
718 62712 : 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 4 : 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 4 : CPASSERT(nkind == 1)
727 4 : number = number + 1
728 4 : map_info%index(number) = number
729 4 : map_info%map_index(number) = number
730 4 : deg_of_freedom(number) = deg_of_freedom(number) + tot_const(nkind)
731 12 : DO kk = point(1, nkind), point(2, nkind)
732 36 : DO jj = 1, 3
733 24 : map_info%p_kin(jj, kk)%point => map_info%s_kin(number)
734 32 : map_info%p_scale(jj, kk)%point => map_info%v_scale(number)
735 : END DO
736 : END DO
737 : ELSE
738 0 : CPABORT("")
739 : END IF
740 136 : IF (nglob_cns /= 0) THEN
741 0 : CPABORT("Molecular thermostats with global constraints are impossible!")
742 : END IF
743 132 : ELSE IF (region == do_region_massive) THEN
744 : ! Massive Region
745 132 : check = (map_info%dis_type == do_thermo_no_communication)
746 132 : CPASSERT(check)
747 8882 : DO ikind = 1, nkind
748 8750 : nmol_local = local_molecules%n_el(ikind)
749 16128 : DO imol_local = 1, nmol_local
750 7246 : icount = icount + 1
751 7246 : imol = local_molecules%list(ikind)%array(imol_local)
752 7246 : molecule => molecule_set(imol)
753 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom, &
754 7246 : first_shell=first_shell, last_shell=last_shell)
755 7246 : IF (shell) THEN
756 904 : first_atom = first_shell
757 904 : last_atom = last_shell
758 : ELSE
759 6342 : IF ((tot_const(icount) > 0) .OR. (nglob_cns /= 0)) THEN
760 0 : CPABORT("Massive thermostats with constraints are impossible!")
761 : END IF
762 : END IF
763 7246 : k = 0
764 30598 : DO ii = point(1, icount), point(2, icount)
765 23352 : ipart = first_atom + k
766 23352 : ielement = locate(massive_atom_list, ipart)
767 23352 : k = k + 1
768 100654 : DO jj = 1, 3
769 70056 : number = number + 1
770 70056 : map_info%index(number) = (ielement - 1)*3 + jj
771 70056 : map_info%map_index(number) = number
772 70056 : map_info%p_kin(jj, ii)%point => map_info%s_kin(number)
773 93408 : map_info%p_scale(jj, ii)%point => map_info%v_scale(number)
774 : END DO
775 : END DO
776 15996 : IF (first_atom + k - 1 /= last_atom) THEN
777 0 : CPABORT("Inconsistent mapping of particles")
778 : END IF
779 : END DO
780 : END DO
781 : ELSE
782 0 : CPABORT("Invalid region!")
783 : END IF
784 :
785 570 : CALL timestop(handle)
786 :
787 570 : 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 570 : 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 570 : POINTER :: colv_list
832 570 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
833 570 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
834 570 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
835 : TYPE(molecule_kind_type), POINTER :: molecule_kind
836 : TYPE(molecule_type), POINTER :: molecule
837 :
838 570 : CALL timeset(routineN, handle)
839 :
840 570 : natoms_local = 0
841 570 : nmol_local = 0
842 570 : nkind = SIZE(molecule_kind_set)
843 570 : 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 30404 : DO ikind = 1, nkind
847 29834 : molecule_kind => molecule_kind_set(ikind)
848 29834 : CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
849 30404 : IF (shell) THEN
850 82 : IF (nshell /= 0) THEN
851 82 : natoms_local = natoms_local + nshell*local_molecules%n_el(ikind)
852 82 : nmol_local = nmol_local + local_molecules%n_el(ikind)
853 : END IF
854 : ELSE
855 29752 : natoms_local = natoms_local + natom*local_molecules%n_el(ikind)
856 29752 : nmol_local = nmol_local + local_molecules%n_el(ikind)
857 : END IF
858 : END DO
859 :
860 570 : CPASSERT(.NOT. ASSOCIATED(const_mol))
861 570 : CPASSERT(.NOT. ASSOCIATED(tot_const))
862 570 : CPASSERT(.NOT. ASSOCIATED(point))
863 570 : IF (dis_type == do_thermo_no_communication) THEN
864 792 : ALLOCATE (const_mol(nmol_local))
865 792 : ALLOCATE (tot_const(nmol_local))
866 792 : ALLOCATE (point(2, nmol_local))
867 :
868 42894 : point(:, :) = 0
869 : atm_offset = 0
870 : icount = 0
871 14926 : DO ikind = 1, nkind
872 14662 : nmol_per_kind = local_molecules%n_el(ikind)
873 14662 : molecule_kind => molecule_kind_set(ikind)
874 : CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
875 14662 : fixd_list=fixd_list, nshell=nshell)
876 14662 : IF (shell) natom = nshell
877 43798 : DO imol_local = 1, nmol_per_kind
878 14210 : icount = icount + 1
879 14210 : point(1, icount) = atm_offset + 1
880 14210 : point(2, icount) = atm_offset + natom
881 14210 : 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 13018 : nfixd = 0
885 13018 : imol = local_molecules%list(ikind)%array(imol_local)
886 13018 : molecule => molecule_set(imol)
887 13018 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
888 13018 : IF (ASSOCIATED(fixd_list)) THEN
889 1766 : DO katom = first_atom, last_atom
890 282600 : DO ilist = 1, SIZE(fixd_list)
891 280834 : IF ((katom == fixd_list(ilist)%fixd) .AND. &
892 1395 : (.NOT. fixd_list(ilist)%restraint%active)) THEN
893 1904 : SELECT CASE (fixd_list(ilist)%itype)
894 : CASE (use_perd_x, use_perd_y, use_perd_z)
895 768 : nfixd = nfixd + 1
896 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
897 256 : nfixd = nfixd + 2
898 : CASE (use_perd_xyz)
899 1136 : nfixd = nfixd + 3
900 : END SELECT
901 : END IF
902 : END DO
903 : END DO
904 : END IF
905 13018 : const_mol(icount) = nc + nfixd
906 13018 : tot_const(icount) = const_mol(icount)
907 : END IF
908 28872 : atm_offset = point(2, icount)
909 : END DO
910 : END DO
911 306 : ELSE IF (dis_type == do_thermo_communication) THEN
912 306 : IF (region == do_region_defined) THEN
913 : ! Setup of the arbitrary region
914 96 : ALLOCATE (tot_const(sum_of_thermostats))
915 32 : ALLOCATE (point(2, 0))
916 32 : ALLOCATE (const_mol(0))
917 32 : atm_offset = 0
918 116 : tot_const = 0
919 32 : const_mol = 0
920 32 : point = 0
921 98 : DO ikind = 1, nkind
922 66 : nmol_per_kind = local_molecules%n_el(ikind)
923 66 : 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 66 : g4x6_list=g4x6_list, nshell=nshell)
927 66 : IF (shell) natom = nshell
928 7450 : DO imol_local = 1, nmol_per_kind
929 7286 : 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 7286 : imol = local_molecules%list(ikind)%array(imol_local)
933 7286 : molecule => molecule_set(imol)
934 7286 : id_region = map_loc_thermo_gen(atm_offset + 1)
935 29178 : 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 7253 : 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 33 : IF (ASSOCIATED(colv_list)) THEN
944 64 : DO i = 1, SIZE(colv_list)
945 64 : IF (.NOT. colv_list(i)%restraint%active) THEN
946 32 : iatm = atm_offset + colv_list(i)%i_atoms(1)
947 64 : DO j = 2, SIZE(colv_list(i)%i_atoms)
948 32 : jatm = atm_offset + colv_list(i)%i_atoms(j)
949 64 : 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 0 : "This is not allowed!")
954 : END IF
955 : END DO
956 32 : id_region = map_loc_thermo_gen(iatm)
957 32 : tot_const(id_region) = tot_const(id_region) + 1
958 : END IF
959 : END DO
960 : END IF
961 33 : IF (ASSOCIATED(g3x3_list)) THEN
962 1 : DO i = 1, SIZE(g3x3_list)
963 0 : IF (.NOT. g3x3_list(i)%restraint%active) THEN
964 0 : iatm = atm_offset + g3x3_list(i)%a
965 0 : jatm = atm_offset + g3x3_list(i)%b
966 0 : 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 0 : "This is not allowed!")
971 : END IF
972 0 : jatm = atm_offset + g3x3_list(i)%c
973 0 : 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 0 : "This is not allowed!")
978 : END IF
979 : END IF
980 0 : id_region = map_loc_thermo_gen(iatm)
981 1 : tot_const(id_region) = tot_const(id_region) + 3
982 : END DO
983 : END IF
984 33 : IF (ASSOCIATED(g4x6_list)) THEN
985 0 : DO i = 1, SIZE(g4x6_list)
986 0 : IF (.NOT. g4x6_list(i)%restraint%active) THEN
987 0 : iatm = atm_offset + g4x6_list(i)%a
988 0 : jatm = atm_offset + g4x6_list(i)%b
989 0 : 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 0 : "This is not allowed!")
994 : END IF
995 0 : jatm = atm_offset + g4x6_list(i)%c
996 0 : 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 0 : "This is not allowed!")
1001 : END IF
1002 0 : jatm = atm_offset + g4x6_list(i)%d
1003 0 : 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 0 : "This is not allowed!")
1008 : END IF
1009 : END IF
1010 0 : id_region = map_loc_thermo_gen(iatm)
1011 0 : 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 7286 : IF (ASSOCIATED(fixd_list)) THEN
1017 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
1018 0 : iatm = 0
1019 0 : DO katom = first_atom, last_atom
1020 0 : iatm = iatm + 1
1021 0 : DO ilist = 1, SIZE(fixd_list)
1022 0 : IF ((katom == fixd_list(ilist)%fixd) .AND. &
1023 0 : (.NOT. fixd_list(ilist)%restraint%active)) THEN
1024 0 : id_region = map_loc_thermo_gen(atm_offset + iatm)
1025 0 : SELECT CASE (fixd_list(ilist)%itype)
1026 : CASE (use_perd_x, use_perd_y, use_perd_z)
1027 0 : tot_const(id_region) = tot_const(id_region) + 1
1028 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
1029 0 : tot_const(id_region) = tot_const(id_region) + 2
1030 : CASE (use_perd_xyz)
1031 0 : 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 7352 : atm_offset = atm_offset + natom
1039 : END DO
1040 : END DO
1041 200 : CALL para_env%sum(tot_const)
1042 : ELSE
1043 822 : ALLOCATE (const_mol(nkind))
1044 822 : ALLOCATE (tot_const(nkind))
1045 822 : ALLOCATE (point(2, nkind))
1046 45592 : point(:, :) = 0
1047 : atm_offset = 0
1048 : ! nc keeps track of all constraints but not fixed ones..
1049 15380 : DO ikind = 1, nkind
1050 15106 : nmol_per_kind = local_molecules%n_el(ikind)
1051 15106 : molecule_kind => molecule_kind_set(ikind)
1052 : CALL get_molecule_kind(molecule_kind, nconstraint=nc, natom=natom, &
1053 15106 : nmolecule=nmolecule, nconstraint_fixd=nfixd, nshell=nshell)
1054 15106 : IF (shell) natom = nshell
1055 15106 : IF (.NOT. shell) THEN
1056 15082 : 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 15082 : tot_const(ikind) = const_mol(ikind)*nmolecule + nfixd
1060 : END IF
1061 15106 : point(1, ikind) = atm_offset + 1
1062 15106 : point(2, ikind) = atm_offset + natom*nmol_per_kind
1063 30486 : atm_offset = point(2, ikind)
1064 : END DO
1065 : END IF
1066 : END IF
1067 570 : IF ((.NOT. simpar%constraint) .OR. shell) THEN
1068 27585 : const_mol = 0
1069 27645 : tot_const = 0
1070 : END IF
1071 :
1072 570 : CALL timestop(handle)
1073 :
1074 570 : 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 570 : 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 570 : INTEGER, DIMENSION(:), POINTER :: array_num_massive_atm, local_atm_list, &
1102 570 : work
1103 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1104 : TYPE(molecule_type), POINTER :: molecule
1105 :
1106 570 : CALL timeset(routineN, handle)
1107 :
1108 570 : num_massive_atm_local = 0
1109 570 : NULLIFY (local_atm_list)
1110 570 : CALL reallocate(local_atm_list, 1, num_massive_atm_local)
1111 :
1112 570 : nkind = SIZE(molecule_kind_set)
1113 30404 : DO ikind = 1, nkind
1114 29834 : nmol_per_kind = local_molecules%n_el(ikind)
1115 64849 : DO imol = 1, nmol_per_kind
1116 34445 : i = local_molecules%list(ikind)%array(imol)
1117 34445 : molecule => molecule_set(i)
1118 34445 : molecule_kind => molecule%molecule_kind
1119 34445 : CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
1120 64279 : IF (region == do_region_massive) THEN
1121 7246 : IF (shell) THEN
1122 904 : natom = nshell
1123 : END IF
1124 7246 : num_massive_atm_local = num_massive_atm_local + natom
1125 7246 : CALL reallocate(local_atm_list, 1, num_massive_atm_local)
1126 7246 : CALL get_molecule(molecule, first_atom=first_atom, first_shell=first_shell)
1127 7246 : IF (shell) THEN
1128 904 : first_atom = first_shell
1129 : END IF
1130 30598 : DO j = 1, natom
1131 30598 : 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 1710 : ALLOCATE (array_num_massive_atm(para_env%num_pe))
1138 1710 : CALL para_env%allgather(num_massive_atm_local, array_num_massive_atm)
1139 :
1140 1710 : num_massive_atm = SUM(array_num_massive_atm)
1141 1272 : ALLOCATE (massive_atom_list(num_massive_atm))
1142 :
1143 570 : offset = 0
1144 1710 : DO iproc = 1, para_env%num_pe
1145 1140 : ncount = array_num_massive_atm(iproc)
1146 2544 : ALLOCATE (work(ncount))
1147 1140 : IF (para_env%mepos == (iproc - 1)) THEN
1148 23922 : DO i = 1, ncount
1149 23922 : work(i) = local_atm_list(i)
1150 : END DO
1151 : ELSE
1152 23922 : work(:) = 0
1153 : END IF
1154 94548 : CALL para_env%bcast(work, iproc - 1)
1155 47844 : DO i = 1, ncount
1156 47844 : massive_atom_list(offset + i) = work(i)
1157 : END DO
1158 1140 : DEALLOCATE (work)
1159 1710 : offset = offset + array_num_massive_atm(iproc)
1160 : END DO
1161 :
1162 : ! Sort atom list
1163 1272 : ALLOCATE (work(num_massive_atm))
1164 570 : CALL sort(massive_atom_list, num_massive_atm, work)
1165 570 : DEALLOCATE (work)
1166 :
1167 570 : DEALLOCATE (local_atm_list)
1168 570 : DEALLOCATE (array_num_massive_atm)
1169 :
1170 570 : CALL timestop(handle)
1171 :
1172 570 : 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 152 : 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 152 : CALL timeset(routineN, handle)
1191 :
1192 456 : ALLOCATE (map_info%s_kin(num_thermo))
1193 456 : ALLOCATE (map_info%v_scale(num_thermo))
1194 456 : ALLOCATE (map_info%p_kin(1, ndeg))
1195 456 : ALLOCATE (map_info%p_scale(1, ndeg))
1196 : ! Allocate the index array
1197 152 : ALLOCATE (map_info%index(1))
1198 152 : ALLOCATE (map_info%map_index(1))
1199 :
1200 : ! Begin the mapping loop
1201 672 : DO i = 1, ndeg
1202 520 : map_info%p_kin(1, i)%point => map_info%s_kin(1)
1203 672 : map_info%p_scale(1, i)%point => map_info%v_scale(1)
1204 : END DO
1205 152 : map_info%index(1) = 1
1206 152 : map_info%map_index(1) = 1
1207 :
1208 152 : CALL timestop(handle)
1209 :
1210 152 : END SUBROUTINE init_baro_map_info
1211 :
1212 : END MODULE thermostat_mapping
1213 :
|
