LCOV - code coverage report
Current view: top level - src/motion - pint_qtb.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 364 456 79.8 %
Date: 2024-12-12 07:20:04 Functions: 10 10 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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 1.15