(git:71c3ab0)
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-2026 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.
497 ! The thermostats energy will rise lineary, because "it will suck up its own mess" (quote from Michele Ceriotti)
498 piglet_therm%thermostat_energy = piglet_therm%thermostat_energy - 0.5_dp*delta_ekin
499
500 CALL timestop(handle)
501
502 END SUBROUTINE pint_piglet_step
503
504! ***************************************************************************
505!> \brief releases the piglet environment
506!> \param piglet_therm piglet data to be released
507!> \author Felix Uhl
508! **************************************************************************************************
509 SUBROUTINE pint_piglet_release(piglet_therm)
510
511 TYPE(piglet_therm_type), INTENT(INOUT) :: piglet_therm
512
513 DEALLOCATE (piglet_therm%a_mat)
514 DEALLOCATE (piglet_therm%c_mat)
515 DEALLOCATE (piglet_therm%gle_t)
516 DEALLOCATE (piglet_therm%gle_s)
517 DEALLOCATE (piglet_therm%smalls)
518 DEALLOCATE (piglet_therm%temp1)
519 DEALLOCATE (piglet_therm%temp2)
520 DEALLOCATE (piglet_therm%sqrtmass)
521
522 END SUBROUTINE pint_piglet_release
523
524! ***************************************************************************
525!> \brief adjust the unit of A MAT for piglet
526!> \param a_mat ...
527!> \param p ...
528!> \param nsp1 ...
529!> \param myunit ...
530!> \param msg ...
531!> \author Felix Uhl
532! **************************************************************************************************
533 SUBROUTINE a_mat_to_cp2k(a_mat, p, nsp1, myunit, msg)
534
535 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: a_mat
536 INTEGER, INTENT(IN) :: p, nsp1
537 CHARACTER(LEN=20), INTENT(IN) :: myunit
538 CHARACTER(default_string_length), INTENT(OUT) :: msg
539
540 CHARACTER(LEN=20) :: isunit
541 INTEGER :: i, imode, j
542
543 msg = ""
544 SELECT CASE (trim(myunit))
545 CASE ("femtoseconds^-1")
546 isunit = "fs^-1"
547 CASE ("picoseconds^-1")
548 isunit = "ps^-1"
549 CASE ("seconds^-1")
550 isunit = "s^-1"
551 CASE ("atomic time units^-1")
552 RETURN
553 CASE DEFAULT
554 msg = "Unknown unit of A-Matrices for PIGLET. Assuming a.u."
555 cpwarn(msg)
556 RETURN
557 END SELECT
558
559 DO imode = 1, p
560 DO j = 1, nsp1
561 DO i = 1, nsp1
562 a_mat(i, j, imode) = cp_unit_to_cp2k(a_mat(i, j, imode), trim(isunit))
563 END DO
564 END DO
565 END DO
566
567 END SUBROUTINE a_mat_to_cp2k
568! ***************************************************************************
569!> \brief adjust the unit of C MAT for piglet
570!> \param c_mat ...
571!> \param p ...
572!> \param nsp1 ...
573!> \param myunit ...
574!> \param msg ...
575!> \author Felix Uhl
576! **************************************************************************************************
577 SUBROUTINE c_mat_to_cp2k(c_mat, p, nsp1, myunit, msg)
578
579 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: c_mat
580 INTEGER, INTENT(IN) :: p, nsp1
581 CHARACTER(LEN=20), INTENT(IN) :: myunit
582 CHARACTER(default_string_length), INTENT(OUT) :: msg
583
584 CHARACTER(LEN=20) :: isunit
585 INTEGER :: i, imode, j
586
587 msg = ""
588 SELECT CASE (trim(myunit))
589 CASE ("eV")
590 isunit = "eV"
591 CASE ("K")
592 isunit = "K_e"
593 CASE ("atomic energy units ")
594 RETURN
595 CASE DEFAULT
596 msg = "Unknown unit of C-Matrices for PIGLET. Assuming a.u."
597 cpwarn(msg)
598 RETURN
599 END SELECT
600
601 DO imode = 1, p
602 DO j = 1, nsp1
603 DO i = 1, nsp1
604 c_mat(i, j, imode) = cp_unit_to_cp2k(c_mat(i, j, imode), trim(isunit))
605 END DO
606 END DO
607 END DO
608
609 END SUBROUTINE c_mat_to_cp2k
610
611! ***************************************************************************
612!> \brief checks if the matrices are suited for the target temperature
613!> \param line ...
614!> \param propagator ...
615!> \param targettemp ...
616!> \author Felix Uhl
617! **************************************************************************************************
618 SUBROUTINE check_temperature(line, propagator, targettemp)
619
620 CHARACTER(len=*), INTENT(IN) :: line
621 INTEGER, INTENT(IN) :: propagator
622 REAL(kind=dp), INTENT(IN) :: targettemp
623
624 CHARACTER(len=default_string_length) :: msg
625 INTEGER :: posnumber
626 REAL(kind=dp) :: convttemp, deviation, matrixtemp, ttemp
627
628 deviation = 100.0d0
629 posnumber = index(line, "T=") + 2
630 IF (propagator == propagator_rpmd) ttemp = targettemp
631 !Get the matrix temperature
632 READ (line(posnumber:), *) matrixtemp
633 msg = ""
634 IF (index(line, "K") /= 0) THEN
635 convttemp = cp_unit_from_cp2k(ttemp, "K")
636 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
637 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
638 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
639 cpwarn(msg)
640 END IF
641 ELSE IF (index(line, "eV") /= 0) THEN
642 convttemp = cp_unit_from_cp2k(ttemp, "K")/11604.505_dp
643 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
644 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
645 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
646 cpwarn(msg)
647 END IF
648 ELSE IF (index(line, "atomic energy units") /= 0) THEN
649 convttemp = ttemp
650 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
651 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
652 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
653 cpwarn(msg)
654 END IF
655 ELSE
656 WRITE (unit=msg, fmt=*) "Unknown PIGLET matrix temperature. Assuming a.u."
657 cpwarn(msg)
658 convttemp = ttemp
659 IF (abs(convttemp - matrixtemp) > convttemp/deviation) THEN
660 WRITE (unit=msg, fmt=*) "PIGLET Simulation temperature (", &
661 convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
662 cpwarn(msg)
663 END IF
664 END IF
665
666 END SUBROUTINE check_temperature
667
668! ***************************************************************************
669!> \brief returns the piglet kinetic energy contribution
670!> \param pint_env ...
671!> \author Felix Uhl
672! **************************************************************************************************
673 ELEMENTAL SUBROUTINE pint_calc_piglet_energy(pint_env)
674 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
675
676 IF (ASSOCIATED(pint_env%piglet_therm)) THEN
677 pint_env%e_piglet = pint_env%piglet_therm%thermostat_energy
678 ELSE
679 pint_env%e_piglet = 0.0d0
680 END IF
681
682 END SUBROUTINE pint_calc_piglet_energy
683
684! ***************************************************************************
685!> \brief calculates S from S*TRANSPOSED(S) by diagonalizing
686!> if S*TRANSPOSED(S) is a positive definite matrix
687!> \param n order of input matrix
688!> \param SST matrix to be decomposed
689!> \param S result matrix
690!> \author Felix Uhl
691! **************************************************************************************************
692 SUBROUTINE sqrt_pos_def_mat(n, SST, S)
693
694 INTEGER, INTENT(IN) :: n
695 REAL(kind=dp), DIMENSION(n, n), INTENT(IN) :: sst
696 REAL(kind=dp), DIMENSION(n, n), INTENT(OUT) :: s
697
698 INTEGER :: i, info, lwork
699 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
700 REAL(kind=dp), DIMENSION(1) :: tmplwork
701 REAL(kind=dp), DIMENSION(n) :: eigval
702 REAL(kind=dp), DIMENSION(n, n) :: a, tmpmatrix
703
704! order of input matrix
705! matrix to be decomposed
706! result matrix
707! Variables for eigenvalue/vector computation
708! store matrix here to pass to lapack routine DSYEVR
709! Array to contain the eigenvalues
710! size of temporary real work array
711! temporary real work array
712! information about success
713! counter
714
715 eigval(:) = 0.0_dp
716 a(:, :) = 0.0_dp
717 a(:, :) = sst(:, :)
718
719 !first call to figure out how big work array needs to be
720 CALL dsyev('V', & ! Compute eigenvalues and eigenvectors
721 'U', & ! Store upper triagonal matrix
722 n, & ! order of matrix A to calculate the eigenvalues/vectors from
723 a, & ! Matrix to calculate the eigenvalues/vectors from
724 n, & ! leading order of matrix A
725 eigval, & ! Array to contain the eigenvalues
726 tmplwork, & ! temporary real work array
727 -1, & ! size of temporary real work array
728 info) ! information about success
729
730 lwork = int(tmplwork(1) + 0.5_dp)
731 ALLOCATE (work(lwork))
732 work(:) = 0.0_dp
733
734 CALL dsyev('V', & ! Compute eigenvalues and eigenvectors
735 'U', & ! Store upper triagonal matrix
736 n, & ! order of matrix A to calculate the eigenvalues/vectors from
737 a, & ! Matrix to calculate the eigenvalues/vectors from
738 n, & ! leading order of matrix A
739 eigval, & ! Array to contain the eigenvalues
740 work, & ! temporary real work array
741 lwork, & ! size of temporary real work array
742 info) ! information about success
743 DEALLOCATE (work)
744 ! A-matrix now contains the eigenvectors
745
746 s(:, :) = 0.0_dp
747 DO i = 1, n
748 ! In case that numerics made some eigenvalues negative
749 IF (eigval(i) > 0.0_dp) THEN
750 s(i, i) = sqrt(eigval(i))
751 END IF
752 END DO
753
754 tmpmatrix(:, :) = 0.0_dp
755 ! Transform matrix back
756 !tmpmatrix = A*S
757 CALL dgemm('N', & ! A-matrix is not transposed
758 'N', & ! S-matrix is not transposed
759 n, & ! number of rows of A-matrix
760 n, & ! number of columns of S-matrix
761 n, & ! number of columns of A-matrix and number of rows of S-matrix
762 1.0d0, & ! scaling factor of A-matrix
763 a, & ! A-matrix
764 n, & ! leading dimension of A-matrix as declared
765 s, & ! S-matrix
766 n, & ! leading dimension of S-matrix as declared
767 0.0d0, & ! scaling of tmpmatrix as additive
768 tmpmatrix, & ! result matrix: tmpmatrix
769 n) ! leading dimension of tmpmatrix
770 !S = tmpmatrix*TRANSPOSE(A) = A*S*TRANSPOSE(A)
771 CALL dgemm('N', & ! tmpmatrix not transposed
772 'T', & ! A-matrix is transposed
773 n, & ! number of rows of tmpmatrix
774 n, & ! number of columns of A-matrix
775 n, & ! number of columns of tmpmatrix and rows of A-matrix
776 1.0d0, & ! scaling factor of tmpmatrix
777 tmpmatrix, & ! tmpmatrix
778 n, & ! leading dimension of tmpmatrix as declared
779 a, & ! A-matrix
780 n, & ! leading dimension of A-matrix as declared
781 0.0d0, & ! scaling of S-matrix as additive
782 s, & ! result matrix: S-matrix
783 n) ! leading dimension of S-matrix as declared
784
785 END SUBROUTINE sqrt_pos_def_mat
786
787END 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:225
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:311
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:122
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:1251
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:1222
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:18
input_constants::propagator_rpmd
integer, parameter, public propagator_rpmd
Definition input_constants.F:315
input_constants::matrix_init_cholesky
integer, parameter, public matrix_init_cholesky
Definition input_constants.F:324
input_constants::matrix_init_diagonal
integer, parameter, public matrix_init_diagonal
Definition input_constants.F:324
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: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
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:565
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:674
pint_piglet::pint_piglet_release
subroutine, public pint_piglet_release(piglet_therm)
releases the piglet environment
Definition pint_piglet.F:510
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:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
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