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 : !> \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 396 : 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 396 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
149 : REAL(KIND=dp) :: fac
150 : TYPE(map_info_type), POINTER :: map_info
151 :
152 396 : CALL timeset(routineN, handle)
153 :
154 396 : NULLIFY (massive_atom_list, deg_of_freedom)
155 :
156 396 : 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 396 : simpar, sum_of_thermostats, gci)
167 :
168 : ! Sum up the number of degrees of freedom on each thermostat.
169 : ! first: initialize the target
170 396 : map_info => nhc%map_info
171 21415 : map_info%s_kin = 0.0_dp
172 1584 : DO i = 1, 3
173 112137 : DO j = 1, natoms_local
174 111741 : 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 860 : 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 396 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
184 204 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
185 204 : IF (fac == 0.0_dp) THEN
186 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
187 : END IF
188 204 : nhc%nvt(1, 1)%nkt = simpar%temp_ext*fac
189 204 : 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 1282 : DO i = 2, nhc%nhc_len
202 42088 : nhc%nvt(i, :)%nkt = simpar%temp_ext
203 42484 : nhc%nvt(i, :)%degrees_of_freedom = 1
204 : END DO
205 396 : DEALLOCATE (deg_of_freedom)
206 396 : DEALLOCATE (massive_atom_list)
207 :
208 : ! Let's clean the arrays
209 21415 : map_info%s_kin = 0.0_dp
210 21811 : map_info%v_scale = 0.0_dp
211 : END SELECT
212 :
213 396 : CALL timestop(handle)
214 :
215 396 : 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)) &
337 0 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
338 : END DO
339 : END DO
340 :
341 : ! if thermostats are replicated but molecules distributed, we have to
342 : ! sum s_kin over all processors
343 0 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
344 :
345 : ! We know the total number of system thermostats.
346 0 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
347 0 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
348 0 : IF (fac == 0.0_dp) THEN
349 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
350 : END IF
351 0 : nhc%nvt(1, 1)%nkt = simpar%temp_slow*fac
352 0 : nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac)
353 : ELSE
354 0 : DO i = 1, nhc%loc_num_nhc
355 0 : imap = map_info%map_index(i)
356 0 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
357 0 : nhc%nvt(1, i)%nkt = simpar%temp_slow*fac
358 0 : nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac)
359 : END DO
360 : END IF
361 :
362 : ! Getting the number of degrees of freedom times k_B T for the rest
363 : ! of the chain
364 0 : DO i = 2, nhc%nhc_len
365 0 : nhc%nvt(i, :)%nkt = simpar%temp_slow
366 0 : nhc%nvt(i, :)%degrees_of_freedom = 1
367 : END DO
368 0 : DEALLOCATE (deg_of_freedom)
369 0 : DEALLOCATE (massive_atom_list)
370 :
371 : ! Let's clean the arrays
372 0 : map_info%s_kin = 0.0_dp
373 0 : map_info%v_scale = 0.0_dp
374 : END SELECT
375 :
376 0 : CALL timestop(handle)
377 :
378 0 : END SUBROUTINE nhc_to_particle_mapping_slow
379 :
380 : ! **************************************************************************************************
381 : !> \brief Creates the thermostatting maps
382 : !> \param thermostat_info ...
383 : !> \param simpar ...
384 : !> \param local_molecules ...
385 : !> \param molecule_set ...
386 : !> \param molecule_kind_set ...
387 : !> \param nhc ...
388 : !> \param para_env ...
389 : !> \param gci ...
390 : !> \par History
391 : !> \author CJM
392 : ! **************************************************************************************************
393 0 : SUBROUTINE nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, &
394 : molecule_set, molecule_kind_set, nhc, para_env, gci)
395 :
396 : TYPE(thermostat_info_type), POINTER :: thermostat_info
397 : TYPE(simpar_type), POINTER :: simpar
398 : TYPE(distribution_1d_type), POINTER :: local_molecules
399 : TYPE(molecule_type), POINTER :: molecule_set(:)
400 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
401 : TYPE(lnhc_parameters_type), POINTER :: nhc
402 : TYPE(mp_para_env_type), POINTER :: para_env
403 : TYPE(global_constraint_type), POINTER :: gci
404 :
405 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping_fast'
406 :
407 : INTEGER :: handle, i, imap, j, natoms_local, &
408 : sum_of_thermostats
409 0 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
410 : REAL(KIND=dp) :: fac
411 : TYPE(map_info_type), POINTER :: map_info
412 :
413 0 : CALL timeset(routineN, handle)
414 :
415 0 : NULLIFY (massive_atom_list, deg_of_freedom)
416 :
417 0 : SELECT CASE (simpar%ensemble)
418 : CASE DEFAULT
419 0 : CPABORT('Unknown ensemble!')
420 : CASE (nvt_adiabatic_ensemble)
421 : CALL setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, &
422 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
423 0 : simpar, sum_of_thermostats, gci)
424 :
425 : ! Sum up the number of degrees of freedom on each thermostat.
426 : ! first: initialize the target
427 0 : map_info => nhc%map_info
428 0 : map_info%s_kin = 0.0_dp
429 0 : DO i = 1, 3
430 0 : DO j = 1, natoms_local
431 0 : IF (ASSOCIATED(map_info%p_kin(i, j)%point)) &
432 0 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
433 : END DO
434 : END DO
435 :
436 : ! if thermostats are replicated but molecules distributed, we have to
437 : ! sum s_kin over all processors
438 0 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
439 :
440 : ! We know the total number of system thermostats.
441 0 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
442 0 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
443 0 : IF (fac == 0.0_dp) THEN
444 0 : CPABORT('Zero degrees of freedom. Nothing to thermalize!')
445 : END IF
446 0 : nhc%nvt(1, 1)%nkt = simpar%temp_fast*fac
447 0 : nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac)
448 : ELSE
449 0 : DO i = 1, nhc%loc_num_nhc
450 0 : imap = map_info%map_index(i)
451 0 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
452 0 : nhc%nvt(1, i)%nkt = simpar%temp_fast*fac
453 0 : nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac)
454 : END DO
455 : END IF
456 :
457 : ! Getting the number of degrees of freedom times k_B T for the rest
458 : ! of the chain
459 0 : DO i = 2, nhc%nhc_len
460 0 : nhc%nvt(i, :)%nkt = simpar%temp_fast
461 0 : nhc%nvt(i, :)%degrees_of_freedom = 1
462 : END DO
463 0 : DEALLOCATE (deg_of_freedom)
464 0 : DEALLOCATE (massive_atom_list)
465 :
466 : ! Let's clean the arrays
467 0 : map_info%s_kin = 0.0_dp
468 0 : map_info%v_scale = 0.0_dp
469 : END SELECT
470 :
471 0 : CALL timestop(handle)
472 :
473 0 : END SUBROUTINE nhc_to_particle_mapping_fast
474 :
475 : ! **************************************************************************************************
476 : !> \brief Main general setup for Nose-Hoover thermostats
477 : !> \param nhc ...
478 : !> \param thermostat_info ...
479 : !> \param deg_of_freedom ...
480 : !> \param massive_atom_list ...
481 : !> \param molecule_kind_set ...
482 : !> \param local_molecules ...
483 : !> \param molecule_set ...
484 : !> \param para_env ...
485 : !> \param natoms_local ...
486 : !> \param simpar ...
487 : !> \param sum_of_thermostats ...
488 : !> \param gci ...
489 : !> \param shell ...
490 : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
491 : ! **************************************************************************************************
492 436 : SUBROUTINE setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, &
493 : massive_atom_list, molecule_kind_set, local_molecules, molecule_set, &
494 : para_env, natoms_local, simpar, sum_of_thermostats, gci, shell)
495 :
496 : TYPE(lnhc_parameters_type), POINTER :: nhc
497 : TYPE(thermostat_info_type), POINTER :: thermostat_info
498 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
499 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
500 : TYPE(distribution_1d_type), POINTER :: local_molecules
501 : TYPE(molecule_type), POINTER :: molecule_set(:)
502 : TYPE(mp_para_env_type), POINTER :: para_env
503 : INTEGER, INTENT(OUT) :: natoms_local
504 : TYPE(simpar_type), POINTER :: simpar
505 : INTEGER, INTENT(OUT) :: sum_of_thermostats
506 : TYPE(global_constraint_type), POINTER :: gci
507 : LOGICAL, INTENT(IN), OPTIONAL :: shell
508 :
509 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_nhc_thermostat'
510 :
511 : INTEGER :: handle, nkind, number, region
512 : LOGICAL :: do_shell
513 : TYPE(map_info_type), POINTER :: map_info
514 :
515 436 : CALL timeset(routineN, handle)
516 :
517 436 : do_shell = .FALSE.
518 436 : IF (PRESENT(shell)) do_shell = shell
519 436 : map_info => nhc%map_info
520 :
521 436 : nkind = SIZE(molecule_kind_set)
522 436 : sum_of_thermostats = thermostat_info%sum_of_thermostats
523 436 : map_info%dis_type = thermostat_info%dis_type
524 436 : number = thermostat_info%number_of_thermostats
525 436 : region = nhc%region
526 :
527 : CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
528 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
529 : simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, &
530 436 : sum_of_thermostats)
531 :
532 109622 : ALLOCATE (nhc%nvt(nhc%nhc_len, number))
533 :
534 : ! Now that we know how many there are stick this into nhc%nkt
535 : ! (number of degrees of freedom times k_B T for the first thermostat
536 : ! on the chain)
537 436 : nhc%loc_num_nhc = number
538 436 : nhc%glob_num_nhc = sum_of_thermostats
539 :
540 436 : CALL timestop(handle)
541 :
542 872 : END SUBROUTINE setup_nhc_thermostat
543 :
544 : ! **************************************************************************************************
545 : !> \brief ...
546 : !> \param thermostat_info ...
547 : !> \param simpar ...
548 : !> \param local_molecules ...
549 : !> \param molecule_set ...
550 : !> \param molecule_kind_set ...
551 : !> \param nhc ...
552 : !> \param para_env ...
553 : !> \param gci ...
554 : ! **************************************************************************************************
555 40 : SUBROUTINE nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, &
556 : molecule_set, molecule_kind_set, nhc, para_env, gci)
557 :
558 : TYPE(thermostat_info_type), POINTER :: thermostat_info
559 : TYPE(simpar_type), POINTER :: simpar
560 : TYPE(distribution_1d_type), POINTER :: local_molecules
561 : TYPE(molecule_type), POINTER :: molecule_set(:)
562 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
563 : TYPE(lnhc_parameters_type), POINTER :: nhc
564 : TYPE(mp_para_env_type), POINTER :: para_env
565 : TYPE(global_constraint_type), POINTER :: gci
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_shell_mapping'
568 :
569 : INTEGER :: handle, i, imap, j, nshell_local, &
570 : sum_of_thermostats
571 40 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_shell_list
572 : TYPE(map_info_type), POINTER :: map_info
573 :
574 40 : CALL timeset(routineN, handle)
575 :
576 40 : NULLIFY (massive_shell_list, deg_of_freedom)
577 :
578 40 : SELECT CASE (simpar%ensemble)
579 : CASE DEFAULT
580 0 : CPABORT('Unknown ensemble!')
581 : CASE (isokin_ensemble, nph_uniaxial_ensemble, &
582 : nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
583 0 : CPABORT('Never reach this point!')
584 : CASE (nve_ensemble, nvt_ensemble, npe_f_ensemble, npe_i_ensemble, npt_i_ensemble, npt_f_ensemble, &
585 : npt_ia_ensemble)
586 :
587 : CALL setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, massive_shell_list, &
588 : molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, &
589 40 : simpar, sum_of_thermostats, gci, shell=.TRUE.)
590 :
591 40 : map_info => nhc%map_info
592 : ! Sum up the number of degrees of freedom on each thermostat.
593 : ! first: initialize the target, via p_kin init s_kin
594 4178 : map_info%s_kin = 0.0_dp
595 1960 : DO j = 1, nshell_local
596 7720 : DO i = 1, 3
597 7680 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
598 : END DO
599 : END DO
600 :
601 : ! If thermostats are replicated but molecules distributed, we have to
602 : ! sum s_kin over all processors
603 60 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
604 :
605 : ! Now that we know how many there are stick this into nhc%nkt
606 : ! (number of degrees of freedom times k_B T )
607 4178 : DO i = 1, nhc%loc_num_nhc
608 4138 : imap = map_info%map_index(i)
609 4138 : nhc%nvt(1, i)%nkt = simpar%temp_sh_ext*map_info%s_kin(imap)
610 4178 : nhc%nvt(1, i)%degrees_of_freedom = INT(map_info%s_kin(imap))
611 : END DO
612 :
613 : ! Getting the number of degrees of freedom times k_B T for the rest of the chain
614 210 : DO i = 2, nhc%nhc_len
615 16540 : nhc%nvt(i, :)%nkt = simpar%temp_sh_ext
616 16580 : nhc%nvt(i, :)%degrees_of_freedom = 1
617 : END DO
618 40 : DEALLOCATE (deg_of_freedom)
619 40 : DEALLOCATE (massive_shell_list)
620 :
621 : ! Let's clean the arrays
622 4178 : map_info%s_kin = 0.0_dp
623 4218 : map_info%v_scale = 0.0_dp
624 : END SELECT
625 :
626 40 : CALL timestop(handle)
627 :
628 40 : END SUBROUTINE nhc_to_shell_mapping
629 :
630 : END MODULE extended_system_mapping
|
