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 : * 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 : */
62 : struct 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 :
79 : static 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 :
83 : void 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 :
87 : static 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 :
91 : static void setup_E_array(struct mmd_data *data, int L_A, int L_B);
92 :
93 : static void setup_R_array(struct mmd_data *data, int L_A, int L_B);
94 :
95 : static 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 : */
101 0 : void 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 0 : assert((ecp_power == 0) || (ecp_power == 1) || (ecp_power == 2));
105 :
106 0 : int size_A = libgrpp_get_shell_size(shell_A);
107 0 : int size_B = libgrpp_get_shell_size(shell_B);
108 :
109 0 : double *buf = calloc(size_A * size_B, sizeof(double));
110 :
111 0 : memset(rpp_matrix, 0, size_A * size_B * sizeof(double));
112 :
113 : /*
114 : * loop over primitives in contractions
115 : */
116 0 : for (int i = 0; i < shell_A->num_primitives; i++) {
117 0 : double coef_A_i = shell_A->coeffs[i];
118 0 : if (fabs(coef_A_i) < LIBGRPP_ZERO_THRESH) {
119 0 : continue;
120 : }
121 :
122 0 : for (int j = 0; j < shell_B->num_primitives; j++) {
123 0 : double coef_B_j = shell_B->coeffs[j];
124 0 : if (fabs(coef_B_j) < LIBGRPP_ZERO_THRESH) {
125 0 : continue;
126 : }
127 :
128 0 : if (ecp_power == 2) {
129 0 : evaluate_rpp_type1_mmd_n2_primitive_shell_pair(
130 0 : shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
131 : alpha_C, buf);
132 0 : } else if (ecp_power == 1) {
133 0 : libgrpp_evaluate_rpp_type1_mmd_n1_primitive_shell_pair(
134 0 : shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
135 : alpha_C, buf);
136 0 : } else if (ecp_power == 0) {
137 0 : evaluate_rpp_type1_mmd_n0_primitive_shell_pair(
138 0 : shell_A, shell_A->alpha[i], shell_B, shell_B->alpha[j], origin_C,
139 : alpha_C, buf);
140 : }
141 :
142 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf, rpp_matrix);
143 : }
144 : }
145 :
146 0 : free(buf);
147 0 : }
148 :
149 : /**
150 : * Integrals of the operator: e^{-ar^2}/r^2
151 : */
152 0 : static 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 0 : double a = alpha_A;
156 0 : double b = alpha_B;
157 0 : double c = rpp_alpha;
158 :
159 0 : double *A = shell_A->origin;
160 0 : double *B = shell_B->origin;
161 0 : double *C = rpp_origin;
162 :
163 0 : double q = a + b + c;
164 0 : double mu_AB = a * b / q;
165 0 : double mu_AC = a * c / q;
166 0 : double mu_BC = b * c / q;
167 :
168 0 : struct mmd_data data;
169 :
170 0 : for (int i = 0; i < 3; i++) {
171 0 : data.A[i] = A[i];
172 0 : data.B[i] = B[i];
173 0 : data.C[i] = C[i];
174 0 : data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
175 :
176 0 : double X_AB = A[i] - B[i];
177 0 : double X_AC = A[i] - C[i];
178 0 : double X_BC = B[i] - C[i];
179 :
180 0 : double K_ab = exp(-mu_AB * X_AB * X_AB);
181 0 : double K_ac = exp(-mu_AC * X_AC * X_AC);
182 0 : double K_bc = exp(-mu_BC * X_BC * X_BC);
183 :
184 0 : data.K_abc[i] = K_ab * K_ac * K_bc;
185 : }
186 :
187 0 : data.q = q;
188 0 : data.R_QC_2 = distance_squared(data.Q, data.C);
189 0 : int L_A = shell_A->L;
190 0 : int L_B = shell_B->L;
191 0 : int Nmax = L_A + L_B;
192 0 : libgrpp_gfun_values(q * data.R_QC_2, Nmax, data.boys);
193 :
194 : /*
195 : * setup E array
196 : */
197 0 : setup_E_array(&data, L_A, L_B);
198 :
199 : /*
200 : * setup R array
201 : */
202 0 : setup_G_array(&data, L_A, L_B);
203 :
204 : /*
205 : * loop over cartesian functions inside the shells
206 : */
207 0 : int size_A = libgrpp_get_shell_size(shell_A);
208 0 : int size_B = libgrpp_get_shell_size(shell_B);
209 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
210 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
211 :
212 0 : for (int m = 0; m < size_A; m++) {
213 0 : for (int n = 0; n < size_B; n++) {
214 0 : int n_A = shell_A->cart_list[3 * m + 0];
215 0 : int l_A = shell_A->cart_list[3 * m + 1];
216 0 : int m_A = shell_A->cart_list[3 * m + 2];
217 0 : int n_B = shell_B->cart_list[3 * n + 0];
218 0 : int l_B = shell_B->cart_list[3 * n + 1];
219 0 : int m_B = shell_B->cart_list[3 * n + 2];
220 :
221 0 : double s = 0.0;
222 :
223 0 : for (int t = 0; t <= n_A + n_B; t++) {
224 0 : double Ex = data.E[0][n_A][n_B][t];
225 :
226 0 : for (int u = 0; u <= l_A + l_B; u++) {
227 0 : double Ey = data.E[1][l_A][l_B][u];
228 :
229 0 : for (int v = 0; v <= m_A + m_B; v++) {
230 0 : double Ez = data.E[2][m_A][m_B][v];
231 :
232 0 : double R_tuv = data.R[t][u][v][0];
233 :
234 0 : s += Ex * Ey * Ez * R_tuv;
235 : }
236 : }
237 : }
238 :
239 0 : s *= N_A * N_B * 2.0 * pow(M_PI, 1.5) / sqrt(q);
240 0 : rpp_matrix[m * size_B + n] = s;
241 : }
242 : }
243 0 : }
244 :
245 : /**
246 : * Integrals of the operator: e^{-ar^2}/r
247 : */
248 0 : void 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 0 : double a = alpha_A;
252 0 : double b = alpha_B;
253 0 : double c = rpp_alpha;
254 :
255 0 : double *A = shell_A->origin;
256 0 : double *B = shell_B->origin;
257 0 : double *C = rpp_origin;
258 :
259 0 : double q = a + b + c;
260 0 : double mu_AB = a * b / q;
261 0 : double mu_AC = a * c / q;
262 0 : double mu_BC = b * c / q;
263 :
264 0 : struct mmd_data data;
265 :
266 0 : for (int i = 0; i < 3; i++) {
267 0 : data.A[i] = A[i];
268 0 : data.B[i] = B[i];
269 0 : data.C[i] = C[i];
270 0 : data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
271 :
272 0 : double X_AB = A[i] - B[i];
273 0 : double X_AC = A[i] - C[i];
274 0 : double X_BC = B[i] - C[i];
275 :
276 0 : double K_ab = exp(-mu_AB * X_AB * X_AB);
277 0 : double K_ac = exp(-mu_AC * X_AC * X_AC);
278 0 : double K_bc = exp(-mu_BC * X_BC * X_BC);
279 :
280 0 : data.K_abc[i] = K_ab * K_ac * K_bc;
281 : }
282 :
283 0 : data.q = q;
284 0 : data.R_QC_2 = distance_squared(data.Q, data.C);
285 0 : int L_A = shell_A->L;
286 0 : int L_B = shell_B->L;
287 0 : int Nmax = L_A + L_B;
288 :
289 0 : libgrpp_boys_values(q * data.R_QC_2, Nmax, data.boys);
290 :
291 : /*
292 : * setup E array
293 : */
294 0 : setup_E_array(&data, L_A, L_B);
295 :
296 : /*
297 : * setup R array
298 : */
299 0 : setup_R_array(&data, L_A, L_B);
300 :
301 : /*
302 : * loop over cartesian functions inside the shells
303 : */
304 0 : int size_A = libgrpp_get_shell_size(shell_A);
305 0 : int size_B = libgrpp_get_shell_size(shell_B);
306 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
307 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
308 :
309 0 : for (int m = 0; m < size_A; m++) {
310 0 : for (int n = 0; n < size_B; n++) {
311 0 : int n_A = shell_A->cart_list[3 * m + 0];
312 0 : int l_A = shell_A->cart_list[3 * m + 1];
313 0 : int m_A = shell_A->cart_list[3 * m + 2];
314 0 : int n_B = shell_B->cart_list[3 * n + 0];
315 0 : int l_B = shell_B->cart_list[3 * n + 1];
316 0 : int m_B = shell_B->cart_list[3 * n + 2];
317 :
318 0 : double s = 0.0;
319 :
320 0 : for (int t = 0; t <= n_A + n_B; t++) {
321 0 : double Ex = data.E[0][n_A][n_B][t];
322 :
323 0 : for (int u = 0; u <= l_A + l_B; u++) {
324 0 : double Ey = data.E[1][l_A][l_B][u];
325 :
326 0 : for (int v = 0; v <= m_A + m_B; v++) {
327 0 : double Ez = data.E[2][m_A][m_B][v];
328 :
329 0 : double R_tuv = data.R[t][u][v][0];
330 :
331 0 : s += Ex * Ey * Ez * R_tuv;
332 : }
333 : }
334 : }
335 :
336 0 : s *= N_A * N_B * 2.0 * M_PI / q;
337 0 : rpp_matrix[m * size_B + n] = s;
338 : }
339 : }
340 0 : }
341 :
342 : /**
343 : * Integrals of the operator: e^{-ar^2}
344 : */
345 0 : static 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 0 : double a = alpha_A;
349 0 : double b = alpha_B;
350 0 : double c = rpp_alpha;
351 :
352 0 : double *A = shell_A->origin;
353 0 : double *B = shell_B->origin;
354 0 : double *C = rpp_origin;
355 :
356 0 : double q = a + b + c;
357 0 : double mu_AB = a * b / q;
358 0 : double mu_AC = a * c / q;
359 0 : double mu_BC = b * c / q;
360 :
361 0 : struct mmd_data data;
362 :
363 0 : for (int i = 0; i < 3; i++) {
364 0 : data.A[i] = A[i];
365 0 : data.B[i] = B[i];
366 0 : data.C[i] = C[i];
367 0 : data.Q[i] = (a * A[i] + b * B[i] + c * C[i]) / q;
368 :
369 0 : double X_AB = A[i] - B[i];
370 0 : double X_AC = A[i] - C[i];
371 0 : double X_BC = B[i] - C[i];
372 :
373 0 : double K_ab = exp(-mu_AB * X_AB * X_AB);
374 0 : double K_ac = exp(-mu_AC * X_AC * X_AC);
375 0 : double K_bc = exp(-mu_BC * X_BC * X_BC);
376 :
377 0 : data.K_abc[i] = K_ab * K_ac * K_bc;
378 : }
379 :
380 0 : data.q = q;
381 0 : data.R_QC_2 = distance_squared(data.Q, data.C);
382 :
383 0 : int L_A = shell_A->L;
384 0 : int L_B = shell_B->L;
385 :
386 : /*
387 : * setup E array
388 : */
389 0 : setup_E_array(&data, L_A, L_B);
390 :
391 : /*
392 : * loop over cartesian functions inside the shells
393 : */
394 0 : int size_A = libgrpp_get_shell_size(shell_A);
395 0 : int size_B = libgrpp_get_shell_size(shell_B);
396 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
397 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
398 :
399 0 : for (int m = 0; m < size_A; m++) {
400 0 : for (int n = 0; n < size_B; n++) {
401 0 : int n_A = shell_A->cart_list[3 * m + 0];
402 0 : int l_A = shell_A->cart_list[3 * m + 1];
403 0 : int m_A = shell_A->cart_list[3 * m + 2];
404 0 : int n_B = shell_B->cart_list[3 * n + 0];
405 0 : int l_B = shell_B->cart_list[3 * n + 1];
406 0 : int m_B = shell_B->cart_list[3 * n + 2];
407 :
408 0 : double E_ij_0 = data.E[0][n_A][n_B][0];
409 0 : double E_kl_0 = data.E[1][l_A][l_B][0];
410 0 : double E_mn_0 = data.E[2][m_A][m_B][0];
411 :
412 0 : double s = N_A * N_B * E_ij_0 * E_kl_0 * E_mn_0 * pow(M_PI / q, 1.5);
413 :
414 0 : rpp_matrix[m * size_B + n] = s;
415 : }
416 : }
417 0 : }
418 :
419 : /**
420 : * Calculation of the R_{tuv}^N auxiliary integrals for the 1/r operator
421 : */
422 0 : static void setup_R_array(struct mmd_data *data, int L_A, int L_B) {
423 0 : double q = data->q;
424 0 : double X_QC = data->Q[0] - data->C[0];
425 0 : double Y_QC = data->Q[1] - data->C[1];
426 0 : double Z_QC = data->Q[2] - data->C[2];
427 :
428 0 : int nmax = L_A + L_B;
429 0 : int L_SUM = L_A + L_B;
430 :
431 0 : for (int t = 0; t <= L_SUM; t++) {
432 0 : for (int u = 0; u <= L_SUM; u++) {
433 0 : for (int v = 0; v <= L_SUM; v++) {
434 :
435 0 : if (t + u + v > L_SUM) {
436 0 : continue;
437 : }
438 :
439 0 : for (int n = 0; n <= nmax - (t + u + v); n++) {
440 :
441 0 : double val = 0.0;
442 :
443 0 : if (t + u + v == 0) {
444 0 : val = pow(-2.0 * q, n) * data->boys[n];
445 0 : } else if (t + u == 0) {
446 0 : if (v > 1) {
447 0 : val += (v - 1) * data->R[t][u][v - 2][n + 1];
448 : }
449 0 : val += Z_QC * data->R[t][u][v - 1][n + 1];
450 0 : } else if (t == 0) {
451 0 : if (u > 1) {
452 0 : val += (u - 1) * data->R[t][u - 2][v][n + 1];
453 : }
454 0 : val += Y_QC * data->R[t][u - 1][v][n + 1];
455 : } else {
456 0 : if (t > 1) {
457 0 : val += (t - 1) * data->R[t - 2][u][v][n + 1];
458 : }
459 0 : val += X_QC * data->R[t - 1][u][v][n + 1];
460 : }
461 :
462 0 : data->R[t][u][v][n] = val;
463 : }
464 : }
465 : }
466 : }
467 0 : }
468 :
469 : /**
470 : * Calculation of the R_{tuv}^N auxiliary integrals for the 1/r^2 operator
471 : */
472 0 : static void setup_G_array(struct mmd_data *data, int L_A, int L_B) {
473 0 : double q = data->q;
474 0 : double X_QC = data->Q[0] - data->C[0];
475 0 : double Y_QC = data->Q[1] - data->C[1];
476 0 : double Z_QC = data->Q[2] - data->C[2];
477 :
478 0 : int nmax = L_A + L_B;
479 0 : int L_SUM = L_A + L_B;
480 :
481 0 : for (int t = 0; t <= L_SUM; t++) {
482 0 : for (int u = 0; u <= L_SUM; u++) {
483 0 : for (int v = 0; v <= L_SUM; v++) {
484 :
485 0 : if (t + u + v > L_SUM) {
486 0 : continue;
487 : }
488 :
489 0 : for (int n = 0; n <= nmax - (t + u + v); n++) {
490 :
491 0 : double val = 0.0;
492 :
493 0 : if (t + u + v == 0) {
494 0 : val = pow(2.0 * q, n) * data->boys[n];
495 0 : } else if (t + u == 0) {
496 0 : if (v > 1) {
497 0 : val += (v - 1) * (data->R[t][u][v - 2][n + 1] -
498 0 : 2.0 * q * data->R[t][u][v - 2][n]);
499 : }
500 0 : val += Z_QC * (data->R[t][u][v - 1][n + 1] -
501 0 : 2.0 * q * data->R[t][u][v - 1][n]);
502 0 : } else if (t == 0) {
503 0 : if (u > 1) {
504 0 : val += (u - 1) * (data->R[t][u - 2][v][n + 1] -
505 0 : 2.0 * q * data->R[t][u - 2][v][n]);
506 : }
507 0 : val += Y_QC * (data->R[t][u - 1][v][n + 1] -
508 0 : 2.0 * q * data->R[t][u - 1][v][n]);
509 : } else {
510 0 : if (t > 1) {
511 0 : val += (t - 1) * (data->R[t - 2][u][v][n + 1] -
512 0 : 2.0 * q * data->R[t - 2][u][v][n]);
513 : }
514 0 : val += X_QC * (data->R[t - 1][u][v][n + 1] -
515 0 : 2.0 * q * data->R[t - 1][u][v][n]);
516 : }
517 :
518 0 : data->R[t][u][v][n] = val;
519 : }
520 : }
521 : }
522 : }
523 0 : }
524 :
525 : /**
526 : * Calculates E^{ij}_t coefficients in the MMD scheme.
527 : */
528 0 : static void setup_E_array(struct mmd_data *data, int L_A, int L_B) {
529 0 : double q = data->q;
530 :
531 0 : for (int coord = 0; coord < 3; coord++) {
532 0 : for (int i = 0; i <= L_A; i++) {
533 0 : for (int j = 0; j <= L_B; j++) {
534 :
535 0 : for (int t = 0; t <= i + j; t++) {
536 :
537 0 : if (t == 0) {
538 0 : if (i == 0 && j == 0) {
539 0 : data->E[coord][0][0][0] = data->K_abc[coord];
540 0 : continue;
541 0 : } else if (i == 1 && j == 0) {
542 0 : double X_QA = data->Q[coord] - data->A[coord];
543 0 : data->E[coord][1][0][0] = X_QA * data->E[coord][0][0][0];
544 0 : continue;
545 0 : } else if (i == 0 && j == 1) {
546 0 : double X_QB = data->Q[coord] - data->B[coord];
547 0 : data->E[coord][0][1][0] = X_QB * data->E[coord][0][0][0];
548 0 : continue;
549 0 : } else if (i == 0) {
550 0 : double X_QB = data->Q[coord] - data->B[coord];
551 0 : data->E[coord][i][j][0] = X_QB * data->E[coord][i][j - 1][0] +
552 0 : data->E[coord][i][j - 1][1];
553 0 : continue;
554 0 : } else {
555 0 : double X_QA = data->Q[coord] - data->A[coord];
556 0 : data->E[coord][i][j][0] = X_QA * data->E[coord][i - 1][j][0] +
557 0 : data->E[coord][i - 1][j][1];
558 0 : continue;
559 : }
560 : } else {
561 0 : double E_ijt = 0.0;
562 0 : double factor = 1.0 / (2.0 * q * t);
563 :
564 0 : if (i > 0) {
565 0 : E_ijt += factor * i * data->E[coord][i - 1][j][t - 1];
566 : }
567 0 : if (j > 0) {
568 0 : E_ijt += factor * j * data->E[coord][i][j - 1][t - 1];
569 : }
570 :
571 0 : data->E[coord][i][j][t] = E_ijt;
572 : }
573 : }
574 : }
575 : }
576 : }
577 0 : }
|
