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 <string.h>
18 : #ifndef M_PI
19 : #define M_PI 3.14159265358979323846
20 : #endif
21 :
22 : #include "grpp_angular_integrals.h"
23 : #include "grpp_binomial.h"
24 : #include "grpp_radial_type2_integral.h"
25 : #include "grpp_spherical_harmonics.h"
26 : #include "grpp_utils.h"
27 : #include "libgrpp.h"
28 :
29 : #define LMAX (2 * LIBGRPP_MAX_BASIS_L + LIBGRPP_MAX_RPP_L)
30 :
31 : static double type2_angular_sum(int L, int lambda_1, int a, int b, int c,
32 : int lambda_2, int d, int e, int f,
33 : double *rsh_values_kA, double *rsh_values_kB);
34 :
35 : /**
36 : * Evaluation of type 2 RPP integrals (scalar-relativistic semilocal RPP with
37 : * L-projectors).
38 : */
39 194176 : void libgrpp_type2_integrals(libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
40 : double *rpp_origin, libgrpp_potential_t *potential,
41 : double *matrix) {
42 194176 : assert(libgrpp_is_initialized());
43 :
44 194176 : int size_A = libgrpp_get_shell_size(shell_A);
45 194176 : int size_B = libgrpp_get_shell_size(shell_B);
46 :
47 194176 : memset(matrix, 0, size_A * size_B * sizeof(double));
48 :
49 194176 : int L = potential->L;
50 194176 : int L_A =
51 194176 : shell_A->cart_list[0] + shell_A->cart_list[1] + shell_A->cart_list[2];
52 194176 : int L_B =
53 194176 : shell_B->cart_list[0] + shell_B->cart_list[1] + shell_B->cart_list[2];
54 :
55 194176 : double *A = shell_A->origin;
56 194176 : double *B = shell_B->origin;
57 194176 : double *C = rpp_origin;
58 :
59 194176 : double CA_x = C[0] - A[0];
60 194176 : double CA_y = C[1] - A[1];
61 194176 : double CA_z = C[2] - A[2];
62 194176 : double CB_x = C[0] - B[0];
63 194176 : double CB_y = C[1] - B[1];
64 194176 : double CB_z = C[2] - B[2];
65 194176 : double CA_2 = CA_x * CA_x + CA_y * CA_y + CA_z * CA_z;
66 194176 : double CB_2 = CB_x * CB_x + CB_y * CB_y + CB_z * CB_z;
67 :
68 194176 : double alpha_A = shell_A->alpha[0];
69 194176 : double alpha_B = shell_B->alpha[0];
70 194176 : double kA_x = -2.0 * (alpha_A * CA_x);
71 194176 : double kA_y = -2.0 * (alpha_A * CA_y);
72 194176 : double kA_z = -2.0 * (alpha_A * CA_z);
73 194176 : double kB_x = -2.0 * (alpha_B * CB_x);
74 194176 : double kB_y = -2.0 * (alpha_B * CB_y);
75 194176 : double kB_z = -2.0 * (alpha_B * CB_z);
76 194176 : double kA_vec[3];
77 194176 : kA_vec[0] = kA_x;
78 194176 : kA_vec[1] = kA_y;
79 194176 : kA_vec[2] = kA_z;
80 194176 : double kB_vec[3];
81 194176 : kB_vec[0] = kB_x;
82 194176 : kB_vec[1] = kB_y;
83 194176 : kB_vec[2] = kB_z;
84 :
85 194176 : int lambda1_max = L + L_A;
86 194176 : int lambda2_max = L + L_B;
87 194176 : int N_max = L_A + L_B;
88 :
89 : /*
90 : * for further evaluation of angular integrals
91 : */
92 194176 : int lmax = int_max3(lambda1_max, lambda2_max, L);
93 : // create_real_spherical_harmonic_coeffs_tables(lmax);
94 :
95 : /*
96 : * pre-compute type 2 radial integrals
97 : */
98 194176 : radial_type2_table_t *radial_table = libgrpp_tabulate_radial_type2_integrals(
99 : lambda1_max, lambda2_max, N_max, CA_2, CB_2, potential, shell_A, shell_B);
100 :
101 : /*
102 : * pre-calculate values of real spherical harmonics for different L
103 : */
104 194176 : double rsh_values_kA[3 * LMAX][6 * LMAX];
105 194176 : double rsh_values_kB[3 * LMAX][6 * LMAX];
106 :
107 986076 : for (int lambda = 0; lambda <= lmax; lambda++) {
108 791900 : libgrpp_evaluate_real_spherical_harmonics_array(lambda, kA_vec,
109 791900 : rsh_values_kA[lambda]);
110 791900 : libgrpp_evaluate_real_spherical_harmonics_array(lambda, kB_vec,
111 791900 : rsh_values_kB[lambda]);
112 : }
113 :
114 : /*
115 : * main loop
116 : * over shell pairs
117 : */
118 809928 : for (int icart = 0; icart < size_A; icart++) {
119 2619457 : for (int jcart = 0; jcart < size_B; jcart++) {
120 :
121 2003705 : double gamma_AB = 0.0;
122 :
123 2003705 : int n_A = shell_A->cart_list[3 * icart + 0];
124 2003705 : int l_A = shell_A->cart_list[3 * icart + 1];
125 2003705 : int m_A = shell_A->cart_list[3 * icart + 2];
126 2003705 : int n_B = shell_B->cart_list[3 * jcart + 0];
127 2003705 : int l_B = shell_B->cart_list[3 * jcart + 1];
128 2003705 : int m_B = shell_B->cart_list[3 * jcart + 2];
129 :
130 5029730 : for (int a = 0; a <= n_A; a++) {
131 :
132 3026025 : double C_nA_a = libgrpp_binomial(n_A, a);
133 3026025 : double pow_CA_x = pow(CA_x, n_A - a);
134 :
135 7328497 : for (int b = 0; b <= l_A; b++) {
136 :
137 4302472 : double C_lA_b = libgrpp_binomial(l_A, b);
138 4302472 : double pow_CA_y = pow(CA_y, l_A - b);
139 :
140 10160005 : for (int c = 0; c <= m_A; c++) {
141 :
142 5857533 : double C_mA_c = libgrpp_binomial(m_A, c);
143 5857533 : double pow_CA_z = pow(CA_z, m_A - c);
144 :
145 14773946 : for (int d = 0; d <= n_B; d++) {
146 :
147 8916413 : double C_nB_d = libgrpp_binomial(n_B, d);
148 8916413 : double pow_CB_x = pow(CB_x, n_B - d);
149 :
150 21676494 : for (int e = 0; e <= l_B; e++) {
151 :
152 12760081 : double C_lB_e = libgrpp_binomial(l_B, e);
153 12760081 : double pow_CB_y = pow(CB_y, l_B - e);
154 :
155 30227492 : for (int f = 0; f <= m_B; f++) {
156 :
157 17467411 : double C_mB_f = libgrpp_binomial(m_B, f);
158 17467411 : double pow_CB_z = pow(CB_z, m_B - f);
159 :
160 17467411 : double factor = C_nA_a * C_lA_b * C_mA_c * C_nB_d * C_lB_e *
161 17467411 : C_mB_f * pow_CA_x * pow_CA_y * pow_CA_z *
162 17467411 : pow_CB_x * pow_CB_y * pow_CB_z;
163 :
164 17467411 : if (fabs(factor) < 1e-13) {
165 13020436 : continue;
166 : }
167 :
168 4446975 : int N = a + b + c + d + e + f;
169 4446975 : double sum_omega_Q = 0.0;
170 :
171 4446975 : int lambda1_lower = int_max2(L - a - b - c, 0);
172 4446975 : int lambda2_lower = int_max2(L - d - e - f, 0);
173 4446975 : int lambda1_upper = L + a + b + c;
174 4446975 : int lambda2_upper = L + d + e + f;
175 :
176 17827196 : for (int lambda_1 = lambda1_lower; lambda_1 <= lambda1_upper;
177 13380221 : lambda_1++) {
178 13380221 : if ((L + a + b + c - lambda_1) % 2 != 0) {
179 4890205 : continue;
180 : }
181 :
182 : for (int lambda_2 = lambda2_lower;
183 35120135 : lambda_2 <= lambda2_upper; lambda_2++) {
184 26630119 : if ((L + d + e + f - lambda_2) % 2 != 0) {
185 9770400 : continue;
186 : }
187 :
188 16859719 : double QN = libgrpp_get_radial_type2_integral(
189 : radial_table, lambda_1, lambda_2, N);
190 16859719 : if (fabs(QN) < 1e-16) {
191 12324324 : continue;
192 : }
193 :
194 9070790 : double sum_angular = type2_angular_sum(
195 : L, lambda_1, a, b, c, lambda_2, d, e, f,
196 4535395 : rsh_values_kA[lambda_1], rsh_values_kB[lambda_2]);
197 :
198 4535395 : sum_omega_Q += QN * sum_angular;
199 : } // loop over lambda_2
200 : } // loop over lambda_1
201 :
202 4446975 : gamma_AB += factor * sum_omega_Q;
203 : }
204 : }
205 : }
206 : }
207 : }
208 : }
209 :
210 2003705 : gamma_AB *= 16 * M_PI * M_PI;
211 :
212 2003705 : matrix[icart * size_B + jcart] = gamma_AB;
213 : }
214 : }
215 :
216 194176 : libgrpp_delete_radial_type2_integrals(radial_table);
217 194176 : }
218 :
219 : /*
220 : * Sum of products of type 2 angular integrals
221 : * (McMurchie, Davidson, 1981, formulas (23) and (24))
222 : */
223 4535395 : static double type2_angular_sum(int L, int lambda_1, int a, int b, int c,
224 : int lambda_2, int d, int e, int f,
225 : double *rsh_values_kA, double *rsh_values_kB) {
226 4535395 : double sum_angular = 0.0;
227 :
228 : /*
229 : * contract tensors with angular integrals
230 : */
231 24625074 : for (int m = -L; m <= L; m++) {
232 20089679 : double omega_1 =
233 20089679 : libgrpp_angular_type2_integral(lambda_1, L, m, a, b, c, rsh_values_kA);
234 20089679 : if (fabs(omega_1) < 1e-16) {
235 13234441 : continue;
236 : }
237 :
238 6855238 : double omega_2 =
239 6855238 : libgrpp_angular_type2_integral(lambda_2, L, m, d, e, f, rsh_values_kB);
240 :
241 6855238 : sum_angular += omega_1 * omega_2;
242 : }
243 :
244 4535395 : return sum_angular;
245 : }
|
