LCOV - code coverage report
Current view: top level - src/motion - pint_qtb.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 79.8 % 456 364
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 10 10

            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
        

Generated by: LCOV version 2.0-1