From 469a92f17a8139dfc2ca74c99abe8776dca5fd0e Mon Sep 17 00:00:00 2001 From: Patrick Seewald Date: Tue, 14 May 2019 12:54:09 +0200 Subject: [PATCH] MME integrals: implement new potentials --- src/cp_eri_mme_interface.F | 135 +++++++---- src/eri_mme/eri_mme_error_control.F | 191 ++++++++++----- src/eri_mme/eri_mme_gaussian.F | 71 +++++- src/eri_mme/eri_mme_integrate.F | 18 +- src/eri_mme/eri_mme_lattice_summation.F | 28 ++- src/eri_mme/eri_mme_test.F | 52 +++- src/eri_mme/eri_mme_types.F | 228 +++++++++++++----- src/mp2_eri.F | 34 +-- src/qmmm_init.F | 6 +- .../test_eri_mme_accuracy_longrange.inp | 27 +++ .../LIBTEST/test_eri_mme_accuracy_yukawa.inp | 27 +++ 11 files changed, 598 insertions(+), 219 deletions(-) create mode 100644 tests/LIBTEST/test_eri_mme_accuracy_longrange.inp create mode 100644 tests/LIBTEST/test_eri_mme_accuracy_yukawa.inp diff --git a/src/cp_eri_mme_interface.F b/src/cp_eri_mme_interface.F index cca82cd3c7..2524a0541a 100644 --- a/src/cp_eri_mme_interface.F +++ b/src/cp_eri_mme_interface.F @@ -23,11 +23,14 @@ MODULE cp_eri_mme_interface USE cp_para_types, ONLY: cp_para_env_type USE eri_mme_test, ONLY: eri_mme_2c_perf_acc_test,& eri_mme_3c_perf_acc_test - USE eri_mme_types, ONLY: eri_mme_init,& + USE eri_mme_types, ONLY: eri_mme_coulomb,& + eri_mme_init,& + eri_mme_longrange,& eri_mme_param,& eri_mme_print_grid_info,& eri_mme_release,& - eri_mme_set_params + eri_mme_set_params,& + eri_mme_yukawa USE input_keyword_types, ONLY: keyword_create,& keyword_release,& keyword_type @@ -46,6 +49,7 @@ MODULE cp_eri_mme_interface USE orbital_pointers, ONLY: init_orbital_pointers USE qs_kind_types, ONLY: get_qs_kind,& qs_kind_type + USE string_utilities, ONLY: s2a #include "./base/base_uses.f90" IMPLICIT NONE @@ -79,10 +83,6 @@ MODULE cp_eri_mme_interface INTEGER :: G_count_2c, R_count_2c INTEGER :: GG_count_3c, GR_count_3c, RR_count_3c REAL(KIND=dp), DIMENSION(3, 3) :: hmat - REAL(KIND=dp) :: zet_err_minimax - REAL(KIND=dp) :: zet_err_cutoff - INTEGER :: l_err_cutoff - INTEGER :: l_max LOGICAL :: do_calib END TYPE cp_eri_mme_param @@ -346,8 +346,11 @@ CONTAINS !> \param basis_type_1 ... !> \param basis_type_2 ... !> \param para_env ... +!> \param potential ... +!> \param pot_par ... ! ************************************************************************************************** - SUBROUTINE eri_mme_set_params_from_basis(param, cell, qs_kind_set, basis_type_1, basis_type_2, para_env) + SUBROUTINE eri_mme_set_params_from_basis(param, cell, qs_kind_set, basis_type_1, basis_type_2, para_env, & + potential, pot_par) TYPE(cp_eri_mme_param), INTENT(INOUT) :: param TYPE(cell_type), INTENT(IN) :: cell TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), & @@ -355,21 +358,24 @@ CONTAINS CHARACTER(len=*), INTENT(IN) :: basis_type_1 CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type_2 TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_set_params_from_basis', & routineP = moduleN//':'//routineN - INTEGER :: handle, l_err_cutoff, l_max - REAL(KIND=dp) :: zet_err_cutoff, zet_err_minimax + INTEGER :: handle, l_max, l_max_zet + REAL(KIND=dp) :: zet_max, zet_min CALL timeset(routineN, handle) CALL error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, & - zet_err_minimax, zet_err_cutoff, l_err_cutoff, l_max) + zet_min, zet_max, l_max_zet, l_max) CALL eri_mme_set_params_custom(param, cell%hmat, cell%orthorhombic, & - zet_err_minimax, zet_err_cutoff, l_err_cutoff, & - l_max, para_env) + zet_min, zet_max, l_max_zet, & + l_max, para_env, & + potential, pot_par) CALL timestop(handle) END SUBROUTINE eri_mme_set_params_from_basis @@ -379,19 +385,24 @@ CONTAINS !> \param param ... !> \param hmat ... !> \param is_ortho ... -!> \param zet_err_minimax ... -!> \param zet_err_cutoff ... -!> \param l_err_cutoff ... +!> \param zet_min ... +!> \param zet_max ... +!> \param l_max_zet ... !> \param l_max ... !> \param para_env ... +!> \param potential ... +!> \param pot_par ... ! ************************************************************************************************** - SUBROUTINE eri_mme_set_params_custom(param, hmat, is_ortho, zet_err_minimax, zet_err_cutoff, l_err_cutoff, l_max, para_env) + SUBROUTINE eri_mme_set_params_custom(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, & + potential, pot_par) TYPE(cp_eri_mme_param), INTENT(INOUT) :: param REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat LOGICAL, INTENT(IN) :: is_ortho - REAL(KIND=dp), INTENT(IN) :: zet_err_minimax, zet_err_cutoff - INTEGER, INTENT(IN) :: l_err_cutoff, l_max + REAL(KIND=dp), INTENT(IN) :: zet_min, zet_max + INTEGER, INTENT(IN) :: l_max_zet, l_max TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par REAL(KIND=dp), PARAMETER :: eps_changed = 1.0E-14_dp @@ -401,10 +412,10 @@ CONTAINS ELSE ! only calibrate cutoff if parameters (cell, basis coefficients) have changed IF (ALL(ABS(param%hmat-hmat) < eps_changed) .AND. & - ABS(param%zet_err_minimax-zet_err_minimax) < eps_changed .AND. & - ABS(param%zet_err_cutoff-zet_err_cutoff) < eps_changed .AND. & - param%l_err_cutoff == l_err_cutoff .AND. & - param%l_max == l_max) THEN + ABS(param%par%zet_min-zet_min) < eps_changed .AND. & + ABS(param%par%zet_max-zet_max) < eps_changed .AND. & + param%par%l_max_zet == l_max_zet .AND. & + param%par%l_max == l_max) THEN param%par%do_calib_cutoff = .FALSE. ELSE param%par%do_calib_cutoff = .TRUE. @@ -414,14 +425,10 @@ CONTAINS param%par%do_calib_cutoff = .FALSE. ENDIF - CALL eri_mme_set_params(param%par, hmat, is_ortho, zet_err_minimax, zet_err_cutoff, l_err_cutoff, l_max, para_env) + CALL eri_mme_set_params(param%par, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, & + potential, pot_par) param%hmat(:, :) = hmat(:, :) - param%zet_err_minimax = zet_err_minimax - param%zet_err_cutoff = zet_err_cutoff - param%l_err_cutoff = l_err_cutoff - - param%l_max = l_max CALL eri_mme_print_info(param) END SUBROUTINE eri_mme_set_params_custom @@ -431,20 +438,20 @@ CONTAINS !> \param qs_kind_set ... !> \param basis_type_1 ... !> \param basis_type_2 ... -!> \param zet_mm Smallest exponent, used to estimate error due to minimax approx. -!> \param zet_c contains max. exponent and the largest exponent for max. l, -!> used to estimate cutoff error -!> \param l_c contains the largest l for max. exponent and max. l, -!> used to estimate cutoff error +!> \param zet_min Smallest exponent, used to estimate error due to minimax approx. +!> \param zet_max contains max. exponent, +!> used to estimate cutoff error +!> \param l_max_zet contains the largest l for max. exponent, +!> used to estimate cutoff error !> \param l_max ... ! ************************************************************************************************** - SUBROUTINE error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, zet_mm, zet_c, l_c, l_max) + SUBROUTINE error_est_pgf_params_from_basis(qs_kind_set, basis_type_1, basis_type_2, zet_min, zet_max, l_max_zet, l_max) TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), & POINTER :: qs_kind_set CHARACTER(len=*), INTENT(IN) :: basis_type_1 CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type_2 - REAL(KIND=dp), INTENT(OUT) :: zet_mm, zet_c - INTEGER, INTENT(OUT) :: l_c, l_max + REAL(KIND=dp), INTENT(OUT) :: zet_min, zet_max + INTEGER, INTENT(OUT) :: l_max_zet, l_max CHARACTER(LEN=*), PARAMETER :: routineN = 'error_est_pgf_params_from_basis', & routineP = moduleN//':'//routineN @@ -461,7 +468,7 @@ CONTAINS l_m = 0 zet_m = 0.0_dp l_zet = -1 - zet_mm = -1.0_dp + zet_min = -1.0_dp nkind = SIZE(qs_kind_set) nbasis = MERGE(2, 1, PRESENT(basis_type_2)) @@ -483,10 +490,10 @@ CONTAINS l_m = MAX(l_m, MAXVAL(basis_set%lmax(:))) DO iset = 1, nset zet_m = MAX(zet_m, MAXVAL(basis_set%zet(1:npgf(iset), iset))) - IF (zet_mm .LT. 0.0_dp) THEN - zet_mm = MINVAL(basis_set%zet(1:npgf(iset), iset)) + IF (zet_min .LT. 0.0_dp) THEN + zet_min = MINVAL(basis_set%zet(1:npgf(iset), iset)) ELSE - zet_mm = MIN(zet_mm, MINVAL(basis_set%zet(1:npgf(iset), iset))) + zet_min = MIN(zet_min, MINVAL(basis_set%zet(1:npgf(iset), iset))) ENDIF ENDDO ENDDO @@ -515,11 +522,11 @@ CONTAINS CPASSERT(l_zet .GE. 0) - zet_c = zet_m + zet_max = zet_m ! l + 1 because we may calculate forces ! this is probably a safe choice also for the case that forces are not needed - l_c = l_zet+1 + l_max_zet = l_zet+1 l_max = l_m+1 CALL timestop(handle) @@ -541,6 +548,18 @@ CONTAINS logger => param%logger unit_nr = param%par%unit_nr + + IF (unit_nr > 0) THEN + SELECT CASE (param%par%potential) + CASE (eri_mme_coulomb) + WRITE (unit_nr, '(/T2, A)') "ERI_MME| Potential: Coulomb" + CASE (eri_mme_yukawa) + WRITE (unit_nr, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", param%par%pot_par + CASE (eri_mme_longrange) + WRITE (unit_nr, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", param%par%pot_par + END SELECT + ENDIF + IF (unit_nr > 0) THEN WRITE (unit_nr, '(/T2, A, T71, F10.1)') "ERI_MME| Cutoff for ERIs [a.u.]:", param%par%cutoff WRITE (unit_nr, '(/T2, A, T78, I3/)') "ERI_MME| Number of terms in minimax approximation:", param%par%n_minimax @@ -674,6 +693,22 @@ CONTAINS default_i_val=1) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="POTENTIAL", & + description="Operator to test", & + default_i_val=eri_mme_coulomb, & + enum_i_vals=(/eri_mme_coulomb, eri_mme_yukawa, eri_mme_longrange/), & + enum_c_vals=s2a("COULOMB", "YUKAWA", "LONGRANGE"), & + enum_desc=s2a("1/r", "exp(-a*r)/r", "erf(a*r)/r")) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="POTENTIAL_PARAM", & + description="Parameter 'a' for chosen potential, ignored for Coulomb", & + default_r_val=0.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_eri_mme_test_section ! ************************************************************************************************** @@ -736,9 +771,10 @@ CONTAINS INTEGER :: count_r, G_count, GG_count, GR_count, i, & ix, iy, iz, l_max, min_nR, nR, nR_xyz, & - nrep, nsample, nzet, R_count, RR_count + nrep, nsample, nzet, potential, & + R_count, RR_count LOGICAL :: test_2c, test_3c, test_accuracy - REAL(KIND=dp) :: zet_fac, zetmax, zetmin + REAL(KIND=dp) :: pot_par, zet_fac, zetmax, zetmin REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zet REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rabc REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par @@ -762,6 +798,8 @@ CONTAINS CALL section_vals_val_get(eri_mme_test_section, "ZET_MAX", r_val=zetmax) CALL section_vals_val_get(eri_mme_test_section, "NZET", i_val=nzet) CALL section_vals_val_get(eri_mme_test_section, "NSAMPLE_3C", i_val=nsample) + CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL", i_val=potential) + CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL_PARAM", r_val=pot_par) IF (nzet .LE. 0) & CPABORT("Number of exponents NZET must be greater than 0.") @@ -807,16 +845,19 @@ CONTAINS ENDDO ! initialize MME method - CALL cp_eri_mme_set_params(param, box%hmat, box%orthorhombic, MINVAL(zet), MAXVAL(zet), l_max, l_max, para_env) + CALL cp_eri_mme_set_params(param, box%hmat, box%orthorhombic, MINVAL(zet), MAXVAL(zet), l_max, l_max, para_env, & + potential, pot_par) IF (iw > 0) WRITE (iw, '(T2, A, T61, I20)') "ERI_MME| Number of atomic distances:", nR G_count = 0; R_count = 0 GG_count = 0; GR_count = 0; RR_count = 0 - IF (test_2c) CALL eri_mme_2c_perf_acc_test(param%par, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, G_count, R_count) + IF (test_2c) CALL eri_mme_2c_perf_acc_test(param%par, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, & + potential=potential, pot_par=pot_par, G_count=G_count, R_count=R_count) IF (test_3c) CALL eri_mme_3c_perf_acc_test(param%par, l_max, zet, rabc, nrep, nsample, & - para_env, iw, GG_count, GR_count, RR_count) + para_env, iw, potential=potential, pot_par=pot_par, & + GG_count=GG_count, GR_count=GR_count, RR_count=RR_count) CALL cp_eri_mme_update_local_counts(param, para_env, G_count, R_count, GG_count, GR_count, RR_count) CALL cp_eri_mme_finalize(param) CALL cell_release(box) diff --git a/src/eri_mme/eri_mme_error_control.F b/src/eri_mme/eri_mme_error_control.F index b85cea6977..6dbf010299 100644 --- a/src/eri_mme/eri_mme_error_control.F +++ b/src/eri_mme/eri_mme_error_control.F @@ -31,7 +31,7 @@ MODULE eri_mme_error_control CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_error_control' - PUBLIC :: calibrate_cutoff, cutoff_minimax_error + PUBLIC :: calibrate_cutoff, cutoff_minimax_error, minimax_error, cutoff_error CONTAINS ! ************************************************************************************************** @@ -41,10 +41,10 @@ CONTAINS !> \param h_inv ... !> \param G_min ... !> \param vol ... -!> \param zet_mm Exponent to estimate minimax error -!> \param l_mm Total ang. mom. quantum number to estimate minimax error -!> \param zet_c Max. exponents to estimate cutoff error -!> \param l_c Max. total ang. mom. quantum numbers to estimate cutoff error +!> \param zet_min Minimum exponent +!> \param l_mm Total ang. mom. quantum number +!> \param zet_max Max. exponents to estimate cutoff error +!> \param l_max_zet Max. total ang. mom. quantum numbers to estimate cutoff error !> \param n_minimax Number of terms in minimax approximation !> \param cutoff_l Initial guess of lower bound for cutoff !> \param cutoff_r Initial guess of upper bound for cutoff @@ -59,16 +59,16 @@ CONTAINS !> \param print_calib ... !> \param unit_nr ... ! ************************************************************************************************** - SUBROUTINE calibrate_cutoff(hmat, h_inv, G_min, vol, zet_mm, l_mm, zet_c, l_c, & + SUBROUTINE calibrate_cutoff(hmat, h_inv, G_min, vol, zet_min, l_mm, zet_max, l_max_zet, & n_minimax, cutoff_l, cutoff_r, tol, delta, & cutoff, err_mm, err_c, C_mm, para_env, print_calib, unit_nr) REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv REAL(KIND=dp), INTENT(IN) :: G_min REAL(KIND=dp) :: vol - REAL(KIND=dp), INTENT(IN) :: zet_mm + REAL(KIND=dp), INTENT(IN) :: zet_min INTEGER, INTENT(IN) :: l_mm - REAL(KIND=dp), INTENT(IN) :: zet_c - INTEGER, INTENT(IN) :: l_c, n_minimax + REAL(KIND=dp), INTENT(IN) :: zet_max + INTEGER, INTENT(IN) :: l_max_zet, n_minimax REAL(KIND=dp), INTENT(IN) :: cutoff_l, cutoff_r, tol, delta REAL(KIND=dp), INTENT(OUT) :: cutoff, err_mm, err_c, C_mm TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env @@ -87,9 +87,9 @@ CONTAINS do_print = unit_nr > 0 .AND. print_calib IF (do_print) THEN WRITE (unit_nr, '(/T2, A)') "ERI_MME| Basis set parameters for estimating minimax error" - WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_mm, l_mm + WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_min, l_mm WRITE (unit_nr, '(T2, A)') "ERI_MME| Basis set parameters for estimating cutoff error" - WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_c, l_c + WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME| exp, l:", zet_max, l_max_zet ENDIF max_iter = 100 @@ -124,7 +124,7 @@ CONTAINS ! approx.) is hit DO i = 1, 2 - CALL cutoff_minimax_error(cutoff_lr(i), hmat, h_inv, vol, G_min, zet_mm, l_mm, zet_c, l_c, & + CALL cutoff_minimax_error(cutoff_lr(i), hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, & n_minimax, minimax_aw, delta_mm(i), delta_c(i), C_mm, para_env) ENDDO @@ -151,7 +151,7 @@ CONTAINS "Maximum number of iterations in bisection to determine cutoff has been exceeded") cutoff_mid = 0.5_dp*(cutoff_lr(1)+cutoff_lr(2)) - CALL cutoff_minimax_error(cutoff_mid, hmat, h_inv, vol, G_min, zet_mm, l_mm, zet_c, l_c, & + CALL cutoff_minimax_error(cutoff_mid, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, & n_minimax, minimax_aw, delta_mm_mid, delta_c_mid, C_mm, para_env) IF (do_print) WRITE (unit_nr, '(T11, I2, F11.1, F11.1, F11.1, 3X, ES9.2, 3X, ES9.2, 3X, ES9.2)') & iter2, cutoff_lr(1), cutoff_lr(2), cutoff_mid, & @@ -189,10 +189,10 @@ CONTAINS !> \param h_inv ... !> \param vol ... !> \param G_min ... -!> \param zet_mm Exponent of P to estimate minimax error +!> \param zet_min Exponent of P to estimate minimax error !> \param l_mm total ang. mom. quantum number of P to estimate minimax error -!> \param zet_ctff Max. exponents of P to estimate cutoff error -!> \param l_ctff Max. total ang. mom. quantum numbers of P to estimate cutoff error +!> \param zet_max Max. exponents of P to estimate cutoff error +!> \param l_max_zet Max. total ang. mom. quantum numbers of P to estimate cutoff error !> \param n_minimax Number of terms in minimax approximation !> \param minimax_aw Minimax coefficients !> \param err_mm Minimax error @@ -201,14 +201,14 @@ CONTAINS !> minimax approx. !> \param para_env ... ! ************************************************************************************************** - SUBROUTINE cutoff_minimax_error(cutoff, hmat, h_inv, vol, G_min, zet_mm, l_mm, zet_ctff, l_ctff, & + SUBROUTINE cutoff_minimax_error(cutoff, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, & n_minimax, minimax_aw, err_mm, err_ctff, C_mm, para_env) REAL(KIND=dp), INTENT(IN) :: cutoff REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv - REAL(KIND=dp), INTENT(IN) :: vol, G_min, zet_mm + REAL(KIND=dp), INTENT(IN) :: vol, G_min, zet_min INTEGER, INTENT(IN) :: l_mm - REAL(KIND=dp), INTENT(IN) :: zet_ctff - INTEGER, INTENT(IN) :: l_ctff, n_minimax + REAL(KIND=dp), INTENT(IN) :: zet_max + INTEGER, INTENT(IN) :: l_max_zet, n_minimax REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw REAL(KIND=dp), INTENT(OUT) :: err_mm, err_ctff, C_mm TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env @@ -216,30 +216,103 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'cutoff_minimax_error', & routineP = moduleN//':'//routineN - INTEGER :: i_aw, i_xyz, iG, iter, max_iter, nG - REAL(KIND=dp) :: C, dG, E_mm, eps_zet, err0, err1, err_c, err_ctff_curr, err_ctff_prev, & - err_d, G, G_1, G_c, gr, prod_mm_k, zet_a, zet_b, zet_c, zet_ctff_tmp, zet_d, zet_div + REAL(KIND=dp) :: delta_mm + + CALL minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, & + n_minimax, minimax_aw, err_mm, delta_mm) + CALL cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, & + n_minimax, minimax_aw, err_ctff, C_mm, para_env) + + END SUBROUTINE cutoff_minimax_error + +! ************************************************************************************************** +!> \brief Minimax error, simple analytical formula +!> Note minimax error may blow up for small exponents. This is also observed numerically, +!> but in this case, error estimate is no upper bound. +!> \param cutoff ... +!> \param hmat ... +!> \param vol ... +!> \param G_min ... +!> \param zet_min Exponent of P to estimate minimax error +!> \param l_mm total ang. mom. quantum number of P to estimate minimax error +!> \param n_minimax Number of terms in minimax approximation +!> \param minimax_aw Minimax coefficients +!> \param err_mm Minimax error +!> \param delta_mm ... +!> \param potential ... +!> \param pot_par ... +! ************************************************************************************************** + SUBROUTINE minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, & + n_minimax, minimax_aw, err_mm, delta_mm, potential, pot_par) + REAL(KIND=dp), INTENT(IN) :: cutoff + REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat + REAL(KIND=dp), INTENT(IN) :: vol, G_min, zet_min + INTEGER, INTENT(IN) :: l_mm, n_minimax + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw + REAL(KIND=dp), INTENT(OUT) :: err_mm, delta_mm + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par + + CHARACTER(LEN=*), PARAMETER :: routineN = 'minimax_error', routineP = moduleN//':'//routineN + + INTEGER :: i_xyz + REAL(KIND=dp) :: prod_mm_k + + CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw(:), & + potential=potential, pot_par=pot_par, err_minimax=delta_mm) + + prod_mm_k = 1.0_dp + DO i_xyz = 1, 3 + prod_mm_k = prod_mm_k*(ABS(hmat(i_xyz, i_xyz))/twopi+ & + MERGE(SQRT(2.0_dp/(zet_min*pi))*EXP(-1.0_dp), 0.0_dp, l_mm .GT. 0)) + ENDDO + err_mm = 32*pi**4/vol*delta_mm*prod_mm_k + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an +!> upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C). +!> +!> Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1. +!> The error is calculated for all l up to l_max and golden section search algorithm is +!> applied to find the exponent that maximizes cutoff error. +!> \param cutoff ... +!> \param h_inv ... +!> \param G_min ... +!> \param zet_max Max. exponents of P to estimate cutoff error +!> \param l_max_zet Max. total ang. mom. quantum numbers of P to estimate cutoff error +!> \param n_minimax Number of terms in minimax approximation +!> \param minimax_aw Minimax coefficients +!> \param err_ctff Cutoff error +!> \param C_mm Scaling constant to generalize AM-GM upper bound estimate to +!> minimax approx. +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, & + n_minimax, minimax_aw, err_ctff, C_mm, para_env) + REAL(KIND=dp), INTENT(IN) :: cutoff + REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv + REAL(KIND=dp), INTENT(IN) :: G_min, zet_max + INTEGER, INTENT(IN) :: l_max_zet, n_minimax + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw + REAL(KIND=dp), INTENT(OUT) :: err_ctff, C_mm + TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'cutoff_error', routineP = moduleN//':'//routineN + + INTEGER :: i_aw, iG, iter, max_iter, nG + REAL(KIND=dp) :: C, dG, eps_zet, err0, err1, err_c, err_ctff_curr, err_ctff_prev, err_d, G, & + G_1, G_c, gr, zet_a, zet_b, zet_c, zet_d, zet_div, zet_max_tmp ! parameters for finding exponent maximizing cutoff error eps_zet = 1.0E-05_dp ! tolerance for exponent zet_div = 2.0_dp ! sampling constant for finding initial values of exponents max_iter = 100 ! maximum number of iterations in golden section search - - zet_ctff_tmp = zet_ctff - - CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw(:), E_mm) G_c = SQRT(2.0*cutoff) - ! 1) Minimax error, simple analytical formula - ! Note minimax error may blow up for small exponents. This is also observed numerically, - ! but in this case, error estimate is no upper bound. - prod_mm_k = 1.0_dp - DO i_xyz = 1, 3 - prod_mm_k = prod_mm_k*(ABS(hmat(i_xyz, i_xyz))/twopi+ & - MERGE(SQRT(2.0_dp/(zet_mm*pi))*EXP(-1.0_dp), 0.0_dp, l_mm .GT. 0)) - ENDDO - err_mm = 32*pi**4/vol*E_mm*prod_mm_k + zet_max_tmp = zet_max ! 2) Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an ! upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C). @@ -277,15 +350,15 @@ CONTAINS CALL cp_abort(__LOCATION__, "Maximum number of iterations for finding "// & "exponent maximizing cutoff error has been exceeded.") - CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_ctff, zet_ctff_tmp, C, err_ctff_curr, para_env) + CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max_tmp, C, err_ctff_curr, para_env) IF (err_ctff_prev .GE. err_ctff_curr) THEN - zet_a = zet_ctff_tmp - zet_b = MIN(zet_ctff_tmp*zet_div**2, zet_ctff) + zet_a = zet_max_tmp + zet_b = MIN(zet_max_tmp*zet_div**2, zet_max) EXIT ELSE err_ctff_prev = err_ctff_curr ENDIF - zet_ctff_tmp = zet_ctff_tmp/zet_div + zet_max_tmp = zet_max_tmp/zet_div ENDDO ! Golden section search @@ -293,13 +366,13 @@ CONTAINS zet_d = zet_a+gr*(zet_b-zet_a) DO iter = 1, max_iter+1 IF (ABS(zet_c-zet_d) .LT. eps_zet*(zet_a+zet_b)) THEN - CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_ctff, zet_a, C, err0, para_env) - CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_ctff, zet_b, C, err1, para_env) + CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_a, C, err0, para_env) + CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_b, C, err1, para_env) err_ctff_curr = MAX(err0, err1) EXIT ENDIF - CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_ctff, zet_c, C, err_c, para_env) - CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_ctff, zet_d, C, err_d, para_env) + CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_c, C, err_c, para_env) + CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_d, C, err_d, para_env) IF (err_c .GT. err_d) THEN zet_b = zet_d zet_d = zet_c @@ -312,28 +385,28 @@ CONTAINS ENDDO err_ctff = err_ctff_curr - END SUBROUTINE cutoff_minimax_error + END SUBROUTINE ! ************************************************************************************************** !> \brief Calculate cutoff error estimate by using C_mm/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 !> as an upper bound for 1/G^2 (and its minimax approximation) for |G| > G_c. -!> Error is referring to a basis function P with fixed exponent zet_c and -!> max. angular momentum l_c. +!> Error is referring to a basis function P with fixed exponent zet_max and +!> max. angular momentum l_max_zet. !> \param cutoff ... !> \param h_inv ... !> \param G_min ... -!> \param l_c ... -!> \param zet_c ... +!> \param l_max_zet ... +!> \param zet_max ... !> \param C_mm ... !> \param err_c ... !> \param para_env ... ! ************************************************************************************************** - SUBROUTINE cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_c, zet_c, C_mm, err_c, para_env) + SUBROUTINE cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max, C_mm, err_c, para_env) REAL(KIND=dp), INTENT(IN) :: cutoff REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv REAL(KIND=dp), INTENT(IN) :: G_min - INTEGER, INTENT(IN) :: l_c - REAL(KIND=dp), INTENT(IN) :: zet_c, C_mm + INTEGER, INTENT(IN) :: l_max_zet + REAL(KIND=dp), INTENT(IN) :: zet_max, C_mm REAL(KIND=dp), INTENT(OUT) :: err_c TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env @@ -353,13 +426,13 @@ CONTAINS G_res = 0.5_dp*G_min ! resolution for screening err_c = 0.0_dp - alpha_G = 1.0_dp/(2.0_dp*zet_c) - prefactor = 1.0_dp/zet_c + alpha_G = 1.0_dp/(2.0_dp*zet_max) + prefactor = 1.0_dp/zet_max - ALLOCATE (S_G_l(0:2*l_c, 3)) - ALLOCATE (S_G_u(0:2*l_c, 3)) + ALLOCATE (S_G_l(0:2*l_max_zet, 3)) + ALLOCATE (S_G_u(0:2*l_max_zet, 3)) - G_rad = exp_radius(2*l_c, alpha_G, eps_G, prefactor, epsin=G_res) + G_rad = exp_radius(2*l_max_zet, alpha_G, eps_G, prefactor, epsin=G_res) ! Parallelization of sum over G vectors my_p = para_env%mepos ! mpi rank @@ -419,7 +492,7 @@ CONTAINS S_G_u = S_G_u*2.0_dp ! to include negative values of G - DO l = 0, l_c + DO l = 0, l_max_zet DO ax = 0, l DO ay = 0, l-ax az = l-ax-ay @@ -434,7 +507,7 @@ CONTAINS S_G_l(2*ax, 1)*S_G_u(2*ay, 2)*S_G_l(2*az, 3)+ & S_G_l(2*ax, 1)*S_G_l(2*ay, 2)*S_G_u(2*az, 3) - err_c_l = 4.0_dp*pi**4*hermite_gauss_norm(zet_c, [ax, ay, az])**2* & + err_c_l = 4.0_dp*pi**4*hermite_gauss_norm(zet_max, [ax, ay, az])**2* & C_mm/3.0_dp*sum_G_diff err_c = MAX(err_c, err_c_l) diff --git a/src/eri_mme/eri_mme_gaussian.F b/src/eri_mme/eri_mme_gaussian.F index 126ef0074f..ac825ea8be 100644 --- a/src/eri_mme/eri_mme_gaussian.F +++ b/src/eri_mme/eri_mme_gaussian.F @@ -24,6 +24,10 @@ MODULE eri_mme_gaussian CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_gaussian' + INTEGER, PARAMETER, PUBLIC :: eri_mme_coulomb = 1, & + eri_mme_yukawa = 2, & + eri_mme_longrange = 3 + PUBLIC :: & create_gaussian_overlap_dist_to_hermite, & create_hermite_to_cartesian, & @@ -83,18 +87,41 @@ CONTAINS !> \param cutoff Plane Wave cutoff !> \param G_min Minimum absolute value of G !> \param minimax_aw Minimax coefficients a_i, w_i +!> \param potential potential to use. Accepts the following values: +!> 1: coulomb potential V(r)=1/r +!> 2: yukawa potential V(r)=e(-a*r)/r +!> 3: long-range coulomb erf(a*r)/r +!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r !> \param err_minimax Maximum error MAX (|1/G^2-\sum_i w_i exp(-a_i G^2)|) ! ************************************************************************************************** - SUBROUTINE get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw, err_minimax) + SUBROUTINE get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw, potential, pot_par, err_minimax) INTEGER, INTENT(IN) :: n_minimax REAL(KIND=dp), INTENT(IN) :: cutoff, G_min REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: minimax_aw + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par REAL(KIND=dp), INTENT(OUT), OPTIONAL :: err_minimax CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_coeff_v_gspace', & routineP = moduleN//':'//routineN + INTEGER :: potential_prv REAL(KIND=dp) :: dG, G_max, minimax_Rc + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a, w + + IF (PRESENT(potential)) THEN + potential_prv = potential + ELSE + potential_prv = eri_mme_coulomb + ENDIF + + IF (potential_prv > 3) THEN + CPABORT("unknown potential") + ENDIF + + IF ((potential_prv >= 2) .AND. .NOT. PRESENT(pot_par)) THEN + CPABORT("potential parameter pot_par required for yukawa or long-range Coulomb") + ENDIF dG = 1.0E-3 ! Resolution in G to determine error of minimax approximation @@ -103,16 +130,50 @@ CONTAINS ! Minimax approx. needs to be valid in range [G_min, G_max] ! 1) compute minimax coefficients - G_max = SQRT(3.0_dp*2.0_dp*cutoff) - minimax_Rc = (G_max/G_min)**2 G_max = SQRT(3.0_dp*2.0_dp*cutoff) CPASSERT(G_max .GT. G_min) + IF (potential_prv == eri_mme_coulomb .OR. potential_prv == eri_mme_longrange) THEN + minimax_Rc = (G_max/G_min)**2 + ELSEIF (potential_prv == eri_mme_yukawa) THEN + minimax_Rc = (G_max**2+pot_par**2)/(G_min**2+pot_par**2) + ENDIF CALL get_exp_minimax_coeff(n_minimax, minimax_Rc, minimax_aw, err_minimax) - minimax_aw = minimax_aw/G_min**2 - IF (PRESENT(err_minimax)) err_minimax = err_minimax/G_min**2 + ALLOCATE (a(n_minimax)); ALLOCATE (w(n_minimax)) + a(:) = minimax_aw(:n_minimax) + w(:) = minimax_aw(n_minimax+1:) + SELECT CASE (potential_prv) + ! Scale minimax coefficients to incorporate different Fourier transforms + CASE (eri_mme_coulomb) + ! FT = 1/G**2 + a(:) = a/G_min**2 + w(:) = w/G_min**2 + CASE (eri_mme_yukawa) + ! FT = 1/(G**2 + pot_par**2) + w(:) = w*EXP((-a*pot_par**2)/(G_min**2+pot_par**2))/(G_min**2+pot_par**2) + a(:) = a/(G_min**2+pot_par**2) + CASE (eri_mme_longrange) + ! FT = exp(-(G/pot_par)**2)/G**2 + ! approximating 1/G**2 as for Coulomb: + a(:) = a/G_min**2 + w(:) = w/G_min**2 + ! incorporate exponential factor: + a(:) = a+1.0_dp/pot_par**2 + END SELECT + minimax_aw = [a(:), w(:)] + + IF (PRESENT(err_minimax)) THEN + IF (potential_prv == eri_mme_coulomb) THEN + err_minimax = err_minimax/G_min**2 + ELSEIF (potential_prv == eri_mme_yukawa) THEN + err_minimax = err_minimax/(G_min**2+pot_par**2) + ELSEIF (potential_prv == eri_mme_longrange) THEN + err_minimax = err_minimax/G_min**2 ! approx. of Coulomb + err_minimax = err_minimax*EXP(-G_min**2/pot_par**2) ! exponential factor + ENDIF + ENDIF END SUBROUTINE get_minimax_coeff_v_gspace diff --git a/src/eri_mme/eri_mme_integrate.F b/src/eri_mme/eri_mme_integrate.F index f2765070d9..5b039e6cf1 100644 --- a/src/eri_mme/eri_mme_integrate.F +++ b/src/eri_mme/eri_mme_integrate.F @@ -61,12 +61,12 @@ CONTAINS !> \param o2 ... !> \param G_count ... !> \param R_count ... -!> \param normalize calculate integrals w.r.t. normalized Hermite-Gaussians -!> \param exact_method create reference values by exact method (no minimax approx. -!> and sum is converged to much higher precision) +!> \param normalize calculate integrals w.r.t. normalized Hermite-Gaussians +!> \param potential use exact potential instead of minimax approx. (for testing only) +!> \param pot_par potential parameter ! ************************************************************************************************** SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, & - hab, o1, o2, G_count, R_count, normalize, exact_method) + hab, o1, o2, G_count, R_count, normalize, potential, pot_par) TYPE(eri_mme_param), INTENT(IN) :: param INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max REAL(KIND=dp), INTENT(IN) :: zeta, zetb @@ -74,7 +74,9 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: hab INTEGER, INTENT(IN) :: o1, o2 INTEGER, INTENT(INOUT), OPTIONAL :: G_count, R_count - LOGICAL, INTENT(IN), OPTIONAL :: normalize, exact_method + LOGICAL, INTENT(IN), OPTIONAL :: normalize + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_2c_integrate', & routineP = moduleN//':'//routineN @@ -117,8 +119,8 @@ CONTAINS l_max = la_max+lb_max - IF (PRESENT(exact_method)) THEN - exact = exact_method + IF (PRESENT(potential)) THEN + exact = .TRUE. ELSE exact = .FALSE. ENDIF @@ -151,7 +153,7 @@ CONTAINS G_rad = exp_radius(l_max, alpha_G, 0.01*param%sum_precision, 1.0_dp, epsin=G_res) G_bounds(:) = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi)) - CALL pgf_sum_2c_gspace_3d(S_G_no, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol, coulomb=.TRUE.) + CALL pgf_sum_2c_gspace_3d(S_G_no, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol, potential, pot_par) ELSE DO i_aw = 1, n_aw diff --git a/src/eri_mme/eri_mme_lattice_summation.F b/src/eri_mme/eri_mme_lattice_summation.F index 1884c0b404..61068a31f2 100644 --- a/src/eri_mme/eri_mme_lattice_summation.F +++ b/src/eri_mme/eri_mme_lattice_summation.F @@ -16,7 +16,8 @@ MODULE eri_mme_lattice_summation USE ao_util, ONLY: exp_radius USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite,& create_hermite_to_cartesian, & - exp_radius_fast + exp_radius_fast, eri_mme_coulomb, eri_mme_yukawa,& + eri_mme_longrange USE kinds, ONLY: dp,& int_8 USE mathconstants, ONLY: gaussi,& @@ -415,7 +416,7 @@ CONTAINS !> \param coulomb ... !> \note MMME Method is not very efficient for non-orthorhombic cells ! ************************************************************************************************** - PURE SUBROUTINE pgf_sum_2c_gspace_3d(S_G, l_max, R, alpha, h_inv, G_c, G_rad, vol, coulomb) + PURE SUBROUTINE pgf_sum_2c_gspace_3d(S_G, l_max, R, alpha, h_inv, G_c, G_rad, vol, potential, pot_par) REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: S_G INTEGER, INTENT(IN) :: l_max REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R @@ -423,23 +424,18 @@ CONTAINS REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_c REAL(KIND=dp), INTENT(IN) :: G_rad, vol - LOGICAL, INTENT(IN), OPTIONAL :: coulomb + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par COMPLEX(KIND=dp) :: exp_tot INTEGER, DIMENSION(3) :: l_xyz INTEGER :: gx, gy, gz, k, l, lco, lx, ly, lz COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: Ig - LOGICAL :: use_coulomb REAL(KIND=dp) :: G_rads_sq, G_sq REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l REAL(KIND=dp), DIMENSION(3, 3) :: ht - IF (PRESENT(coulomb)) THEN - use_coulomb = coulomb - ELSE - use_coulomb = .FALSE. - ENDIF ht = twopi*TRANSPOSE(h_inv) Ig(:) = 0.0_dp @@ -457,13 +453,21 @@ CONTAINS G_sq = G(1)**2+G(2)**2+G(3)**2 IF (G_sq .GT. G_rads_sq) CYCLE - IF (use_coulomb) THEN + IF (PRESENT(potential)) THEN IF (gx .EQ. 0 .AND. gy .EQ. 0 .AND. gz .EQ. 0) CYCLE ENDIF exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R)) ! cost: 2 exp_w flops - IF (use_coulomb) THEN - exp_tot = exp_tot/G_sq + IF (PRESENT(potential)) THEN + SELECT CASE (potential) + CASE (eri_mme_coulomb) + exp_tot = exp_tot/G_sq + CASE (eri_mme_yukawa) + exp_tot = exp_tot/(G_sq + pot_par**2) + !exp_tot = exp_tot/G_sq + CASE (eri_mme_longrange) + exp_tot = exp_tot/G_sq*EXP(-G_sq/pot_par**2) + END SELECT ENDIF DO k = 1, 3 G_pow_l(k, 0) = 1.0_dp diff --git a/src/eri_mme/eri_mme_test.F b/src/eri_mme/eri_mme_test.F index 8b2740251c..05b3880fb6 100644 --- a/src/eri_mme/eri_mme_test.F +++ b/src/eri_mme/eri_mme_test.F @@ -16,7 +16,11 @@ MODULE eri_mme_test create_hermite_to_cartesian USE eri_mme_integrate, ONLY: eri_mme_2c_integrate,& eri_mme_3c_integrate - USE eri_mme_types, ONLY: eri_mme_param + USE eri_mme_types, ONLY: eri_mme_coulomb,& + eri_mme_longrange,& + eri_mme_param,& + eri_mme_set_potential,& + eri_mme_yukawa USE kinds, ONLY: dp USE mathconstants, ONLY: twopi USE message_passing, ONLY: mp_sum @@ -45,19 +49,24 @@ CONTAINS !> \param test_accuracy ... !> \param para_env ... !> \param iw ... +!> \param potential ... +!> \param pot_par ... !> \param G_count ... !> \param R_count ... ! ************************************************************************************************** - SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, G_count, R_count) - TYPE(eri_mme_param), INTENT(IN) :: param + SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, & + para_env, iw, potential, pot_par, G_count, R_count) + TYPE(eri_mme_param), INTENT(INOUT) :: param INTEGER, INTENT(IN) :: l_max REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & INTENT(IN) :: zet REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rabc INTEGER, INTENT(IN) :: nrep - LOGICAL, INTENT(IN) :: test_accuracy + LOGICAL, INTENT(INOUT) :: test_accuracy TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env INTEGER, INTENT(IN) :: iw + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par INTEGER, INTENT(OUT), OPTIONAL :: G_count, R_count CHARACTER(len=*), PARAMETER :: routineN = 'eri_mme_2c_perf_acc_test', & @@ -76,6 +85,10 @@ CONTAINS nzet = SIZE(zet) nR = SIZE(rabc, 2) + IF (PRESENT(potential)) THEN + CALL eri_mme_set_potential(param, potential, pot_par) + ENDIF + ! Calculate reference values (Exact expression in G space converged to high precision) IF (test_accuracy) THEN ht = twopi*TRANSPOSE(param%h_inv) @@ -87,7 +100,7 @@ CONTAINS DO iab = 1, nR CALL eri_mme_2c_integrate(param, 0, l_max, 0, l_max, zet(izet), zet(izet), rabc(:, iab), & I_ref(:, :, iab, izet), 0, 0, & - normalize=.TRUE., exact_method=.TRUE.) + normalize=.TRUE., potential=potential, pot_par=pot_par) ENDDO ENDDO @@ -121,7 +134,7 @@ CONTAINS ENDIF IF (iw > 0) THEN - WRITE (iw, '(/T2, A)') "ERI_MME| Test results for 2c cpu time" + WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time" WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy" DO l = 0, l_max @@ -165,12 +178,15 @@ CONTAINS !> \param nsample ... !> \param para_env ... !> \param iw ... +!> \param potential ... +!> \param pot_par ... !> \param GG_count ... !> \param GR_count ... !> \param RR_count ... ! ************************************************************************************************** - SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, para_env, iw, GG_count, GR_count, RR_count) - TYPE(eri_mme_param), INTENT(IN) :: param + SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, & + para_env, iw, potential, pot_par, GG_count, GR_count, RR_count) + TYPE(eri_mme_param), INTENT(INOUT) :: param INTEGER, INTENT(IN) :: l_max REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & INTENT(IN) :: zet @@ -180,6 +196,8 @@ CONTAINS INTEGER, INTENT(IN), OPTIONAL :: nsample TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env INTEGER, INTENT(IN) :: iw + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par INTEGER, INTENT(OUT), OPTIONAL :: GG_count, GR_count, RR_count CHARACTER(len=*), PARAMETER :: routineN = 'eri_mme_3c_perf_acc_test', & @@ -200,6 +218,10 @@ CONTAINS nzet = SIZE(zet) nR = SIZE(rabc, 2) + IF (PRESENT(potential)) THEN + CALL eri_mme_set_potential(param, potential, pot_par) + ENDIF + IF (param%debug) THEN DO izeta = 1, nzet DO izetb = 1, nzet @@ -216,7 +238,19 @@ CONTAINS ENDIF IF (iw > 0) THEN - WRITE (iw, '(/T2, A)') "ERI_MME| Test results for 3c cpu time" + IF (PRESENT(potential)) THEN + SELECT CASE (potential) + CASE (eri_mme_coulomb) + WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb" + CASE (eri_mme_yukawa) + WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par + CASE (eri_mme_longrange) + WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par + END SELECT + ELSE + WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb" + ENDIF + WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time" WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time" ENDIF diff --git a/src/eri_mme/eri_mme_types.F b/src/eri_mme/eri_mme_types.F index 27ca477e92..bb39885dd0 100644 --- a/src/eri_mme/eri_mme_types.F +++ b/src/eri_mme/eri_mme_types.F @@ -15,8 +15,12 @@ MODULE eri_mme_types USE cp_para_types, ONLY: cp_para_env_type USE eri_mme_error_control, ONLY: calibrate_cutoff,& - cutoff_minimax_error - USE eri_mme_gaussian, ONLY: get_minimax_coeff_v_gspace + cutoff_minimax_error,& + minimax_error + USE eri_mme_gaussian, ONLY: eri_mme_coulomb,& + eri_mme_longrange,& + eri_mme_yukawa,& + get_minimax_coeff_v_gspace USE eri_mme_util, ONLY: G_abs_min,& R_abs_min USE kinds, ONLY: dp @@ -40,7 +44,11 @@ MODULE eri_mme_types eri_mme_release, & eri_mme_set_params, & eri_mme_print_grid_info, & - get_minimax_from_cutoff + get_minimax_from_cutoff, & + eri_mme_coulomb, & + eri_mme_yukawa, & + eri_mme_longrange, & + eri_mme_set_potential TYPE minimax_grid REAL(KIND=dp) :: cutoff @@ -62,6 +70,7 @@ MODULE eri_mme_types REAL(KIND=dp) :: cutoff_min, cutoff_max, cutoff_delta, & cutoff_eps REAL(KIND=dp) :: err_mm, err_c + REAL(KIND=dp) :: mm_delta REAL(KIND=dp) :: G_min, R_min LOGICAL :: is_valid LOGICAL :: debug @@ -73,6 +82,10 @@ MODULE eri_mme_types INTEGER :: n_grids TYPE(minimax_grid), DIMENSION(:), & ALLOCATABLE :: minimax_grid + REAL(KIND=dp) :: zet_max, zet_min + INTEGER :: l_max, l_mm, l_max_zet + INTEGER :: potential + REAL(KIND=dp) :: pot_par END TYPE eri_mme_param CONTAINS @@ -130,6 +143,8 @@ CONTAINS param%debug_nsum = debug_nsum param%print_calib = print_calib param%unit_nr = unit_nr + param%err_mm = -1.0_dp + param%err_c = -1.0_dp param%is_valid = .FALSE. ALLOCATE (param%minimax_grid(param%n_grids)) @@ -141,32 +156,37 @@ CONTAINS !> \param param ... !> \param hmat ... !> \param is_ortho ... -!> \param zet_err_minimax Exponent used to estimate error of minimax approximation. -!> \param zet_err_cutoff Exponents used to estimate error of finite cutoff. -!> \param l_err_cutoff Total ang. mom. quantum numbers l to be combined with exponents in -!> zet_err_cutoff. +!> \param zet_min Exponent used to estimate error of minimax approximation. +!> \param zet_max Exponent used to estimate error of finite cutoff. +!> \param l_max_zet Total ang. mom. quantum numbers l to be combined with exponents in +!> zet_max. !> \param l_max Maximum total angular momentum quantum number !> \param para_env ... +!> \param potential potential to use. Accepts the following values: +!> 1: coulomb potential V(r)=1/r +!> 2: yukawa potential V(r)=e(-a*r)/r +!> 3: long-range coulomb erf(a*r)/r +!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r ! ************************************************************************************************** - SUBROUTINE eri_mme_set_params(param, hmat, is_ortho, zet_err_minimax, zet_err_cutoff, l_err_cutoff, l_max, para_env) + SUBROUTINE eri_mme_set_params(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, & + potential, pot_par) TYPE(eri_mme_param), INTENT(INOUT) :: param REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat LOGICAL, INTENT(IN) :: is_ortho - REAL(KIND=dp), INTENT(IN) :: zet_err_minimax, zet_err_cutoff - INTEGER, INTENT(IN) :: l_err_cutoff, l_max - TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env + REAL(KIND=dp), INTENT(IN) :: zet_min, zet_max + INTEGER, INTENT(IN) :: l_max_zet, l_max + TYPE(cp_para_env_type), INTENT(IN), OPTIONAL, & + POINTER :: para_env + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_set_params', & routineP = moduleN//':'//routineN - INTEGER :: handle, i_grid, l_c, l_mm, n_grids, & - n_minimax, n_mm + INTEGER :: handle, l_mm, n_grids LOGICAL :: s_only - REAL(KIND=dp) :: cutoff, cutoff_delta, cutoff_max, cutoff_min, cutoff_rel, E_mm, err_c, & - err_c_prev, err_mm, err_mm_prev, err_target, exp_max, exp_min - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw, minimax_aw_prev - - !REAL(KIND=dp) :: zet_c + REAL(KIND=dp) :: cutoff + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw CALL timeset(routineN, handle) @@ -175,9 +195,8 @@ CONTAINS CALL init_orbital_pointers(3*l_max) ! allow for orbital pointers of combined index - ! l values for minimax error estimate (l_mm) and for cutoff error estimate (l_c) + ! l values for minimax error estimate (l_mm) and for cutoff error estimate (l_max_zet) l_mm = MERGE(0, 1, s_only) - l_c = l_err_cutoff ! cell info ! Note: we recompute basic quantities from hmat to avoid dependency on cp2k cell type @@ -190,6 +209,12 @@ CONTAINS param%G_min = G_abs_min(param%h_inv) param%R_min = R_abs_min(param%hmat) + ! Minimum and maximum exponents + param%zet_max = zet_max + param%zet_min = zet_min + param%l_max_zet = l_max_zet + param%l_mm = l_mm + ! cutoff calibration not yet implemented for general cell IF (.NOT. param%is_ortho) THEN param%do_calib_cutoff = .FALSE. @@ -199,9 +224,11 @@ CONTAINS n_grids = param%n_grids ! Cutoff calibration and error estimate for orthorhombic cell + ! Here we assume Coulomb potential which should give an upper bound error also for the other + ! potentials IF (param%do_calib_cutoff) THEN CALL calibrate_cutoff(param%hmat, param%h_inv, param%G_min, param%vol, & - zet_err_minimax, l_mm, zet_err_cutoff, l_c, param%n_minimax, & + zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, & param%cutoff_min, param%cutoff_max, param%cutoff_eps, & param%cutoff_delta, cutoff, param%err_mm, param%err_c, & param%C_mm, para_env, param%print_calib, param%unit_nr) @@ -210,44 +237,126 @@ CONTAINS ELSE IF (param%do_error_est) THEN ALLOCATE (minimax_aw(2*param%n_minimax)) CALL cutoff_minimax_error(param%cutoff, param%hmat, param%h_inv, param%vol, param%G_min, & - zet_err_minimax, l_mm, zet_err_cutoff, l_c, param%n_minimax, & + zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, & minimax_aw, param%err_mm, param%err_c, param%C_mm, para_env) DEALLOCATE (minimax_aw) + ENDIF + + param%is_valid = .TRUE. + + CALL eri_mme_set_potential(param, potential=potential, pot_par=pot_par) + + CALL timestop(handle) + END SUBROUTINE eri_mme_set_params + +! ************************************************************************************************** +!> \brief ... +!> \param param ... +!> \param potential potential to use. Accepts the following values: +!> 1: coulomb potential V(r)=1/r +!> 2: yukawa potential V(r)=e(-a*r)/r +!> 3: long-range coulomb erf(a*r)/r +!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r +! ************************************************************************************************** + SUBROUTINE eri_mme_set_potential(param, potential, pot_par) + TYPE(eri_mme_param), INTENT(INOUT) :: param + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par + + REAL(KIND=dp) :: cutoff_max, cutoff_min, cutoff_rel + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw + + CPASSERT(param%is_valid) + + IF (PRESENT(potential)) THEN + param%potential = potential ELSE - param%err_mm = -1.0_dp - param%err_c = -1.0_dp + param%potential = eri_mme_coulomb + ENDIF + + IF (PRESENT(pot_par)) THEN + param%pot_par = pot_par + ELSE + param%pot_par = 0.0_dp ENDIF ALLOCATE (minimax_aw(2*param%n_minimax)) - CALL get_minimax_coeff_v_gspace(param%n_minimax, param%cutoff, param%G_min, minimax_aw, E_mm) + + CALL minimax_error(param%cutoff, param%hmat, param%vol, param%G_min, param%zet_min, param%l_mm, & + param%n_minimax, minimax_aw, param%err_mm, param%mm_delta, potential=potential, pot_par=pot_par) DEALLOCATE (minimax_aw) - err_target = E_mm - - exp_max = zet_err_cutoff - exp_min = zet_err_minimax - - CPASSERT(exp_max+1.0E-12 .GT. exp_min) - CPASSERT(n_grids .GE. 1) + CPASSERT(param%zet_max+1.0E-12 .GT. param%zet_min) + CPASSERT(param%n_grids .GE. 1) cutoff_max = param%cutoff - cutoff_rel = param%cutoff/exp_max - cutoff_min = exp_min*cutoff_rel - cutoff_delta = (cutoff_max/cutoff_min)**(1.0_dp/(n_grids)) + cutoff_rel = param%cutoff/param%zet_max + cutoff_min = param%zet_min*cutoff_rel - n_minimax = param%n_minimax + CALL eri_mme_destroy_minimax_grids(param%minimax_grid) + ALLOCATE (param%minimax_grid(param%n_grids)) + + CALL eri_mme_create_minimax_grids(param%n_grids, param%minimax_grid, param%n_minimax, & + cutoff_max, cutoff_min, param%G_min, & + param%mm_delta, potential=potential, pot_par=pot_par) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param n_grids ... +!> \param minimax_grids ... +!> \param n_minimax ... +!> \param cutoff_max ... +!> \param cutoff_min ... +!> \param G_min ... +!> \param target_error ... +!> \param potential ... +!> \param pot_par ... +! ************************************************************************************************** + SUBROUTINE eri_mme_create_minimax_grids(n_grids, minimax_grids, n_minimax, & + cutoff_max, cutoff_min, G_min, & + target_error, potential, pot_par) + INTEGER, INTENT(IN) :: n_grids + TYPE(minimax_grid), DIMENSION(n_grids), & + INTENT(OUT) :: minimax_grids + INTEGER, INTENT(IN) :: n_minimax + REAL(KIND=dp), INTENT(IN) :: cutoff_max, cutoff_min, G_min, & + target_error + INTEGER, INTENT(IN), OPTIONAL :: potential + REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par + + INTEGER :: i_grid, n_mm + REAL(KIND=dp) :: cutoff, cutoff_delta, err_c, err_c_prev, & + err_mm, err_mm_prev + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw, minimax_aw_prev + + cutoff_delta = (cutoff_max/cutoff_min)**(1.0_dp/(n_grids)) cutoff = cutoff_max - DO i_grid = n_grids, 1, -1 + + ALLOCATE (minimax_aw(2*n_minimax)) + ! for first grid (for max. cutoff) always use default n_minimax + CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw, err_minimax=err_mm, & + potential=potential, pot_par=pot_par) + CPASSERT(err_mm .LT. 1.1_dp*target_error+1.0E-12) + CALL create_minimax_grid(minimax_grids(n_grids), cutoff, n_minimax, minimax_aw, err_mm) + DEALLOCATE (minimax_aw) + + DO i_grid = n_grids-1, 1, -1 DO n_mm = n_minimax, 1, -1 ALLOCATE (minimax_aw(2*n_mm)) - CALL get_minimax_coeff_v_gspace(n_mm, cutoff, param%G_min, minimax_aw, err_mm) - IF (err_mm .GT. 1.1_dp*err_target) THEN + CALL get_minimax_coeff_v_gspace(n_mm, cutoff, G_min, minimax_aw, err_minimax=err_mm, & + potential=potential, pot_par=pot_par) + + IF (err_mm .GT. 1.1_dp*target_error) THEN CPASSERT(n_mm .NE. n_minimax) - CALL create_minimax_grid(param%minimax_grid(i_grid), cutoff, n_mm+1, minimax_aw_prev, err_mm_prev) + CALL create_minimax_grid(minimax_grids(i_grid), cutoff, n_mm+1, minimax_aw_prev, err_mm_prev) + DEALLOCATE (minimax_aw) EXIT ENDIF + IF (ALLOCATED(minimax_aw_prev)) DEALLOCATE (minimax_aw_prev) ALLOCATE (minimax_aw_prev(2*n_mm)) minimax_aw_prev(:) = minimax_aw(:) @@ -257,11 +366,27 @@ CONTAINS ENDDO cutoff = cutoff/cutoff_delta ENDDO + END SUBROUTINE - param%is_valid = .TRUE. +! ************************************************************************************************** +!> \brief ... +!> \param minimax_grids ... +! ************************************************************************************************** + SUBROUTINE eri_mme_destroy_minimax_grids(minimax_grids) + TYPE(minimax_grid), ALLOCATABLE, DIMENSION(:), & + INTENT(INOUT) :: minimax_grids - CALL timestop(handle) - END SUBROUTINE eri_mme_set_params + INTEGER :: igrid + + IF (ALLOCATED(minimax_grids)) THEN + DO igrid = 1, SIZE(minimax_grids) + IF (ASSOCIATED(minimax_grids(igrid)%minimax_aw)) THEN + DEALLOCATE (minimax_grids(igrid)%minimax_aw) + ENDIF + ENDDO + DEALLOCATE (minimax_grids) + ENDIF + END SUBROUTINE ! ************************************************************************************************** !> \brief ... @@ -320,16 +445,6 @@ CONTAINS ENDIF END SUBROUTINE -! ************************************************************************************************** -!> \brief ... -!> \param grid ... -! ************************************************************************************************** - SUBROUTINE destroy_minimax_grid(grid) - TYPE(minimax_grid), INTENT(INOUT) :: grid - - DEALLOCATE (grid%minimax_aw) - END SUBROUTINE - ! ************************************************************************************************** !> \brief ... !> \param grid ... @@ -357,13 +472,8 @@ CONTAINS SUBROUTINE eri_mme_release(param) TYPE(eri_mme_param), INTENT(INOUT) :: param - INTEGER :: igrid - IF (ALLOCATED(param%minimax_grid)) THEN - DO igrid = 1, param%n_grids - CALL destroy_minimax_grid(param%minimax_grid(igrid)) - ENDDO - DEALLOCATE (param%minimax_grid) + CALL eri_mme_destroy_minimax_grids(param%minimax_grid) ENDIF END SUBROUTINE eri_mme_release diff --git a/src/mp2_eri.F b/src/mp2_eri.F index ffc00189aa..be0381f67e 100644 --- a/src/mp2_eri.F +++ b/src/mp2_eri.F @@ -1051,23 +1051,23 @@ $: 'habdc=habdc, &'*abdc !> \note ! ************************************************************************************************** SUBROUTINE integrate_set_3c(param, & - la_min, la_max, lb_min, lb_max, lc_min, lc_max, & - npgfa, npgfb, npgfc, & - zeta, zetb, zetc, & - ra, rb, rc, & - habc, & - n_habc_a, n_habc_b, n_habc_c, & - offset_habc_a, offset_habc_b, offset_habc_c, & - offset_set_a, offset_set_b, offset_set_c, & - sphi_a, sphi_b, sphi_c, & - sgfa, sgfb, sgfc, & - nsgfa, nsgfb, nsgfc, & - eri_method, & - pabc, & - force_a, force_b, force_c, & - do_symmetric, on_diagonal, & - hdabc, hadbc, habdc, & - GG_count, GR_count, RR_count) + la_min, la_max, lb_min, lb_max, lc_min, lc_max, & + npgfa, npgfb, npgfc, & + zeta, zetb, zetc, & + ra, rb, rc, & + habc, & + n_habc_a, n_habc_b, n_habc_c, & + offset_habc_a, offset_habc_b, offset_habc_c, & + offset_set_a, offset_set_b, offset_set_c, & + sphi_a, sphi_b, sphi_c, & + sgfa, sgfb, sgfc, & + nsgfa, nsgfb, nsgfc, & + eri_method, & + pabc, & + force_a, force_b, force_c, & + do_symmetric, on_diagonal, & + hdabc, hadbc, habdc, & + GG_count, GR_count, RR_count) TYPE(eri_mme_param), INTENT(INOUT) :: param INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, & diff --git a/src/qmmm_init.F b/src/qmmm_init.F index 5c72266a77..7a28fbaf1e 100644 --- a/src/qmmm_init.F +++ b/src/qmmm_init.F @@ -843,9 +843,9 @@ CONTAINS CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param) CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, & hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, & - zet_err_minimax=qmmm_env%image_charge_pot%eta, & - zet_err_cutoff=qmmm_env%image_charge_pot%eta, & - l_err_cutoff=0, & + zet_min=qmmm_env%image_charge_pot%eta, & + zet_max=qmmm_env%image_charge_pot%eta, & + l_max_zet=0, & l_max=0, & para_env=para_env) diff --git a/tests/LIBTEST/test_eri_mme_accuracy_longrange.inp b/tests/LIBTEST/test_eri_mme_accuracy_longrange.inp new file mode 100644 index 0000000000..4d25c40169 --- /dev/null +++ b/tests/LIBTEST/test_eri_mme_accuracy_longrange.inp @@ -0,0 +1,27 @@ +&GLOBAL + PROJECT test_eri_mme + PRINT_LEVEL MEDIUM + PROGRAM_NAME TEST + RUN_TYPE NONE +&END GLOBAL +&TEST + &ERI_MME_TEST + ABC 5.0 4.0 6.0 + MIN_NPOS 8 + NREP 1 + CHECK_2C_ACCURACY .TRUE. + LMAX 5 + ZET_MIN 0.001 + ZET_MAX 1.0 + NZET 4 + TEST_3C .FALSE. + POTENTIAL LONGRANGE + POTENTIAL_PARAM 0.25 + &ERI_MME + N_MINIMAX 15 + &CUTOFF_CALIB + EPS 0.001 + &END CUTOFF_CALIB + &END ERI_MME + &END ERI_MME_TEST +&END TEST diff --git a/tests/LIBTEST/test_eri_mme_accuracy_yukawa.inp b/tests/LIBTEST/test_eri_mme_accuracy_yukawa.inp new file mode 100644 index 0000000000..818eb75248 --- /dev/null +++ b/tests/LIBTEST/test_eri_mme_accuracy_yukawa.inp @@ -0,0 +1,27 @@ +&GLOBAL + PROJECT test_eri_mme + PRINT_LEVEL MEDIUM + PROGRAM_NAME TEST + RUN_TYPE NONE +&END GLOBAL +&TEST + &ERI_MME_TEST + ABC 5.0 4.0 6.0 + MIN_NPOS 8 + NREP 1 + CHECK_2C_ACCURACY .TRUE. + LMAX 5 + ZET_MIN 0.001 + ZET_MAX 1.0 + NZET 4 + TEST_3C .FALSE. + POTENTIAL YUKAWA + POTENTIAL_PARAM 0.5 + &ERI_MME + N_MINIMAX 15 + &CUTOFF_CALIB + EPS 0.001 + &END CUTOFF_CALIB + &END ERI_MME + &END ERI_MME_TEST +&END TEST