Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: MIT */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : /*
9 : * libgrpp - a library for the evaluation of integrals over
10 : * generalized relativistic pseudopotentials.
11 : *
12 : * Copyright (C) 2021-2023 Alexander Oleynichenko
13 : */
14 :
15 : /*
16 : * Screening of radial integrals.
17 : *
18 : * The technique of screening is adopted from:
19 : * R. A. Shaw, J. G. Hill. Prescreening and efficiency in the evaluation
20 : * of integrals over ab initio effective core potentials.
21 : * J. Chem. Phys. 147, 074108 (2017). doi: 10.1063/1.4986887
22 : * (see also Supplementary Material for this article).
23 : *
24 : * Note that in this publication the transcendental equation (2) for
25 : * type 2 integrals is not correct.
26 : */
27 :
28 : #include <math.h>
29 : #include <stdlib.h>
30 : #ifndef M_PI
31 : #define M_PI 3.14159265358979323846
32 : #endif
33 :
34 : #include "grpp_factorial.h"
35 : #include "grpp_screening.h"
36 : #include "grpp_specfunc.h"
37 : #include "libgrpp.h"
38 :
39 : /*
40 : * functions defined below in the file
41 : */
42 :
43 : static int screening_radial_type1_integral_primitive(
44 : int lambda, int n, double CA_2, double CB_2, double alpha_A, double alpha_B,
45 : double k, double eta, double *screened_value);
46 :
47 : static double screening_type1_equation_for_maximum(double r, int n, int lambda,
48 : double p, double k);
49 :
50 : static int screening_radial_type2_integral_primitive(
51 : int lambda1, int lambda2, int n, double CA_2, double CB_2, double alpha_A,
52 : double alpha_B, double k1, double k2, double eta, double *screened_value);
53 :
54 : static double screening_type2_equation_for_maximum(double r, int n, int lambda1,
55 : int lambda2, double p,
56 : double k1, double k2);
57 :
58 : // static double analytic_one_center_rpp_integral_primitive(int L, double
59 : // alpha1,
60 : // double alpha2, int
61 : // n, double zeta);
62 :
63 : /**
64 : * screening for the type 1 radial integrals
65 : * for the pair of contracted gaussian functions.
66 : */
67 0 : int libgrpp_screening_radial_type1(int lambda, int n, double CA_2, double CB_2,
68 : double alpha_A, double alpha_B, double k,
69 : double prefactor,
70 : libgrpp_potential_t *potential,
71 : double *screened_value) {
72 0 : *screened_value = 0.0;
73 :
74 0 : if (lambda >= 1 && fabs(k) <= LIBGRPP_ZERO_THRESH) {
75 : return EXIT_SUCCESS;
76 : }
77 :
78 : /*
79 : * loop over RPP primitives
80 : */
81 0 : for (int iprim = 0; iprim < potential->num_primitives; iprim++) {
82 0 : double eta = potential->alpha[iprim];
83 0 : int ni = n + potential->powers[iprim];
84 0 : double coef = potential->coeffs[iprim];
85 :
86 0 : double val_i = 0.0;
87 0 : int err_code = screening_radial_type1_integral_primitive(
88 : lambda, ni, CA_2, CB_2, alpha_A, alpha_B, k, eta, &val_i);
89 0 : if (err_code == EXIT_FAILURE) {
90 0 : return EXIT_FAILURE;
91 : }
92 :
93 0 : *screened_value += prefactor * coef * val_i;
94 : }
95 :
96 : return EXIT_SUCCESS;
97 : }
98 :
99 : /**
100 : * screening for the type 1 radial integrals
101 : * for the pair of primitive gaussian functions.
102 : */
103 0 : static int screening_radial_type1_integral_primitive(
104 : int lambda, int n, double CA_2, double CB_2, double alpha_A, double alpha_B,
105 : double k, double eta, double *screened_value) {
106 0 : double p = alpha_A + alpha_B + eta;
107 0 : double CA = sqrt(CA_2);
108 0 : double CB = sqrt(CB_2);
109 :
110 : /*
111 : * find position of the maximum of the integrand
112 : */
113 0 : const double tol = 1e-2;
114 0 : double r0 = (alpha_A * CA + alpha_B * CB) / p;
115 0 : double r0_prev = 0.0;
116 :
117 0 : int nsteps = 0;
118 0 : do {
119 0 : nsteps++;
120 0 : if (nsteps == 10) {
121 0 : *screened_value = 0.0;
122 0 : return EXIT_FAILURE;
123 : }
124 :
125 0 : r0_prev = r0;
126 0 : r0 = screening_type1_equation_for_maximum(r0, n, lambda, p, k);
127 0 : } while (fabs(r0 - r0_prev) > tol);
128 :
129 : /*
130 : * envelope function for the integrand
131 : */
132 0 : *screened_value =
133 0 : sqrt(M_PI / p) * pow(r0, n) *
134 0 : libgrpp_modified_bessel_scaled(lambda, k * r0) *
135 0 : exp(-p * r0 * r0 - alpha_A * CA_2 - alpha_B * CB_2 + k * r0) * 0.5 *
136 0 : (1 + erf(sqrt(p) * r0));
137 :
138 0 : return EXIT_SUCCESS;
139 : }
140 :
141 : /**
142 : * ratio of two modified scaled Bessel functions guarded against divide by zero
143 : */
144 55895616 : static double modified_bessel_scaled_ratio(int n, double x) {
145 55895616 : double numerator, denominator;
146 55895616 : if (n == 0) {
147 15676608 : numerator = libgrpp_modified_bessel_scaled(1, x);
148 15676608 : denominator = libgrpp_modified_bessel_scaled(0, x);
149 : } else {
150 40219008 : numerator = libgrpp_modified_bessel_scaled(n - 1, x);
151 40219008 : denominator = libgrpp_modified_bessel_scaled(n, x);
152 : }
153 55895616 : if (denominator == 0.0) {
154 57504 : return (2 * n + 1) / x; // asymptote for x->0
155 : } else {
156 55838112 : return numerator / denominator;
157 : }
158 : }
159 :
160 : /**
161 : * transcendental equation for finding maximum of the type 1 integrand
162 : */
163 0 : static double screening_type1_equation_for_maximum(double r, int n, int lambda,
164 : double p, double k) {
165 0 : double k_r = k * r;
166 0 : double K_ratio = modified_bessel_scaled_ratio(lambda, k_r);
167 0 : double a = n + K_ratio * k_r;
168 0 : if (lambda > 0) {
169 0 : a = a - lambda - 1;
170 : }
171 :
172 0 : return sqrt(a / (2.0 * p));
173 : }
174 :
175 : /**
176 : * screening for the type 2 radial integrals
177 : * for the pair of contracted gaussian functions.
178 : */
179 8041609 : int libgrpp_screening_radial_type2(int lambda1, int lambda2, int n, double CA_2,
180 : double CB_2, libgrpp_shell_t *bra,
181 : libgrpp_shell_t *ket,
182 : libgrpp_potential_t *potential,
183 : double *screened_value) {
184 8041609 : *screened_value = 0.0;
185 :
186 8041609 : double CA = sqrt(CA_2);
187 8041609 : double CB = sqrt(CB_2);
188 :
189 : /*
190 : * loop over 'bra' contracted function
191 : */
192 13108478 : for (int i = 0; i < bra->num_primitives; i++) {
193 8041609 : double alpha_A = bra->alpha[i];
194 8041609 : double coef_i = bra->coeffs[i];
195 8041609 : double k1 = 2 * alpha_A * CA;
196 :
197 : /*
198 : * loop over 'ket' contracted function
199 : */
200 13108478 : for (int j = 0; j < ket->num_primitives; j++) {
201 8041609 : double alpha_B = ket->alpha[j];
202 8041609 : double coef_j = ket->coeffs[j];
203 8041609 : double k2 = 2 * alpha_B * CB;
204 :
205 : /*
206 : * loop over RPP primitives
207 : */
208 17178776 : for (int k = 0; k < potential->num_primitives; k++) {
209 12111907 : double eta = potential->alpha[k];
210 12111907 : int ni = n + potential->powers[k];
211 12111907 : double coef_k = potential->coeffs[k];
212 :
213 12111907 : double val_ijk = 0.0;
214 12111907 : int err_code = screening_radial_type2_integral_primitive(
215 : lambda1, lambda2, ni, CA_2, CB_2, alpha_A, alpha_B, k1, k2, eta,
216 : &val_ijk);
217 12111907 : if (err_code == EXIT_FAILURE) {
218 2974740 : return EXIT_FAILURE;
219 : }
220 :
221 9137167 : *screened_value += coef_i * coef_j * coef_k * val_ijk;
222 : }
223 : }
224 : }
225 :
226 : return EXIT_SUCCESS;
227 : }
228 :
229 : /**
230 : * Analytically evaluates Gaussian integral:
231 : * \int_0^\infty r^n e^(-a r^2) dr
232 : */
233 504052 : double libgrpp_gaussian_integral(int n, double a) {
234 504052 : if (n % 2 == 0) {
235 293255 : int k = n / 2;
236 293255 : return libgrpp_double_factorial(2 * k - 1) / (pow(2.0, k + 1) * pow(a, k)) *
237 293255 : sqrt(M_PI / a);
238 : } else {
239 210797 : int k = (n - 1) / 2;
240 210797 : return libgrpp_factorial(k) / (2.0 * pow(a, k + 1));
241 : }
242 : }
243 :
244 : /**
245 : * screening for the type 2 radial integrals
246 : * for the pair of primitive gaussian functions.
247 : */
248 12111907 : static int screening_radial_type2_integral_primitive(
249 : int lambda1, int lambda2, int n, double CA_2, double CB_2, double alpha_A,
250 : double alpha_B, double k1, double k2, double eta, double *screened_value) {
251 12111907 : *screened_value = 0.0;
252 :
253 12111907 : if (lambda1 >= 1 && fabs(k1) <= LIBGRPP_ZERO_THRESH) {
254 : return EXIT_SUCCESS;
255 : }
256 10396807 : if (lambda2 >= 1 && fabs(k2) <= LIBGRPP_ZERO_THRESH) {
257 : return EXIT_SUCCESS;
258 : }
259 :
260 8145963 : double p = alpha_A + alpha_B + eta;
261 8145963 : double CA = sqrt(CA_2);
262 8145963 : double CB = sqrt(CB_2);
263 :
264 : /*
265 : * special case:
266 : * lambda1 = lambda2 = 0,
267 : * k1 = k2 = 0.
268 : * => M_0(0) = 1
269 : * => we have one-center integral which can be evaluated analytically
270 : */
271 8145963 : if (lambda1 == 0 && lambda2 == 0) {
272 847945 : if (fabs(k1) <= LIBGRPP_ZERO_THRESH && fabs(k2) <= LIBGRPP_ZERO_THRESH) {
273 0 : *screened_value = exp(-alpha_A * CA * CA - alpha_B * CB * CB) *
274 0 : libgrpp_gaussian_integral(n, p);
275 0 : return EXIT_SUCCESS;
276 : }
277 : }
278 :
279 : /*
280 : * find position of the maximum of the integrand
281 : */
282 8145963 : const double tol = 1e-2;
283 8145963 : double r0 = (alpha_A * CA + alpha_B * CB) / p;
284 8145963 : double r0_prev = 0.0;
285 8145963 : int nsteps = 0;
286 :
287 30922548 : do {
288 30922548 : nsteps++;
289 30922548 : if (nsteps == 5) {
290 2974740 : *screened_value = 0.0;
291 2974740 : return EXIT_FAILURE;
292 : }
293 27947808 : r0_prev = r0;
294 27947808 : r0 = screening_type2_equation_for_maximum(r0, n, lambda1, lambda2, p, k1,
295 : k2);
296 27947808 : } while (fabs(r0 - r0_prev) > tol);
297 :
298 : /*
299 : * envelope function for the integrand
300 : */
301 15513669 : *screened_value = sqrt(M_PI / p) * pow(r0, n) *
302 5171223 : libgrpp_modified_bessel_scaled(lambda1, k1 * r0) *
303 5171223 : libgrpp_modified_bessel_scaled(lambda2, k2 * r0) *
304 5171223 : exp(-eta * r0 * r0 - alpha_A * (r0 - CA) * (r0 - CA) -
305 5171223 : alpha_B * (r0 - CB) * (r0 - CB)) *
306 5171223 : 0.5 * (1 + erf(sqrt(p) * r0));
307 :
308 5171223 : return EXIT_SUCCESS;
309 : }
310 :
311 : /**
312 : * transcendental equation for finding maximum of the type 2 integrand
313 : */
314 27947808 : static double screening_type2_equation_for_maximum(double r, int n, int lambda1,
315 : int lambda2, double p,
316 : double k1, double k2) {
317 27947808 : double k1_r = k1 * r;
318 27947808 : double K1_ratio = modified_bessel_scaled_ratio(lambda1, k1_r);
319 :
320 27947808 : double k2_r = k2 * r;
321 27947808 : double K2_ratio = modified_bessel_scaled_ratio(lambda2, k2_r);
322 :
323 27947808 : double a = K1_ratio * k1_r + K2_ratio * k2_r + n;
324 :
325 27947808 : if (lambda1 > 0) {
326 20272638 : a = a - lambda1 - 1;
327 : }
328 27947808 : if (lambda2 > 0) {
329 19946370 : a = a - lambda2 - 1;
330 : }
331 :
332 27947808 : return sqrt(a / (2.0 * p));
333 : }
|
