(git:b279b6b)
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 ! **************************************************************************************************
14 MODULE 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 
82 CONTAINS
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 
786 END 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_cholesky_stab
subroutine, public gle_cholesky_stab(SST, S, n)
...
Definition: gle_system_dynamics.F:278
gle_system_dynamics::gle_matrix_exp
subroutine, public gle_matrix_exp(M, n, j, k, EM)
...
Definition: gle_system_dynamics.F:231
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