(git:374b731)
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-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! **************************************************************************************************
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) &
130 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
131 "MASS section (or none) in the BAROSTAT section")
132 restart = explicit .AND. restart
133 END IF
134
135 IF (restart) THEN
136 CALL section_vals_val_get(work_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
137 ind = 0
138 DO i = 1, SIZE(npt_info, 1)
139 DO j = 1, SIZE(npt_info, 2)
140 ind = ind + 1
141 npt_info(i, j)%v = buffer(ind)
142 END DO
143 END DO
144 CALL section_vals_val_get(work_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
145 ind = 0
146 DO i = 1, SIZE(npt_info, 1)
147 DO j = 1, SIZE(npt_info, 2)
148 ind = ind + 1
149 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 globenv)
156 END IF
157
158 END IF
159
160 CALL timestop(handle)
161
162 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 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 CALL timeset(routinen, handle)
190
191 restart = .false.
192
193 CALL nhc_to_barostat_mapping(simpar, nhc)
194
195 ! Set up the Yoshida weights
196 IF (nhc%nyosh > 0) THEN
197 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
198 CALL set_yoshida_coef(nhc, simpar%dt)
199 END IF
200
201 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
202
203 IF (.NOT. restart) THEN
204 ! Initializing thermostat forces and velocities for the Nose-Hoover
205 ! Chain variables
206 SELECT CASE (simpar%ensemble)
207 CASE DEFAULT
208 temp = simpar%temp_baro_ext
209 END SELECT
210 IF (nhc%nhc_len /= 0) THEN
211 CALL init_nhc_variables(nhc, temp, para_env, globenv)
212 END IF
213 END IF
214
215 CALL init_nhc_forces(nhc)
216
217 CALL timestop(handle)
218
219 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 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 CALL timeset(routinen, handle)
258
259 restart = .false.
260 ! fire up the thermostats, if not NVE
261
262 CALL nhc_to_particle_mapping_slow(thermostat_info, simpar, local_molecules, &
263 molecule, molecule_kind_set, nhc, para_env, gci)
264
265 ! Set up the Yoshida weights
266 IF (nhc%nyosh > 0) THEN
267 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
268 CALL set_yoshida_coef(nhc, simpar%dt)
269 END IF
270
271 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
272
273 IF (.NOT. restart) THEN
274 ! Initializing thermostat forces and velocities for the Nose-Hoover
275 ! Chain variables
276 IF (nhc%nhc_len /= 0) THEN
277 CALL init_nhc_variables(nhc, simpar%temp_slow, para_env, globenv)
278 END IF
279 END IF
280
281 CALL init_nhc_forces(nhc)
282
283 CALL timestop(handle)
284
285 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 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 CALL timeset(routinen, handle)
324
325 restart = .false.
326 ! fire up the thermostats, if not NVE
327
328 CALL nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, &
329 molecule, molecule_kind_set, nhc, para_env, gci)
330
331 ! Set up the Yoshida weights
332 IF (nhc%nyosh > 0) THEN
333 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
334 CALL set_yoshida_coef(nhc, simpar%dt)
335 END IF
336
337 CALL restart_nose(nhc, nose_section, save_mem, restart, "", "", para_env)
338
339 IF (.NOT. restart) THEN
340 ! Initializing thermostat forces and velocities for the Nose-Hoover
341 ! Chain variables
342 IF (nhc%nhc_len /= 0) THEN
343 CALL init_nhc_variables(nhc, simpar%temp_fast, para_env, globenv)
344 END IF
345 END IF
346
347 CALL init_nhc_forces(nhc)
348
349 CALL timestop(handle)
350
351 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 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 CALL timeset(routinen, handle)
392
393 restart = .false.
394 ! fire up the thermostats, if not NVE
395
396 CALL nhc_to_particle_mapping(thermostat_info, simpar, local_molecules, &
397 molecule, molecule_kind_set, nhc, para_env, gci)
398
399 ! Set up the Yoshida weights
400 IF (nhc%nyosh > 0) THEN
401 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
402 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 "PARTICLE", para_env)
407
408 IF (.NOT. restart) THEN
409 ! Initializing thermostat forces and velocities for the Nose-Hoover
410 ! Chain variables
411 IF (nhc%nhc_len /= 0) THEN
412 CALL init_nhc_variables(nhc, simpar%temp_ext, para_env, globenv)
413 END IF
414 END IF
415
416 CALL init_nhc_forces(nhc)
417
418 CALL timestop(handle)
419
420 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 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 CALL timeset(routinen, handle)
461
462 CALL nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, &
463 molecule, molecule_kind_set, nhc, para_env, gci)
464
465 restart = .false.
466 ! Set up the Yoshida weights
467 IF (nhc%nyosh > 0) THEN
468 ALLOCATE (nhc%dt_yosh(1:nhc%nyosh))
469 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 "SHELL", para_env)
474
475 IF (.NOT. restart) THEN
476 ! Initialize thermostat forces and velocities
477 ! Chain variables
478 IF (nhc%nhc_len /= 0) THEN
479 CALL init_nhc_variables(nhc, simpar%temp_sh_ext, para_env, globenv)
480 END IF
481 END IF
482
483 CALL init_nhc_forces(nhc)
484
485 CALL timestop(handle)
486
487 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 SUBROUTINE set_yoshida_coef(nhc, dt)
499
500 TYPE(lnhc_parameters_type), POINTER :: nhc
501 REAL(kind=dp), INTENT(IN) :: dt
502
503 REAL(kind=dp), DIMENSION(nhc%nyosh) :: yosh_wt
504
505 SELECT CASE (nhc%nyosh)
506 CASE DEFAULT
507 cpabort('Value not available.')
508 CASE (1)
509 yosh_wt(1) = 1.0_dp
510 CASE (3)
511 yosh_wt(1) = 1.0_dp/(2.0_dp - (2.0_dp)**(1.0_dp/3.0_dp))
512 yosh_wt(2) = 1.0_dp - 2.0_dp*yosh_wt(1)
513 yosh_wt(3) = yosh_wt(1)
514 CASE (5)
515 yosh_wt(1) = 1.0_dp/(4.0_dp - (4.0_dp)**(1.0_dp/3.0_dp))
516 yosh_wt(2) = yosh_wt(1)
517 yosh_wt(4) = yosh_wt(1)
518 yosh_wt(5) = yosh_wt(1)
519 yosh_wt(3) = 1.0_dp - 4.0_dp*yosh_wt(1)
520 CASE (7)
521 yosh_wt(1) = .78451361047756_dp
522 yosh_wt(2) = .235573213359357_dp
523 yosh_wt(3) = -1.17767998417887_dp
524 yosh_wt(4) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + yosh_wt(3))
525 yosh_wt(5) = yosh_wt(3)
526 yosh_wt(6) = yosh_wt(2)
527 yosh_wt(7) = yosh_wt(1)
528 CASE (9)
529 yosh_wt(1) = 0.192_dp
530 yosh_wt(2) = 0.554910818409783619692725006662999_dp
531 yosh_wt(3) = 0.124659619941888644216504240951585_dp
532 yosh_wt(4) = -0.843182063596933505315033808282941_dp
533 yosh_wt(5) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
534 yosh_wt(3) + yosh_wt(4))
535 yosh_wt(6) = yosh_wt(4)
536 yosh_wt(7) = yosh_wt(3)
537 yosh_wt(8) = yosh_wt(2)
538 yosh_wt(9) = yosh_wt(1)
539 CASE (15)
540 yosh_wt(1) = 0.102799849391985_dp
541 yosh_wt(2) = -0.196061023297549e1_dp
542 yosh_wt(3) = 0.193813913762276e1_dp
543 yosh_wt(4) = -0.158240635368243_dp
544 yosh_wt(5) = -0.144485223686048e1_dp
545 yosh_wt(6) = 0.253693336566229_dp
546 yosh_wt(7) = 0.914844246229740_dp
547 yosh_wt(8) = 1.0_dp - 2.0_dp*(yosh_wt(1) + yosh_wt(2) + &
548 yosh_wt(3) + yosh_wt(4) + yosh_wt(5) + yosh_wt(6) + yosh_wt(7))
549 yosh_wt(9) = yosh_wt(7)
550 yosh_wt(10) = yosh_wt(6)
551 yosh_wt(11) = yosh_wt(5)
552 yosh_wt(12) = yosh_wt(4)
553 yosh_wt(13) = yosh_wt(3)
554 yosh_wt(14) = yosh_wt(2)
555 yosh_wt(15) = yosh_wt(1)
556 END SELECT
557 nhc%dt_yosh = dt*yosh_wt/real(nhc%nc, kind=dp)
558
559 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: buffer
591 TYPE(map_info_type), POINTER :: map_info
592 TYPE(section_vals_type), POINTER :: work_section
593
594 CALL timeset(routinen, handle)
595
596 NULLIFY (buffer)
597 NULLIFY (work_section)
598
599 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 restart, para_env)
605
606 ELSE
607
608 ! Read the default restart file in ASCII format
609
610 explicit = .false.
611 restart = .false.
612
613 IF (ASSOCIATED(nose_section)) THEN
614 work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
615 CALL section_vals_get(work_section, explicit=explicit)
616 restart = explicit
617 work_section => section_vals_get_subs_vals(nose_section, "COORD")
618 CALL section_vals_get(work_section, explicit=explicit)
619 IF (.NOT. restart .AND. explicit) &
620 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
621 "COORD and MASS and FORCE section (or none) in the NOSE section")
622 restart = explicit .AND. restart
623 work_section => section_vals_get_subs_vals(nose_section, "MASS")
624 CALL section_vals_get(work_section, explicit=explicit)
625 IF (.NOT. restart .AND. explicit) &
626 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
627 "COORD and MASS and FORCE section (or none) in the NOSE section")
628 restart = explicit .AND. restart
629 work_section => section_vals_get_subs_vals(nose_section, "FORCE")
630 CALL section_vals_get(work_section, explicit=explicit)
631 IF (.NOT. restart .AND. explicit) &
632 CALL cp_abort(__location__, "You need to define both VELOCITY and "// &
633 "COORD and MASS and FORCE section (or none) in the NOSE section")
634 restart = explicit .AND. restart
635 END IF
636
637 IF (restart) THEN
638 map_info => nhc%map_info
639 CALL section_vals_val_get(nose_section, "COORD%_DEFAULT_KEYWORD_", r_vals=buffer)
640 DO i = 1, SIZE(nhc%nvt, 2)
641 ind = map_info%index(i)
642 ind = (ind - 1)*nhc%nhc_len
643 DO j = 1, SIZE(nhc%nvt, 1)
644 ind = ind + 1
645 nhc%nvt(j, i)%eta = buffer(ind)
646 END DO
647 END DO
648 CALL section_vals_val_get(nose_section, "VELOCITY%_DEFAULT_KEYWORD_", r_vals=buffer)
649 DO i = 1, SIZE(nhc%nvt, 2)
650 ind = map_info%index(i)
651 ind = (ind - 1)*nhc%nhc_len
652 DO j = 1, SIZE(nhc%nvt, 1)
653 ind = ind + 1
654 nhc%nvt(j, i)%v = buffer(ind)
655 END DO
656 END DO
657 CALL section_vals_val_get(nose_section, "MASS%_DEFAULT_KEYWORD_", r_vals=buffer)
658 DO i = 1, SIZE(nhc%nvt, 2)
659 ind = map_info%index(i)
660 ind = (ind - 1)*nhc%nhc_len
661 DO j = 1, SIZE(nhc%nvt, 1)
662 ind = ind + 1
663 nhc%nvt(j, i)%mass = buffer(ind)
664 END DO
665 END DO
666 CALL section_vals_val_get(nose_section, "FORCE%_DEFAULT_KEYWORD_", r_vals=buffer)
667 DO i = 1, SIZE(nhc%nvt, 2)
668 ind = map_info%index(i)
669 ind = (ind - 1)*nhc%nhc_len
670 DO j = 1, SIZE(nhc%nvt, 1)
671 ind = ind + 1
672 nhc%nvt(j, i)%f = buffer(ind)
673 END DO
674 END DO
675 END IF
676
677 IF (save_mem) THEN
678 NULLIFY (work_section)
679 work_section => section_vals_get_subs_vals(nose_section, "COORD")
680 CALL section_vals_remove_values(work_section)
681 NULLIFY (work_section)
682 work_section => section_vals_get_subs_vals(nose_section, "VELOCITY")
683 CALL section_vals_remove_values(work_section)
684 NULLIFY (work_section)
685 work_section => section_vals_get_subs_vals(nose_section, "FORCE")
686 CALL section_vals_remove_values(work_section)
687 NULLIFY (work_section)
688 work_section => section_vals_get_subs_vals(nose_section, "MASS")
689 CALL section_vals_remove_values(work_section)
690 END IF
691
692 END IF
693
694 CALL timestop(handle)
695
696 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: array_of_rn
719 TYPE(map_info_type), POINTER :: map_info
720
721 CALL timeset(routinen, handle)
722
723 tot_rn = 0
724
725 ! first initializing the mass of the nhc variables
726 nhc%nvt(:, :)%mass = nhc%nvt(:, :)%nkt*nhc%tau_nhc**2
727 nhc%nvt(:, :)%eta = 0._dp
728 nhc%nvt(:, :)%v = 0._dp
729 nhc%nvt(:, :)%f = 0._dp
730
731 map_info => nhc%map_info
732 SELECT CASE (map_info%dis_type)
733 CASE (do_thermo_only_master) ! for NPT
734 CASE DEFAULT
735 tot_rn = nhc%glob_num_nhc*nhc%nhc_len
736
737 ALLOCATE (array_of_rn(tot_rn))
738 array_of_rn(:) = 0.0_dp
739 END SELECT
740
741 SELECT CASE (map_info%dis_type)
742 CASE (do_thermo_only_master) ! for NPT
743 ! Map deterministically determined random number to nhc % v
744 DO i = 1, nhc%loc_num_nhc
745 DO j = 1, nhc%nhc_len
746 nhc%nvt(j, i)%v = globenv%gaussian_rng_stream%next()
747 END DO
748 END DO
749
750 akin = 0.0_dp
751 DO i = 1, nhc%loc_num_nhc
752 DO j = 1, nhc%nhc_len
753 akin = akin + 0.5_dp*(nhc%nvt(j, i)%mass* &
754 nhc%nvt(j, i)%v* &
755 nhc%nvt(j, i)%v)
756 END DO
757 END DO
758 number = nhc%loc_num_nhc
759
760 ! scale velocities to get the correct initial temperature
761 temp = 2.0_dp*akin/real(number, kind=dp)
762 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
763 DO i = 1, nhc%loc_num_nhc
764 DO j = 1, nhc%nhc_len
765 nhc%nvt(j, i)%v = temp*nhc%nvt(j, i)%v
766 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 DO i = 1, nhc%loc_num_nhc
772 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 nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
775 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
776 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 DO i = 1, tot_rn
783 array_of_rn(i) = globenv%gaussian_rng_stream%next()
784 END DO
785 ! Map deterministically determined random number to nhc % v
786 DO i = 1, nhc%loc_num_nhc
787 icount = map_info%index(i)
788 icount = (icount - 1)*nhc%nhc_len
789 DO j = 1, nhc%nhc_len
790 icount = icount + 1
791 nhc%nvt(j, i)%v = array_of_rn(icount)
792 ! WRITE ( *, * ) 'VEL', para_env%mepos, i,j, nhc%nvt(j,i)%v
793 nhc%nvt(j, i)%eta = 0.0_dp
794 END DO
795 END DO
796 DEALLOCATE (array_of_rn)
797
798 number = nhc%glob_num_nhc
799 CALL get_nhc_energies(nhc, dum, akin, para_env)
800
801 ! scale velocities to get the correct initial temperature
802 temp = 2.0_dp*akin/real(number, kind=dp)
803 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
804 DO i = 1, nhc%loc_num_nhc
805 DO j = 1, nhc%nhc_len
806 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 DO i = 1, nhc%loc_num_nhc
812 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 nhc%nvt(j - 1, i)%v - nhc%nvt(j, i)%nkt
815 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
816 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 CALL timestop(handle)
824
825 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 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 CALL timeset(routinen, handle)
856
857 temp = 0.0_dp
858
859 ! first initializing the mass of the nhc variables
860
861 npt(:, :)%eps = 0.0_dp
862 npt(:, :)%v = 0.0_dp
863 npt(:, :)%f = 0.0_dp
864 SELECT CASE (ensemble)
865 CASE (npt_i_ensemble, npt_ia_ensemble)
866 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2
867 CASE (npt_f_ensemble)
868 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 npt(:, :)%mass = cmass
871 CASE (npe_f_ensemble)
872 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2/3.0_dp
873 CASE (npe_i_ensemble)
874 npt(:, :)%mass = real(nfree + 3, kind=dp)*temp_ext*tau_cell**2
875 END SELECT
876 ! initializing velocities
877 DO i = 1, SIZE(npt, 1)
878 DO j = i, SIZE(npt, 2)
879 v = globenv%gaussian_rng_stream%next()
880 ! Symmetrizing the initial barostat velocities to ensure
881 ! no rotation of the cell under NPT_F
882 npt(j, i)%v = v
883 npt(i, j)%v = v
884 END DO
885 END DO
886
887 akin = 0.0_dp
888 DO i = 1, SIZE(npt, 1)
889 DO j = 1, SIZE(npt, 2)
890 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 number = SIZE(npt, 1)*SIZE(npt, 2)
895
896 ! scale velocities to get the correct initial temperature
897 IF (number /= 0) THEN
898 temp = 2.0_dp*akin/real(number, kind=dp)
899 IF (temp > 0.0_dp) temp = sqrt(temp_ext/temp)
900 END IF
901 DO i = 1, SIZE(npt, 1)
902 DO j = i, SIZE(npt, 2)
903 npt(j, i)%v = temp*npt(j, i)%v
904 npt(i, j)%v = npt(j, i)%v
905 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 npt(:, :)%v = 0.0_dp
918 END SELECT
919
920 CALL timestop(handle)
921
922 END SUBROUTINE init_barostat_variables
923
924! **************************************************************************************************
925!> \brief Assigns extended parameters from the restart file.
926!> \param nhc ...
927!> \author CJM
928! **************************************************************************************************
929 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 CALL timeset(routinen, handle)
938
939 cpassert(ASSOCIATED(nhc))
940 ! assign the forces
941 DO i = 1, SIZE(nhc%nvt, 2)
942 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 nhc%nvt(j, i)%nkt
946 IF (nhc%nvt(j, i)%mass > 0.0_dp) THEN
947 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 CALL timestop(handle)
953
954 END SUBROUTINE init_nhc_forces
955
956END 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:175
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:239
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:441
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:305
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:372
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:395
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:557
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:17
input_constants::do_thermo_only_master
integer, parameter, public do_thermo_only_master
Definition input_constants.F:296
input_constants::nph_uniaxial_ensemble
integer, parameter, public nph_uniaxial_ensemble
Definition input_constants.F:95
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:95
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:95
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:95
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:95
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:95
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:95
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:677
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:889
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:724
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:697
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:1040
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:1379
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
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:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::molecule_type
Definition molecule_types.F:110
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