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 : * Calculation of nuclear attraction integrals.
17 : *
18 : * For the point charge nuclear model the recursive Obara-Saika scheme is used
19 : * to calculate nuclear attraction integrals. For details, see
20 : * T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-Structure Theory,
21 : * John Wiley & Sons Ltd, 2000.
22 : * Chapter 9.10.1, "The Obara-Saika scheme for one-electron Coulomb integrals"
23 : *
24 : * For the other three models,
25 : * - uniformly charged ball
26 : * - Gaussian nucleus
27 : * - Fermi nucleus,
28 : * the scheme is actually the same as for the type 1 (radially-local) ECP
29 : * integrals. Electrostatic potential V(r) induced by the finite nuclear charge
30 : * distribution is integrated numerically on a radial grid.
31 : */
32 :
33 : #include <assert.h>
34 : #include <math.h>
35 : #include <stdio.h>
36 : #include <stdlib.h>
37 : #include <string.h>
38 :
39 : #include "grpp_norm_gaussian.h"
40 : #include "grpp_nuclear_models.h"
41 : #include "grpp_utils.h"
42 : #include "libgrpp.h"
43 :
44 : #ifndef M_PI
45 : #define M_PI 3.14159265358979323846
46 : #endif
47 :
48 : void libgrpp_evaluate_radially_local_potential_integral_primitive_gaussians(
49 : double *A, int n_cart_A, int *cart_list_A, double alpha_A, double *B,
50 : int n_cart_B, int *cart_list_B, double alpha_B, double *C,
51 : double (*potential)(double r, void *params), void *potential_params,
52 : double *matrix);
53 :
54 : void libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
55 : libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
56 : double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix);
57 :
58 : static double wrapper_coulomb_potential_point(double r, void *params);
59 :
60 : static double wrapper_coulomb_potential_ball(double r, void *params);
61 :
62 : static double wrapper_coulomb_potential_gaussian(double r, void *params);
63 :
64 : static double wrapper_coulomb_potential_fermi(double r, void *params);
65 :
66 : static double wrapper_coulomb_potential_fermi_bubble(double r, void *params);
67 :
68 : /**
69 : * Calculates nuclear attraction integral between two shells
70 : * represented by contracted Gaussian functions.
71 : *
72 : * nuclear model should be one of:
73 : * LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE
74 : * LIBGRPP_NUCLEAR_MODEL_CHARGED_BALL
75 : * LIBGRPP_NUCLEAR_MODEL_GAUSSIAN
76 : * LIBGRPP_NUCLEAR_MODEL_FERMI
77 : * LIBGRPP_NUCLEAR_MODEL_FERMI_BUBBLE
78 : * LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE_NUMERICAL
79 : */
80 0 : void libgrpp_nuclear_attraction_integrals(libgrpp_shell_t *shell_A,
81 : libgrpp_shell_t *shell_B,
82 : double *charge_origin, int charge,
83 : int nuclear_model,
84 : double *model_params,
85 : double *coulomb_matrix) {
86 0 : assert(libgrpp_is_initialized());
87 :
88 0 : int size_A = libgrpp_get_shell_size(shell_A);
89 0 : int size_B = libgrpp_get_shell_size(shell_B);
90 :
91 0 : double *buf = calloc(size_A * size_B, sizeof(double));
92 :
93 0 : memset(coulomb_matrix, 0, size_A * size_B * sizeof(double));
94 :
95 : // loop over primitives in contractions
96 0 : for (int i = 0; i < shell_A->num_primitives; i++) {
97 0 : double coef_A_i = shell_A->coeffs[i];
98 :
99 0 : for (int j = 0; j < shell_B->num_primitives; j++) {
100 0 : double coef_B_j = shell_B->coeffs[j];
101 :
102 0 : if (fabs(coef_A_i * coef_B_j) < 1e-13) {
103 0 : continue;
104 : }
105 :
106 0 : if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE) {
107 :
108 : // use code for RPP type-1 integrals with RPP exponent = 0.0
109 0 : libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
110 0 : shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j],
111 : charge_origin, 0.0, buf);
112 :
113 0 : for (int k = 0; k < size_A * size_B; k++) {
114 0 : buf[k] *= (-1) * charge;
115 : }
116 0 : } else if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_CHARGED_BALL ||
117 : nuclear_model == LIBGRPP_NUCLEAR_MODEL_GAUSSIAN ||
118 : nuclear_model == LIBGRPP_NUCLEAR_MODEL_FERMI ||
119 0 : nuclear_model == LIBGRPP_NUCLEAR_MODEL_FERMI_BUBBLE ||
120 : nuclear_model ==
121 : LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE_NUMERICAL) {
122 :
123 0 : double params[10];
124 0 : params[0] = charge;
125 :
126 : /*
127 : * choose nuclear model
128 : */
129 0 : double (*electrostatic_potential_fun)(double, void *) = NULL;
130 :
131 0 : if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE_NUMERICAL) {
132 : // printf("charge distribution: point\n");
133 : electrostatic_potential_fun = wrapper_coulomb_potential_point;
134 0 : } else if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_CHARGED_BALL) {
135 0 : params[1] = model_params[0]; // R_rms
136 0 : electrostatic_potential_fun = wrapper_coulomb_potential_ball;
137 0 : } else if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_GAUSSIAN) {
138 0 : params[1] = model_params[0]; // R_rms
139 0 : electrostatic_potential_fun = wrapper_coulomb_potential_gaussian;
140 0 : } else if (nuclear_model == LIBGRPP_NUCLEAR_MODEL_FERMI) {
141 0 : params[1] = model_params[0]; // c
142 0 : params[2] = model_params[1]; // a
143 0 : electrostatic_potential_fun = wrapper_coulomb_potential_fermi;
144 : } else {
145 0 : params[1] = model_params[0]; // c
146 0 : params[2] = model_params[1]; // a
147 0 : params[3] = model_params[2]; // k
148 0 : electrostatic_potential_fun = wrapper_coulomb_potential_fermi_bubble;
149 : }
150 :
151 : /*
152 : * calculate integrals for the shell pair
153 : */
154 0 : libgrpp_evaluate_radially_local_potential_integral_primitive_gaussians(
155 0 : shell_A->origin, size_A, shell_A->cart_list, shell_A->alpha[i],
156 0 : shell_B->origin, size_B, shell_B->cart_list, shell_B->alpha[j],
157 : charge_origin, electrostatic_potential_fun, params, buf);
158 : } else {
159 0 : printf("LIBGRPP: unknown finite nuclear charge distribution model!\n");
160 0 : exit(0);
161 : }
162 :
163 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf, coulomb_matrix);
164 : }
165 : }
166 :
167 0 : free(buf);
168 0 : }
169 :
170 : /**
171 : * Calculates nuclear attraction integral between two shells
172 : * represented by contracted Gaussian functions for the electrostatic potential
173 : * generated by the point charge.
174 : */
175 0 : void libgrpp_nuclear_attraction_integrals_point_charge(libgrpp_shell_t *shell_A,
176 : libgrpp_shell_t *shell_B,
177 : double *charge_origin,
178 : int charge,
179 : double *coulomb_matrix) {
180 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, charge,
181 : LIBGRPP_NUCLEAR_MODEL_POINT_CHARGE, NULL,
182 : coulomb_matrix);
183 0 : }
184 :
185 : /**
186 : * Calculates nuclear attraction integral between two shells
187 : * represented by contracted Gaussian functions for the electrostatic potential
188 : * generated by the charged ball.
189 : *
190 : * r_rms stands for the root mean square radius (in bohrs)
191 : */
192 0 : void libgrpp_nuclear_attraction_integrals_charged_ball(libgrpp_shell_t *shell_A,
193 : libgrpp_shell_t *shell_B,
194 : double *charge_origin,
195 : int charge, double r_rms,
196 : double *coulomb_matrix) {
197 0 : double params[10];
198 0 : params[0] = charge;
199 0 : params[1] = r_rms;
200 :
201 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, charge,
202 : LIBGRPP_NUCLEAR_MODEL_CHARGED_BALL,
203 : params, coulomb_matrix);
204 0 : }
205 :
206 : /**
207 : * Calculates nuclear attraction integral between two shells
208 : * represented by contracted Gaussian functions for the electrostatic potential
209 : * generated by the Gaussian distribution.
210 : *
211 : * r_rms stands for the root mean square radius (in bohrs)
212 : */
213 0 : void libgrpp_nuclear_attraction_integrals_gaussian_model(
214 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *charge_origin,
215 : int charge, double r_rms, double *coulomb_matrix) {
216 0 : double params[10];
217 0 : params[0] = charge;
218 0 : params[1] = r_rms;
219 :
220 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, charge,
221 : LIBGRPP_NUCLEAR_MODEL_GAUSSIAN, params,
222 : coulomb_matrix);
223 0 : }
224 :
225 : /**
226 : * Calculates nuclear attraction integral between two shells
227 : * represented by contracted Gaussian functions for the electrostatic potential
228 : * generated by the Fermi distribution.
229 : *
230 : * Model parameters 'c' and 'a' must be given in bohrs.
231 : */
232 0 : void libgrpp_nuclear_attraction_integrals_fermi_model(
233 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *charge_origin,
234 : int charge, double fermi_param_c, double fermi_param_a,
235 : double *coulomb_matrix) {
236 0 : double params[10];
237 0 : params[0] = charge;
238 0 : params[1] = fermi_param_c;
239 0 : params[2] = fermi_param_a;
240 :
241 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, charge,
242 : LIBGRPP_NUCLEAR_MODEL_FERMI, params,
243 : coulomb_matrix);
244 0 : }
245 :
246 : /**
247 : * Calculates nuclear attraction integral between two shells
248 : * represented by contracted Gaussian functions for the electrostatic potential
249 : * generated by the "Fermi + bubble" distribution.
250 : *
251 : * Model parameters 'c' and 'a' must be given in bohrs.
252 : * The 'k' constant is dimensionless.
253 : */
254 0 : void libgrpp_nuclear_attraction_integrals_fermi_bubble_model(
255 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *charge_origin,
256 : int charge, double param_c, double param_a, double param_k,
257 : double *coulomb_matrix) {
258 0 : double params[10];
259 0 : params[0] = charge;
260 0 : params[1] = param_c;
261 0 : params[2] = param_a;
262 0 : params[3] = param_k;
263 :
264 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, charge,
265 : LIBGRPP_NUCLEAR_MODEL_FERMI_BUBBLE,
266 : params, coulomb_matrix);
267 0 : }
268 :
269 : /**
270 : * wrappers for charge distribution functions.
271 : * are used to provide a unified interface to radially-local potentials.
272 : * the 'params' argument is unpacked, then the specific routines are invoked.
273 : */
274 :
275 0 : static double wrapper_coulomb_potential_point(double r, void *params) {
276 0 : double Z = ((double *)params)[0];
277 :
278 0 : return libgrpp_coulomb_potential_point(r, Z);
279 : }
280 :
281 0 : static double wrapper_coulomb_potential_ball(double r, void *params) {
282 0 : double Z = ((double *)params)[0];
283 0 : double R_rms = ((double *)params)[1];
284 :
285 0 : return libgrpp_coulomb_potential_ball(r, Z, R_rms);
286 : }
287 :
288 0 : static double wrapper_coulomb_potential_gaussian(double r, void *params) {
289 0 : double Z = ((double *)params)[0];
290 0 : double R_rms = ((double *)params)[1];
291 :
292 0 : return libgrpp_coulomb_potential_gaussian(r, Z, R_rms);
293 : }
294 :
295 0 : static double wrapper_coulomb_potential_fermi(double r, void *params) {
296 0 : double Z = ((double *)params)[0];
297 0 : double c = ((double *)params)[1];
298 0 : double a = ((double *)params)[2];
299 :
300 0 : return libgrpp_coulomb_potential_fermi(r, Z, c, a);
301 : }
302 :
303 0 : static double wrapper_coulomb_potential_fermi_bubble(double r, void *params) {
304 0 : double Z = ((double *)params)[0];
305 0 : double c = ((double *)params)[1];
306 0 : double a = ((double *)params)[2];
307 0 : double k = ((double *)params)[3];
308 :
309 0 : return libgrpp_coulomb_potential_fermi_bubble(r, Z, c, a, k);
310 : }
|
