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