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 194176 : 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 194176 : radial_type2_table_t *table;
106 194176 : table = (radial_type2_table_t *)calloc(1, sizeof(radial_type2_table_t));
107 194176 : table->lambda1_max = lambda1_max;
108 194176 : table->lambda2_max = lambda2_max;
109 194176 : table->n_max = n_max;
110 194176 : table->radial_integrals = (double *)calloc(
111 194176 : (lambda1_max + 1) * (lambda2_max + 1) * (n_max + 1), sizeof(double));
112 :
113 : /*
114 : * the special case of one-center RPP integrals
115 : */
116 194176 : if (CA_2 < 1e-14 && CB_2 < 1e-14) {
117 :
118 75144 : for (int i = 0; i < bra->num_primitives; i++) {
119 37572 : double alpha_A = bra->alpha[i];
120 75144 : double coef_i =
121 37572 : bra->coeffs[i] * libgrpp_gaussian_norm_factor(bra->L, 0, 0, alpha_A);
122 :
123 75144 : for (int j = 0; j < ket->num_primitives; j++) {
124 37572 : double alpha_B = ket->alpha[j];
125 75144 : double coef_j = ket->coeffs[j] *
126 37572 : libgrpp_gaussian_norm_factor(ket->L, 0, 0, alpha_B);
127 :
128 201824 : for (int k = 0; k < potential->num_primitives; k++) {
129 164252 : double eta = potential->alpha[k];
130 164252 : int ni = potential->powers[k];
131 164252 : double coef_k = potential->coeffs[k];
132 :
133 164252 : double p = alpha_A + alpha_B + eta;
134 164252 : double factor = coef_i * coef_j * coef_k;
135 :
136 668304 : for (int n = 0; n <= n_max; n++) {
137 :
138 504052 : double val_ijk = libgrpp_gaussian_integral(ni + n, p);
139 504052 : table->radial_integrals[n] += factor * val_ijk;
140 504052 : ;
141 : }
142 : }
143 : }
144 : }
145 :
146 : return table;
147 : }
148 :
149 : /*
150 : * for numerical integration on the grid
151 : */
152 156604 : radial_type2_params_t params;
153 156604 : params.CA = sqrt(CA_2);
154 156604 : params.CB = sqrt(CB_2);
155 156604 : params.potential = libgrpp_shrink_potential(potential);
156 :
157 156604 : params.bra = libgrpp_shell_deep_copy(bra);
158 156604 : libgrpp_shell_shrink(params.bra);
159 156604 : libgrpp_shell_mult_normcoef(params.bra);
160 :
161 156604 : params.ket = libgrpp_shell_deep_copy(ket);
162 156604 : libgrpp_shell_shrink(params.ket);
163 156604 : libgrpp_shell_mult_normcoef(params.ket);
164 :
165 : /*
166 : * create radial grid
167 : */
168 156604 : radial_type2_grid_t *grid =
169 156604 : 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 742629 : for (int lambda_1 = 0; lambda_1 <= lambda1_max; lambda_1++) {
175 3106590 : for (int lambda_2 = 0; lambda_2 <= lambda2_max; lambda_2++) {
176 10562174 : for (int n = 0; n <= n_max; n++) {
177 :
178 8041609 : int converged;
179 8041609 : double Q = calculate_radial_type2_integral(grid, n, lambda_1, lambda_2,
180 : 1e-16, &converged);
181 :
182 : // int dim1 = lambda1_max + 1;
183 8041609 : int dim2 = lambda2_max + 1;
184 8041609 : int dimn = n_max + 1;
185 8041609 : table->radial_integrals[dim2 * dimn * lambda_1 + dimn * lambda_2 + n] =
186 : Q;
187 : }
188 : }
189 : }
190 :
191 : /*
192 : * clean-up
193 : */
194 156604 : libgrpp_delete_potential(params.potential);
195 156604 : libgrpp_delete_shell(params.bra);
196 156604 : libgrpp_delete_shell(params.ket);
197 156604 : delete_radial_type2_grid(grid);
198 :
199 156604 : return table;
200 : }
201 :
202 : /**
203 : * destructor for the table of radial type 2 integrals
204 : */
205 194176 : void libgrpp_delete_radial_type2_integrals(radial_type2_table_t *table) {
206 194176 : free(table->radial_integrals);
207 194176 : free(table);
208 194176 : }
209 :
210 : /**
211 : * Returns radial integral at complex index (lambda1,lambda2,N)
212 : */
213 16859719 : 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 16859719 : int lambda2_max = table->lambda2_max;
217 16859719 : int n_max = table->n_max;
218 : // int dim1 = lambda1_max + 1;
219 16859719 : int dim2 = lambda2_max + 1;
220 16859719 : int dimn = n_max + 1;
221 :
222 16859719 : double Q =
223 16859719 : table->radial_integrals[dim2 * dimn * lambda1 + dimn * lambda2 + n];
224 :
225 16859719 : 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 8041609 : 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 8041609 : int nr = MIN_GRID;
237 :
238 8041609 : *converged = 0;
239 8041609 : double prev_sum = 0.0;
240 8041609 : double sum = 0.0;
241 :
242 8041609 : double *w = grid->w;
243 : // double *r = grid->r;
244 8041609 : double *pot_values = grid->rpp_values;
245 8041609 : double *F1 = grid->F1[lambda1];
246 8041609 : double *F2 = grid->F2[lambda2];
247 8041609 : double *r_N = grid->r_N[n];
248 :
249 : /*
250 : * first step: integral screening
251 : */
252 8041609 : double CA = grid->params->CA;
253 8041609 : double CB = grid->params->CB;
254 :
255 8041609 : double screened = 0.0;
256 8041609 : 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 8041609 : if (screen_success == EXIT_SUCCESS && fabs(screened) < tolerance) {
261 3877417 : *converged = 1;
262 3877417 : return screened;
263 : }
264 :
265 : /*
266 : * second step: calculation on the smallest possible grid
267 : */
268 133254144 : for (int i = 0; i < nr; i++) {
269 129089952 : 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 4708977 : do {
276 4708977 : int idx = nr;
277 4708977 : nr = 2 * nr + 1;
278 4708977 : if (nr > MAX_GRID) {
279 : break;
280 : }
281 :
282 4705102 : prev_sum = sum;
283 4705102 : sum = 0.5 * sum;
284 :
285 4705102 : expand_radial_type2_grid(grid, nr);
286 :
287 226072526 : for (int i = idx; i < nr; i++) {
288 221367424 : sum += w[i] * pot_values[i] * F1[i] * F2[i] * r_N[i];
289 : }
290 :
291 4705102 : if (screen_success == EXIT_SUCCESS &&
292 1549693 : (fabs(sum) / fabs(screened) < 0.001)) {
293 0 : *converged = 0;
294 0 : continue;
295 : }
296 :
297 4705102 : *converged = fabs(sum - prev_sum) <= tolerance;
298 :
299 4705102 : } while (!(*converged));
300 :
301 : return sum;
302 : }
303 :
304 : /**
305 : * Numerical integration on the Log3 grid
306 : */
307 : static radial_type2_grid_t *
308 156604 : create_radial_type2_grid(int lambda1_max, int lambda2_max, int n_max,
309 : radial_type2_params_t *params) {
310 156604 : radial_type2_grid_t *grid =
311 156604 : (radial_type2_grid_t *)calloc(1, sizeof(radial_type2_grid_t));
312 :
313 156604 : grid->nr = MIN_GRID;
314 156604 : grid->n_max = n_max;
315 156604 : grid->lambda1_max = lambda1_max;
316 156604 : grid->lambda2_max = lambda2_max;
317 156604 : grid->params = params;
318 :
319 156604 : grid->r = alloc_zeros_1d(MAX_GRID);
320 156604 : grid->w = alloc_zeros_1d(MAX_GRID);
321 156604 : grid->rpp_values = alloc_zeros_1d(MAX_GRID);
322 156604 : grid->F1 = alloc_zeros_2d(lambda1_max + 1, MAX_GRID);
323 156604 : grid->F2 = alloc_zeros_2d(lambda2_max + 1, MAX_GRID);
324 156604 : grid->r_N = alloc_zeros_2d(n_max + 1, MAX_GRID);
325 :
326 : // vectors 'k': k = - 2 * alpha * |CA|
327 156604 : double *bra_k = alloc_zeros_1d(params->bra->num_primitives);
328 156604 : double *ket_k = alloc_zeros_1d(params->ket->num_primitives);
329 156604 : calc_k_values(params->bra->num_primitives, params->bra->alpha, params->CA,
330 : bra_k);
331 156604 : calc_k_values(params->ket->num_primitives, params->ket->alpha, params->CB,
332 : ket_k);
333 :
334 : // initial set of pre-calculated points
335 156604 : int nr = grid->nr;
336 156604 : const double R = 5.0;
337 156604 : const double R3 = R * R * R;
338 :
339 5011328 : for (int i = 1; i <= nr; i++) {
340 4854724 : double xi = i / (nr + 1.0);
341 4854724 : double xi3 = xi * xi * xi;
342 4854724 : double ln_xi = log(1 - xi3);
343 4854724 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
344 4854724 : double ri = -R * ln_xi;
345 :
346 4854724 : grid->r[i - 1] = ri;
347 4854724 : grid->w[i - 1] = wi;
348 9709448 : grid->rpp_values[i - 1] =
349 4854724 : libgrpp_potential_value(grid->params->potential, ri);
350 18623436 : for (int n = 0; n <= n_max; n++) {
351 13768712 : grid->r_N[n][i - 1] = pow(ri, n);
352 : }
353 23021499 : for (int lambda1 = 0; lambda1 <= lambda1_max; lambda1++) {
354 18166775 : grid->F1[lambda1][i - 1] = radial_type2_integrand_fun_contracted(
355 : ri, lambda1, bra_k, params->CA, params->bra);
356 : }
357 23110407 : for (int lambda2 = 0; lambda2 <= lambda2_max; lambda2++) {
358 18255683 : grid->F2[lambda2][i - 1] = radial_type2_integrand_fun_contracted(
359 : ri, lambda2, ket_k, params->CB, params->ket);
360 : }
361 : }
362 :
363 156604 : free(bra_k);
364 156604 : free(ket_k);
365 :
366 156604 : return grid;
367 : }
368 :
369 : /**
370 : * constructs new radial grid points
371 : */
372 4705102 : static void expand_radial_type2_grid(radial_type2_grid_t *grid, int nr) {
373 4705102 : const double R = 5.0;
374 4705102 : const double R3 = R * R * R;
375 :
376 4705102 : if (nr > MAX_GRID) {
377 : return;
378 : }
379 :
380 4705102 : if (nr <= grid->nr) { // nothing to do
381 : return;
382 : }
383 :
384 209777 : radial_type2_params_t *params = grid->params;
385 :
386 : // vectors 'k': k = - 2 * alpha * |CA|
387 209777 : double *bra_k = alloc_zeros_1d(params->bra->num_primitives);
388 209777 : double *ket_k = alloc_zeros_1d(params->ket->num_primitives);
389 209777 : calc_k_values(params->bra->num_primitives, params->bra->alpha, params->CA,
390 : bra_k);
391 209777 : calc_k_values(params->ket->num_primitives, params->ket->alpha, params->CB,
392 : ket_k);
393 :
394 : // additional set of grid points
395 209777 : int idx = grid->nr;
396 34271601 : for (int i = 1; i <= nr; i += 2) {
397 34061824 : double xi = i / (nr + 1.0);
398 34061824 : double xi3 = xi * xi * xi;
399 34061824 : double ln_xi = log(1 - xi3);
400 34061824 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
401 34061824 : double ri = -R * ln_xi;
402 :
403 34061824 : grid->r[idx] = ri;
404 34061824 : grid->w[idx] = wi;
405 68123648 : grid->rpp_values[idx] =
406 34061824 : libgrpp_potential_value(grid->params->potential, ri);
407 :
408 119638432 : for (int n = 0; n <= grid->n_max; n++) {
409 85576608 : grid->r_N[n][idx] = pow(ri, n);
410 : }
411 :
412 127136160 : for (int lambda1 = 0; lambda1 <= grid->lambda1_max; lambda1++) {
413 93074336 : grid->F1[lambda1][idx] = radial_type2_integrand_fun_contracted(
414 93074336 : ri, lambda1, bra_k, grid->params->CA, params->bra);
415 : }
416 132910336 : for (int lambda2 = 0; lambda2 <= grid->lambda2_max; lambda2++) {
417 98848512 : grid->F2[lambda2][idx] = radial_type2_integrand_fun_contracted(
418 98848512 : ri, lambda2, ket_k, grid->params->CB, params->ket);
419 : }
420 :
421 34061824 : idx++;
422 : }
423 :
424 209777 : grid->nr = nr;
425 :
426 209777 : free(bra_k);
427 209777 : free(ket_k);
428 : }
429 :
430 : /**
431 : * deallocates memory used for the radial grid
432 : */
433 156604 : static void delete_radial_type2_grid(radial_type2_grid_t *grid) {
434 156604 : free(grid->r);
435 156604 : free(grid->w);
436 156604 : free(grid->rpp_values);
437 156604 : free_2d(grid->F1, grid->lambda1_max + 1);
438 156604 : free_2d(grid->F2, grid->lambda2_max + 1);
439 156604 : free_2d(grid->r_N, grid->n_max + 1);
440 156604 : free(grid);
441 156604 : }
442 :
443 : /**
444 : * Calculate the value of the integrand function
445 : */
446 : static double
447 228345306 : radial_type2_integrand_fun_contracted(double r, int lambda, double *k,
448 : double CA, libgrpp_shell_t *gauss_fun) {
449 228345306 : double F = 0.0;
450 228345306 : double r_CA_2 = (r - CA) * (r - CA);
451 :
452 228345306 : int nprim = gauss_fun->num_primitives;
453 228345306 : double *alpha = gauss_fun->alpha;
454 228345306 : double *coeffs = gauss_fun->coeffs;
455 :
456 456690612 : for (int i = 0; i < nprim; i++) {
457 228345306 : double power = -alpha[i] * r_CA_2;
458 228345306 : F += coeffs[i] * exp(power) *
459 228345306 : libgrpp_modified_bessel_scaled(lambda, k[i] * r);
460 : }
461 :
462 228345306 : return F;
463 : }
464 :
465 732762 : static void calc_k_values(int nprim, const double *alpha, double CA,
466 : double *k) {
467 1465524 : for (int i = 0; i < nprim; i++) {
468 732762 : k[i] = 2.0 * alpha[i] * CA;
469 : }
470 732762 : }
|
