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 : * Evaluation of type 2 radial integrals.
17 : *
18 : * The procedure in general follows that described in:
19 : * R. Flores-Moreno et al. Half-numerical evaluation of pseudopotential
20 : * integrals. J. Comp. Chem. 27, 1009 (2006) (see formulas (12) and (13) for
21 : * radial integrals invoking contracted Gaussian functions and RPPs)
22 : *
23 : * The Log3 integration scheme used here is detailed in:
24 : * C.-K. Skylaris et al. An efficient method for calculating effective core
25 : * potential integrals which involve projection operators. Chem. Phys. Lett.
26 : * 296, 445 (1998)
27 : */
28 :
29 : #include <math.h>
30 : #include <stdlib.h>
31 : #ifndef M_PI
32 : #define M_PI 3.14159265358979323846
33 : #endif
34 :
35 : #include "grpp_radial_type2_integral.h"
36 :
37 : #include "grpp_norm_gaussian.h"
38 : #include "grpp_screening.h"
39 : #include "grpp_specfunc.h"
40 : #include "grpp_utils.h"
41 :
42 : #define MIN_GRID 31
43 : #define MAX_GRID 10000
44 :
45 : typedef struct {
46 : double CA;
47 : double CB;
48 : libgrpp_potential_t *potential;
49 : libgrpp_shell_t *bra;
50 : libgrpp_shell_t *ket;
51 : } radial_type2_params_t;
52 :
53 : /**
54 : * RPP and radial contracted Gaussians are pre-calculated on a grid,
55 : * and then combined into radial integrals
56 : */
57 : typedef struct {
58 : int nr;
59 : int n_max;
60 : int lambda1_max;
61 : int lambda2_max;
62 : double *r;
63 : double *w;
64 : double *rpp_values;
65 : double **r_N;
66 : double **F1;
67 : double **F2;
68 : radial_type2_params_t *params;
69 : } radial_type2_grid_t;
70 :
71 : /**
72 : * pre-definitions of the functions used below
73 : */
74 :
75 : static double calculate_radial_type2_integral(radial_type2_grid_t *grid, int n,
76 : int lambda1, int lambda2,
77 : double tolerance, int *converged);
78 :
79 : static radial_type2_grid_t *
80 : create_radial_type2_grid(int lambda1_max, int lambda2_max, int n_max,
81 : radial_type2_params_t *params);
82 :
83 : static void delete_radial_type2_grid(radial_type2_grid_t *grid);
84 :
85 : static double radial_type2_integrand_fun_contracted(double r, int lambda,
86 : double *k, double CA,
87 : libgrpp_shell_t *gauss_fun);
88 :
89 : static void expand_radial_type2_grid(radial_type2_grid_t *grid, int nr);
90 :
91 : static void calc_k_values(int nprim, const double *alpha, double CA, double *k);
92 :
93 : double libgrpp_gaussian_integral(int n, double a);
94 :
95 : /**
96 : * Creates table with pre-calculated radial type 2 integrals.
97 : */
98 149632 : radial_type2_table_t *libgrpp_tabulate_radial_type2_integrals(
99 : int lambda1_max, int lambda2_max, int n_max, double CA_2, double CB_2,
100 : libgrpp_potential_t *potential, libgrpp_shell_t *bra,
101 : libgrpp_shell_t *ket) {
102 : /*
103 : * create empty table containing pre-tabulated radial type 2 integrals
104 : */
105 149632 : radial_type2_table_t *table;
106 149632 : table = (radial_type2_table_t *)calloc(1, sizeof(radial_type2_table_t));
107 149632 : table->lambda1_max = lambda1_max;
108 149632 : table->lambda2_max = lambda2_max;
109 149632 : table->n_max = n_max;
110 149632 : table->radial_integrals = (double *)calloc(
111 149632 : (lambda1_max + 1) * (lambda2_max + 1) * (n_max + 1), sizeof(double));
112 :
113 : /*
114 : * the special case of one-center RPP integrals
115 : */
116 149632 : if (CA_2 < 1e-14 && CB_2 < 1e-14) {
117 :
118 23850 : for (int i = 0; i < bra->num_primitives; i++) {
119 11925 : double alpha_A = bra->alpha[i];
120 23850 : double coef_i =
121 11925 : bra->coeffs[i] * libgrpp_gaussian_norm_factor(bra->L, 0, 0, alpha_A);
122 :
123 23850 : for (int j = 0; j < ket->num_primitives; j++) {
124 11925 : double alpha_B = ket->alpha[j];
125 23850 : double coef_j = ket->coeffs[j] *
126 11925 : libgrpp_gaussian_norm_factor(ket->L, 0, 0, alpha_B);
127 :
128 39393 : for (int k = 0; k < potential->num_primitives; k++) {
129 27468 : double eta = potential->alpha[k];
130 27468 : int ni = potential->powers[k];
131 27468 : double coef_k = potential->coeffs[k];
132 :
133 27468 : double p = alpha_A + alpha_B + eta;
134 27468 : double factor = coef_i * coef_j * coef_k;
135 :
136 99792 : for (int n = 0; n <= n_max; n++) {
137 :
138 72324 : double val_ijk = libgrpp_gaussian_integral(ni + n, p);
139 72324 : table->radial_integrals[n] += factor * val_ijk;
140 72324 : ;
141 : }
142 : }
143 : }
144 : }
145 :
146 : return table;
147 : }
148 :
149 : /*
150 : * for numerical integration on the grid
151 : */
152 137707 : radial_type2_params_t params;
153 137707 : params.CA = sqrt(CA_2);
154 137707 : params.CB = sqrt(CB_2);
155 137707 : params.potential = libgrpp_shrink_potential(potential);
156 :
157 137707 : params.bra = libgrpp_shell_deep_copy(bra);
158 137707 : libgrpp_shell_shrink(params.bra);
159 137707 : libgrpp_shell_mult_normcoef(params.bra);
160 :
161 137707 : params.ket = libgrpp_shell_deep_copy(ket);
162 137707 : libgrpp_shell_shrink(params.ket);
163 137707 : libgrpp_shell_mult_normcoef(params.ket);
164 :
165 : /*
166 : * create radial grid
167 : */
168 137707 : radial_type2_grid_t *grid =
169 137707 : create_radial_type2_grid(lambda1_max, lambda2_max, n_max, ¶ms);
170 :
171 : /*
172 : * calculate radial integrals and store them into the table
173 : */
174 673203 : for (int lambda_1 = 0; lambda_1 <= lambda1_max; lambda_1++) {
175 2902163 : for (int lambda_2 = 0; lambda_2 <= lambda2_max; lambda_2++) {
176 9964982 : for (int n = 0; n <= n_max; n++) {
177 :
178 7598315 : int converged;
179 7598315 : double Q = calculate_radial_type2_integral(grid, n, lambda_1, lambda_2,
180 : 1e-16, &converged);
181 :
182 : // int dim1 = lambda1_max + 1;
183 7598315 : int dim2 = lambda2_max + 1;
184 7598315 : int dimn = n_max + 1;
185 7598315 : table->radial_integrals[dim2 * dimn * lambda_1 + dimn * lambda_2 + n] =
186 : Q;
187 : }
188 : }
189 : }
190 :
191 : /*
192 : * clean-up
193 : */
194 137707 : libgrpp_delete_potential(params.potential);
195 137707 : libgrpp_delete_shell(params.bra);
196 137707 : libgrpp_delete_shell(params.ket);
197 137707 : delete_radial_type2_grid(grid);
198 :
199 137707 : return table;
200 : }
201 :
202 : /**
203 : * destructor for the table of radial type 2 integrals
204 : */
205 149632 : void libgrpp_delete_radial_type2_integrals(radial_type2_table_t *table) {
206 149632 : free(table->radial_integrals);
207 149632 : free(table);
208 149632 : }
209 :
210 : /**
211 : * Returns radial integral at complex index (lambda1,lambda2,N)
212 : */
213 14370774 : double libgrpp_get_radial_type2_integral(radial_type2_table_t *table,
214 : int lambda1, int lambda2, int n) {
215 : // int lambda1_max = table->lambda1_max;
216 14370774 : int lambda2_max = table->lambda2_max;
217 14370774 : int n_max = table->n_max;
218 : // int dim1 = lambda1_max + 1;
219 14370774 : int dim2 = lambda2_max + 1;
220 14370774 : int dimn = n_max + 1;
221 :
222 14370774 : double Q =
223 14370774 : table->radial_integrals[dim2 * dimn * lambda1 + dimn * lambda2 + n];
224 :
225 14370774 : return Q;
226 : }
227 :
228 : /**
229 : * calculates type 2 radial integral T^N_{lambda1,lambda2}
230 : * for the two given contracted gaussian functions and the contracted potential
231 : */
232 7598315 : static double calculate_radial_type2_integral(radial_type2_grid_t *grid, int n,
233 : int lambda1, int lambda2,
234 : double tolerance,
235 : int *converged) {
236 7598315 : int nr = MIN_GRID;
237 :
238 7598315 : *converged = 0;
239 7598315 : double prev_sum = 0.0;
240 7598315 : double sum = 0.0;
241 :
242 7598315 : double *w = grid->w;
243 : // double *r = grid->r;
244 7598315 : double *pot_values = grid->rpp_values;
245 7598315 : double *F1 = grid->F1[lambda1];
246 7598315 : double *F2 = grid->F2[lambda2];
247 7598315 : double *r_N = grid->r_N[n];
248 :
249 : /*
250 : * first step: integral screening
251 : */
252 7598315 : double CA = grid->params->CA;
253 7598315 : double CB = grid->params->CB;
254 :
255 7598315 : double screened = 0.0;
256 7598315 : int screen_success = libgrpp_screening_radial_type2(
257 : lambda1, lambda2, n, CA * CA, CB * CB, grid->params->bra,
258 : grid->params->ket, grid->params->potential, &screened);
259 :
260 7598315 : if (screen_success == EXIT_SUCCESS && fabs(screened) < tolerance) {
261 3594781 : *converged = 1;
262 3594781 : return screened;
263 : }
264 :
265 : /*
266 : * second step: calculation on the smallest possible grid
267 : */
268 128113088 : for (int i = 0; i < nr; i++) {
269 124109554 : sum += w[i] * pot_values[i] * F1[i] * F2[i] * r_N[i];
270 : }
271 :
272 : /*
273 : * third step: adaptive integration, refinement of the result
274 : */
275 4350771 : do {
276 4350771 : int idx = nr;
277 4350771 : nr = 2 * nr + 1;
278 4350771 : if (nr > MAX_GRID) {
279 : break;
280 : }
281 :
282 4350170 : prev_sum = sum;
283 4350170 : sum = 0.5 * sum;
284 :
285 4350170 : expand_radial_type2_grid(grid, nr);
286 :
287 163783066 : for (int i = idx; i < nr; i++) {
288 159432896 : sum += w[i] * pot_values[i] * F1[i] * F2[i] * r_N[i];
289 : }
290 :
291 4350170 : if (screen_success == EXIT_SUCCESS &&
292 1194761 : (fabs(sum) / fabs(screened) < 0.001)) {
293 0 : *converged = 0;
294 0 : continue;
295 : }
296 :
297 4350170 : *converged = fabs(sum - prev_sum) <= tolerance;
298 :
299 4350170 : } while (!(*converged));
300 :
301 : return sum;
302 : }
303 :
304 : /**
305 : * Numerical integration on the Log3 grid
306 : */
307 : static radial_type2_grid_t *
308 137707 : create_radial_type2_grid(int lambda1_max, int lambda2_max, int n_max,
309 : radial_type2_params_t *params) {
310 137707 : radial_type2_grid_t *grid =
311 137707 : (radial_type2_grid_t *)calloc(1, sizeof(radial_type2_grid_t));
312 :
313 137707 : grid->nr = MIN_GRID;
314 137707 : grid->n_max = n_max;
315 137707 : grid->lambda1_max = lambda1_max;
316 137707 : grid->lambda2_max = lambda2_max;
317 137707 : grid->params = params;
318 :
319 137707 : grid->r = alloc_zeros_1d(MAX_GRID);
320 137707 : grid->w = alloc_zeros_1d(MAX_GRID);
321 137707 : grid->rpp_values = alloc_zeros_1d(MAX_GRID);
322 137707 : grid->F1 = alloc_zeros_2d(lambda1_max + 1, MAX_GRID);
323 137707 : grid->F2 = alloc_zeros_2d(lambda2_max + 1, MAX_GRID);
324 137707 : grid->r_N = alloc_zeros_2d(n_max + 1, MAX_GRID);
325 :
326 : // vectors 'k': k = - 2 * alpha * |CA|
327 137707 : double *bra_k = alloc_zeros_1d(params->bra->num_primitives);
328 137707 : double *ket_k = alloc_zeros_1d(params->ket->num_primitives);
329 137707 : calc_k_values(params->bra->num_primitives, params->bra->alpha, params->CA,
330 : bra_k);
331 137707 : calc_k_values(params->ket->num_primitives, params->ket->alpha, params->CB,
332 : ket_k);
333 :
334 : // initial set of pre-calculated points
335 137707 : int nr = grid->nr;
336 137707 : const double R = 5.0;
337 137707 : const double R3 = R * R * R;
338 :
339 4406624 : for (int i = 1; i <= nr; i++) {
340 4268917 : double xi = i / (nr + 1.0);
341 4268917 : double xi3 = xi * xi * xi;
342 4268917 : double ln_xi = log(1 - xi3);
343 4268917 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
344 4268917 : double ri = -R * ln_xi;
345 :
346 4268917 : grid->r[i - 1] = ri;
347 4268917 : grid->w[i - 1] = wi;
348 8537834 : grid->rpp_values[i - 1] =
349 4268917 : libgrpp_potential_value(grid->params->potential, ri);
350 16590084 : for (int n = 0; n <= n_max; n++) {
351 12321167 : grid->r_N[n][i - 1] = pow(ri, n);
352 : }
353 20869293 : for (int lambda1 = 0; lambda1 <= lambda1_max; lambda1++) {
354 16600376 : grid->F1[lambda1][i - 1] = radial_type2_integrand_fun_contracted(
355 : ri, lambda1, bra_k, params->CA, params->bra);
356 : }
357 20886033 : for (int lambda2 = 0; lambda2 <= lambda2_max; lambda2++) {
358 16617116 : grid->F2[lambda2][i - 1] = radial_type2_integrand_fun_contracted(
359 : ri, lambda2, ket_k, params->CB, params->ket);
360 : }
361 : }
362 :
363 137707 : free(bra_k);
364 137707 : free(ket_k);
365 :
366 137707 : return grid;
367 : }
368 :
369 : /**
370 : * constructs new radial grid points
371 : */
372 4350170 : static void expand_radial_type2_grid(radial_type2_grid_t *grid, int nr) {
373 4350170 : const double R = 5.0;
374 4350170 : const double R3 = R * R * R;
375 :
376 4350170 : if (nr > MAX_GRID) {
377 : return;
378 : }
379 :
380 4350170 : if (nr <= grid->nr) { // nothing to do
381 : return;
382 : }
383 :
384 159854 : radial_type2_params_t *params = grid->params;
385 :
386 : // vectors 'k': k = - 2 * alpha * |CA|
387 159854 : double *bra_k = alloc_zeros_1d(params->bra->num_primitives);
388 159854 : double *ket_k = alloc_zeros_1d(params->ket->num_primitives);
389 159854 : calc_k_values(params->bra->num_primitives, params->bra->alpha, params->CA,
390 : bra_k);
391 159854 : calc_k_values(params->ket->num_primitives, params->ket->alpha, params->CB,
392 : ket_k);
393 :
394 : // additional set of grid points
395 159854 : int idx = grid->nr;
396 11165774 : for (int i = 1; i <= nr; i += 2) {
397 11005920 : double xi = i / (nr + 1.0);
398 11005920 : double xi3 = xi * xi * xi;
399 11005920 : double ln_xi = log(1 - xi3);
400 11005920 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
401 11005920 : double ri = -R * ln_xi;
402 :
403 11005920 : grid->r[idx] = ri;
404 11005920 : grid->w[idx] = wi;
405 22011840 : grid->rpp_values[idx] =
406 11005920 : libgrpp_potential_value(grid->params->potential, ri);
407 :
408 40524288 : for (int n = 0; n <= grid->n_max; n++) {
409 29518368 : grid->r_N[n][idx] = pow(ri, n);
410 : }
411 :
412 45512640 : for (int lambda1 = 0; lambda1 <= grid->lambda1_max; lambda1++) {
413 34506720 : grid->F1[lambda1][idx] = radial_type2_integrand_fun_contracted(
414 34506720 : ri, lambda1, bra_k, grid->params->CA, params->bra);
415 : }
416 46238336 : for (int lambda2 = 0; lambda2 <= grid->lambda2_max; lambda2++) {
417 35232416 : grid->F2[lambda2][idx] = radial_type2_integrand_fun_contracted(
418 35232416 : ri, lambda2, ket_k, grid->params->CB, params->ket);
419 : }
420 :
421 11005920 : idx++;
422 : }
423 :
424 159854 : grid->nr = nr;
425 :
426 159854 : free(bra_k);
427 159854 : free(ket_k);
428 : }
429 :
430 : /**
431 : * deallocates memory used for the radial grid
432 : */
433 137707 : static void delete_radial_type2_grid(radial_type2_grid_t *grid) {
434 137707 : free(grid->r);
435 137707 : free(grid->w);
436 137707 : free(grid->rpp_values);
437 137707 : free_2d(grid->F1, grid->lambda1_max + 1);
438 137707 : free_2d(grid->F2, grid->lambda2_max + 1);
439 137707 : free_2d(grid->r_N, grid->n_max + 1);
440 137707 : free(grid);
441 137707 : }
442 :
443 : /**
444 : * Calculate the value of the integrand function
445 : */
446 : static double
447 102956628 : radial_type2_integrand_fun_contracted(double r, int lambda, double *k,
448 : double CA, libgrpp_shell_t *gauss_fun) {
449 102956628 : double F = 0.0;
450 102956628 : double r_CA_2 = (r - CA) * (r - CA);
451 :
452 102956628 : int nprim = gauss_fun->num_primitives;
453 102956628 : double *alpha = gauss_fun->alpha;
454 102956628 : double *coeffs = gauss_fun->coeffs;
455 :
456 205913256 : for (int i = 0; i < nprim; i++) {
457 102956628 : double power = -alpha[i] * r_CA_2;
458 102956628 : F += coeffs[i] * exp(power) *
459 102956628 : libgrpp_modified_bessel_scaled(lambda, k[i] * r);
460 : }
461 :
462 102956628 : return F;
463 : }
464 :
465 595122 : static void calc_k_values(int nprim, const double *alpha, double CA,
466 : double *k) {
467 1190244 : for (int i = 0; i < nprim; i++) {
468 595122 : k[i] = 2.0 * alpha[i] * CA;
469 : }
470 595122 : }
|
