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 : * Wrappers for the LIBGRPP subroutines to be used from Fortran projects.
17 : *
18 : * C99 Fortran
19 : * --------------------
20 : * int32_t integer(4)
21 : * double real(8)
22 : */
23 :
24 : #include <math.h>
25 : #include <stdint.h>
26 : #include <stdio.h>
27 : #include <stdlib.h>
28 :
29 : #include "grpp_factorial.h"
30 : #include "libgrpp.h"
31 :
32 : /*
33 : * Fine-tuning of the LIBGRPP internal parameters.
34 : */
35 :
36 0 : void libgrpp_set_default_parameters_() { libgrpp_set_default_parameters(); }
37 :
38 0 : void libgrpp_set_radial_tolerance_(const double *tolerance) {
39 0 : libgrpp_set_radial_tolerance(*tolerance);
40 0 : }
41 :
42 0 : void libgrpp_set_angular_screening_tolerance_(const double *tolerance) {
43 0 : libgrpp_set_angular_screening_tolerance(*tolerance);
44 0 : }
45 :
46 0 : void libgrpp_set_modified_bessel_tolerance_(const double *tolerance) {
47 0 : libgrpp_set_modified_bessel_tolerance(*tolerance);
48 0 : }
49 :
50 0 : void libgrpp_set_cartesian_order_(const int32_t *order) {
51 0 : libgrpp_set_cartesian_order(*order);
52 0 : }
53 :
54 : /*
55 : * initialization and finalization
56 : */
57 12600 : void libgrpp_init_() { libgrpp_init(); }
58 :
59 0 : void libgrpp_finalize_() { libgrpp_finalize(); }
60 :
61 : /*
62 : * Type 1 RPP integrals (local term)
63 : */
64 :
65 0 : void libgrpp_type1_integrals_(
66 : // contracted Gaussian A
67 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
68 : double *alpha_A,
69 : // contracted Gaussian B
70 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
71 : double *alpha_B,
72 : // pseudopotential
73 : double *rpp_origin, int32_t *rpp_num_primitives, int32_t *rpp_powers,
74 : double *rpp_coeffs, double *rpp_alpha,
75 : // answer
76 : double *matrix) {
77 0 : int *pot_powers_int = (int *)calloc(*rpp_num_primitives, sizeof(int));
78 :
79 0 : for (int i = 0; i < *rpp_num_primitives; i++) {
80 0 : pot_powers_int[i] = rpp_powers[i];
81 : }
82 :
83 0 : libgrpp_potential_t *pot = libgrpp_new_potential(
84 : 0, 0, *rpp_num_primitives, pot_powers_int, rpp_coeffs, rpp_alpha);
85 0 : libgrpp_shell_t *shell_A =
86 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
87 0 : libgrpp_shell_t *shell_B =
88 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
89 :
90 0 : libgrpp_type1_integrals(shell_A, shell_B, rpp_origin, pot, matrix);
91 :
92 0 : libgrpp_delete_shell(shell_A);
93 0 : libgrpp_delete_shell(shell_B);
94 0 : libgrpp_delete_potential(pot);
95 0 : free(pot_powers_int);
96 0 : }
97 :
98 : /*
99 : * Type 2 RPP integrals (semilocal terms with projectors)
100 : */
101 :
102 194176 : void libgrpp_type2_integrals_(
103 : // contracted Gaussian A
104 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
105 : double *alpha_A, double *origin_B,
106 : // contracted Gaussian B
107 : int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B, double *alpha_B,
108 : // pseudopotential
109 : double *pot_origin, int32_t *pot_L, int32_t *pot_num_primitives,
110 : int32_t *pot_powers, double *pot_coeffs, double *pot_alpha,
111 : // answer
112 : double *matrix) {
113 194176 : int *pot_powers_int = (int *)calloc(*pot_num_primitives, sizeof(int));
114 :
115 719265 : for (int i = 0; i < *pot_num_primitives; i++) {
116 525089 : pot_powers_int[i] = pot_powers[i];
117 : }
118 :
119 194176 : libgrpp_potential_t *pot = libgrpp_new_potential(
120 : *pot_L, 0, *pot_num_primitives, pot_powers_int, pot_coeffs, pot_alpha);
121 194176 : libgrpp_shell_t *shell_A =
122 194176 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
123 194176 : libgrpp_shell_t *shell_B =
124 194176 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
125 :
126 194176 : libgrpp_type2_integrals(shell_A, shell_B, pot_origin, pot, matrix);
127 :
128 194176 : libgrpp_delete_shell(shell_A);
129 194176 : libgrpp_delete_shell(shell_B);
130 194176 : libgrpp_delete_potential(pot);
131 194176 : free(pot_powers_int);
132 194176 : }
133 :
134 : /*
135 : * Effective spin-orbit operator ("Type 3") RPP integrals (semilocal terms with
136 : * projectors)
137 : */
138 :
139 0 : void libgrpp_spin_orbit_integrals_(
140 : // contracted Gaussian A
141 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
142 : double *alpha_A,
143 : // contracted Gaussian B
144 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
145 : double *alpha_B,
146 : // pseudopotential
147 : double *pot_origin, int32_t *pot_L, int32_t *pot_num_primitives,
148 : int32_t *pot_powers, double *pot_coeffs, double *pot_alpha,
149 : // answer
150 : double *so_x_matrix, double *so_y_matrix, double *so_z_matrix) {
151 0 : int *pot_powers_int = (int *)calloc(*pot_num_primitives, sizeof(int));
152 :
153 0 : for (int i = 0; i < *pot_num_primitives; i++) {
154 0 : pot_powers_int[i] = pot_powers[i];
155 : }
156 :
157 : /*
158 : * construct RPP structure
159 : */
160 0 : libgrpp_potential_t *pot = libgrpp_new_potential(
161 : *pot_L, 0, *pot_num_primitives, pot_powers_int, pot_coeffs, pot_alpha);
162 :
163 : /*
164 : * construct shells
165 : */
166 0 : libgrpp_shell_t *shell_A =
167 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
168 0 : libgrpp_shell_t *shell_B =
169 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
170 :
171 0 : libgrpp_spin_orbit_integrals(shell_A, shell_B, pot_origin, pot, so_x_matrix,
172 : so_y_matrix, so_z_matrix);
173 :
174 0 : libgrpp_delete_shell(shell_A);
175 0 : libgrpp_delete_shell(shell_B);
176 0 : libgrpp_delete_potential(pot);
177 0 : free(pot_powers_int);
178 0 : }
179 :
180 : /**
181 : * Outercore RPP integrals (non-local terms with projectors onto outercore
182 : * spinors)
183 : *
184 : * Part 1: integration of the first non-local term:
185 : * U*|nlj><nlj| + |nlj><nlj|*U
186 : */
187 0 : void libgrpp_outercore_potential_integrals_part_1_(
188 : // contracted Gaussian A
189 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
190 : double *alpha_A, double *origin_B,
191 : // contracted Gaussian B
192 : int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B, double *alpha_B,
193 : // pseudopotential for the outercore shell LJ
194 : double *pot_origin, int32_t *pot_L, int32_t *pot_J,
195 : int32_t *pot_num_primitives, int32_t *pot_powers, double *pot_coeffs,
196 : double *pot_alpha,
197 : // expansion of the outercore shell LJ
198 : int32_t *oc_shell_num_primitives, double *oc_shell_coeffs,
199 : double *oc_shell_alpha,
200 : // answer
201 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
202 : double *so_z_matrix) {
203 0 : void libgrpp_outercore_potential_integrals_part_1(
204 : libgrpp_shell_t * shell_A, libgrpp_shell_t * shell_B, double *C,
205 : libgrpp_potential_t *oc_potential, libgrpp_shell_t *oc_shell,
206 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
207 : double *so_z_matrix);
208 :
209 : /*
210 : * array conversion: Fortran -> C
211 : */
212 0 : int *pot_powers_int = (int *)calloc(*pot_num_primitives, sizeof(int));
213 0 : for (int i = 0; i < *pot_num_primitives; i++) {
214 0 : pot_powers_int[i] = pot_powers[i];
215 : }
216 :
217 : /*
218 : * construct shells
219 : */
220 0 : libgrpp_shell_t *shell_A =
221 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
222 0 : libgrpp_shell_t *shell_B =
223 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
224 :
225 : /*
226 : * construct pseudopotential for the given L,J numbers
227 : */
228 0 : libgrpp_potential_t *oc_potential =
229 0 : libgrpp_new_potential(*pot_L, *pot_J, *pot_num_primitives, pot_powers_int,
230 : pot_coeffs, pot_alpha);
231 :
232 : /*
233 : * construct outercore shell associated with the pseudopotential
234 : */
235 0 : libgrpp_shell_t *oc_shell =
236 0 : libgrpp_new_shell(pot_origin, *pot_L, *oc_shell_num_primitives,
237 : oc_shell_coeffs, oc_shell_alpha);
238 :
239 : /*
240 : * evaluate RPP integrals
241 : */
242 0 : libgrpp_outercore_potential_integrals_part_1(
243 : shell_A, shell_B, pot_origin, oc_potential, oc_shell, arep_matrix,
244 : so_x_matrix, so_y_matrix, so_z_matrix);
245 :
246 : /*
247 : * clean-up
248 : */
249 0 : libgrpp_delete_shell(shell_A);
250 0 : libgrpp_delete_shell(shell_B);
251 0 : libgrpp_delete_potential(oc_potential);
252 0 : libgrpp_delete_shell(oc_shell);
253 0 : free(pot_powers_int);
254 0 : }
255 :
256 : /**
257 : * Outercore RPP integrals (non-local terms with projectors onto outercore
258 : * spinors)
259 : *
260 : * Part 2: integration of the second non-local term:
261 : * |nlj><nlj| U |n'lj><n'lj|
262 : */
263 0 : void libgrpp_outercore_potential_integrals_part_2_(
264 : // contracted Gaussian A
265 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
266 : double *alpha_A,
267 : // contracted Gaussian B
268 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
269 : double *alpha_B,
270 : // origin of the RPP
271 : double *pot_origin,
272 : // outercore shell 1:
273 : int32_t *oc_shell_1_L, int32_t *oc_shell_1_J, int32_t *pot1_num_primitives,
274 : int32_t *pot1_powers, double *pot1_coeffs, double *pot1_alpha,
275 : int32_t *oc_shell_1_num_primitives, double *oc_shell_1_coeffs,
276 : double *oc_shell_1_alpha,
277 : // outercore shell 2:
278 : int32_t *oc_shell_2_L, int32_t *oc_shell_2_J, int32_t *pot2_num_primitives,
279 : int32_t *pot2_powers, double *pot2_coeffs, double *pot2_alpha,
280 : int32_t *oc_shell_2_num_primitives, double *oc_shell_2_coeffs,
281 : double *oc_shell_2_alpha,
282 : // answer:
283 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
284 : double *so_z_matrix) {
285 0 : void libgrpp_outercore_potential_integrals_part_2(
286 : libgrpp_shell_t * shell_A, libgrpp_shell_t * shell_B, double *C,
287 : libgrpp_potential_t *oc_potential_1, libgrpp_shell_t *oc_shell_1,
288 : libgrpp_potential_t *oc_potential_2, libgrpp_shell_t *oc_shell_2,
289 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
290 : double *so_z_matrix);
291 :
292 : /*
293 : * array conversion: Fortran -> C
294 : */
295 0 : int *pot1_powers_int = (int *)calloc(*pot1_num_primitives, sizeof(int));
296 0 : int *pot2_powers_int = (int *)calloc(*pot2_num_primitives, sizeof(int));
297 :
298 0 : for (int i = 0; i < *pot1_num_primitives; i++) {
299 0 : pot1_powers_int[i] = pot1_powers[i];
300 : }
301 0 : for (int i = 0; i < *pot2_num_primitives; i++) {
302 0 : pot2_powers_int[i] = pot2_powers[i];
303 : }
304 :
305 : /*
306 : * construct shells
307 : */
308 0 : libgrpp_shell_t *shell_A =
309 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
310 0 : libgrpp_shell_t *shell_B =
311 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
312 :
313 : /*
314 : * the first outercore pseudopotential U_{n,L,J} and the corresponding
315 : * outercore shell
316 : */
317 0 : libgrpp_potential_t *oc_potential_1 =
318 0 : libgrpp_new_potential(*oc_shell_1_L, *oc_shell_1_J, *pot1_num_primitives,
319 : pot1_powers_int, pot1_coeffs, pot1_alpha);
320 0 : libgrpp_shell_t *oc_shell_1 =
321 0 : libgrpp_new_shell(pot_origin, *oc_shell_1_L, *oc_shell_1_num_primitives,
322 : oc_shell_1_coeffs, oc_shell_1_alpha);
323 :
324 : /*
325 : * the second outercore pseudopotential U_{n',L',J'} and the corresponding
326 : * outercore shell
327 : */
328 0 : libgrpp_potential_t *oc_potential_2 =
329 0 : libgrpp_new_potential(*oc_shell_2_L, *oc_shell_2_J, *pot2_num_primitives,
330 : pot2_powers_int, pot2_coeffs, pot2_alpha);
331 0 : libgrpp_shell_t *oc_shell_2 =
332 0 : libgrpp_new_shell(pot_origin, *oc_shell_2_L, *oc_shell_2_num_primitives,
333 : oc_shell_2_coeffs, oc_shell_2_alpha);
334 :
335 : /*
336 : * evaluate integrals
337 : */
338 0 : libgrpp_outercore_potential_integrals_part_2(
339 : shell_A, shell_B, pot_origin, oc_potential_1, oc_shell_1, oc_potential_2,
340 : oc_shell_2, arep_matrix, so_x_matrix, so_y_matrix, so_z_matrix);
341 :
342 : /*
343 : * clean-up
344 : */
345 0 : libgrpp_delete_shell(shell_A);
346 0 : libgrpp_delete_shell(shell_B);
347 0 : libgrpp_delete_potential(oc_potential_1);
348 0 : libgrpp_delete_potential(oc_potential_2);
349 0 : libgrpp_delete_shell(oc_shell_1);
350 0 : libgrpp_delete_shell(oc_shell_2);
351 0 : free(pot1_powers_int);
352 0 : free(pot2_powers_int);
353 0 : }
354 :
355 : /**
356 : * Analytic calculation of gradients of LOCAL potential integrals for a given
357 : * shell pair with respect to the point 'point_3d'.
358 : */
359 0 : void libgrpp_type1_integrals_gradient_(
360 : // contracted Gaussian A
361 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
362 : double *alpha_A,
363 : // contracted Gaussian B
364 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
365 : double *alpha_B,
366 : // pseudopotential
367 : double *rpp_origin, int32_t *rpp_num_primitives, int32_t *rpp_powers,
368 : double *rpp_coeffs, double *rpp_alpha,
369 : // differentiation wrt the 3d point (x,y,z)
370 : double *point_3d,
371 : // answer: matrices d<Int>/dx, d<Int>/dy, d<Int>/dZ
372 : double *grad_arep_x, double *grad_arep_y, double *grad_arep_z) {
373 0 : int *pot_powers_int = (int *)calloc(*rpp_num_primitives, sizeof(int));
374 0 : double *grad_array[3];
375 0 : grad_array[0] = grad_arep_x;
376 0 : grad_array[1] = grad_arep_y;
377 0 : grad_array[2] = grad_arep_z;
378 :
379 0 : for (int i = 0; i < *rpp_num_primitives; i++) {
380 0 : pot_powers_int[i] = rpp_powers[i];
381 : }
382 :
383 0 : libgrpp_potential_t *pot = libgrpp_new_potential(
384 : 0, 0, *rpp_num_primitives, pot_powers_int, rpp_coeffs, rpp_alpha);
385 0 : libgrpp_shell_t *shell_A =
386 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
387 0 : libgrpp_shell_t *shell_B =
388 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
389 :
390 0 : libgrpp_type1_integrals_gradient(shell_A, shell_B, rpp_origin, pot, point_3d,
391 : grad_array);
392 :
393 0 : libgrpp_delete_shell(shell_A);
394 0 : libgrpp_delete_shell(shell_B);
395 0 : libgrpp_delete_potential(pot);
396 0 : free(pot_powers_int);
397 0 : }
398 :
399 : /**
400 : * Analytic calculation of gradients of SEMI-LOCAL potential integrals for a
401 : * given shell pair with respect to the point 'point_3d'.
402 : */
403 0 : void libgrpp_type2_integrals_gradient_(
404 : // contracted Gaussian A
405 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
406 : double *alpha_A, double *origin_B,
407 : // contracted Gaussian B
408 : int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B, double *alpha_B,
409 : // pseudopotential
410 : double *pot_origin, int32_t *pot_L, int32_t *pot_num_primitives,
411 : int32_t *pot_powers, double *pot_coeffs, double *pot_alpha,
412 : // differentiation wrt the 3d point (x,y,z)
413 : double *point_3d,
414 : // answer: matrices d<Int>/dx, d<Int>/dy, d<Int>/dZ
415 : double *grad_arep_x, double *grad_arep_y, double *grad_arep_z) {
416 0 : int *pot_powers_int = (int *)calloc(*pot_num_primitives, sizeof(int));
417 0 : for (int i = 0; i < *pot_num_primitives; i++) {
418 0 : pot_powers_int[i] = pot_powers[i];
419 : }
420 :
421 0 : double *grad_array[3];
422 0 : grad_array[0] = grad_arep_x;
423 0 : grad_array[1] = grad_arep_y;
424 0 : grad_array[2] = grad_arep_z;
425 :
426 0 : libgrpp_potential_t *pot = libgrpp_new_potential(
427 : *pot_L, 0, *pot_num_primitives, pot_powers_int, pot_coeffs, pot_alpha);
428 0 : libgrpp_shell_t *shell_A =
429 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
430 0 : libgrpp_shell_t *shell_B =
431 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
432 :
433 0 : libgrpp_type2_integrals_gradient(shell_A, shell_B, pot_origin, pot, point_3d,
434 : grad_array);
435 :
436 0 : libgrpp_delete_shell(shell_A);
437 0 : libgrpp_delete_shell(shell_B);
438 0 : libgrpp_delete_potential(pot);
439 0 : free(pot_powers_int);
440 0 : }
441 :
442 : /**
443 : * Analytic calculation of gradients of integrals over the effective spin-orbit
444 : * operator (potential) for a given shell pair (with respect to the point
445 : * 'point_3d').
446 : */
447 0 : void libgrpp_spin_orbit_integrals_gradient_(
448 : // contracted Gaussian A
449 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
450 : double *alpha_A,
451 : // contracted Gaussian B
452 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
453 : double *alpha_B,
454 : // pseudopotential
455 : double *pot_origin, int32_t *pot_L, int32_t *pot_num_primitives,
456 : int32_t *pot_powers, double *pot_coeffs, double *pot_alpha,
457 : // differentiation wrt the 3d point (x,y,z)
458 : double *point_3d,
459 : // answer: 9 matrices
460 : // d<SO-x>/dx, d<SO-x>/dy, d<SO-x>/dZ
461 : double *grad_sox_x, double *grad_sox_y, double *grad_sox_z,
462 : // d<SO-y>/dx, d<SO-y>/dy, d<SO-y>/dZ
463 : double *grad_soy_x, double *grad_soy_y, double *grad_soy_z,
464 : // d<SO-z>/dx, d<SO-z>/dy, d<SO-z>/dZ
465 : double *grad_soz_x, double *grad_soz_y, double *grad_soz_z) {
466 0 : int *pot_powers_int = (int *)calloc(*pot_num_primitives, sizeof(int));
467 :
468 0 : double *grad_array_SO_x[3];
469 0 : grad_array_SO_x[0] = grad_sox_x;
470 0 : grad_array_SO_x[1] = grad_sox_y;
471 0 : grad_array_SO_x[2] = grad_sox_z;
472 :
473 0 : double *grad_array_SO_y[3];
474 0 : grad_array_SO_y[0] = grad_soy_x;
475 0 : grad_array_SO_y[1] = grad_soy_y;
476 0 : grad_array_SO_y[2] = grad_soy_z;
477 :
478 0 : double *grad_array_SO_z[3];
479 0 : grad_array_SO_z[0] = grad_soz_x;
480 0 : grad_array_SO_z[1] = grad_soz_y;
481 0 : grad_array_SO_z[2] = grad_soz_z;
482 :
483 0 : for (int i = 0; i < *pot_num_primitives; i++) {
484 0 : pot_powers_int[i] = pot_powers[i];
485 : }
486 :
487 : /*
488 : * construct RPP structure
489 : */
490 0 : libgrpp_potential_t *pot = libgrpp_new_potential(
491 : *pot_L, 0, *pot_num_primitives, pot_powers_int, pot_coeffs, pot_alpha);
492 :
493 : /*
494 : * construct shells
495 : */
496 0 : libgrpp_shell_t *shell_A =
497 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
498 0 : libgrpp_shell_t *shell_B =
499 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
500 :
501 0 : libgrpp_spin_orbit_integrals_gradient(shell_A, shell_B, pot_origin, pot,
502 : point_3d, grad_array_SO_x,
503 : grad_array_SO_y, grad_array_SO_z);
504 :
505 0 : libgrpp_delete_shell(shell_A);
506 0 : libgrpp_delete_shell(shell_B);
507 0 : libgrpp_delete_potential(pot);
508 0 : free(pot_powers_int);
509 0 : }
510 :
511 : /**
512 : * Overlap integrals between two contracted Gaussians with given cartesian parts
513 : * x^n y^l z^m (auxiliary function)
514 : */
515 :
516 0 : void evaluate_overlap_integral_contracted_(
517 : // contracted Gaussian A
518 : double *origin_A, int32_t *n_A, int32_t *l_A, int32_t *m_A,
519 : int32_t *num_primitives_A, double *coeffs_A, double *alpha_A,
520 : // contracted Gaussian B
521 : double *origin_B, int32_t *n_B, int32_t *l_B, int32_t *m_B,
522 : int32_t *num_primitives_B, double *coeffs_B, double *alpha_B,
523 : // answer
524 : double *overlap_integral) {
525 0 : void libgrpp_overlap_integrals(libgrpp_shell_t * shell_A,
526 : libgrpp_shell_t * shell_B,
527 : double *overlap_matrix);
528 :
529 0 : libgrpp_shell_t *shell_A = libgrpp_new_shell(
530 0 : origin_A, *n_A + *l_A + *m_A, *num_primitives_A, coeffs_A, alpha_A);
531 0 : libgrpp_shell_t *shell_B = libgrpp_new_shell(
532 0 : origin_B, *n_B + *l_B + *m_B, *num_primitives_B, coeffs_B, alpha_B);
533 :
534 0 : shell_A->cart_size = 1;
535 0 : shell_A->cart_list[0] = *n_A;
536 0 : shell_A->cart_list[1] = *l_A;
537 0 : shell_A->cart_list[2] = *m_A;
538 :
539 0 : shell_B->cart_size = 1;
540 0 : shell_B->cart_list[0] = *n_B;
541 0 : shell_B->cart_list[1] = *l_B;
542 0 : shell_B->cart_list[2] = *m_B;
543 :
544 0 : libgrpp_overlap_integrals(shell_A, shell_B, overlap_integral);
545 :
546 0 : libgrpp_delete_shell(shell_A);
547 0 : libgrpp_delete_shell(shell_B);
548 0 : }
549 :
550 : /*
551 : * calculates normalization factor for the given contracted Gaussians
552 : * (auxiliary function)
553 : */
554 0 : void radial_gto_norm_factor_(int32_t *L, int32_t *num_primitives,
555 : double *coeffs, double *alpha, double *norm) {
556 0 : *norm = 0.0;
557 0 : double S = 0.0;
558 0 : double origin[] = {0, 0, 0};
559 0 : int zero = 0;
560 :
561 0 : evaluate_overlap_integral_contracted_(origin, L, &zero, &zero, num_primitives,
562 : coeffs, alpha, origin, L, &zero, &zero,
563 : num_primitives, coeffs, alpha, &S);
564 :
565 0 : *norm = sqrt(libgrpp_double_factorial(2 * (*L) - 1)) / sqrt(S);
566 0 : }
567 :
568 : /*
569 : * overlap integrals (for the shell pair)
570 : */
571 :
572 0 : void libgrpp_overlap_integrals_(
573 : // contracted Gaussian A
574 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
575 : double *alpha_A,
576 : // contracted Gaussian B
577 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
578 : double *alpha_B,
579 : // answer
580 : double *matrix) {
581 0 : libgrpp_shell_t *shell_A =
582 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
583 0 : libgrpp_shell_t *shell_B =
584 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
585 :
586 0 : libgrpp_overlap_integrals(shell_A, shell_B, matrix);
587 :
588 0 : libgrpp_delete_shell(shell_A);
589 0 : libgrpp_delete_shell(shell_B);
590 0 : }
591 :
592 : /*
593 : * kinetic-energy integrals (for the shell pair)
594 : */
595 :
596 0 : void libgrpp_kinetic_energy_integrals_(
597 : // contracted Gaussian A
598 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
599 : double *alpha_A,
600 : // contracted Gaussian B
601 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
602 : double *alpha_B,
603 : // answer
604 : double *matrix) {
605 0 : libgrpp_shell_t *shell_A =
606 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
607 0 : libgrpp_shell_t *shell_B =
608 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
609 :
610 0 : libgrpp_kinetic_energy_integrals(shell_A, shell_B, matrix);
611 :
612 0 : libgrpp_delete_shell(shell_A);
613 0 : libgrpp_delete_shell(shell_B);
614 0 : }
615 :
616 : /*
617 : * momentum operator integrals (for the shell pair)
618 : */
619 :
620 0 : void libgrpp_momentum_integrals_(
621 : // contracted Gaussian A
622 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
623 : double *alpha_A,
624 : // contracted Gaussian B
625 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
626 : double *alpha_B,
627 : // answer
628 : double *matrix_x, double *matrix_y, double *matrix_z) {
629 0 : libgrpp_shell_t *shell_A =
630 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
631 0 : libgrpp_shell_t *shell_B =
632 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
633 :
634 0 : libgrpp_momentum_integrals(shell_A, shell_B, matrix_x, matrix_y, matrix_z);
635 :
636 0 : libgrpp_delete_shell(shell_A);
637 0 : libgrpp_delete_shell(shell_B);
638 0 : }
639 :
640 : /*
641 : * nuclear attraction integrals
642 : */
643 :
644 0 : void libgrpp_nuclear_attraction_integrals_(
645 : // contracted Gaussian A
646 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
647 : double *alpha_A,
648 : // contracted Gaussian B
649 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
650 : double *alpha_B,
651 : // potential definition
652 : double *charge_origin, int32_t *charge, int32_t *nuclear_model,
653 : double *model_params,
654 : // answer
655 : double *matrix) {
656 0 : libgrpp_shell_t *shell_A =
657 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
658 0 : libgrpp_shell_t *shell_B =
659 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
660 :
661 0 : libgrpp_nuclear_attraction_integrals(shell_A, shell_B, charge_origin, *charge,
662 : *nuclear_model, model_params, matrix);
663 :
664 0 : libgrpp_delete_shell(shell_A);
665 0 : libgrpp_delete_shell(shell_B);
666 0 : }
667 :
668 0 : void libgrpp_nuclear_attraction_integrals_point_charge_(
669 : // contracted Gaussian A
670 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
671 : double *alpha_A,
672 : // contracted Gaussian B
673 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
674 : double *alpha_B,
675 : // potential definition
676 : double *charge_origin, int32_t *charge,
677 : // answer
678 : double *matrix) {
679 0 : libgrpp_shell_t *shell_A =
680 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
681 0 : libgrpp_shell_t *shell_B =
682 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
683 :
684 0 : libgrpp_nuclear_attraction_integrals_point_charge(
685 : shell_A, shell_B, charge_origin, *charge, matrix);
686 :
687 0 : libgrpp_delete_shell(shell_A);
688 0 : libgrpp_delete_shell(shell_B);
689 0 : }
690 :
691 0 : void libgrpp_nuclear_attraction_integrals_charged_ball_(
692 : // contracted Gaussian A
693 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
694 : double *alpha_A,
695 : // contracted Gaussian B
696 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
697 : double *alpha_B,
698 : // potential definition
699 : double *charge_origin, int32_t *charge, double *r_rms,
700 : // answer
701 : double *matrix) {
702 0 : libgrpp_shell_t *shell_A =
703 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
704 0 : libgrpp_shell_t *shell_B =
705 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
706 :
707 0 : libgrpp_nuclear_attraction_integrals_charged_ball(
708 : shell_A, shell_B, charge_origin, *charge, *r_rms, matrix);
709 :
710 0 : libgrpp_delete_shell(shell_A);
711 0 : libgrpp_delete_shell(shell_B);
712 0 : }
713 :
714 0 : void libgrpp_nuclear_attraction_integrals_gaussian_model_(
715 : // contracted Gaussian A
716 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
717 : double *alpha_A,
718 : // contracted Gaussian B
719 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
720 : double *alpha_B,
721 : // potential definition
722 : double *charge_origin, int32_t *charge, double *r_rms,
723 : // answer
724 : double *matrix) {
725 0 : libgrpp_shell_t *shell_A =
726 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
727 0 : libgrpp_shell_t *shell_B =
728 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
729 :
730 0 : libgrpp_nuclear_attraction_integrals_gaussian_model(
731 : shell_A, shell_B, charge_origin, *charge, *r_rms, matrix);
732 :
733 0 : libgrpp_delete_shell(shell_A);
734 0 : libgrpp_delete_shell(shell_B);
735 0 : }
736 :
737 0 : void libgrpp_nuclear_attraction_integrals_fermi_model_(
738 : // contracted Gaussian A
739 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
740 : double *alpha_A,
741 : // contracted Gaussian B
742 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
743 : double *alpha_B,
744 : // potential definition
745 : double *charge_origin, int32_t *charge, double *fermi_param_c,
746 : double *fermi_param_a,
747 : // answer
748 : double *matrix) {
749 0 : libgrpp_shell_t *shell_A =
750 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
751 0 : libgrpp_shell_t *shell_B =
752 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
753 :
754 0 : libgrpp_nuclear_attraction_integrals_fermi_model(
755 : shell_A, shell_B, charge_origin, *charge, *fermi_param_c, *fermi_param_a,
756 : matrix);
757 :
758 0 : libgrpp_delete_shell(shell_A);
759 0 : libgrpp_delete_shell(shell_B);
760 0 : }
761 :
762 0 : void libgrpp_nuclear_attraction_integrals_fermi_bubble_model_(
763 : // contracted Gaussian A
764 : double *origin_A, int32_t *L_A, int32_t *num_primitives_A, double *coeffs_A,
765 : double *alpha_A,
766 : // contracted Gaussian B
767 : double *origin_B, int32_t *L_B, int32_t *num_primitives_B, double *coeffs_B,
768 : double *alpha_B,
769 : // potential definition
770 : double *charge_origin, int32_t *charge, double *fermi_param_c,
771 : double *fermi_param_a, double *param_k,
772 : // answer
773 : double *matrix) {
774 0 : libgrpp_shell_t *shell_A =
775 0 : libgrpp_new_shell(origin_A, *L_A, *num_primitives_A, coeffs_A, alpha_A);
776 0 : libgrpp_shell_t *shell_B =
777 0 : libgrpp_new_shell(origin_B, *L_B, *num_primitives_B, coeffs_B, alpha_B);
778 :
779 0 : libgrpp_nuclear_attraction_integrals_fermi_bubble_model(
780 : shell_A, shell_B, charge_origin, *charge, *fermi_param_c, *fermi_param_a,
781 : *param_k, matrix);
782 :
783 0 : libgrpp_delete_shell(shell_A);
784 0 : libgrpp_delete_shell(shell_B);
785 0 : }
786 :
787 : /*
788 : * Fortran interface to the nuclear models
789 : */
790 :
791 0 : void libgrpp_estimate_nuclear_rms_radius_johnson_1985_(int32_t *A,
792 : double *R_rms) {
793 0 : *R_rms = libgrpp_estimate_nuclear_rms_radius_johnson_1985(*A);
794 0 : }
795 :
796 0 : void libgrpp_estimate_nuclear_rms_radius_golovko_2008_(int32_t *A,
797 : double *R_rms) {
798 0 : *R_rms = libgrpp_estimate_nuclear_rms_radius_golovko_2008(*A);
799 0 : }
800 :
801 0 : void libgrpp_estimate_fermi_model_parameters_(double *R_rms, double *c,
802 : double *a, int32_t *err_code) {
803 0 : *err_code = (int32_t)libgrpp_estimate_fermi_model_parameters(*R_rms, c, a);
804 0 : }
805 :
806 0 : void libgrpp_charge_density_ball_(double *r, double *Z, double *R_rms,
807 : double *rho) {
808 0 : *rho = libgrpp_charge_density_ball(*r, *Z, *R_rms);
809 0 : }
810 :
811 0 : void libgrpp_charge_density_gaussian_(double *r, double *Z, double *R_rms,
812 : double *rho) {
813 0 : *rho = libgrpp_charge_density_gaussian(*r, *Z, *R_rms);
814 0 : }
815 :
816 0 : void libgrpp_charge_density_fermi_(double *r, double *Z, double *c, double *a,
817 : double *rho) {
818 0 : *rho = libgrpp_charge_density_fermi(*r, *Z, *c, *a);
819 0 : }
820 :
821 0 : void libgrpp_charge_density_fermi_bubble_(double *r, double *Z, double *c,
822 : double *a, double *k, double *rho) {
823 0 : *rho = libgrpp_charge_density_fermi_bubble(*r, *Z, *c, *a, *k);
824 0 : }
825 :
826 0 : void libgrpp_coulomb_potential_point_(double *r, double *Z, double *potential) {
827 0 : *potential = libgrpp_coulomb_potential_point(*r, *Z);
828 0 : }
829 :
830 0 : void libgrpp_coulomb_potential_ball_(double *r, double *Z, double *R_rms,
831 : double *potential) {
832 0 : *potential = libgrpp_coulomb_potential_ball(*r, *Z, *R_rms);
833 0 : }
834 :
835 0 : void libgrpp_coulomb_potential_gaussian_(double *r, double *Z, double *R_rms,
836 : double *potential) {
837 0 : *potential = libgrpp_coulomb_potential_gaussian(*r, *Z, *R_rms);
838 0 : }
839 :
840 0 : void libgrpp_coulomb_potential_fermi_(double *r, double *Z, double *c,
841 : double *a, double *potential) {
842 0 : *potential = libgrpp_coulomb_potential_fermi(*r, *Z, *c, *a);
843 0 : }
844 :
845 0 : void libgrpp_coulomb_potential_fermi_bubble_(double *r, double *Z, double *c,
846 : double *a, double *k,
847 : double *potential) {
848 0 : *potential = libgrpp_coulomb_potential_fermi_bubble(*r, *Z, *c, *a, *k);
849 0 : }
850 :
851 0 : void libgrpp_rms_radius_fermi_(int32_t *Z, double *c, double *a,
852 : double *r_rms) {
853 0 : *r_rms = libgrpp_rms_radius_fermi(*Z, *c, *a);
854 0 : }
855 :
856 0 : void libgrpp_rms_radius_fermi_bubble_(int32_t *Z, double *c, double *a,
857 : double *k, double *r_rms) {
858 0 : *r_rms = libgrpp_rms_radius_fermi_bubble(*Z, *c, *a, *k);
859 0 : }
|
