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 <assert.h>
16 : #include <math.h>
17 : #include <stdlib.h>
18 : #include <string.h>
19 :
20 : #ifndef M_PI
21 : #define M_PI 3.14159265358979323846
22 : #endif
23 :
24 : #include "grpp_angular_integrals.h"
25 : #include "grpp_binomial.h"
26 : #include "grpp_lmatrix.h"
27 : #include "grpp_radial_type2_integral.h"
28 : #include "grpp_spherical_harmonics.h"
29 : #include "grpp_utils.h"
30 : #include "libgrpp.h"
31 :
32 : #define LMAX (2 * LIBGRPP_MAX_BASIS_L + LIBGRPP_MAX_RPP_L)
33 :
34 : static void type3_angular_sum(int L, double *Lx_matrix, double *Ly_matrix,
35 : double *Lz_matrix, int lambda_1, int a, int b,
36 : int c, double *rsh_values_kA, int lambda_2, int d,
37 : int e, int f, double *rsh_values_kB,
38 : double *sum_angular_x, double *sum_angular_y,
39 : double *sum_angular_z);
40 :
41 : /**
42 : * Evaluation of spin-orbit ("type 3") RPP integrals.
43 : *
44 : * The theoretical outline is given in the paper:
45 : * R. M. Pitzer, N. W. Winter. Spin-orbit (core) and core potential integrals.
46 : * Int. J. Quantum Chem. 40(6), 773 (1991). doi: 10.1002/qua.560400606
47 : * However, the formula on page 776 of Pitzer & Winter is not reproduced in the
48 : * code exactly.
49 : */
50 0 : void libgrpp_spin_orbit_integrals(libgrpp_shell_t *shell_A,
51 : libgrpp_shell_t *shell_B, double *rpp_origin,
52 : libgrpp_potential_t *potential,
53 : double *so_x_matrix, double *so_y_matrix,
54 : double *so_z_matrix) {
55 0 : assert(libgrpp_is_initialized());
56 :
57 0 : int size_A = libgrpp_get_shell_size(shell_A);
58 0 : int size_B = libgrpp_get_shell_size(shell_B);
59 :
60 0 : memset(so_x_matrix, 0, size_A * size_B * sizeof(double));
61 0 : memset(so_y_matrix, 0, size_A * size_B * sizeof(double));
62 0 : memset(so_z_matrix, 0, size_A * size_B * sizeof(double));
63 :
64 0 : int L = potential->L;
65 0 : int L_A =
66 0 : shell_A->cart_list[0] + shell_A->cart_list[1] + shell_A->cart_list[2];
67 0 : int L_B =
68 0 : shell_B->cart_list[0] + shell_B->cart_list[1] + shell_B->cart_list[2];
69 :
70 0 : double *A = shell_A->origin;
71 0 : double *B = shell_B->origin;
72 0 : double *C = rpp_origin;
73 :
74 0 : double CA_x = C[0] - A[0];
75 0 : double CA_y = C[1] - A[1];
76 0 : double CA_z = C[2] - A[2];
77 0 : double CB_x = C[0] - B[0];
78 0 : double CB_y = C[1] - B[1];
79 0 : double CB_z = C[2] - B[2];
80 0 : double CA_2 = CA_x * CA_x + CA_y * CA_y + CA_z * CA_z;
81 0 : double CB_2 = CB_x * CB_x + CB_y * CB_y + CB_z * CB_z;
82 :
83 0 : double alpha_A = shell_A->alpha[0];
84 0 : double alpha_B = shell_B->alpha[0];
85 0 : double kA_x = -2.0 * (alpha_A * CA_x);
86 0 : double kA_y = -2.0 * (alpha_A * CA_y);
87 0 : double kA_z = -2.0 * (alpha_A * CA_z);
88 0 : double kB_x = -2.0 * (alpha_B * CB_x);
89 0 : double kB_y = -2.0 * (alpha_B * CB_y);
90 0 : double kB_z = -2.0 * (alpha_B * CB_z);
91 0 : double kA_vec[3];
92 0 : kA_vec[0] = kA_x;
93 0 : kA_vec[1] = kA_y;
94 0 : kA_vec[2] = kA_z;
95 0 : double kB_vec[3];
96 0 : kB_vec[0] = kB_x;
97 0 : kB_vec[1] = kB_y;
98 0 : kB_vec[2] = kB_z;
99 :
100 0 : int lambda1_max = L + L_A;
101 0 : int lambda2_max = L + L_B;
102 0 : int N_max = L_A + L_B; // + n_RPP;
103 :
104 : /*
105 : * pre-compute matrices of the Lx, Ly, Lz operators
106 : */
107 0 : double *Lx_matrix = calloc((2 * L + 1) * (2 * L + 1), sizeof(double));
108 0 : double *Ly_matrix = calloc((2 * L + 1) * (2 * L + 1), sizeof(double));
109 0 : double *Lz_matrix = calloc((2 * L + 1) * (2 * L + 1), sizeof(double));
110 0 : libgrpp_construct_angular_momentum_matrices_rsh(L, Lx_matrix, Ly_matrix,
111 : Lz_matrix);
112 :
113 : /*
114 : * for further evaluation of angular integrals
115 : */
116 0 : int lmax = int_max3(lambda1_max, lambda2_max, L);
117 : // create_real_spherical_harmonic_coeffs_tables(lmax);
118 :
119 : /*
120 : * pre-calculate values of real spherical harmonics for different L
121 : */
122 0 : double rsh_values_kA[LMAX][2 * LMAX + 1];
123 0 : double rsh_values_kB[LMAX][2 * LMAX + 1];
124 :
125 0 : for (int lambda = 0; lambda <= lmax; lambda++) {
126 0 : libgrpp_evaluate_real_spherical_harmonics_array(lambda, kA_vec,
127 0 : rsh_values_kA[lambda]);
128 0 : libgrpp_evaluate_real_spherical_harmonics_array(lambda, kB_vec,
129 0 : rsh_values_kB[lambda]);
130 : }
131 :
132 : /*
133 : * pre-compute radial integrals
134 : */
135 0 : radial_type2_table_t *radial_table = libgrpp_tabulate_radial_type2_integrals(
136 : lambda1_max, lambda2_max, N_max, CA_2, CB_2, potential, shell_A, shell_B);
137 :
138 : /*
139 : * loop over shell pairs
140 : */
141 0 : for (int icart = 0; icart < size_A; icart++) {
142 0 : for (int jcart = 0; jcart < size_B; jcart++) {
143 :
144 0 : double SO_x = 0.0;
145 0 : double SO_y = 0.0;
146 0 : double SO_z = 0.0;
147 :
148 0 : int n_A = shell_A->cart_list[3 * icart + 0];
149 0 : int l_A = shell_A->cart_list[3 * icart + 1];
150 0 : int m_A = shell_A->cart_list[3 * icart + 2];
151 0 : int n_B = shell_B->cart_list[3 * jcart + 0];
152 0 : int l_B = shell_B->cart_list[3 * jcart + 1];
153 0 : int m_B = shell_B->cart_list[3 * jcart + 2];
154 :
155 0 : for (int a = 0; a <= n_A; a++) {
156 :
157 0 : double C_nA_a = libgrpp_binomial(n_A, a);
158 0 : double pow_CA_x = pow(CA_x, n_A - a);
159 :
160 0 : for (int b = 0; b <= l_A; b++) {
161 :
162 0 : double C_lA_b = libgrpp_binomial(l_A, b);
163 0 : double pow_CA_y = pow(CA_y, l_A - b);
164 :
165 0 : for (int c = 0; c <= m_A; c++) {
166 :
167 0 : double C_mA_c = libgrpp_binomial(m_A, c);
168 0 : double pow_CA_z = pow(CA_z, m_A - c);
169 :
170 0 : for (int d = 0; d <= n_B; d++) {
171 :
172 0 : double C_nB_d = libgrpp_binomial(n_B, d);
173 0 : double pow_CB_x = pow(CB_x, n_B - d);
174 :
175 0 : for (int e = 0; e <= l_B; e++) {
176 :
177 0 : double C_lB_e = libgrpp_binomial(l_B, e);
178 0 : double pow_CB_y = pow(CB_y, l_B - e);
179 :
180 0 : for (int f = 0; f <= m_B; f++) {
181 :
182 0 : double C_mB_f = libgrpp_binomial(m_B, f);
183 0 : double pow_CB_z = pow(CB_z, m_B - f);
184 :
185 0 : int N = a + b + c + d + e + f;
186 0 : double factor = C_nA_a * C_lA_b * C_mA_c * C_nB_d * C_lB_e *
187 0 : C_mB_f * pow_CA_x * pow_CA_y * pow_CA_z *
188 0 : pow_CB_x * pow_CB_y * pow_CB_z;
189 :
190 0 : if (fabs(factor) < LIBGRPP_ZERO_THRESH) {
191 0 : continue;
192 : }
193 :
194 : /*
195 : * contraction of radial integrals with angular integrals
196 : */
197 :
198 0 : double sum_omega_Q_x = 0.0;
199 0 : double sum_omega_Q_y = 0.0;
200 0 : double sum_omega_Q_z = 0.0;
201 :
202 0 : int lambda1_lower = int_max2(L - a - b - c, 0);
203 0 : int lambda2_lower = int_max2(L - d - e - f, 0);
204 0 : int lambda1_upper = L + a + b + c;
205 0 : int lambda2_upper = L + d + e + f;
206 :
207 0 : for (int lambda_1 = lambda1_lower; lambda_1 <= lambda1_upper;
208 0 : lambda_1++) {
209 0 : if ((L + a + b + c - lambda_1) % 2 != 0) {
210 0 : continue;
211 : }
212 :
213 : for (int lambda_2 = lambda2_lower;
214 0 : lambda_2 <= lambda2_upper; lambda_2++) {
215 0 : if ((L + d + e + f - lambda_2) % 2 != 0) {
216 0 : continue;
217 : }
218 :
219 0 : double QN = libgrpp_get_radial_type2_integral(
220 : radial_table, lambda_1, lambda_2, N);
221 0 : if (fabs(QN) < LIBGRPP_ZERO_THRESH) {
222 0 : continue;
223 : }
224 :
225 0 : double sum_angular_x, sum_angular_y, sum_angular_z;
226 0 : type3_angular_sum(
227 : L, Lx_matrix, Ly_matrix, Lz_matrix, lambda_1, a, b, c,
228 0 : rsh_values_kA[lambda_1], lambda_2, d, e, f,
229 0 : rsh_values_kB[lambda_2], &sum_angular_x,
230 : &sum_angular_y, &sum_angular_z);
231 :
232 0 : sum_omega_Q_x += QN * sum_angular_x;
233 0 : sum_omega_Q_y += QN * sum_angular_y;
234 0 : sum_omega_Q_z += QN * sum_angular_z;
235 : }
236 : }
237 :
238 0 : SO_x += factor * sum_omega_Q_x;
239 0 : SO_y += factor * sum_omega_Q_y;
240 0 : SO_z += factor * sum_omega_Q_z;
241 : }
242 : }
243 : }
244 : }
245 : }
246 : }
247 :
248 0 : so_x_matrix[icart * size_B + jcart] = SO_x * (16.0 * M_PI * M_PI);
249 0 : so_y_matrix[icart * size_B + jcart] = SO_y * (16.0 * M_PI * M_PI);
250 0 : so_z_matrix[icart * size_B + jcart] = SO_z * (16.0 * M_PI * M_PI);
251 : }
252 : }
253 :
254 0 : libgrpp_delete_radial_type2_integrals(radial_table);
255 0 : free(Lx_matrix);
256 0 : free(Ly_matrix);
257 0 : free(Lz_matrix);
258 0 : }
259 :
260 : /*
261 : * Double sum of products of type 2 angular integrals
262 : * (Pitzer, Winter, 1991, formula on the top of the page 776)
263 : */
264 0 : static void type3_angular_sum(int L, double *Lx_matrix, double *Ly_matrix,
265 : double *Lz_matrix, int lambda_1, int a, int b,
266 : int c, double *rsh_values_kA, int lambda_2, int d,
267 : int e, int f, double *rsh_values_kB,
268 : double *sum_angular_x, double *sum_angular_y,
269 : double *sum_angular_z) {
270 0 : *sum_angular_x = 0.0;
271 0 : *sum_angular_y = 0.0;
272 0 : *sum_angular_z = 0.0;
273 :
274 : /*
275 : * contract tensors with angular integrals
276 : */
277 0 : for (int m1 = -L; m1 <= L; m1++) {
278 0 : for (int m2 = -L; m2 <= L; m2++) {
279 :
280 0 : double lx = Lx_matrix[(2 * L + 1) * (m1 + L) + (m2 + L)];
281 0 : double ly = Ly_matrix[(2 * L + 1) * (m1 + L) + (m2 + L)];
282 0 : double lz = Lz_matrix[(2 * L + 1) * (m1 + L) + (m2 + L)];
283 0 : if (fabs(lx) < LIBGRPP_ZERO_THRESH && fabs(ly) < LIBGRPP_ZERO_THRESH &&
284 0 : fabs(lz) < LIBGRPP_ZERO_THRESH) {
285 0 : continue;
286 : }
287 :
288 0 : double omega_1 = libgrpp_angular_type2_integral(lambda_1, L, m1, a, b, c,
289 : rsh_values_kA);
290 0 : if (fabs(omega_1) < LIBGRPP_ZERO_THRESH) {
291 0 : continue;
292 : }
293 :
294 0 : double omega_2 = libgrpp_angular_type2_integral(lambda_2, L, m2, d, e, f,
295 : rsh_values_kB);
296 :
297 0 : *sum_angular_x += omega_1 * omega_2 * lx;
298 0 : *sum_angular_y += omega_1 * omega_2 * ly;
299 0 : *sum_angular_z += omega_1 * omega_2 * lz;
300 : }
301 : }
302 0 : }
|
