(git:ed6f26b)
Loading...
Searching...
No Matches
grpp_type1_mcmurchie_davidson.c
Go to the documentation of this file.
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 * Implementation of the McMurchie-Davidson (MMD) recurrence relations
17 * for radially local RPP integrals.
18 *
19 * Operators to be integrated are:
20 * n = 2: e^{-ar^2}
21 * n = 1: e^{-ar^2}/r^1
22 * n = 0: e^{-ar^2}/r^2
23 *
24 * General description of the MMD scheme is given in:
25 * (1) T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-Structure
26 * Theory, John Wiley & Sons Ltd, 2000. (2) L. E. McMurchie, E. R. Davidson,
27 * One- and two-electron integrals over cartesian gaussian functions.
28 * J. Comput. Phys. 26(2), 218 (1978).
29 * doi: 10.1016/0021-9991(78)90092-X
30 *
31 * More on the evaluation of integrals over the 1/r^2 operator:
32 * (1) J. 0. Jensen, A. H. Cameri, C. P. Vlahacos, D. Zeroka, H. F. Hameka, C.
33 * N. Merrow, Evaluation of one-electron integrals for arbitrary operators V(r)
34 * over cartesian Gaussians: Application to inverse-square distance and Yukawa
35 * operators. J. Comput. Chem. 14(8), 986 (1993). doi: 10.1002/jcc.540140814 (2)
36 * B. Gao, A. J. Thorvaldsen, K. Ruud, GEN1INT: A unified procedure for the
37 * evaluation of one-electron integrals over Gaussian basis functions and their
38 * geometric derivatives. Int. J. Quantum Chem. 111(4), 858 (2011).
39 * doi: 10.1002/qua.22886
40 */
41#include <assert.h>
42#include <math.h>
43#include <stdlib.h>
44#include <string.h>
45
46#include "grpp_type1_mcmurchie_davidson.h"
47
48#ifndef M_PI
49#define M_PI 3.14159265358979323846
50#endif
51
52#include "grpp_norm_gaussian.h"
53#include "grpp_specfunc.h"
54#include "grpp_utils.h"
55#include "libgrpp.h"
56
57#define LMAX LIBGRPP_MAX_BASIS_L
58
59/*
60 * temporary data for the McMurchie-Davidson algorithm
61 */
62struct mmd_data {
63 double A[3];
64 double B[3];
65 double C[3];
66 double Q[3];
67 double q;
68 double K_abc[3];
69 double R_QC_2;
70 double boys[100];
71 double E[3][LMAX][LMAX][2 * LMAX];
72 double R[2 * LMAX][2 * LMAX][2 * LMAX][2 * LMAX];
73};
74
75/*
76 * functions used below in the file
77 */
78
79static void evaluate_rpp_type1_mmd_n2_primitive_shell_pair(
80 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
81 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix);
82
83void libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
84 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
85 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix);
86
87static void evaluate_rpp_type1_mmd_n0_primitive_shell_pair(
88 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
89 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix);
90
91static void setup_E_array(struct mmd_data *data, int L_A, int L_B);
92
93static void setup_R_array(struct mmd_data *data, int L_A, int L_B);
94
95static void setup_G_array(struct mmd_data *data, int L_A, int L_B);
96
97/**
98 * General interface for the McMurchie-Davidson algorithm for integrals
99 * over the radially local RPP operator.
100 */
101void libgrpp_type1_integrals_mcmurchie_davidson_1978(
102 libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *origin_C,
103 double alpha_C, int ecp_power, double *rpp_matrix) {
104 assert((ecp_power == 0) || (ecp_power == 1) || (ecp_power == 2));
105
106 int size_A = libgrpp_get_shell_size(shell_A);
107 int size_B = libgrpp_get_shell_size(shell_B);
108
109 double *buf = calloc(size_A * size_B, sizeof(double));
110
111 memset(rpp_matrix, 0, size_A * size_B * sizeof(double));
112
113 /*
114 * loop over primitives in contractions
115 */
116 for (int i = 0; i < shell_A->num_primitives; i++) {
117 double coef_A_i = shell_A->coeffs[i];
118 if (fabs(coef_A_i) < LIBGRPP_ZERO_THRESH) {
119 continue;
120 }
121
122 for (int j = 0; j < shell_B->num_primitives; j++) {
123 double coef_B_j = shell_B->coeffs[j];
124 if (fabs(coef_B_j) < LIBGRPP_ZERO_THRESH) {
125 continue;
126 }
127
128 if (ecp_power == 2) {
129 evaluate_rpp_type1_mmd_n2_primitive_shell_pair(
130 shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
131 alpha_C, buf);
132 } else if (ecp_power == 1) {
133 libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
134 shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
135 alpha_C, buf);
136 } else if (ecp_power == 0) {
137 evaluate_rpp_type1_mmd_n0_primitive_shell_pair(
138 shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
139 alpha_C, buf);
140 }
141
142 libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf, rpp_matrix);
143 }
144 }
145
146 free(buf);
147}
148
149/**
150 * Integrals of the operator: e^{-ar^2}/r^2
151 */
152static void evaluate_rpp_type1_mmd_n0_primitive_shell_pair(
153 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
154 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix) {
155 double a = alpha_A;
156 double b = alpha_B;
157 double c = rpp_alpha;
158
159 double *A = shell_A->origin;
160 double *B = shell_B->origin;
161 double *C = rpp_origin;
162
163 double q = a + b + c;
164 double mu_AB = a * b / q;
165 double mu_AC = a * c / q;
166 double mu_BC = b * c / q;
167
168 struct mmd_data data;
169
170 for (int i = 0; i < 3; i++) {
171 data.A[i] = A[i];
172 data.B[i] = B[i];
173 data.C[i] = C[i];
174 data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
175
176 double X_AB = A[i] - B[i];
177 double X_AC = A[i] - C[i];
178 double X_BC = B[i] - C[i];
179
180 double K_ab = exp(-mu_AB * X_AB * X_AB);
181 double K_ac = exp(-mu_AC * X_AC * X_AC);
182 double K_bc = exp(-mu_BC * X_BC * X_BC);
183
184 data.K_abc[i] = K_ab * K_ac * K_bc;
185 }
186
187 data.q = q;
188 data.R_QC_2 = distance_squared(data.Q, data.C);
189 int L_A = shell_A->L;
190 int L_B = shell_B->L;
191 int Nmax = L_A + L_B;
192 libgrpp_gfun_values(q * data.R_QC_2, Nmax, data.boys);
193
194 /*
195 * setup E array
196 */
197 setup_E_array(&data, L_A, L_B);
198
199 /*
200 * setup R array
201 */
202 setup_G_array(&data, L_A, L_B);
203
204 /*
205 * loop over cartesian functions inside the shells
206 */
207 int size_A = libgrpp_get_shell_size(shell_A);
208 int size_B = libgrpp_get_shell_size(shell_B);
209 double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
210 double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
211
212 for (int m = 0; m < size_A; m++) {
213 for (int n = 0; n < size_B; n++) {
214 int n_A = shell_A->cart_list[3 * m + 0];
215 int l_A = shell_A->cart_list[3 * m + 1];
216 int m_A = shell_A->cart_list[3 * m + 2];
217 int n_B = shell_B->cart_list[3 * n + 0];
218 int l_B = shell_B->cart_list[3 * n + 1];
219 int m_B = shell_B->cart_list[3 * n + 2];
220
221 double s = 0.0;
222
223 for (int t = 0; t <= n_A + n_B; t++) {
224 double Ex = data.E[0][n_A][n_B][t];
225
226 for (int u = 0; u <= l_A + l_B; u++) {
227 double Ey = data.E[1][l_A][l_B][u];
228
229 for (int v = 0; v <= m_A + m_B; v++) {
230 double Ez = data.E[2][m_A][m_B][v];
231
232 double R_tuv = data.R[t][u][v][0];
233
234 s += Ex * Ey * Ez * R_tuv;
235 }
236 }
237 }
238
239 s *= N_A * N_B * 2.0 * pow(M_PI, 1.5) / sqrt(q);
240 rpp_matrix[m * size_B + n] = s;
241 }
242 }
243}
244
245/**
246 * Integrals of the operator: e^{-ar^2}/r
247 */
248void libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
249 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
250 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix) {
251 double a = alpha_A;
252 double b = alpha_B;
253 double c = rpp_alpha;
254
255 double *A = shell_A->origin;
256 double *B = shell_B->origin;
257 double *C = rpp_origin;
258
259 double q = a + b + c;
260 double mu_AB = a * b / q;
261 double mu_AC = a * c / q;
262 double mu_BC = b * c / q;
263
264 struct mmd_data data;
265
266 for (int i = 0; i < 3; i++) {
267 data.A[i] = A[i];
268 data.B[i] = B[i];
269 data.C[i] = C[i];
270 data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
271
272 double X_AB = A[i] - B[i];
273 double X_AC = A[i] - C[i];
274 double X_BC = B[i] - C[i];
275
276 double K_ab = exp(-mu_AB * X_AB * X_AB);
277 double K_ac = exp(-mu_AC * X_AC * X_AC);
278 double K_bc = exp(-mu_BC * X_BC * X_BC);
279
280 data.K_abc[i] = K_ab * K_ac * K_bc;
281 }
282
283 data.q = q;
284 data.R_QC_2 = distance_squared(data.Q, data.C);
285 int L_A = shell_A->L;
286 int L_B = shell_B->L;
287 int Nmax = L_A + L_B;
288
289 libgrpp_boys_values(q * data.R_QC_2, Nmax, data.boys);
290
291 /*
292 * setup E array
293 */
294 setup_E_array(&data, L_A, L_B);
295
296 /*
297 * setup R array
298 */
299 setup_R_array(&data, L_A, L_B);
300
301 /*
302 * loop over cartesian functions inside the shells
303 */
304 int size_A = libgrpp_get_shell_size(shell_A);
305 int size_B = libgrpp_get_shell_size(shell_B);
306 double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
307 double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
308
309 for (int m = 0; m < size_A; m++) {
310 for (int n = 0; n < size_B; n++) {
311 int n_A = shell_A->cart_list[3 * m + 0];
312 int l_A = shell_A->cart_list[3 * m + 1];
313 int m_A = shell_A->cart_list[3 * m + 2];
314 int n_B = shell_B->cart_list[3 * n + 0];
315 int l_B = shell_B->cart_list[3 * n + 1];
316 int m_B = shell_B->cart_list[3 * n + 2];
317
318 double s = 0.0;
319
320 for (int t = 0; t <= n_A + n_B; t++) {
321 double Ex = data.E[0][n_A][n_B][t];
322
323 for (int u = 0; u <= l_A + l_B; u++) {
324 double Ey = data.E[1][l_A][l_B][u];
325
326 for (int v = 0; v <= m_A + m_B; v++) {
327 double Ez = data.E[2][m_A][m_B][v];
328
329 double R_tuv = data.R[t][u][v][0];
330
331 s += Ex * Ey * Ez * R_tuv;
332 }
333 }
334 }
335
336 s *= N_A * N_B * 2.0 * M_PI / q;
337 rpp_matrix[m * size_B + n] = s;
338 }
339 }
340}
341
342/**
343 * Integrals of the operator: e^{-ar^2}
344 */
345static void evaluate_rpp_type1_mmd_n2_primitive_shell_pair(
346 libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B,
347 double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix) {
348 double a = alpha_A;
349 double b = alpha_B;
350 double c = rpp_alpha;
351
352 double *A = shell_A->origin;
353 double *B = shell_B->origin;
354 double *C = rpp_origin;
355
356 double q = a + b + c;
357 double mu_AB = a * b / q;
358 double mu_AC = a * c / q;
359 double mu_BC = b * c / q;
360
361 struct mmd_data data;
362
363 for (int i = 0; i < 3; i++) {
364 data.A[i] = A[i];
365 data.B[i] = B[i];
366 data.C[i] = C[i];
367 data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
368
369 double X_AB = A[i] - B[i];
370 double X_AC = A[i] - C[i];
371 double X_BC = B[i] - C[i];
372
373 double K_ab = exp(-mu_AB * X_AB * X_AB);
374 double K_ac = exp(-mu_AC * X_AC * X_AC);
375 double K_bc = exp(-mu_BC * X_BC * X_BC);
376
377 data.K_abc[i] = K_ab * K_ac * K_bc;
378 }
379
380 data.q = q;
381 data.R_QC_2 = distance_squared(data.Q, data.C);
382
383 int L_A = shell_A->L;
384 int L_B = shell_B->L;
385
386 /*
387 * setup E array
388 */
389 setup_E_array(&data, L_A, L_B);
390
391 /*
392 * loop over cartesian functions inside the shells
393 */
394 int size_A = libgrpp_get_shell_size(shell_A);
395 int size_B = libgrpp_get_shell_size(shell_B);
396 double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
397 double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
398
399 for (int m = 0; m < size_A; m++) {
400 for (int n = 0; n < size_B; n++) {
401 int n_A = shell_A->cart_list[3 * m + 0];
402 int l_A = shell_A->cart_list[3 * m + 1];
403 int m_A = shell_A->cart_list[3 * m + 2];
404 int n_B = shell_B->cart_list[3 * n + 0];
405 int l_B = shell_B->cart_list[3 * n + 1];
406 int m_B = shell_B->cart_list[3 * n + 2];
407
408 double E_ij_0 = data.E[0][n_A][n_B][0];
409 double E_kl_0 = data.E[1][l_A][l_B][0];
410 double E_mn_0 = data.E[2][m_A][m_B][0];
411
412 double s = N_A * N_B * E_ij_0 * E_kl_0 * E_mn_0 * pow(M_PI / q, 1.5);
413
414 rpp_matrix[m * size_B + n] = s;
415 }
416 }
417}
418
419/**
420 * Calculation of the R_{tuv}^N auxiliary integrals for the 1/r operator
421 */
422static void setup_R_array(struct mmd_data *data, int L_A, int L_B) {
423 double q = data->q;
424 double X_QC = data->Q[0] - data->C[0];
425 double Y_QC = data->Q[1] - data->C[1];
426 double Z_QC = data->Q[2] - data->C[2];
427
428 int nmax = L_A + L_B;
429 int L_SUM = L_A + L_B;
430
431 for (int t = 0; t <= L_SUM; t++) {
432 for (int u = 0; u <= L_SUM; u++) {
433 for (int v = 0; v <= L_SUM; v++) {
434
435 if (t + u + v > L_SUM) {
436 continue;
437 }
438
439 for (int n = 0; n <= nmax - (t + u + v); n++) {
440
441 double val = 0.0;
442
443 if (t + u + v == 0) {
444 val = pow(-2.0 * q, n) * data->boys[n];
445 } else if (t + u == 0) {
446 if (v > 1) {
447 val += (v - 1) * data->R[t][u][v - 2][n + 1];
448 }
449 val += Z_QC * data->R[t][u][v - 1][n + 1];
450 } else if (t == 0) {
451 if (u > 1) {
452 val += (u - 1) * data->R[t][u - 2][v][n + 1];
453 }
454 val += Y_QC * data->R[t][u - 1][v][n + 1];
455 } else {
456 if (t > 1) {
457 val += (t - 1) * data->R[t - 2][u][v][n + 1];
458 }
459 val += X_QC * data->R[t - 1][u][v][n + 1];
460 }
461
462 data->R[t][u][v][n] = val;
463 }
464 }
465 }
466 }
467}
468
469/**
470 * Calculation of the R_{tuv}^N auxiliary integrals for the 1/r^2 operator
471 */
472static void setup_G_array(struct mmd_data *data, int L_A, int L_B) {
473 double q = data->q;
474 double X_QC = data->Q[0] - data->C[0];
475 double Y_QC = data->Q[1] - data->C[1];
476 double Z_QC = data->Q[2] - data->C[2];
477
478 int nmax = L_A + L_B;
479 int L_SUM = L_A + L_B;
480
481 for (int t = 0; t <= L_SUM; t++) {
482 for (int u = 0; u <= L_SUM; u++) {
483 for (int v = 0; v <= L_SUM; v++) {
484
485 if (t + u + v > L_SUM) {
486 continue;
487 }
488
489 for (int n = 0; n <= nmax - (t + u + v); n++) {
490
491 double val = 0.0;
492
493 if (t + u + v == 0) {
494 val = pow(2.0 * q, n) * data->boys[n];
495 } else if (t + u == 0) {
496 if (v > 1) {
497 val += (v - 1) * (data->R[t][u][v - 2][n + 1] -
498 2.0 * q * data->R[t][u][v - 2][n]);
499 }
500 val += Z_QC * (data->R[t][u][v - 1][n + 1] -
501 2.0 * q * data->R[t][u][v - 1][n]);
502 } else if (t == 0) {
503 if (u > 1) {
504 val += (u - 1) * (data->R[t][u - 2][v][n + 1] -
505 2.0 * q * data->R[t][u - 2][v][n]);
506 }
507 val += Y_QC * (data->R[t][u - 1][v][n + 1] -
508 2.0 * q * data->R[t][u - 1][v][n]);
509 } else {
510 if (t > 1) {
511 val += (t - 1) * (data->R[t - 2][u][v][n + 1] -
512 2.0 * q * data->R[t - 2][u][v][n]);
513 }
514 val += X_QC * (data->R[t - 1][u][v][n + 1] -
515 2.0 * q * data->R[t - 1][u][v][n]);
516 }
517
518 data->R[t][u][v][n] = val;
519 }
520 }
521 }
522 }
523}
524
525/**
526 * Calculates E^{ij}_t coefficients in the MMD scheme.
527 */
528static void setup_E_array(struct mmd_data *data, int L_A, int L_B) {
529 double q = data->q;
530
531 for (int coord = 0; coord < 3; coord++) {
532 for (int i = 0; i <= L_A; i++) {
533 for (int j = 0; j <= L_B; j++) {
534
535 for (int t = 0; t <= i + j; t++) {
536
537 if (t == 0) {
538 if (i == 0 && j == 0) {
539 data->E[coord][0][0][0] = data->K_abc[coord];
540 continue;
541 } else if (i == 1 && j == 0) {
542 double X_QA = data->Q[coord] - data->A[coord];
543 data->E[coord][1][0][0] = X_QA * data->E[coord][0][0][0];
544 continue;
545 } else if (i == 0 && j == 1) {
546 double X_QB = data->Q[coord] - data->B[coord];
547 data->E[coord][0][1][0] = X_QB * data->E[coord][0][0][0];
548 continue;
549 } else if (i == 0) {
550 double X_QB = data->Q[coord] - data->B[coord];
551 data->E[coord][i][j][0] = X_QB * data->E[coord][i][j - 1][0] +
552 data->E[coord][i][j - 1][1];
553 continue;
554 } else {
555 double X_QA = data->Q[coord] - data->A[coord];
556 data->E[coord][i][j][0] = X_QA * data->E[coord][i - 1][j][0] +
557 data->E[coord][i - 1][j][1];
558 continue;
559 }
560 } else {
561 double E_ijt = 0.0;
562 double factor = 1.0 / (2.0 * q * t);
563
564 if (i > 0) {
565 E_ijt += factor * i * data->E[coord][i - 1][j][t - 1];
566 }
567 if (j > 0) {
568 E_ijt += factor * j * data->E[coord][i][j - 1][t - 1];
569 }
570
571 data->E[coord][i][j][t] = E_ijt;
572 }
573 }
574 }
575 }
576 }
577}
i
static void const int const int i
Definition grid_cpu_collint.h:38
M_PI
#define M_PI
Definition grpp_angular_integrals.c:31
libgrpp_gaussian_norm_factor
double libgrpp_gaussian_norm_factor(int n, int l, int m, double alpha)
Definition grpp_norm_gaussian.c:27
grpp_norm_gaussian.h
libgrpp_get_shell_size
int libgrpp_get_shell_size(libgrpp_shell_t *shell)
Definition grpp_shell.c:98
grpp_specfunc.h
libgrpp_boys_values
void libgrpp_boys_values(double x, int nmax, double *b)
Definition grpp_specfunc_boys.c:53
libgrpp_gfun_values
void libgrpp_gfun_values(double x, int nmax, double *g)
Definition grpp_specfunc_gfun.c:45
setup_G_array
static void setup_G_array(struct mmd_data *data, int L_A, int L_B)
Definition grpp_type1_mcmurchie_davidson.c:472
evaluate_rpp_type1_mmd_n0_primitive_shell_pair
static void evaluate_rpp_type1_mmd_n0_primitive_shell_pair(libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B, double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix)
Definition grpp_type1_mcmurchie_davidson.c:152
LMAX
#define LMAX
Definition grpp_type1_mcmurchie_davidson.c:57
setup_E_array
static void setup_E_array(struct mmd_data *data, int L_A, int L_B)
Definition grpp_type1_mcmurchie_davidson.c:528
libgrpp_type1_integrals_mcmurchie_davidson_1978
void libgrpp_type1_integrals_mcmurchie_davidson_1978(libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *origin_C, double alpha_C, int ecp_power, double *rpp_matrix)
Definition grpp_type1_mcmurchie_davidson.c:101
libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair
void libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B, double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix)
Definition grpp_type1_mcmurchie_davidson.c:248
evaluate_rpp_type1_mmd_n2_primitive_shell_pair
static void evaluate_rpp_type1_mmd_n2_primitive_shell_pair(libgrpp_shell_t *shell_A, double alpha_A, libgrpp_shell_t *shell_B, double alpha_B, double *rpp_origin, double rpp_alpha, double *rpp_matrix)
Definition grpp_type1_mcmurchie_davidson.c:345
setup_R_array
static void setup_R_array(struct mmd_data *data, int L_A, int L_B)
Definition grpp_type1_mcmurchie_davidson.c:422
grpp_type1_mcmurchie_davidson.h
distance_squared
double distance_squared(double *A, double *B)
Definition grpp_utils.c:68
libgrpp_daxpy
void libgrpp_daxpy(int n, double a, double *x, double *y)
Definition grpp_utils.c:46
grpp_utils.h
LIBGRPP_ZERO_THRESH
#define LIBGRPP_ZERO_THRESH
Definition libgrpp_types.h:56
libgrpp_shell_t
Definition libgrpp_types.h:20
libgrpp_shell_t::num_primitives
int num_primitives
Definition libgrpp_types.h:24
libgrpp_shell_t::cart_list
int * cart_list
Definition libgrpp_types.h:23
libgrpp_shell_t::origin
double origin[3]
Definition libgrpp_types.h:27
libgrpp_shell_t::coeffs
double * coeffs
Definition libgrpp_types.h:25
libgrpp_shell_t::L
int L
Definition libgrpp_types.h:21
libgrpp_shell_t::alpha
double * alpha
Definition libgrpp_types.h:26
mmd_data
Definition grpp_type1_mcmurchie_davidson.c:62
mmd_data::Q
double Q[3]
Definition grpp_type1_mcmurchie_davidson.c:66
mmd_data::A
double A[3]
Definition grpp_type1_mcmurchie_davidson.c:63
mmd_data::E
double E[3][LMAX][LMAX][2 *LMAX]
Definition grpp_type1_mcmurchie_davidson.c:71
mmd_data::R
double R[2 *LMAX][2 *LMAX][2 *LMAX][2 *LMAX]
Definition grpp_type1_mcmurchie_davidson.c:72
mmd_data::K_abc
double K_abc[3]
Definition grpp_type1_mcmurchie_davidson.c:68
mmd_data::R_QC_2
double R_QC_2
Definition grpp_type1_mcmurchie_davidson.c:69
mmd_data::boys
double boys[100]
Definition grpp_type1_mcmurchie_davidson.c:70
mmd_data::B
double B[3]
Definition grpp_type1_mcmurchie_davidson.c:64
mmd_data::q
double q
Definition grpp_type1_mcmurchie_davidson.c:67
mmd_data::C
double C[3]
Definition grpp_type1_mcmurchie_davidson.c:65