1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculates 2-center integrals for different r12 operators comparing the Solid harmonic 8!> Gaussian integral scheme to the Obara-Saika (OS) scheme 9!> \author Dorothea Golze [05.2016] 10! ************************************************************************************************** 11MODULE shg_integrals_test 12 13 USE basis_set_types, ONLY: allocate_gto_basis_set,& 14 deallocate_gto_basis_set,& 15 gto_basis_set_type,& 16 init_orb_basis_set,& 17 read_gto_basis_set 18 USE constants_operator, ONLY: operator_coulomb,& 19 operator_gauss,& 20 operator_verf,& 21 operator_verfc,& 22 operator_vgauss 23 USE cp_log_handling, ONLY: cp_to_string 24 USE generic_os_integrals, ONLY: int_operators_r12_ab_os,& 25 int_overlap_ab_os,& 26 int_overlap_aba_os,& 27 int_overlap_abb_os,& 28 int_ra2m_ab_os 29 USE generic_shg_integrals, ONLY: int_operators_r12_ab_shg,& 30 int_overlap_ab_shg,& 31 int_overlap_aba_shg,& 32 int_overlap_abb_shg,& 33 int_ra2m_ab_shg 34 USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,& 35 contraction_matrix_shg_mix,& 36 contraction_matrix_shg_rx2m,& 37 get_clebsch_gordon_coefficients 38 USE input_cp2k_subsys, ONLY: create_basis_section 39 USE input_keyword_types, ONLY: keyword_create,& 40 keyword_release,& 41 keyword_type 42 USE input_section_types, ONLY: & 43 section_add_keyword, section_add_subsection, section_create, section_release, & 44 section_type, section_vals_get, section_vals_get_subs_vals, section_vals_type, & 45 section_vals_val_get 46 USE input_val_types, ONLY: real_t 47 USE kinds, ONLY: default_string_length,& 48 dp 49 USE orbital_pointers, ONLY: init_orbital_pointers 50 USE orbital_transformation_matrices, ONLY: init_spherical_harmonics 51#include "./base/base_uses.f90" 52 53 IMPLICIT NONE 54 55 PRIVATE 56 57! ************************************************************************************************** 58 59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'shg_integrals_test' 60 61 PUBLIC :: create_shg_integrals_test_section, shg_integrals_perf_acc_test 62 63CONTAINS 64 65! ************************************************************************************************** 66!> \brief Create input section for unit testing 67!> \param section ... 68! ************************************************************************************************** 69 SUBROUTINE create_shg_integrals_test_section(section) 70 TYPE(section_type), INTENT(INOUT), POINTER :: section 71 72 CHARACTER(len=*), PARAMETER :: routineN = 'create_shg_integrals_test_section', & 73 routineP = moduleN//':'//routineN 74 75 TYPE(keyword_type), POINTER :: keyword 76 TYPE(section_type), POINTER :: subsection 77 78 NULLIFY (keyword, subsection) 79 80 CPASSERT(.NOT. ASSOCIATED(section)) 81 CALL section_create(section, __LOCATION__, name="SHG_INTEGRALS_TEST", & 82 description="Parameters for testing the SHG 2-center integrals for "// & 83 "different r12 operators. Test w.r.t. performance and accurarcy.", & 84 n_keywords=4, n_subsections=1) 85 86 CALL create_basis_section(subsection) 87 CALL section_add_subsection(section, subsection) 88 CALL section_release(subsection) 89 90 CALL keyword_create(keyword, __LOCATION__, & 91 name="_SECTION_PARAMETERS_", & 92 description="Controls the activation the SHG integral test. ", & 93 default_l_val=.FALSE., & 94 lone_keyword_l_val=.TRUE.) 95 CALL section_add_keyword(section, keyword) 96 CALL keyword_release(keyword) 97 98 CALL keyword_create(keyword, __LOCATION__, name="ABC", & 99 description="Specify the lengths of the cell vectors A, B, and C. ", & 100 usage="ABC 10.000 10.000 10.000", unit_str="angstrom", & 101 n_var=3, type_of_var=real_t) 102 CALL section_add_keyword(section, keyword) 103 CALL keyword_release(keyword) 104 105 CALL keyword_create(keyword, __LOCATION__, name="NAB_MIN", & 106 description="Minimum number of atomic distances to consider. ", & 107 default_i_val=8) 108 CALL section_add_keyword(section, keyword) 109 CALL keyword_release(keyword) 110 111 CALL keyword_create(keyword, __LOCATION__, name="NREP", & 112 description="Number of repeated calculation of each integral. ", & 113 default_i_val=1) 114 CALL section_add_keyword(section, keyword) 115 CALL keyword_release(keyword) 116 117 CALL keyword_create(keyword, __LOCATION__, name="CHECK_ACCURACY", & 118 description="Causes abortion when SHG and OS integrals differ "// & 119 "more what's given by ACCURACY_LEVEL.", & 120 default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) 121 CALL section_add_keyword(section, keyword) 122 CALL keyword_release(keyword) 123 124 CALL keyword_create(keyword, __LOCATION__, name="ACCURACY_LEVEL", & 125 description="Level of accuracy for comparison of SHG and OS "// & 126 "integrals.", & 127 default_r_val=1.0E-8_dp) 128 CALL section_add_keyword(section, keyword) 129 CALL keyword_release(keyword) 130 131 CALL keyword_create(keyword, __LOCATION__, name="CALCULATE_DERIVATIVES", & 132 description="Calculates also the derivatives of the integrals.", & 133 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 134 CALL section_add_keyword(section, keyword) 135 CALL keyword_release(keyword) 136 137 CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP", & 138 description="Calculates the integrals (a|b).", & 139 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 140 CALL section_add_keyword(section, keyword) 141 CALL keyword_release(keyword) 142 143 CALL keyword_release(keyword) 144 CALL keyword_create(keyword, __LOCATION__, name="TEST_COULOMB", & 145 description="Calculates the integrals (a|1/r12|b).", & 146 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 147 CALL section_add_keyword(section, keyword) 148 CALL keyword_release(keyword) 149 150 CALL keyword_create(keyword, __LOCATION__, name="TEST_VERF", & 151 description="Calculates the integrals (a|erf(omega*r12)/r12|b).", & 152 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 153 CALL section_add_keyword(section, keyword) 154 CALL keyword_release(keyword) 155 156 CALL keyword_create(keyword, __LOCATION__, name="TEST_VERFC", & 157 description="Calculates the integrals (a|erfc(omega*r12)/r12|b).", & 158 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 159 CALL section_add_keyword(section, keyword) 160 CALL keyword_release(keyword) 161 162 CALL keyword_create(keyword, __LOCATION__, name="TEST_VGAUSS", & 163 description="Calculates the integrals (a|exp(omega*r12^2)/r12|b).", & 164 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 165 CALL section_add_keyword(section, keyword) 166 CALL keyword_release(keyword) 167 168 CALL keyword_create(keyword, __LOCATION__, name="TEST_GAUSS", & 169 description="Calculates the integrals (a|exp(omega*r12^2)|b).", & 170 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 171 CALL section_add_keyword(section, keyword) 172 CALL keyword_release(keyword) 173 174 CALL keyword_create(keyword, __LOCATION__, name="TEST_RA2M", & 175 description="Calculates the integrals (a|(r-Ra)^(2m)|b).", & 176 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 177 CALL section_add_keyword(section, keyword) 178 CALL keyword_release(keyword) 179 180 CALL keyword_create(keyword, __LOCATION__, name="M", & 181 description="Exponent in integral (a|(r-Ra)^(2m)|b).", & 182 default_i_val=1) 183 CALL section_add_keyword(section, keyword) 184 CALL keyword_release(keyword) 185 186 CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP_ABA", & 187 description="Calculates the integrals (a|b|b).", & 188 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 189 CALL section_add_keyword(section, keyword) 190 CALL keyword_release(keyword) 191 192 CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP_ABB", & 193 description="Calculates the integrals (a|b|b).", & 194 default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) 195 CALL section_add_keyword(section, keyword) 196 CALL keyword_release(keyword) 197 198 END SUBROUTINE create_shg_integrals_test_section 199 200! ************************************************************************************************** 201!> \brief Unit test for performance and accuracy of the SHG integrals 202!> \param iw output unit 203!> \param shg_integrals_test_section ... 204! ************************************************************************************************** 205 SUBROUTINE shg_integrals_perf_acc_test(iw, shg_integrals_test_section) 206 INTEGER, INTENT(IN) :: iw 207 TYPE(section_vals_type), INTENT(INOUT), POINTER :: shg_integrals_test_section 208 209 CHARACTER(len=*), PARAMETER :: routineN = 'shg_integrals_perf_acc_test', & 210 routineP = moduleN//':'//routineN 211 212 CHARACTER(LEN=default_string_length) :: basis_type 213 INTEGER :: count_ab, handle, iab, jab, kab, lamax, & 214 lbmax, lcamax, lcbmax, lmax, nab, & 215 nab_min, nab_xyz, nrep, nrep_bas 216 LOGICAL :: acc_check, calc_derivatives, & 217 test_overlap_aba, test_overlap_abb 218 REAL(KIND=dp) :: acc_param 219 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rab 220 REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par 221 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg 222 TYPE(gto_basis_set_type), POINTER :: fba, fbb, oba, obb 223 TYPE(section_vals_type), POINTER :: basis_section 224 225 CALL timeset(routineN, handle) 226 NULLIFY (oba, obb, fba, fbb, basis_section, cell_par) 227 CALL section_vals_val_get(shg_integrals_test_section, "ABC", r_vals=cell_par) 228 CALL section_vals_val_get(shg_integrals_test_section, "NAB_MIN", i_val=nab_min) 229 CALL section_vals_val_get(shg_integrals_test_section, "NREP", i_val=nrep) 230 CALL section_vals_val_get(shg_integrals_test_section, "CHECK_ACCURACY", l_val=acc_check) 231 CALL section_vals_val_get(shg_integrals_test_section, "ACCURACY_LEVEL", r_val=acc_param) 232 CALL section_vals_val_get(shg_integrals_test_section, "CALCULATE_DERIVATIVES", l_val=calc_derivatives) 233 CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", l_val=test_overlap_aba) 234 CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", l_val=test_overlap_abb) 235 236 !*** Read the basis set information 237 basis_section => section_vals_get_subs_vals(shg_integrals_test_section, "BASIS") 238 CALL section_vals_get(basis_section, n_repetition=nrep_bas) 239 IF (.NOT. (nrep_bas == 2 .OR. nrep_bas == 3)) THEN 240 CALL cp_abort(__LOCATION__, & 241 "Provide basis sets") 242 END IF 243 CALL allocate_gto_basis_set(oba) 244 CALL read_gto_basis_set(TRIM("A"), basis_type, oba, basis_section, irep=1) 245 lamax = MAXVAL(oba%lmax) 246 CALL allocate_gto_basis_set(obb) 247 CALL read_gto_basis_set(TRIM("B"), basis_type, obb, basis_section, irep=2) 248 lbmax = MAXVAL(obb%lmax) 249 lmax = MAX(lamax, lbmax) 250 IF (test_overlap_aba) THEN 251 CALL allocate_gto_basis_set(fba) 252 CALL read_gto_basis_set(TRIM("CA"), basis_type, fba, basis_section, irep=3) 253 lcamax = MAXVAL(fba%lmax) 254 lmax = MAX(lamax + lcamax, lbmax) 255 ENDIF 256 IF (test_overlap_abb) THEN 257 CALL allocate_gto_basis_set(fbb) 258 CALL read_gto_basis_set(TRIM("CB"), basis_type, fbb, basis_section, irep=3) 259 lcbmax = MAXVAL(fbb%lmax) 260 lmax = MAX(lamax, lbmax + lcbmax) 261 ENDIF 262 IF (test_overlap_aba .AND. test_overlap_abb) THEN 263 lmax = MAX(MAX(lamax + lcamax, lbmax), MAX(lamax, lbmax + lcbmax)) 264 ENDIF 265 !*** Initialize basis set information 266 CALL init_orbital_pointers(lmax + 1) 267 CALL init_spherical_harmonics(lmax, output_unit=-100) 268 oba%norm_type = 2 269 CALL init_orb_basis_set(oba) 270 obb%norm_type = 2 271 CALL init_orb_basis_set(obb) 272 IF (test_overlap_aba) THEN 273 fba%norm_type = 2 274 CALL init_orb_basis_set(fba) 275 ENDIF 276 IF (test_overlap_abb) THEN 277 fbb%norm_type = 2 278 CALL init_orb_basis_set(fbb) 279 ENDIF 280 ! if shg integrals are later actually used in the code, contraction_matrix_shg should be 281 ! moved to init_orb_basis_set and scon_shg should become an element of gto_basis_set_type 282 CALL contraction_matrix_shg(oba, scona_shg) 283 CALL contraction_matrix_shg(obb, sconb_shg) 284 285 !*** Create range of rab (atomic distances) to be tested 286 nab_xyz = CEILING(REAL(nab_min, KIND=dp)**(1.0_dp/3.0_dp) - 1.0E-06) 287 nab = nab_xyz**3 288 289 ALLOCATE (rab(3, nab)) 290 count_ab = 0 291 DO iab = 1, nab_xyz 292 DO jab = 1, nab_xyz 293 DO kab = 1, nab_xyz 294 count_ab = count_ab + 1 295 rab(:, count_ab) = [iab*ABS(cell_par(1)), jab*ABS(cell_par(2)), kab*ABS(cell_par(3))]/nab_xyz 296 ENDDO 297 ENDDO 298 ENDDO 299 300 !*** Calculate the SHG integrals 301 302 CALL test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 303 shg_integrals_test_section, acc_check, & 304 acc_param, calc_derivatives, iw) 305 306 CALL test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 307 shg_integrals_test_section, acc_check, & 308 acc_param, calc_derivatives, iw) 309 CALL test_shg_ra2m_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 310 shg_integrals_test_section, acc_check, & 311 acc_param, calc_derivatives, iw) 312 313 CALL test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scona_shg, sconb_shg, & 314 shg_integrals_test_section, acc_check, & 315 acc_param, calc_derivatives, iw) 316 317 DEALLOCATE (scona_shg, sconb_shg, rab) 318 CALL deallocate_gto_basis_set(oba) 319 CALL deallocate_gto_basis_set(obb) 320 IF (test_overlap_aba) CALL deallocate_gto_basis_set(fba) 321 IF (test_overlap_abb) CALL deallocate_gto_basis_set(fbb) 322 323 CALL timestop(handle) 324 325 END SUBROUTINE shg_integrals_perf_acc_test 326 327! ************************************************************************************************** 328!> \brief tests two-center integrals of the type [a|O(r12)|b] 329!> \param oba basis set on a 330!> \param obb basis set on b 331!> \param rab distance between a and b 332!> \param nrep ... 333!> \param scona_shg SHG contraction matrix for a 334!> \param sconb_shg SHG contraction matrix for b 335!> \param shg_integrals_test_section ... 336!> \param acc_check if accuracy is checked 337!> \param acc_param accuracy level, if deviation larger abort 338!> \param calc_derivatives ... 339!> \param iw ... 340! ************************************************************************************************** 341 SUBROUTINE test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 342 shg_integrals_test_section, acc_check, & 343 acc_param, calc_derivatives, iw) 344 TYPE(gto_basis_set_type), POINTER :: oba, obb 345 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab 346 INTEGER, INTENT(IN) :: nrep 347 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg 348 TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section 349 LOGICAL, INTENT(IN) :: acc_check 350 REAL(KIND=dp), INTENT(IN) :: acc_param 351 LOGICAL, INTENT(IN) :: calc_derivatives 352 INTEGER, INTENT(IN) :: iw 353 354 CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_operator12_integrals', & 355 routineP = moduleN//':'//routineN 356 357 INTEGER :: iab, irep, nab, nfa, nfb 358 LOGICAL :: test_any, test_coulomb, test_gauss, & 359 test_verf, test_verfc, test_vgauss 360 REAL(KIND=dp) :: ddmax_coulomb, ddmax_gauss, ddmax_verf, ddmax_verfc, ddmax_vgauss, ddtemp, & 361 dmax_coulomb, dmax_gauss, dmax_verf, dmax_verfc, dmax_vgauss, dtemp, omega 362 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vab_os, vab_shg 363 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dvab_os, dvab_shg 364 365 CALL section_vals_val_get(shg_integrals_test_section, "TEST_COULOMB", l_val=test_coulomb) 366 CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERF", l_val=test_verf) 367 CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERFC", l_val=test_verfc) 368 CALL section_vals_val_get(shg_integrals_test_section, "TEST_VGAUSS", l_val=test_vgauss) 369 CALL section_vals_val_get(shg_integrals_test_section, "TEST_GAUSS", l_val=test_gauss) 370 371 test_any = (test_coulomb .OR. test_verf .OR. test_verfc .OR. test_vgauss .OR. test_gauss) 372 373 IF (test_any) THEN 374 nfa = oba%nsgf 375 nfb = obb%nsgf 376 ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3)) 377 ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3)) 378 omega = 2.3_dp 379 dmax_coulomb = 0.0_dp 380 ddmax_coulomb = 0.0_dp 381 dmax_verf = 0.0_dp 382 ddmax_verf = 0.0_dp 383 dmax_verfc = 0.0_dp 384 ddmax_verfc = 0.0_dp 385 dmax_vgauss = 0.0_dp 386 ddmax_vgauss = 0.0_dp 387 dmax_gauss = 0.0_dp 388 ddmax_gauss = 0.0_dp 389 390 nab = SIZE(rab, 2) 391 DO irep = 1, nrep 392 DO iab = 1, nab 393 !*** Coulomb: (a|1/r12|b) 394 IF (test_coulomb) THEN 395 CALL int_operators_r12_ab_shg(operator_coulomb, vab_shg, dvab_shg, rab(:, iab), & 396 oba, obb, scona_shg, sconb_shg, & 397 calculate_forces=calc_derivatives) 398 CALL int_operators_r12_ab_os(operator_coulomb, vab_os, dvab_os, rab(:, iab), & 399 oba, obb, calculate_forces=calc_derivatives) 400 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 401 dmax_coulomb = MAX(dmax_coulomb, dtemp) 402 ddmax_coulomb = MAX(ddmax_coulomb, ddtemp) 403 ENDIF 404 !*** verf: (a|erf(omega*r12)/r12|b) 405 IF (test_verf) THEN 406 CALL int_operators_r12_ab_shg(operator_verf, vab_shg, dvab_shg, rab(:, iab), & 407 oba, obb, scona_shg, sconb_shg, omega, & 408 calc_derivatives) 409 CALL int_operators_r12_ab_os(operator_verf, vab_os, dvab_os, rab(:, iab), & 410 oba, obb, omega=omega, calculate_forces=calc_derivatives) 411 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 412 dmax_verf = MAX(dmax_verf, dtemp) 413 ddmax_verf = MAX(ddmax_verf, ddtemp) 414 ENDIF 415 !*** verfc: (a|erfc(omega*r12)/r12|b) 416 IF (test_verfc) THEN 417 CALL int_operators_r12_ab_shg(operator_verfc, vab_shg, dvab_shg, rab(:, iab), & 418 oba, obb, scona_shg, sconb_shg, omega, & 419 calc_derivatives) 420 CALL int_operators_r12_ab_os(operator_verfc, vab_os, dvab_os, rab(:, iab), & 421 oba, obb, omega=omega, calculate_forces=calc_derivatives) 422 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 423 dmax_verfc = MAX(dmax_verfc, dtemp) 424 ddmax_verfc = MAX(ddmax_verfc, ddtemp) 425 ENDIF 426 !*** vgauss: (a|exp(omega*r12^2)/r12|b) 427 IF (test_vgauss) THEN 428 CALL int_operators_r12_ab_shg(operator_vgauss, vab_shg, dvab_shg, rab(:, iab), & 429 oba, obb, scona_shg, sconb_shg, omega, & 430 calc_derivatives) 431 CALL int_operators_r12_ab_os(operator_vgauss, vab_os, dvab_os, rab(:, iab), & 432 oba, obb, omega=omega, calculate_forces=calc_derivatives) 433 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 434 dmax_vgauss = MAX(dmax_vgauss, dtemp) 435 ddmax_vgauss = MAX(ddmax_vgauss, ddtemp) 436 ENDIF 437 !*** gauss: (a|exp(omega*r12^2)|b) 438 IF (test_gauss) THEN 439 CALL int_operators_r12_ab_shg(operator_gauss, vab_shg, dvab_shg, rab(:, iab), & 440 oba, obb, scona_shg, sconb_shg, omega, & 441 calc_derivatives) 442 CALL int_operators_r12_ab_os(operator_gauss, vab_os, dvab_os, rab(:, iab), & 443 oba, obb, omega=omega, calculate_forces=calc_derivatives) 444 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 445 dmax_gauss = MAX(dmax_gauss, dtemp) 446 ddmax_gauss = MAX(ddmax_gauss, ddtemp) 447 ENDIF 448 ENDDO 449 ENDDO 450 451 IF (iw > 0) THEN 452 WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER SHG and OS INTEGRALS:" 453 WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" 454 IF (test_coulomb) THEN 455 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|1/r12|b]", & 456 dmax_coulomb, ddmax_coulomb 457 ENDIF 458 IF (test_verf) THEN 459 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erf(omega*r12)/r12|b]", & 460 dmax_verf, ddmax_verf 461 ENDIF 462 IF (test_verfc) THEN 463 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erfc(omega*r12)/r12|b]", & 464 dmax_verfc, ddmax_verfc 465 ENDIF 466 IF (test_vgauss) THEN 467 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)/r12|b]", & 468 dmax_vgauss, ddmax_vgauss 469 ENDIF 470 IF (test_gauss) THEN 471 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)|b]", & 472 dmax_gauss, ddmax_gauss 473 ENDIF 474 475 IF (acc_check) THEN 476 IF ((dmax_coulomb >= acc_param) .OR. (ddmax_coulomb >= acc_param)) THEN 477 CPABORT("[a|1/r12|b]: Dev. larger than"//cp_to_string(acc_param)) 478 ENDIF 479 IF ((dmax_verf >= acc_param) .OR. (ddmax_verf >= acc_param)) THEN 480 CPABORT("[a|erf(omega*r12)/r12|b]: Dev. larger than"//cp_to_string(acc_param)) 481 ENDIF 482 IF ((dmax_verfc >= acc_param) .OR. (ddmax_verfc >= acc_param)) THEN 483 CPABORT("[a|erfc(omega*r12)/r12|b]: Dev. larger than"//cp_to_string(acc_param)) 484 ENDIF 485 IF ((dmax_vgauss >= acc_param) .OR. (ddmax_vgauss >= acc_param)) THEN 486 CPABORT("[a|exp(-omega*r12^2)/r12|b]: Dev. larger than"//cp_to_string(acc_param)) 487 ENDIF 488 IF ((dmax_gauss >= acc_param) .OR. (ddmax_gauss >= acc_param)) THEN 489 CPABORT("[a|exp(-omega*r12^2)|b]: Dev. larger than"//cp_to_string(acc_param)) 490 ENDIF 491 ENDIF 492 ENDIF 493 DEALLOCATE (vab_shg, vab_os, dvab_shg, dvab_os) 494 ENDIF 495 496 END SUBROUTINE test_shg_operator12_integrals 497 498! ************************************************************************************************** 499!> \brief tests two center overlap integrals [a|b] 500!> \param oba ... 501!> \param obb ... 502!> \param rab ... 503!> \param nrep ... 504!> \param scona_shg ... 505!> \param sconb_shg ... 506!> \param shg_integrals_test_section ... 507!> \param acc_check ... 508!> \param acc_param ... 509!> \param calc_derivatives ... 510!> \param iw ... 511! ************************************************************************************************** 512 SUBROUTINE test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 513 shg_integrals_test_section, acc_check, & 514 acc_param, calc_derivatives, iw) 515 TYPE(gto_basis_set_type), POINTER :: oba, obb 516 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab 517 INTEGER, INTENT(IN) :: nrep 518 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg 519 TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section 520 LOGICAL, INTENT(IN) :: acc_check 521 REAL(KIND=dp), INTENT(IN) :: acc_param 522 LOGICAL, INTENT(IN) :: calc_derivatives 523 INTEGER, INTENT(IN) :: iw 524 525 CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_overlap_integrals', & 526 routineP = moduleN//':'//routineN 527 528 INTEGER :: iab, irep, nab, nfa, nfb 529 LOGICAL :: test_overlap 530 REAL(KIND=dp) :: ddmax_overlap, ddtemp, dmax_overlap, & 531 dtemp, dummy 532 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab_os, sab_shg 533 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsab_os, dsab_shg 534 535 CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP", & 536 l_val=test_overlap) 537 IF (test_overlap) THEN 538 !effectively switch off screening; makes no sense for the tests 539 oba%set_radius(:) = 1.0E+09_dp 540 obb%set_radius(:) = 1.0E+09_dp 541 oba%pgf_radius(:, :) = 1.0E+09_dp 542 obb%pgf_radius(:, :) = 1.0E+09_dp 543 nfa = oba%nsgf 544 nfb = obb%nsgf 545 dummy = 0.0_dp 546 dmax_overlap = 0.0_dp 547 ddmax_overlap = 0.0_dp 548 ALLOCATE (sab_shg(nfa, nfb), dsab_shg(nfa, nfb, 3)) 549 ALLOCATE (sab_os(nfa, nfb), dsab_os(nfa, nfb, 3)) 550 nab = SIZE(rab, 2) 551 DO irep = 1, nrep 552 DO iab = 1, nab 553 CALL int_overlap_ab_shg(sab_shg, dsab_shg, rab(:, iab), oba, obb, & 554 scona_shg, sconb_shg, calc_derivatives) 555 CALL int_overlap_ab_os(sab_os, dsab_os, rab(:, iab), oba, obb, & 556 calc_derivatives, debug=.FALSE., dmax=dummy) 557 CALL calculate_deviation_ab(sab_shg, sab_os, dsab_shg, dsab_os, dtemp, ddtemp) 558 dmax_overlap = MAX(dmax_overlap, dtemp) 559 ddmax_overlap = MAX(ddmax_overlap, ddtemp) 560 ENDDO 561 ENDDO 562 563 IF (iw > 0) THEN 564 WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER OVERLAP SHG and OS INTEGRALS:" 565 WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" 566 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b]", & 567 dmax_overlap, ddmax_overlap 568 ENDIF 569 IF (acc_check) THEN 570 IF ((dmax_overlap >= acc_param) .OR. (ddmax_overlap >= acc_param)) THEN 571 CPABORT("[a|b]: Deviation larger than"//cp_to_string(acc_param)) 572 ENDIF 573 ENDIF 574 DEALLOCATE (sab_shg, sab_os, dsab_shg, dsab_os) 575 ENDIF 576 577 END SUBROUTINE test_shg_overlap_integrals 578 579! ************************************************************************************************** 580!> \brief tests two-center integrals of the type [a|(r-Ra)^(2m)|b] 581!> \param oba ... 582!> \param obb ... 583!> \param rab ... 584!> \param nrep ... 585!> \param scona_shg ... 586!> \param sconb_shg ... 587!> \param shg_integrals_test_section ... 588!> \param acc_check ... 589!> \param acc_param ... 590!> \param calc_derivatives ... 591!> \param iw ... 592! ************************************************************************************************** 593 SUBROUTINE test_shg_ra2m_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, & 594 shg_integrals_test_section, acc_check, & 595 acc_param, calc_derivatives, iw) 596 TYPE(gto_basis_set_type), POINTER :: oba, obb 597 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab 598 INTEGER, INTENT(IN) :: nrep 599 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg 600 TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section 601 LOGICAL, INTENT(IN) :: acc_check 602 REAL(KIND=dp), INTENT(IN) :: acc_param 603 LOGICAL, INTENT(IN) :: calc_derivatives 604 INTEGER, INTENT(IN) :: iw 605 606 CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_ra2m_integrals', & 607 routineP = moduleN//':'//routineN 608 609 INTEGER :: iab, irep, m, nab, nfa, nfb 610 LOGICAL :: test_ra2m 611 REAL(KIND=dp) :: ddmax, ddtemp, dmax, dtemp 612 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vab_os, vab_shg 613 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dvab_os, dvab_shg 614 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: scon_ra2m 615 616 CALL section_vals_val_get(shg_integrals_test_section, "TEST_RA2M", & 617 l_val=test_ra2m) 618 IF (test_ra2m) THEN 619 CALL section_vals_val_get(shg_integrals_test_section, "M", & 620 i_val=m) 621 nfa = oba%nsgf 622 nfb = obb%nsgf 623 dmax = 0.0_dp 624 ddmax = 0.0_dp 625 CALL contraction_matrix_shg_rx2m(oba, m, scona_shg, scon_ra2m) 626 ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3)) 627 ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3)) 628 nab = SIZE(rab, 2) 629 DO irep = 1, nrep 630 DO iab = 1, nab 631 CALL int_ra2m_ab_shg(vab_shg, dvab_shg, rab(:, iab), oba, obb, & 632 scon_ra2m, sconb_shg, m, calc_derivatives) 633 CALL int_ra2m_ab_os(vab_os, dvab_os, rab(:, iab), oba, obb, m, calc_derivatives) 634 CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp) 635 dmax = MAX(dmax, dtemp) 636 ddmax = MAX(ddmax, ddtemp) 637 ENDDO 638 ENDDO 639 IF (iw > 0) THEN 640 WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER RA2m SHG and OS INTEGRALS:" 641 WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" 642 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|(r-Ra)^(2m)|b]", & 643 dmax, ddmax 644 ENDIF 645 IF (acc_check) THEN 646 IF ((dmax >= acc_param) .OR. (ddmax >= acc_param)) THEN 647 CPABORT("[a|ra^(2m)|b]: Deviation larger than"//cp_to_string(acc_param)) 648 ENDIF 649 ENDIF 650 DEALLOCATE (scon_ra2m) 651 DEALLOCATE (vab_shg, vab_os, dvab_shg, dvab_os) 652 ENDIF 653 END SUBROUTINE test_shg_ra2m_integrals 654 655! ************************************************************************************************** 656!> \brief test overlap integrals [a|b|a] and [a|b|b] 657!> \param oba ... 658!> \param obb ... 659!> \param fba ... 660!> \param fbb ... 661!> \param rab ... 662!> \param nrep ... 663!> \param scon_oba ... 664!> \param scon_obb ... 665!> \param shg_integrals_test_section ... 666!> \param acc_check ... 667!> \param acc_param ... 668!> \param calc_derivatives ... 669!> \param iw ... 670! ************************************************************************************************** 671 SUBROUTINE test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scon_oba, scon_obb, & 672 shg_integrals_test_section, acc_check, & 673 acc_param, calc_derivatives, iw) 674 TYPE(gto_basis_set_type), POINTER :: oba, obb, fba, fbb 675 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab 676 INTEGER, INTENT(IN) :: nrep 677 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_oba, scon_obb 678 TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section 679 LOGICAL, INTENT(IN) :: acc_check 680 REAL(KIND=dp), INTENT(IN) :: acc_param 681 LOGICAL, INTENT(IN) :: calc_derivatives 682 INTEGER, INTENT(IN) :: iw 683 684 CHARACTER(len=*), PARAMETER :: routineN = 'test_shg_overlap_aba_integrals', & 685 routineP = moduleN//':'//routineN 686 687 INTEGER :: iab, irep, la_max, lb_max, lbb_max, & 688 maxl_orb, maxl_ri, nab, nba, nbb, nfa, & 689 nfb 690 INTEGER, DIMENSION(:, :), POINTER :: ncg_none0 691 INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list, fba_index, fbb_index, & 692 oba_index, obb_index 693 LOGICAL :: test_overlap_aba, test_overlap_abb 694 REAL(KIND=dp) :: ddmax_overlap_aba, ddmax_overlap_abb, & 695 ddtemp, dmax_overlap_aba, & 696 dmax_overlap_abb, dtemp, dummy 697 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: saba_os, saba_shg, sabb_os, sabb_shg 698 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dsaba_os, dsaba_shg, dsabb_os, dsabb_shg 699 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cg_coeff 700 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix 701 702 CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", & 703 l_val=test_overlap_aba) 704 CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", & 705 l_val=test_overlap_abb) 706 IF (test_overlap_aba .OR. test_overlap_abb) THEN 707 !effectively switch off screening; makes no sense for the tests 708 oba%set_radius(:) = 1.0E+09_dp 709 obb%set_radius(:) = 1.0E+09_dp 710 oba%pgf_radius(:, :) = 1.0E+09_dp 711 obb%pgf_radius(:, :) = 1.0E+09_dp 712 nba = oba%nsgf 713 nbb = obb%nsgf 714 maxl_orb = MAX(MAXVAL(oba%lmax), MAXVAL(obb%lmax)) 715 la_max = MAXVAL(oba%lmax) 716 lb_max = MAXVAL(obb%lmax) 717 IF (test_overlap_aba) THEN 718 fba%set_radius(:) = 1.0E+09_dp 719 fba%pgf_radius(:, :) = 1.0E+09_dp 720 nfa = fba%nsgf 721 maxl_ri = MAXVAL(fba%lmax) + 1 ! + 1 to avoid fail for l=0 722 ALLOCATE (saba_shg(nba, nbb, nfa), dsaba_shg(nba, nbb, nfa, 3)) 723 ALLOCATE (saba_os(nba, nbb, nfa), dsaba_os(nba, nbb, nfa, 3)) 724 CALL contraction_matrix_shg_mix(oba, fba, oba_index, fba_index, scona_mix) 725 ENDIF 726 IF (test_overlap_abb) THEN 727 fbb%set_radius(:) = 1.0E+09_dp 728 fbb%pgf_radius(:, :) = 1.0E+09_dp 729 nfb = fbb%nsgf 730 maxl_ri = MAXVAL(fbb%lmax) + 1 731 lbb_max = MAXVAL(obb%lmax) + MAXVAL(fbb%lmax) 732 ALLOCATE (sabb_shg(nba, nbb, nfb), dsabb_shg(nba, nbb, nfb, 3)) 733 ALLOCATE (sabb_os(nba, nbb, nfb), dsabb_os(nba, nbb, nfb, 3)) 734 CALL contraction_matrix_shg_mix(obb, fbb, obb_index, fbb_index, sconb_mix) 735 ENDIF 736 dummy = 0.0_dp 737 dmax_overlap_aba = 0.0_dp 738 ddmax_overlap_aba = 0.0_dp 739 dmax_overlap_abb = 0.0_dp 740 ddmax_overlap_abb = 0.0_dp 741 CALL get_clebsch_gordon_coefficients(cg_coeff, cg_none0_list, ncg_none0, maxl_orb, maxl_ri) 742 nab = SIZE(rab, 2) 743 IF (test_overlap_aba) THEN 744 DO irep = 1, nrep 745 DO iab = 1, nab 746 CALL int_overlap_aba_shg(saba_shg, dsaba_shg, rab(:, iab), oba, obb, fba, & 747 scon_obb, scona_mix, oba_index, fba_index, & 748 cg_coeff, cg_none0_list, ncg_none0, & 749 calc_derivatives) 750 CALL int_overlap_aba_os(saba_os, dsaba_os, rab(:, iab), oba, obb, fba, & 751 calc_derivatives, debug=.FALSE., dmax=dummy) 752 CALL calculate_deviation_abx(saba_shg, saba_os, dsaba_shg, dsaba_os, dtemp, ddtemp) 753 dmax_overlap_aba = MAX(dmax_overlap_aba, dtemp) 754 ddmax_overlap_aba = MAX(ddmax_overlap_aba, ddtemp) 755 ENDDO 756 ENDDO 757 DEALLOCATE (oba_index, fba_index, scona_mix) 758 DEALLOCATE (saba_shg, saba_os, dsaba_shg, dsaba_os) 759 ENDIF 760 IF (test_overlap_abb) THEN 761 DO irep = 1, nrep 762 DO iab = 1, nab 763 CALL int_overlap_abb_shg(sabb_shg, dsabb_shg, rab(:, iab), oba, obb, fbb, & 764 scon_oba, sconb_mix, obb_index, fbb_index, & 765 cg_coeff, cg_none0_list, ncg_none0, & 766 calc_derivatives) 767 CALL int_overlap_abb_os(sabb_os, dsabb_os, rab(:, iab), oba, obb, fbb, & 768 calc_derivatives, debug=.FALSE., dmax=dummy) 769 CALL calculate_deviation_abx(sabb_shg, sabb_os, dsabb_shg, dsabb_os, dtemp, ddtemp) 770 dmax_overlap_abb = MAX(dmax_overlap_abb, dtemp) 771 ddmax_overlap_abb = MAX(ddmax_overlap_abb, ddtemp) 772 ENDDO 773 ENDDO 774 DEALLOCATE (obb_index, fbb_index, sconb_mix) 775 DEALLOCATE (sabb_shg, sabb_os, dsabb_shg, dsabb_os) 776 ENDIF 777 IF (iw > 0) THEN 778 WRITE (iw, FMT="(/,T2,A)") "TEST INFO [a|b|x] OVERLAP SHG and OS INTEGRALS:" 779 WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives" 780 IF (test_overlap_aba) THEN 781 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|a]", & 782 dmax_overlap_aba, ddmax_overlap_aba 783 ENDIF 784 IF (test_overlap_abb) THEN 785 WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|b]", & 786 dmax_overlap_abb, ddmax_overlap_abb 787 ENDIF 788 ENDIF 789 IF (acc_check) THEN 790 IF ((dmax_overlap_aba >= acc_param) .OR. (ddmax_overlap_aba >= acc_param)) THEN 791 CPABORT("[a|b|a]: Dev. larger than"//cp_to_string(acc_param)) 792 ENDIF 793 IF ((dmax_overlap_abb >= acc_param) .OR. (ddmax_overlap_abb >= acc_param)) THEN 794 CPABORT("[a|b|b]: Dev. larger than"//cp_to_string(acc_param)) 795 ENDIF 796 ENDIF 797 DEALLOCATE (cg_coeff, cg_none0_list, ncg_none0) 798 ENDIF 799 800 END SUBROUTINE test_shg_overlap_aba_integrals 801 802! ************************************************************************************************** 803!> \brief Calculation of the deviation between SHG and OS integrals 804!> \param vab_shg integral matrix obtained from the SHG scheme 805!> \param vab_os integral matrix obtained from the OS scheme 806!> \param dvab_shg derivative of the integrals, SHG 807!> \param dvab_os derivative of the integrals, OS 808!> \param dmax maximal deviation of vab matrices 809!> \param ddmax maximal deviation of dvab matrices 810! ************************************************************************************************** 811 SUBROUTINE calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax) 812 813 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vab_shg, vab_os 814 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dvab_shg, dvab_os 815 REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax 816 817 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation_ab', & 818 routineP = moduleN//':'//routineN 819 820 INTEGER :: i, j, k 821 REAL(KIND=dp) :: diff 822 823 dmax = 0.0_dp 824 ddmax = 0.0_dp 825 826 ! integrals vab 827 DO j = 1, SIZE(vab_shg, 2) 828 DO i = 1, SIZE(vab_shg, 1) 829 diff = ABS(vab_shg(i, j) - vab_os(i, j)) 830 dmax = MAX(dmax, diff) 831 ENDDO 832 ENDDO 833 834 ! derivatives dvab 835 DO k = 1, 3 836 DO j = 1, SIZE(dvab_shg, 2) 837 DO i = 1, SIZE(dvab_shg, 1) 838 diff = ABS(dvab_shg(i, j, k) - dvab_os(i, j, k)) 839 ddmax = MAX(ddmax, diff) 840 ENDDO 841 ENDDO 842 ENDDO 843 844 END SUBROUTINE calculate_deviation_ab 845 846! ************************************************************************************************** 847!> \brief Calculation of the deviation between SHG and OS integrals 848!> \param vab_shg integral matrix obtained from the SHG scheme 849!> \param vab_os integral matrix obtained from the OS scheme 850!> \param dvab_shg derivative of the integrals, SHG 851!> \param dvab_os derivative of the integrals, OS 852!> \param dmax maximal deviation of vab matrices 853!> \param ddmax maximal deviation of dvab matrices 854! ************************************************************************************************** 855 SUBROUTINE calculate_deviation_abx(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax) 856 857 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: vab_shg, vab_os 858 REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dvab_shg, dvab_os 859 REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax 860 861 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation_abx', & 862 routineP = moduleN//':'//routineN 863 864 INTEGER :: i, j, k, l 865 REAL(KIND=dp) :: diff 866 867 dmax = 0.0_dp 868 ddmax = 0.0_dp 869 870 ! integrals vab 871 DO k = 1, SIZE(vab_shg, 3) 872 DO j = 1, SIZE(vab_shg, 2) 873 DO i = 1, SIZE(vab_shg, 1) 874 diff = ABS(vab_shg(i, j, k) - vab_os(i, j, k)) 875 dmax = MAX(dmax, diff) 876 ENDDO 877 ENDDO 878 ENDDO 879 880 ! derivatives dvab 881 DO l = 1, 3 882 DO k = 1, SIZE(dvab_shg, 3) 883 DO j = 1, SIZE(dvab_shg, 2) 884 DO i = 1, SIZE(dvab_shg, 1) 885 diff = ABS(dvab_shg(i, j, k, l) - dvab_os(i, j, k, l)) 886 ddmax = MAX(ddmax, diff) 887 ENDDO 888 ENDDO 889 ENDDO 890 ENDDO 891 892 END SUBROUTINE calculate_deviation_abx 893 894END MODULE shg_integrals_test 895