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 link - please recompile to access GXAC.")
      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