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 : #include "libgrpp.h"
16 :
17 : #include <stdio.h>
18 : #include <stdlib.h>
19 : #include <string.h>
20 :
21 : #include "grpp_diff_gaussian.h"
22 : #include "grpp_utils.h"
23 :
24 : void grpp_gradient_diff_bra_contribution(
25 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
26 : libgrpp_grpp_t *grpp_operator, double *grpp_origin, double **grad_arep,
27 : double **grad_so_x, double **grad_so_y, double **grad_so_z, double factor);
28 :
29 : void grpp_gradient_diff_bra_grpp_integrals(
30 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
31 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
32 : double **arep_matrix_down, double **so_x_matrix_down,
33 : double **so_y_matrix_down, double **so_z_matrix_down,
34 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
35 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up);
36 :
37 : void grpp_gradient_diff_ket_contribution(
38 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
39 : libgrpp_grpp_t *grpp_operator, double *grpp_origin, double **grad_arep,
40 : double **grad_so_x, double **grad_so_y, double **grad_so_z, double factor);
41 :
42 : void grpp_gradient_diff_ket_grpp_integrals(
43 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
44 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
45 : double **arep_matrix_down, double **so_x_matrix_down,
46 : double **so_y_matrix_down, double **so_z_matrix_down,
47 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
48 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up);
49 :
50 : void grpp_gradient_contribution(libgrpp_shell_t *shell_A,
51 : libgrpp_shell_t *shell_B,
52 : libgrpp_grpp_t *grpp_operator,
53 : double *grpp_origin, double **grad_arep,
54 : double **grad_so_x, double **grad_so_y,
55 : double **grad_so_z, int diff_bra,
56 : double factor);
57 :
58 : void grpp_gradient_diff_gaussian(
59 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
60 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
61 : double **arep_matrix_down, double **so_x_matrix_down,
62 : double **so_y_matrix_down, double **so_z_matrix_down,
63 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
64 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up,
65 : int diff_bra);
66 :
67 : extern int libgrpp_nlm_to_linear(int *nlm);
68 :
69 : double **libgrpp_alloc_gradients(libgrpp_shell_t *bra, libgrpp_shell_t *ket);
70 :
71 : void libgrpp_dealloc_gradients(double **grad);
72 :
73 : /**
74 : * Analytic calculation of gradients of LOCAL potential integrals for a given
75 : * shell pair with respect to the point 'point_3d'.
76 : */
77 0 : void libgrpp_type1_integrals_gradient(libgrpp_shell_t *shell_A,
78 : libgrpp_shell_t *shell_B,
79 : double *grpp_origin,
80 : libgrpp_potential_t *potential,
81 : double *point_3d, double **grad_arep) {
82 0 : libgrpp_grpp_t *grpp_operator = libgrpp_new_grpp();
83 0 : libgrpp_grpp_set_local_potential(grpp_operator, potential);
84 :
85 : /*
86 : * these arrays are not actually used.
87 : * they are needed only in order to use the
88 : * libgrpp_full_grpp_integrals_gradient() routine.
89 : */
90 0 : double **stub_grad_so_x = libgrpp_alloc_gradients(shell_A, shell_B);
91 0 : double **stub_grad_so_y = libgrpp_alloc_gradients(shell_A, shell_B);
92 0 : double **stub_grad_so_z = libgrpp_alloc_gradients(shell_A, shell_B);
93 :
94 0 : libgrpp_full_grpp_integrals_gradient(
95 : shell_A, shell_B, grpp_operator, grpp_origin, point_3d, grad_arep,
96 : stub_grad_so_x, stub_grad_so_y, stub_grad_so_z);
97 :
98 0 : libgrpp_dealloc_gradients(stub_grad_so_x);
99 0 : libgrpp_dealloc_gradients(stub_grad_so_y);
100 0 : libgrpp_dealloc_gradients(stub_grad_so_z);
101 :
102 0 : grpp_operator->U_L = NULL;
103 0 : libgrpp_delete_grpp(grpp_operator);
104 0 : }
105 :
106 : /**
107 : * Analytic calculation of gradients of SEMI-LOCAL potential integrals for a
108 : * given shell pair with respect to the point 'point_3d'.
109 : */
110 0 : void libgrpp_type2_integrals_gradient(libgrpp_shell_t *shell_A,
111 : libgrpp_shell_t *shell_B,
112 : double *grpp_origin,
113 : libgrpp_potential_t *potential,
114 : double *point_3d, double **grad_arep) {
115 0 : libgrpp_grpp_t *grpp_operator = libgrpp_new_grpp();
116 0 : libgrpp_grpp_add_averaged_potential(grpp_operator, potential);
117 :
118 : /*
119 : * these arrays are not actually used.
120 : * they are needed only in order to use the
121 : * libgrpp_full_grpp_integrals_gradient() routine.
122 : */
123 0 : double **stub_grad_so_x = libgrpp_alloc_gradients(shell_A, shell_B);
124 0 : double **stub_grad_so_y = libgrpp_alloc_gradients(shell_A, shell_B);
125 0 : double **stub_grad_so_z = libgrpp_alloc_gradients(shell_A, shell_B);
126 :
127 0 : libgrpp_full_grpp_integrals_gradient(
128 : shell_A, shell_B, grpp_operator, grpp_origin, point_3d, grad_arep,
129 : stub_grad_so_x, stub_grad_so_y, stub_grad_so_z);
130 :
131 0 : libgrpp_dealloc_gradients(stub_grad_so_x);
132 0 : libgrpp_dealloc_gradients(stub_grad_so_y);
133 0 : libgrpp_dealloc_gradients(stub_grad_so_z);
134 :
135 0 : grpp_operator->n_arep = 0;
136 0 : libgrpp_delete_grpp(grpp_operator);
137 0 : }
138 :
139 : /**
140 : * Analytic calculation of gradients of integrals over the effective spin-orbit
141 : * operator (potential) for a given shell pair (with respect to the point
142 : * 'point_3d').
143 : */
144 0 : void libgrpp_spin_orbit_integrals_gradient(
145 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *grpp_origin,
146 : libgrpp_potential_t *potential, double *point_3d, double **grad_so_x,
147 : double **grad_so_y, double **grad_so_z) {
148 0 : libgrpp_grpp_t *grpp_operator = libgrpp_new_grpp();
149 0 : libgrpp_grpp_add_spin_orbit_potential(grpp_operator, potential);
150 :
151 : /*
152 : * this array is not actually used and is needed only in order
153 : * to use the libgrpp_full_grpp_integrals_gradient() routine.
154 : */
155 0 : double **stub_grad_arep = libgrpp_alloc_gradients(shell_A, shell_B);
156 :
157 0 : libgrpp_full_grpp_integrals_gradient(shell_A, shell_B, grpp_operator,
158 : grpp_origin, point_3d, stub_grad_arep,
159 : grad_so_x, grad_so_y, grad_so_z);
160 :
161 : /*
162 : * inside the libgrpp_full_grpp_integrals_gradient() function
163 : * the SO potential was scaled by 2/(2L+1). Thus the result has to be
164 : * re-scaled by (2L+1)/2 to get rid of any problems with pre-factor
165 : */
166 0 : int L = potential->L;
167 0 : int buf_size = shell_A->cart_size * shell_B->cart_size;
168 0 : for (int icoord = 0; icoord < 3; icoord++) {
169 0 : for (int i = 0; i < buf_size; i++) {
170 0 : grad_so_x[icoord][i] *= (2.0 * L + 1.0) / 2.0;
171 0 : grad_so_y[icoord][i] *= (2.0 * L + 1.0) / 2.0;
172 0 : grad_so_z[icoord][i] *= (2.0 * L + 1.0) / 2.0;
173 : }
174 : }
175 :
176 0 : libgrpp_dealloc_gradients(stub_grad_arep);
177 :
178 0 : grpp_operator->n_esop = 0;
179 0 : libgrpp_delete_grpp(grpp_operator);
180 0 : }
181 :
182 0 : void libgrpp_outercore_potential_integrals_gradient(
183 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *rpp_origin,
184 : int num_oc_shells, libgrpp_potential_t **oc_potentials,
185 : libgrpp_shell_t **oc_shells, double *point_3d, double **grad_arep,
186 : double **grad_so_x, double **grad_so_y, double **grad_so_z) {
187 0 : libgrpp_grpp_t *grpp_operator = libgrpp_new_grpp();
188 0 : for (int ioc = 0; ioc < num_oc_shells; ioc++) {
189 0 : libgrpp_grpp_add_outercore_potential(grpp_operator, oc_potentials[ioc],
190 0 : oc_shells[ioc]);
191 : }
192 :
193 0 : libgrpp_full_grpp_integrals_gradient(shell_A, shell_B, grpp_operator,
194 : rpp_origin, point_3d, grad_arep,
195 : grad_so_x, grad_so_y, grad_so_z);
196 :
197 0 : grpp_operator->n_oc_shells = 0;
198 0 : libgrpp_delete_grpp(grpp_operator);
199 0 : }
200 :
201 : /**
202 : * Analytic calculation of gradients of GRPP integrals for a given shell pair
203 : * with respect to the point 'point_3d'.
204 : * (for the full GRPP operator which includes local, semi-local and non-local
205 : * parts)
206 : */
207 0 : void libgrpp_full_grpp_integrals_gradient(
208 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
209 : libgrpp_grpp_t *grpp_operator, double *grpp_origin, double *point_3d,
210 : double **grad_arep, double **grad_so_x, double **grad_so_y,
211 : double **grad_so_z) {
212 0 : int cart_size_A = shell_A->cart_size;
213 0 : int cart_size_B = shell_B->cart_size;
214 0 : int buf_size = cart_size_A * cart_size_B;
215 :
216 : /*
217 : * initialization: set all gradients to zero
218 : */
219 0 : for (int icoord = 0; icoord < 3; icoord++) {
220 0 : memset(grad_arep[icoord], 0, sizeof(double) * buf_size);
221 0 : memset(grad_so_x[icoord], 0, sizeof(double) * buf_size);
222 0 : memset(grad_so_y[icoord], 0, sizeof(double) * buf_size);
223 0 : memset(grad_so_z[icoord], 0, sizeof(double) * buf_size);
224 : }
225 :
226 : /*
227 : * d<AAA>/d... = 0
228 : */
229 0 : if (points_are_equal(shell_A->origin, grpp_origin) &&
230 0 : points_are_equal(shell_B->origin, grpp_origin)) {
231 : return;
232 : }
233 :
234 : /*
235 : * d<ACB>/dD = 0
236 : */
237 0 : if (!points_are_equal(shell_A->origin, point_3d) &&
238 0 : !points_are_equal(shell_B->origin, point_3d) &&
239 0 : !points_are_equal(grpp_origin, point_3d)) {
240 : return;
241 : }
242 :
243 0 : double *A = shell_A->origin;
244 0 : double *B = shell_B->origin;
245 0 : double *C = grpp_origin;
246 0 : double *D = point_3d;
247 :
248 0 : const int diff_bra = 1;
249 0 : const int diff_ket = 0;
250 :
251 : /*
252 : * Type ACB
253 : */
254 0 : if (!points_are_equal(A, C) && !points_are_equal(C, B) &&
255 0 : !points_are_equal(A, B)) {
256 0 : if (points_are_equal(A, D)) {
257 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
258 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
259 : diff_bra, +1.0);
260 : }
261 0 : if (points_are_equal(B, D)) {
262 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
263 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
264 : diff_ket, +1.0);
265 : }
266 0 : if (points_are_equal(C, D)) {
267 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
268 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
269 : diff_bra, -1.0);
270 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
271 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
272 : diff_ket, -1.0);
273 : }
274 : }
275 :
276 : /*
277 : * Type ACA
278 : */
279 0 : if (points_are_equal(A, B) && !points_are_equal(A, C)) {
280 0 : if (points_are_equal(A, D)) {
281 0 : grpp_gradient_diff_bra_contribution(shell_A, shell_B, grpp_operator,
282 : grpp_origin, grad_arep, grad_so_x,
283 : grad_so_y, grad_so_z, +1.0);
284 0 : grpp_gradient_diff_ket_contribution(shell_A, shell_B, grpp_operator,
285 : grpp_origin, grad_arep, grad_so_x,
286 : grad_so_y, grad_so_z, +1.0);
287 : } else {
288 0 : grpp_gradient_diff_bra_contribution(shell_A, shell_B, grpp_operator,
289 : grpp_origin, grad_arep, grad_so_x,
290 : grad_so_y, grad_so_z, -1.0);
291 0 : grpp_gradient_diff_ket_contribution(shell_A, shell_B, grpp_operator,
292 : grpp_origin, grad_arep, grad_so_x,
293 : grad_so_y, grad_so_z, -1.0);
294 : }
295 : }
296 :
297 : /*
298 : * Type ACC
299 : */
300 0 : if (!points_are_equal(A, C) && points_are_equal(C, B)) {
301 0 : if (points_are_equal(A, D)) {
302 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
303 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
304 : diff_bra, +1.0);
305 : } else {
306 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
307 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
308 : diff_bra, -1.0);
309 : }
310 : }
311 :
312 : /*
313 : * Type CCB
314 : */
315 0 : if (points_are_equal(A, C) && !points_are_equal(C, B)) {
316 0 : if (points_are_equal(B, D)) {
317 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
318 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
319 : diff_ket, +1.0);
320 : } else {
321 0 : grpp_gradient_contribution(shell_A, shell_B, grpp_operator, grpp_origin,
322 : grad_arep, grad_so_x, grad_so_y, grad_so_z,
323 : diff_ket, -1.0);
324 : }
325 : }
326 : }
327 :
328 : /**
329 : * Calculates contribution to gradients arising from the < df/dA | V | g > term:
330 : *
331 : * grad += factor * < df/dA | V | g >
332 : *
333 : * (bra basis function is differentiated).
334 : */
335 0 : void grpp_gradient_diff_bra_contribution(
336 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
337 : libgrpp_grpp_t *grpp_operator, double *grpp_origin, double **grad_arep,
338 : double **grad_so_x, double **grad_so_y, double **grad_so_z, double factor) {
339 : /*
340 : * calculate integrals < df/dA | V | B >
341 : */
342 0 : double *arep_matrix_down = NULL;
343 0 : double *so_x_matrix_down = NULL;
344 0 : double *so_y_matrix_down = NULL;
345 0 : double *so_z_matrix_down = NULL;
346 0 : double *arep_matrix_up = NULL;
347 0 : double *so_x_matrix_up = NULL;
348 0 : double *so_y_matrix_up = NULL;
349 0 : double *so_z_matrix_up = NULL;
350 0 : int cart_size_down = 0;
351 0 : int cart_size_up = 0;
352 :
353 0 : grpp_gradient_diff_bra_grpp_integrals(
354 : shell_A, shell_B, grpp_operator, grpp_origin, &arep_matrix_down,
355 : &so_x_matrix_down, &so_y_matrix_down, &so_z_matrix_down, &arep_matrix_up,
356 : &so_x_matrix_up, &so_y_matrix_up, &so_z_matrix_up, &cart_size_down,
357 : &cart_size_up);
358 :
359 : /*
360 : * construct contributions to gradients:
361 : * d<A|V|B>/dA += < df/dA | V | B >
362 : */
363 0 : for (int icoord = 0; icoord < 3; icoord++) {
364 0 : for (int i = 0; i < shell_A->cart_size; i++) {
365 0 : for (int j = 0; j < shell_B->cart_size; j++) {
366 :
367 0 : int bra_nlm[3];
368 0 : bra_nlm[0] = shell_A->cart_list[3 * i + 0];
369 0 : bra_nlm[1] = shell_A->cart_list[3 * i + 1];
370 0 : bra_nlm[2] = shell_A->cart_list[3 * i + 2];
371 :
372 0 : int ket_nlm[3];
373 0 : ket_nlm[0] = shell_B->cart_list[3 * j + 0];
374 0 : ket_nlm[1] = shell_B->cart_list[3 * j + 1];
375 0 : ket_nlm[2] = shell_B->cart_list[3 * j + 2];
376 :
377 0 : int index = i * shell_B->cart_size + j;
378 :
379 : /*
380 : * contribution from the L-1 gaussian
381 : */
382 0 : if (shell_A->L > 0) {
383 0 : bra_nlm[icoord] -= 1;
384 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
385 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
386 0 : bra_nlm[icoord] += 1;
387 :
388 0 : grad_arep[icoord][index] -=
389 0 : factor * bra_nlm[icoord] *
390 0 : arep_matrix_down[shell_B->cart_size * bra_index + ket_index];
391 0 : grad_so_x[icoord][index] -=
392 0 : factor * bra_nlm[icoord] *
393 0 : so_x_matrix_down[shell_B->cart_size * bra_index + ket_index];
394 0 : grad_so_y[icoord][index] -=
395 0 : factor * bra_nlm[icoord] *
396 0 : so_y_matrix_down[shell_B->cart_size * bra_index + ket_index];
397 0 : grad_so_z[icoord][index] -=
398 0 : factor * bra_nlm[icoord] *
399 0 : so_z_matrix_down[shell_B->cart_size * bra_index + ket_index];
400 : }
401 :
402 : /*
403 : * contribution from the L+1 gaussian
404 : */
405 0 : bra_nlm[icoord] += 1;
406 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
407 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
408 0 : bra_nlm[icoord] -= 1;
409 :
410 0 : grad_arep[icoord][index] +=
411 0 : factor * arep_matrix_up[shell_B->cart_size * bra_index + ket_index];
412 0 : grad_so_x[icoord][index] +=
413 0 : factor * so_x_matrix_up[shell_B->cart_size * bra_index + ket_index];
414 0 : grad_so_y[icoord][index] +=
415 0 : factor * so_y_matrix_up[shell_B->cart_size * bra_index + ket_index];
416 0 : grad_so_z[icoord][index] +=
417 0 : factor * so_z_matrix_up[shell_B->cart_size * bra_index + ket_index];
418 : }
419 : }
420 : }
421 :
422 0 : if (arep_matrix_down) {
423 0 : free(arep_matrix_down);
424 0 : free(so_x_matrix_down);
425 0 : free(so_y_matrix_down);
426 0 : free(so_z_matrix_down);
427 : }
428 0 : free(arep_matrix_up);
429 0 : free(so_x_matrix_up);
430 0 : free(so_y_matrix_up);
431 0 : free(so_z_matrix_up);
432 0 : }
433 :
434 : /**
435 : * To assemble the contribution < df/dA | V | g > to gradients, one have to
436 : * differentiate a Gaussian function. Such a differentiation yields two
437 : * Gaussians with angular momenta L-1 ("down") and L+1 ("up"): dG/dA -> G(L-1)
438 : * and G(L+1)
439 : *
440 : * This function constructs overlap matrices with these "downgraded" and
441 : * "upgraded" Gaussian functions: < G(L-1) | V | G' > and < G(L+1) | V | G' >
442 : *
443 : */
444 0 : void grpp_gradient_diff_bra_grpp_integrals(
445 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
446 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
447 : double **arep_matrix_down, double **so_x_matrix_down,
448 : double **so_y_matrix_down, double **so_z_matrix_down,
449 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
450 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up) {
451 : /*
452 : * differentiation of contracted Gaussian functions
453 : */
454 0 : libgrpp_shell_t *shell_A_down = NULL;
455 0 : libgrpp_shell_t *shell_A_up = NULL;
456 0 : libgrpp_differentiate_shell(shell_A, &shell_A_down, &shell_A_up);
457 :
458 0 : *cart_size_down = 0;
459 0 : if (shell_A_down != NULL) {
460 0 : *cart_size_down = shell_A_down->cart_size;
461 : }
462 0 : *cart_size_up = shell_A_up->cart_size;
463 :
464 : /*
465 : * matrix < L-1 | V | L >
466 : */
467 0 : if (shell_A_down != NULL) {
468 0 : size_t mat_size_down = shell_A_down->cart_size * shell_B->cart_size;
469 0 : *arep_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
470 0 : *so_x_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
471 0 : *so_y_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
472 0 : *so_z_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
473 :
474 0 : libgrpp_full_grpp_integrals(
475 : shell_A_down, shell_B, grpp_operator, grpp_origin, *arep_matrix_down,
476 : *so_x_matrix_down, *so_y_matrix_down, *so_z_matrix_down);
477 : } else {
478 0 : *arep_matrix_down = NULL;
479 0 : *so_x_matrix_down = NULL;
480 0 : *so_y_matrix_down = NULL;
481 0 : *so_z_matrix_down = NULL;
482 : }
483 :
484 : /*
485 : * matrix < L+1 | V | L >
486 : */
487 0 : size_t mat_size_up = shell_A_up->cart_size * shell_B->cart_size;
488 0 : *arep_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
489 0 : *so_x_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
490 0 : *so_y_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
491 0 : *so_z_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
492 :
493 0 : libgrpp_full_grpp_integrals(shell_A_up, shell_B, grpp_operator, grpp_origin,
494 : *arep_matrix_up, *so_x_matrix_up, *so_y_matrix_up,
495 : *so_z_matrix_up);
496 :
497 : /*
498 : * clean up
499 : */
500 0 : if (shell_A_down) {
501 0 : libgrpp_delete_shell(shell_A_down);
502 : }
503 0 : libgrpp_delete_shell(shell_A_up);
504 0 : }
505 :
506 : /**
507 : * Calculates contribution to gradients arising from the < df/dA | V | g > term:
508 : *
509 : * grad += factor * < f | V | dg/dA >
510 : *
511 : * (bra basis function is differentiated).
512 : */
513 0 : void grpp_gradient_diff_ket_contribution(
514 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
515 : libgrpp_grpp_t *grpp_operator, double *grpp_origin, double **grad_arep,
516 : double **grad_so_x, double **grad_so_y, double **grad_so_z, double factor) {
517 : /*
518 : * calculate integrals < df/dA | V | B >
519 : */
520 0 : double *arep_matrix_down = NULL;
521 0 : double *so_x_matrix_down = NULL;
522 0 : double *so_y_matrix_down = NULL;
523 0 : double *so_z_matrix_down = NULL;
524 0 : double *arep_matrix_up = NULL;
525 0 : double *so_x_matrix_up = NULL;
526 0 : double *so_y_matrix_up = NULL;
527 0 : double *so_z_matrix_up = NULL;
528 0 : int cart_size_down = 0;
529 0 : int cart_size_up = 0;
530 :
531 0 : grpp_gradient_diff_ket_grpp_integrals(
532 : shell_A, shell_B, grpp_operator, grpp_origin, &arep_matrix_down,
533 : &so_x_matrix_down, &so_y_matrix_down, &so_z_matrix_down, &arep_matrix_up,
534 : &so_x_matrix_up, &so_y_matrix_up, &so_z_matrix_up, &cart_size_down,
535 : &cart_size_up);
536 :
537 : /*
538 : * construct contributions to gradients:
539 : * d<A|B>/dA += < df/dA | V | B >
540 : */
541 0 : for (int icoord = 0; icoord < 3; icoord++) {
542 0 : for (int i = 0; i < shell_A->cart_size; i++) {
543 0 : for (int j = 0; j < shell_B->cart_size; j++) {
544 :
545 0 : int bra_nlm[3];
546 0 : bra_nlm[0] = shell_A->cart_list[3 * i + 0];
547 0 : bra_nlm[1] = shell_A->cart_list[3 * i + 1];
548 0 : bra_nlm[2] = shell_A->cart_list[3 * i + 2];
549 :
550 0 : int ket_nlm[3];
551 0 : ket_nlm[0] = shell_B->cart_list[3 * j + 0];
552 0 : ket_nlm[1] = shell_B->cart_list[3 * j + 1];
553 0 : ket_nlm[2] = shell_B->cart_list[3 * j + 2];
554 :
555 0 : int index = i * shell_B->cart_size + j;
556 :
557 : /*
558 : * contribution from the L-1 gaussian
559 : */
560 0 : if (shell_B->L > 0) {
561 0 : ket_nlm[icoord] -= 1;
562 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
563 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
564 0 : ket_nlm[icoord] += 1;
565 :
566 0 : grad_arep[icoord][index] -=
567 0 : factor * ket_nlm[icoord] *
568 0 : arep_matrix_down[cart_size_down * bra_index + ket_index];
569 0 : grad_so_x[icoord][index] -=
570 0 : factor * ket_nlm[icoord] *
571 0 : so_x_matrix_down[cart_size_down * bra_index + ket_index];
572 0 : grad_so_y[icoord][index] -=
573 0 : factor * ket_nlm[icoord] *
574 0 : so_y_matrix_down[cart_size_down * bra_index + ket_index];
575 0 : grad_so_z[icoord][index] -=
576 0 : factor * ket_nlm[icoord] *
577 0 : so_z_matrix_down[cart_size_down * bra_index + ket_index];
578 : }
579 :
580 : /*
581 : * contribution from the L+1 gaussian
582 : */
583 0 : ket_nlm[icoord] += 1;
584 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
585 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
586 0 : ket_nlm[icoord] -= 1;
587 :
588 0 : grad_arep[icoord][index] +=
589 0 : factor * arep_matrix_up[cart_size_up * bra_index + ket_index];
590 0 : grad_so_x[icoord][index] +=
591 0 : factor * so_x_matrix_up[cart_size_up * bra_index + ket_index];
592 0 : grad_so_y[icoord][index] +=
593 0 : factor * so_y_matrix_up[cart_size_up * bra_index + ket_index];
594 0 : grad_so_z[icoord][index] +=
595 0 : factor * so_z_matrix_up[cart_size_up * bra_index + ket_index];
596 : }
597 : }
598 : }
599 :
600 0 : if (arep_matrix_down) {
601 0 : free(arep_matrix_down);
602 0 : free(so_x_matrix_down);
603 0 : free(so_y_matrix_down);
604 0 : free(so_z_matrix_down);
605 : }
606 0 : free(arep_matrix_up);
607 0 : free(so_x_matrix_up);
608 0 : free(so_y_matrix_up);
609 0 : free(so_z_matrix_up);
610 0 : }
611 :
612 : /**
613 : * To assemble the contribution < df/dA | V | g > to gradients, one have to
614 : * differentiate Gaussian function. Such a differentiation yields two Gaussians
615 : * with angular momenta L-1 ("down") and L+1 ("up"): dG/dA -> G(L-1) and G(L+1)
616 : *
617 : * This function constructs matrices with these "downgraded" and "upgraded"
618 : * Gaussian functions:
619 : * < G(L-1) | V | G' > and < G(L+1) | V | G' >
620 : *
621 : */
622 0 : void grpp_gradient_diff_ket_grpp_integrals(
623 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
624 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
625 : double **arep_matrix_down, double **so_x_matrix_down,
626 : double **so_y_matrix_down, double **so_z_matrix_down,
627 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
628 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up) {
629 : /*
630 : * differentiation of contracted Gaussian functions
631 : */
632 0 : libgrpp_shell_t *shell_B_down = NULL;
633 0 : libgrpp_shell_t *shell_B_up = NULL;
634 0 : libgrpp_differentiate_shell(shell_B, &shell_B_down, &shell_B_up);
635 :
636 0 : *cart_size_down = 0;
637 0 : if (shell_B_down != NULL) {
638 0 : *cart_size_down = shell_B_down->cart_size;
639 : }
640 0 : *cart_size_up = shell_B_up->cart_size;
641 :
642 : /*
643 : * matrix < L-1 | L>
644 : */
645 0 : if (shell_B_down != NULL) {
646 0 : size_t mat_size_down = shell_A->cart_size * shell_B_down->cart_size;
647 0 : *arep_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
648 0 : *so_x_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
649 0 : *so_y_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
650 0 : *so_z_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
651 :
652 0 : libgrpp_full_grpp_integrals(
653 : // evaluate_grpp_integrals_shell_pair(
654 : shell_A, shell_B_down, grpp_operator, grpp_origin, *arep_matrix_down,
655 : *so_x_matrix_down, *so_y_matrix_down, *so_z_matrix_down);
656 : } else {
657 0 : *arep_matrix_down = NULL;
658 0 : *so_x_matrix_down = NULL;
659 0 : *so_y_matrix_down = NULL;
660 0 : *so_z_matrix_down = NULL;
661 : }
662 :
663 : /*
664 : * matrix < L+1 | L>
665 : */
666 0 : size_t mat_size_up = shell_A->cart_size * shell_B_up->cart_size;
667 0 : *arep_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
668 0 : *so_x_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
669 0 : *so_y_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
670 0 : *so_z_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
671 :
672 0 : libgrpp_full_grpp_integrals(
673 : // evaluate_grpp_integrals_shell_pair(
674 : shell_A, shell_B_up, grpp_operator, grpp_origin, *arep_matrix_up,
675 : *so_x_matrix_up, *so_y_matrix_up, *so_z_matrix_up);
676 :
677 : /*
678 : * clean up
679 : */
680 0 : if (shell_B_down) {
681 0 : libgrpp_delete_shell(shell_B_down);
682 : }
683 0 : libgrpp_delete_shell(shell_B_up);
684 0 : }
685 :
686 0 : void grpp_gradient_contribution(libgrpp_shell_t *shell_A,
687 : libgrpp_shell_t *shell_B,
688 : libgrpp_grpp_t *grpp_operator,
689 : double *grpp_origin, double **grad_arep,
690 : double **grad_so_x, double **grad_so_y,
691 : double **grad_so_z, int diff_bra,
692 : double factor) {
693 : // int diff_ket = 0;
694 :
695 0 : if (diff_bra == 0) {
696 : diff_bra = 0;
697 : // diff_ket = 1;
698 : } else {
699 0 : diff_bra = 1;
700 : // diff_ket = 0;
701 : }
702 :
703 : /*
704 : * calculate overlap integrals < df/dA | V | B >
705 : */
706 0 : double *arep_matrix_down = NULL;
707 0 : double *so_x_matrix_down = NULL;
708 0 : double *so_y_matrix_down = NULL;
709 0 : double *so_z_matrix_down = NULL;
710 0 : double *arep_matrix_up = NULL;
711 0 : double *so_x_matrix_up = NULL;
712 0 : double *so_y_matrix_up = NULL;
713 0 : double *so_z_matrix_up = NULL;
714 0 : int cart_size_down = 0;
715 0 : int cart_size_up = 0;
716 :
717 0 : grpp_gradient_diff_gaussian(
718 : shell_A, shell_B, grpp_operator, grpp_origin, &arep_matrix_down,
719 : &so_x_matrix_down, &so_y_matrix_down, &so_z_matrix_down, &arep_matrix_up,
720 : &so_x_matrix_up, &so_y_matrix_up, &so_z_matrix_up, &cart_size_down,
721 : &cart_size_up, diff_bra);
722 :
723 : /*
724 : * construct contributions to gradients:
725 : * d<A|U|B>/dA += < df/dA | U | B >
726 : */
727 0 : for (int icoord = 0; icoord < 3; icoord++) {
728 0 : for (int i = 0; i < shell_A->cart_size; i++) {
729 0 : for (int j = 0; j < shell_B->cart_size; j++) {
730 :
731 0 : int bra_nlm[3];
732 0 : bra_nlm[0] = shell_A->cart_list[3 * i + 0];
733 0 : bra_nlm[1] = shell_A->cart_list[3 * i + 1];
734 0 : bra_nlm[2] = shell_A->cart_list[3 * i + 2];
735 :
736 0 : int ket_nlm[3];
737 0 : ket_nlm[0] = shell_B->cart_list[3 * j + 0];
738 0 : ket_nlm[1] = shell_B->cart_list[3 * j + 1];
739 0 : ket_nlm[2] = shell_B->cart_list[3 * j + 2];
740 :
741 0 : int index = i * shell_B->cart_size + j;
742 :
743 0 : int *diff_nlm = diff_bra ? bra_nlm : ket_nlm;
744 :
745 : /*
746 : * contribution from the L-1 gaussian
747 : */
748 0 : if (cart_size_down > 0) {
749 0 : diff_nlm[icoord] -= 1;
750 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
751 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
752 0 : diff_nlm[icoord] += 1;
753 :
754 0 : int n = diff_nlm[icoord];
755 0 : int row_len = diff_bra ? shell_B->cart_size : cart_size_down;
756 0 : int index_down = row_len * bra_index + ket_index;
757 :
758 0 : grad_arep[icoord][index] -= factor * n * arep_matrix_down[index_down];
759 0 : grad_so_x[icoord][index] -= factor * n * so_x_matrix_down[index_down];
760 0 : grad_so_y[icoord][index] -= factor * n * so_y_matrix_down[index_down];
761 0 : grad_so_z[icoord][index] -= factor * n * so_z_matrix_down[index_down];
762 : }
763 :
764 : /*
765 : * contribution from the L+1 gaussian
766 : */
767 0 : diff_nlm[icoord] += 1;
768 0 : int bra_index = libgrpp_nlm_to_linear(bra_nlm);
769 0 : int ket_index = libgrpp_nlm_to_linear(ket_nlm);
770 0 : diff_nlm[icoord] -= 1;
771 :
772 0 : int row_len = diff_bra ? shell_B->cart_size : cart_size_up;
773 0 : int index_up = row_len * bra_index + ket_index;
774 :
775 0 : grad_arep[icoord][index] += factor * arep_matrix_up[index_up];
776 0 : grad_so_x[icoord][index] += factor * so_x_matrix_up[index_up];
777 0 : grad_so_y[icoord][index] += factor * so_y_matrix_up[index_up];
778 0 : grad_so_z[icoord][index] += factor * so_z_matrix_up[index_up];
779 : }
780 : }
781 : }
782 :
783 0 : if (arep_matrix_down) {
784 0 : free(arep_matrix_down);
785 0 : free(so_x_matrix_down);
786 0 : free(so_y_matrix_down);
787 0 : free(so_z_matrix_down);
788 : }
789 0 : free(arep_matrix_up);
790 0 : free(so_x_matrix_up);
791 0 : free(so_y_matrix_up);
792 0 : free(so_z_matrix_up);
793 0 : }
794 :
795 : /**
796 : * To assemble the contribution < df/dA | V | g > to gradients, one have to
797 : * differentiate Gaussian function. Such a differentiation yields two Gaussians
798 : * with angular momenta L-1 ("down") and L+1 ("up"): dG/dA -> G(L-1) and G(L+1)
799 : *
800 : * This function constructs matrices with these "downgraded" and "upgraded"
801 : * Gaussian functions:
802 : * < G(L-1) | V | G' > and < G(L+1) | V | G' >
803 : *
804 : */
805 0 : void grpp_gradient_diff_gaussian(
806 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
807 : libgrpp_grpp_t *grpp_operator, double *grpp_origin,
808 : double **arep_matrix_down, double **so_x_matrix_down,
809 : double **so_y_matrix_down, double **so_z_matrix_down,
810 : double **arep_matrix_up, double **so_x_matrix_up, double **so_y_matrix_up,
811 : double **so_z_matrix_up, int *cart_size_down, int *cart_size_up,
812 : int diff_bra) {
813 : // int diff_ket = 0;
814 :
815 0 : if (diff_bra == 0) {
816 : diff_bra = 0;
817 : // diff_ket = 1;
818 : } else {
819 0 : diff_bra = 1;
820 : // diff_ket = 0;
821 : }
822 :
823 : /*
824 : * which shell should be differentiated, bra or ket
825 : */
826 0 : libgrpp_shell_t *const_shell = NULL;
827 0 : libgrpp_shell_t *diff_shell = NULL;
828 0 : if (diff_bra) {
829 : diff_shell = shell_A;
830 : const_shell = shell_B;
831 : } else {
832 : diff_shell = shell_B;
833 : const_shell = shell_A;
834 : }
835 :
836 : /*
837 : * differentiation of contracted Gaussian functions
838 : */
839 0 : libgrpp_shell_t *shell_down = NULL;
840 0 : libgrpp_shell_t *shell_up = NULL;
841 0 : libgrpp_differentiate_shell(diff_shell, &shell_down, &shell_up);
842 :
843 0 : *cart_size_down = 0;
844 0 : if (shell_down != NULL) {
845 0 : *cart_size_down = shell_down->cart_size;
846 : }
847 0 : *cart_size_up = shell_up->cart_size;
848 :
849 : /*
850 : * GRPP matrix:
851 : * < L-1 | U | L > or < L | U | L-1 >
852 : */
853 0 : if (shell_down != NULL) {
854 0 : size_t mat_size_down = const_shell->cart_size * shell_down->cart_size;
855 0 : *arep_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
856 0 : *so_x_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
857 0 : *so_y_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
858 0 : *so_z_matrix_down = (double *)calloc(mat_size_down, sizeof(double));
859 :
860 0 : libgrpp_full_grpp_integrals(
861 : diff_bra ? shell_down : shell_A, diff_bra ? shell_B : shell_down,
862 : grpp_operator, grpp_origin, *arep_matrix_down, *so_x_matrix_down,
863 : *so_y_matrix_down, *so_z_matrix_down);
864 : } else {
865 0 : *arep_matrix_down = NULL;
866 0 : *so_x_matrix_down = NULL;
867 0 : *so_y_matrix_down = NULL;
868 0 : *so_z_matrix_down = NULL;
869 : }
870 :
871 : /*
872 : * GRPP matrix:
873 : * < L+1 | U | L > or < L | U | L+1 >
874 : */
875 0 : size_t mat_size_up = const_shell->cart_size * shell_up->cart_size;
876 0 : *arep_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
877 0 : *so_x_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
878 0 : *so_y_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
879 0 : *so_z_matrix_up = (double *)calloc(mat_size_up, sizeof(double));
880 :
881 0 : libgrpp_full_grpp_integrals(diff_bra ? shell_up : shell_A,
882 : diff_bra ? shell_B : shell_up, grpp_operator,
883 : grpp_origin, *arep_matrix_up, *so_x_matrix_up,
884 : *so_y_matrix_up, *so_z_matrix_up);
885 :
886 : /*
887 : * clean up
888 : */
889 0 : if (shell_down) {
890 0 : libgrpp_delete_shell(shell_down);
891 : }
892 0 : libgrpp_delete_shell(shell_up);
893 0 : }
894 :
895 : /**
896 : * Allocates memory for gradients for a given shell pair.
897 : */
898 0 : double **libgrpp_alloc_gradients(libgrpp_shell_t *bra, libgrpp_shell_t *ket) {
899 0 : size_t size = bra->cart_size * ket->cart_size;
900 :
901 0 : double **grad = (double **)calloc(3, sizeof(double *));
902 0 : for (int icoord = 0; icoord < 3; icoord++) {
903 0 : grad[icoord] = (double *)calloc(size, sizeof(double));
904 : }
905 :
906 0 : return grad;
907 : }
908 :
909 : /**
910 : * Deallocates arrays containing gradients of AO integrals.
911 : */
912 0 : void libgrpp_dealloc_gradients(double **grad) {
913 0 : free(grad[0]);
914 0 : free(grad[1]);
915 0 : free(grad[2]);
916 0 : free(grad);
917 0 : }
|
