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 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              :    ! Abort and inform that GreenX was not included in the compilation
      16              :    ! Ideally, this will be avoided in testing by the conditional tests
      17              :    CPABORT("CP2K not compiled with GreenX library.")
      18              : #else
      19              :    USE kinds, ONLY: dp
      20              :    USE gx_ac, ONLY: create_thiele_pade, &
      21              :                     evaluate_thiele_pade_at, &
      22              :                     free_params, &
      23              :                     params
      24              : 
      25              :    ! Create the dataset containing the fitting data
      26              :    ! Two Lorentzian peaks with some overlap
      27              :    COMPLEX(kind=dp)                   :: damp_one = (2, 0), &
      28              :                                          damp_two = (4, 0), &
      29              :                                          center_one = (2, 0), &
      30              :                                          center_two = (8, 0), &
      31              :                                          amp_one = (10, 0), &
      32              :                                          amp_two = (10, 0), &
      33              :                                          min_source = (0, 0), &
      34              :                                          min_fit = (0, 0), &
      35              :                                          max_source = (10, 0), &
      36              :                                          max_fit = (10, 0)
      37              :    INTEGER, PARAMETER                 :: n_source = 20, &
      38              :                                          n_fit = 100, &
      39              :                                          n_param = 10
      40              :    INTEGER                            :: i
      41              :    COMPLEX(kind=dp)                   :: d_source, &
      42              :                                          d_fit
      43              :    COMPLEX(kind=dp), DIMENSION(n_source) :: x_source, &
      44              :                                             y_source
      45              :    COMPLEX(kind=dp), DIMENSION(n_fit)    :: x_fit, &
      46              :                                             y_fit
      47            2 :    TYPE(params)                       :: fit_params
      48              : 
      49            2 :    d_source = (max_source - min_source)/CMPLX(n_source - 1, kind=dp)
      50            2 :    d_fit = (max_fit - min_fit)/CMPLX(n_fit - 1, kind=dp)
      51              : 
      52            2 :    PRINT '(A12)', "#Source data"
      53              : 
      54           42 :    DO i = 1, n_source
      55           40 :       x_source(i) = min_source + CMPLX(i - 1, 0.0, kind=dp)*d_source
      56           40 :       y_source(i) = amp_one/(damp_one*damp_one + (x_source(i) - center_one)*(x_source(i) - center_one))
      57           40 :       y_source(i) = y_source(i) + amp_two/(damp_two*damp_two + (x_source(i) - center_two)*(x_source(i) - center_two))
      58           42 :       PRINT '(E20.8E3,E20.8E3)', REAL(x_source(i), kind=dp), REAL(y_source(i), kind=dp)
      59              :    END DO
      60              : 
      61            2 :    PRINT '(A9)', "#Fit data"
      62              : 
      63              :    ! Fit points created
      64              :    ! Now do the actual fitting
      65            2 :    fit_params = create_thiele_pade(n_param, x_source, y_source)
      66              : 
      67              :    ! Create the evaluation grid
      68          202 :    DO i = 1, n_fit
      69          202 :       x_fit(i) = min_fit + d_fit*CMPLX(i - 1, 0, kind=dp)
      70              :    END DO
      71              : 
      72            2 :    y_fit(1:n_fit) = evaluate_thiele_pade_at(fit_params, x_fit)
      73              : 
      74          202 :    DO i = 1, n_fit
      75          202 :       PRINT '(E20.8E3,E20.8E3)', REAL(x_fit(i), kind=dp), REAL(y_fit(i), kind=dp)
      76              :    END DO
      77              : 
      78            2 :    CALL free_params(fit_params)
      79              : #endif
      80            2 : END PROGRAM gx_ac_unittest