Line data Source code
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 : ! **************************************************************************************************
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 54 : 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 2 : 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 2 : piglet_therm%nsp1 = piglet_therm%nsp1 + 1
101 2 : p = pint_env%p
102 2 : piglet_therm%p = pint_env%p
103 2 : ndim = pint_env%ndim
104 2 : piglet_therm%ndim = pint_env%ndim
105 2 : NULLIFY (piglet_therm%a_mat)
106 10 : ALLOCATE (piglet_therm%a_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
107 2186 : piglet_therm%a_mat(:, :, :) = 0.0_dp
108 2 : NULLIFY (piglet_therm%c_mat)
109 8 : ALLOCATE (piglet_therm%c_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
110 2186 : piglet_therm%c_mat(:, :, :) = 0.0_dp
111 2 : NULLIFY (piglet_therm%gle_t)
112 8 : ALLOCATE (piglet_therm%gle_t(piglet_therm%nsp1, piglet_therm%nsp1, p))
113 2186 : piglet_therm%gle_t(:, :, :) = 0.0_dp
114 2 : NULLIFY (piglet_therm%gle_s)
115 8 : ALLOCATE (piglet_therm%gle_s(piglet_therm%nsp1, piglet_therm%nsp1, p))
116 2186 : piglet_therm%gle_s(:, :, :) = 0.0_dp
117 2 : NULLIFY (piglet_therm%smalls)
118 8 : ALLOCATE (piglet_therm%smalls(piglet_therm%nsp1, ndim*p))
119 1105922 : piglet_therm%smalls(:, :) = 0.0_dp
120 2 : NULLIFY (piglet_therm%temp1)
121 8 : ALLOCATE (piglet_therm%temp1(piglet_therm%nsp1, ndim))
122 92162 : piglet_therm%temp1(:, :) = 0.0_dp
123 2 : NULLIFY (piglet_therm%temp2)
124 6 : ALLOCATE (piglet_therm%temp2(piglet_therm%nsp1, ndim))
125 92162 : piglet_therm%temp2(:, :) = 0.0_dp
126 2 : NULLIFY (piglet_therm%sqrtmass)
127 8 : ALLOCATE (piglet_therm%sqrtmass(piglet_therm%p, piglet_therm%ndim))
128 119810 : piglet_therm%sqrtmass(:, :) = 0.0_dp
129 :
130 2 : 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 2 : 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 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: smallstmp
158 : REAL(KIND=dp), &
159 4 : DIMENSION(piglet_therm%nsp1, piglet_therm%nsp1) :: Mtmp, tmpmatrix
160 : TYPE(section_vals_type), POINTER :: rng_section, smalls_section
161 :
162 2 : p = piglet_therm%p
163 2 : pint_env%e_piglet = 0.0_dp
164 2 : pint_env%piglet_therm%thermostat_energy = 0.0_dp
165 2 : CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=piglet_therm%thermostat_energy)
166 2 : CALL section_vals_val_get(section, "SMATRIX_INIT", i_val=matrix_init)
167 : !Read the matices
168 :
169 2 : IF (pint_env%propagator%prop_kind /= propagator_rpmd) THEN
170 0 : CPABORT("PIGLET is designed to work with the RPMD propagator")
171 : END IF
172 : ! Select algorithm for S-matrix initialization
173 2 : IF (matrix_init == matrix_init_cholesky) THEN
174 0 : IF (para_env%is_source()) THEN
175 0 : CALL pint_write_line("PIGLET| Initalizing S-matrices using cholesky decomposition.")
176 : END IF
177 2 : ELSE IF (matrix_init == matrix_init_diagonal) THEN
178 2 : IF (para_env%is_source()) THEN
179 1 : CALL pint_write_line("PIGLET| Initalizing S-matrices using full diagonalization.")
180 : END IF
181 : ELSE
182 0 : CPWARN("No PIGLET init algorithm selected. Selecting cholesky decomposition")
183 0 : matrix_init = matrix_init_cholesky
184 : END IF
185 :
186 2 : IF (para_env%is_source()) THEN
187 : ! Read input matrices
188 1 : WRITE (default_format, '(A,I10,A)') "(A", default_string_length, ")"
189 : CALL section_vals_val_get(section, "MATRICES_FILE_NAME", &
190 1 : c_val=matrices_file_name)
191 1 : CALL pint_write_line("PIGLET| Reading PIGLET matrices from file: ")
192 1 : CALL pint_write_line("PIGLET| "//TRIM(matrices_file_name))
193 : CALL open_file(file_name=TRIM(matrices_file_name), &
194 1 : file_action="READ", file_status="OLD", unit_number=input_unit)
195 1 : read_err = 0
196 1 : msg = ""
197 20 : DO WHILE (read_err == 0)
198 20 : READ (input_unit, default_format, iostat=read_err) temp_input
199 20 : IF (read_err /= 0) EXIT
200 : !Parse comment section
201 20 : IF (INDEX(temp_input, "#") /= 0) THEN
202 16 : IF (INDEX(temp_input, "T=") /= 0) THEN
203 : CALL check_temperature(line=temp_input, &
204 : propagator=pint_env%propagator%prop_kind, &
205 1 : targettemp=pint_env%kT*pint_env%propagator%temp_sim2phys)
206 15 : ELSE IF (INDEX(temp_input, "A MATRIX") /= 0) THEN
207 1 : obrac = INDEX(temp_input, "(") + 1
208 1 : cbrac = INDEX(temp_input, ")") - 1
209 1 : read_unit = temp_input(obrac:cbrac)
210 13 : DO imode = 1, p
211 12 : READ (input_unit, default_format) temp_input
212 121 : DO i = 1, piglet_therm%nsp1
213 : READ (input_unit, *, iostat=read_err) &
214 1080 : (piglet_therm%a_mat(i, j, imode), j=1, piglet_therm%nsp1)
215 120 : IF (read_err /= 0) THEN
216 0 : WRITE (UNIT=msg, FMT=*) "Invalid PIGLET A-matrix Nr.", i - 1
217 0 : CPABORT(msg)
218 0 : EXIT
219 : END IF
220 : END DO
221 : END DO
222 : !convert to cp2k units
223 1 : IF (read_err == 0) THEN
224 : CALL a_mat_to_cp2k(piglet_therm%a_mat, p, &
225 1 : piglet_therm%nsp1, read_unit, msg)
226 : END IF
227 14 : ELSE IF (INDEX(temp_input, "C MATRIX") /= 0) THEN
228 1 : obrac = INDEX(temp_input, "(") + 1
229 1 : cbrac = INDEX(temp_input, ")") - 1
230 1 : read_unit = temp_input(obrac:cbrac)
231 13 : DO imode = 1, p
232 12 : READ (input_unit, default_format) temp_input
233 121 : DO i = 1, piglet_therm%nsp1
234 : READ (input_unit, *, iostat=read_err) &
235 1080 : (piglet_therm%c_mat(i, j, imode), j=1, piglet_therm%nsp1)
236 120 : IF (read_err /= 0) THEN
237 0 : WRITE (UNIT=msg, FMT=*) "Invalid PIGLET C-matrix Nr.", i - 1
238 0 : CPABORT(msg)
239 0 : EXIT
240 : END IF
241 : END DO
242 : END DO
243 1 : IF (read_err == 0) THEN
244 : CALL c_mat_to_cp2k(piglet_therm%c_mat, p, &
245 1 : piglet_therm%nsp1, read_unit, msg)
246 : END IF
247 : END IF
248 : END IF
249 : END DO
250 1 : 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 4370 : para_env%source)
255 : CALL para_env%bcast(piglet_therm%c_mat, &
256 4370 : para_env%source)
257 :
258 : !prepare Random number generator
259 2 : NULLIFY (rng_section)
260 : rng_section => section_vals_get_subs_vals(section, &
261 2 : subsection_name="RNG_INIT")
262 2 : CALL section_vals_get(rng_section, explicit=explicit)
263 2 : IF (explicit) THEN
264 : CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
265 0 : i_rep_val=1, c_val=rng_record)
266 0 : piglet_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
267 : ELSE
268 18 : 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 2 : seed=initial_seed)
273 : END IF
274 :
275 : !Compute the T and S matrices on every mpi process
276 26 : 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 2184 : 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 24 : 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 24 : piglet_therm%nsp1) ! leading dimension of Mtmp
314 : ! C - T*C*TRANSPOSE(T):
315 2184 : Mtmp(:, :) = piglet_therm%c_mat(:, :, i) - Mtmp(:, :)
316 :
317 26 : 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 0 : piglet_therm%nsp1) ! Size of the matrix
322 24 : 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 24 : 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 1105922 : piglet_therm%smalls = 0.0_dp
335 2 : NULLIFY (smalls_section)
336 2 : smalls_section => section_vals_get_subs_vals(section, subsection_name="EXTRA_DOF")
337 2 : CALL section_vals_get(smalls_section, explicit=explicit)
338 2 : IF (explicit) THEN
339 0 : NULLIFY (smallstmp)
340 : CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
341 0 : n_rep_val=ns)
342 : CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
343 0 : r_vals=smallstmp)
344 0 : i = 1
345 0 : DO isp = 2, piglet_therm%nsp1
346 0 : DO ibead = 1, piglet_therm%p*piglet_therm%ndim
347 0 : piglet_therm%smalls(isp, ibead) = smallstmp(i)
348 0 : i = i + 1
349 : END DO
350 : END DO
351 : ELSE
352 26 : DO ibead = 1, piglet_therm%p
353 24 : IF (matrix_init == matrix_init_cholesky) THEN
354 : CALL gle_cholesky_stab(piglet_therm%c_mat(:, :, ibead), & ! Matrix to decompose
355 : Mtmp, & ! Result
356 0 : piglet_therm%nsp1) ! Size of Matrix
357 24 : 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 24 : Mtmp) ! result
362 : END IF
363 : ! Fill a vector with random numbers
364 110616 : DO idim = 1, piglet_therm%ndim
365 1105944 : DO j = 1, piglet_therm%nsp1
366 1105920 : 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 24 : piglet_therm%nsp1) ! leading dimension of result matrix
383 :
384 110618 : DO idim = 1, piglet_therm%ndim
385 110592 : j = (idim - 1)*piglet_therm%p + ibead
386 1105944 : DO i = 1, piglet_therm%nsp1
387 1105920 : 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 9218 : DO idim = 1, piglet_therm%ndim
395 119810 : DO ibead = 1, piglet_therm%p
396 119808 : piglet_therm%sqrtmass(ibead, idim) = SQRT(pint_env%mass_fict(ibead, idim))
397 : END DO
398 : END DO
399 :
400 2 : 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 20 : 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 20 : CALL timeset(routineN, handle)
423 20 : nsp1 = piglet_therm%nsp1
424 20 : ndim = piglet_therm%ndim
425 20 : p = piglet_therm%p
426 20 : 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 260 : 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 1106160 : DO idim = 1, ndim
435 1106160 : 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 1106160 : DO idim = 1, ndim
439 9953520 : DO i = 2, nsp1
440 9953280 : 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 2400 : DO j = 1, nsp1
446 9955680 : DO idim = 1, ndim
447 9955440 : 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 240 : 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 240 : 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 260 : nsp1) ! Leading dimension of result vector
483 : END DO
484 :
485 : ! Copy the mass scales momenta to the outgoing velocities
486 20 : delta_ekin = 0.0_dp
487 92180 : DO idim = 1, ndim
488 1198100 : DO ibead = 1, p
489 1105920 : 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 1198080 : 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 20 : piglet_therm%thermostat_energy = piglet_therm%thermostat_energy - 0.5_dp*delta_ekin
499 :
500 20 : CALL timestop(handle)
501 :
502 20 : 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 2 : SUBROUTINE pint_piglet_release(piglet_therm)
510 :
511 : TYPE(piglet_therm_type), INTENT(INOUT) :: piglet_therm
512 :
513 2 : DEALLOCATE (piglet_therm%a_mat)
514 2 : DEALLOCATE (piglet_therm%c_mat)
515 2 : DEALLOCATE (piglet_therm%gle_t)
516 2 : DEALLOCATE (piglet_therm%gle_s)
517 2 : DEALLOCATE (piglet_therm%smalls)
518 2 : DEALLOCATE (piglet_therm%temp1)
519 2 : DEALLOCATE (piglet_therm%temp2)
520 2 : DEALLOCATE (piglet_therm%sqrtmass)
521 :
522 2 : 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 1 : 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 1 : msg = ""
544 1 : SELECT CASE (TRIM(myunit))
545 : CASE ("femtoseconds^-1")
546 0 : isunit = "fs^-1"
547 : CASE ("picoseconds^-1")
548 1 : isunit = "ps^-1"
549 : CASE ("seconds^-1")
550 0 : isunit = "s^-1"
551 : CASE ("atomic time units^-1")
552 0 : RETURN
553 : CASE DEFAULT
554 0 : msg = "Unknown unit of A-Matrices for PIGLET. Assuming a.u."
555 0 : CPWARN(msg)
556 1 : RETURN
557 : END SELECT
558 :
559 13 : DO imode = 1, p
560 121 : DO j = 1, nsp1
561 1092 : DO i = 1, nsp1
562 1080 : 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 1 : 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 1 : msg = ""
588 1 : SELECT CASE (TRIM(myunit))
589 : CASE ("eV")
590 0 : isunit = "eV"
591 : CASE ("K")
592 1 : isunit = "K_e"
593 : CASE ("atomic energy units ")
594 0 : RETURN
595 : CASE DEFAULT
596 0 : msg = "Unknown unit of C-Matrices for PIGLET. Assuming a.u."
597 0 : CPWARN(msg)
598 1 : RETURN
599 : END SELECT
600 :
601 13 : DO imode = 1, p
602 121 : DO j = 1, nsp1
603 1092 : DO i = 1, nsp1
604 1080 : 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 1 : 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 1 : deviation = 100.0d0
629 1 : posnumber = INDEX(line, "T=") + 2
630 1 : IF (propagator == propagator_rpmd) ttemp = targettemp
631 : !Get the matrix temperature
632 1 : READ (line(posnumber:), *) matrixtemp
633 1 : msg = ""
634 1 : IF (INDEX(line, "K") /= 0) THEN
635 1 : convttemp = cp_unit_from_cp2k(ttemp, "K")
636 1 : IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
637 0 : WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
638 0 : convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
639 0 : CPWARN(msg)
640 : END IF
641 0 : ELSE IF (INDEX(line, "eV") /= 0) THEN
642 0 : convttemp = cp_unit_from_cp2k(ttemp, "K")/11604.505_dp
643 0 : IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
644 0 : WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
645 0 : convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
646 0 : CPWARN(msg)
647 : END IF
648 0 : ELSE IF (INDEX(line, "atomic energy units") /= 0) THEN
649 0 : convttemp = ttemp
650 0 : IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
651 0 : WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
652 0 : convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
653 0 : CPWARN(msg)
654 : END IF
655 : ELSE
656 0 : WRITE (UNIT=msg, FMT=*) "Unknown PIGLET matrix temperature. Assuming a.u."
657 0 : CPWARN(msg)
658 0 : convttemp = ttemp
659 0 : IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
660 0 : WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
661 0 : convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
662 0 : CPWARN(msg)
663 : END IF
664 : END IF
665 :
666 1 : END SUBROUTINE check_temperature
667 :
668 : ! ***************************************************************************
669 : !> \brief returns the piglet kinetic energy contribution
670 : !> \param pint_env ...
671 : !> \author Felix Uhl
672 : ! **************************************************************************************************
673 12 : ELEMENTAL SUBROUTINE pint_calc_piglet_energy(pint_env)
674 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
675 :
676 12 : IF (ASSOCIATED(pint_env%piglet_therm)) THEN
677 12 : pint_env%e_piglet = pint_env%piglet_therm%thermostat_energy
678 : ELSE
679 0 : pint_env%e_piglet = 0.0d0
680 : END IF
681 :
682 12 : 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 48 : 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 48 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
700 : REAL(KIND=dp), DIMENSION(1) :: tmplwork
701 96 : REAL(KIND=dp), DIMENSION(n) :: eigval
702 48 : 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 480 : eigval(:) = 0.0_dp
716 4368 : A(:, :) = 0.0_dp
717 4368 : 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 48 : info) ! information about success
729 :
730 48 : lwork = INT(tmplwork(1) + 0.5_dp)
731 144 : ALLOCATE (work(lwork))
732 48 : 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 48 : info) ! information about success
743 48 : DEALLOCATE (work)
744 : ! A-matrix now contains the eigenvectors
745 :
746 4368 : S(:, :) = 0.0_dp
747 480 : DO i = 1, n
748 : ! In case that numerics made some eigenvalues negative
749 480 : IF (eigval(i) > 0.0_dp) THEN
750 432 : S(i, i) = SQRT(eigval(i))
751 : END IF
752 : END DO
753 :
754 4368 : 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 48 : 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 48 : n) ! leading dimension of S-matrix as declared
784 :
785 48 : END SUBROUTINE sqrt_pos_def_mat
786 :
787 : END MODULE pint_piglet
|
