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