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 QTB thermostat to PI runs.
10 : !> Based on the PILE implementation from Felix Uhl (pint_pile.F)
11 : !> \author Fabien Brieuc
12 : !> \par History
13 : !> 02.2018 created [Fabien Brieuc]
14 : ! **************************************************************************************************
15 : MODULE pint_qtb
16 : USE cp_files, ONLY: open_file
17 : USE cp_log_handling, ONLY: cp_get_default_logger,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: debug_print_level,&
20 : silent_print_level
21 : USE fft_tools, ONLY: FWFFT,&
22 : fft_1d_many,&
23 : fft_alloc,&
24 : fft_dealloc
25 : USE input_constants, ONLY: propagator_rpmd
26 : USE input_section_types, ONLY: section_vals_get,&
27 : section_vals_get_subs_vals,&
28 : section_vals_type,&
29 : section_vals_val_get
30 : USE kinds, ONLY: dp
31 : USE mathconstants, ONLY: pi,&
32 : twopi
33 : USE message_passing, ONLY: mp_para_env_type
34 : USE parallel_rng_types, ONLY: GAUSSIAN,&
35 : rng_record_length,&
36 : rng_stream_type,&
37 : rng_stream_type_from_record
38 : USE pint_io, ONLY: pint_write_line
39 : USE pint_types, ONLY: normalmode_env_type,&
40 : pint_env_type,&
41 : qtb_therm_type
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : PUBLIC :: pint_qtb_step, &
49 : pint_qtb_init, &
50 : pint_qtb_release, &
51 : pint_calc_qtb_energy
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_qtb'
54 :
55 : CONTAINS
56 :
57 : ! ***************************************************************************
58 : !> \brief initializes the data for a QTB run
59 : !> \brief ...
60 : !> \param qtb_therm ...
61 : !> \param pint_env ...
62 : !> \param normalmode_env ...
63 : !> \param section ...
64 : ! **************************************************************************************************
65 6 : SUBROUTINE pint_qtb_init(qtb_therm, pint_env, normalmode_env, section)
66 : TYPE(qtb_therm_type), POINTER :: qtb_therm
67 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
68 : TYPE(normalmode_env_type), POINTER :: normalmode_env
69 : TYPE(section_vals_type), POINTER :: section
70 :
71 : CHARACTER(LEN=rng_record_length) :: rng_record
72 : INTEGER :: i, j, p
73 : LOGICAL :: restart
74 : REAL(KIND=dp) :: dti2, ex
75 : REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed
76 : TYPE(section_vals_type), POINTER :: rng_section
77 :
78 6 : IF (pint_env%propagator%prop_kind /= propagator_rpmd) THEN
79 0 : CPABORT("QTB is designed to work with the RPMD propagator only")
80 : END IF
81 :
82 6 : pint_env%e_qtb = 0.0_dp
83 150 : ALLOCATE (qtb_therm)
84 : qtb_therm%thermostat_energy = 0.0_dp
85 :
86 : !Get input parameters
87 6 : CALL section_vals_val_get(section, "TAU", r_val=qtb_therm%tau)
88 6 : CALL section_vals_val_get(section, "LAMBDA", r_val=qtb_therm%lamb)
89 6 : CALL section_vals_val_get(section, "TAUCUT", r_val=qtb_therm%taucut)
90 6 : CALL section_vals_val_get(section, "LAMBCUT", r_val=qtb_therm%lambcut)
91 6 : CALL section_vals_val_get(section, "FP", i_val=qtb_therm%fp)
92 6 : CALL section_vals_val_get(section, "NF", i_val=qtb_therm%nf)
93 6 : CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=qtb_therm%thermostat_energy)
94 :
95 6 : p = pint_env%p
96 6 : dti2 = 0.5_dp*pint_env%dt
97 18 : ALLOCATE (qtb_therm%c1(p))
98 12 : ALLOCATE (qtb_therm%c2(p))
99 12 : ALLOCATE (qtb_therm%g_fric(p))
100 24 : ALLOCATE (qtb_therm%massfact(p, pint_env%ndim))
101 :
102 : !Initialize everything
103 6 : qtb_therm%g_fric(1) = 1.0_dp/qtb_therm%tau
104 24 : DO i = 2, p
105 : qtb_therm%g_fric(i) = SQRT((1.d0/qtb_therm%tau)**2 + (qtb_therm%lamb)**2* &
106 24 : normalmode_env%lambda(i))
107 : END DO
108 30 : DO i = 1, p
109 24 : ex = -dti2*qtb_therm%g_fric(i)
110 24 : qtb_therm%c1(i) = EXP(ex)
111 24 : ex = qtb_therm%c1(i)*qtb_therm%c1(i)
112 30 : qtb_therm%c2(i) = SQRT(1.0_dp - ex)
113 : END DO
114 27654 : DO j = 1, pint_env%ndim
115 138246 : DO i = 1, pint_env%p
116 138240 : qtb_therm%massfact(i, j) = SQRT(1.0_dp/pint_env%mass_fict(i, j))
117 : END DO
118 : END DO
119 :
120 : !prepare Random number generator
121 6 : NULLIFY (rng_section)
122 : rng_section => section_vals_get_subs_vals(section, &
123 6 : subsection_name="RNG_INIT")
124 6 : CALL section_vals_get(rng_section, explicit=restart)
125 6 : IF (restart) THEN
126 : CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
127 2 : i_rep_val=1, c_val=rng_record)
128 2 : qtb_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
129 : ELSE
130 36 : initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
131 : qtb_therm%gaussian_rng_stream = rng_stream_type( &
132 : name="qtb_rng_gaussian", distribution_type=GAUSSIAN, &
133 : extended_precision=.TRUE., &
134 4 : seed=initial_seed)
135 : END IF
136 :
137 : !Initialization of the QTB random forces
138 6 : CALL pint_qtb_forces_init(pint_env, normalmode_env, qtb_therm, restart)
139 :
140 6 : END SUBROUTINE pint_qtb_init
141 :
142 : ! **************************************************************************************************
143 : !> \brief ...
144 : !> \param vold ...
145 : !> \param vnew ...
146 : !> \param p ...
147 : !> \param ndim ...
148 : !> \param masses ...
149 : !> \param qtb_therm ...
150 : ! **************************************************************************************************
151 60 : SUBROUTINE pint_qtb_step(vold, vnew, p, ndim, masses, qtb_therm)
152 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vold, vnew
153 : INTEGER, INTENT(IN) :: p, ndim
154 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: masses
155 : TYPE(qtb_therm_type), POINTER :: qtb_therm
156 :
157 : CHARACTER(len=*), PARAMETER :: routineN = 'pint_qtb_step'
158 :
159 : INTEGER :: handle, i, ibead, idim
160 : REAL(KIND=dp) :: delta_ekin
161 :
162 60 : CALL timeset(routineN, handle)
163 60 : delta_ekin = 0.0_dp
164 :
165 : !update random forces
166 300 : DO ibead = 1, p
167 240 : qtb_therm%cpt(ibead) = qtb_therm%cpt(ibead) + 1
168 : !new random forces at every qtb_therm%step
169 300 : IF (qtb_therm%cpt(ibead) == 2*qtb_therm%step(ibead)) THEN
170 0 : IF (ibead == 1) THEN
171 : !update the rng status
172 0 : DO i = 1, qtb_therm%nf - 1
173 0 : qtb_therm%rng_status(i) = qtb_therm%rng_status(i + 1)
174 : END DO
175 0 : CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(qtb_therm%nf))
176 : END IF
177 0 : DO idim = 1, ndim
178 : !update random numbers
179 0 : DO i = 1, qtb_therm%nf - 1
180 0 : qtb_therm%r(i, ibead, idim) = qtb_therm%r(i + 1, ibead, idim)
181 : END DO
182 0 : qtb_therm%r(qtb_therm%nf, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
183 : !compute new random force through the convolution product
184 0 : qtb_therm%rf(ibead, idim) = 0.0_dp
185 0 : DO i = 1, qtb_therm%nf
186 : qtb_therm%rf(ibead, idim) = qtb_therm%rf(ibead, idim) + &
187 0 : qtb_therm%h(i, ibead)*qtb_therm%r(i, ibead, idim)
188 : END DO
189 : END DO
190 0 : qtb_therm%cpt(ibead) = 0
191 : END IF
192 : END DO
193 :
194 : !perform MD step
195 276540 : DO idim = 1, ndim
196 1382460 : DO ibead = 1, p
197 : vnew(ibead, idim) = qtb_therm%c1(ibead)*vold(ibead, idim) + &
198 : qtb_therm%massfact(ibead, idim)*qtb_therm%c2(ibead)* &
199 1105920 : qtb_therm%rf(ibead, idim)
200 : delta_ekin = delta_ekin + masses(ibead, idim)*( &
201 : vnew(ibead, idim)*vnew(ibead, idim) - &
202 1382400 : vold(ibead, idim)*vold(ibead, idim))
203 : END DO
204 : END DO
205 :
206 60 : qtb_therm%thermostat_energy = qtb_therm%thermostat_energy - 0.5_dp*delta_ekin
207 :
208 60 : CALL timestop(handle)
209 60 : END SUBROUTINE pint_qtb_step
210 :
211 : ! ***************************************************************************
212 : !> \brief releases the qtb environment
213 : !> \param qtb_therm qtb data to be released
214 : ! **************************************************************************************************
215 6 : SUBROUTINE pint_qtb_release(qtb_therm)
216 :
217 : TYPE(qtb_therm_type), INTENT(INOUT) :: qtb_therm
218 :
219 6 : DEALLOCATE (qtb_therm%c1)
220 6 : DEALLOCATE (qtb_therm%c2)
221 6 : DEALLOCATE (qtb_therm%g_fric)
222 6 : DEALLOCATE (qtb_therm%massfact)
223 6 : DEALLOCATE (qtb_therm%rf)
224 6 : DEALLOCATE (qtb_therm%h)
225 6 : DEALLOCATE (qtb_therm%r)
226 6 : DEALLOCATE (qtb_therm%cpt)
227 6 : DEALLOCATE (qtb_therm%step)
228 6 : DEALLOCATE (qtb_therm%rng_status)
229 :
230 6 : END SUBROUTINE pint_qtb_release
231 :
232 : ! ***************************************************************************
233 : !> \brief returns the qtb kinetic energy contribution
234 : !> \param pint_env ...
235 : ! **************************************************************************************************
236 36 : SUBROUTINE pint_calc_qtb_energy(pint_env)
237 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
238 :
239 36 : IF (ASSOCIATED(pint_env%qtb_therm)) THEN
240 36 : pint_env%e_qtb = pint_env%qtb_therm%thermostat_energy
241 : END IF
242 :
243 36 : END SUBROUTINE pint_calc_qtb_energy
244 :
245 : ! ***************************************************************************
246 : !> \brief initialize the QTB random forces
247 : !> \param pint_env ...
248 : !> \param normalmode_env ...
249 : !> \param qtb_therm ...
250 : !> \param restart ...
251 : !> \author Fabien Brieuc
252 : ! **************************************************************************************************
253 6 : SUBROUTINE pint_qtb_forces_init(pint_env, normalmode_env, qtb_therm, restart)
254 : TYPE(pint_env_type), INTENT(IN) :: pint_env
255 : TYPE(normalmode_env_type), POINTER :: normalmode_env
256 : TYPE(qtb_therm_type), POINTER :: qtb_therm
257 : LOGICAL :: restart
258 :
259 : CHARACTER(len=*), PARAMETER :: routineN = 'pint_qtb_forces_init'
260 :
261 : COMPLEX(KIND=dp) :: tmp1
262 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:), &
263 6 : POINTER :: filter_in, filter_out
264 : INTEGER :: handle, i, ibead, idim, log_unit, ndim, &
265 : nf, p, print_level, step
266 : REAL(KIND=dp) :: aa, bb, correct, dt, dw, fcut, h, kT, &
267 : tmp, w
268 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fp
269 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: fp1
270 : TYPE(cp_logger_type), POINTER :: logger
271 : TYPE(mp_para_env_type), POINTER :: para_env
272 :
273 6 : CALL timeset(routineN, handle)
274 :
275 6 : p = pint_env%p
276 6 : ndim = pint_env%ndim
277 6 : dt = pint_env%dt
278 6 : IF (MOD(qtb_therm%nf, 2) /= 0) qtb_therm%nf = qtb_therm%nf + 1
279 6 : nf = qtb_therm%nf
280 :
281 6 : para_env => pint_env%logger%para_env
282 :
283 18 : ALLOCATE (qtb_therm%rng_status(nf))
284 24 : ALLOCATE (qtb_therm%h(nf, p))
285 18 : ALLOCATE (qtb_therm%step(p))
286 :
287 : !initialize random forces on ionode only
288 6 : IF (para_env%is_source()) THEN
289 :
290 3 : NULLIFY (logger)
291 3 : logger => cp_get_default_logger()
292 3 : print_level = logger%iter_info%print_level
293 :
294 : !physical temperature (T) not the simulation one (TxP)
295 3 : kT = pint_env%kT*pint_env%propagator%temp_sim2phys
296 :
297 9 : ALLOCATE (fp(nf/2))
298 6 : CALL fft_alloc(filter_in, [nf])
299 6 : CALL fft_alloc(filter_out, [nf])
300 :
301 3 : IF (print_level == debug_print_level) THEN
302 : !create log file if print_level is debug
303 : CALL open_file(file_name=TRIM(logger%iter_info%project_name)//".qtbLog", &
304 0 : file_action="WRITE", file_status="UNKNOWN", unit_number=log_unit)
305 0 : WRITE (log_unit, '(A)') ' # Log file for the QTB random forces generation'
306 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
307 0 : WRITE (log_unit, '(A,I5)') ' # Number of beads P = ', p
308 0 : WRITE (log_unit, '(A,I6)') ' # Number of dimension 3*N = ', ndim
309 0 : WRITE (log_unit, '(A,I4)') ' # Number of filter parameters Nf=', nf
310 : END IF
311 :
312 15 : DO ibead = 1, p
313 : !fcut is adapted to the NM freq.
314 : !Note that lambda is the angular free ring freq. squared
315 : fcut = SQRT((1.d0/qtb_therm%taucut)**2 + (qtb_therm%lambcut)**2* &
316 12 : normalmode_env%lambda(ibead))
317 12 : fcut = fcut/twopi
318 : !new random forces are drawn every step
319 12 : qtb_therm%step(ibead) = NINT(1.0_dp/(2.0_dp*fcut*dt))
320 12 : IF (qtb_therm%step(ibead) == 0) qtb_therm%step(ibead) = 1
321 12 : step = qtb_therm%step(ibead)
322 : !effective timestep h = step*dt = 1/(2*fcut)
323 12 : h = step*dt
324 : !angular freq. step - dw = 2*pi/(nf*h) = 2*wcut/nf
325 12 : dw = twopi/(nf*h)
326 :
327 : !generate f_P function
328 12 : IF (qtb_therm%fp == 0) THEN
329 4 : CALL pint_qtb_computefp0(pint_env, fp, fp1, dw, aa, bb, log_unit, ibead, print_level)
330 : ELSE
331 8 : CALL pint_qtb_computefp1(pint_env, fp, fp1, dw, aa, bb, log_unit, ibead, print_level)
332 : END IF
333 780 : fp = p*kT*fp ! fp is now in cp2k energy units
334 :
335 12 : IF (print_level == debug_print_level) THEN
336 0 : WRITE (log_unit, '(A,I4,A)') ' # -------- NM ', ibead, ' --------'
337 0 : WRITE (log_unit, '(A,I4,A)') ' # New random forces every ', step, ' MD steps'
338 0 : WRITE (log_unit, '(A,ES13.3,A)') ' # Angular cutoff freq. = ', twopi*fcut*4.1341e4_dp, ' rad/ps'
339 0 : WRITE (log_unit, '(A,ES13.3,A)') ' # Free ring polymer angular freq.= ', &
340 0 : SQRT(normalmode_env%lambda(ibead))*4.1341e4_dp, ' rad/ps'
341 0 : WRITE (log_unit, '(A,ES13.3,A)') ' # Friction coeff. = ', qtb_therm%g_fric(ibead)*4.1341e4_dp, ' THz'
342 0 : WRITE (log_unit, '(A,ES13.3,A)') ' # Angular frequency step dw = ', dw*4.1341e4_dp, ' rad/ps'
343 : END IF
344 :
345 : !compute the filter in Fourier space
346 12 : IF (p == 1) THEN
347 0 : filter_in(1) = SQRT(kT)*(1.0_dp, 0.0_dp)
348 12 : ELSE IF (qtb_therm%fp == 1 .AND. ibead == 1) THEN
349 2 : filter_in(1) = SQRT(p*kT)*(1.0_dp, 0.0_dp)
350 : ELSE
351 10 : filter_in(1) = SQRT(p*kT*fp1(1))*(1.0_dp, 0.0_dp)
352 : END IF
353 780 : DO i = 1, nf/2
354 768 : w = i*dw
355 768 : tmp = 0.5_dp*w*h
356 768 : correct = SIN(tmp)/tmp
357 768 : filter_in(i + 1) = SQRT(fp(i))/correct*(1.0_dp, 0.0_dp)
358 780 : filter_in(nf - i + 1) = CONJG(filter_in(i + 1))
359 : END DO
360 :
361 : !compute the filter in time space - FFT
362 12 : CALL pint_qtb_fft(filter_in, filter_out, nf)
363 : !reordering + normalisation
364 : !normalisation : 1/nf comes from the DFT, 1/sqrt(step) is to
365 : !take into account the effective timestep h = step*dt and
366 : !1/sqrt(2.0_dp) is to take into account the fact that the
367 : !same random force is used for the two thermostat "half-steps"
368 780 : DO i = 1, nf/2
369 768 : tmp1 = filter_out(i)/(nf*SQRT(2.0_dp*step))
370 768 : filter_out(i) = filter_out(nf/2 + i)/(nf*SQRT(2.0_dp*step))
371 780 : filter_out(nf/2 + i) = tmp1
372 : END DO
373 :
374 1551 : DO i = 1, nf
375 1548 : qtb_therm%h(i, ibead) = REAL(filter_out(i), dp)
376 : END DO
377 : END DO
378 :
379 3 : CALL fft_dealloc(filter_in)
380 3 : CALL fft_dealloc(filter_out)
381 3 : DEALLOCATE (fp)
382 3 : IF (p > 1) DEALLOCATE (fp1)
383 : END IF
384 :
385 6198 : CALL para_env%bcast(qtb_therm%h)
386 54 : CALL para_env%bcast(qtb_therm%step)
387 :
388 30 : ALLOCATE (qtb_therm%r(nf, p, ndim))
389 12 : ALLOCATE (qtb_therm%cpt(p))
390 24 : ALLOCATE (qtb_therm%rf(p, ndim))
391 :
392 6 : IF (restart) THEN
393 2 : CALL pint_qtb_restart(pint_env, qtb_therm)
394 : ELSE
395 : !update the rng status
396 516 : DO i = 1, qtb_therm%nf
397 516 : CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(i))
398 : END DO
399 : !if no restart then initialize random numbers from scratch
400 20 : qtb_therm%cpt = 0
401 18436 : DO idim = 1, ndim
402 92164 : DO ibead = 1, p
403 9529344 : DO i = 1, nf
404 9510912 : qtb_therm%r(i, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
405 : END DO
406 : END DO
407 : END DO
408 : END IF
409 :
410 : !compute the first random forces
411 27654 : DO idim = 1, ndim
412 138246 : DO ibead = 1, p
413 110592 : qtb_therm%rf(ibead, idim) = 0.0_dp
414 14294016 : DO i = 1, nf
415 : qtb_therm%rf(ibead, idim) = qtb_therm%rf(ibead, idim) + &
416 14266368 : qtb_therm%h(i, ibead)*qtb_therm%r(i, ibead, idim)
417 : END DO
418 : END DO
419 : END DO
420 :
421 6 : CALL timestop(handle)
422 12 : END SUBROUTINE pint_qtb_forces_init
423 :
424 : ! ***************************************************************************
425 : !> \brief control the generation of the first random forces in the case
426 : !> of a restart
427 : !> \param pint_env ...
428 : !> \param qtb_therm ...
429 : !> \author Fabien Brieuc
430 : ! **************************************************************************************************
431 2 : SUBROUTINE pint_qtb_restart(pint_env, qtb_therm)
432 : TYPE(pint_env_type), INTENT(IN) :: pint_env
433 : TYPE(qtb_therm_type), POINTER :: qtb_therm
434 :
435 : INTEGER :: begin, i, ibead, idim, istep
436 :
437 : begin = pint_env%first_step - MOD(pint_env%first_step, qtb_therm%step(1)) - &
438 2 : (qtb_therm%nf - 1)*qtb_therm%step(1)
439 :
440 2 : IF (begin <= 0) THEN
441 10 : qtb_therm%cpt = 0
442 : !update the rng status
443 258 : DO i = 1, qtb_therm%nf
444 258 : CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(i))
445 : END DO
446 : !first random numbers initialized from scratch
447 9218 : DO idim = 1, pint_env%ndim
448 46082 : DO ibead = 1, pint_env%p
449 4764672 : DO i = 1, qtb_therm%nf
450 4755456 : qtb_therm%r(i, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
451 : END DO
452 : END DO
453 : END DO
454 : begin = 1
455 : ELSE
456 0 : qtb_therm%cpt(1) = 2*(qtb_therm%step(1) - 1)
457 0 : DO ibead = 2, pint_env%p
458 0 : qtb_therm%cpt(ibead) = 2*MOD(begin - 1, qtb_therm%step(ibead))
459 : END DO
460 : END IF
461 :
462 : !from istep = 1,2*(the last previous MD step - begin) because
463 : !the thermostat step is called two times per MD step
464 : !DO istep = 2*begin, 2*pint_env%first_step
465 22 : DO istep = 1, 2*(pint_env%first_step - begin + 1)
466 102 : DO ibead = 1, pint_env%p
467 80 : qtb_therm%cpt(ibead) = qtb_therm%cpt(ibead) + 1
468 : !new random forces at every qtb_therm%step
469 100 : IF (qtb_therm%cpt(ibead) == 2*qtb_therm%step(ibead)) THEN
470 0 : IF (ibead == 1) THEN
471 : !update the rng status
472 0 : DO i = 1, qtb_therm%nf - 1
473 0 : qtb_therm%rng_status(i) = qtb_therm%rng_status(i + 1)
474 : END DO
475 0 : CALL qtb_therm%gaussian_rng_stream%dump(qtb_therm%rng_status(qtb_therm%nf))
476 : END IF
477 0 : DO idim = 1, pint_env%ndim
478 : !update random numbers
479 0 : DO i = 1, qtb_therm%nf - 1
480 0 : qtb_therm%r(i, ibead, idim) = qtb_therm%r(i + 1, ibead, idim)
481 : END DO
482 0 : qtb_therm%r(qtb_therm%nf, ibead, idim) = qtb_therm%gaussian_rng_stream%next()
483 : END DO
484 0 : qtb_therm%cpt(ibead) = 0
485 : END IF
486 : END DO
487 : END DO
488 :
489 2 : END SUBROUTINE pint_qtb_restart
490 :
491 : ! ***************************************************************************
492 : !> \brief compute the f_P^(0) function necessary for coupling QTB with PIMD
493 : !> \param pint_env ...
494 : !> \param fp stores the computed function on the grid used for the generation
495 : !> of the filter h
496 : !> \param fp1 stores the computed function on an larger and finer grid
497 : !> \param dw angular frequency step
498 : !> \param aa ...
499 : !> \param bb ...
500 : !> \param log_unit ...
501 : !> \param ibead ...
502 : !> \param print_level ...
503 : !> \author Fabien Brieuc
504 : ! **************************************************************************************************
505 4 : SUBROUTINE pint_qtb_computefp0(pint_env, fp, fp1, dw, aa, bb, log_unit, ibead, print_level)
506 : TYPE(pint_env_type), INTENT(IN) :: pint_env
507 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: fp
508 : REAL(KIND=dp), DIMENSION(:), POINTER :: fp1
509 : REAL(KIND=dp), INTENT(IN) :: dw, aa, bb
510 : INTEGER, INTENT(IN) :: log_unit, ibead, print_level
511 :
512 : CHARACTER(len=200) :: line
513 : INTEGER :: i, j, k, n, niter, nx, p
514 4 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: kk
515 : REAL(KIND=dp) :: dx, dx1, err, fprev, hbokT, malpha, op, &
516 : r2, tmp, w, x1, xmax, xmin, xx
517 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: h, x, x2
518 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fpxk, xk, xk2
519 :
520 4 : n = SIZE(fp)
521 4 : p = pint_env%p
522 :
523 : !using the physical temperature (T) not the simulation one (TxP)
524 4 : hbokT = 1.0_dp/(pint_env%kT*pint_env%propagator%temp_sim2phys)
525 :
526 : !P = 1 : standard QTB
527 : !fp = theta(w, T) / kT
528 4 : IF (p == 1) THEN
529 0 : DO j = 1, n
530 0 : w = j*dw
531 0 : tmp = hbokT*w
532 0 : fp(j) = tmp*(0.5_dp + 1.0_dp/(EXP(tmp) - 1.0_dp))
533 : END DO
534 :
535 0 : IF (print_level == debug_print_level) THEN
536 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
537 0 : WRITE (log_unit, '(A)') ' # computed fp^(0) function'
538 0 : WRITE (log_unit, '(A)') ' # i, w(a.u.), fp'
539 0 : DO j = 1, n
540 0 : WRITE (log_unit, *) j, j*dw, j*0.5_dp*hbokt*dw, fp(j)
541 : END DO
542 : END IF
543 : ! P > 1: QTB-PIMD
544 : ELSE
545 : !**** initialization ****
546 4 : dx1 = 0.5_dp*hbokt*dw
547 4 : xmin = 1.0e-7_dp !these values allows for an acceptable
548 4 : dx = 0.05_dp !ratio between accuracy, computing time and
549 4 : xmax = 10000.0_dp !memory requirement - tested for P up to 1024
550 : nx = INT((xmax - xmin)/dx) + 1
551 4 : nx = nx + nx/5 !add 20% points to avoid any problems at the end
552 : !of the interval (probably unnecessary)
553 4 : IF (ibead == 1) THEN
554 1 : op = 1.0_dp/p
555 1 : malpha = op !mixing parameter alpha = 1/P
556 1 : niter = 30 !30 iterations are enough to converge
557 :
558 1 : IF (print_level == debug_print_level) THEN
559 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
560 0 : WRITE (log_unit, '(A)') ' # computing fp^(0) function'
561 0 : WRITE (log_unit, '(A)') ' # parameters used:'
562 0 : WRITE (log_unit, '(A,ES13.3)') ' # dx = ', dx
563 0 : WRITE (log_unit, '(A,ES13.3)') ' # xmin = ', xmin
564 0 : WRITE (log_unit, '(A,ES13.3)') ' # xmax = ', xmax
565 0 : WRITE (log_unit, '(A,I8,I8)') ' # nx, n = ', nx, n
566 : END IF
567 :
568 1 : ALLOCATE (x(nx))
569 1 : ALLOCATE (x2(nx))
570 1 : ALLOCATE (h(nx))
571 1 : ALLOCATE (fp1(nx))
572 4 : ALLOCATE (xk(p - 1, nx))
573 3 : ALLOCATE (xk2(p - 1, nx))
574 4 : ALLOCATE (kk(p - 1, nx))
575 3 : ALLOCATE (fpxk(p - 1, nx))
576 :
577 : ! initialize fp(x)
578 : ! fp1 = fp(x) = h(x/P)
579 : ! fpxk = fp(xk) = h(xk/P)
580 240001 : DO j = 1, nx
581 240000 : x(j) = xmin + (j - 1)*dx
582 240000 : x2(j) = x(j)**2
583 240000 : h(j) = x(j)/TANH(x(j))
584 240000 : IF (x(j) <= 1.0e-10_dp) h(j) = 1.0_dp
585 240000 : fp1(j) = op*x(j)/TANH(x(j)*op)
586 240000 : IF (x(j)*op <= 1.0e-10_dp) fp1(j) = 1.0_dp
587 960001 : DO k = 1, p - 1
588 720000 : xk2(k, j) = x2(j) + (p*SIN(k*pi*op))**2
589 720000 : xk(k, j) = SQRT(xk2(k, j))
590 720000 : kk(k, j) = NINT((xk(k, j) - xmin)/dx) + 1
591 720000 : fpxk(k, j) = xk(k, j)*op/TANH(xk(k, j)*op)
592 960000 : IF (xk(k, j)*op <= 1.0e-10_dp) fpxk(k, j) = 1.0_dp
593 : END DO
594 : END DO
595 :
596 : ! **** resolution ****
597 : ! compute fp(x)
598 31 : DO i = 1, niter
599 30 : err = 0.0_dp
600 7200030 : DO j = 1, nx
601 : tmp = 0.0_dp
602 28800000 : DO k = 1, p - 1
603 28800000 : tmp = tmp + fpxk(k, j)*x2(j)/xk2(k, j)
604 : END DO
605 7200000 : fprev = fp1(j)
606 7200000 : fp1(j) = malpha*(h(j) - tmp) + (1.0_dp - malpha)*fp1(j)
607 7200030 : IF (j <= n) err = err + ABS(1.0_dp - fp1(j)/fprev) ! compute "errors"
608 : END DO
609 30 : err = err/n
610 :
611 : ! Linear regression on the last 20% of the F_P function
612 30 : CALL pint_qtb_linreg(fp1(8*nx/10:nx), x(8*nx/10:nx), aa, bb, r2, log_unit, print_level)
613 :
614 : ! compute the new F_P(xk*sqrt(P))
615 : ! through linear interpolation
616 : ! or linear extrapolation if outside of the range
617 7200031 : DO j = 1, nx
618 28800030 : DO k = 1, p - 1
619 28800000 : IF (kk(k, j) < nx) THEN
620 : fpxk(k, j) = fp1(kk(k, j)) + (fp1(kk(k, j) + 1) - fp1(kk(k, j)))/dx* &
621 21599910 : (xk(k, j) - x(kk(k, j)))
622 : ELSE
623 90 : fpxk(k, j) = aa*xk(k, j) + bb
624 : END IF
625 : END DO
626 : END DO
627 : END DO
628 :
629 1 : IF (print_level == debug_print_level) THEN
630 : ! **** tests ****
631 0 : WRITE (log_unit, '(A,ES9.3)') ' # average error during computation: ', err
632 0 : WRITE (log_unit, '(A,ES9.3)') ' # slope of F_P at high freq. - theoretical: ', op
633 0 : WRITE (log_unit, '(A,ES9.3)') ' # slope of F_P at high freq. - calculated: ', aa
634 0 : WRITE (log_unit, '(A,F6.3)') ' # F_P at zero freq. - theoretical: ', 1.0_dp
635 0 : WRITE (log_unit, '(A,F6.3)') ' # F_P at zero freq. - calculated: ', fp1(1)
636 1 : ELSE IF (print_level > silent_print_level) THEN
637 1 : CALL pint_write_line("QTB| Initialization of random forces using fP0 function")
638 1 : CALL pint_write_line("QTB| Computation of fP0 function")
639 1 : WRITE (line, '(A,ES9.3)') 'QTB| average error ', err
640 1 : CALL pint_write_line(TRIM(line))
641 1 : WRITE (line, '(A,ES9.3)') 'QTB| slope at high frequency - theoretical: ', op
642 1 : CALL pint_write_line(TRIM(line))
643 1 : WRITE (line, '(A,ES9.3)') 'QTB| slope at high frequency - calculated: ', aa
644 1 : CALL pint_write_line(TRIM(line))
645 1 : WRITE (line, '(A,F6.3)') 'QTB| value at zero frequency - theoretical: ', 1.0_dp
646 1 : CALL pint_write_line(TRIM(line))
647 1 : WRITE (line, '(A,F6.3)') 'QTB| value at zero frequency - calculated: ', fp1(1)
648 1 : CALL pint_write_line(TRIM(line))
649 : END IF
650 :
651 1 : IF (print_level == debug_print_level) THEN
652 : ! **** write solution ****
653 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
654 0 : WRITE (log_unit, '(A)') ' # computed fp function'
655 0 : WRITE (log_unit, '(A)') ' # i, w(a.u.), x, fp'
656 0 : DO j = 1, nx
657 0 : WRITE (log_unit, *) j, j*dw, xmin + (j - 1)*dx, fp1(j)
658 : END DO
659 : END IF
660 :
661 1 : DEALLOCATE (x)
662 1 : DEALLOCATE (x2)
663 1 : DEALLOCATE (h)
664 1 : DEALLOCATE (xk)
665 1 : DEALLOCATE (xk2)
666 1 : DEALLOCATE (kk)
667 1 : DEALLOCATE (fpxk)
668 : END IF
669 :
670 : ! compute values of fP on the grid points for the current NM
671 : ! through linear interpolation / regression
672 260 : DO j = 1, n
673 256 : x1 = j*dx1
674 256 : k = NINT((x1 - xmin)/dx) + 1
675 260 : IF (k > nx) THEN
676 0 : fp(j) = aa*x1 + bb
677 256 : ELSE IF (k <= 0) THEN
678 0 : CALL pint_write_line("QTB| error in fp computation x < xmin")
679 0 : CPABORT("Error in fp computation (x < xmin) in initialization of QTB random forces")
680 : ELSE
681 256 : xx = xmin + (k - 1)*dx
682 256 : IF (x1 > xx) THEN
683 117 : fp(j) = fp1(k) + (fp1(k + 1) - fp1(k))/dx*(x1 - xx)
684 : ELSE
685 139 : fp(j) = fp1(k) + (fp1(k) - fp1(k - 1))/dx*(x1 - xx)
686 : END IF
687 : END IF
688 : END DO
689 :
690 : END IF
691 :
692 4 : END SUBROUTINE pint_qtb_computefp0
693 :
694 : ! ***************************************************************************
695 : !> \brief compute the f_P^(1) function necessary for coupling QTB with PIMD
696 : !> \param pint_env ...
697 : !> \param fp stores the computed function on the grid used for the generation
698 : !> of the filter h
699 : !> \param fp1 stores the computed function on an larger and finer grid
700 : !> \param dw angular frequency step
701 : !> \param aa ...
702 : !> \param bb ...
703 : !> \param log_unit ...
704 : !> \param ibead ...
705 : !> \param print_level ...
706 : !> \author Fabien Brieuc
707 : ! **************************************************************************************************
708 8 : SUBROUTINE pint_qtb_computefp1(pint_env, fp, fp1, dw, aa, bb, log_unit, ibead, print_level)
709 : TYPE(pint_env_type), INTENT(IN) :: pint_env
710 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: fp
711 : REAL(KIND=dp), DIMENSION(:), POINTER :: fp1
712 : REAL(KIND=dp) :: dw, aa, bb
713 : INTEGER, INTENT(IN) :: log_unit, ibead, print_level
714 :
715 : CHARACTER(len=200) :: line
716 : INTEGER :: i, j, k, n, niter, nx, p
717 8 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: kk
718 : REAL(KIND=dp) :: dx, dx1, err, fprev, hbokT, malpha, op, &
719 : op1, r2, tmp, tmp1, xmax, xmin, xx
720 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: h, x, x2
721 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fpxk, xk, xk2
722 :
723 8 : n = SIZE(fp)
724 8 : p = pint_env%p
725 :
726 : !using the physical temperature (T) not the simulation one (TxP)
727 8 : hbokT = 1.0_dp/(pint_env%kT*pint_env%propagator%temp_sim2phys)
728 :
729 : !Centroid NM (ibead=1) : classical
730 : !fp = 1
731 8 : IF (ibead == 1) THEN
732 130 : DO j = 1, n
733 130 : fp(j) = 1.0_dp
734 : END DO
735 : ELSE
736 : !**** initialization ****
737 6 : dx1 = 0.5_dp*hbokt*dw
738 6 : xmin = 1.0e-3_dp !these values allows for an acceptable
739 6 : dx = 0.05_dp !ratio between accuracy, computing time and
740 6 : xmax = 10000.0_dp !memory requirement - tested for P up to 1024
741 : nx = INT((xmax - xmin)/dx) + 1
742 6 : nx = nx + nx/5 !add 20% points to avoid problem at the end
743 : !of the interval (probably unnecessary)
744 6 : op = 1.0_dp/p
745 6 : IF (ibead == 2) THEN
746 2 : op1 = 1.0_dp/(p - 1)
747 2 : malpha = op !mixing parameter alpha = 1/P
748 2 : niter = 40 !40 iterations are enough to converge
749 :
750 2 : IF (print_level == debug_print_level) THEN
751 : ! **** write solution ****
752 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
753 0 : WRITE (log_unit, '(A)') ' # computing fp^(1) function'
754 0 : WRITE (log_unit, '(A)') ' # parameters used:'
755 0 : WRITE (log_unit, '(A,ES13.3)') ' # dx = ', dx
756 0 : WRITE (log_unit, '(A,ES13.3)') ' # xmin = ', xmin
757 0 : WRITE (log_unit, '(A,ES13.3)') ' # xmax = ', xmax
758 0 : WRITE (log_unit, '(A,I8,I8)') ' # nx, n = ', nx, n
759 : END IF
760 :
761 2 : ALLOCATE (x(nx))
762 2 : ALLOCATE (x2(nx))
763 2 : ALLOCATE (h(nx))
764 2 : ALLOCATE (fp1(nx))
765 10 : ALLOCATE (xk(p - 1, nx))
766 6 : ALLOCATE (xk2(p - 1, nx))
767 8 : ALLOCATE (kk(p - 1, nx))
768 6 : ALLOCATE (fpxk(p - 1, nx))
769 :
770 : ! initialize F_P(x) = f_P(x_1)
771 : ! fp1 = fp(x) = h(x/(P-1))
772 : ! fpxk = fp(xk) = h(xk/(P-1))
773 480002 : DO j = 1, nx
774 480000 : x(j) = xmin + (j - 1)*dx
775 480000 : x2(j) = x(j)**2
776 480000 : h(j) = x(j)/TANH(x(j))
777 480000 : IF (x(j) <= 1.0e-10_dp) h(j) = 1.0_dp
778 480000 : fp1(j) = op1*x(j)/TANH(x(j)*op1)
779 480000 : IF (x(j)*op1 <= 1.0e-10_dp) fp1(j) = 1.0_dp
780 1920002 : DO k = 1, p - 1
781 1440000 : xk2(k, j) = x2(j) + (p*SIN(k*pi*op))**2
782 1440000 : xk(k, j) = SQRT(xk2(k, j) - (p*SIN(pi*op))**2)
783 1440000 : kk(k, j) = NINT((xk(k, j) - xmin)/dx) + 1
784 1440000 : fpxk(k, j) = xk(k, j)*op1/TANH(xk(k, j)*op1)
785 1920000 : IF (xk(k, j)*op1 <= 1.0e-10_dp) fpxk(k, j) = 1.0_dp
786 : END DO
787 : END DO
788 :
789 : ! **** resolution ****
790 : ! compute fp(x)
791 82 : DO i = 1, niter
792 80 : err = 0.0_dp
793 19200080 : DO j = 1, nx
794 : tmp = 0.0_dp
795 57600000 : DO k = 2, p - 1
796 57600000 : tmp = tmp + fpxk(k, j)*x2(j)/xk2(k, j)
797 : END DO
798 19200000 : fprev = fp1(j)
799 19200000 : tmp1 = 1.0_dp + (p*SIN(pi*op)/x(j))**2
800 19200000 : fp1(j) = malpha*tmp1*(h(j) - 1.0_dp - tmp) + (1.0_dp - malpha)*fp1(j)
801 19200080 : IF (j <= n) err = err + ABS(1.0_dp - fp1(j)/fprev) ! compute "errors"
802 : END DO
803 80 : err = err/n
804 :
805 : ! Linear regression on the last 20% of the F_P function
806 80 : CALL pint_qtb_linreg(fp1(8*nx/10:nx), x(8*nx/10:nx), aa, bb, r2, log_unit, print_level)
807 :
808 : ! compute the new F_P(xk*sqrt(P))
809 : ! through linear interpolation
810 : ! or linear extrapolation if outside of the range
811 19200082 : DO j = 1, nx
812 76800080 : DO k = 1, p - 1
813 76800000 : IF (kk(k, j) < nx) THEN
814 : fpxk(k, j) = fp1(kk(k, j)) + (fp1(kk(k, j) + 1) - fp1(kk(k, j)))/dx* &
815 57599760 : (xk(k, j) - x(kk(k, j)))
816 : ELSE
817 240 : fpxk(k, j) = aa*xk(k, j) + bb
818 : END IF
819 : END DO
820 : END DO
821 : END DO
822 :
823 2 : IF (print_level == debug_print_level) THEN
824 : ! **** tests ****
825 0 : WRITE (log_unit, '(A,ES9.3)') ' # average error during computation: ', err
826 0 : WRITE (log_unit, '(A,ES9.3)') ' # slope of F_P at high freq. - theoretical: ', op1
827 0 : WRITE (log_unit, '(A,ES9.3)') ' # slope of F_P at high freq. - calculated: ', aa
828 2 : ELSE IF (print_level > silent_print_level) THEN
829 2 : CALL pint_write_line("QTB| Initialization of random forces using fP1 function")
830 2 : CALL pint_write_line("QTB| Computation of fP1 function")
831 2 : WRITE (line, '(A,ES9.3)') 'QTB| average error ', err
832 2 : CALL pint_write_line(TRIM(line))
833 2 : WRITE (line, '(A,ES9.3)') 'QTB| slope at high frequency - theoretical: ', op1
834 2 : CALL pint_write_line(TRIM(line))
835 2 : WRITE (line, '(A,ES9.3)') 'QTB| slope at high frequency - calculated: ', aa
836 2 : CALL pint_write_line(TRIM(line))
837 : END IF
838 :
839 2 : IF (print_level == debug_print_level) THEN
840 : ! **** write solution ****
841 0 : WRITE (log_unit, '(A)') ' # ------------------------------------------------'
842 0 : WRITE (log_unit, '(A)') ' # computed fp function'
843 0 : WRITE (log_unit, '(A)') ' # i, w(a.u.), x, fp'
844 0 : DO j = 1, nx
845 0 : WRITE (log_unit, *) j, j*dw, xmin + (j - 1)*dx, fp1(j)
846 : END DO
847 : END IF
848 :
849 2 : DEALLOCATE (x2)
850 2 : DEALLOCATE (h)
851 2 : DEALLOCATE (xk)
852 2 : DEALLOCATE (xk2)
853 2 : DEALLOCATE (kk)
854 2 : DEALLOCATE (fpxk)
855 : END IF
856 :
857 : ! compute values of fP on the grid points for the current NM
858 : ! trough linear interpolation / regression
859 390 : DO j = 1, n
860 384 : tmp = (j*dx1)**2 - (p*SIN(pi*op))**2
861 390 : IF (tmp < 0.d0) THEN
862 72 : fp(j) = fp1(1)
863 : ELSE
864 312 : tmp = SQRT(tmp)
865 312 : k = NINT((tmp - xmin)/dx) + 1
866 312 : IF (k > nx) THEN
867 0 : fp(j) = aa*tmp + bb
868 312 : ELSE IF (k <= 0) THEN
869 0 : fp(j) = fp1(1)
870 : ELSE
871 312 : xx = xmin + (k - 1)*dx
872 312 : IF (tmp > xx) THEN
873 168 : fp(j) = fp1(k) + (fp1(k + 1) - fp1(k))/dx*(tmp - xx)
874 : ELSE
875 144 : fp(j) = fp1(k) + (fp1(k) - fp1(k - 1))/dx*(tmp - xx)
876 : END IF
877 : END IF
878 : END IF
879 : END DO
880 :
881 : END IF
882 :
883 8 : END SUBROUTINE pint_qtb_computefp1
884 :
885 : ! ***************************************************************************
886 : !> \brief perform a simple linear regression - y(x) = a*x + b
887 : !> \param y ...
888 : !> \param x ...
889 : !> \param a ...
890 : !> \param b ...
891 : !> \param r2 ...
892 : !> \param log_unit ...
893 : !> \param print_level ...
894 : !> \author Fabien Brieuc
895 : ! **************************************************************************************************
896 110 : SUBROUTINE pint_qtb_linreg(y, x, a, b, r2, log_unit, print_level)
897 : REAL(KIND=dp), DIMENSION(:) :: y, x
898 : REAL(KIND=dp) :: a, b, r2
899 : INTEGER :: log_unit, print_level
900 :
901 : CHARACTER(len=200) :: line
902 : INTEGER :: i, n
903 : REAL(KIND=dp) :: xav, xvar, xycov, yav, yvar
904 :
905 110 : n = SIZE(y)
906 :
907 110 : xav = 0.0_dp
908 110 : yav = 0.0_dp
909 110 : xycov = 0.0_dp
910 110 : xvar = 0.0_dp
911 110 : yvar = 0.0_dp
912 :
913 5280220 : DO i = 1, n
914 5280110 : xav = xav + x(i)
915 5280110 : yav = yav + y(i)
916 5280110 : xycov = xycov + x(i)*y(i)
917 5280110 : xvar = xvar + x(i)**2
918 5280220 : yvar = yvar + y(i)**2
919 : END DO
920 :
921 110 : xav = xav/n
922 110 : yav = yav/n
923 110 : xycov = xycov/n
924 110 : xycov = xycov - xav*yav
925 110 : xvar = xvar/n
926 110 : xvar = xvar - xav**2
927 110 : yvar = yvar/n
928 110 : yvar = yvar - yav**2
929 :
930 110 : a = xycov/xvar
931 110 : b = yav - a*xav
932 :
933 110 : r2 = xycov/SQRT(xvar*yvar)
934 :
935 110 : IF (r2 < 0.9_dp) THEN
936 0 : IF (print_level == debug_print_level) THEN
937 0 : WRITE (log_unit, '(A, E10.3)') '# possible error during linear regression: r^2 = ', r2
938 0 : ELSE IF (print_level > silent_print_level) THEN
939 0 : WRITE (line, '(A,E10.3)') 'QTB| possible error during linear regression: r^2 = ', r2
940 0 : CALL pint_write_line(TRIM(line))
941 : END IF
942 : END IF
943 :
944 110 : END SUBROUTINE pint_qtb_linreg
945 :
946 : ! **************************************************************************************************
947 : !> \brief ...
948 : !> \param z_in ...
949 : !> \param z_out ...
950 : !> \param n ...
951 : ! **************************************************************************************************
952 12 : SUBROUTINE pint_qtb_fft(z_in, z_out, n)
953 :
954 : INTEGER :: n
955 : COMPLEX(KIND=dp), DIMENSION(n) :: z_out, z_in
956 :
957 : INTEGER :: stat
958 :
959 12 : CALL fft_1d_many(FWFFT, n, 1, .FALSE., .FALSE., n, n, z_in, z_out, 1.0_dp, stat)
960 12 : END SUBROUTINE pint_qtb_fft
961 :
962 : END MODULE pint_qtb
|
