LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_test.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 62.4 % 141 88
Test Date: 2025-07-25 12:55:17 Functions: 66.7 % 3 2

            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              : !> \brief Methods for testing / debugging.
       9              : !> \par History
      10              : !>       2015 09 created
      11              : !> \author Patrick Seewald
      12              : ! **************************************************************************************************
      13              : 
      14              : MODULE eri_mme_test
      15              : 
      16              :    USE eri_mme_gaussian,                ONLY: create_gaussian_overlap_dist_to_hermite,&
      17              :                                               create_hermite_to_cartesian
      18              :    USE eri_mme_integrate,               ONLY: eri_mme_2c_integrate,&
      19              :                                               eri_mme_3c_integrate
      20              :    USE eri_mme_types,                   ONLY: eri_mme_coulomb,&
      21              :                                               eri_mme_longrange,&
      22              :                                               eri_mme_param,&
      23              :                                               eri_mme_set_potential,&
      24              :                                               eri_mme_yukawa
      25              :    USE kinds,                           ONLY: dp
      26              :    USE mathconstants,                   ONLY: twopi
      27              :    USE message_passing,                 ONLY: mp_para_env_type
      28              :    USE orbital_pointers,                ONLY: ncoset
      29              : #include "../base/base_uses.f90"
      30              : 
      31              :    IMPLICIT NONE
      32              : 
      33              :    PRIVATE
      34              : 
      35              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      36              : 
      37              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_test'
      38              : 
      39              :    PUBLIC :: eri_mme_2c_perf_acc_test, &
      40              :              eri_mme_3c_perf_acc_test
      41              : 
      42              : CONTAINS
      43              : ! **************************************************************************************************
      44              : !> \brief Unit test for performance and accuracy
      45              : !> \param param ...
      46              : !> \param l_max ...
      47              : !> \param zet ...
      48              : !> \param rabc ...
      49              : !> \param nrep ...
      50              : !> \param test_accuracy ...
      51              : !> \param para_env ...
      52              : !> \param iw ...
      53              : !> \param potential ...
      54              : !> \param pot_par ...
      55              : !> \param G_count ...
      56              : !> \param R_count ...
      57              : ! **************************************************************************************************
      58            8 :    SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, &
      59              :                                        para_env, iw, potential, pot_par, G_count, R_count)
      60              :       TYPE(eri_mme_param), INTENT(INOUT)                 :: param
      61              :       INTEGER, INTENT(IN)                                :: l_max
      62              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
      63              :          INTENT(IN)                                      :: zet
      64              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: rabc
      65              :       INTEGER, INTENT(IN)                                :: nrep
      66              :       LOGICAL, INTENT(INOUT)                             :: test_accuracy
      67              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
      68              :       INTEGER, INTENT(IN)                                :: iw
      69              :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
      70              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
      71              :       INTEGER, INTENT(OUT), OPTIONAL                     :: G_count, R_count
      72              : 
      73              :       INTEGER                                            :: iab, irep, izet, l, nR, nzet
      74              :       LOGICAL                                            :: acc_check
      75              :       REAL(KIND=dp)                                      :: acc, t0, t1
      76            8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: time
      77            8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: I_diff, I_ref, I_test
      78              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
      79              : 
      80            8 :       IF (PRESENT(G_count)) G_count = 0
      81            8 :       IF (PRESENT(R_count)) R_count = 0
      82              : 
      83            8 :       nzet = SIZE(zet)
      84            8 :       nR = SIZE(rabc, 2)
      85              : 
      86            8 :       IF (PRESENT(potential)) THEN
      87            8 :          CALL eri_mme_set_potential(param, potential, pot_par)
      88              :       END IF
      89              : 
      90              :       ! Calculate reference values (Exact expression in G space converged to high precision)
      91            8 :       IF (test_accuracy) THEN
      92           78 :          ht = twopi*TRANSPOSE(param%h_inv)
      93              : 
      94           36 :          ALLOCATE (I_ref(ncoset(l_max), ncoset(l_max), nR, nzet))
      95       613086 :          I_ref(:, :, :, :) = 0.0_dp
      96              : 
      97           30 :          DO izet = 1, nzet
      98          222 :             DO iab = 1, nR
      99              :                CALL eri_mme_2c_integrate(param, 0, l_max, 0, l_max, zet(izet), zet(izet), rabc(:, iab), &
     100              :                                          I_ref(:, :, iab, izet), 0, 0, &
     101          216 :                                          normalize=.TRUE., potential=potential, pot_par=pot_par)
     102              : 
     103              :             END DO
     104              :          END DO
     105              :       END IF
     106              : 
     107              :       ! test performance and accuracy of MME method
     108           48 :       ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), nR, nzet))
     109           40 :       ALLOCATE (I_diff(ncoset(l_max), ncoset(l_max), nR, nzet))
     110              : 
     111           32 :       ALLOCATE (time(0:l_max, nzet))
     112           56 :       DO l = 0, l_max
     113          320 :          DO izet = 1, nzet
     114          264 :             CALL CPU_TIME(t0)
     115          528 :             DO irep = 1, nrep
     116         2640 :                DO iab = 1, nR
     117              :                   CALL eri_mme_2c_integrate(param, 0, l, 0, l, zet(izet), zet(izet), rabc(:, iab), &
     118              :                                             I_test(:, :, iab, izet), 0, 0, &
     119              :                                             G_count=G_count, R_count=R_count, &
     120         2376 :                                             normalize=.TRUE.)
     121              :                END DO
     122              :             END DO
     123          264 :             CALL CPU_TIME(t1)
     124          312 :             time(l, izet) = t1 - t0
     125              :          END DO
     126              :       END DO
     127              : 
     128            8 :       CALL para_env%sum(time)
     129              : 
     130            8 :       IF (test_accuracy) THEN
     131       613086 :          I_diff(:, :, :, :) = ABS(I_test - I_ref)
     132              :       END IF
     133              : 
     134            8 :       IF (iw > 0) THEN
     135            4 :          WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time"
     136            4 :          WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy"
     137              : 
     138           28 :          DO l = 0, l_max
     139          160 :             DO izet = 1, nzet
     140          132 :                IF (test_accuracy) THEN
     141        83976 :                   acc = MAXVAL(I_diff(ncoset(l - 1) + 1:ncoset(l), ncoset(l - 1) + 1:ncoset(l), :, izet))
     142              :                ELSE
     143           60 :                   acc = 0.0_dp
     144              :                END IF
     145              : 
     146              :                WRITE (iw, '(T11, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
     147          156 :                   l, zet(izet), time(l, izet)/nrep, acc
     148              :             END DO
     149              :          END DO
     150              : 
     151            4 :          IF (test_accuracy) THEN
     152            3 :             WRITE (iw, '(/T2, A, 47X, ES9.2)') "ERI_MME| Maximum error:", &
     153       306546 :                MAXVAL(I_diff)
     154              : 
     155            3 :             IF (param%is_ortho) THEN
     156       306543 :                acc_check = param%err_mm + param%err_c .GE. MAXVAL(I_diff)
     157              :             ELSE
     158              :                acc_check = .TRUE.
     159              :             END IF
     160              : 
     161            3 :             IF (.NOT. acc_check) &
     162            0 :                CPABORT("Actual error greater than upper bound estimate.")
     163              : 
     164              :          END IF
     165              :       END IF
     166              : 
     167            8 :    END SUBROUTINE eri_mme_2c_perf_acc_test
     168              : 
     169              : ! **************************************************************************************************
     170              : !> \brief ...
     171              : !> \param param ...
     172              : !> \param l_max ...
     173              : !> \param zet ...
     174              : !> \param rabc ...
     175              : !> \param nrep ...
     176              : !> \param nsample ...
     177              : !> \param para_env ...
     178              : !> \param iw ...
     179              : !> \param potential ...
     180              : !> \param pot_par ...
     181              : !> \param GG_count ...
     182              : !> \param GR_count ...
     183              : !> \param RR_count ...
     184              : ! **************************************************************************************************
     185            2 :    SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, &
     186              :                                        para_env, iw, potential, pot_par, GG_count, GR_count, RR_count)
     187              :       TYPE(eri_mme_param), INTENT(INOUT)                 :: param
     188              :       INTEGER, INTENT(IN)                                :: l_max
     189              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     190              :          INTENT(IN)                                      :: zet
     191              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     192              :          INTENT(IN)                                      :: rabc
     193              :       INTEGER, INTENT(IN)                                :: nrep
     194              :       INTEGER, INTENT(IN), OPTIONAL                      :: nsample
     195              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     196              :       INTEGER, INTENT(IN)                                :: iw
     197              :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
     198              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
     199              :       INTEGER, INTENT(OUT), OPTIONAL                     :: GG_count, GR_count, RR_count
     200              : 
     201              :       INTEGER                                            :: ira, irb, irc, irep, ixyz, izeta, izetb, &
     202              :                                                             izetc, la, lb, lc, nintg, nR, ns, nzet
     203              :       REAL(KIND=dp)                                      :: t0, t1, time
     204            2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: I_test
     205              : 
     206            2 :       IF (PRESENT(GG_count)) GG_count = 0
     207            2 :       IF (PRESENT(GR_count)) GR_count = 0
     208            2 :       IF (PRESENT(RR_count)) RR_count = 0
     209              : 
     210            2 :       ns = 1
     211            2 :       IF (PRESENT(nsample)) ns = nsample
     212              : 
     213            2 :       nzet = SIZE(zet)
     214            2 :       nR = SIZE(rabc, 2)
     215              : 
     216            2 :       IF (PRESENT(potential)) THEN
     217            2 :          CALL eri_mme_set_potential(param, potential, pot_par)
     218              :       END IF
     219              : 
     220            2 :       IF (param%debug) THEN
     221            0 :          DO izeta = 1, nzet
     222            0 :          DO izetb = 1, nzet
     223            0 :             DO ira = 1, nR
     224            0 :             DO irb = 1, nR
     225            0 :                DO ixyz = 1, 3
     226              :                   CALL overlap_dist_expansion_test(l_max, l_max, zet(izeta), zet(izetb), &
     227            0 :                                                    rabc(ixyz, ira), rabc(ixyz, irb), 0.0_dp, param%debug_delta)
     228              :                END DO
     229              :             END DO
     230              :             END DO
     231              :          END DO
     232              :          END DO
     233              :       END IF
     234              : 
     235            2 :       IF (iw > 0) THEN
     236            1 :          IF (PRESENT(potential)) THEN
     237            2 :             SELECT CASE (potential)
     238              :             CASE (eri_mme_coulomb)
     239            1 :                WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
     240              :             CASE (eri_mme_yukawa)
     241            0 :                WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par
     242              :             CASE (eri_mme_longrange)
     243            1 :                WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par
     244              :             END SELECT
     245              :          ELSE
     246            0 :             WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
     247              :          END IF
     248            1 :          WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time"
     249            1 :          WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time"
     250              :       END IF
     251              : 
     252           10 :       ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), ncoset(l_max)))
     253              : 
     254            2 :       nintg = 0
     255           14 :       DO la = 0, l_max
     256           86 :       DO lb = 0, l_max
     257          516 :       DO lc = 0, l_max
     258         4824 :          DO izeta = 1, nzet
     259        47952 :          DO izetb = 1, nzet
     260       479520 :          DO izetc = 1, nzet
     261       432000 :             nintg = nintg + 1
     262       475200 :             IF (MOD(nintg, ns) .EQ. 0) THEN
     263      2145708 :                I_test(:, :, :) = 0.0_dp
     264           12 :                CALL CPU_TIME(t0)
     265           24 :                DO irep = 1, nrep
     266          120 :                   DO ira = 1, nR
     267          876 :                   DO irb = 1, nR
     268         7008 :                   DO irc = 1, nR
     269              :                      CALL eri_mme_3c_integrate(param, 0, la, 0, lb, 0, lc, zet(izeta), zet(izetb), zet(izetc), &
     270              :                                                rabc(:, ira), rabc(:, irb), rabc(:, irc), I_test, 0, 0, 0, &
     271         6912 :                                                GG_count, GR_count, RR_count)
     272              :                   END DO
     273              :                   END DO
     274              :                   END DO
     275              :                END DO
     276           12 :                CALL CPU_TIME(t1)
     277           12 :                time = t1 - t0
     278           12 :                CALL para_env%sum(time)
     279           12 :                IF (iw > 0) THEN
     280              :                   WRITE (iw, '(T11, I1, 1X, I1, 1X, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
     281            6 :                      la, lb, lc, zet(izeta), zet(izetb), zet(izetc), time/nrep
     282              :                END IF
     283              :             END IF
     284              :          END DO
     285              :          END DO
     286              :          END DO
     287              :       END DO
     288              :       END DO
     289              :       END DO
     290              : 
     291            2 :    END SUBROUTINE eri_mme_3c_perf_acc_test
     292              : 
     293              : ! **************************************************************************************************
     294              : !> \brief check that expanding an overlap distribution of cartesian/hermite Gaussians into a
     295              : !>        lin combi of single cartesian/hermite Gaussians is correct.
     296              : !> \param l_max ...
     297              : !> \param m_max ...
     298              : !> \param zeta ...
     299              : !> \param zetb ...
     300              : !> \param R1 ...
     301              : !> \param R2 ...
     302              : !> \param r ...
     303              : !> \param tolerance ...
     304              : !> \note STATUS: tested
     305              : ! **************************************************************************************************
     306            0 :    SUBROUTINE overlap_dist_expansion_test(l_max, m_max, zeta, zetb, R1, R2, r, tolerance)
     307              :       INTEGER, INTENT(IN)                                :: l_max, m_max
     308              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, R1, R2, r, tolerance
     309              : 
     310              :       INTEGER                                            :: l, m, t
     311              :       REAL(KIND=dp)                                      :: C_prod_err, H_prod_err, Rp, zetp
     312            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: C1, C2, C_ol, H1, H2, H_ol
     313            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: C_prod_ref, C_prod_test, H_prod_ref, &
     314            0 :                                                             H_prod_test, h_to_c_1, h_to_c_2, &
     315            0 :                                                             h_to_c_ol
     316            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: E_C, E_H
     317              : 
     318            0 :       zetp = zeta + zetb
     319            0 :       Rp = (zeta*R1 + zetb*R2)/zetp
     320            0 :       ALLOCATE (C1(0:l_max), H1(0:l_max))
     321            0 :       ALLOCATE (C2(0:m_max), H2(0:m_max))
     322            0 :       ALLOCATE (C_ol(0:l_max + m_max))
     323            0 :       ALLOCATE (H_ol(0:l_max + m_max))
     324            0 :       ALLOCATE (C_prod_ref(0:l_max, 0:m_max))
     325            0 :       ALLOCATE (C_prod_test(0:l_max, 0:m_max))
     326            0 :       ALLOCATE (H_prod_ref(0:l_max, 0:m_max))
     327            0 :       ALLOCATE (H_prod_test(0:l_max, 0:m_max))
     328              : 
     329            0 :       ALLOCATE (E_C(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
     330            0 :       ALLOCATE (E_H(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
     331            0 :       CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 1, E_C)
     332            0 :       CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 2, E_H)
     333            0 :       CALL create_hermite_to_cartesian(zetp, l_max + m_max, h_to_c_ol)
     334              :       CALL create_hermite_to_cartesian(zeta, l_max, h_to_c_1)
     335              :       CALL create_hermite_to_cartesian(zetb, m_max, h_to_c_2)
     336              : 
     337            0 :       DO t = 0, l_max + m_max
     338            0 :          C_ol(t) = (r - Rp)**t*EXP(-zetp*(r - Rp)**2)
     339              :       END DO
     340              : 
     341            0 :       DO l = 0, l_max
     342            0 :          C1(l) = (r - R1)**l*EXP(-zeta*(r - R1)**2)
     343              :       END DO
     344            0 :       DO m = 0, m_max
     345            0 :          C2(m) = (r - R2)**m*EXP(-zetb*(r - R2)**2)
     346              :       END DO
     347              : 
     348            0 :       H1(:) = MATMUL(TRANSPOSE(h_to_c_1(0:, 0:)), C1)
     349            0 :       H2(:) = MATMUL(TRANSPOSE(h_to_c_2(0:, 0:)), C2)
     350            0 :       H_ol(:) = MATMUL(TRANSPOSE(h_to_c_ol(0:, 0:)), C_ol)
     351              : 
     352            0 :       DO m = 0, m_max
     353            0 :          DO l = 0, l_max
     354            0 :             C_prod_ref(l, m) = C1(l)*C2(m)
     355            0 :             H_prod_ref(l, m) = H1(l)*H2(m)
     356            0 :             C_prod_test(l, m) = 0.0_dp
     357            0 :             H_prod_test(l, m) = 0.0_dp
     358            0 :             DO t = 0, l + m
     359            0 :                C_prod_test(l, m) = C_prod_test(l, m) + E_C(t, l, m)*H_ol(t)
     360            0 :                H_prod_test(l, m) = H_prod_test(l, m) + E_H(t, l, m)*H_ol(t)
     361              :             END DO
     362              :          END DO
     363              :       END DO
     364              : 
     365            0 :       C_prod_err = MAXVAL(ABS(C_prod_test - C_prod_ref)/(0.5_dp*(ABS(C_prod_test) + ABS(C_prod_ref)) + 1.0_dp))
     366            0 :       H_prod_err = MAXVAL(ABS(H_prod_test - H_prod_ref)/(0.5_dp*(ABS(C_prod_test) + ABS(C_prod_ref)) + 1.0_dp))
     367              : 
     368            0 :       CPASSERT(C_prod_err .LE. tolerance)
     369            0 :       CPASSERT(H_prod_err .LE. tolerance)
     370              :       MARK_USED(tolerance)
     371              : 
     372            0 :    END SUBROUTINE overlap_dist_expansion_test
     373              : 
     374              : END MODULE eri_mme_test
        

Generated by: LCOV version 2.0-1