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