comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

MME integrals: implement new potentials

Browse source
This commit is contained in:
Patrick Seewald 2019-05-14 12:54:09 +02:00
parent da5dcd41b5
commit 469a92f17a
11 changed files with 598 additions and 219 deletions
Download patch file Download diff file Expand all files Collapse all files

135
src/cp_eri_mme_interface.F
View file

@ -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)

191
src/eri_mme/eri_mme_error_control.F
View file

@ -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)

71
src/eri_mme/eri_mme_gaussian.F
View file

@ -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

18
src/eri_mme/eri_mme_integrate.F
View file

@ -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

28
src/eri_mme/eri_mme_lattice_summation.F
View file

@ -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

52
src/eri_mme/eri_mme_test.F
View file

@ -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

228
src/eri_mme/eri_mme_types.F
View file

@ -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

34
src/mp2_eri.F
View file

@ -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, &

6
src/qmmm_init.F
View file

@ -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)

27
tests/LIBTEST/test_eri_mme_accuracy_longrange.inp Normal file
View file

@ -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

27
tests/LIBTEST/test_eri_mme_accuracy_yukawa.inp Normal file
View file

@ -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