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