Line data    Source code

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