(git:71c3ab0)
Loading...
Searching...
No Matches
extended_system_init.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
62CONTAINS
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 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 REAL(kind=dp), DIMENSION(:), POINTER :: buffer
87 TYPE(section_vals_type), POINTER :: work_section2
88
89 CALL timeset(routinen, handle)
90
91 NULLIFY (work_section2)
92
93 explicit = .false.
94 restart = .false.
95
96 cpassert(.NOT. ASSOCIATED(npt_info))
97
98 ! first allocating the npt_info_type if requested
99 SELECT CASE (simpar%ensemble)
100 CASE (npt_i_ensemble, npe_i_ensemble, npt_ia_ensemble)
101 ALLOCATE (npt_info(1, 1))
102 npt_info(:, :)%eps = log(cell%deth)/3.0_dp
103 temp = simpar%temp_baro_ext
104
105 CASE (npt_f_ensemble, npe_f_ensemble)
106 ALLOCATE (npt_info(3, 3))
107 temp = simpar%temp_baro_ext
108
109 CASE (nph_uniaxial_ensemble)
110 ALLOCATE (npt_info(1, 1))
111 temp = simpar%temp_baro_ext
112
113 CASE (nph_uniaxial_damped_ensemble)
114 ALLOCATE (npt_info(1, 1))
115 temp = simpar%temp_baro_ext
116
117 CASE DEFAULT
118 ! Do nothing..
119 NULLIFY (npt_info)
120 END SELECT
121
122 IF (ASSOCIATED(npt_info)) THEN
123 IF (ASSOCIATED(work_section)) THEN
124 work_section2 => section_vals_get_subs_vals(work_section, "VELOCITY")
125 CALL section_vals_get(work_section2, explicit=explicit)
126 restart = explicit
127 work_section2 => section_vals_get_subs_vals(work_section, "MASS")
128 CALL section_vals_get(work_section2, explicit=explicit)
129 IF (restart .NEQV. explicit) THEN
130 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
131 "MASS section (or none) in the BAROSTAT section")
132 END IF
133 restart = explicit .AND. restart
134 END IF
135
136 IF (restart) THEN
137 CALL section_vals_val_get(work_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
138 ind = 0
139 DO i = 1, SIZE(npt_info, 1)
140 DO j = 1, SIZE(npt_info, 2)
141 ind = ind + 1
142 npt_info(i, j)%v = buffer(ind)
143 END DO
144 END DO
145 CALL section_vals_val_get(work_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
146 ind = 0
147 DO i = 1, SIZE(npt_info, 1)
148 DO j = 1, SIZE(npt_info, 2)
149 ind = ind + 1
150 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 globenv)
157 END IF
158
159 END IF
160
161 CALL timestop(handle)
162
163 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 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 CALL timeset(routinen, handle)
191
192 restart = .false.
193
194 CALL nhc_to_barostat_mapping(simpar, nhc)
195
196 ! Set up the Yoshida weights
197 IF (nhc%nyosh > 0) THEN
198 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
199 CALL set_yoshida_coef(nhc, simpar%dt)
200 END IF
201
202 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
203
204 IF (.NOT. restart) THEN
205 ! Initializing thermostat forces and velocities for the Nose-Hoover
206 ! Chain variables
207 SELECT CASE (simpar%ensemble)
208 CASE DEFAULT
209 temp = simpar%temp_baro_ext
210 END SELECT
211 IF (nhc%nhc_len /= 0) THEN
212 CALL init_nhc_variables(nhc, temp, para_env, globenv)
213 END IF
214 END IF
215
216 CALL init_nhc_forces(nhc)
217
218 CALL timestop(handle)
219
220 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 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 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 molecule, molecule_kind_set, nhc, para_env, gci)
265
266 ! Set up the Yoshida weights
267 IF (nhc%nyosh > 0) THEN
268 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
269 CALL set_yoshida_coef(nhc, simpar%dt)
270 END IF
271
272 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
273
274 IF (.NOT. restart) THEN
275 ! Initializing thermostat forces and velocities for the Nose-Hoover
276 ! Chain variables
277 IF (nhc%nhc_len /= 0) THEN
278 CALL init_nhc_variables(nhc, simpar%temp_slow, para_env, globenv)
279 END IF
280 END IF
281
282 CALL init_nhc_forces(nhc)
283
284 CALL timestop(handle)
285
286 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 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 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 molecule, molecule_kind_set, nhc, para_env, gci)
331
332 ! Set up the Yoshida weights
333 IF (nhc%nyosh > 0) THEN
334 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
335 CALL set_yoshida_coef(nhc, simpar%dt)
336 END IF
337
338 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
339
340 IF (.NOT. restart) THEN
341 ! Initializing thermostat forces and velocities for the Nose-Hoover
342 ! Chain variables
343 IF (nhc%nhc_len /= 0) THEN
344 CALL init_nhc_variables(nhc, simpar%temp_fast, para_env, globenv)
345 END IF
346 END IF
347
348 CALL init_nhc_forces(nhc)
349
350 CALL timestop(handle)
351
352 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 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 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 molecule, molecule_kind_set, nhc, para_env, gci)
399
400 ! Set up the Yoshida weights
401 IF (nhc%nyosh > 0) THEN
402 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
403 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 "PARTICLE", para_env)
408
409 IF (.NOT. restart) THEN
410 ! Initializing thermostat forces and velocities for the Nose-Hoover
411 ! Chain variables
412 IF (nhc%nhc_len /= 0) THEN
413 CALL init_nhc_variables(nhc, simpar%temp_ext, para_env, globenv)
414 END IF
415 END IF
416
417 CALL init_nhc_forces(nhc)
418
419 CALL timestop(handle)
420
421 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 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 CALL timeset(routinen, handle)
462
463 CALL nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, &
464 molecule, molecule_kind_set, nhc, para_env, gci)
465
466 restart = .false.
467 ! Set up the Yoshida weights
468 IF (nhc%nyosh > 0) THEN
469 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
470 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 "SHELL", para_env)
475
476 IF (.NOT. restart) THEN
477 ! Initialize thermostat forces and velocities
478 ! Chain variables
479 IF (nhc%nhc_len /= 0) THEN
480 CALL init_nhc_variables(nhc, simpar%temp_sh_ext, para_env, globenv)
481 END IF
482 END IF
483
484 CALL init_nhc_forces(nhc)
485
486 CALL timestop(handle)
487
488 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 SUBROUTINE set_yoshida_coef(nhc, dt)
500
501 TYPE(lnhc_parameters_type), POINTER :: nhc
502 REAL(kind=dp), INTENT(IN) :: dt
503
504 REAL(kind=dp), DIMENSION(nhc%nyosh) :: yosh_wt
505
506 SELECT CASE (nhc%nyosh)
507 CASE DEFAULT
508 cpabort('Value not available.')
509 CASE (1)
510 yosh_wt(1) = 1.0_dp
511 CASE (3)
512 yosh_wt(1) = 1.0_dp/(2.0_dp - (2.0_dp)**(1.0_dp/3.0_dp))
513 yosh_wt(2) = 1.0_dp - 2.0_dp*yosh_wt(1)
514 yosh_wt(3) = yosh_wt(1)
515 CASE (5)
516 yosh_wt(1) = 1.0_dp/(4.0_dp - (4.0_dp)**(1.0_dp/3.0_dp))
517 yosh_wt(2) = yosh_wt(1)
518 yosh_wt(4) = yosh_wt(1)
519 yosh_wt(5) = yosh_wt(1)
520 yosh_wt(3) = 1.0_dp - 4.0_dp*yosh_wt(1)
521 CASE (7)
522 yosh_wt(1) = .78451361047756_dp
523 yosh_wt(2) = .235573213359357_dp
524 yosh_wt(3) = -1.17767998417887_dp
525 yosh_wt(4) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + yosh_wt(3))
526 yosh_wt(5) = yosh_wt(3)
527 yosh_wt(6) = yosh_wt(2)
528 yosh_wt(7) = yosh_wt(1)
529 CASE (9)
530 yosh_wt(1) = 0.192_dp
531 yosh_wt(2) = 0.554910818409783619692725006662999_dp
532 yosh_wt(3) = 0.124659619941888644216504240951585_dp
533 yosh_wt(4) = -0.843182063596933505315033808282941_dp
534 yosh_wt(5) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
535 yosh_wt(3) + yosh_wt(4))
536 yosh_wt(6) = yosh_wt(4)
537 yosh_wt(7) = yosh_wt(3)
538 yosh_wt(8) = yosh_wt(2)
539 yosh_wt(9) = yosh_wt(1)
540 CASE (15)
541 yosh_wt(1) = 0.102799849391985_dp
542 yosh_wt(2) = -0.196061023297549e1_dp
543 yosh_wt(3) = 0.193813913762276e1_dp
544 yosh_wt(4) = -0.158240635368243_dp
545 yosh_wt(5) = -0.144485223686048e1_dp
546 yosh_wt(6) = 0.253693336566229_dp
547 yosh_wt(7) = 0.914844246229740_dp
548 yosh_wt(8) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
549 yosh_wt(3) + yosh_wt(4) + yosh_wt(5) + yosh_wt(6) + yosh_wt(7))
550 yosh_wt(9) = yosh_wt(7)
551 yosh_wt(10) = yosh_wt(6)
552 yosh_wt(11) = yosh_wt(5)
553 yosh_wt(12) = yosh_wt(4)
554 yosh_wt(13) = yosh_wt(3)
555 yosh_wt(14) = yosh_wt(2)
556 yosh_wt(15) = yosh_wt(1)
557 END SELECT
558 nhc%dt_yosh = dt*yosh_wt/real(nhc%nc, kind=dp)
559
560 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: buffer
592 TYPE(map_info_type), POINTER :: map_info
593 TYPE(section_vals_type), POINTER :: work_section
594
595 CALL timeset(routinen, handle)
596
597 NULLIFY (buffer)
598 NULLIFY (work_section)
599
600 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 restart, para_env)
606
607 ELSE
608
609 ! Read the default restart file in ASCII format
610
611 explicit = .false.
612 restart = .false.
613
614 IF (ASSOCIATED(nose_section)) THEN
615 work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
616 CALL section_vals_get(work_section, explicit=explicit)
617 restart = explicit
618 work_section => section_vals_get_subs_vals(nose_section, "COORD")
619 CALL section_vals_get(work_section, explicit=explicit)
620 IF (.NOT. restart .AND. explicit) THEN
621 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
622 "COORD and MASS and FORCE section (or none) in the NOSE section")
623 END IF
624 restart = explicit .AND. restart
625 work_section => section_vals_get_subs_vals(nose_section, "MASS")
626 CALL section_vals_get(work_section, explicit=explicit)
627 IF (.NOT. restart .AND. explicit) THEN
628 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
629 "COORD and MASS and FORCE section (or none) in the NOSE section")
630 END IF
631 restart = explicit .AND. restart
632 work_section => section_vals_get_subs_vals(nose_section, "FORCE")
633 CALL section_vals_get(work_section, explicit=explicit)
634 IF (.NOT. restart .AND. explicit) THEN
635 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
636 "COORD and MASS and FORCE section (or none) in the NOSE section")
637 END IF
638 restart = explicit .AND. restart
639 END IF
640
641 IF (restart) THEN
642 map_info => nhc%map_info
643 CALL section_vals_val_get(nose_section, "COORD%_DEFAULT_KEYWORD_", r_vals=buffer)
644 DO i = 1, SIZE(nhc%nvt, 2)
645 ind = map_info%index(i)
646 ind = (ind - 1)*nhc%nhc_len
647 DO j = 1, SIZE(nhc%nvt, 1)
648 ind = ind + 1
649 nhc%nvt(j, i)%eta = buffer(ind)
650 END DO
651 END DO
652 CALL section_vals_val_get(nose_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
653 DO i = 1, SIZE(nhc%nvt, 2)
654 ind = map_info%index(i)
655 ind = (ind - 1)*nhc%nhc_len
656 DO j = 1, SIZE(nhc%nvt, 1)
657 ind = ind + 1
658 nhc%nvt(j, i)%v = buffer(ind)
659 END DO
660 END DO
661 CALL section_vals_val_get(nose_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
662 DO i = 1, SIZE(nhc%nvt, 2)
663 ind = map_info%index(i)
664 ind = (ind - 1)*nhc%nhc_len
665 DO j = 1, SIZE(nhc%nvt, 1)
666 ind = ind + 1
667 nhc%nvt(j, i)%mass = buffer(ind)
668 END DO
669 END DO
670 CALL section_vals_val_get(nose_section, "FORCE%_DEFAULT_KEYWORD_", r_vals=buffer)
671 DO i = 1, SIZE(nhc%nvt, 2)
672 ind = map_info%index(i)
673 ind = (ind - 1)*nhc%nhc_len
674 DO j = 1, SIZE(nhc%nvt, 1)
675 ind = ind + 1
676 nhc%nvt(j, i)%f = buffer(ind)
677 END DO
678 END DO
679 END IF
680
681 IF (save_mem) THEN
682 NULLIFY (work_section)
683 work_section => section_vals_get_subs_vals(nose_section, "COORD")
684 CALL section_vals_remove_values(work_section)
685 NULLIFY (work_section)
686 work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
687 CALL section_vals_remove_values(work_section)
688 NULLIFY (work_section)
689 work_section => section_vals_get_subs_vals(nose_section, "FORCE")
690 CALL section_vals_remove_values(work_section)
691 NULLIFY (work_section)
692 work_section => section_vals_get_subs_vals(nose_section, "MASS")
693 CALL section_vals_remove_values(work_section)
694 END IF
695
696 END IF
697
698 CALL timestop(handle)
699
700 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: array_of_rn
723 TYPE(map_info_type), POINTER :: map_info
724
725 CALL timeset(routinen, handle)
726
727 tot_rn = 0
728
729 ! first initializing the mass of the nhc variables
730 nhc%nvt(:, :)%mass = nhc%nvt(:, :)%nkt*nhc%tau_nhc**2
731 nhc%nvt(:, :)%eta = 0._dp
732 nhc%nvt(:, :)%v = 0._dp
733 nhc%nvt(:, :)%f = 0._dp
734
735 map_info => nhc%map_info
736 SELECT CASE (map_info%dis_type)
737 CASE (do_thermo_only_master) ! for NPT
738 CASE DEFAULT
739 tot_rn = nhc%glob_num_nhc*nhc%nhc_len
740
741 ALLOCATE (array_of_rn(tot_rn))
742 array_of_rn(:) = 0.0_dp
743 END SELECT
744
745 SELECT CASE (map_info%dis_type)
746 CASE (do_thermo_only_master) ! for NPT
747 ! Map deterministically determined random number to nhc % v
748 DO i = 1, nhc%loc_num_nhc
749 DO j = 1, nhc%nhc_len
750 nhc%nvt(j, i)%v = globenv%gaussian_rng_stream%next()
751 END DO
752 END DO
753
754 akin = 0.0_dp
755 DO i = 1, nhc%loc_num_nhc
756 DO j = 1, nhc%nhc_len
757 akin = akin + 0.5_dp*(nhc%nvt(j, i)%mass* &
758 nhc%nvt(j, i)%v* &
759 nhc%nvt(j, i)%v)
760 END DO
761 END DO
762 number = nhc%loc_num_nhc
763
764 ! scale velocities to get the correct initial temperature
765 temp = 2.0_dp*akin/real(number, kind=dp)
766 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
767 DO i = 1, nhc%loc_num_nhc
768 DO j = 1, nhc%nhc_len
769 nhc%nvt(j, i)%v = temp*nhc%nvt(j, i)%v
770 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 DO i = 1, nhc%loc_num_nhc
776 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 nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
779 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
780 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 DO i = 1, tot_rn
787 array_of_rn(i) = globenv%gaussian_rng_stream%next()
788 END DO
789 ! Map deterministically determined random number to nhc % v
790 DO i = 1, nhc%loc_num_nhc
791 icount = map_info%index(i)
792 icount = (icount - 1)*nhc%nhc_len
793 DO j = 1, nhc%nhc_len
794 icount = icount + 1
795 nhc%nvt(j, i)%v = array_of_rn(icount)
796 ! WRITE ( *, * ) 'VEL', para_env%mepos, i,j, nhc%nvt(j,i)%v
797 nhc%nvt(j, i)%eta = 0.0_dp
798 END DO
799 END DO
800 DEALLOCATE (array_of_rn)
801
802 number = nhc%glob_num_nhc
803 CALL get_nhc_energies(nhc, dum, akin, para_env)
804
805 ! scale velocities to get the correct initial temperature
806 temp = 2.0_dp*akin/real(number, kind=dp)
807 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
808 DO i = 1, nhc%loc_num_nhc
809 DO j = 1, nhc%nhc_len
810 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 DO i = 1, nhc%loc_num_nhc
816 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 nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
819 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
820 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 CALL timestop(handle)
828
829 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 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 CALL timeset(routinen, handle)
860
861 temp = 0.0_dp
862
863 ! first initializing the mass of the nhc variables
864
865 npt(:, :)%eps = 0.0_dp
866 npt(:, :)%v = 0.0_dp
867 npt(:, :)%f = 0.0_dp
868 SELECT CASE (ensemble)
869 CASE (npt_i_ensemble, npt_ia_ensemble)
870 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2
871 CASE (npt_f_ensemble)
872 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 npt(:, :)%mass = cmass
875 CASE (npe_f_ensemble)
876 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2/3.0_dp
877 CASE (npe_i_ensemble)
878 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2
879 END SELECT
880 ! initializing velocities
881 DO i = 1, SIZE(npt, 1)
882 DO j = i, SIZE(npt, 2)
883 v = globenv%gaussian_rng_stream%next()
884 ! Symmetrizing the initial barostat velocities to ensure
885 ! no rotation of the cell under NPT_F
886 npt(j, i)%v = v
887 npt(i, j)%v = v
888 END DO
889 END DO
890
891 akin = 0.0_dp
892 DO i = 1, SIZE(npt, 1)
893 DO j = 1, SIZE(npt, 2)
894 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 number = SIZE(npt, 1)*SIZE(npt, 2)
899
900 ! scale velocities to get the correct initial temperature
901 IF (number /= 0) THEN
902 temp = 2.0_dp*akin/real(number, kind=dp)
903 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
904 END IF
905 DO i = 1, SIZE(npt, 1)
906 DO j = i, SIZE(npt, 2)
907 npt(j, i)%v = temp*npt(j, i)%v
908 npt(i, j)%v = npt(j, i)%v
909 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 npt(:, :)%v = 0.0_dp
922 END SELECT
923
924 CALL timestop(handle)
925
926 END SUBROUTINE init_barostat_variables
927
928! **************************************************************************************************
929!> \brief Assigns extended parameters from the restart file.
930!> \param nhc ...
931!> \author CJM
932! **************************************************************************************************
933 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 CALL timeset(routinen, handle)
942
943 cpassert(ASSOCIATED(nhc))
944 ! assign the forces
945 DO i = 1, SIZE(nhc%nvt, 2)
946 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 nhc%nvt(j, i)%nkt
950 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
951 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 CALL timestop(handle)
957
958 END SUBROUTINE init_nhc_forces
959
960END MODULE extended_system_init
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
extended_system_init
Definition extended_system_init.F:15
extended_system_init::initialize_npt
subroutine, public initialize_npt(simpar, globenv, npt_info, cell, work_section)
...
Definition extended_system_init.F:74
extended_system_init::initialize_nhc_baro
subroutine, public initialize_nhc_baro(simpar, para_env, globenv, nhc, nose_section, save_mem)
fire up the thermostats, if NPT
Definition extended_system_init.F:176
extended_system_init::initialize_nhc_slow
subroutine, public initialize_nhc_slow(thermostat_info, simpar, local_molecules, molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, gci, save_mem)
...
Definition extended_system_init.F:240
extended_system_init::initialize_nhc_shell
subroutine, public initialize_nhc_shell(thermostat_info, simpar, local_molecules, molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, gci, save_mem, binary_restart_file_name)
...
Definition extended_system_init.F:442
extended_system_init::initialize_nhc_fast
subroutine, public initialize_nhc_fast(thermostat_info, simpar, local_molecules, molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, gci, save_mem)
...
Definition extended_system_init.F:306
extended_system_init::initialize_nhc_part
subroutine, public initialize_nhc_part(thermostat_info, simpar, local_molecules, molecule, molecule_kind_set, para_env, globenv, nhc, nose_section, gci, save_mem, binary_restart_file_name)
...
Definition extended_system_init.F:373
extended_system_mapping
Definition extended_system_mapping.F:15
extended_system_mapping::nhc_to_particle_mapping_fast
subroutine, public nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, molecule_set, molecule_kind_set, nhc, para_env, gci)
Creates the thermostatting maps.
Definition extended_system_mapping.F:396
extended_system_mapping::nhc_to_particle_mapping_slow
subroutine, public nhc_to_particle_mapping_slow(thermostat_info, simpar, local_molecules, molecule_set, molecule_kind_set, nhc, para_env, gci)
Creates the thermostatting maps.
Definition extended_system_mapping.F:300
extended_system_mapping::nhc_to_barostat_mapping
subroutine, public nhc_to_barostat_mapping(simpar, nhc)
Creates the thermostatting for the barostat.
Definition extended_system_mapping.F:60
extended_system_mapping::nhc_to_shell_mapping
subroutine, public nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, molecule_set, molecule_kind_set, nhc, para_env, gci)
...
Definition extended_system_mapping.F:559
extended_system_mapping::nhc_to_particle_mapping
subroutine, public nhc_to_particle_mapping(thermostat_info, simpar, local_molecules, molecule_set, molecule_kind_set, nhc, para_env, gci)
Creates the thermostatting maps.
Definition extended_system_mapping.F:134
extended_system_types
Lumps all possible extended system variables into one type for easy access and passing.
Definition extended_system_types.F:17
extended_system_types::debug_isotropic_limit
logical, parameter, public debug_isotropic_limit
Definition extended_system_types.F:38
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_thermo_only_master
integer, parameter, public do_thermo_only_master
Definition input_constants.F:306
input_constants::nph_uniaxial_ensemble
integer, parameter, public nph_uniaxial_ensemble
Definition input_constants.F:102
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:102
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:102
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:102
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:102
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:102
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:102
input_cp2k_binary_restarts
Routines to read the binary restart file of CP2K.
Definition input_cp2k_binary_restarts.F:15
input_cp2k_binary_restarts::read_binary_thermostats_nose
subroutine, public read_binary_thermostats_nose(prefix, nhc, binary_restart_file_name, restart, para_env)
Read the input section &THERMOSTAT for Nose thermostats from an external file written in binary forma...
Definition input_cp2k_binary_restarts.F:694
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_remove_values
subroutine, public section_vals_remove_values(section_vals)
removes the values of a repetition of the section
Definition input_section_types.F:901
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
thermostat_types
Thermostat structure: module containing thermostat available for MD.
Definition thermostat_types.F:12
thermostat_utils
Utilities for thermostats.
Definition thermostat_utils.F:12
thermostat_utils::get_nhc_energies
subroutine, public get_nhc_energies(nhc, nhc_pot, nhc_kin, para_env, array_kin, array_pot)
Calculates kinetic energy and potential energy of the nhc variables.
Definition thermostat_utils.F:1462
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
extended_system_types::lnhc_parameters_type
Definition extended_system_types.F:75
extended_system_types::map_info_type
Definition extended_system_types.F:65
extended_system_types::npt_info_type
Definition extended_system_types.F:42
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:151
molecule_types::global_constraint_type
Definition molecule_types.F:134
molecule_types::molecule_type
Definition molecule_types.F:159
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
thermostat_types::thermostat_info_type
Definition thermostat_types.F:90