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 : * Integration of the non-local terms in the GRPP operator.
17 : * These integrals are reduced to the type 1 integrals and overlap integrals.
18 : */
19 :
20 : #include <assert.h>
21 : #include <math.h>
22 : #include <stdio.h>
23 : #include <stdlib.h>
24 : #include <string.h>
25 :
26 : #ifndef M_PI
27 : #define M_PI 3.14159265358979323846
28 : #endif
29 :
30 : #include "grpp_factorial.h"
31 : #include "grpp_lmatrix.h"
32 : #include "grpp_overlap.h"
33 : #include "grpp_spherical_harmonics.h"
34 : #include "grpp_utils.h"
35 : #include "libgrpp.h"
36 :
37 : /*
38 : * pre-definitions of function used below
39 : */
40 :
41 : void libgrpp_outercore_potential_integrals_part_1(
42 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *C,
43 : libgrpp_potential_t *oc_potential, libgrpp_shell_t *oc_shell,
44 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
45 : double *so_z_matrix);
46 :
47 : void libgrpp_outercore_potential_integrals_part_2(
48 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *C,
49 : libgrpp_potential_t *oc_potential_1, libgrpp_shell_t *oc_shell_1,
50 : libgrpp_potential_t *oc_potential_2, libgrpp_shell_t *oc_shell_2,
51 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
52 : double *so_z_matrix);
53 :
54 : static double calculate_delta_integral(libgrpp_potential_t *oc_pot_1,
55 : libgrpp_shell_t *oc_shell_1,
56 : libgrpp_potential_t *oc_pot_2,
57 : libgrpp_shell_t *oc_shell_2);
58 :
59 : static void transform_to_sph_basis_ket(int dim_bra, int dim_ket_cart,
60 : int dim_ket_sph, double *A_in,
61 : double *A_out, double *S_lm_coef);
62 :
63 : static void transform_to_sph_basis_bra(int dim_bra_cart, int dim_bra_sph,
64 : int dim_ket, double *A_in, double *A_out,
65 : double *S_lm_coef);
66 :
67 : static double ang_norm_factor(int lx, int ly, int lz);
68 :
69 : static double analytic_one_center_rpp_integral_contracted(
70 : libgrpp_shell_t *bra, libgrpp_shell_t *ket, libgrpp_potential_t *pot);
71 :
72 : static double analytic_one_center_rpp_integral_primitive(int L, double alpha1,
73 : double alpha2, int n,
74 : double zeta);
75 :
76 : static double radial_gto_norm_factor(int L, double alpha);
77 :
78 : /**
79 : * Calculates non-local contributions to the scalar-relativistic ECP and
80 : * effective spin-orbit interaction matrices from the outercore (OC) potentials.
81 : */
82 0 : void libgrpp_outercore_potential_integrals(
83 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *rpp_origin,
84 : int num_oc_shells, libgrpp_potential_t **oc_potentials,
85 : libgrpp_shell_t **oc_shells, double *arep, double *esop_x, double *esop_y,
86 : double *esop_z) {
87 0 : assert(libgrpp_is_initialized());
88 :
89 : // clear output matrices
90 0 : int size_A = libgrpp_get_shell_size(shell_A);
91 0 : int size_B = libgrpp_get_shell_size(shell_B);
92 0 : memset(arep, 0, size_A * size_B * sizeof(double));
93 0 : memset(esop_x, 0, size_A * size_B * sizeof(double));
94 0 : memset(esop_y, 0, size_A * size_B * sizeof(double));
95 0 : memset(esop_z, 0, size_A * size_B * sizeof(double));
96 :
97 : // bind outercore shells of the GRPP to the RPP center
98 0 : for (int ioc = 0; ioc < num_oc_shells; ioc++) {
99 0 : oc_shells[ioc]->origin[0] = rpp_origin[0];
100 0 : oc_shells[ioc]->origin[1] = rpp_origin[1];
101 0 : oc_shells[ioc]->origin[2] = rpp_origin[2];
102 : }
103 :
104 : // 1. the U * |nlj><nlj| + |nlj><nlj| * U part
105 0 : for (int ioc = 0; ioc < num_oc_shells; ioc++) {
106 0 : libgrpp_potential_t *pot = oc_potentials[ioc];
107 0 : libgrpp_shell_t *nlj = oc_shells[ioc];
108 :
109 0 : libgrpp_outercore_potential_integrals_part_1(
110 : shell_A, shell_B, rpp_origin, pot, nlj, arep, esop_x, esop_y, esop_z);
111 : }
112 :
113 : // 2. the |nlj><nlj| U |n'lj><n'lj| part
114 0 : for (int ioc = 0; ioc < num_oc_shells; ioc++) {
115 0 : for (int joc = 0; joc < num_oc_shells; joc++) {
116 :
117 0 : libgrpp_potential_t *pot_1 = oc_potentials[ioc];
118 0 : libgrpp_potential_t *pot_2 = oc_potentials[joc];
119 0 : libgrpp_shell_t *nlj_1 = oc_shells[ioc];
120 0 : libgrpp_shell_t *nlj_2 = oc_shells[joc];
121 :
122 0 : libgrpp_outercore_potential_integrals_part_2(
123 : shell_A, shell_B, rpp_origin, pot_1, nlj_1, pot_2, nlj_2, arep,
124 : esop_x, esop_y, esop_z);
125 : }
126 : }
127 0 : }
128 :
129 : /**
130 : * integration of the non-local outercore potential:
131 : * the U*|nlj><nlj| + |nlj><nlj|*U part
132 : */
133 0 : void libgrpp_outercore_potential_integrals_part_1(
134 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *C,
135 : libgrpp_potential_t *oc_potential, libgrpp_shell_t *oc_shell,
136 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
137 : double *so_z_matrix) {
138 0 : int L = oc_shell->L;
139 0 : double J = oc_potential->J / 2.0;
140 0 : int size_A = libgrpp_get_shell_size(shell_A);
141 0 : int size_B = libgrpp_get_shell_size(shell_B);
142 0 : int size_nlj = libgrpp_get_shell_size(oc_shell);
143 :
144 : /*
145 : * foolproof: bind outercore shells to the ECP center
146 : */
147 0 : oc_shell->origin[0] = C[0];
148 0 : oc_shell->origin[1] = C[1];
149 0 : oc_shell->origin[2] = C[2];
150 :
151 : /*
152 : * Transformation coefficients: Cartesian -> Real spherical
153 : * Rows: real spherical harmonics S_lm
154 : * Columns: cartesians x^r y^s z^t
155 : */
156 0 : double *S_lm_coef = (double *)calloc((2 * L + 1) * size_nlj, sizeof(double));
157 0 : for (int m = -L; m <= +L; m++) {
158 0 : for (int icart = 0; icart < size_nlj; icart++) {
159 0 : int r = oc_shell->cart_list[3 * icart + 0];
160 0 : int s = oc_shell->cart_list[3 * icart + 1];
161 0 : int t = oc_shell->cart_list[3 * icart + 2];
162 0 : double u = libgrpp_spherical_to_cartesian_coef(L, m, r, s) /
163 0 : ang_norm_factor(r, s, t);
164 0 : S_lm_coef[size_nlj * (m + L) + icart] = u;
165 : }
166 : }
167 :
168 : /*
169 : * Overlap integrals of the A and B shells with the outercore shell |nlj>:
170 : * <A|nljm> and <nljm|B>
171 : *
172 : * Integrals are evaluated first for the cartesian components of |nlj>, and
173 : * then transformed to the basis of real spherical components |nljm>.
174 : * Resulting integrals are stored in the 'S_a_nljm' and 'S_nljm_b' arrays.
175 : */
176 0 : double *S_nljm_b_cart =
177 0 : alloc_zeros_1d(size_nlj * size_B); // in cart_list basis
178 0 : double *S_a_nljm_cart = alloc_zeros_1d(size_A * size_nlj);
179 0 : double *S_nljm_b = alloc_zeros_1d((2 * L + 1) * size_B); // in spherical basis
180 0 : double *S_a_nljm = alloc_zeros_1d(size_A * (2 * L + 1));
181 0 : libgrpp_overlap_integrals(oc_shell, shell_B, S_nljm_b_cart); // <nljm|B>
182 0 : libgrpp_overlap_integrals(shell_A, oc_shell, S_a_nljm_cart); // <A|nljm>
183 0 : transform_to_sph_basis_ket(size_A, size_nlj, 2 * L + 1, S_a_nljm_cart,
184 : S_a_nljm, S_lm_coef);
185 0 : transform_to_sph_basis_bra(size_nlj, 2 * L + 1, size_B, S_nljm_b_cart,
186 : S_nljm_b, S_lm_coef);
187 0 : free(S_nljm_b_cart);
188 0 : free(S_a_nljm_cart);
189 :
190 : /*
191 : * ECP type-2 (semilocal) integrals of the A and B shells with the outercore
192 : * shell |nlj>: <A|U(r)P_L|nljm> and <nljm|U(r)P_L|B>
193 : *
194 : * Integrals are evaluated first for the cartesian components of |nlj>, and
195 : * then transformed to the basis of real spherical components |nljm>.
196 : * Resulting integrals are stored in the 'U_a_nljm' and 'U_nljm_b' arrays.
197 : */
198 0 : double *U_nljm_b_cart =
199 0 : alloc_zeros_1d(size_nlj * size_B); // in cart_list basis
200 0 : double *U_a_nljm_cart = alloc_zeros_1d(size_A * size_nlj);
201 0 : double *U_nljm_b = alloc_zeros_1d((2 * L + 1) * size_B); // in spherical basis
202 0 : double *U_a_nljm = alloc_zeros_1d(size_A * (2 * L + 1));
203 0 : libgrpp_type1_integrals(oc_shell, shell_B, C, oc_potential,
204 : U_nljm_b_cart); // <nljm|U(r)P_L|B>
205 0 : libgrpp_type1_integrals(shell_A, oc_shell, C, oc_potential,
206 : U_a_nljm_cart); // <A|U(r)P_L|nljm>
207 0 : transform_to_sph_basis_ket(size_A, size_nlj, 2 * L + 1, U_a_nljm_cart,
208 : U_a_nljm, S_lm_coef);
209 0 : transform_to_sph_basis_bra(size_nlj, 2 * L + 1, size_B, U_nljm_b_cart,
210 : U_nljm_b, S_lm_coef);
211 0 : free(U_nljm_b_cart);
212 0 : free(U_a_nljm_cart);
213 :
214 : /*
215 : * Construct outercore AREP matrix elements
216 : * < a | P_nlj U(r) P_L + U(r) P_nlj P_L | b >
217 : */
218 0 : double arep_factor =
219 0 : (J < L) ? (L / (2.0 * L + 1)) : ((L + 1) / (2.0 * L + 1));
220 0 : double *buf = alloc_zeros_1d(size_A * size_B);
221 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, S_a_nljm, U_nljm_b, buf);
222 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, U_a_nljm, S_nljm_b, buf);
223 0 : libgrpp_daxpy(size_A * size_B, arep_factor, buf, arep_matrix);
224 0 : free(buf);
225 :
226 : /*
227 : * Construct outercore effective SO potential matrix elements
228 : * < a | P_nlj U(r) L P_L + U(r) P_nlj L P_L | b >
229 : */
230 0 : double **L_matrices = alloc_zeros_2d(3, (2 * L + 1) * (2 * L + 1));
231 0 : libgrpp_construct_angular_momentum_matrices_rsh(L, L_matrices[0],
232 : L_matrices[1], L_matrices[2]);
233 :
234 0 : double **so_buf = alloc_zeros_2d(3, size_A * size_B);
235 0 : buf = alloc_zeros_1d((2 * L + 1) * int_max2(size_A, size_B));
236 :
237 0 : for (int icoord = 0; icoord < 3; icoord++) {
238 : // U(r) P_nlj
239 0 : memset(buf, 0, (2 * L + 1) * size_B * sizeof(double));
240 0 : libgrpp_multiply_matrices(2 * L + 1, size_B, 2 * L + 1, L_matrices[icoord],
241 : S_nljm_b, buf);
242 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, U_a_nljm, buf,
243 0 : so_buf[icoord]);
244 :
245 : // P_nlj U(r)
246 0 : memset(buf, 0, (2 * L + 1) * size_A * sizeof(double));
247 0 : libgrpp_multiply_matrices(size_A, 2 * L + 1, 2 * L + 1, S_a_nljm,
248 : L_matrices[icoord], buf);
249 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, buf, U_nljm_b,
250 : so_buf[icoord]);
251 : }
252 :
253 0 : free(buf);
254 :
255 0 : double esop_factor = (J < L) ? (-2.0 / (2 * L + 1)) : (+2.0 / (2 * L + 1));
256 0 : libgrpp_daxpy(size_A * size_B, esop_factor, so_buf[0], so_x_matrix);
257 0 : libgrpp_daxpy(size_A * size_B, esop_factor, so_buf[1], so_y_matrix);
258 0 : libgrpp_daxpy(size_A * size_B, esop_factor, so_buf[2], so_z_matrix);
259 :
260 : /*
261 : * Cleanup
262 : */
263 0 : free(S_lm_coef);
264 0 : free(S_a_nljm);
265 0 : free(S_nljm_b);
266 0 : free(U_a_nljm);
267 0 : free(U_nljm_b);
268 0 : free_2d(L_matrices, 3);
269 0 : free_2d(so_buf, 3);
270 0 : }
271 :
272 : /**
273 : * integration of the non-local outercore potential:
274 : * the |nlj><nlj| U |n'lj><n'lj| part
275 : */
276 0 : void libgrpp_outercore_potential_integrals_part_2(
277 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double *C,
278 : libgrpp_potential_t *oc_potential_1, libgrpp_shell_t *oc_shell_1,
279 : libgrpp_potential_t *oc_potential_2, libgrpp_shell_t *oc_shell_2,
280 : double *arep_matrix, double *so_x_matrix, double *so_y_matrix,
281 : double *so_z_matrix) {
282 0 : int size_A = libgrpp_get_shell_size(shell_A);
283 0 : int size_B = libgrpp_get_shell_size(shell_B);
284 :
285 0 : if (oc_potential_1->L != oc_potential_2->L) {
286 : return;
287 : }
288 0 : if (oc_potential_1->J != oc_potential_2->J) {
289 : return;
290 : }
291 :
292 : /*
293 : * foolproof: bind outercore shells to the ECP center
294 : */
295 0 : oc_shell_1->origin[0] = C[0];
296 0 : oc_shell_1->origin[1] = C[1];
297 0 : oc_shell_1->origin[2] = C[2];
298 0 : oc_shell_2->origin[0] = C[0];
299 0 : oc_shell_2->origin[1] = C[1];
300 0 : oc_shell_2->origin[2] = C[2];
301 :
302 : /*
303 : * just to be consistent with the MOLGEP code by A. Titov, N. Mosyagin, A.
304 : * Petrov: off-diagonal elements <n'lj|U|n''lj> = 0
305 : */
306 : /*if (!ecps_are_equal(oc_potential_1, oc_potential_2)) {
307 : return;
308 : }*/
309 :
310 0 : int L = oc_potential_1->L;
311 0 : double J = oc_potential_1->J / 2.0;
312 0 : int size_nlj = libgrpp_get_shell_size(oc_shell_1);
313 :
314 0 : double delta = calculate_delta_integral(oc_potential_1, oc_shell_1,
315 : oc_potential_2, oc_shell_2);
316 :
317 : /*
318 : * Transformation coefficients: Cartesian -> Real spherical
319 : * Rows: real spherical harmonics S_lm
320 : * Columns: cartesians x^r y^s z^t
321 : */
322 0 : double *S_lm_coef = (double *)calloc((2 * L + 1) * size_nlj, sizeof(double));
323 0 : for (int m = -L; m <= +L; m++) {
324 0 : for (int icart = 0; icart < size_nlj; icart++) {
325 0 : int r = oc_shell_1->cart_list[3 * icart + 0];
326 0 : int s = oc_shell_1->cart_list[3 * icart + 1];
327 0 : int t = oc_shell_1->cart_list[3 * icart + 2];
328 0 : double u = libgrpp_spherical_to_cartesian_coef(L, m, r, s) /
329 0 : ang_norm_factor(r, s, t);
330 0 : S_lm_coef[size_nlj * (m + L) + icart] = u;
331 : }
332 : }
333 :
334 : /*
335 : * Overlap integrals of the A and B shells with the outercore shells |nlj> and
336 : * |n'lj>: <A|nljm> and <n'ljm|B>
337 : *
338 : * Integrals are evaluated first for the cartesian components of |nlj>, and
339 : * then transformed to the basis of real spherical components |nljm>.
340 : * Resulting integrals are stored in the 'S_a_nljm1' and 'S_nljm2_b' arrays.
341 : */
342 0 : double *S_nljm2_b_cart =
343 0 : alloc_zeros_1d(size_nlj * size_B); // in cart_list basis
344 0 : double *S_a_nljm1_cart = alloc_zeros_1d(size_A * size_nlj);
345 0 : double *S_nljm2_b =
346 0 : alloc_zeros_1d((2 * L + 1) * size_B); // in spherical basis
347 0 : double *S_a_nljm1 = alloc_zeros_1d(size_A * (2 * L + 1));
348 0 : libgrpp_overlap_integrals(oc_shell_2, shell_B, S_nljm2_b_cart); // <nljm|B>
349 0 : libgrpp_overlap_integrals(shell_A, oc_shell_1, S_a_nljm1_cart); // <A|nljm>
350 0 : transform_to_sph_basis_ket(size_A, size_nlj, 2 * L + 1, S_a_nljm1_cart,
351 : S_a_nljm1, S_lm_coef);
352 0 : transform_to_sph_basis_bra(size_nlj, 2 * L + 1, size_B, S_nljm2_b_cart,
353 : S_nljm2_b, S_lm_coef);
354 0 : free(S_nljm2_b_cart);
355 0 : free(S_a_nljm1_cart);
356 :
357 : /*
358 : * Construct outercore AREP matrix elements
359 : * < a | P_nlj U(r) P_n'lj | b >
360 : */
361 0 : double arep_factor =
362 0 : (J < L) ? (L / (2.0 * L + 1)) : ((L + 1) / (2.0 * L + 1));
363 0 : double *buf = alloc_zeros_1d(size_A * size_B);
364 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, S_a_nljm1, S_nljm2_b,
365 : buf);
366 0 : libgrpp_daxpy(size_A * size_B, (-1.0) * delta * arep_factor, buf,
367 : arep_matrix);
368 0 : free(buf);
369 :
370 : /*
371 : * Construct outercore effective SO potential matrix elements
372 : * < a | P_nlj U(r) L P_L + U(r) P_nlj L P_L | b >
373 : */
374 0 : double **L_matrices = alloc_zeros_2d(3, (2 * L + 1) * (2 * L + 1));
375 0 : libgrpp_construct_angular_momentum_matrices_rsh(L, L_matrices[0],
376 : L_matrices[1], L_matrices[2]);
377 :
378 0 : double **so_buf = alloc_zeros_2d(3, size_A * size_B);
379 0 : buf = alloc_zeros_1d((2 * L + 1) * size_B);
380 0 : for (int icoord = 0; icoord < 3; icoord++) {
381 0 : memset(buf, 0, (2 * L + 1) * size_B * sizeof(double));
382 0 : libgrpp_multiply_matrices(2 * L + 1, size_B, 2 * L + 1, L_matrices[icoord],
383 : S_nljm2_b, buf);
384 0 : libgrpp_multiply_matrices(size_A, size_B, 2 * L + 1, S_a_nljm1, buf,
385 0 : so_buf[icoord]);
386 : }
387 0 : free(buf);
388 :
389 0 : double esop_factor = (J < L) ? (-2.0 / (2 * L + 1)) : (+2.0 / (2 * L + 1));
390 0 : libgrpp_daxpy(size_A * size_B, (-1.0) * delta * esop_factor, so_buf[0],
391 : so_x_matrix);
392 0 : libgrpp_daxpy(size_A * size_B, (-1.0) * delta * esop_factor, so_buf[1],
393 : so_y_matrix);
394 0 : libgrpp_daxpy(size_A * size_B, (-1.0) * delta * esop_factor, so_buf[2],
395 : so_z_matrix);
396 :
397 0 : free_2d(so_buf, 3);
398 0 : free_2d(L_matrices, 3);
399 :
400 0 : free(S_nljm2_b);
401 0 : free(S_a_nljm1);
402 0 : free(S_lm_coef);
403 : }
404 :
405 : /**
406 : * Calculation of radial "delta" integrals.
407 : * Is performed analytically.
408 : */
409 0 : static double calculate_delta_integral(libgrpp_potential_t *oc_pot_1,
410 : libgrpp_shell_t *oc_shell_1,
411 : libgrpp_potential_t *oc_pot_2,
412 : libgrpp_shell_t *oc_shell_2) {
413 : // both shells must have equal L,J quantum numbers, otherwise
414 : // the < nlj | U | n'l'j' > integral is strictly zero
415 0 : if (oc_pot_1->L != oc_pot_2->L || oc_pot_1->J != oc_pot_2->J) {
416 : return 0.0;
417 : }
418 :
419 0 : double U1 = analytic_one_center_rpp_integral_contracted(oc_shell_1,
420 : oc_shell_2, oc_pot_1);
421 0 : double U2 = analytic_one_center_rpp_integral_contracted(oc_shell_1,
422 : oc_shell_2, oc_pot_2);
423 :
424 0 : return 0.5 * (U1 + U2);
425 : }
426 :
427 : /**
428 : * analytic formula for one-center RPP integral between two contracted gaussian
429 : * functions.
430 : */
431 0 : static double analytic_one_center_rpp_integral_contracted(
432 : libgrpp_shell_t *bra, libgrpp_shell_t *ket, libgrpp_potential_t *pot) {
433 0 : double sum = 0.0;
434 0 : int L = pot->L;
435 :
436 0 : for (int i = 0; i < bra->num_primitives; i++) {
437 0 : double coef_i = bra->coeffs[i];
438 0 : double alpha_i = bra->alpha[i];
439 0 : double N_i = radial_gto_norm_factor(L, alpha_i);
440 :
441 0 : for (int j = 0; j < ket->num_primitives; j++) {
442 0 : double coef_j = ket->coeffs[j];
443 0 : double alpha_j = ket->alpha[j];
444 0 : double N_j = radial_gto_norm_factor(L, alpha_j);
445 0 : double factor = N_i * N_j * coef_i * coef_j;
446 :
447 0 : for (int k = 0; k < pot->num_primitives; k++) {
448 0 : double coef_k = pot->coeffs[k];
449 0 : double zeta = pot->alpha[k];
450 0 : int n_rpp = pot->powers[k];
451 :
452 0 : double u_ijk = analytic_one_center_rpp_integral_primitive(
453 : L, alpha_i, alpha_j, n_rpp, zeta);
454 :
455 0 : sum += factor * coef_k * u_ijk;
456 : }
457 : }
458 : }
459 :
460 0 : return sum;
461 : }
462 :
463 : /**
464 : * analytic formula for one-center RPP integral between two gaussian primitives.
465 : * normalization factors are omitted here.
466 : */
467 0 : static double analytic_one_center_rpp_integral_primitive(int L, double alpha1,
468 : double alpha2, int n,
469 : double zeta) {
470 0 : double a = alpha1 + alpha2 + zeta;
471 :
472 0 : if (n % 2 == 0) { // even n
473 0 : int k = L + n / 2;
474 0 : return libgrpp_double_factorial(2 * k - 1) / (pow(2.0, k + 1) * pow(a, k)) *
475 0 : sqrt(M_PI / a);
476 : } else { // odd n
477 0 : int k = L + (n - 1) / 2;
478 0 : return libgrpp_factorial(k) / (2.0 * pow(a, k + 1));
479 : }
480 : }
481 :
482 : /**
483 : * calculate normalization factor for the radial Gaussian-type orbital:
484 : * G(r) = N * r^L * exp(-alpha * r^2)
485 : */
486 0 : static double radial_gto_norm_factor(int L, double alpha) {
487 : // pre-tabulated factors for calculation of normalization constants
488 : // (for each value of L)
489 0 : static const double factors[] = {
490 : 2.5264751109842587, 2.9173221708553032, 2.6093322745198853,
491 : 1.9724697960897537, 1.3149798640598356, 7.9296269381073192e-1,
492 : 4.3985656185609934e-1, 2.2714095183849672e-1, 1.1017954545099481e-1,
493 : 5.0553842554329785e-2, 2.2063505731056757e-2, 9.2011179391124215e-3,
494 : 3.6804471756449694e-3, 1.4166047783978804e-3, 5.2611380677564405e-4,
495 : 1.8898565833279173e-4, 6.5796360823633550e-5, 2.2243229718298637e-5,
496 : 7.3135288801774484e-6, 2.3422037547660024e-6};
497 :
498 0 : return factors[L] * pow(alpha, 0.75 + L / 2.0);
499 : }
500 :
501 : /**
502 : * Transforms matrix from the basis of unitary sphere polynomials
503 : * to the basis of real spherical harmonics S_lm
504 : * (separately for 'bra' and 'ket' vectors)
505 : */
506 0 : static void transform_to_sph_basis_ket(int dim_bra, int dim_ket_cart,
507 : int dim_ket_sph, double *A_in,
508 : double *A_out, double *S_lm_coef) {
509 0 : for (int i = 0; i < dim_bra; i++) {
510 0 : for (int j = 0; j < dim_ket_sph; j++) {
511 : double s = 0.0;
512 :
513 0 : for (int icart = 0; icart < dim_ket_cart; icart++) {
514 0 : double u_nljm = S_lm_coef[dim_ket_cart * j + icart];
515 0 : s += u_nljm * A_in[i * dim_ket_cart + icart];
516 : }
517 :
518 0 : A_out[i * dim_ket_sph + j] = s;
519 : }
520 : }
521 0 : }
522 :
523 0 : static void transform_to_sph_basis_bra(int dim_bra_cart, int dim_bra_sph,
524 : int dim_ket, double *A_in, double *A_out,
525 : double *S_lm_coef) {
526 0 : for (int j = 0; j < dim_ket; j++) {
527 0 : for (int i = 0; i < dim_bra_sph; i++) {
528 :
529 : double s = 0.0;
530 :
531 0 : for (int icart = 0; icart < dim_bra_cart; icart++) {
532 0 : double u_nljm = S_lm_coef[dim_bra_cart * i + icart];
533 0 : s += u_nljm * A_in[icart * dim_ket + j];
534 : }
535 :
536 0 : A_out[i * dim_ket + j] = s;
537 : }
538 : }
539 0 : }
540 :
541 0 : static double ang_norm_factor(int lx, int ly, int lz) {
542 0 : int L = lx + ly + lz;
543 :
544 0 : return 1.0 / (2.0 * sqrt(M_PI)) *
545 0 : sqrt(libgrpp_double_factorial(2 * L + 1)
546 : /*(double_factorial(2 * lx - 1) * double_factorial(2 * ly - 1) *
547 : double_factorial(2 * lz - 1))*/
548 : );
549 : }
|
