LCOV - code coverage report
Current view: top level - src - mp2_grids.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.1 % 542 499
Test Date: 2025-07-25 12:55:17 Functions: 94.1 % 17 16

            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 Routines to calculate frequency and time grids (integration points and weights)
      10              : !>        for correlation methods
      11              : !> \par History
      12              : !>      05.2019 Refactored from rpa_ri_gpw [Frederick Stein]
      13              : ! **************************************************************************************************
      14              : MODULE mp2_grids
      15              :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      16              :                                               cp_fm_type
      17              :    USE greenx_interface,                ONLY: greenx_get_minimax_grid
      18              :    USE input_section_types,             ONLY: section_vals_type,&
      19              :                                               section_vals_val_set
      20              :    USE kinds,                           ONLY: dp
      21              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      22              :                                               kpoint_env_type,&
      23              :                                               kpoint_type
      24              :    USE machine,                         ONLY: m_flush
      25              :    USE mathconstants,                   ONLY: pi
      26              :    USE message_passing,                 ONLY: mp_para_env_release,&
      27              :                                               mp_para_env_type
      28              :    USE minimax_exp,                     ONLY: get_exp_minimax_coeff
      29              :    USE minimax_exp_gw,                  ONLY: get_exp_minimax_coeff_gw
      30              :    USE minimax_rpa,                     ONLY: get_rpa_minimax_coeff,&
      31              :                                               get_rpa_minimax_coeff_larger_grid
      32              :    USE qs_environment_types,            ONLY: get_qs_env,&
      33              :                                               qs_environment_type
      34              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      35              :                                               mo_set_type
      36              : #include "./base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              : 
      40              :    PRIVATE
      41              : 
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_grids'
      43              : 
      44              :    PUBLIC :: get_minimax_grid, get_clenshaw_grid, test_least_square_ft, get_l_sq_wghts_cos_tf_t_to_w, &
      45              :              get_l_sq_wghts_cos_tf_w_to_t, get_l_sq_wghts_sin_tf_t_to_w
      46              : 
      47              : CONTAINS
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief ...
      51              : !> \param para_env ...
      52              : !> \param unit_nr ...
      53              : !> \param homo ...
      54              : !> \param Eigenval ...
      55              : !> \param num_integ_points ...
      56              : !> \param do_im_time ...
      57              : !> \param do_ri_sos_laplace_mp2 ...
      58              : !> \param do_print ...
      59              : !> \param tau_tj ...
      60              : !> \param tau_wj ...
      61              : !> \param qs_env ...
      62              : !> \param do_gw_im_time ...
      63              : !> \param do_kpoints_cubic_RPA ...
      64              : !> \param e_fermi ...
      65              : !> \param tj ...
      66              : !> \param wj ...
      67              : !> \param weights_cos_tf_t_to_w ...
      68              : !> \param weights_cos_tf_w_to_t ...
      69              : !> \param weights_sin_tf_t_to_w ...
      70              : !> \param regularization ...
      71              : ! **************************************************************************************************
      72          204 :    SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, &
      73              :                                do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, &
      74              :                                do_kpoints_cubic_RPA, e_fermi, tj, wj, weights_cos_tf_t_to_w, &
      75              :                                weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
      76              : 
      77              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
      78              :       INTEGER, INTENT(IN)                                :: unit_nr
      79              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
      80              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Eigenval
      81              :       INTEGER, INTENT(IN)                                :: num_integ_points
      82              :       LOGICAL, INTENT(IN)                                :: do_im_time, do_ri_sos_laplace_mp2, &
      83              :                                                             do_print
      84              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
      85              :          INTENT(OUT)                                     :: tau_tj, tau_wj
      86              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      87              :       LOGICAL, INTENT(IN)                                :: do_gw_im_time, do_kpoints_cubic_RPA
      88              :       REAL(KIND=dp), INTENT(OUT)                         :: e_fermi
      89              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
      90              :          INTENT(OUT)                                     :: tj, wj
      91              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
      92              :          INTENT(OUT)                                     :: weights_cos_tf_t_to_w, &
      93              :                                                             weights_cos_tf_w_to_t, &
      94              :                                                             weights_sin_tf_t_to_w
      95              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: regularization
      96              : 
      97              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_minimax_grid'
      98              :       INTEGER, PARAMETER                                 :: num_points_per_magnitude = 200
      99              : 
     100              :       INTEGER                                            :: handle, ierr, jquad, nspins
     101              :       LOGICAL                                            :: my_do_kpoints
     102              :       REAL(KIND=dp)                                      :: E_Range, Emax, Emin, max_error_min, &
     103              :                                                             my_regularization, scaling
     104          204 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: x_tw
     105              : 
     106          204 :       CALL timeset(routineN, handle)
     107              : 
     108              :       CALL determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
     109          204 :                                   do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
     110              : 
     111              :       CALL greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, &
     112              :                                    tau_tj, tau_wj, qs_env%mp2_env%ri_g0w0%regularization_minimax, &
     113              :                                    tj, wj, weights_cos_tf_t_to_w, &
     114          204 :                                    weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
     115              : 
     116              : ! Shortcut if Greenx was available and successful
     117          204 :       IF (ierr == 0) THEN
     118           74 :          CALL timestop(handle)
     119              :          RETURN
     120              :       END IF
     121              : 
     122              :       ! Test for spin unrestricted
     123          130 :       nspins = SIZE(homo)
     124              : 
     125              :       ! Test whether all necessary variables are available
     126          130 :       my_do_kpoints = .FALSE.
     127          130 :       IF (.NOT. do_ri_sos_laplace_mp2) THEN
     128          130 :          my_do_kpoints = do_kpoints_cubic_RPA
     129              :       END IF
     130              : 
     131              :       my_regularization = 0.0_dp
     132          130 :       IF (PRESENT(regularization)) THEN
     133          130 :          my_regularization = regularization
     134              : 
     135          130 :          IF (num_integ_points > 20 .AND. e_range < 100.0_dp) THEN
     136            0 :             IF (unit_nr > 0) &
     137              :                CALL cp_warn(__LOCATION__, &
     138              :                             "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// &
     139              :                             "That may lead to numerical "// &
     140              :                             "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// &
     141            0 :                             "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
     142              :          END IF
     143              : 
     144          130 :          IF (.NOT. do_ri_sos_laplace_mp2) THEN
     145          216 :             ALLOCATE (x_tw(2*num_integ_points))
     146          620 :             x_tw = 0.0_dp
     147           72 :             ierr = 0
     148           72 :             IF (num_integ_points .LE. 20) THEN
     149           72 :                CALL get_rpa_minimax_coeff(num_integ_points, e_range, x_tw, ierr)
     150              :             ELSE
     151            0 :                CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, e_range, x_tw)
     152              :             END IF
     153              : 
     154          216 :             ALLOCATE (tj(num_integ_points))
     155          346 :             tj = 0.0_dp
     156              : 
     157          144 :             ALLOCATE (wj(num_integ_points))
     158          346 :             wj = 0.0_dp
     159              : 
     160          346 :             DO jquad = 1, num_integ_points
     161          274 :                tj(jquad) = x_tw(jquad)
     162          346 :                wj(jquad) = x_tw(jquad + num_integ_points)
     163              :             END DO
     164              : 
     165              :             ! for the smaller grids, the factor of 4 is included in get_rpa_minimax_coeff for wj
     166           72 :             IF (num_integ_points >= 26) THEN
     167            0 :                wj(:) = wj(:)*4.0_dp
     168              :             END IF
     169              : 
     170           72 :             DEALLOCATE (x_tw)
     171              : 
     172           72 :             IF (unit_nr > 0 .AND. do_print) THEN
     173              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     174           35 :                   "MINIMAX_INFO| Number of integration points:", num_integ_points
     175              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
     176           35 :                   "MINIMAX_INFO| Gap for the minimax approximation:", Emin
     177              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
     178           35 :                   "MINIMAX_INFO| Range for the minimax approximation:", e_range
     179           35 :                WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
     180          167 :                DO jquad = 1, num_integ_points
     181          167 :                   WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
     182              :                END DO
     183           35 :                CALL m_flush(unit_nr)
     184              :             END IF
     185              : 
     186              :             ! scale the minimax parameters
     187          346 :             tj(:) = tj(:)*Emin
     188          346 :             wj(:) = wj(:)*Emin
     189              :          END IF
     190              : 
     191              :          ! set up the minimax time grid
     192          130 :          IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
     193              : 
     194          300 :             ALLOCATE (x_tw(2*num_integ_points))
     195          836 :             x_tw = 0.0_dp
     196              : 
     197          100 :             IF (num_integ_points .LE. 20) THEN
     198          100 :                CALL get_exp_minimax_coeff(num_integ_points, e_range, x_tw)
     199              :             ELSE
     200            0 :                CALL get_exp_minimax_coeff_gw(num_integ_points, e_range, x_tw)
     201              :             END IF
     202              : 
     203              :             ! For RPA we include already a factor of two (see later steps)
     204          100 :             scaling = 2.0_dp
     205          100 :             IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp
     206              : 
     207          300 :             ALLOCATE (tau_tj(num_integ_points))
     208          468 :             tau_tj = 0.0_dp
     209              : 
     210          200 :             ALLOCATE (tau_wj(num_integ_points))
     211          468 :             tau_wj = 0.0_dp
     212              : 
     213          468 :             DO jquad = 1, num_integ_points
     214          368 :                tau_tj(jquad) = x_tw(jquad)/scaling
     215          468 :                tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
     216              :             END DO
     217              : 
     218          100 :             DEALLOCATE (x_tw)
     219              : 
     220          100 :             IF (unit_nr > 0 .AND. do_print) THEN
     221              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
     222           49 :                   "MINIMAX_INFO| Range for the minimax approximation:", e_range
     223              :                ! For testing the gap
     224              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
     225           49 :                   "MINIMAX_INFO| Gap:", Emin
     226              :                WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
     227           49 :                   "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
     228          228 :                DO jquad = 1, num_integ_points
     229          228 :                   WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
     230              :                END DO
     231           49 :                CALL m_flush(unit_nr)
     232              :             END IF
     233              : 
     234              :             ! scale grid from [1,R] to [Emin,Emax]
     235          468 :             tau_tj(:) = tau_tj(:)/Emin
     236          468 :             tau_wj(:) = tau_wj(:)/Emin
     237              : 
     238          100 :             IF (.NOT. do_ri_sos_laplace_mp2) THEN
     239          168 :                ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
     240          678 :                weights_cos_tf_t_to_w = 0.0_dp
     241              : 
     242              :                CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, &
     243              :                                                  Emin, Emax, max_error_min, num_points_per_magnitude, &
     244           42 :                                                  my_regularization)
     245              : 
     246              :                ! get the weights for the cosine transform W^c(iw) -> W^c(it)
     247          126 :                ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
     248          678 :                weights_cos_tf_w_to_t = 0.0_dp
     249              : 
     250              :                CALL get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, tj, &
     251              :                                                  Emin, Emax, max_error_min, num_points_per_magnitude, &
     252           42 :                                                  my_regularization)
     253              : 
     254           42 :                IF (do_gw_im_time) THEN
     255              : 
     256              :                   ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71)
     257           30 :                   ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points))
     258          310 :                   weights_sin_tf_t_to_w = 0.0_dp
     259              : 
     260              :                   CALL get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, tj, &
     261              :                                                     Emin, Emax, max_error_min, num_points_per_magnitude, &
     262           10 :                                                     my_regularization)
     263              : 
     264           10 :                   IF (unit_nr > 0) THEN
     265              :                      WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
     266            5 :                         "MINIMAX_INFO| Maximum deviation of the imag. time fit:", max_error_min
     267              :                   END IF
     268              :                END IF
     269              : 
     270              :             END IF
     271              : 
     272              :          END IF
     273              :       END IF
     274              : 
     275          130 :       CALL timestop(handle)
     276              : 
     277          130 :    END SUBROUTINE get_minimax_grid
     278              : 
     279              : ! **************************************************************************************************
     280              : !> \brief ...
     281              : !> \param para_env ...
     282              : !> \param para_env_RPA ...
     283              : !> \param unit_nr ...
     284              : !> \param homo ...
     285              : !> \param virtual ...
     286              : !> \param Eigenval ...
     287              : !> \param num_integ_points ...
     288              : !> \param num_integ_group ...
     289              : !> \param color_rpa_group ...
     290              : !> \param fm_mat_S ...
     291              : !> \param my_do_gw ...
     292              : !> \param ext_scaling ...
     293              : !> \param a_scaling ...
     294              : !> \param tj ...
     295              : !> \param wj ...
     296              : ! **************************************************************************************************
     297           96 :    SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
     298           96 :                                 num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, &
     299              :                                 ext_scaling, a_scaling, tj, wj)
     300              : 
     301              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env, para_env_RPA
     302              :       INTEGER, INTENT(IN)                                :: unit_nr
     303              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo, virtual
     304              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Eigenval
     305              :       INTEGER, INTENT(IN)                                :: num_integ_points, num_integ_group, &
     306              :                                                             color_rpa_group
     307              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_mat_S
     308              :       LOGICAL, INTENT(IN)                                :: my_do_gw
     309              :       REAL(KIND=dp), INTENT(IN)                          :: ext_scaling
     310              :       REAL(KIND=dp), INTENT(OUT)                         :: a_scaling
     311              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     312              :          INTENT(OUT)                                     :: tj, wj
     313              : 
     314              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_clenshaw_grid'
     315              : 
     316              :       INTEGER                                            :: handle, jquad, nspins
     317              :       LOGICAL                                            :: my_open_shell
     318              : 
     319           96 :       CALL timeset(routineN, handle)
     320              : 
     321           96 :       nspins = SIZE(homo)
     322           96 :       my_open_shell = (nspins == 2)
     323              : 
     324              :       ! Now, start to prepare the different grid
     325          288 :       ALLOCATE (tj(num_integ_points))
     326         3766 :       tj = 0.0_dp
     327              : 
     328          192 :       ALLOCATE (wj(num_integ_points))
     329         3766 :       wj = 0.0_dp
     330              : 
     331         3670 :       DO jquad = 1, num_integ_points - 1
     332         3574 :          tj(jquad) = jquad*pi/(2.0_dp*num_integ_points)
     333         3670 :          wj(jquad) = pi/(num_integ_points*SIN(tj(jquad))**2)
     334              :       END DO
     335           96 :       tj(num_integ_points) = pi/2.0_dp
     336           96 :       wj(num_integ_points) = pi/(2.0_dp*num_integ_points*SIN(tj(num_integ_points))**2)
     337              : 
     338           96 :       IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
     339           58 :          a_scaling = ext_scaling
     340              :       ELSE
     341              :          CALL calc_scaling_factor(a_scaling, para_env, para_env_RPA, homo, virtual, Eigenval, &
     342              :                                   num_integ_points, num_integ_group, color_rpa_group, &
     343           38 :                                   tj, wj, fm_mat_S)
     344              :       END IF
     345              : 
     346           96 :       IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
     347              : 
     348         3766 :       wj(:) = wj(:)*a_scaling
     349              : 
     350           96 :       CALL timestop(handle)
     351              : 
     352           96 :    END SUBROUTINE get_clenshaw_grid
     353              : 
     354              : ! **************************************************************************************************
     355              : !> \brief ...
     356              : !> \param a_scaling_ext ...
     357              : !> \param para_env ...
     358              : !> \param para_env_RPA ...
     359              : !> \param homo ...
     360              : !> \param virtual ...
     361              : !> \param Eigenval ...
     362              : !> \param num_integ_points ...
     363              : !> \param num_integ_group ...
     364              : !> \param color_rpa_group ...
     365              : !> \param tj_ext ...
     366              : !> \param wj_ext ...
     367              : !> \param fm_mat_S ...
     368              : ! **************************************************************************************************
     369           38 :    SUBROUTINE calc_scaling_factor(a_scaling_ext, para_env, para_env_RPA, homo, virtual, Eigenval, &
     370              :                                   num_integ_points, num_integ_group, color_rpa_group, &
     371           38 :                                   tj_ext, wj_ext, fm_mat_S)
     372              :       REAL(KIND=dp), INTENT(OUT)                         :: a_scaling_ext
     373              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env, para_env_RPA
     374              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo, virtual
     375              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Eigenval
     376              :       INTEGER, INTENT(IN)                                :: num_integ_points, num_integ_group, &
     377              :                                                             color_rpa_group
     378              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     379              :          INTENT(IN)                                      :: tj_ext, wj_ext
     380              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_mat_S
     381              : 
     382              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_scaling_factor'
     383              : 
     384              :       INTEGER                                            :: handle, icycle, jquad, ncol_local, &
     385              :                                                             ncol_local_beta, nspins
     386              :       LOGICAL                                            :: my_open_shell
     387              :       REAL(KIND=dp) :: a_high, a_low, a_scaling, conv_param, eps, first_deriv, left_term, &
     388              :          right_term, right_term_ref, right_term_ref_beta, step
     389           38 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cottj, D_ia, D_ia_beta, iaia_RI, &
     390           38 :                                                             iaia_RI_beta, M_ia, M_ia_beta
     391              :       TYPE(mp_para_env_type), POINTER                    :: para_env_col, para_env_col_beta
     392              : 
     393           38 :       CALL timeset(routineN, handle)
     394              : 
     395           38 :       nspins = SIZE(homo)
     396           38 :       my_open_shell = (nspins == 2)
     397              : 
     398           38 :       eps = 1.0E-10_dp
     399              : 
     400          114 :       ALLOCATE (cottj(num_integ_points))
     401              : 
     402              :       ! calculate the cotangent of the abscissa tj
     403          488 :       DO jquad = 1, num_integ_points
     404          488 :          cottj(jquad) = 1.0_dp/TAN(tj_ext(jquad))
     405              :       END DO
     406              : 
     407              :       CALL calc_ia_ia_integrals(para_env_RPA, homo(1), virtual(1), ncol_local, right_term_ref, Eigenval(:, 1, 1), &
     408           38 :                                 D_ia, iaia_RI, M_ia, fm_mat_S(1), para_env_col)
     409              : 
     410              :       ! In the open shell case do point 1-2-3 for the beta spin
     411           38 :       IF (my_open_shell) THEN
     412              :          CALL calc_ia_ia_integrals(para_env_RPA, homo(2), virtual(2), ncol_local_beta, right_term_ref_beta, Eigenval(:, 1, 2), &
     413            8 :                                    D_ia_beta, iaia_RI_beta, M_ia_beta, fm_mat_S(2), para_env_col_beta)
     414              : 
     415            8 :          right_term_ref = right_term_ref + right_term_ref_beta
     416              :       END IF
     417              : 
     418              :       ! bcast the result
     419           38 :       IF (para_env%mepos == 0) THEN
     420           19 :          CALL para_env%bcast(right_term_ref, 0)
     421              :       ELSE
     422           19 :          right_term_ref = 0.0_dp
     423           19 :          CALL para_env%bcast(right_term_ref, 0)
     424              :       END IF
     425              : 
     426              :       ! 5) start iteration for solving the non-linear equation by bisection
     427              :       ! find limit, here step=0.5 seems a good compromise
     428           38 :       conv_param = 100.0_dp*EPSILON(right_term_ref)
     429           38 :       step = 0.5_dp
     430           38 :       a_low = 0.0_dp
     431           38 :       a_high = step
     432           38 :       right_term = -right_term_ref
     433          100 :       DO icycle = 1, num_integ_points*2
     434           92 :          a_scaling = a_high
     435              : 
     436              :          CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
     437              :                                 M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
     438              :                                 ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
     439           92 :                                 para_env, para_env_col, para_env_col_beta)
     440           92 :          left_term = left_term/4.0_dp/pi*a_scaling
     441              : 
     442           92 :          IF (ABS(left_term) > ABS(right_term) .OR. ABS(left_term + right_term) <= conv_param) EXIT
     443           62 :          a_low = a_high
     444          100 :          a_high = a_high + step
     445              : 
     446              :       END DO
     447              : 
     448           38 :       IF (ABS(left_term + right_term) >= conv_param) THEN
     449           32 :          IF (a_scaling >= 2*num_integ_points*step) THEN
     450           10 :             a_scaling = 1.0_dp
     451              :          ELSE
     452              : 
     453          340 :             DO icycle = 1, num_integ_points*2
     454          336 :                a_scaling = (a_low + a_high)/2.0_dp
     455              : 
     456              :                CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
     457              :                                       M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
     458              :                                       ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
     459          336 :                                       para_env, para_env_col, para_env_col_beta)
     460          336 :                left_term = left_term/4.0_dp/pi*a_scaling
     461              : 
     462          336 :                IF (ABS(left_term) > ABS(right_term)) THEN
     463              :                   a_high = a_scaling
     464              :                ELSE
     465          186 :                   a_low = a_scaling
     466              :                END IF
     467              : 
     468          340 :                IF (ABS(a_high - a_low) < 1.0e-5_dp) EXIT
     469              : 
     470              :             END DO
     471              : 
     472              :          END IF
     473              :       END IF
     474              : 
     475           38 :       a_scaling_ext = a_scaling
     476           38 :       CALL para_env%bcast(a_scaling_ext, 0)
     477              : 
     478           38 :       DEALLOCATE (cottj)
     479           38 :       DEALLOCATE (iaia_RI)
     480           38 :       DEALLOCATE (D_ia)
     481           38 :       DEALLOCATE (M_ia)
     482           38 :       CALL mp_para_env_release(para_env_col)
     483              : 
     484           38 :       IF (my_open_shell) THEN
     485            8 :          DEALLOCATE (iaia_RI_beta)
     486            8 :          DEALLOCATE (D_ia_beta)
     487            8 :          DEALLOCATE (M_ia_beta)
     488            8 :          CALL mp_para_env_release(para_env_col_beta)
     489              :       END IF
     490              : 
     491           38 :       CALL timestop(handle)
     492              : 
     493           76 :    END SUBROUTINE calc_scaling_factor
     494              : 
     495              : ! **************************************************************************************************
     496              : !> \brief ...
     497              : !> \param para_env_RPA ...
     498              : !> \param homo ...
     499              : !> \param virtual ...
     500              : !> \param ncol_local ...
     501              : !> \param right_term_ref ...
     502              : !> \param Eigenval ...
     503              : !> \param D_ia ...
     504              : !> \param iaia_RI ...
     505              : !> \param M_ia ...
     506              : !> \param fm_mat_S ...
     507              : !> \param para_env_col ...
     508              : ! **************************************************************************************************
     509           46 :    SUBROUTINE calc_ia_ia_integrals(para_env_RPA, homo, virtual, ncol_local, right_term_ref, Eigenval, &
     510              :                                    D_ia, iaia_RI, M_ia, fm_mat_S, para_env_col)
     511              : 
     512              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env_RPA
     513              :       INTEGER, INTENT(IN)                                :: homo, virtual
     514              :       INTEGER, INTENT(OUT)                               :: ncol_local
     515              :       REAL(KIND=dp), INTENT(OUT)                         :: right_term_ref
     516              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
     517              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     518              :          INTENT(OUT)                                     :: D_ia, iaia_RI, M_ia
     519              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S
     520              :       TYPE(mp_para_env_type), POINTER                    :: para_env_col
     521              : 
     522              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ia_ia_integrals'
     523              : 
     524              :       INTEGER                                            :: avirt, color_col, color_row, handle, &
     525              :                                                             i_global, iiB, iocc, nrow_local
     526           46 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     527              :       REAL(KIND=dp)                                      :: eigen_diff
     528              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: iaia_RI_dp
     529              :       TYPE(mp_para_env_type), POINTER                    :: para_env_row
     530              : 
     531           46 :       CALL timeset(routineN, handle)
     532              : 
     533              :       ! calculate the (ia|ia) RI integrals
     534              :       ! ----------------------------------
     535              :       ! 1) get info fm_mat_S
     536              :       CALL cp_fm_get_info(matrix=fm_mat_S, &
     537              :                           nrow_local=nrow_local, &
     538              :                           ncol_local=ncol_local, &
     539              :                           row_indices=row_indices, &
     540           46 :                           col_indices=col_indices)
     541              : 
     542              :       ! allocate the local buffer of iaia_RI integrals (dp kind)
     543          136 :       ALLOCATE (iaia_RI_dp(ncol_local))
     544         2764 :       iaia_RI_dp = 0.0_dp
     545              : 
     546              :       ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
     547         2764 :       DO iiB = 1, ncol_local
     548       188006 :          iaia_RI_dp(iiB) = iaia_RI_dp(iiB) + DOT_PRODUCT(fm_mat_S%local_data(:, iiB), fm_mat_S%local_data(:, iiB))
     549              :       END DO
     550              : 
     551              :       ! 3) sum the result with the processes of the RPA_group having the same columns
     552              :       !          _______ia______               _
     553              :       !         |   |   |   |   |             | |
     554              :       !     --> | 1 | 5 | 9 | 13|   SUM -->   | |
     555              :       !         |___|__ |___|___|             |_|
     556              :       !         |   |   |   |   |             | |
     557              :       !     --> | 2 | 6 | 10| 14|   SUM -->   | |
     558              :       !       K |___|___|___|___|             |_|   (ia|ia)_RI
     559              :       !         |   |   |   |   |             | |
     560              :       !     --> | 3 | 7 | 11| 15|   SUM -->   | |
     561              :       !         |___|___|___|___|             |_|
     562              :       !         |   |   |   |   |             | |
     563              :       !     --> | 4 | 8 | 12| 16|   SUM -->   | |
     564              :       !         |___|___|___|___|             |_|
     565              :       !
     566              : 
     567           46 :       color_col = fm_mat_S%matrix_struct%context%mepos(2)
     568           46 :       ALLOCATE (para_env_col)
     569           46 :       CALL para_env_col%from_split(para_env_RPA, color_col)
     570              : 
     571           46 :       CALL para_env_col%sum(iaia_RI_dp)
     572              : 
     573              :       ! convert the iaia_RI_dp into double-double precision
     574          136 :       ALLOCATE (iaia_RI(ncol_local))
     575         2764 :       DO iiB = 1, ncol_local
     576         2764 :          iaia_RI(iiB) = iaia_RI_dp(iiB)
     577              :       END DO
     578           46 :       DEALLOCATE (iaia_RI_dp)
     579              : 
     580              :       ! 4) calculate the right hand term, D_ia is the matrix containing the
     581              :       ! orbital energy differences, M_ia is the diagonal of the full RPA 'excitation'
     582              :       ! matrix
     583          136 :       ALLOCATE (D_ia(ncol_local))
     584              : 
     585           90 :       ALLOCATE (M_ia(ncol_local))
     586              : 
     587         2764 :       DO iiB = 1, ncol_local
     588         2718 :          i_global = col_indices(iiB)
     589              : 
     590         2718 :          iocc = MAX(1, i_global - 1)/virtual + 1
     591         2718 :          avirt = i_global - (iocc - 1)*virtual
     592         2718 :          eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
     593              : 
     594         2764 :          D_ia(iiB) = eigen_diff
     595              :       END DO
     596              : 
     597         2764 :       DO iiB = 1, ncol_local
     598         2764 :          M_ia(iiB) = D_ia(iiB)*D_ia(iiB) + 2.0_dp*D_ia(iiB)*iaia_RI(iiB)
     599              :       END DO
     600              : 
     601           46 :       right_term_ref = 0.0_dp
     602         2764 :       DO iiB = 1, ncol_local
     603         2764 :          right_term_ref = right_term_ref + (SQRT(M_ia(iiB)) - D_ia(iiB) - iaia_RI(iiB))
     604              :       END DO
     605           46 :       right_term_ref = right_term_ref/2.0_dp
     606              : 
     607              :       ! sum the result with the processes of the RPA_group having the same row
     608           46 :       color_row = fm_mat_S%matrix_struct%context%mepos(1)
     609           46 :       ALLOCATE (para_env_row)
     610           46 :       CALL para_env_row%from_split(para_env_RPA, color_row)
     611              : 
     612              :       ! allocate communication array for rows
     613           46 :       CALL para_env_row%sum(right_term_ref)
     614              : 
     615           46 :       CALL mp_para_env_release(para_env_row)
     616              : 
     617           46 :       CALL timestop(handle)
     618              : 
     619           46 :    END SUBROUTINE calc_ia_ia_integrals
     620              : 
     621              : ! **************************************************************************************************
     622              : !> \brief ...
     623              : !> \param a_scaling ...
     624              : !> \param left_term ...
     625              : !> \param first_deriv ...
     626              : !> \param num_integ_points ...
     627              : !> \param my_open_shell ...
     628              : !> \param M_ia ...
     629              : !> \param cottj ...
     630              : !> \param wj ...
     631              : !> \param D_ia ...
     632              : !> \param D_ia_beta ...
     633              : !> \param M_ia_beta ...
     634              : !> \param ncol_local ...
     635              : !> \param ncol_local_beta ...
     636              : !> \param num_integ_group ...
     637              : !> \param color_rpa_group ...
     638              : !> \param para_env ...
     639              : !> \param para_env_col ...
     640              : !> \param para_env_col_beta ...
     641              : ! **************************************************************************************************
     642          428 :    SUBROUTINE calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
     643              :                                 M_ia, cottj, wj, D_ia, D_ia_beta, M_ia_beta, &
     644              :                                 ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
     645              :                                 para_env, para_env_col, para_env_col_beta)
     646              :       REAL(KIND=dp), INTENT(IN)                          :: a_scaling
     647              :       REAL(KIND=dp), INTENT(INOUT)                       :: left_term, first_deriv
     648              :       INTEGER, INTENT(IN)                                :: num_integ_points
     649              :       LOGICAL, INTENT(IN)                                :: my_open_shell
     650              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     651              :          INTENT(IN)                                      :: M_ia, cottj, wj, D_ia, D_ia_beta, &
     652              :                                                             M_ia_beta
     653              :       INTEGER, INTENT(IN)                                :: ncol_local, ncol_local_beta, &
     654              :                                                             num_integ_group, color_rpa_group
     655              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env, para_env_col
     656              :       TYPE(mp_para_env_type), POINTER                    :: para_env_col_beta
     657              : 
     658              :       INTEGER                                            :: iiB, jquad
     659              :       REAL(KIND=dp)                                      :: first_deriv_beta, left_term_beta, omega
     660              : 
     661          428 :       left_term = 0.0_dp
     662          428 :       first_deriv = 0.0_dp
     663          428 :       left_term_beta = 0.0_dp
     664          428 :       first_deriv_beta = 0.0_dp
     665         4452 :       DO jquad = 1, num_integ_points
     666              :          ! parallelize over integration points
     667         4024 :          IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
     668         2212 :          omega = a_scaling*cottj(jquad)
     669              : 
     670       162484 :          DO iiB = 1, ncol_local
     671              :             ! parallelize over ia elements in the para_env_row group
     672       160272 :             IF (MODULO(iiB, para_env_col%num_pe) /= para_env_col%mepos) CYCLE
     673              :             ! calculate left_term
     674              :             left_term = left_term + wj(jquad)* &
     675              :                         (LOG(1.0_dp + (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2)) - &
     676       145072 :                          (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2))
     677              :             first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
     678       162484 :                           ((-M_ia(iiB) + D_ia(iiB)**2)**2/((omega**2 + D_ia(iiB)**2)**2*(omega**2 + M_ia(iiB))))
     679              :          END DO
     680              : 
     681         2640 :          IF (my_open_shell) THEN
     682        14490 :             DO iiB = 1, ncol_local_beta
     683              :                ! parallelize over ia elements in the para_env_row group
     684        14140 :                IF (MODULO(iiB, para_env_col_beta%num_pe) /= para_env_col_beta%mepos) CYCLE
     685              :                ! calculate left_term
     686              :                left_term_beta = left_term_beta + wj(jquad)* &
     687              :                                 (LOG(1.0_dp + (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2)) - &
     688        14140 :                                  (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2))
     689              :                first_deriv_beta = &
     690              :                   first_deriv_beta + wj(jquad)*cottj(jquad)**2* &
     691        14490 :                   ((-M_ia_beta(iiB) + D_ia_beta(iiB)**2)**2/((omega**2 + D_ia_beta(iiB)**2)**2*(omega**2 + M_ia_beta(iiB))))
     692              :             END DO
     693              :          END IF
     694              : 
     695              :       END DO
     696              : 
     697              :       ! sum the contribution from all proc, starting form the row group
     698          428 :       CALL para_env%sum(left_term)
     699          428 :       CALL para_env%sum(first_deriv)
     700              : 
     701          428 :       IF (my_open_shell) THEN
     702           70 :          CALL para_env%sum(left_term_beta)
     703           70 :          CALL para_env%sum(first_deriv_beta)
     704              : 
     705           70 :          left_term = left_term + left_term_beta
     706           70 :          first_deriv = first_deriv + first_deriv_beta
     707              :       END IF
     708              : 
     709          428 :    END SUBROUTINE calculate_objfunc
     710              : 
     711              : ! **************************************************************************************************
     712              : !> \brief Calculate integration weights for the tau grid (in dependency of the omega node)
     713              : !> \param num_integ_points ...
     714              : !> \param tau_tj ...
     715              : !> \param weights_cos_tf_t_to_w ...
     716              : !> \param omega_tj ...
     717              : !> \param E_min ...
     718              : !> \param E_max ...
     719              : !> \param max_error ...
     720              : !> \param num_points_per_magnitude ...
     721              : !> \param regularization ...
     722              : ! **************************************************************************************************
     723           70 :    SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, &
     724              :                                            E_min, E_max, max_error, num_points_per_magnitude, &
     725              :                                            regularization)
     726              : 
     727              :       INTEGER, INTENT(IN)                                :: num_integ_points
     728              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     729              :          INTENT(IN)                                      :: tau_tj
     730              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     731              :          INTENT(INOUT)                                   :: weights_cos_tf_t_to_w
     732              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     733              :          INTENT(IN)                                      :: omega_tj
     734              :       REAL(KIND=dp), INTENT(IN)                          :: E_min, E_max
     735              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_error
     736              :       INTEGER, INTENT(IN)                                :: num_points_per_magnitude
     737              :       REAL(KIND=dp), INTENT(IN)                          :: regularization
     738              : 
     739              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_t_to_w'
     740              : 
     741              :       INTEGER                                            :: handle, iii, info, jjj, jquad, lwork, &
     742              :                                                             num_x_nodes
     743           70 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
     744              :       REAL(KIND=dp)                                      :: multiplicator, omega
     745           70 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: sing_values, tau_wj_work, vec_UTy, work, &
     746              :                                                             x_values, y_values
     747           70 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: mat_A, mat_SinvVSinvSigma, &
     748           70 :                                                             mat_SinvVSinvT, mat_U
     749              : 
     750           70 :       CALL timeset(routineN, handle)
     751              : 
     752              :       ! take num_points_per_magnitude points per 10-interval
     753           70 :       num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
     754              : 
     755              :       ! take at least as many x points as integration points to have clear
     756              :       ! input for the singular value decomposition
     757           70 :       num_x_nodes = MAX(num_x_nodes, num_integ_points)
     758              : 
     759          210 :       ALLOCATE (x_values(num_x_nodes))
     760        24070 :       x_values = 0.0_dp
     761          140 :       ALLOCATE (y_values(num_x_nodes))
     762        24070 :       y_values = 0.0_dp
     763          280 :       ALLOCATE (mat_A(num_x_nodes, num_integ_points))
     764       150976 :       mat_A = 0.0_dp
     765          210 :       ALLOCATE (tau_wj_work(num_integ_points))
     766          576 :       tau_wj_work = 0.0_dp
     767          140 :       ALLOCATE (sing_values(num_integ_points))
     768          576 :       sing_values = 0.0_dp
     769          280 :       ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
     770      9144070 :       mat_U = 0.0_dp
     771          210 :       ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
     772              : 
     773       150976 :       mat_SinvVSinvT = 0.0_dp
     774              :       ! double the value nessary for 'A' to achieve good performance
     775           70 :       lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
     776          210 :       ALLOCATE (work(lwork))
     777       105230 :       work = 0.0_dp
     778          210 :       ALLOCATE (iwork(8*num_integ_points))
     779         4118 :       iwork = 0
     780          210 :       ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
     781       174470 :       mat_SinvVSinvSigma = 0.0_dp
     782          140 :       ALLOCATE (vec_UTy(num_x_nodes))
     783        24070 :       vec_UTy = 0.0_dp
     784              : 
     785           70 :       max_error = 0.0_dp
     786              : 
     787              :       ! loop over all omega frequency points
     788          576 :       DO jquad = 1, num_integ_points
     789              : 
     790              :          ! set the x-values logarithmically in the interval [Emin,Emax]
     791          506 :          multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
     792       150906 :          DO iii = 1, num_x_nodes
     793       150906 :             x_values(iii) = E_min*multiplicator**(iii - 1)
     794              :          END DO
     795              : 
     796          506 :          omega = omega_tj(jquad)
     797              : 
     798              :          ! y=2x/(x^2+omega_k^2)
     799       150906 :          DO iii = 1, num_x_nodes
     800       150906 :             y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2)
     801              :          END DO
     802              : 
     803              :          ! calculate mat_A
     804         6892 :          DO jjj = 1, num_integ_points
     805      1590892 :             DO iii = 1, num_x_nodes
     806      1590386 :                mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
     807              :             END DO
     808              :          END DO
     809              : 
     810              :          ! Singular value decomposition of mat_A
     811              :          CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
     812          506 :                      mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
     813              : 
     814          506 :          CPASSERT(info == 0)
     815              : 
     816              :          ! integration weights = V Sigma U^T y
     817              :          ! 1) V*Sigma
     818         6892 :          DO jjj = 1, num_integ_points
     819       114582 :             DO iii = 1, num_integ_points
     820              : !               mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
     821              :                mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
     822       114076 :                                               /(regularization**2 + sing_values(jjj)**2)
     823              :             END DO
     824              :          END DO
     825              : 
     826              :          ! 2) U^T y
     827              :          CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
     828          506 :                     0.0_dp, vec_UTy, num_x_nodes)
     829              : 
     830              :          ! 3) (V*Sigma) * (U^T y)
     831              :          CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
     832          506 :                     num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
     833              : 
     834         6892 :          weights_cos_tf_t_to_w(jquad, :) = tau_wj_work(:)
     835              : 
     836              :          CALL calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, &
     837          576 :                                                       y_values, num_integ_points, num_x_nodes)
     838              : 
     839              :       END DO ! jquad
     840              : 
     841            0 :       DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, sing_values, mat_U, mat_SinvVSinvT, &
     842           70 :                   work, iwork, mat_SinvVSinvSigma, vec_UTy)
     843              : 
     844           70 :       CALL timestop(handle)
     845              : 
     846           70 :    END SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w
     847              : 
     848              : ! **************************************************************************************************
     849              : !> \brief Calculate integration weights for the tau grid (in dependency of the omega node)
     850              : !> \param num_integ_points ...
     851              : !> \param tau_tj ...
     852              : !> \param weights_sin_tf_t_to_w ...
     853              : !> \param omega_tj ...
     854              : !> \param E_min ...
     855              : !> \param E_max ...
     856              : !> \param max_error ...
     857              : !> \param num_points_per_magnitude ...
     858              : !> \param regularization ...
     859              : ! **************************************************************************************************
     860           38 :    SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, &
     861              :                                            E_min, E_max, max_error, num_points_per_magnitude, regularization)
     862              : 
     863              :       INTEGER, INTENT(IN)                                :: num_integ_points
     864              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     865              :          INTENT(IN)                                      :: tau_tj
     866              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     867              :          INTENT(INOUT)                                   :: weights_sin_tf_t_to_w
     868              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     869              :          INTENT(IN)                                      :: omega_tj
     870              :       REAL(KIND=dp), INTENT(IN)                          :: E_min, E_max
     871              :       REAL(KIND=dp), INTENT(OUT)                         :: max_error
     872              :       INTEGER, INTENT(IN)                                :: num_points_per_magnitude
     873              :       REAL(KIND=dp), INTENT(IN)                          :: regularization
     874              : 
     875              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_sin_tf_t_to_w'
     876              : 
     877              :       INTEGER                                            :: handle, iii, info, jjj, jquad, lwork, &
     878              :                                                             num_x_nodes
     879           38 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
     880              :       REAL(KIND=dp)                                      :: chi2_min_jquad, multiplicator, omega
     881           38 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: sing_values, tau_wj_work, vec_UTy, work, &
     882           38 :                                                             work_array, x_values, y_values
     883           38 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: mat_A, mat_SinvVSinvSigma, &
     884           38 :                                                             mat_SinvVSinvT, mat_U
     885              : 
     886           38 :       CALL timeset(routineN, handle)
     887              : 
     888              :       ! take num_points_per_magnitude points per 10-interval
     889           38 :       num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
     890              : 
     891              :       ! take at least as many x points as integration points to have clear
     892              :       ! input for the singular value decomposition
     893           38 :       num_x_nodes = MAX(num_x_nodes, num_integ_points)
     894              : 
     895          114 :       ALLOCATE (x_values(num_x_nodes))
     896        12838 :       x_values = 0.0_dp
     897           76 :       ALLOCATE (y_values(num_x_nodes))
     898        12838 :       y_values = 0.0_dp
     899          152 :       ALLOCATE (mat_A(num_x_nodes, num_integ_points))
     900       119656 :       mat_A = 0.0_dp
     901          114 :       ALLOCATE (tau_wj_work(num_integ_points))
     902          456 :       tau_wj_work = 0.0_dp
     903          114 :       ALLOCATE (work_array(2*num_integ_points))
     904          874 :       work_array = 0.0_dp
     905           76 :       ALLOCATE (sing_values(num_integ_points))
     906          456 :       sing_values = 0.0_dp
     907          152 :       ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
     908      4972838 :       mat_U = 0.0_dp
     909          114 :       ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
     910              : 
     911       119656 :       mat_SinvVSinvT = 0.0_dp
     912              :       ! double the value nessary for 'A' to achieve good performance
     913           38 :       lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
     914          114 :       ALLOCATE (work(lwork))
     915        79758 :       work = 0.0_dp
     916          114 :       ALLOCATE (iwork(8*num_integ_points))
     917         3382 :       iwork = 0
     918          114 :       ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
     919       132038 :       mat_SinvVSinvSigma = 0.0_dp
     920           76 :       ALLOCATE (vec_UTy(num_x_nodes))
     921        12838 :       vec_UTy = 0.0_dp
     922              : 
     923           38 :       max_error = 0.0_dp
     924              : 
     925              :       ! loop over all omega frequency points
     926          456 :       DO jquad = 1, num_integ_points
     927              : 
     928          418 :          chi2_min_jquad = 100.0_dp
     929              : 
     930              :          ! set the x-values logarithmically in the interval [Emin,Emax]
     931          418 :          multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
     932       119618 :          DO iii = 1, num_x_nodes
     933       119618 :             x_values(iii) = E_min*multiplicator**(iii - 1)
     934              :          END DO
     935              : 
     936          418 :          omega = omega_tj(jquad)
     937              : 
     938              :          ! y=2x/(x^2+omega_k^2)
     939       119618 :          DO iii = 1, num_x_nodes
     940              : !            y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2)
     941       119618 :             y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2)
     942              :          END DO
     943              : 
     944              :          ! calculate mat_A
     945         6556 :          DO jjj = 1, num_integ_points
     946      1501756 :             DO iii = 1, num_x_nodes
     947      1501338 :                mat_A(iii, jjj) = SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
     948              :             END DO
     949              :          END DO
     950              : 
     951              :          ! Singular value decomposition of mat_A
     952              :          CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
     953          418 :                      mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
     954              : 
     955          418 :          CPASSERT(info == 0)
     956              : 
     957              :          ! integration weights = V Sigma U^T y
     958              :          ! 1) V*Sigma
     959         6556 :          DO jjj = 1, num_integ_points
     960       113534 :             DO iii = 1, num_integ_points
     961              : !               mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
     962              :                mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
     963       113116 :                                               /(regularization**2 + sing_values(jjj)**2)
     964              :             END DO
     965              :          END DO
     966              : 
     967              :          ! 2) U^T y
     968              :          CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
     969          418 :                     0.0_dp, vec_UTy, num_x_nodes)
     970              : 
     971              :          ! 3) (V*Sigma) * (U^T y)
     972              :          CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
     973          418 :                     num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
     974              : 
     975         6556 :          weights_sin_tf_t_to_w(jquad, :) = tau_wj_work(:)
     976              : 
     977              :          CALL calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, &
     978          456 :                                                     y_values, num_integ_points, num_x_nodes)
     979              : 
     980              :       END DO ! jquad
     981              : 
     982            0 :       DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
     983           38 :                   work, iwork, mat_SinvVSinvSigma, vec_UTy)
     984              : 
     985           38 :       CALL timestop(handle)
     986              : 
     987           38 :    END SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w
     988              : 
     989              : ! **************************************************************************************************
     990              : !> \brief ...
     991              : !> \param max_error ...
     992              : !> \param omega ...
     993              : !> \param tau_tj ...
     994              : !> \param tau_wj_work ...
     995              : !> \param x_values ...
     996              : !> \param y_values ...
     997              : !> \param num_integ_points ...
     998              : !> \param num_x_nodes ...
     999              : ! **************************************************************************************************
    1000          506 :    PURE SUBROUTINE calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, &
    1001              :                                                            y_values, num_integ_points, num_x_nodes)
    1002              : 
    1003              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_error, omega
    1004              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1005              :          INTENT(IN)                                      :: tau_tj, tau_wj_work, x_values, y_values
    1006              :       INTEGER, INTENT(IN)                                :: num_integ_points, num_x_nodes
    1007              : 
    1008              :       INTEGER                                            :: kkk
    1009              :       REAL(KIND=dp)                                      :: func_val, func_val_temp, max_error_tmp
    1010              : 
    1011          506 :       max_error_tmp = 0.0_dp
    1012              : 
    1013       150906 :       DO kkk = 1, num_x_nodes
    1014              : 
    1015              :          func_val = 0.0_dp
    1016              : 
    1017       150400 :          CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
    1018              : 
    1019       150906 :          IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
    1020         2494 :             max_error_tmp = ABS(y_values(kkk) - func_val)
    1021         2494 :             func_val_temp = func_val
    1022              :          END IF
    1023              : 
    1024              :       END DO
    1025              : 
    1026          506 :       IF (max_error_tmp > max_error) THEN
    1027              : 
    1028          144 :          max_error = max_error_tmp
    1029              : 
    1030              :       END IF
    1031              : 
    1032          506 :    END SUBROUTINE calc_max_error_fit_tau_grid_with_cosine
    1033              : 
    1034              : ! **************************************************************************************************
    1035              : !> \brief Evaluate fit function when calculating tau grid for cosine transform
    1036              : !> \param func_val ...
    1037              : !> \param x_value ...
    1038              : !> \param num_integ_points ...
    1039              : !> \param tau_tj ...
    1040              : !> \param tau_wj_work ...
    1041              : !> \param omega ...
    1042              : ! **************************************************************************************************
    1043       150400 :    PURE SUBROUTINE eval_fit_func_tau_grid_cosine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega)
    1044              : 
    1045              :       REAL(KIND=dp), INTENT(OUT)                         :: func_val
    1046              :       REAL(KIND=dp), INTENT(IN)                          :: x_value
    1047              :       INTEGER, INTENT(IN)                                :: num_integ_points
    1048              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1049              :          INTENT(IN)                                      :: tau_tj, tau_wj_work
    1050              :       REAL(KIND=dp), INTENT(IN)                          :: omega
    1051              : 
    1052              :       INTEGER                                            :: iii
    1053              : 
    1054       150400 :       func_val = 0.0_dp
    1055              : 
    1056      1734400 :       DO iii = 1, num_integ_points
    1057              : 
    1058              :          ! calculate value of the fit function
    1059      1734400 :          func_val = func_val + tau_wj_work(iii)*COS(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii))
    1060              : 
    1061              :       END DO
    1062              : 
    1063       150400 :    END SUBROUTINE eval_fit_func_tau_grid_cosine
    1064              : 
    1065              : ! **************************************************************************************************
    1066              : !> \brief Evaluate fit function when calculating tau grid for sine transform
    1067              : !> \param func_val ...
    1068              : !> \param x_value ...
    1069              : !> \param num_integ_points ...
    1070              : !> \param tau_tj ...
    1071              : !> \param tau_wj_work ...
    1072              : !> \param omega ...
    1073              : ! **************************************************************************************************
    1074       119200 :    PURE SUBROUTINE eval_fit_func_tau_grid_sine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega)
    1075              : 
    1076              :       REAL(KIND=dp), INTENT(INOUT)                       :: func_val
    1077              :       REAL(KIND=dp), INTENT(IN)                          :: x_value
    1078              :       INTEGER, INTENT(in)                                :: num_integ_points
    1079              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1080              :          INTENT(IN)                                      :: tau_tj, tau_wj_work
    1081              :       REAL(KIND=dp), INTENT(IN)                          :: omega
    1082              : 
    1083              :       INTEGER                                            :: iii
    1084              : 
    1085       119200 :       func_val = 0.0_dp
    1086              : 
    1087      1614400 :       DO iii = 1, num_integ_points
    1088              : 
    1089              :          ! calculate value of the fit function
    1090      1614400 :          func_val = func_val + tau_wj_work(iii)*SIN(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii))
    1091              : 
    1092              :       END DO
    1093              : 
    1094       119200 :    END SUBROUTINE eval_fit_func_tau_grid_sine
    1095              : 
    1096              : ! **************************************************************************************************
    1097              : !> \brief ...
    1098              : !> \param max_error ...
    1099              : !> \param omega ...
    1100              : !> \param tau_tj ...
    1101              : !> \param tau_wj_work ...
    1102              : !> \param x_values ...
    1103              : !> \param y_values ...
    1104              : !> \param num_integ_points ...
    1105              : !> \param num_x_nodes ...
    1106              : ! **************************************************************************************************
    1107          418 :    PURE SUBROUTINE calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, &
    1108              :                                                          y_values, num_integ_points, num_x_nodes)
    1109              : 
    1110              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_error, omega
    1111              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1112              :          INTENT(IN)                                      :: tau_tj, tau_wj_work, x_values, y_values
    1113              :       INTEGER, INTENT(IN)                                :: num_integ_points, num_x_nodes
    1114              : 
    1115              :       INTEGER                                            :: kkk
    1116              :       REAL(KIND=dp)                                      :: func_val, func_val_temp, max_error_tmp
    1117              : 
    1118          418 :       max_error_tmp = 0.0_dp
    1119              : 
    1120       119618 :       DO kkk = 1, num_x_nodes
    1121              : 
    1122       119200 :          func_val = 0.0_dp
    1123              : 
    1124       119200 :          CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
    1125              : 
    1126       119618 :          IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
    1127         2410 :             max_error_tmp = ABS(y_values(kkk) - func_val)
    1128         2410 :             func_val_temp = func_val
    1129              :          END IF
    1130              : 
    1131              :       END DO
    1132              : 
    1133          418 :       IF (max_error_tmp > max_error) THEN
    1134              : 
    1135           42 :          max_error = max_error_tmp
    1136              : 
    1137              :       END IF
    1138              : 
    1139          418 :    END SUBROUTINE calc_max_error_fit_tau_grid_with_sine
    1140              : 
    1141              : ! **************************************************************************************************
    1142              : !> \brief test the singular value decomposition for the computation of integration weights for the
    1143              : !>         Fourier transform between time and frequency grid in cubic-scaling RPA
    1144              : !> \param nR ...
    1145              : !> \param iw ...
    1146              : ! **************************************************************************************************
    1147            0 :    SUBROUTINE test_least_square_ft(nR, iw)
    1148              :       INTEGER, INTENT(IN)                                :: nR, iw
    1149              : 
    1150              :       INTEGER                                            :: ierr, iR, jquad, num_integ_points
    1151              :       REAL(KIND=dp)                                      :: max_error, multiplicator, Rc, Rc_max
    1152            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: tau_tj, tau_wj, tj, wj, x_tw
    1153            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: weights_cos_tf_t_to_w
    1154              : 
    1155            0 :       Rc_max = 1.0E+7
    1156              : 
    1157            0 :       multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp))
    1158              : 
    1159            0 :       DO num_integ_points = 1, 20
    1160              : 
    1161            0 :          ALLOCATE (x_tw(2*num_integ_points))
    1162            0 :          x_tw = 0.0_dp
    1163            0 :          ALLOCATE (tau_tj(num_integ_points))
    1164            0 :          tau_tj = 0.0_dp
    1165            0 :          ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
    1166            0 :          weights_cos_tf_t_to_w = 0.0_dp
    1167            0 :          ALLOCATE (tau_wj(num_integ_points))
    1168            0 :          tau_wj = 0.0_dp
    1169            0 :          ALLOCATE (tj(num_integ_points))
    1170            0 :          tj = 0.0_dp
    1171            0 :          ALLOCATE (wj(num_integ_points))
    1172            0 :          wj = 0.0_dp
    1173              : 
    1174            0 :          DO iR = 0, nR - 1
    1175              : 
    1176            0 :             Rc = 2.0_dp*multiplicator**iR
    1177              : 
    1178            0 :             ierr = 0
    1179            0 :             CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.)
    1180              : 
    1181            0 :             DO jquad = 1, num_integ_points
    1182            0 :                tj(jquad) = x_tw(jquad)
    1183            0 :                wj(jquad) = x_tw(jquad + num_integ_points)
    1184              :             END DO
    1185              : 
    1186            0 :             x_tw = 0.0_dp
    1187              : 
    1188            0 :             CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw)
    1189              : 
    1190            0 :             DO jquad = 1, num_integ_points
    1191            0 :                tau_tj(jquad) = x_tw(jquad)/2.0_dp
    1192            0 :                tau_wj(jquad) = x_tw(jquad + num_integ_points)/2.0_dp
    1193              :             END DO
    1194              : 
    1195              :             CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, &
    1196              :                                               weights_cos_tf_t_to_w, tj, &
    1197            0 :                                               1.0_dp, Rc, max_error, 200, 0.0_dp)
    1198              : 
    1199            0 :             IF (iw > 0) THEN
    1200            0 :                WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error
    1201              :             END IF
    1202              : 
    1203              :          END DO
    1204              : 
    1205            0 :          DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj)
    1206              : 
    1207              :       END DO
    1208              : 
    1209            0 :    END SUBROUTINE test_least_square_ft
    1210              : 
    1211              : ! **************************************************************************************************
    1212              : !> \brief ...
    1213              : !> \param num_integ_points ...
    1214              : !> \param tau_tj ...
    1215              : !> \param weights_cos_tf_w_to_t ...
    1216              : !> \param omega_tj ...
    1217              : !> \param E_min ...
    1218              : !> \param E_max ...
    1219              : !> \param max_error ...
    1220              : !> \param num_points_per_magnitude ...
    1221              : !> \param regularization ...
    1222              : ! **************************************************************************************************
    1223           70 :    SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, &
    1224              :                                            E_min, E_max, max_error, num_points_per_magnitude, regularization)
    1225              : 
    1226              :       INTEGER, INTENT(IN)                                :: num_integ_points
    1227              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1228              :          INTENT(IN)                                      :: tau_tj
    1229              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
    1230              :          INTENT(INOUT)                                   :: weights_cos_tf_w_to_t
    1231              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1232              :          INTENT(IN)                                      :: omega_tj
    1233              :       REAL(KIND=dp), INTENT(IN)                          :: E_min, E_max
    1234              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_error
    1235              :       INTEGER, INTENT(IN)                                :: num_points_per_magnitude
    1236              :       REAL(KIND=dp), INTENT(IN)                          :: regularization
    1237              : 
    1238              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_w_to_t'
    1239              : 
    1240              :       INTEGER                                            :: handle, iii, info, jjj, jquad, lwork, &
    1241              :                                                             num_x_nodes
    1242           70 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
    1243              :       REAL(KIND=dp)                                      :: chi2_min_jquad, multiplicator, omega, &
    1244              :                                                             tau, x_value
    1245           70 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: omega_wj_work, sing_values, vec_UTy, &
    1246           70 :                                                             work, work_array, x_values, y_values
    1247           70 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: mat_A, mat_SinvVSinvSigma, &
    1248           70 :                                                             mat_SinvVSinvT, mat_U
    1249              : 
    1250           70 :       CALL timeset(routineN, handle)
    1251              : 
    1252              :       ! take num_points_per_magnitude points per 10-interval
    1253           70 :       num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
    1254              : 
    1255              :       ! take at least as many x points as integration points to have clear
    1256              :       ! input for the singular value decomposition
    1257           70 :       num_x_nodes = MAX(num_x_nodes, num_integ_points)
    1258              : 
    1259          210 :       ALLOCATE (x_values(num_x_nodes))
    1260        24070 :       x_values = 0.0_dp
    1261          140 :       ALLOCATE (y_values(num_x_nodes))
    1262        24070 :       y_values = 0.0_dp
    1263          280 :       ALLOCATE (mat_A(num_x_nodes, num_integ_points))
    1264       150976 :       mat_A = 0.0_dp
    1265          210 :       ALLOCATE (omega_wj_work(num_integ_points))
    1266          576 :       omega_wj_work = 0.0_dp
    1267          210 :       ALLOCATE (work_array(2*num_integ_points))
    1268         1082 :       work_array = 0.0_dp
    1269          140 :       ALLOCATE (sing_values(num_integ_points))
    1270          576 :       sing_values = 0.0_dp
    1271          280 :       ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
    1272      9144070 :       mat_U = 0.0_dp
    1273          210 :       ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
    1274              : 
    1275       150976 :       mat_SinvVSinvT = 0.0_dp
    1276              :       ! double the value nessary for 'A' to achieve good performance
    1277           70 :       lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
    1278          210 :       ALLOCATE (work(lwork))
    1279       105230 :       work = 0.0_dp
    1280          210 :       ALLOCATE (iwork(8*num_integ_points))
    1281         4118 :       iwork = 0
    1282          210 :       ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
    1283       174470 :       mat_SinvVSinvSigma = 0.0_dp
    1284          140 :       ALLOCATE (vec_UTy(num_x_nodes))
    1285        24070 :       vec_UTy = 0.0_dp
    1286              : 
    1287              :       ! set the x-values logarithmically in the interval [Emin,Emax]
    1288           70 :       multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
    1289        24070 :       DO iii = 1, num_x_nodes
    1290        24070 :          x_values(iii) = E_min*multiplicator**(iii - 1)
    1291              :       END DO
    1292              : 
    1293           70 :       max_error = 0.0_dp
    1294              : 
    1295              :       ! loop over all tau time points
    1296          576 :       DO jquad = 1, num_integ_points
    1297              : 
    1298          506 :          chi2_min_jquad = 100.0_dp
    1299              : 
    1300          506 :          tau = tau_tj(jquad)
    1301              : 
    1302              :          ! y=exp(-x*|tau_k|)
    1303       150906 :          DO iii = 1, num_x_nodes
    1304       150906 :             y_values(iii) = EXP(-x_values(iii)*tau)
    1305              :          END DO
    1306              : 
    1307              :          ! calculate mat_A
    1308         6892 :          DO jjj = 1, num_integ_points
    1309      1590892 :             DO iii = 1, num_x_nodes
    1310      1584000 :                omega = omega_tj(jjj)
    1311      1584000 :                x_value = x_values(iii)
    1312      1590386 :                mat_A(iii, jjj) = COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
    1313              :             END DO
    1314              :          END DO
    1315              : 
    1316              :          ! Singular value decomposition of mat_A
    1317              :          CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
    1318          506 :                      mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
    1319              : 
    1320          506 :          CPASSERT(info == 0)
    1321              : 
    1322              :          ! integration weights = V Sigma U^T y
    1323              :          ! 1) V*Sigma
    1324         6892 :          DO jjj = 1, num_integ_points
    1325       114582 :             DO iii = 1, num_integ_points
    1326              : !               mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
    1327              :                mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
    1328       114076 :                                               /(regularization**2 + sing_values(jjj)**2)
    1329              :             END DO
    1330              :          END DO
    1331              : 
    1332              :          ! 2) U^T y
    1333              :          CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
    1334          506 :                     0.0_dp, vec_UTy, num_x_nodes)
    1335              : 
    1336              :          ! 3) (V*Sigma) * (U^T y)
    1337              :          CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
    1338          506 :                     num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points)
    1339              : 
    1340         6892 :          weights_cos_tf_w_to_t(jquad, :) = omega_wj_work(:)
    1341              : 
    1342              :          CALL calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, &
    1343          576 :                                                         y_values, num_integ_points, num_x_nodes)
    1344              : 
    1345              :       END DO ! jquad
    1346              : 
    1347            0 :       DEALLOCATE (x_values, y_values, mat_A, omega_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
    1348           70 :                   work, iwork, mat_SinvVSinvSigma, vec_UTy)
    1349              : 
    1350           70 :       CALL timestop(handle)
    1351              : 
    1352           70 :    END SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t
    1353              : 
    1354              : ! **************************************************************************************************
    1355              : !> \brief ...
    1356              : !> \param max_error ...
    1357              : !> \param tau ...
    1358              : !> \param omega_tj ...
    1359              : !> \param omega_wj_work ...
    1360              : !> \param x_values ...
    1361              : !> \param y_values ...
    1362              : !> \param num_integ_points ...
    1363              : !> \param num_x_nodes ...
    1364              : ! **************************************************************************************************
    1365          506 :    SUBROUTINE calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, &
    1366              :                                                         y_values, num_integ_points, num_x_nodes)
    1367              : 
    1368              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_error
    1369              :       REAL(KIND=dp), INTENT(IN)                          :: tau
    1370              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1371              :          INTENT(IN)                                      :: omega_tj, omega_wj_work, x_values, &
    1372              :                                                             y_values
    1373              :       INTEGER, INTENT(IN)                                :: num_integ_points, num_x_nodes
    1374              : 
    1375              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_max_error_fit_omega_grid_with_cosine'
    1376              : 
    1377              :       INTEGER                                            :: handle, kkk
    1378              :       REAL(KIND=dp)                                      :: func_val, func_val_temp, max_error_tmp
    1379              : 
    1380          506 :       CALL timeset(routineN, handle)
    1381              : 
    1382          506 :       max_error_tmp = 0.0_dp
    1383              : 
    1384       150906 :       DO kkk = 1, num_x_nodes
    1385              : 
    1386              :          func_val = 0.0_dp
    1387              : 
    1388       150400 :          CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau)
    1389              : 
    1390       150906 :          IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
    1391         4244 :             max_error_tmp = ABS(y_values(kkk) - func_val)
    1392         4244 :             func_val_temp = func_val
    1393              :          END IF
    1394              : 
    1395              :       END DO
    1396              : 
    1397          506 :       IF (max_error_tmp > max_error) THEN
    1398              : 
    1399          132 :          max_error = max_error_tmp
    1400              : 
    1401              :       END IF
    1402              : 
    1403          506 :       CALL timestop(handle)
    1404              : 
    1405          506 :    END SUBROUTINE calc_max_error_fit_omega_grid_with_cosine
    1406              : 
    1407              : ! **************************************************************************************************
    1408              : !> \brief ...
    1409              : !> \param func_val ...
    1410              : !> \param x_value ...
    1411              : !> \param num_integ_points ...
    1412              : !> \param omega_tj ...
    1413              : !> \param omega_wj_work ...
    1414              : !> \param tau ...
    1415              : ! **************************************************************************************************
    1416       150400 :    PURE SUBROUTINE eval_fit_func_omega_grid_cosine(func_val, x_value, num_integ_points, omega_tj, omega_wj_work, tau)
    1417              :       REAL(KIND=dp), INTENT(OUT)                         :: func_val
    1418              :       REAL(KIND=dp), INTENT(IN)                          :: x_value
    1419              :       INTEGER, INTENT(IN)                                :: num_integ_points
    1420              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1421              :          INTENT(IN)                                      :: omega_tj, omega_wj_work
    1422              :       REAL(KIND=dp), INTENT(IN)                          :: tau
    1423              : 
    1424              :       INTEGER                                            :: iii
    1425              :       REAL(KIND=dp)                                      :: omega
    1426              : 
    1427       150400 :       func_val = 0.0_dp
    1428              : 
    1429      1734400 :       DO iii = 1, num_integ_points
    1430              : 
    1431              :          ! calculate value of the fit function
    1432      1584000 :          omega = omega_tj(iii)
    1433      1734400 :          func_val = func_val + omega_wj_work(iii)*COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
    1434              : 
    1435              :       END DO
    1436              : 
    1437       150400 :    END SUBROUTINE eval_fit_func_omega_grid_cosine
    1438              : 
    1439              : ! **************************************************************************************************
    1440              : !> \brief ...
    1441              : !> \param qs_env ...
    1442              : !> \param para_env ...
    1443              : !> \param gap ...
    1444              : !> \param max_eig_diff ...
    1445              : !> \param e_fermi ...
    1446              : ! **************************************************************************************************
    1447            8 :    SUBROUTINE gap_and_max_eig_diff_kpoints(qs_env, para_env, gap, max_eig_diff, e_fermi)
    1448              : 
    1449              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1450              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1451              :       REAL(KIND=dp), INTENT(OUT)                         :: gap, max_eig_diff, e_fermi
    1452              : 
    1453              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'gap_and_max_eig_diff_kpoints'
    1454              : 
    1455              :       INTEGER                                            :: handle, homo, ikpgr, ispin, kplocal, &
    1456              :                                                             nmo, nspin
    1457              :       INTEGER, DIMENSION(2)                              :: kp_range
    1458              :       REAL(KIND=dp)                                      :: e_homo, e_homo_temp, e_lumo, e_lumo_temp
    1459              :       REAL(KIND=dp), DIMENSION(3)                        :: tmp
    1460            4 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvalues
    1461              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1462              :       TYPE(kpoint_type), POINTER                         :: kpoint
    1463              :       TYPE(mo_set_type), POINTER                         :: mo_set
    1464              : 
    1465            4 :       CALL timeset(routineN, handle)
    1466              : 
    1467              :       CALL get_qs_env(qs_env, &
    1468            4 :                       kpoints=kpoint)
    1469              : 
    1470            4 :       mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
    1471            4 :       CALL get_mo_set(mo_set, nmo=nmo)
    1472              : 
    1473            4 :       CALL get_kpoint_info(kpoint, kp_range=kp_range)
    1474            4 :       kplocal = kp_range(2) - kp_range(1) + 1
    1475              : 
    1476            4 :       gap = 1000.0_dp
    1477            4 :       max_eig_diff = 0.0_dp
    1478            4 :       e_homo = -1000.0_dp
    1479            4 :       e_lumo = 1000.0_dp
    1480              : 
    1481           12 :       DO ikpgr = 1, kplocal
    1482            8 :          kp => kpoint%kp_env(ikpgr)%kpoint_env
    1483            8 :          nspin = SIZE(kp%mos, 2)
    1484           20 :          DO ispin = 1, nspin
    1485            8 :             mo_set => kp%mos(1, ispin)
    1486            8 :             CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
    1487            8 :             e_homo_temp = eigenvalues(homo)
    1488            8 :             e_lumo_temp = eigenvalues(homo + 1)
    1489              : 
    1490            8 :             IF (e_homo_temp > e_homo) e_homo = e_homo_temp
    1491            8 :             IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp
    1492           24 :             IF (eigenvalues(nmo) - eigenvalues(1) > max_eig_diff) max_eig_diff = eigenvalues(nmo) - eigenvalues(1)
    1493              : 
    1494              :          END DO
    1495              :       END DO
    1496              : 
    1497              :       ! Collect all three numbers in an array
    1498              :       ! Reverse sign of lumo to reduce number of MPI calls
    1499            4 :       tmp(1) = e_homo
    1500            4 :       tmp(2) = -e_lumo
    1501            4 :       tmp(3) = max_eig_diff
    1502            4 :       CALL para_env%max(tmp)
    1503              : 
    1504            4 :       gap = -tmp(2) - tmp(1)
    1505            4 :       e_fermi = (tmp(1) - tmp(2))*0.5_dp
    1506            4 :       max_eig_diff = tmp(3)
    1507              : 
    1508            4 :       CALL timestop(handle)
    1509              : 
    1510            4 :    END SUBROUTINE
    1511              : 
    1512              : ! **************************************************************************************************
    1513              : !> \brief returns minimal and maximal energy values for the E_range for the minimax grid selection
    1514              : !> \param qs_env ...
    1515              : !> \param para_env ...
    1516              : !> \param homo index of the homo level for the respective spin channel
    1517              : !> \param Eigenval eigenvalues
    1518              : !> \param do_ri_sos_laplace_mp2 flag for SOS-MP2
    1519              : !> \param do_kpoints_cubic_RPA flag for cubic-scaling RPA with k-points
    1520              : !> \param Emin minimal eigenvalue difference (gap of the system)
    1521              : !> \param Emax maximal eigenvalue difference
    1522              : !> \param e_range ...
    1523              : !> \param e_fermi Fermi level
    1524              : ! **************************************************************************************************
    1525          204 :    SUBROUTINE determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
    1526              :                                      do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
    1527              : 
    1528              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1529              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1530              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
    1531              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Eigenval
    1532              :       LOGICAL, INTENT(IN)                                :: do_ri_sos_laplace_mp2, &
    1533              :                                                             do_kpoints_cubic_RPA
    1534              :       REAL(KIND=dp), INTENT(OUT)                         :: Emin, Emax, e_range, e_fermi
    1535              : 
    1536              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'determine_energy_range'
    1537              : 
    1538              :       INTEGER                                            :: handle, ispin, nspins
    1539              :       LOGICAL                                            :: my_do_kpoints
    1540              :       TYPE(section_vals_type), POINTER                   :: input
    1541              : 
    1542          204 :       CALL timeset(routineN, handle)
    1543              :       ! Test for spin unrestricted
    1544          204 :       nspins = SIZE(homo)
    1545              : 
    1546              :       ! Test whether all necessary variables are available
    1547          204 :       my_do_kpoints = .FALSE.
    1548          204 :       IF (.NOT. do_ri_sos_laplace_mp2) THEN
    1549          146 :          my_do_kpoints = do_kpoints_cubic_RPA
    1550              :       END IF
    1551              : 
    1552          146 :       IF (my_do_kpoints) THEN
    1553            4 :          CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
    1554            4 :          E_Range = Emax/Emin
    1555              :       ELSE
    1556          200 :          IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
    1557          152 :             Emin = HUGE(dp)
    1558          152 :             Emax = 0.0_dp
    1559          340 :             DO ispin = 1, nspins
    1560          340 :                IF (homo(ispin) > 0) THEN
    1561          184 :                   Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin))
    1562        14732 :                   Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin)))
    1563              :                END IF
    1564              :             END DO
    1565          152 :             E_Range = Emax/Emin
    1566          152 :             qs_env%mp2_env%e_range = e_range
    1567          152 :             qs_env%mp2_env%e_gap = Emin
    1568              : 
    1569          152 :             CALL get_qs_env(qs_env, input=input)
    1570          152 :             CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=e_range)
    1571          152 :             CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
    1572              :          ELSE
    1573           48 :             E_range = qs_env%mp2_env%E_range
    1574           48 :             Emin = qs_env%mp2_env%E_gap
    1575           48 :             Emax = Emin*E_range
    1576              :          END IF
    1577              :       END IF
    1578              : 
    1579              :       ! When we perform SOS-MP2, we need an additional factor of 2 for the energies (compare with mp2_laplace.F)
    1580              :       ! We do not need weights etc. for the cosine transform
    1581              :       ! We do not scale Emax because it is not needed for SOS-MP2
    1582          204 :       IF (do_ri_sos_laplace_mp2) THEN
    1583           58 :          Emin = Emin*2.0_dp
    1584           58 :          Emax = Emax*2.0_dp
    1585              :       END IF
    1586              : 
    1587          204 :       CALL timestop(handle)
    1588          204 :    END SUBROUTINE
    1589              : 
    1590              : END MODULE mp2_grids
        

Generated by: LCOV version 2.0-1