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 1 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.
21 : * J. Comp. Chem. 27, 1009 (2006)
22 : * (see formulas (12) and (13) for radial integrals invoking contracted Gaussian
23 : functions and RPPs).
24 :
25 : * In contrast to type 2 integrals, the special case of type 1 integrals Bessel
26 : functions
27 : * cannot be factorized, and one cannot use contracted Gaussians directly
28 : * (and we have to use primitive Gaussians instead).
29 : * However, the RPP radial function still can be used as a whole in the
30 : integrand.
31 : *
32 : * The Log3 integration scheme used here is detailed in:
33 : * C.-K. Skylaris et al. An efficient method for calculating effective core
34 : potential integrals
35 : * which involve projection operators.
36 : * Chem. Phys. Lett. 296, 445 (1998)
37 : */
38 : #include <math.h>
39 : #include <stdlib.h>
40 : #ifndef M_PI
41 : #define M_PI 3.14159265358979323846
42 : #endif
43 :
44 : #include "grpp_radial_type1_integral.h"
45 : #include "libgrpp.h"
46 :
47 : #include "grpp_specfunc.h"
48 : #include "grpp_utils.h"
49 :
50 : #define MIN_GRID 2047
51 : #define MAX_GRID 10000
52 :
53 : typedef struct {
54 : double k;
55 : double alpha_A;
56 : double alpha_B;
57 : double CA_2;
58 : double CB_2;
59 : double prefactor;
60 : double (*potential)(double r, void *params);
61 : void *potential_params;
62 : } radial_type1_params_t;
63 :
64 : typedef struct {
65 : int nr;
66 : int n_max;
67 : int lambda_max;
68 : double *r;
69 : double *w;
70 : double *pot_values;
71 : double *gto_values;
72 : double **r_N;
73 : double **mod_bessel;
74 : radial_type1_params_t *params;
75 : } radial_type1_grid_t;
76 :
77 : static radial_type1_grid_t *
78 : create_radial_type1_grid(int lambda_max, int n_max,
79 : radial_type1_params_t *params);
80 :
81 : static void expand_radial_type1_grid(radial_type1_grid_t *grid, int nr);
82 :
83 : static void delete_radial_type1_grid(radial_type1_grid_t *grid);
84 :
85 : static double calculate_radial_type1_integral(radial_type1_grid_t *grid, int n,
86 : int lambda, double tolerance,
87 : int *converged);
88 :
89 0 : radial_type1_table_t *libgrpp_tabulate_radial_type1_integrals(
90 : int lambda_max, int n_max, double CA_2, double CB_2, double alpha_A,
91 : double alpha_B, double k, double prefactor,
92 : double (*potential)(double r, void *params), void *potential_params) {
93 0 : radial_type1_table_t *table;
94 0 : double const tolerance = libgrpp_params.radial_tolerance;
95 :
96 0 : table = (radial_type1_table_t *)calloc(1, sizeof(radial_type1_table_t));
97 0 : table->lambda_max = lambda_max;
98 0 : table->n_max = n_max;
99 0 : table->radial_integrals =
100 0 : (double *)calloc((lambda_max + 1) * (n_max + 1), sizeof(double));
101 :
102 0 : radial_type1_params_t params;
103 0 : params.CA_2 = CA_2;
104 0 : params.CB_2 = CB_2;
105 0 : params.alpha_A = alpha_A;
106 0 : params.alpha_B = alpha_B;
107 0 : params.k = k;
108 0 : params.prefactor = prefactor;
109 0 : params.potential = potential;
110 0 : params.potential_params = potential_params;
111 :
112 0 : radial_type1_grid_t *grid =
113 0 : create_radial_type1_grid(lambda_max, n_max, ¶ms);
114 :
115 0 : for (int lambda = 0; lambda <= lambda_max; lambda++) {
116 0 : for (int n = 0; n <= n_max; n++) {
117 :
118 0 : int converged;
119 0 : double Q = calculate_radial_type1_integral(grid, n, lambda, tolerance,
120 : &converged);
121 :
122 0 : table->radial_integrals[lambda * (lambda_max + 1) + n] = Q;
123 : }
124 : }
125 :
126 0 : delete_radial_type1_grid(grid);
127 :
128 0 : return table;
129 : }
130 :
131 0 : void libgrpp_delete_radial_type1_integrals(radial_type1_table_t *table) {
132 0 : free(table->radial_integrals);
133 0 : free(table);
134 0 : }
135 :
136 0 : double libgrpp_get_radial_type1_integral(radial_type1_table_t *table,
137 : int lambda, int n) {
138 0 : int lambda_max = table->lambda_max;
139 0 : return table->radial_integrals[lambda * (lambda_max + 1) + n];
140 : }
141 :
142 0 : static double radial_type1_integrand_fun(double r,
143 : radial_type1_params_t *params) {
144 0 : double alpha_A = params->alpha_A;
145 0 : double alpha_B = params->alpha_B;
146 0 : double k = params->k;
147 0 : double CA_2 = params->CA_2;
148 0 : double CB_2 = params->CB_2;
149 0 : double prefactor = params->prefactor;
150 :
151 0 : double power = k * r - (alpha_A + alpha_B) * r * r - alpha_A * CA_2 -
152 0 : alpha_B * CB_2; // + N * log(r);
153 :
154 0 : return prefactor * exp(power);
155 : }
156 :
157 : static radial_type1_grid_t *
158 0 : create_radial_type1_grid(int lambda_max, int n_max,
159 : radial_type1_params_t *params) {
160 0 : radial_type1_grid_t *grid =
161 0 : (radial_type1_grid_t *)calloc(1, sizeof(radial_type1_grid_t));
162 :
163 0 : grid->nr = MIN_GRID;
164 0 : grid->n_max = n_max;
165 0 : grid->lambda_max = lambda_max;
166 0 : grid->params = params;
167 :
168 0 : grid->r = (double *)calloc(MAX_GRID, sizeof(double));
169 0 : grid->w = (double *)calloc(MAX_GRID, sizeof(double));
170 0 : grid->pot_values = (double *)calloc(MAX_GRID, sizeof(double));
171 0 : grid->gto_values = (double *)calloc(MAX_GRID, sizeof(double));
172 0 : grid->r_N = alloc_zeros_2d(n_max + 1, MAX_GRID);
173 0 : grid->mod_bessel = alloc_zeros_2d(lambda_max + 1, MAX_GRID);
174 :
175 : // initial set of pre-calculated points
176 0 : int nr = grid->nr;
177 0 : const double R = 5.0;
178 0 : const double R3 = R * R * R;
179 :
180 0 : for (int i = 1; i <= nr; i++) {
181 0 : double xi = i / (nr + 1.0);
182 0 : double xi3 = xi * xi * xi;
183 0 : double ln_xi = log(1 - xi3);
184 0 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
185 0 : double ri = -R * ln_xi;
186 :
187 0 : grid->r[i - 1] = ri;
188 0 : grid->w[i - 1] = wi;
189 0 : grid->pot_values[i - 1] = params->potential(ri, params->potential_params);
190 0 : grid->gto_values[i - 1] = radial_type1_integrand_fun(ri, params);
191 :
192 0 : for (int lambda = 0; lambda <= lambda_max; lambda++) {
193 0 : grid->mod_bessel[lambda][i - 1] =
194 0 : libgrpp_modified_bessel_scaled(lambda, ri * params->k);
195 : }
196 :
197 0 : for (int n = 0; n <= n_max; n++) {
198 0 : grid->r_N[n][i - 1] = pow(ri, n);
199 : }
200 : }
201 :
202 0 : return grid;
203 : }
204 :
205 0 : static void delete_radial_type1_grid(radial_type1_grid_t *grid) {
206 0 : free(grid->r);
207 0 : free(grid->w);
208 0 : free(grid->pot_values);
209 0 : free(grid->gto_values);
210 0 : free_2d(grid->r_N, grid->n_max + 1);
211 0 : free_2d(grid->mod_bessel, grid->lambda_max + 1);
212 0 : free(grid);
213 0 : }
214 :
215 0 : static void expand_radial_type1_grid(radial_type1_grid_t *grid, int nr) {
216 0 : const double R = 5.0;
217 0 : const double R3 = R * R * R;
218 :
219 0 : if (nr > MAX_GRID) {
220 : return;
221 : }
222 :
223 0 : if (nr <= grid->nr) { // nothing to do
224 : return;
225 : }
226 :
227 : int idx = grid->nr;
228 0 : for (int i = 1; i <= nr; i += 2) {
229 0 : double xi = i / (nr + 1.0);
230 0 : double xi3 = xi * xi * xi;
231 0 : double ln_xi = log(1 - xi3);
232 0 : double wi = 3 * R3 * xi * xi * ln_xi * ln_xi / ((1 - xi3) * (nr + 1.0));
233 0 : double ri = -R * ln_xi;
234 :
235 0 : grid->r[idx] = ri;
236 0 : grid->w[idx] = wi;
237 0 : grid->pot_values[idx] =
238 0 : grid->params->potential(ri, grid->params->potential_params);
239 0 : grid->gto_values[idx] = radial_type1_integrand_fun(ri, grid->params);
240 :
241 0 : for (int lambda = 0; lambda <= grid->lambda_max; lambda++) {
242 0 : double kr = grid->params->k * ri;
243 0 : double bessel = libgrpp_modified_bessel_scaled(lambda, kr);
244 0 : grid->mod_bessel[lambda][idx] = bessel;
245 : }
246 :
247 0 : for (int n = 0; n <= grid->n_max; n++) {
248 0 : grid->r_N[n][idx] = pow(ri, n);
249 : }
250 0 : idx++;
251 : }
252 :
253 0 : grid->nr = nr;
254 : }
255 :
256 0 : static double calculate_radial_type1_integral(radial_type1_grid_t *grid, int n,
257 : int lambda, double tolerance,
258 : int *converged) {
259 0 : int nr = MIN_GRID;
260 :
261 0 : *converged = 0;
262 0 : double prev_sum = 0.0;
263 0 : double sum = 0.0;
264 :
265 0 : double *w = grid->w;
266 : // double *r = grid->r;
267 0 : double *pot_values = grid->pot_values;
268 0 : double *gto_values = grid->gto_values;
269 0 : double *r_N = grid->r_N[n];
270 0 : double *mod_bessel = grid->mod_bessel[lambda];
271 :
272 : /*
273 : * first step: screening of an integral
274 : */
275 : /*double screened = 0.0;
276 : int screened_success = screening_radial_type1(
277 : lambda,
278 : n,
279 : grid->params->CA_2,
280 : grid->params->CB_2,
281 : grid->params->alpha_A,
282 : grid->params->alpha_B,
283 : grid->params->k,
284 : grid->params->prefactor,
285 : grid->params->potential_params,
286 : &screened
287 : );
288 :
289 : if (screened_success == EXIT_SUCCESS && fabs(screened) < tolerance) {
290 : *converged = 1;
291 : return screened;
292 : }*/
293 :
294 : /*
295 : * second step: calculation on the smallest possible grid
296 : */
297 0 : for (int i = 0; i < nr; i++) {
298 0 : sum += w[i] * pot_values[i] * gto_values[i] * r_N[i] * mod_bessel[i];
299 : }
300 :
301 : /*
302 : * third step: adaptive integration, refinement of the result
303 : */
304 0 : do {
305 0 : int idx = nr;
306 0 : nr = 2 * nr + 1;
307 :
308 0 : if (nr > MAX_GRID) {
309 : break;
310 : }
311 :
312 0 : prev_sum = sum;
313 0 : sum = 0.5 * sum;
314 :
315 0 : expand_radial_type1_grid(grid, nr);
316 :
317 0 : for (int i = idx; i < nr; i++) {
318 0 : sum += w[i] * pot_values[i] * gto_values[i] * r_N[i] * mod_bessel[i];
319 : }
320 :
321 : /*if (screened_success == EXIT_SUCCESS && (fabs(sum) / fabs(screened) <
322 : 0.001)) { *converged = 0; continue;
323 : }*/
324 :
325 0 : *converged = fabs(sum - prev_sum) <= tolerance;
326 :
327 0 : } while (!(*converged));
328 :
329 0 : return sum;
330 : }
|
