Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief GreenX Analytic continuation unit test
      10              : !> \author Stepan Marek
      11              : ! **************************************************************************************************
      12            2 : PROGRAM gx_ac_unittest
      13              : #include "base/base_uses.f90"
      14              : #if defined(__GREENX)
      15              :    USE kinds, ONLY: dp
      16              :    USE gx_ac, ONLY: create_thiele_pade, &
      17              :                     evaluate_thiele_pade_at, &
      18              :                     free_params, &
      19              :                     params
      20              : #endif
      21              : 
      22              :    IMPLICIT NONE
      23              : 
      24              : #if !defined(__GREENX)
      25              :    ! Abort and inform that GreenX was not included in the compilation
      26              :    ! Ideally, this will be avoided in testing by the conditional tests
      27              :    CPABORT("CP2K not compiled with GreenX library.")
      28              : #else
      29              :    ! Create the dataset containing the fitting data
      30              :    ! Two Lorentzian peaks with some overlap
      31              :    COMPLEX(kind=dp)                   :: damp_one = (2, 0), &
      32              :                                          damp_two = (4, 0), &
      33              :                                          center_one = (2, 0), &
      34              :                                          center_two = (8, 0), &
      35              :                                          amp_one = (10, 0), &
      36              :                                          amp_two = (10, 0), &
      37              :                                          min_source = (0, 0), &
      38              :                                          min_fit = (0, 0), &
      39              :                                          max_source = (10, 0), &
      40              :                                          max_fit = (10, 0)
      41              :    INTEGER, PARAMETER                 :: n_source = 20, &
      42              :                                          n_fit = 100, &
      43              :                                          n_param = 10
      44              :    INTEGER                            :: i
      45              :    COMPLEX(kind=dp)                   :: d_source, &
      46              :                                          d_fit
      47              :    COMPLEX(kind=dp), DIMENSION(n_source) :: x_source, &
      48              :                                             y_source
      49              :    COMPLEX(kind=dp), DIMENSION(n_fit)    :: x_fit, &
      50              :                                             y_fit
      51            2 :    TYPE(params)                       :: fit_params
      52              : 
      53            2 :    d_source = (max_source - min_source)/CMPLX(n_source - 1, kind=dp)
      54            2 :    d_fit = (max_fit - min_fit)/CMPLX(n_fit - 1, kind=dp)
      55              : 
      56            2 :    PRINT '(A12)', "#Source data"
      57              : 
      58           42 :    DO i = 1, n_source
      59           40 :       x_source(i) = min_source + CMPLX(i - 1, 0.0, kind=dp)*d_source
      60           40 :       y_source(i) = amp_one/(damp_one*damp_one + (x_source(i) - center_one)*(x_source(i) - center_one))
      61           40 :       y_source(i) = y_source(i) + amp_two/(damp_two*damp_two + (x_source(i) - center_two)*(x_source(i) - center_two))
      62           42 :       PRINT '(E20.8E3,E20.8E3)', REAL(x_source(i), kind=dp), REAL(y_source(i), kind=dp)
      63              :    END DO
      64              : 
      65            2 :    PRINT '(A9)', "#Fit data"
      66              : 
      67              :    ! Fit points created
      68              :    ! Now do the actual fitting
      69            2 :    fit_params = create_thiele_pade(n_param, x_source, y_source)
      70              : 
      71              :    ! Create the evaluation grid
      72          202 :    DO i = 1, n_fit
      73          202 :       x_fit(i) = min_fit + d_fit*CMPLX(i - 1, 0, kind=dp)
      74              :    END DO
      75              : 
      76            2 :    y_fit(1:n_fit) = evaluate_thiele_pade_at(fit_params, x_fit)
      77              : 
      78          202 :    DO i = 1, n_fit
      79          202 :       PRINT '(E20.8E3,E20.8E3)', REAL(x_fit(i), kind=dp), REAL(y_fit(i), kind=dp)
      80              :    END DO
      81              : 
      82            2 :    CALL free_params(fit_params)
      83              : #endif
      84            2 : END PROGRAM gx_ac_unittest