(git:374b731)
Loading...
Searching...
No Matches
pint_piglet.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!> \brief Methods to apply the piglet thermostat to PI runs.
10!> \author Felix Uhl
11!> \par History
12!> 10.2014 created [Felix Uhl]
13! **************************************************************************************************
14MODULE pint_piglet
15 USE cp_files, ONLY: close_file,&
16 open_file
17 USE cp_units, ONLY: cp_unit_from_cp2k,&
18 cp_unit_to_cp2k
19 USE gle_system_dynamics, ONLY: gle_cholesky_stab,&
20 gle_matrix_exp
21 USE input_constants, ONLY: matrix_init_cholesky,&
22 matrix_init_diagonal,&
23 propagator_rpmd
24 USE input_section_types, ONLY: section_vals_get,&
25 section_vals_get_subs_vals,&
26 section_vals_type,&
27 section_vals_val_get
28 USE kinds, ONLY: default_path_length,&
29 default_string_length,&
30 dp
31 USE message_passing, ONLY: mp_para_env_type
32 USE parallel_rng_types, ONLY: gaussian,&
33 rng_record_length,&
34 rng_stream_type,&
35 rng_stream_type_from_record
36 USE pint_io, ONLY: pint_write_line
37 USE pint_types, ONLY: piglet_therm_type,&
38 pint_env_type
39#include "../base/base_uses.f90"
40
41 IMPLICIT NONE
42
43 PRIVATE
44
45 PUBLIC :: pint_piglet_create, &
46 pint_piglet_init, &
47 pint_piglet_step, &
48 pint_piglet_release, &
49 pint_calc_piglet_energy
50
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_piglet'
52
53!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
54! !
55! _._ _..._ .-', _.._ _ !
56! '-. ` ' /-._.-' ',/ ) !
57! ) \ '. !
58! / _ _ | \ !
59! | o o / | !
60! \ .-. ; !
61! '-('' ).-' ,' ; !
62! '-; | .' !
63! \ \ / !
64! | 7 .__ _.-\ \ !
65! | | | ``/ /` / !
66! /,_| | /,_/ / !
67! /,_/ '`-' !
68! !
69! ( ( ( !
70! )\ ) )\ ) ( )\ ) * ) !
71! (()/((()/( )\ ) (()/( ( ` ) /( !
72! /(_))/(_))(()/( /(_)) )\ ( )(_)) !
73! (_)) (_)) /(_))_ (_)) ((_) (_(_()) !
74! | _ \|_ _| (_)) __|| | | __||_ _| !
75! | _/ | | | (_ || |__ | _| | | !
76! |_| |___| \___||____||___| |_| !
77! !
78! Make Quantum Mechanics Hot !
79! !
80!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
81
82CONTAINS
83
84! ***************************************************************************
85!> \brief creates the data structure for a piglet thermostating in PI runs
86!> \param piglet_therm ...
87!> \param pint_env ...
88!> \param section ...
89!> \author Felix Uhl
90! **************************************************************************************************
91 SUBROUTINE pint_piglet_create(piglet_therm, pint_env, section)
92 TYPE(piglet_therm_type), INTENT(OUT) :: piglet_therm
93 TYPE(pint_env_type), INTENT(IN) :: pint_env
94 TYPE(section_vals_type), POINTER :: section
95
96 INTEGER :: ndim, p
97
98 CALL section_vals_val_get(section, "NEXTRA_DOF", i_val=piglet_therm%nsp1)
99 !add real degree of freedom to ns to reach nsp1
100 piglet_therm%nsp1 = piglet_therm%nsp1 + 1
101 p = pint_env%p
102 piglet_therm%p = pint_env%p
103 ndim = pint_env%ndim
104 piglet_therm%ndim = pint_env%ndim
105 NULLIFY (piglet_therm%a_mat)
106 ALLOCATE (piglet_therm%a_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
107 piglet_therm%a_mat(:, :, :) = 0.0_dp
108 NULLIFY (piglet_therm%c_mat)
109 ALLOCATE (piglet_therm%c_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
110 piglet_therm%c_mat(:, :, :) = 0.0_dp
111 NULLIFY (piglet_therm%gle_t)
112 ALLOCATE (piglet_therm%gle_t(piglet_therm%nsp1, piglet_therm%nsp1, p))
113 piglet_therm%gle_t(:, :, :) = 0.0_dp
114 NULLIFY (piglet_therm%gle_s)
115 ALLOCATE (piglet_therm%gle_s(piglet_therm%nsp1, piglet_therm%nsp1, p))
116 piglet_therm%gle_s(:, :, :) = 0.0_dp
117 NULLIFY (piglet_therm%smalls)
118 ALLOCATE (piglet_therm%smalls(piglet_therm%nsp1, ndim*p))
119 piglet_therm%smalls(:, :) = 0.0_dp
120 NULLIFY (piglet_therm%temp1)
121 ALLOCATE (piglet_therm%temp1(piglet_therm%nsp1, ndim))
122 piglet_therm%temp1(:, :) = 0.0_dp
123 NULLIFY (piglet_therm%temp2)
124 ALLOCATE (piglet_therm%temp2(piglet_therm%nsp1, ndim))
125 piglet_therm%temp2(:, :) = 0.0_dp
126 NULLIFY (piglet_therm%sqrtmass)
127 ALLOCATE (piglet_therm%sqrtmass(piglet_therm%p, piglet_therm%ndim))
128 piglet_therm%sqrtmass(:, :) = 0.0_dp
129
130 END SUBROUTINE pint_piglet_create
131
132! ***************************************************************************
133!> \brief initializes the data for a piglet run
134!> \param piglet_therm ...
135!> \param pint_env ...
136!> \param section ...
137!> \param dt ...
138!> \param para_env ...
139!> \author Felix Uhl
140! **************************************************************************************************
141 SUBROUTINE pint_piglet_init(piglet_therm, pint_env, section, dt, para_env)
142 TYPE(piglet_therm_type), INTENT(INOUT) :: piglet_therm
143 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
144 TYPE(section_vals_type), POINTER :: section
145 REAL(kind=dp), INTENT(IN) :: dt
146 TYPE(mp_para_env_type), POINTER :: para_env
147
148 CHARACTER(len=20) :: default_format, read_unit
149 CHARACTER(len=default_path_length) :: matrices_file_name
150 CHARACTER(len=default_string_length) :: msg, temp_input
151 CHARACTER(LEN=rng_record_length) :: rng_record
152 INTEGER :: cbrac, i, ibead, idim, imode, &
153 input_unit, isp, j, matrix_init, ns, &
154 obrac, p, read_err
155 LOGICAL :: explicit
156 REAL(kind=dp), DIMENSION(3, 2) :: initial_seed
157 REAL(kind=dp), DIMENSION(:), POINTER :: smallstmp
158 REAL(kind=dp), &
159 DIMENSION(piglet_therm%nsp1, piglet_therm%nsp1) :: mtmp, tmpmatrix
160 TYPE(section_vals_type), POINTER :: rng_section, smalls_section
161
162 p = piglet_therm%p
163 pint_env%e_piglet = 0.0_dp
164 pint_env%piglet_therm%thermostat_energy = 0.0_dp
165 CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=piglet_therm%thermostat_energy)
166 CALL section_vals_val_get(section, "SMATRIX_INIT", i_val=matrix_init)
167 !Read the matices
168
169 IF (pint_env%propagator%prop_kind /= propagator_rpmd) THEN
170 cpabort("PIGLET is designed to work with the RPMD propagator")
171 END IF
172 ! Select algorithm for S-matrix initialization
173 IF (matrix_init == matrix_init_cholesky) THEN
174 IF (para_env%is_source()) THEN
175 CALL pint_write_line("PIGLET| Initalizing S-matrices using cholesky decomposition.")
176 END IF
177 ELSE IF (matrix_init == matrix_init_diagonal) THEN
178 IF (para_env%is_source()) THEN
179 CALL pint_write_line("PIGLET| Initalizing S-matrices using full diagonalization.")
180 END IF
181 ELSE
182 cpwarn("No PIGLET init algorithm selected. Selecting cholesky decomposition")
183 matrix_init = matrix_init_cholesky
184 END IF
185
186 IF (para_env%is_source()) THEN
187 ! Read input matrices
188 WRITE (default_format, '(A,I10,A)') "(A", default_string_length, ")"
189 CALL section_vals_val_get(section, "MATRICES_FILE_NAME", &
190 c_val=matrices_file_name)
191 CALL pint_write_line("PIGLET| Reading PIGLET matrices from file: ")
192 CALL pint_write_line("PIGLET| "//trim(matrices_file_name))
193 CALL open_file(file_name=trim(matrices_file_name), &
194 file_action="READ", file_status="OLD", unit_number=input_unit)
195 read_err = 0
196 msg = ""
197 DO WHILE (read_err == 0)
198 READ (input_unit, default_format, iostat=read_err) temp_input
199 IF (read_err /= 0) EXIT
200 !Parse comment section
201 IF (index(temp_input, "#") /= 0) THEN
202 IF (index(temp_input, "T=") /= 0) THEN
203 CALL check_temperature(line=temp_input, &
204 propagator=pint_env%propagator%prop_kind, &
205 targettemp=pint_env%kT*pint_env%propagator%temp_sim2phys)
206 ELSE IF (index(temp_input, "A MATRIX") /= 0) THEN
207 obrac = index(temp_input, "(") + 1
208 cbrac = index(temp_input, ")") - 1
209 read_unit = temp_input(obrac:cbrac)
210 DO imode = 1, p
211 READ (input_unit, default_format) temp_input
212 DO i = 1, piglet_therm%nsp1
213 READ (input_unit, *, iostat=read_err) &
214 (piglet_therm%a_mat(i, j, imode), j=1, piglet_therm%nsp1)
215 IF (read_err /= 0) THEN
216 WRITE (unit=msg, fmt=*) "Invalid PIGLET A-matrix Nr.", i - 1
217 cpabort(msg)
218 EXIT
219 END IF
220 END DO
221 END DO
222 !convert to cp2k units
223 IF (read_err == 0) THEN
224 CALL a_mat_to_cp2k(piglet_therm%a_mat, p, &
225 piglet_therm%nsp1, read_unit, msg)
226 END IF
227 ELSE IF (index(temp_input, "C MATRIX") /= 0) THEN
228 obrac = index(temp_input, "(") + 1
229 cbrac = index(temp_input, ")") - 1
230 read_unit = temp_input(obrac:cbrac)
231 DO imode = 1, p
232 READ (input_unit, default_format) temp_input
233 DO i = 1, piglet_therm%nsp1
234 READ (input_unit, *, iostat=read_err) &
235 (piglet_therm%c_mat(i, j, imode), j=1, piglet_therm%nsp1)
236 IF (read_err /= 0) THEN
237 WRITE (unit=msg, fmt=*) "Invalid PIGLET C-matrix Nr.", i - 1
238 cpabort(msg)
239 EXIT
240 END IF
241 END DO
242 END DO
243 IF (read_err == 0) THEN
244 CALL c_mat_to_cp2k(piglet_therm%c_mat, p, &
245 piglet_therm%nsp1, read_unit, msg)
246 END IF
247 END IF
248 END IF
249 END DO
250 CALL close_file(unit_number=input_unit)
251 END IF
252 ! communicate A and C matrix to other nodes
253 CALL para_env%bcast(piglet_therm%a_mat, &
254 para_env%source)
255 CALL para_env%bcast(piglet_therm%c_mat, &
256 para_env%source)
257
258 !prepare Random number generator
259 NULLIFY (rng_section)
260 rng_section => section_vals_get_subs_vals(section, &
261 subsection_name="RNG_INIT")
262 CALL section_vals_get(rng_section, explicit=explicit)
263 IF (explicit) THEN
264 CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
265 i_rep_val=1, c_val=rng_record)
266 piglet_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
267 ELSE
268 initial_seed(:, :) = real(pint_env%thermostat_rng_seed, dp)
269 piglet_therm%gaussian_rng_stream = rng_stream_type( &
270 name="piglet_rng_gaussian", distribution_type=gaussian, &
271 extended_precision=.true., &
272 seed=initial_seed)
273 END IF
274
275 !Compute the T and S matrices on every mpi process
276 DO i = 1, p
277 ! T = EXP(-A_mat*dt) = EXP(-A_mat*dt/2)
278 ! Values for j and k = 15 are way to high, but to be sure.
279 ! (its only executed once anyway)
280 CALL gle_matrix_exp(m=(-0.5_dp*dt)*piglet_therm%a_mat(:, :, i), & !dt scaled A matrix
281 n=piglet_therm%nsp1, & !size of matrix
282 j=15, & !truncation for taylor expansion
283 k=15, & !scaling parameter for faster convergence of taylor expansion
284 em=piglet_therm%gle_t(:, :, i)) !output T matrices
285 ! S*TRANSPOSE(S) = C-T*C*TRANSPOSE(T)
286 ! T*C:
287 CALL dgemm('N', & ! T-matrix is not transposed
288 'N', & ! C-matrix is not transposed
289 piglet_therm%nsp1, & ! number of rows of T-matrix
290 piglet_therm%nsp1, & ! number of columns of C-matrix
291 piglet_therm%nsp1, & ! number of columns of T-matrix and number of rows of C-matrix
292 1.0d0, & ! scaling factor alpha
293 piglet_therm%gle_t(:, :, i), & ! T-matrix
294 piglet_therm%nsp1, & ! leading dimension of T-matrix as declared
295 piglet_therm%c_mat(:, :, i), & ! C-matrix
296 piglet_therm%nsp1, & ! leading dimension of C-matrix as declared
297 0.0d0, & ! scaling of tmpmatrix as additive
298 tmpmatrix, & ! result matrix: tmpmatrix
299 piglet_therm%nsp1) ! leading dimension of tmpmatrix
300 ! T*C*TRANSPOSE(T):
301 CALL dgemm('N', & ! tmpmatrix is not transposed
302 'T', & ! T-matrix is transposed
303 piglet_therm%nsp1, & ! number of rows of tmpmatrix
304 piglet_therm%nsp1, & ! number of columns of T-matrix
305 piglet_therm%nsp1, & ! number of columns of tmpmatrix and number of rows of T-matrix
306 1.0d0, & ! scaling factor alpha
307 tmpmatrix, & ! tmpmatrix
308 piglet_therm%nsp1, & ! leading dimension of tmpmatrix as declared
309 piglet_therm%gle_t(:, :, i), & ! T-matrix
310 piglet_therm%nsp1, & ! leading dimension of T-matrix as declared
311 0.0d0, & ! scaling of Mtmp as additive
312 mtmp, & ! result matrix: Mtmp
313 piglet_therm%nsp1) ! leading dimension of Mtmp
314 ! C - T*C*TRANSPOSE(T):
315 mtmp(:, :) = piglet_therm%c_mat(:, :, i) - mtmp(:, :)
316
317 IF (matrix_init == matrix_init_cholesky) THEN
318 ! Get S by cholesky decomposition of Mtmp
319 CALL gle_cholesky_stab(mtmp, & ! Matrix to decompose
320 piglet_therm%gle_s(:, :, i), & ! result
321 piglet_therm%nsp1) ! Size of the matrix
322 ELSE IF (matrix_init == matrix_init_diagonal) THEN
323 ! Get S by full diagonalization of MTmp matrix
324 CALL sqrt_pos_def_mat(piglet_therm%nsp1, & ! Size of the matrix
325 mtmp, & ! matrix to decompose
326 piglet_therm%gle_s(:, :, i)) ! result
327 END IF
328
329 END DO
330
331 ! Initialize extra degrees of freedom for Markovian Dynamics
332 ! as a cholesky decomposition of C-matrix multiplied by a random number vector
333 ! Or from restart
334 piglet_therm%smalls = 0.0_dp
335 NULLIFY (smalls_section)
336 smalls_section => section_vals_get_subs_vals(section, subsection_name="EXTRA_DOF")
337 CALL section_vals_get(smalls_section, explicit=explicit)
338 IF (explicit) THEN
339 NULLIFY (smallstmp)
340 CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
341 n_rep_val=ns)
342 CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
343 r_vals=smallstmp)
344 i = 1
345 DO isp = 2, piglet_therm%nsp1
346 DO ibead = 1, piglet_therm%p*piglet_therm%ndim
347 piglet_therm%smalls(isp, ibead) = smallstmp(i)
348 i = i + 1
349 END DO
350 END DO
351 ELSE
352 DO ibead = 1, piglet_therm%p
353 IF (matrix_init == matrix_init_cholesky) THEN
354 CALL gle_cholesky_stab(piglet_therm%c_mat(:, :, ibead), & ! Matrix to decompose
355 mtmp, & ! Result
356 piglet_therm%nsp1) ! Size of Matrix
357 ELSE IF (matrix_init == matrix_init_diagonal) THEN
358 ! Get S by full diagonalization of c_mat matrix
359 CALL sqrt_pos_def_mat(piglet_therm%nsp1, & ! Size of the matrix
360 piglet_therm%c_mat(:, :, ibead), & ! matrix to decompose
361 mtmp) ! result
362 END IF
363 ! Fill a vector with random numbers
364 DO idim = 1, piglet_therm%ndim
365 DO j = 1, piglet_therm%nsp1
366 piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
367 !piglet_therm%temp2(j,idim) = 1.0_dp
368 END DO
369 END DO
370 CALL dgemm("N", & ! Matrix Mtmp is not transposed
371 "N", & ! Matrix temp2 is not transposed
372 piglet_therm%nsp1, & ! Number of rows of matrix Mtmp
373 piglet_therm%ndim, & ! Number of columns of matrix temp2
374 piglet_therm%nsp1, & ! Number of columns of matrix Mtmp
375 1.0_dp, & !scaling of Mtmp
376 mtmp(1, 1), & ! Matrix Mtmp
377 piglet_therm%nsp1, & ! leading dimension of Mtmp
378 piglet_therm%temp2, & ! temp2 matrix
379 piglet_therm%nsp1, & ! leading dimension of temp2
380 0.0_dp, & ! scaling of added matrix smalls
381 piglet_therm%temp1, & ! result matrix
382 piglet_therm%nsp1) ! leading dimension of result matrix
383
384 DO idim = 1, piglet_therm%ndim
385 j = (idim - 1)*piglet_therm%p + ibead
386 DO i = 1, piglet_therm%nsp1
387 piglet_therm%smalls(i, j) = piglet_therm%temp1(i, idim)
388 END DO
389 END DO
390 END DO
391 END IF
392
393 !Fill the array for the sqrt of the masses
394 DO idim = 1, piglet_therm%ndim
395 DO ibead = 1, piglet_therm%p
396 piglet_therm%sqrtmass(ibead, idim) = sqrt(pint_env%mass_fict(ibead, idim))
397 END DO
398 END DO
399
400 END SUBROUTINE pint_piglet_init
401
402! **************************************************************************************************
403!> \brief ...
404!> \param vold ...
405!> \param vnew ...
406!> \param first_mode ...
407!> \param masses ...
408!> \param piglet_therm ...
409! **************************************************************************************************
410 SUBROUTINE pint_piglet_step(vold, vnew, first_mode, masses, piglet_therm)
411 REAL(kind=dp), DIMENSION(:, :), POINTER :: vold, vnew
412 INTEGER, INTENT(IN) :: first_mode
413 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: masses
414 TYPE(piglet_therm_type), POINTER :: piglet_therm
415
416 CHARACTER(len=*), PARAMETER :: routinen = 'pint_piglet_step'
417
418 INTEGER :: handle, i, ibead, idim, j, ndim, ndimp, &
419 nsp1, p
420 REAL(kind=dp) :: delta_ekin
421
422 CALL timeset(routinen, handle)
423 nsp1 = piglet_therm%nsp1
424 ndim = piglet_therm%ndim
425 p = piglet_therm%p
426 ndimp = ndim*p
427 ! perform the following operation for all 3N*P
428 ! smalls = gle_t*smalls + gle_s*rand_mat
429 ! Copy mass scaled momenta to temp1 matrix
430 ! p/sqrt(m) we use velocities so v*sqrt(m)
431 DO ibead = first_mode, p
432 ! Copy mass scaled momenta to temp1 matrix
433 ! p/sqrt(m) we use velocities so v*sqrt(m)
434 DO idim = 1, ndim
435 piglet_therm%temp1(1, idim) = vold(ibead, idim)*piglet_therm%sqrtmass(ibead, idim)
436 END DO
437 ! copy the extra degrees of freedom to the temp1 matrix
438 DO idim = 1, ndim
439 DO i = 2, nsp1
440 piglet_therm%temp1(i, idim) = piglet_therm%smalls(i, (ibead - 1)*ndim + idim)
441 END DO
442 END DO
443
444 !fill temp2 with gaussian random noise
445 DO j = 1, nsp1
446 DO idim = 1, ndim
447 piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
448 !piglet_therm%temp2(j,idim) = 1.0_dp
449 END DO
450 END DO
451
452 i = (ibead - 1)*piglet_therm%ndim + 1
453 !smalls(:,i) = 1*S*temp2 + 0 * smalls
454 CALL dgemm("N", & ! S-matrix should not be transposed
455 "N", & ! tmp2 matrix shoud not be transposed
456 nsp1, & ! Number of rows of S-Matrix
457 ndim, & ! Number of columns of temp2 vector
458 nsp1, & ! Number of Columns of S-Matrix
459 1.0_dp, & ! Scaling of S-Matrix
460 piglet_therm%gle_s(:, :, ibead), & ! S-matrix
461 nsp1, & ! Leading dimension of S-matrix
462 piglet_therm%temp2, & ! temp2 vector
463 nsp1, & ! Leading dimension of temp2
464 0.0_dp, & ! scaling factor of added smalls vector
465 piglet_therm%smalls(:, i), & ! result vector
466 nsp1) ! Leading dimension of result vector
467
468 ! Now add the product of T-matrix * old smalls vectors
469 ! smalls (:,i) = 1*T*temp1 + 1*smalls
470 CALL dgemm("N", & ! T-matrix should not be transposed
471 "N", & ! temp1 matrix shoud not be transposed
472 nsp1, & ! Number of rows of T-Matrix
473 ndim, & ! Number of columns of temp1 vector
474 nsp1, & ! Number of Columns of T-Matrix
475 1.0_dp, & ! Scaling of T-Matrix
476 piglet_therm%gle_t(:, :, ibead), & ! T-matrix
477 nsp1, & ! Leading dimension of T-matrix
478 piglet_therm%temp1, & ! temp1 vector
479 nsp1, & ! Leading dimension of temp1
480 1.0_dp, & ! scaling factor of added smalls vector
481 piglet_therm%smalls(:, i), & ! result vector
482 nsp1) ! Leading dimension of result vector
483 END DO
484
485 ! Copy the mass scales momenta to the outgoing velocities
486 delta_ekin = 0.0_dp
487 DO idim = 1, ndim
488 DO ibead = 1, p
489 vnew(ibead, idim) = piglet_therm%smalls(1, (ibead - 1)*ndim + idim)/piglet_therm%sqrtmass(ibead, idim)
490 delta_ekin = delta_ekin + masses(ibead, idim)*( &
491 vnew(ibead, idim)*vnew(ibead, idim) - &
492 vold(ibead, idim)*vold(ibead, idim))
493 END DO
494 END DO
495
496 ! the piglet is such a strong thermostat, that it messes up the "exact" integration. The thermostats energy will rise lineary, because "it will suck up its own mess" (quote from Michele Ceriotti)
497 piglet_therm%thermostat_energy = piglet_therm%thermostat_energy - 0.5_dp*delta_ekin
498
499 CALL timestop(handle)
500
501 END SUBROUTINE pint_piglet_step
502
503! ***************************************************************************
504!> \brief releases the piglet environment
505!> \param piglet_therm piglet data to be released
506!> \author Felix Uhl
507! **************************************************************************************************
508 SUBROUTINE pint_piglet_release(piglet_therm)
509
510 TYPE(piglet_therm_type), INTENT(INOUT) :: piglet_therm
511
512 DEALLOCATE (piglet_therm%a_mat)
513 DEALLOCATE (piglet_therm%c_mat)
514 DEALLOCATE (piglet_therm%gle_t)
515 DEALLOCATE (piglet_therm%gle_s)
516 DEALLOCATE (piglet_therm%smalls)
517 DEALLOCATE (piglet_therm%temp1)
518 DEALLOCATE (piglet_therm%temp2)
519 DEALLOCATE (piglet_therm%sqrtmass)
520
521 END SUBROUTINE pint_piglet_release
522
523! ***************************************************************************
524!> \brief adjust the unit of A MAT for piglet
525!> \param a_mat ...
526!> \param p ...
527!> \param nsp1 ...
528!> \param myunit ...
529!> \param msg ...
530!> \author Felix Uhl
531! **************************************************************************************************
532 SUBROUTINE a_mat_to_cp2k(a_mat, p, nsp1, myunit, msg)
533
534 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: a_mat
535 INTEGER, INTENT(IN) :: p, nsp1
536 CHARACTER(LEN=20), INTENT(IN) :: myunit
537 CHARACTER(default_string_length), INTENT(OUT) :: msg
538
539 CHARACTER(LEN=20) :: isunit
540 INTEGER :: i, imode, j
541
542 msg = ""
543 SELECT CASE (trim(myunit))
544 CASE ("femtoseconds^-1")
545 isunit = "fs^-1"
546 CASE ("picoseconds^-1")
547 isunit = "ps^-1"
548 CASE ("seconds^-1")
549 isunit = "s^-1"
550 CASE ("atomic time units^-1")
551 RETURN
552 CASE DEFAULT
553 msg = "Unknown unit of A-Matrices for PIGLET. Assuming a.u."
554 cpwarn(msg)
555 RETURN
556 END SELECT
557
558 DO imode = 1, p
559 DO j = 1, nsp1
560 DO i = 1, nsp1
561 a_mat(i, j, imode) = cp_unit_to_cp2k(a_mat(i, j, imode), trim(isunit))
562 END DO
563 END DO
564 END DO
565
566 END SUBROUTINE
567! ***************************************************************************
568!> \brief adjust the unit of C MAT for piglet
569!> \param c_mat ...
570!> \param p ...
571!> \param nsp1 ...
572!> \param myunit ...
573!> \param msg ...
574!> \author Felix Uhl
575! **************************************************************************************************
576 SUBROUTINE c_mat_to_cp2k(c_mat, p, nsp1, myunit, msg)
577
578 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: c_mat
579 INTEGER, INTENT(IN) :: p, nsp1
580 CHARACTER(LEN=20), INTENT(IN) :: myunit
581 CHARACTER(default_string_length), INTENT(OUT) :: msg
582
583 CHARACTER(LEN=20) :: isunit
584 INTEGER :: i, imode, j
585
586 msg = ""
587 SELECT CASE (trim(myunit))
588 CASE ("eV")
589 isunit = "eV"
590 CASE ("K")
591 isunit = "K_e"
592 CASE ("atomic energy units ")
593 RETURN
594 CASE DEFAULT
595 msg = "Unknown unit of C-Matrices for PIGLET. Assuming a.u."
596 cpwarn(msg)
597 RETURN
598 END SELECT
599
600 DO imode = 1, p
601 DO j = 1, nsp1
602 DO i = 1, nsp1
603 c_mat(i, j, imode) = cp_unit_to_cp2k(c_mat(i, j, imode), trim(isunit))
604 END DO
605 END DO
606 END DO
607
608 END SUBROUTINE c_mat_to_cp2k
609
610! ***************************************************************************
611!> \brief checks if the matrices are suited for the target temperature
612!> \param line ...
613!> \param propagator ...
614!> \param targettemp ...
615!> \author Felix Uhl
616! **************************************************************************************************
617 SUBROUTINE check_temperature(line, propagator, targettemp)
618
619 CHARACTER(len=*), INTENT(IN) :: line
620 INTEGER, INTENT(IN) :: propagator
621 REAL(kind=dp), INTENT(IN) :: targettemp
622
623 CHARACTER(len=default_string_length) :: msg
624 INTEGER :: posnumber
625 REAL(kind=dp) :: convttemp, deviation, matrixtemp, ttemp
626
627 deviation = 100.0d0
628 posnumber = index(line, "T=") + 2
629 IF (propagator == propagator_rpmd) ttemp = targettemp
630 !Get the matrix temperature
631 READ (line(posnumber:), *) matrixtemp
632 msg = ""
633 IF (index(line, "K") /= 0) THEN
634 convttemp = cp_unit_from_cp2k(ttemp, "K")
635 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
636 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
637 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
638 cpwarn(msg)
639 END IF
640 ELSE IF (index(line, "eV") /= 0) THEN
641 convttemp = cp_unit_from_cp2k(ttemp, "K")/11604.505_dp
642 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
643 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
644 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
645 cpwarn(msg)
646 END IF
647 ELSE IF (index(line, "atomic energy units") /= 0) THEN
648 convttemp = ttemp
649 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
650 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
651 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
652 cpwarn(msg)
653 END IF
654 ELSE
655 WRITE (unit=msg, fmt=*) "Unknown PIGLET matrix temperature. Assuming a.u."
656 cpwarn(msg)
657 convttemp = ttemp
658 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
659 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
660 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
661 cpwarn(msg)
662 END IF
663 END IF
664
665 END SUBROUTINE check_temperature
666
667! ***************************************************************************
668!> \brief returns the piglet kinetic energy contribution
669!> \param pint_env ...
670!> \author Felix Uhl
671! **************************************************************************************************
672 ELEMENTAL SUBROUTINE pint_calc_piglet_energy(pint_env)
673 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
674
675 IF (ASSOCIATED(pint_env%piglet_therm)) THEN
676 pint_env%e_piglet = pint_env%piglet_therm%thermostat_energy
677 ELSE
678 pint_env%e_piglet = 0.0d0
679 END IF
680
681 END SUBROUTINE pint_calc_piglet_energy
682
683! ***************************************************************************
684!> \brief calculates S from S*TRANSPOSED(S) by diagonalizing
685!> if S*TRANSPOSED(S) is a positive definite matrix
686!> \param n order of input matrix
687!> \param SST matrix to be decomposed
688!> \param S result matrix
689!> \author Felix Uhl
690! **************************************************************************************************
691 SUBROUTINE sqrt_pos_def_mat(n, SST, S)
692
693 INTEGER, INTENT(IN) :: n
694 REAL(kind=dp), DIMENSION(n, n), INTENT(IN) :: sst
695 REAL(kind=dp), DIMENSION(n, n), INTENT(OUT) :: s
696
697 INTEGER :: i, info, lwork
698 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
699 REAL(kind=dp), DIMENSION(1) :: tmplwork
700 REAL(kind=dp), DIMENSION(n) :: eigval
701 REAL(kind=dp), DIMENSION(n, n) :: a, tmpmatrix
702
703! order of input matrix
704! matrix to be decomposed
705! result matrix
706! Variables for eigenvalue/vector computation
707! store matrix here to pass to lapack routine DSYEVR
708! Array to contain the eigenvalues
709! size of temporary real work array
710! temporary real work array
711! information about success
712! counter
713
714 eigval(:) = 0.0_dp
715 a(:, :) = 0.0_dp
716 a(:, :) = sst(:, :)
717
718 !first call to figure out how big work array needs to be
719 CALL dsyev('V', & ! Compute eigenvalues and eigenvectors
720 'U', & ! Store upper triagonal matrix
721 n, & ! order of matrix A to calculate the eigenvalues/vectors from
722 a, & ! Matrix to calculate the eigenvalues/vectors from
723 n, & ! leading order of matrix A
724 eigval, & ! Array to contain the eigenvalues
725 tmplwork, & ! temporary real work array
726 -1, & ! size of temporary real work array
727 info) ! information about success
728
729 lwork = int(tmplwork(1) + 0.5_dp)
730 ALLOCATE (work(lwork))
731 work(:) = 0.0_dp
732
733 CALL dsyev('V', & ! Compute eigenvalues and eigenvectors
734 'U', & ! Store upper triagonal matrix
735 n, & ! order of matrix A to calculate the eigenvalues/vectors from
736 a, & ! Matrix to calculate the eigenvalues/vectors from
737 n, & ! leading order of matrix A
738 eigval, & ! Array to contain the eigenvalues
739 work, & ! temporary real work array
740 lwork, & ! size of temporary real work array
741 info) ! information about success
742 DEALLOCATE (work)
743 ! A-matrix now contains the eigenvectors
744
745 s(:, :) = 0.0_dp
746 DO i = 1, n
747 ! In case that numerics made some eigenvalues negative
748 IF (eigval(i) > 0.0_dp) THEN
749 s(i, i) = sqrt(eigval(i))
750 END IF
751 END DO
752
753 tmpmatrix(:, :) = 0.0_dp
754 ! Transform matrix back
755 !tmpmatrix = A*S
756 CALL dgemm('N', & ! A-matrix is not transposed
757 'N', & ! S-matrix is not transposed
758 n, & ! number of rows of A-matrix
759 n, & ! number of columns of S-matrix
760 n, & ! number of columns of A-matrix and number of rows of S-matrix
761 1.0d0, & ! scaling factor of A-matrix
762 a, & ! A-matrix
763 n, & ! leading dimension of A-matrix as declared
764 s, & ! S-matrix
765 n, & ! leading dimension of S-matrix as declared
766 0.0d0, & ! scaling of tmpmatrix as additive
767 tmpmatrix, & ! result matrix: tmpmatrix
768 n) ! leading dimension of tmpmatrix
769 !S = tmpmatrix*TRANSPOSE(A) = A*S*TRANSPOSE(A)
770 CALL dgemm('N', & ! tmpmatrix not transposed
771 'T', & ! A-matrix is transposed
772 n, & ! number of rows of tmpmatrix
773 n, & ! number of columns of A-matrix
774 n, & ! number of columns of tmpmatrix and rows of A-matrix
775 1.0d0, & ! scaling factor of tmpmatrix
776 tmpmatrix, & ! tmpmatrix
777 n, & ! leading dimension of tmpmatrix as declared
778 a, & ! A-matrix
779 n, & ! leading dimension of A-matrix as declared
780 0.0d0, & ! scaling of S-matrix as additive
781 s, & ! result matrix: S-matrix
782 n) ! leading dimension of S-matrix as declared
783
784 END SUBROUTINE sqrt_pos_def_mat
785
786END MODULE pint_piglet
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:195
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1179
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1150
gle_system_dynamics
Definition gle_system_dynamics.F:12
gle_system_dynamics::gle_matrix_exp
subroutine, public gle_matrix_exp(m, n, j, k, em)
...
Definition gle_system_dynamics.F:231
gle_system_dynamics::gle_cholesky_stab
subroutine, public gle_cholesky_stab(sst, s, n)
...
Definition gle_system_dynamics.F:278
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::propagator_rpmd
integer, parameter, public propagator_rpmd
Definition input_constants.F:305
input_constants::matrix_init_cholesky
integer, parameter, public matrix_init_cholesky
Definition input_constants.F:313
input_constants::matrix_init_diagonal
integer, parameter, public matrix_init_diagonal
Definition input_constants.F:313
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_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
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::rng_stream_type_from_record
type(rng_stream_type) function, public rng_stream_type_from_record(rng_record)
Create a RNG stream from a record given as an internal file (string).
Definition parallel_rng_types.F:564
parallel_rng_types::rng_record_length
integer, parameter, public rng_record_length
Definition parallel_rng_types.F:65
parallel_rng_types::gaussian
integer, parameter, public gaussian
Definition parallel_rng_types.F:73
pint_io
I/O subroutines for pint_env.
Definition pint_io.F:13
pint_io::pint_write_line
subroutine, public pint_write_line(line)
Writes out a line of text to the default output unit.
Definition pint_io.F:76
pint_piglet
Methods to apply the piglet thermostat to PI runs.
Definition pint_piglet.F:14
pint_piglet::pint_calc_piglet_energy
elemental subroutine, public pint_calc_piglet_energy(pint_env)
returns the piglet kinetic energy contribution
Definition pint_piglet.F:673
pint_piglet::pint_piglet_release
subroutine, public pint_piglet_release(piglet_therm)
releases the piglet environment
Definition pint_piglet.F:509
pint_piglet::pint_piglet_create
subroutine, public pint_piglet_create(piglet_therm, pint_env, section)
creates the data structure for a piglet thermostating in PI runs
Definition pint_piglet.F:92
pint_piglet::pint_piglet_step
subroutine, public pint_piglet_step(vold, vnew, first_mode, masses, piglet_therm)
...
Definition pint_piglet.F:411
pint_piglet::pint_piglet_init
subroutine, public pint_piglet_init(piglet_therm, pint_env, section, dt, para_env)
initializes the data for a piglet run
Definition pint_piglet.F:142
pint_types
Definition pint_types.F:8
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
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
pint_types::piglet_therm_type
data to use the piglet thermostat
Definition pint_types.F:232
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112