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 : !> \par History
10 : !> CJM, 20-Feb-01
11 : !> JGH (10-Mar-2001)
12 : !> CJM (10-Apr-2001)
13 : !> \author CJM
14 : ! **************************************************************************************************
15 : MODULE extended_system_mapping
16 :
17 : USE distribution_1d_types, ONLY: distribution_1d_type
18 : USE extended_system_types, ONLY: debug_isotropic_limit,&
19 : lnhc_parameters_type,&
20 : map_info_type
21 : USE input_constants, ONLY: &
22 : do_thermo_communication, do_thermo_no_communication, do_thermo_only_master, &
23 : isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
24 : nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
25 : npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, reftraj_ensemble
26 : USE kinds, ONLY: dp
27 : USE message_passing, ONLY: mp_para_env_type
28 : USE molecule_kind_types, ONLY: molecule_kind_type
29 : USE molecule_types, ONLY: global_constraint_type,&
30 : molecule_type
31 : USE simpar_types, ONLY: simpar_type
32 : USE thermostat_mapping, ONLY: adiabatic_mapping_region,&
33 : init_baro_map_info,&
34 : thermostat_mapping_region
35 : USE thermostat_types, ONLY: thermostat_info_type
36 : #include "../../base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_mapping'
43 :
44 : PUBLIC :: nhc_to_particle_mapping, nhc_to_barostat_mapping, &
45 : nhc_to_shell_mapping, nhc_to_particle_mapping_fast, &
46 : nhc_to_particle_mapping_slow
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Creates the thermostatting for the barostat
52 : !> \param simpar ...
53 : !> \param nhc ...
54 : !> \par History
55 : !> CJM, 20-Feb-01 : nhc structure allocated to zero when not in use
56 : !> JGH (10-Mar-2001) : set nhc variables to zero when not in use
57 : !> \author CJM
58 : ! **************************************************************************************************
59 120 : SUBROUTINE nhc_to_barostat_mapping(simpar, nhc)
60 :
61 : TYPE(simpar_type), POINTER :: simpar
62 : TYPE(lnhc_parameters_type), POINTER :: nhc
63 :
64 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_barostat_mapping'
65 :
66 : INTEGER :: handle, i, number
67 : TYPE(map_info_type), POINTER :: map_info
68 :
69 120 : CALL timeset(routineN, handle)
70 :
71 120 : SELECT CASE (simpar%ensemble)
72 : CASE DEFAULT
73 0 : CPABORT('Never reach this point!')
74 : CASE (npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
75 120 : map_info => nhc%map_info
76 120 : map_info%dis_type = do_thermo_only_master
77 :
78 : ! Counting the total number of thermostats ( 1 for NPT_I, NPT_IA, and NPT_F )
79 120 : nhc%loc_num_nhc = 1
80 120 : nhc%glob_num_nhc = 1
81 120 : IF (simpar%ensemble == npt_f_ensemble) THEN
82 42 : number = 9
83 : ELSE
84 78 : number = 1
85 : END IF
86 :
87 120 : CALL init_baro_map_info(map_info, number, nhc%loc_num_nhc)
88 :
89 968 : ALLOCATE (nhc%nvt(nhc%nhc_len, nhc%loc_num_nhc))
90 : ! Now that we know how many there are stick this into nhc % nkt
91 : ! (number of degrees of freedom times k_B T )
92 240 : DO i = 1, nhc%loc_num_nhc
93 120 : nhc%nvt(1, i)%nkt = simpar%temp_ext*number
94 120 : nhc%nvt(1, i)%degrees_of_freedom = number
95 120 : IF (debug_isotropic_limit) THEN
96 : nhc%nvt(1, i)%nkt = simpar%temp_ext
97 : END IF
98 : END DO
99 :
100 : ! getting the number of degrees of freedom times k_B T for the rest of the chain
101 368 : DO i = 2, nhc%nhc_len
102 616 : nhc%nvt(i, :)%nkt = simpar%temp_ext
103 : END DO
104 :
105 : ! Let's clean the arrays
106 240 : map_info%s_kin = 0.0_dp
107 360 : map_info%v_scale = 0.0_dp
108 : END SELECT
109 :
110 120 : CALL timestop(handle)
111 :
112 120 : END SUBROUTINE nhc_to_barostat_mapping
113 :
114 : ! **************************************************************************************************
115 : !> \brief Creates the thermostatting maps
116 : !> \param thermostat_info ...
117 : !> \param simpar ...
118 : !> \param local_molecules ...
119 : !> \param molecule_set ...
120 : !> \param molecule_kind_set ...
121 : !> \param nhc ...
122 : !> \param para_env ...
123 : !> \param gci ...
124 : !> \par History
125 : !> 29-Nov-00 (JGH) correct counting of DOF if constraints are off
126 : !> CJM, 20-Feb-01 : nhc structure allocated to zero when not in use
127 : !> JGH (10-Mar-2001) : set nhc variables to zero when not in use
128 : !> CJM(10-NOV-2001) : New parallelization with new molecule structures
129 : !> Teodoro Laino 09.2007 [tlaino] - University of Zurich - cleaning and updating
130 : !> \author CJM
131 : ! **************************************************************************************************
132 376 : SUBROUTINE nhc_to_particle_mapping(thermostat_info, simpar, local_molecules, &
133 : molecule_set, molecule_kind_set, nhc, para_env, gci)
134 :
135 : TYPE(thermostat_info_type), POINTER :: thermostat_info
136 : TYPE(simpar_type), POINTER :: simpar
137 : TYPE(distribution_1d_type), POINTER :: local_molecules
138 : TYPE(molecule_type), POINTER :: molecule_set(:)
139 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
140 : TYPE(lnhc_parameters_type), POINTER :: nhc
141 : TYPE(mp_para_env_type), POINTER :: para_env
142 : TYPE(global_constraint_type), POINTER :: gci
143 :
144 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping'
145 :
146 : INTEGER :: handle, i, imap, j, natoms_local, &
147 : sum_of_thermostats
148 376 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
149 : REAL(KIND=dp) :: fac
150 : TYPE(map_info_type), POINTER :: map_info
151 :
152 376 : CALL timeset(routineN, handle)
153 :
154 376 : NULLIFY (massive_atom_list, deg_of_freedom)
155 :
156 376 : SELECT CASE (simpar%ensemble)
157 : CASE DEFAULT
158 0 : CPABORT('Unknown ensemble!')
159 : CASE (nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble, &
160 : nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
161 0 : CPABORT('Never reach this point!')
162 : CASE (nvt_ensemble, npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
163 :
164 : CALL setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, &
165 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
166 376 : simpar, sum_of_thermostats, gci)
167 :
168 : ! Sum up the number of degrees of freedom on each thermostat.
169 : ! first: initialize the target
170 376 : map_info => nhc%map_info
171 21375 : map_info%s_kin = 0.0_dp
172 1504 : DO i = 1, 3
173 109957 : DO j = 1, natoms_local
174 109581 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
175 : END DO
176 : END DO
177 :
178 : ! if thermostats are replicated but molecules distributed, we have to
179 : ! sum s_kin over all processors
180 800 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
181 :
182 : ! We know the total number of system thermostats.
183 376 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
184 184 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
185 184 : IF (fac == 0.0_dp) THEN
186 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
187 : END IF
188 184 : nhc%nvt(1, 1)%nkt = simpar%temp_ext*fac
189 184 : nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac)
190 : ELSE
191 21003 : DO i = 1, nhc%loc_num_nhc
192 20811 : imap = map_info%map_index(i)
193 20811 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
194 20811 : nhc%nvt(1, i)%nkt = simpar%temp_ext*fac
195 21003 : nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac)
196 : END DO
197 : END IF
198 :
199 : ! Getting the number of degrees of freedom times k_B T for the rest
200 : ! of the chain
201 1222 : DO i = 2, nhc%nhc_len
202 42008 : nhc%nvt(i, :)%nkt = simpar%temp_ext
203 42384 : nhc%nvt(i, :)%degrees_of_freedom = 1
204 : END DO
205 376 : DEALLOCATE (deg_of_freedom)
206 376 : DEALLOCATE (massive_atom_list)
207 :
208 : ! Let's clean the arrays
209 21375 : map_info%s_kin = 0.0_dp
210 21751 : map_info%v_scale = 0.0_dp
211 : END SELECT
212 :
213 376 : CALL timestop(handle)
214 :
215 376 : END SUBROUTINE nhc_to_particle_mapping
216 :
217 : ! **************************************************************************************************
218 : !> \brief Main general setup for Adiabatic Nose-Hoover thermostats
219 : !> \param nhc ...
220 : !> \param thermostat_info ...
221 : !> \param deg_of_freedom ...
222 : !> \param massive_atom_list ...
223 : !> \param molecule_kind_set ...
224 : !> \param local_molecules ...
225 : !> \param molecule_set ...
226 : !> \param para_env ...
227 : !> \param natoms_local ...
228 : !> \param simpar ...
229 : !> \param sum_of_thermostats ...
230 : !> \param gci ...
231 : !> \param shell ...
232 : !> \author CJM -PNNL -2011
233 : ! **************************************************************************************************
234 0 : SUBROUTINE setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, &
235 : massive_atom_list, molecule_kind_set, local_molecules, molecule_set, &
236 : para_env, natoms_local, simpar, sum_of_thermostats, gci, shell)
237 :
238 : TYPE(lnhc_parameters_type), POINTER :: nhc
239 : TYPE(thermostat_info_type), POINTER :: thermostat_info
240 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
241 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
242 : TYPE(distribution_1d_type), POINTER :: local_molecules
243 : TYPE(molecule_type), POINTER :: molecule_set(:)
244 : TYPE(mp_para_env_type), POINTER :: para_env
245 : INTEGER, INTENT(OUT) :: natoms_local
246 : TYPE(simpar_type), POINTER :: simpar
247 : INTEGER, INTENT(OUT) :: sum_of_thermostats
248 : TYPE(global_constraint_type), POINTER :: gci
249 : LOGICAL, INTENT(IN), OPTIONAL :: shell
250 :
251 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_adiabatic_thermostat'
252 :
253 : INTEGER :: handle, nkind, number, region
254 : LOGICAL :: do_shell
255 : TYPE(map_info_type), POINTER :: map_info
256 :
257 0 : CALL timeset(routineN, handle)
258 :
259 0 : do_shell = .FALSE.
260 0 : IF (PRESENT(shell)) do_shell = shell
261 0 : map_info => nhc%map_info
262 :
263 0 : nkind = SIZE(molecule_kind_set)
264 0 : sum_of_thermostats = thermostat_info%sum_of_thermostats
265 0 : map_info%dis_type = thermostat_info%dis_type
266 0 : number = thermostat_info%number_of_thermostats
267 0 : region = nhc%region
268 :
269 : CALL adiabatic_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
270 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
271 : simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, &
272 0 : sum_of_thermostats)
273 0 : ALLOCATE (nhc%nvt(nhc%nhc_len, number))
274 :
275 : ! Now that we know how many there are stick this into nhc%nkt
276 : ! (number of degrees of freedom times k_B T for the first thermostat
277 : ! on the chain)
278 0 : nhc%loc_num_nhc = number
279 0 : nhc%glob_num_nhc = sum_of_thermostats
280 :
281 0 : CALL timestop(handle)
282 :
283 0 : END SUBROUTINE setup_adiabatic_thermostat
284 :
285 : ! **************************************************************************************************
286 : !> \brief Creates the thermostatting maps
287 : !> \param thermostat_info ...
288 : !> \param simpar ...
289 : !> \param local_molecules ...
290 : !> \param molecule_set ...
291 : !> \param molecule_kind_set ...
292 : !> \param nhc ...
293 : !> \param para_env ...
294 : !> \param gci ...
295 : !> \par History
296 : !> \author CJM
297 : ! **************************************************************************************************
298 0 : SUBROUTINE nhc_to_particle_mapping_slow(thermostat_info, simpar, local_molecules, &
299 : molecule_set, molecule_kind_set, nhc, para_env, gci)
300 :
301 : TYPE(thermostat_info_type), POINTER :: thermostat_info
302 : TYPE(simpar_type), POINTER :: simpar
303 : TYPE(distribution_1d_type), POINTER :: local_molecules
304 : TYPE(molecule_type), POINTER :: molecule_set(:)
305 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
306 : TYPE(lnhc_parameters_type), POINTER :: nhc
307 : TYPE(mp_para_env_type), POINTER :: para_env
308 : TYPE(global_constraint_type), POINTER :: gci
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping_slow'
311 :
312 : INTEGER :: handle, i, imap, j, natoms_local, &
313 : sum_of_thermostats
314 0 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
315 : REAL(KIND=dp) :: fac
316 : TYPE(map_info_type), POINTER :: map_info
317 :
318 0 : CALL timeset(routineN, handle)
319 :
320 0 : NULLIFY (massive_atom_list, deg_of_freedom)
321 :
322 0 : SELECT CASE (simpar%ensemble)
323 : CASE DEFAULT
324 0 : CPABORT('Unknown ensemble!')
325 : CASE (nvt_adiabatic_ensemble)
326 : CALL setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, &
327 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
328 0 : simpar, sum_of_thermostats, gci)
329 :
330 : ! Sum up the number of degrees of freedom on each thermostat.
331 : ! first: initialize the target
332 0 : map_info => nhc%map_info
333 0 : map_info%s_kin = 0.0_dp
334 0 : DO i = 1, 3
335 0 : DO j = 1, natoms_local
336 0 : IF (ASSOCIATED(map_info%p_kin(i, j)%point)) THEN
337 0 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
338 : END IF
339 : END DO
340 : END DO
341 :
342 : ! if thermostats are replicated but molecules distributed, we have to
343 : ! sum s_kin over all processors
344 0 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
345 :
346 : ! We know the total number of system thermostats.
347 0 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
348 0 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
349 0 : IF (fac == 0.0_dp) THEN
350 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
351 : END IF
352 0 : nhc%nvt(1, 1)%nkt = simpar%temp_slow*fac
353 0 : nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac)
354 : ELSE
355 0 : DO i = 1, nhc%loc_num_nhc
356 0 : imap = map_info%map_index(i)
357 0 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
358 0 : nhc%nvt(1, i)%nkt = simpar%temp_slow*fac
359 0 : nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac)
360 : END DO
361 : END IF
362 :
363 : ! Getting the number of degrees of freedom times k_B T for the rest
364 : ! of the chain
365 0 : DO i = 2, nhc%nhc_len
366 0 : nhc%nvt(i, :)%nkt = simpar%temp_slow
367 0 : nhc%nvt(i, :)%degrees_of_freedom = 1
368 : END DO
369 0 : DEALLOCATE (deg_of_freedom)
370 0 : DEALLOCATE (massive_atom_list)
371 :
372 : ! Let's clean the arrays
373 0 : map_info%s_kin = 0.0_dp
374 0 : map_info%v_scale = 0.0_dp
375 : END SELECT
376 :
377 0 : CALL timestop(handle)
378 :
379 0 : END SUBROUTINE nhc_to_particle_mapping_slow
380 :
381 : ! **************************************************************************************************
382 : !> \brief Creates the thermostatting maps
383 : !> \param thermostat_info ...
384 : !> \param simpar ...
385 : !> \param local_molecules ...
386 : !> \param molecule_set ...
387 : !> \param molecule_kind_set ...
388 : !> \param nhc ...
389 : !> \param para_env ...
390 : !> \param gci ...
391 : !> \par History
392 : !> \author CJM
393 : ! **************************************************************************************************
394 0 : SUBROUTINE nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, &
395 : molecule_set, molecule_kind_set, nhc, para_env, gci)
396 :
397 : TYPE(thermostat_info_type), POINTER :: thermostat_info
398 : TYPE(simpar_type), POINTER :: simpar
399 : TYPE(distribution_1d_type), POINTER :: local_molecules
400 : TYPE(molecule_type), POINTER :: molecule_set(:)
401 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
402 : TYPE(lnhc_parameters_type), POINTER :: nhc
403 : TYPE(mp_para_env_type), POINTER :: para_env
404 : TYPE(global_constraint_type), POINTER :: gci
405 :
406 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping_fast'
407 :
408 : INTEGER :: handle, i, imap, j, natoms_local, &
409 : sum_of_thermostats
410 0 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
411 : REAL(KIND=dp) :: fac
412 : TYPE(map_info_type), POINTER :: map_info
413 :
414 0 : CALL timeset(routineN, handle)
415 :
416 0 : NULLIFY (massive_atom_list, deg_of_freedom)
417 :
418 0 : SELECT CASE (simpar%ensemble)
419 : CASE DEFAULT
420 0 : CPABORT('Unknown ensemble!')
421 : CASE (nvt_adiabatic_ensemble)
422 : CALL setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, &
423 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
424 0 : simpar, sum_of_thermostats, gci)
425 :
426 : ! Sum up the number of degrees of freedom on each thermostat.
427 : ! first: initialize the target
428 0 : map_info => nhc%map_info
429 0 : map_info%s_kin = 0.0_dp
430 0 : DO i = 1, 3
431 0 : DO j = 1, natoms_local
432 0 : IF (ASSOCIATED(map_info%p_kin(i, j)%point)) THEN
433 0 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
434 : END IF
435 : END DO
436 : END DO
437 :
438 : ! if thermostats are replicated but molecules distributed, we have to
439 : ! sum s_kin over all processors
440 0 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
441 :
442 : ! We know the total number of system thermostats.
443 0 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
444 0 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
445 0 : IF (fac == 0.0_dp) THEN
446 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
447 : END IF
448 0 : nhc%nvt(1, 1)%nkt = simpar%temp_fast*fac
449 0 : nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac)
450 : ELSE
451 0 : DO i = 1, nhc%loc_num_nhc
452 0 : imap = map_info%map_index(i)
453 0 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
454 0 : nhc%nvt(1, i)%nkt = simpar%temp_fast*fac
455 0 : nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac)
456 : END DO
457 : END IF
458 :
459 : ! Getting the number of degrees of freedom times k_B T for the rest
460 : ! of the chain
461 0 : DO i = 2, nhc%nhc_len
462 0 : nhc%nvt(i, :)%nkt = simpar%temp_fast
463 0 : nhc%nvt(i, :)%degrees_of_freedom = 1
464 : END DO
465 0 : DEALLOCATE (deg_of_freedom)
466 0 : DEALLOCATE (massive_atom_list)
467 :
468 : ! Let's clean the arrays
469 0 : map_info%s_kin = 0.0_dp
470 0 : map_info%v_scale = 0.0_dp
471 : END SELECT
472 :
473 0 : CALL timestop(handle)
474 :
475 0 : END SUBROUTINE nhc_to_particle_mapping_fast
476 :
477 : ! **************************************************************************************************
478 : !> \brief Main general setup for Nose-Hoover thermostats
479 : !> \param nhc ...
480 : !> \param thermostat_info ...
481 : !> \param deg_of_freedom ...
482 : !> \param massive_atom_list ...
483 : !> \param molecule_kind_set ...
484 : !> \param local_molecules ...
485 : !> \param molecule_set ...
486 : !> \param para_env ...
487 : !> \param natoms_local ...
488 : !> \param simpar ...
489 : !> \param sum_of_thermostats ...
490 : !> \param gci ...
491 : !> \param shell ...
492 : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
493 : ! **************************************************************************************************
494 416 : SUBROUTINE setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, &
495 : massive_atom_list, molecule_kind_set, local_molecules, molecule_set, &
496 : para_env, natoms_local, simpar, sum_of_thermostats, gci, shell)
497 :
498 : TYPE(lnhc_parameters_type), POINTER :: nhc
499 : TYPE(thermostat_info_type), POINTER :: thermostat_info
500 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
501 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
502 : TYPE(distribution_1d_type), POINTER :: local_molecules
503 : TYPE(molecule_type), POINTER :: molecule_set(:)
504 : TYPE(mp_para_env_type), POINTER :: para_env
505 : INTEGER, INTENT(OUT) :: natoms_local
506 : TYPE(simpar_type), POINTER :: simpar
507 : INTEGER, INTENT(OUT) :: sum_of_thermostats
508 : TYPE(global_constraint_type), POINTER :: gci
509 : LOGICAL, INTENT(IN), OPTIONAL :: shell
510 :
511 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_nhc_thermostat'
512 :
513 : INTEGER :: handle, nkind, number, region
514 : LOGICAL :: do_shell
515 : TYPE(map_info_type), POINTER :: map_info
516 :
517 416 : CALL timeset(routineN, handle)
518 :
519 416 : do_shell = .FALSE.
520 416 : IF (PRESENT(shell)) do_shell = shell
521 416 : map_info => nhc%map_info
522 :
523 416 : nkind = SIZE(molecule_kind_set)
524 416 : sum_of_thermostats = thermostat_info%sum_of_thermostats
525 416 : map_info%dis_type = thermostat_info%dis_type
526 416 : number = thermostat_info%number_of_thermostats
527 416 : region = nhc%region
528 :
529 : CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
530 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
531 : simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, &
532 416 : sum_of_thermostats)
533 :
534 109462 : ALLOCATE (nhc%nvt(nhc%nhc_len, number))
535 :
536 : ! Now that we know how many there are stick this into nhc%nkt
537 : ! (number of degrees of freedom times k_B T for the first thermostat
538 : ! on the chain)
539 416 : nhc%loc_num_nhc = number
540 416 : nhc%glob_num_nhc = sum_of_thermostats
541 :
542 416 : CALL timestop(handle)
543 :
544 832 : END SUBROUTINE setup_nhc_thermostat
545 :
546 : ! **************************************************************************************************
547 : !> \brief ...
548 : !> \param thermostat_info ...
549 : !> \param simpar ...
550 : !> \param local_molecules ...
551 : !> \param molecule_set ...
552 : !> \param molecule_kind_set ...
553 : !> \param nhc ...
554 : !> \param para_env ...
555 : !> \param gci ...
556 : ! **************************************************************************************************
557 40 : SUBROUTINE nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, &
558 : molecule_set, molecule_kind_set, nhc, para_env, gci)
559 :
560 : TYPE(thermostat_info_type), POINTER :: thermostat_info
561 : TYPE(simpar_type), POINTER :: simpar
562 : TYPE(distribution_1d_type), POINTER :: local_molecules
563 : TYPE(molecule_type), POINTER :: molecule_set(:)
564 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
565 : TYPE(lnhc_parameters_type), POINTER :: nhc
566 : TYPE(mp_para_env_type), POINTER :: para_env
567 : TYPE(global_constraint_type), POINTER :: gci
568 :
569 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_shell_mapping'
570 :
571 : INTEGER :: handle, i, imap, j, nshell_local, &
572 : sum_of_thermostats
573 40 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_shell_list
574 : TYPE(map_info_type), POINTER :: map_info
575 :
576 40 : CALL timeset(routineN, handle)
577 :
578 40 : NULLIFY (massive_shell_list, deg_of_freedom)
579 :
580 40 : SELECT CASE (simpar%ensemble)
581 : CASE DEFAULT
582 0 : CPABORT('Unknown ensemble!')
583 : CASE (isokin_ensemble, nph_uniaxial_ensemble, &
584 : nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
585 0 : CPABORT('Never reach this point!')
586 : CASE (nve_ensemble, nvt_ensemble, npe_f_ensemble, npe_i_ensemble, npt_i_ensemble, npt_f_ensemble, &
587 : npt_ia_ensemble)
588 :
589 : CALL setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, massive_shell_list, &
590 : molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, &
591 40 : simpar, sum_of_thermostats, gci, shell=.TRUE.)
592 :
593 40 : map_info => nhc%map_info
594 : ! Sum up the number of degrees of freedom on each thermostat.
595 : ! first: initialize the target, via p_kin init s_kin
596 4178 : map_info%s_kin = 0.0_dp
597 1960 : DO j = 1, nshell_local
598 7720 : DO i = 1, 3
599 7680 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
600 : END DO
601 : END DO
602 :
603 : ! If thermostats are replicated but molecules distributed, we have to
604 : ! sum s_kin over all processors
605 60 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
606 :
607 : ! Now that we know how many there are stick this into nhc%nkt
608 : ! (number of degrees of freedom times k_B T )
609 4178 : DO i = 1, nhc%loc_num_nhc
610 4138 : imap = map_info%map_index(i)
611 4138 : nhc%nvt(1, i)%nkt = simpar%temp_sh_ext*map_info%s_kin(imap)
612 4178 : nhc%nvt(1, i)%degrees_of_freedom = INT(map_info%s_kin(imap))
613 : END DO
614 :
615 : ! Getting the number of degrees of freedom times k_B T for the rest of the chain
616 210 : DO i = 2, nhc%nhc_len
617 16540 : nhc%nvt(i, :)%nkt = simpar%temp_sh_ext
618 16580 : nhc%nvt(i, :)%degrees_of_freedom = 1
619 : END DO
620 40 : DEALLOCATE (deg_of_freedom)
621 40 : DEALLOCATE (massive_shell_list)
622 :
623 : ! Let's clean the arrays
624 4178 : map_info%s_kin = 0.0_dp
625 4218 : map_info%v_scale = 0.0_dp
626 : END SELECT
627 :
628 40 : CALL timestop(handle)
629 :
630 40 : END SUBROUTINE nhc_to_shell_mapping
631 :
632 : END MODULE extended_system_mapping
|
