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

Grpp: Fix floating point exception in grpp_screening.c

Browse source
This commit is contained in:
Ole Schütt 2025-06-30 13:54:49 +02:00 • committed by Ole Schütt
parent 12a8a74dd0
commit 082fe5dfbd
3 changed files with 26 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

50
src/grpp/grpp_screening.c
View file

@ -138,21 +138,33 @@ static int screening_radial_type1_integral_primitive(
return EXIT_SUCCESS;
}
/**
* ratio of two modified scaled Bessel functions guarded against divide by zero
*/
static double modified_bessel_scaled_ratio(int n, double x) {
double numerator, denominator;
if (n == 0) {
numerator = libgrpp_modified_bessel_scaled(1, x);
denominator = libgrpp_modified_bessel_scaled(0, x);
} else {
numerator = libgrpp_modified_bessel_scaled(n - 1, x);
denominator = libgrpp_modified_bessel_scaled(n, x);
}
if (denominator == 0.0) {
return (2 * n + 1) / x; // asymptote for x->0
} else {
return numerator / denominator;
}
}
/**
* transcendental equation for finding maximum of the type 1 integrand
*/
static double screening_type1_equation_for_maximum(double r, int n, int lambda,
double p, double k) {
double K_ratio = 0.0;
if (lambda == 0) {
K_ratio = libgrpp_modified_bessel_scaled(1, k * r) /
libgrpp_modified_bessel_scaled(0, k * r);
} else {
K_ratio = libgrpp_modified_bessel_scaled(lambda - 1, k * r) /
libgrpp_modified_bessel_scaled(lambda, k * r);
}
double a = n + K_ratio * k * r;
double k_r = k * r;
double K_ratio = modified_bessel_scaled_ratio(lambda, k_r);
double a = n + K_ratio * k_r;
if (lambda > 0) {
a = a - lambda - 1;
}
@ -302,25 +314,11 @@ static int screening_radial_type2_integral_primitive(
static double screening_type2_equation_for_maximum(double r, int n, int lambda1,
int lambda2, double p,
double k1, double k2) {
double K1_ratio = 0.0;
double k1_r = k1 * r;
if (lambda1 == 0) {
K1_ratio = libgrpp_modified_bessel_scaled(1, k1_r) /
libgrpp_modified_bessel_scaled(0, k1_r);
} else {
K1_ratio = libgrpp_modified_bessel_scaled(lambda1 - 1, k1_r) /
libgrpp_modified_bessel_scaled(lambda1, k1_r);
}
double K1_ratio = modified_bessel_scaled_ratio(lambda1, k1_r);
double K2_ratio = 0.0;
double k2_r = k2 * r;
if (lambda2 == 0) {
K2_ratio = libgrpp_modified_bessel_scaled(1, k2_r) /
libgrpp_modified_bessel_scaled(0, k2_r);
} else {
K2_ratio = libgrpp_modified_bessel_scaled(lambda2 - 1, k2_r) /
libgrpp_modified_bessel_scaled(lambda2, k2_r);
}
double K2_ratio = modified_bessel_scaled_ratio(lambda2, k2_r);
double a = K1_ratio * k1_r + K2_ratio * k2_r + n;

2
tests/QS/regtest-ecp-2/TEST_FILES.toml
View file

@ -1,3 +1,3 @@
"ICl_lanl2dz_gpw.inp" = []
"Rn_stuttgart_gapw.inp" = [{matcher="M031", tol=1.0E-08, ref=3.56470376649E+03}]
#SbH3_def2_gapw.inp 11 1.0E-11 -241.595963788403253
"SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.595963788403253}]

2
tests/QS/regtest-ecp/TEST_FILES.toml
View file

@ -1,5 +1,5 @@
"ICl_lanl2dz_gpw.inp" = [{matcher="M011", tol=1.0E-11, ref=-23.515911094537550}]
"Rn_stuttgart_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-25.632304810373267}]
#SbH3_def2_gapw.inp 11 1.0E-11 -241.414769303399083
"SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.414769303399083}]
"HCl_ccECP.inp" = [{matcher="M011", tol=1.0E-12, ref=-15.457709431781252}]
"HCl_force.inp" = []