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 : * transcendental equation for finding maximum of the type 1 integrand
143 : */
144 0 : static double screening_type1_equation_for_maximum(double r, int n, int lambda,
145 : double p, double k) {
146 0 : double K_ratio = 0.0;
147 0 : if (lambda == 0) {
148 0 : K_ratio = libgrpp_modified_bessel_scaled(1, k * r) /
149 0 : libgrpp_modified_bessel_scaled(0, k * r);
150 : } else {
151 0 : K_ratio = libgrpp_modified_bessel_scaled(lambda - 1, k * r) /
152 0 : libgrpp_modified_bessel_scaled(lambda, k * r);
153 : }
154 :
155 0 : double a = n + K_ratio * k * r;
156 0 : if (lambda > 0) {
157 0 : a = a - lambda - 1;
158 : }
159 :
160 0 : return sqrt(a / (2.0 * p));
161 : }
162 :
163 : /**
164 : * screening for the type 2 radial integrals
165 : * for the pair of contracted gaussian functions.
166 : */
167 7598315 : int libgrpp_screening_radial_type2(int lambda1, int lambda2, int n, double CA_2,
168 : double CB_2, libgrpp_shell_t *bra,
169 : libgrpp_shell_t *ket,
170 : libgrpp_potential_t *potential,
171 : double *screened_value) {
172 7598315 : *screened_value = 0.0;
173 :
174 7598315 : double CA = sqrt(CA_2);
175 7598315 : double CB = sqrt(CB_2);
176 :
177 : /*
178 : * loop over 'bra' contracted function
179 : */
180 12221890 : for (int i = 0; i < bra->num_primitives; i++) {
181 7598315 : double alpha_A = bra->alpha[i];
182 7598315 : double coef_i = bra->coeffs[i];
183 7598315 : double k1 = 2 * alpha_A * CA;
184 :
185 : /*
186 : * loop over 'ket' contracted function
187 : */
188 12221890 : for (int j = 0; j < ket->num_primitives; j++) {
189 7598315 : double alpha_B = ket->alpha[j];
190 7598315 : double coef_j = ket->coeffs[j];
191 7598315 : double k2 = 2 * alpha_B * CB;
192 :
193 : /*
194 : * loop over RPP primitives
195 : */
196 14195744 : for (int k = 0; k < potential->num_primitives; k++) {
197 9572169 : double eta = potential->alpha[k];
198 9572169 : int ni = n + potential->powers[k];
199 9572169 : double coef_k = potential->coeffs[k];
200 :
201 9572169 : double val_ijk = 0.0;
202 9572169 : int err_code = screening_radial_type2_integral_primitive(
203 : lambda1, lambda2, ni, CA_2, CB_2, alpha_A, alpha_B, k1, k2, eta,
204 : &val_ijk);
205 9572169 : if (err_code == EXIT_FAILURE) {
206 2974740 : return EXIT_FAILURE;
207 : }
208 :
209 6597429 : *screened_value += coef_i * coef_j * coef_k * val_ijk;
210 : }
211 : }
212 : }
213 :
214 : return EXIT_SUCCESS;
215 : }
216 :
217 : /**
218 : * Analytically evaluates Gaussian integral:
219 : * \int_0^\infty r^n e^(-a r^2) dr
220 : */
221 72324 : double libgrpp_gaussian_integral(int n, double a) {
222 72324 : if (n % 2 == 0) {
223 42215 : int k = n / 2;
224 42215 : return libgrpp_double_factorial(2 * k - 1) / (pow(2.0, k + 1) * pow(a, k)) *
225 42215 : sqrt(M_PI / a);
226 : } else {
227 30109 : int k = (n - 1) / 2;
228 30109 : return libgrpp_factorial(k) / (2.0 * pow(a, k + 1));
229 : }
230 : }
231 :
232 : /**
233 : * screening for the type 2 radial integrals
234 : * for the pair of primitive gaussian functions.
235 : */
236 9572169 : static int screening_radial_type2_integral_primitive(
237 : int lambda1, int lambda2, int n, double CA_2, double CB_2, double alpha_A,
238 : double alpha_B, double k1, double k2, double eta, double *screened_value) {
239 9572169 : *screened_value = 0.0;
240 :
241 9572169 : if (lambda1 >= 1 && fabs(k1) <= LIBGRPP_ZERO_THRESH) {
242 : return EXIT_SUCCESS;
243 : }
244 8266401 : if (lambda2 >= 1 && fabs(k2) <= LIBGRPP_ZERO_THRESH) {
245 : return EXIT_SUCCESS;
246 : }
247 :
248 6996809 : double p = alpha_A + alpha_B + eta;
249 6996809 : double CA = sqrt(CA_2);
250 6996809 : double CB = sqrt(CB_2);
251 :
252 : /*
253 : * special case:
254 : * lambda1 = lambda2 = 0,
255 : * k1 = k2 = 0.
256 : * => M_0(0) = 1
257 : * => we have one-center integral which can be evaluated analytically
258 : */
259 6996809 : if (lambda1 == 0 && lambda2 == 0) {
260 598905 : if (fabs(k1) <= LIBGRPP_ZERO_THRESH && fabs(k2) <= LIBGRPP_ZERO_THRESH) {
261 0 : *screened_value = exp(-alpha_A * CA * CA - alpha_B * CB * CB) *
262 0 : libgrpp_gaussian_integral(n, p);
263 0 : return EXIT_SUCCESS;
264 : }
265 : }
266 :
267 : /*
268 : * find position of the maximum of the integrand
269 : */
270 6996809 : const double tol = 1e-2;
271 6996809 : double r0 = (alpha_A * CA + alpha_B * CB) / p;
272 6996809 : double r0_prev = 0.0;
273 6996809 : int nsteps = 0;
274 :
275 27987208 : do {
276 27987208 : nsteps++;
277 27987208 : if (nsteps == 5) {
278 2974740 : *screened_value = 0.0;
279 2974740 : return EXIT_FAILURE;
280 : }
281 25012468 : r0_prev = r0;
282 25012468 : r0 = screening_type2_equation_for_maximum(r0, n, lambda1, lambda2, p, k1,
283 : k2);
284 25012468 : } while (fabs(r0 - r0_prev) > tol);
285 :
286 : /*
287 : * envelope function for the integrand
288 : */
289 12066207 : *screened_value = sqrt(M_PI / p) * pow(r0, n) *
290 4022069 : libgrpp_modified_bessel_scaled(lambda1, k1 * r0) *
291 4022069 : libgrpp_modified_bessel_scaled(lambda2, k2 * r0) *
292 4022069 : exp(-eta * r0 * r0 - alpha_A * (r0 - CA) * (r0 - CA) -
293 4022069 : alpha_B * (r0 - CB) * (r0 - CB)) *
294 4022069 : 0.5 * (1 + erf(sqrt(p) * r0));
295 :
296 4022069 : return EXIT_SUCCESS;
297 : }
298 :
299 : /**
300 : * transcendental equation for finding maximum of the type 2 integrand
301 : */
302 25012468 : static double screening_type2_equation_for_maximum(double r, int n, int lambda1,
303 : int lambda2, double p,
304 : double k1, double k2) {
305 25012468 : double K1_ratio = 0.0;
306 25012468 : double k1_r = k1 * r;
307 25012468 : if (lambda1 == 0) {
308 12929966 : K1_ratio = libgrpp_modified_bessel_scaled(1, k1_r) /
309 6464983 : libgrpp_modified_bessel_scaled(0, k1_r);
310 : } else {
311 37094970 : K1_ratio = libgrpp_modified_bessel_scaled(lambda1 - 1, k1_r) /
312 18547485 : libgrpp_modified_bessel_scaled(lambda1, k1_r);
313 : }
314 :
315 25012468 : double K2_ratio = 0.0;
316 25012468 : double k2_r = k2 * r;
317 25012468 : if (lambda2 == 0) {
318 12885116 : K2_ratio = libgrpp_modified_bessel_scaled(1, k2_r) /
319 6442558 : libgrpp_modified_bessel_scaled(0, k2_r);
320 : } else {
321 37139820 : K2_ratio = libgrpp_modified_bessel_scaled(lambda2 - 1, k2_r) /
322 18569910 : libgrpp_modified_bessel_scaled(lambda2, k2_r);
323 : }
324 :
325 25012468 : double a = K1_ratio * k1_r + K2_ratio * k2_r + n;
326 :
327 25012468 : if (lambda1 > 0) {
328 18547485 : a = a - lambda1 - 1;
329 : }
330 25012468 : if (lambda2 > 0) {
331 18569910 : a = a - lambda2 - 1;
332 : }
333 :
334 25012468 : return sqrt(a / (2.0 * p));
335 : }
|
