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, FEB 20 2001: added subroutine initialize_extended_parameters
11 : !> cjm, MAY 03 2001: reorganized and added separtate routines for
12 : !> nhc_part, nhc_baro, nhc_ao, npt
13 : !> \author CJM
14 : ! **************************************************************************************************
15 : MODULE extended_system_init
16 :
17 : USE cell_types, ONLY: cell_type
18 : USE distribution_1d_types, ONLY: distribution_1d_type
19 : USE extended_system_mapping, ONLY: nhc_to_barostat_mapping,&
20 : nhc_to_particle_mapping,&
21 : nhc_to_particle_mapping_fast,&
22 : nhc_to_particle_mapping_slow,&
23 : nhc_to_shell_mapping
24 : USE extended_system_types, ONLY: debug_isotropic_limit,&
25 : lnhc_parameters_type,&
26 : map_info_type,&
27 : npt_info_type
28 : USE global_types, ONLY: global_environment_type
29 : USE input_constants, ONLY: do_thermo_only_master,&
30 : npe_f_ensemble,&
31 : npe_i_ensemble,&
32 : nph_uniaxial_damped_ensemble,&
33 : nph_uniaxial_ensemble,&
34 : npt_f_ensemble,&
35 : npt_i_ensemble,&
36 : npt_ia_ensemble
37 : USE input_cp2k_binary_restarts, ONLY: read_binary_thermostats_nose
38 : USE input_section_types, ONLY: section_vals_get,&
39 : section_vals_get_subs_vals,&
40 : section_vals_remove_values,&
41 : section_vals_type,&
42 : section_vals_val_get
43 : USE kinds, ONLY: dp
44 : USE message_passing, ONLY: mp_para_env_type
45 : USE molecule_kind_types, ONLY: molecule_kind_type
46 : USE molecule_types, ONLY: global_constraint_type,&
47 : molecule_type
48 : USE simpar_types, ONLY: simpar_type
49 : USE thermostat_types, ONLY: thermostat_info_type
50 : USE thermostat_utils, ONLY: get_nhc_energies
51 : #include "../../base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_init'
58 :
59 : PUBLIC :: initialize_nhc_part, initialize_nhc_baro, initialize_npt, &
60 : initialize_nhc_shell, initialize_nhc_slow, initialize_nhc_fast
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief ...
66 : !> \param simpar ...
67 : !> \param globenv ...
68 : !> \param npt_info ...
69 : !> \param cell ...
70 : !> \param work_section ...
71 : !> \author CJM
72 : ! **************************************************************************************************
73 174 : SUBROUTINE initialize_npt(simpar, globenv, npt_info, cell, work_section)
74 :
75 : TYPE(simpar_type), POINTER :: simpar
76 : TYPE(global_environment_type), POINTER :: globenv
77 : TYPE(npt_info_type), DIMENSION(:, :), POINTER :: npt_info
78 : TYPE(cell_type), POINTER :: cell
79 : TYPE(section_vals_type), POINTER :: work_section
80 :
81 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_npt'
82 :
83 : INTEGER :: handle, i, ind, j
84 : LOGICAL :: explicit, restart
85 : REAL(KIND=dp) :: temp
86 174 : REAL(KIND=dp), DIMENSION(:), POINTER :: buffer
87 : TYPE(section_vals_type), POINTER :: work_section2
88 :
89 174 : CALL timeset(routineN, handle)
90 :
91 174 : NULLIFY (work_section2)
92 :
93 : explicit = .FALSE.
94 174 : restart = .FALSE.
95 :
96 174 : CPASSERT(.NOT. ASSOCIATED(npt_info))
97 :
98 : ! first allocating the npt_info_type if requested
99 282 : SELECT CASE (simpar%ensemble)
100 : CASE (npt_i_ensemble, npe_i_ensemble, npt_ia_ensemble)
101 324 : ALLOCATE (npt_info(1, 1))
102 324 : npt_info(:, :)%eps = LOG(cell%deth)/3.0_dp
103 108 : temp = simpar%temp_baro_ext
104 :
105 : CASE (npt_f_ensemble, npe_f_ensemble)
106 780 : ALLOCATE (npt_info(3, 3))
107 60 : temp = simpar%temp_baro_ext
108 :
109 : CASE (nph_uniaxial_ensemble)
110 12 : ALLOCATE (npt_info(1, 1))
111 4 : temp = simpar%temp_baro_ext
112 :
113 : CASE (nph_uniaxial_damped_ensemble)
114 6 : ALLOCATE (npt_info(1, 1))
115 2 : temp = simpar%temp_baro_ext
116 :
117 : CASE DEFAULT
118 : ! Do nothing..
119 174 : NULLIFY (npt_info)
120 : END SELECT
121 :
122 174 : IF (ASSOCIATED(npt_info)) THEN
123 174 : IF (ASSOCIATED(work_section)) THEN
124 174 : work_section2 => section_vals_get_subs_vals(work_section, "VELOCITY")
125 174 : CALL section_vals_get(work_section2, explicit=explicit)
126 174 : restart = explicit
127 174 : work_section2 => section_vals_get_subs_vals(work_section, "MASS")
128 174 : CALL section_vals_get(work_section2, explicit=explicit)
129 174 : IF (restart .NEQV. explicit) THEN
130 : CALL cp_abort(__LOCATION__, "You need to define both VELOCITY and "// &
131 0 : "MASS section (or none) in the BAROSTAT section")
132 : END IF
133 174 : restart = explicit .AND. restart
134 : END IF
135 :
136 : IF (restart) THEN
137 22 : CALL section_vals_val_get(work_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
138 22 : ind = 0
139 72 : DO i = 1, SIZE(npt_info, 1)
140 206 : DO j = 1, SIZE(npt_info, 2)
141 134 : ind = ind + 1
142 184 : npt_info(i, j)%v = buffer(ind)
143 : END DO
144 : END DO
145 22 : CALL section_vals_val_get(work_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
146 22 : ind = 0
147 72 : DO i = 1, SIZE(npt_info, 1)
148 206 : DO j = 1, SIZE(npt_info, 2)
149 134 : ind = ind + 1
150 184 : npt_info(i, j)%mass = buffer(ind)
151 : END DO
152 : END DO
153 : ELSE
154 : CALL init_barostat_variables(npt_info, simpar%tau_cell, temp, &
155 : simpar%nfree, simpar%ensemble, simpar%cmass, &
156 152 : globenv)
157 : END IF
158 :
159 : END IF
160 :
161 174 : CALL timestop(handle)
162 :
163 174 : END SUBROUTINE initialize_npt
164 :
165 : ! **************************************************************************************************
166 : !> \brief fire up the thermostats, if NPT
167 : !> \param simpar ...
168 : !> \param para_env ...
169 : !> \param globenv ...
170 : !> \param nhc ...
171 : !> \param nose_section ...
172 : !> \param save_mem ...
173 : !> \author CJM
174 : ! **************************************************************************************************
175 240 : SUBROUTINE initialize_nhc_baro(simpar, para_env, globenv, nhc, nose_section, save_mem)
176 :
177 : TYPE(simpar_type), POINTER :: simpar
178 : TYPE(mp_para_env_type), POINTER :: para_env
179 : TYPE(global_environment_type), POINTER :: globenv
180 : TYPE(lnhc_parameters_type), POINTER :: nhc
181 : TYPE(section_vals_type), POINTER :: nose_section
182 : LOGICAL, INTENT(IN) :: save_mem
183 :
184 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_baro'
185 :
186 : INTEGER :: handle
187 : LOGICAL :: restart
188 : REAL(KIND=dp) :: temp
189 :
190 120 : CALL timeset(routineN, handle)
191 :
192 : restart = .FALSE.
193 :
194 120 : CALL nhc_to_barostat_mapping(simpar, nhc)
195 :
196 : ! Set up the Yoshida weights
197 120 : IF (nhc%nyosh > 0) THEN
198 360 : ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
199 120 : CALL set_yoshida_coef(nhc, simpar%dt)
200 : END IF
201 :
202 120 : CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
203 :
204 120 : IF (.NOT. restart) THEN
205 : ! Initializing thermostat forces and velocities for the Nose-Hoover
206 : ! Chain variables
207 102 : SELECT CASE (simpar%ensemble)
208 : CASE DEFAULT
209 102 : temp = simpar%temp_baro_ext
210 : END SELECT
211 102 : IF (nhc%nhc_len /= 0) THEN
212 102 : CALL init_nhc_variables(nhc, temp, para_env, globenv)
213 : END IF
214 : END IF
215 :
216 120 : CALL init_nhc_forces(nhc)
217 :
218 120 : CALL timestop(handle)
219 :
220 120 : END SUBROUTINE initialize_nhc_baro
221 :
222 : ! **************************************************************************************************
223 : !> \brief ...
224 : !> \param thermostat_info ...
225 : !> \param simpar ...
226 : !> \param local_molecules ...
227 : !> \param molecule ...
228 : !> \param molecule_kind_set ...
229 : !> \param para_env ...
230 : !> \param globenv ...
231 : !> \param nhc ...
232 : !> \param nose_section ...
233 : !> \param gci ...
234 : !> \param save_mem ...
235 : !> \author CJM
236 : ! **************************************************************************************************
237 0 : SUBROUTINE initialize_nhc_slow(thermostat_info, simpar, local_molecules, &
238 : molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, &
239 : gci, save_mem)
240 :
241 : TYPE(thermostat_info_type), POINTER :: thermostat_info
242 : TYPE(simpar_type), POINTER :: simpar
243 : TYPE(distribution_1d_type), POINTER :: local_molecules
244 : TYPE(molecule_type), POINTER :: molecule(:)
245 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
246 : TYPE(mp_para_env_type), POINTER :: para_env
247 : TYPE(global_environment_type), POINTER :: globenv
248 : TYPE(lnhc_parameters_type), POINTER :: nhc
249 : TYPE(section_vals_type), POINTER :: nose_section
250 : TYPE(global_constraint_type), POINTER :: gci
251 : LOGICAL, INTENT(IN) :: save_mem
252 :
253 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_slow'
254 :
255 : INTEGER :: handle
256 : LOGICAL :: restart
257 :
258 0 : CALL timeset(routineN, handle)
259 :
260 : restart = .FALSE.
261 : ! fire up the thermostats, if not NVE
262 :
263 : CALL nhc_to_particle_mapping_slow(thermostat_info, simpar, local_molecules, &
264 0 : molecule, molecule_kind_set, nhc, para_env, gci)
265 :
266 : ! Set up the Yoshida weights
267 0 : IF (nhc%nyosh > 0) THEN
268 0 : ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
269 0 : CALL set_yoshida_coef(nhc, simpar%dt)
270 : END IF
271 :
272 0 : CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
273 :
274 0 : IF (.NOT. restart) THEN
275 : ! Initializing thermostat forces and velocities for the Nose-Hoover
276 : ! Chain variables
277 0 : IF (nhc%nhc_len /= 0) THEN
278 0 : CALL init_nhc_variables(nhc, simpar%temp_slow, para_env, globenv)
279 : END IF
280 : END IF
281 :
282 0 : CALL init_nhc_forces(nhc)
283 :
284 0 : CALL timestop(handle)
285 :
286 0 : END SUBROUTINE initialize_nhc_slow
287 :
288 : ! **************************************************************************************************
289 : !> \brief ...
290 : !> \param thermostat_info ...
291 : !> \param simpar ...
292 : !> \param local_molecules ...
293 : !> \param molecule ...
294 : !> \param molecule_kind_set ...
295 : !> \param para_env ...
296 : !> \param globenv ...
297 : !> \param nhc ...
298 : !> \param nose_section ...
299 : !> \param gci ...
300 : !> \param save_mem ...
301 : !> \author CJM
302 : ! **************************************************************************************************
303 0 : SUBROUTINE initialize_nhc_fast(thermostat_info, simpar, local_molecules, &
304 : molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, &
305 : gci, save_mem)
306 :
307 : TYPE(thermostat_info_type), POINTER :: thermostat_info
308 : TYPE(simpar_type), POINTER :: simpar
309 : TYPE(distribution_1d_type), POINTER :: local_molecules
310 : TYPE(molecule_type), POINTER :: molecule(:)
311 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
312 : TYPE(mp_para_env_type), POINTER :: para_env
313 : TYPE(global_environment_type), POINTER :: globenv
314 : TYPE(lnhc_parameters_type), POINTER :: nhc
315 : TYPE(section_vals_type), POINTER :: nose_section
316 : TYPE(global_constraint_type), POINTER :: gci
317 : LOGICAL, INTENT(IN) :: save_mem
318 :
319 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_fast'
320 :
321 : INTEGER :: handle
322 : LOGICAL :: restart
323 :
324 0 : CALL timeset(routineN, handle)
325 :
326 : restart = .FALSE.
327 : ! fire up the thermostats, if not NVE
328 :
329 : CALL nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, &
330 0 : molecule, molecule_kind_set, nhc, para_env, gci)
331 :
332 : ! Set up the Yoshida weights
333 0 : IF (nhc%nyosh > 0) THEN
334 0 : ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
335 0 : CALL set_yoshida_coef(nhc, simpar%dt)
336 : END IF
337 :
338 0 : CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
339 :
340 0 : IF (.NOT. restart) THEN
341 : ! Initializing thermostat forces and velocities for the Nose-Hoover
342 : ! Chain variables
343 0 : IF (nhc%nhc_len /= 0) THEN
344 0 : CALL init_nhc_variables(nhc, simpar%temp_fast, para_env, globenv)
345 : END IF
346 : END IF
347 :
348 0 : CALL init_nhc_forces(nhc)
349 :
350 0 : CALL timestop(handle)
351 :
352 0 : END SUBROUTINE initialize_nhc_fast
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param thermostat_info ...
357 : !> \param simpar ...
358 : !> \param local_molecules ...
359 : !> \param molecule ...
360 : !> \param molecule_kind_set ...
361 : !> \param para_env ...
362 : !> \param globenv ...
363 : !> \param nhc ...
364 : !> \param nose_section ...
365 : !> \param gci ...
366 : !> \param save_mem ...
367 : !> \param binary_restart_file_name ...
368 : !> \author CJM
369 : ! **************************************************************************************************
370 752 : SUBROUTINE initialize_nhc_part(thermostat_info, simpar, local_molecules, &
371 : molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, &
372 : gci, save_mem, binary_restart_file_name)
373 :
374 : TYPE(thermostat_info_type), POINTER :: thermostat_info
375 : TYPE(simpar_type), POINTER :: simpar
376 : TYPE(distribution_1d_type), POINTER :: local_molecules
377 : TYPE(molecule_type), POINTER :: molecule(:)
378 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
379 : TYPE(mp_para_env_type), POINTER :: para_env
380 : TYPE(global_environment_type), POINTER :: globenv
381 : TYPE(lnhc_parameters_type), POINTER :: nhc
382 : TYPE(section_vals_type), POINTER :: nose_section
383 : TYPE(global_constraint_type), POINTER :: gci
384 : LOGICAL, INTENT(IN) :: save_mem
385 : CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name
386 :
387 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_part'
388 :
389 : INTEGER :: handle
390 : LOGICAL :: restart
391 :
392 376 : CALL timeset(routineN, handle)
393 :
394 : restart = .FALSE.
395 : ! fire up the thermostats, if not NVE
396 :
397 : CALL nhc_to_particle_mapping(thermostat_info, simpar, local_molecules, &
398 376 : molecule, molecule_kind_set, nhc, para_env, gci)
399 :
400 : ! Set up the Yoshida weights
401 376 : IF (nhc%nyosh > 0) THEN
402 1128 : ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
403 376 : CALL set_yoshida_coef(nhc, simpar%dt)
404 : END IF
405 :
406 : CALL restart_nose(nhc, nose_section, save_mem, restart, binary_restart_file_name, &
407 376 : "PARTICLE", para_env)
408 :
409 376 : IF (.NOT. restart) THEN
410 : ! Initializing thermostat forces and velocities for the Nose-Hoover
411 : ! Chain variables
412 300 : IF (nhc%nhc_len /= 0) THEN
413 300 : CALL init_nhc_variables(nhc, simpar%temp_ext, para_env, globenv)
414 : END IF
415 : END IF
416 :
417 376 : CALL init_nhc_forces(nhc)
418 :
419 376 : CALL timestop(handle)
420 :
421 376 : END SUBROUTINE initialize_nhc_part
422 :
423 : ! **************************************************************************************************
424 : !> \brief ...
425 : !> \param thermostat_info ...
426 : !> \param simpar ...
427 : !> \param local_molecules ...
428 : !> \param molecule ...
429 : !> \param molecule_kind_set ...
430 : !> \param para_env ...
431 : !> \param globenv ...
432 : !> \param nhc ...
433 : !> \param nose_section ...
434 : !> \param gci ...
435 : !> \param save_mem ...
436 : !> \param binary_restart_file_name ...
437 : !> \author MI
438 : ! **************************************************************************************************
439 80 : SUBROUTINE initialize_nhc_shell(thermostat_info, simpar, local_molecules, &
440 : molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, &
441 : gci, save_mem, binary_restart_file_name)
442 :
443 : TYPE(thermostat_info_type), POINTER :: thermostat_info
444 : TYPE(simpar_type), POINTER :: simpar
445 : TYPE(distribution_1d_type), POINTER :: local_molecules
446 : TYPE(molecule_type), POINTER :: molecule(:)
447 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
448 : TYPE(mp_para_env_type), POINTER :: para_env
449 : TYPE(global_environment_type), POINTER :: globenv
450 : TYPE(lnhc_parameters_type), POINTER :: nhc
451 : TYPE(section_vals_type), POINTER :: nose_section
452 : TYPE(global_constraint_type), POINTER :: gci
453 : LOGICAL, INTENT(IN) :: save_mem
454 : CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name
455 :
456 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_nhc_shell'
457 :
458 : INTEGER :: handle
459 : LOGICAL :: restart
460 :
461 40 : CALL timeset(routineN, handle)
462 :
463 : CALL nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, &
464 40 : molecule, molecule_kind_set, nhc, para_env, gci)
465 :
466 : restart = .FALSE.
467 : ! Set up the Yoshida weights
468 40 : IF (nhc%nyosh > 0) THEN
469 120 : ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
470 40 : CALL set_yoshida_coef(nhc, simpar%dt)
471 : END IF
472 :
473 : CALL restart_nose(nhc, nose_section, save_mem, restart, binary_restart_file_name, &
474 40 : "SHELL", para_env)
475 :
476 40 : IF (.NOT. restart) THEN
477 : ! Initialize thermostat forces and velocities
478 : ! Chain variables
479 28 : IF (nhc%nhc_len /= 0) THEN
480 28 : CALL init_nhc_variables(nhc, simpar%temp_sh_ext, para_env, globenv)
481 : END IF
482 : END IF
483 :
484 40 : CALL init_nhc_forces(nhc)
485 :
486 40 : CALL timestop(handle)
487 :
488 40 : END SUBROUTINE initialize_nhc_shell
489 :
490 : ! **************************************************************************************************
491 : !> \brief This lists the coefficients for the Yoshida method (higher
492 : !> order integrator used in NVT)
493 : !> \param nhc ...
494 : !> \param dt ...
495 : !> \date 14-NOV-2000
496 : !> \par History
497 : !> none
498 : ! **************************************************************************************************
499 536 : SUBROUTINE set_yoshida_coef(nhc, dt)
500 :
501 : TYPE(lnhc_parameters_type), POINTER :: nhc
502 : REAL(KIND=dp), INTENT(IN) :: dt
503 :
504 1072 : REAL(KIND=dp), DIMENSION(nhc%nyosh) :: yosh_wt
505 :
506 0 : SELECT CASE (nhc%nyosh)
507 : CASE DEFAULT
508 0 : CPABORT('Value not available.')
509 : CASE (1)
510 0 : yosh_wt(1) = 1.0_dp
511 : CASE (3)
512 536 : yosh_wt(1) = 1.0_dp/(2.0_dp - (2.0_dp)**(1.0_dp/3.0_dp))
513 536 : yosh_wt(2) = 1.0_dp - 2.0_dp*yosh_wt(1)
514 536 : yosh_wt(3) = yosh_wt(1)
515 : CASE (5)
516 0 : yosh_wt(1) = 1.0_dp/(4.0_dp - (4.0_dp)**(1.0_dp/3.0_dp))
517 0 : yosh_wt(2) = yosh_wt(1)
518 0 : yosh_wt(4) = yosh_wt(1)
519 0 : yosh_wt(5) = yosh_wt(1)
520 0 : yosh_wt(3) = 1.0_dp - 4.0_dp*yosh_wt(1)
521 : CASE (7)
522 0 : yosh_wt(1) = .78451361047756_dp
523 0 : yosh_wt(2) = .235573213359357_dp
524 0 : yosh_wt(3) = -1.17767998417887_dp
525 0 : yosh_wt(4) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + yosh_wt(3))
526 0 : yosh_wt(5) = yosh_wt(3)
527 0 : yosh_wt(6) = yosh_wt(2)
528 0 : yosh_wt(7) = yosh_wt(1)
529 : CASE (9)
530 0 : yosh_wt(1) = 0.192_dp
531 0 : yosh_wt(2) = 0.554910818409783619692725006662999_dp
532 0 : yosh_wt(3) = 0.124659619941888644216504240951585_dp
533 0 : yosh_wt(4) = -0.843182063596933505315033808282941_dp
534 : yosh_wt(5) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
535 0 : yosh_wt(3) + yosh_wt(4))
536 0 : yosh_wt(6) = yosh_wt(4)
537 0 : yosh_wt(7) = yosh_wt(3)
538 0 : yosh_wt(8) = yosh_wt(2)
539 0 : yosh_wt(9) = yosh_wt(1)
540 : CASE (15)
541 0 : yosh_wt(1) = 0.102799849391985_dp
542 0 : yosh_wt(2) = -0.196061023297549e1_dp
543 0 : yosh_wt(3) = 0.193813913762276e1_dp
544 0 : yosh_wt(4) = -0.158240635368243_dp
545 0 : yosh_wt(5) = -0.144485223686048e1_dp
546 0 : yosh_wt(6) = 0.253693336566229_dp
547 0 : yosh_wt(7) = 0.914844246229740_dp
548 : yosh_wt(8) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
549 0 : yosh_wt(3) + yosh_wt(4) + yosh_wt(5) + yosh_wt(6) + yosh_wt(7))
550 0 : yosh_wt(9) = yosh_wt(7)
551 0 : yosh_wt(10) = yosh_wt(6)
552 0 : yosh_wt(11) = yosh_wt(5)
553 0 : yosh_wt(12) = yosh_wt(4)
554 0 : yosh_wt(13) = yosh_wt(3)
555 0 : yosh_wt(14) = yosh_wt(2)
556 536 : yosh_wt(15) = yosh_wt(1)
557 : END SELECT
558 2144 : nhc%dt_yosh = dt*yosh_wt/REAL(nhc%nc, KIND=dp)
559 :
560 536 : END SUBROUTINE set_yoshida_coef
561 :
562 : ! **************************************************************************************************
563 : !> \brief read coordinate, velocities, forces and masses of the
564 : !> thermostat from restart file
565 : !> \param nhc ...
566 : !> \param nose_section ...
567 : !> \param save_mem ...
568 : !> \param restart ...
569 : !> \param binary_restart_file_name ...
570 : !> \param thermostat_name ...
571 : !> \param para_env ...
572 : !> \par History
573 : !> 24-07-07 created
574 : !> \author MI
575 : ! **************************************************************************************************
576 536 : SUBROUTINE restart_nose(nhc, nose_section, save_mem, restart, &
577 : binary_restart_file_name, thermostat_name, &
578 : para_env)
579 :
580 : TYPE(lnhc_parameters_type), POINTER :: nhc
581 : TYPE(section_vals_type), POINTER :: nose_section
582 : LOGICAL, INTENT(IN) :: save_mem
583 : LOGICAL, INTENT(OUT) :: restart
584 : CHARACTER(LEN=*), INTENT(IN) :: binary_restart_file_name, thermostat_name
585 : TYPE(mp_para_env_type), POINTER :: para_env
586 :
587 : CHARACTER(len=*), PARAMETER :: routineN = 'restart_nose'
588 :
589 : INTEGER :: handle, i, ind, j
590 : LOGICAL :: explicit
591 536 : REAL(KIND=dp), DIMENSION(:), POINTER :: buffer
592 : TYPE(map_info_type), POINTER :: map_info
593 : TYPE(section_vals_type), POINTER :: work_section
594 :
595 536 : CALL timeset(routineN, handle)
596 :
597 536 : NULLIFY (buffer)
598 536 : NULLIFY (work_section)
599 :
600 536 : IF (LEN_TRIM(binary_restart_file_name) > 0) THEN
601 :
602 : ! Read binary restart file, if available
603 :
604 : CALL read_binary_thermostats_nose(thermostat_name, nhc, binary_restart_file_name, &
605 38 : restart, para_env)
606 :
607 : ELSE
608 :
609 : ! Read the default restart file in ASCII format
610 :
611 : explicit = .FALSE.
612 498 : restart = .FALSE.
613 :
614 498 : IF (ASSOCIATED(nose_section)) THEN
615 498 : work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
616 498 : CALL section_vals_get(work_section, explicit=explicit)
617 498 : restart = explicit
618 498 : work_section => section_vals_get_subs_vals(nose_section, "COORD")
619 498 : CALL section_vals_get(work_section, explicit=explicit)
620 498 : IF (.NOT. restart .AND. explicit) THEN
621 : CALL cp_abort(__LOCATION__, "You need to define both VELOCITY and "// &
622 0 : "COORD and MASS and FORCE section (or none) in the NOSE section")
623 : END IF
624 498 : restart = explicit .AND. restart
625 498 : work_section => section_vals_get_subs_vals(nose_section, "MASS")
626 498 : CALL section_vals_get(work_section, explicit=explicit)
627 498 : IF (.NOT. restart .AND. explicit) THEN
628 : CALL cp_abort(__LOCATION__, "You need to define both VELOCITY and "// &
629 0 : "COORD and MASS and FORCE section (or none) in the NOSE section")
630 : END IF
631 498 : restart = explicit .AND. restart
632 498 : work_section => section_vals_get_subs_vals(nose_section, "FORCE")
633 498 : CALL section_vals_get(work_section, explicit=explicit)
634 498 : IF (.NOT. restart .AND. explicit) THEN
635 : CALL cp_abort(__LOCATION__, "You need to define both VELOCITY and "// &
636 0 : "COORD and MASS and FORCE section (or none) in the NOSE section")
637 : END IF
638 922 : restart = explicit .AND. restart
639 : END IF
640 :
641 498 : IF (restart) THEN
642 74 : map_info => nhc%map_info
643 74 : CALL section_vals_val_get(nose_section, "COORD%_DEFAULT_KEYWORD_", r_vals=buffer)
644 4442 : DO i = 1, SIZE(nhc%nvt, 2)
645 4368 : ind = map_info%index(i)
646 4368 : ind = (ind - 1)*nhc%nhc_len
647 17928 : DO j = 1, SIZE(nhc%nvt, 1)
648 13486 : ind = ind + 1
649 17854 : nhc%nvt(j, i)%eta = buffer(ind)
650 : END DO
651 : END DO
652 74 : CALL section_vals_val_get(nose_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
653 4442 : DO i = 1, SIZE(nhc%nvt, 2)
654 4368 : ind = map_info%index(i)
655 4368 : ind = (ind - 1)*nhc%nhc_len
656 17928 : DO j = 1, SIZE(nhc%nvt, 1)
657 13486 : ind = ind + 1
658 17854 : nhc%nvt(j, i)%v = buffer(ind)
659 : END DO
660 : END DO
661 74 : CALL section_vals_val_get(nose_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
662 4442 : DO i = 1, SIZE(nhc%nvt, 2)
663 4368 : ind = map_info%index(i)
664 4368 : ind = (ind - 1)*nhc%nhc_len
665 17928 : DO j = 1, SIZE(nhc%nvt, 1)
666 13486 : ind = ind + 1
667 17854 : nhc%nvt(j, i)%mass = buffer(ind)
668 : END DO
669 : END DO
670 74 : CALL section_vals_val_get(nose_section, "FORCE%_DEFAULT_KEYWORD_", r_vals=buffer)
671 4442 : DO i = 1, SIZE(nhc%nvt, 2)
672 4368 : ind = map_info%index(i)
673 4368 : ind = (ind - 1)*nhc%nhc_len
674 17928 : DO j = 1, SIZE(nhc%nvt, 1)
675 13486 : ind = ind + 1
676 17854 : nhc%nvt(j, i)%f = buffer(ind)
677 : END DO
678 : END DO
679 : END IF
680 :
681 498 : IF (save_mem) THEN
682 2 : NULLIFY (work_section)
683 2 : work_section => section_vals_get_subs_vals(nose_section, "COORD")
684 2 : CALL section_vals_remove_values(work_section)
685 2 : NULLIFY (work_section)
686 2 : work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
687 2 : CALL section_vals_remove_values(work_section)
688 2 : NULLIFY (work_section)
689 2 : work_section => section_vals_get_subs_vals(nose_section, "FORCE")
690 2 : CALL section_vals_remove_values(work_section)
691 2 : NULLIFY (work_section)
692 2 : work_section => section_vals_get_subs_vals(nose_section, "MASS")
693 2 : CALL section_vals_remove_values(work_section)
694 : END IF
695 :
696 : END IF
697 :
698 536 : CALL timestop(handle)
699 :
700 536 : END SUBROUTINE restart_nose
701 :
702 : ! **************************************************************************************************
703 : !> \brief Initializes the NHC velocities to the Maxwellian distribution
704 : !> \param nhc ...
705 : !> \param temp_ext ...
706 : !> \param para_env ...
707 : !> \param globenv ...
708 : !> \date 14-NOV-2000
709 : !> \par History
710 : !> none
711 : ! **************************************************************************************************
712 430 : SUBROUTINE init_nhc_variables(nhc, temp_ext, para_env, globenv)
713 : TYPE(lnhc_parameters_type), POINTER :: nhc
714 : REAL(KIND=dp), INTENT(IN) :: temp_ext
715 : TYPE(mp_para_env_type), POINTER :: para_env
716 : TYPE(global_environment_type), POINTER :: globenv
717 :
718 : CHARACTER(len=*), PARAMETER :: routineN = 'init_nhc_variables'
719 :
720 : INTEGER :: handle, i, icount, j, number, tot_rn
721 : REAL(KIND=dp) :: akin, dum, temp
722 430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: array_of_rn
723 : TYPE(map_info_type), POINTER :: map_info
724 :
725 430 : CALL timeset(routineN, handle)
726 :
727 430 : tot_rn = 0
728 :
729 : ! first initializing the mass of the nhc variables
730 72430 : nhc%nvt(:, :)%mass = nhc%nvt(:, :)%nkt*nhc%tau_nhc**2
731 72430 : nhc%nvt(:, :)%eta = 0._dp
732 72430 : nhc%nvt(:, :)%v = 0._dp
733 72430 : nhc%nvt(:, :)%f = 0._dp
734 :
735 430 : map_info => nhc%map_info
736 758 : SELECT CASE (map_info%dis_type)
737 : CASE (do_thermo_only_master) ! for NPT
738 : CASE DEFAULT
739 328 : tot_rn = nhc%glob_num_nhc*nhc%nhc_len
740 :
741 984 : ALLOCATE (array_of_rn(tot_rn))
742 758 : array_of_rn(:) = 0.0_dp
743 : END SELECT
744 :
745 102 : SELECT CASE (map_info%dis_type)
746 : CASE (do_thermo_only_master) ! for NPT
747 : ! Map deterministically determined random number to nhc % v
748 204 : DO i = 1, nhc%loc_num_nhc
749 522 : DO j = 1, nhc%nhc_len
750 420 : nhc%nvt(j, i)%v = globenv%gaussian_rng_stream%next()
751 : END DO
752 : END DO
753 :
754 102 : akin = 0.0_dp
755 204 : DO i = 1, nhc%loc_num_nhc
756 522 : DO j = 1, nhc%nhc_len
757 : akin = akin + 0.5_dp*(nhc%nvt(j, i)%mass* &
758 : nhc%nvt(j, i)%v* &
759 420 : nhc%nvt(j, i)%v)
760 : END DO
761 : END DO
762 102 : number = nhc%loc_num_nhc
763 :
764 : ! scale velocities to get the correct initial temperature
765 102 : temp = 2.0_dp*akin/REAL(number, KIND=dp)
766 102 : IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp)
767 204 : DO i = 1, nhc%loc_num_nhc
768 522 : DO j = 1, nhc%nhc_len
769 318 : nhc%nvt(j, i)%v = temp*nhc%nvt(j, i)%v
770 420 : nhc%nvt(j, i)%eta = 0.0_dp
771 : END DO
772 : END DO
773 :
774 : ! initializing all of the forces on the thermostats
775 204 : DO i = 1, nhc%loc_num_nhc
776 420 : DO j = 2, nhc%nhc_len
777 : nhc%nvt(j, i)%f = nhc%nvt(j - 1, i)%mass*nhc%nvt(j - 1, i)%v* &
778 216 : nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
779 318 : IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
780 216 : nhc%nvt(j, i)%f = nhc%nvt(j, i)%f/nhc%nvt(j, i)%mass
781 : END IF
782 : END DO
783 : END DO
784 :
785 : CASE DEFAULT
786 110524 : DO i = 1, tot_rn
787 110524 : array_of_rn(i) = globenv%gaussian_rng_stream%next()
788 : END DO
789 : ! Map deterministically determined random number to nhc % v
790 16503 : DO i = 1, nhc%loc_num_nhc
791 16175 : icount = map_info%index(i)
792 16175 : icount = (icount - 1)*nhc%nhc_len
793 71908 : DO j = 1, nhc%nhc_len
794 55405 : icount = icount + 1
795 55405 : nhc%nvt(j, i)%v = array_of_rn(icount)
796 : ! WRITE ( *, * ) 'VEL', para_env%mepos, i,j, nhc%nvt(j,i)%v
797 71580 : nhc%nvt(j, i)%eta = 0.0_dp
798 : END DO
799 : END DO
800 328 : DEALLOCATE (array_of_rn)
801 :
802 328 : number = nhc%glob_num_nhc
803 328 : CALL get_nhc_energies(nhc, dum, akin, para_env)
804 :
805 : ! scale velocities to get the correct initial temperature
806 328 : temp = 2.0_dp*akin/REAL(number, KIND=dp)
807 328 : IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp)
808 16503 : DO i = 1, nhc%loc_num_nhc
809 71908 : DO j = 1, nhc%nhc_len
810 71580 : nhc%nvt(j, i)%v = temp*nhc%nvt(j, i)%v
811 : END DO
812 : END DO
813 :
814 : ! initializing all of the forces on the thermostats
815 17261 : DO i = 1, nhc%loc_num_nhc
816 55733 : DO j = 2, nhc%nhc_len
817 : nhc%nvt(j, i)%f = nhc%nvt(j - 1, i)%mass*nhc%nvt(j - 1, i)%v* &
818 39230 : nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
819 55405 : IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
820 38654 : nhc%nvt(j, i)%f = nhc%nvt(j, i)%f/nhc%nvt(j, i)%mass
821 : END IF
822 : END DO
823 : END DO
824 :
825 : END SELECT
826 :
827 430 : CALL timestop(handle)
828 :
829 430 : END SUBROUTINE init_nhc_variables
830 :
831 : ! **************************************************************************************************
832 : !> \brief Initializes the barostat velocities to the Maxwellian distribution
833 : !> \param npt ...
834 : !> \param tau_cell ...
835 : !> \param temp_ext ...
836 : !> \param nfree ...
837 : !> \param ensemble ...
838 : !> \param cmass ...
839 : !> \param globenv ...
840 : !> \date 14-NOV-2000
841 : !> \par History
842 : !> none
843 : ! **************************************************************************************************
844 152 : SUBROUTINE init_barostat_variables(npt, tau_cell, temp_ext, nfree, ensemble, &
845 : cmass, globenv)
846 :
847 : TYPE(npt_info_type), DIMENSION(:, :), &
848 : INTENT(INOUT) :: npt
849 : REAL(KIND=dp), INTENT(IN) :: tau_cell, temp_ext
850 : INTEGER, INTENT(IN) :: nfree, ensemble
851 : REAL(KIND=dp), INTENT(IN) :: cmass
852 : TYPE(global_environment_type), POINTER :: globenv
853 :
854 : CHARACTER(len=*), PARAMETER :: routineN = 'init_barostat_variables'
855 :
856 : INTEGER :: handle, i, j, number
857 : REAL(KIND=dp) :: akin, temp, v
858 :
859 152 : CALL timeset(routineN, handle)
860 :
861 152 : temp = 0.0_dp
862 :
863 : ! first initializing the mass of the nhc variables
864 :
865 916 : npt(:, :)%eps = 0.0_dp
866 916 : npt(:, :)%v = 0.0_dp
867 916 : npt(:, :)%f = 0.0_dp
868 250 : SELECT CASE (ensemble)
869 : CASE (npt_i_ensemble, npt_ia_ensemble)
870 294 : npt(:, :)%mass = REAL(nfree + 3, KIND=dp)*temp_ext*tau_cell**2
871 : CASE (npt_f_ensemble)
872 468 : npt(:, :)%mass = REAL(nfree + 3, KIND=dp)*temp_ext*tau_cell**2/3.0_dp
873 : CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
874 18 : npt(:, :)%mass = cmass
875 : CASE (npe_f_ensemble)
876 130 : npt(:, :)%mass = REAL(nfree + 3, KIND=dp)*temp_ext*tau_cell**2/3.0_dp
877 : CASE (npe_i_ensemble)
878 156 : npt(:, :)%mass = REAL(nfree + 3, KIND=dp)*temp_ext*tau_cell**2
879 : END SELECT
880 : ! initializing velocities
881 396 : DO i = 1, SIZE(npt, 1)
882 778 : DO j = i, SIZE(npt, 2)
883 382 : v = globenv%gaussian_rng_stream%next()
884 : ! Symmetrizing the initial barostat velocities to ensure
885 : ! no rotation of the cell under NPT_F
886 382 : npt(j, i)%v = v
887 626 : npt(i, j)%v = v
888 : END DO
889 : END DO
890 :
891 396 : akin = 0.0_dp
892 396 : DO i = 1, SIZE(npt, 1)
893 916 : DO j = 1, SIZE(npt, 2)
894 764 : akin = akin + 0.5_dp*(npt(j, i)%mass*npt(j, i)%v*npt(j, i)%v)
895 : END DO
896 : END DO
897 :
898 152 : number = SIZE(npt, 1)*SIZE(npt, 2)
899 :
900 : ! scale velocities to get the correct initial temperature
901 152 : IF (number /= 0) THEN
902 152 : temp = 2.0_dp*akin/REAL(number, KIND=dp)
903 152 : IF (temp > 0.0_dp) temp = SQRT(temp_ext/temp)
904 : END IF
905 396 : DO i = 1, SIZE(npt, 1)
906 778 : DO j = i, SIZE(npt, 2)
907 382 : npt(j, i)%v = temp*npt(j, i)%v
908 382 : npt(i, j)%v = npt(j, i)%v
909 244 : IF (debug_isotropic_limit) THEN
910 : npt(j, i)%v = 0.0_dp
911 : npt(i, j)%v = 0.0_dp
912 : WRITE (*, *) 'DEBUG ISOTROPIC LIMIT| INITIAL v_eps', npt(j, i)%v
913 : END IF
914 : END DO
915 : END DO
916 :
917 : ! Zero barostat velocities for nph_uniaxial
918 : SELECT CASE (ensemble)
919 : ! Zero barostat velocities for nph_uniaxial
920 : CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
921 164 : npt(:, :)%v = 0.0_dp
922 : END SELECT
923 :
924 152 : CALL timestop(handle)
925 :
926 152 : END SUBROUTINE init_barostat_variables
927 :
928 : ! **************************************************************************************************
929 : !> \brief Assigns extended parameters from the restart file.
930 : !> \param nhc ...
931 : !> \author CJM
932 : ! **************************************************************************************************
933 536 : SUBROUTINE init_nhc_forces(nhc)
934 :
935 : TYPE(lnhc_parameters_type), POINTER :: nhc
936 :
937 : CHARACTER(len=*), PARAMETER :: routineN = 'init_nhc_forces'
938 :
939 : INTEGER :: handle, i, j
940 :
941 536 : CALL timeset(routineN, handle)
942 :
943 536 : CPASSERT(ASSOCIATED(nhc))
944 : ! assign the forces
945 25789 : DO i = 1, SIZE(nhc%nvt, 2)
946 83569 : DO j = 2, SIZE(nhc%nvt, 1)
947 : nhc%nvt(j, i)%f = nhc%nvt(j - 1, i)%mass* &
948 : nhc%nvt(j - 1, i)%v**2 - &
949 57780 : nhc%nvt(j, i)%nkt
950 83033 : IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
951 57204 : nhc%nvt(j, i)%f = nhc%nvt(j, i)%f/nhc%nvt(j, i)%mass
952 : END IF
953 : END DO
954 : END DO
955 :
956 536 : CALL timestop(handle)
957 :
958 536 : END SUBROUTINE init_nhc_forces
959 :
960 : END MODULE extended_system_init
|
