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