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 : * Constructs tables with the expansion coefficients of real spherical harmonics
17 : * in the basis of (cartesian) unitary spherical polynomials.
18 : *
19 : * For more details about the algorithm used, see:
20 : * R. Flores-Moreno et al, J. Comput. Chem. 27, 1009 (2006),
21 : * doi: 10.1002/jcc.20410
22 : */
23 :
24 : #include <math.h>
25 : #include <stdio.h>
26 : #include <stdlib.h>
27 : #include <string.h>
28 :
29 : #ifndef M_PI
30 : #define M_PI 3.1415926535897932384626433
31 : #endif
32 :
33 : #ifndef M_SQRT1_2
34 : #define M_SQRT1_2 0.70710678118654752440
35 : #endif
36 :
37 : #include "grpp_binomial.h"
38 : #include "grpp_factorial.h"
39 : #include "grpp_spherical_harmonics.h"
40 : #include "libgrpp.h"
41 : /*
42 : * Tables with pretabulated expansion coefficients
43 : */
44 :
45 : static rsh_coef_table_t **rsh_coef_tables = NULL;
46 : static int rsh_tables_lmax = -1;
47 :
48 : /*
49 : * Function pre-definitions
50 : */
51 : rsh_coef_table_t *libgrpp_tabulate_real_spherical_harmonic_coeffs(int L);
52 :
53 : static int *generate_cartesian_combinations(int L, int *num);
54 :
55 : /**
56 : * Constructs the set of tables with C_{l,m}^{lx,ly,lz} coefficients
57 : * (up to maximum angular momentum Lmax).
58 : * (pretabulation step)
59 : */
60 18 : void libgrpp_create_real_spherical_harmonic_coeffs_tables(int Lmax) {
61 18 : if (Lmax <= rsh_tables_lmax) {
62 : // nothing to do
63 : } else {
64 : // expand tables: realloc memory and add tables for the highest L values
65 18 : rsh_coef_tables = (rsh_coef_table_t **)realloc(
66 18 : rsh_coef_tables, (Lmax + 1) * sizeof(rsh_coef_table_t *));
67 :
68 756 : for (int L = rsh_tables_lmax + 1; L <= Lmax; L++) {
69 738 : rsh_coef_tables[L] = libgrpp_tabulate_real_spherical_harmonic_coeffs(L);
70 : }
71 18 : rsh_tables_lmax = Lmax;
72 : }
73 18 : }
74 :
75 : /**
76 : * Calculates all C_{l,m}^{lx,ly,lz} coefficients for the given angular momentum
77 : * L.
78 : */
79 738 : rsh_coef_table_t *libgrpp_tabulate_real_spherical_harmonic_coeffs(int L) {
80 738 : int ncart = (L + 1) * (L + 2) / 2;
81 :
82 738 : rsh_coef_table_t *coef_table =
83 738 : (rsh_coef_table_t *)calloc(1, sizeof(rsh_coef_table_t));
84 738 : coef_table->n_cart_comb = ncart;
85 738 : coef_table->cartesian_comb = generate_cartesian_combinations(L, &ncart);
86 738 : coef_table->coeffs = (double *)calloc((2 * L + 1) * ncart, sizeof(double));
87 :
88 30996 : for (int m = -L; m <= L; m++) {
89 13580676 : for (int icomb = 0; icomb < ncart; icomb++) {
90 13550418 : int lx = coef_table->cartesian_comb[3 * icomb];
91 13550418 : int ly = coef_table->cartesian_comb[3 * icomb + 1];
92 : // int lz = coef_table->cartesian_comb[3 * icomb + 2];
93 13550418 : double u_lm_lx_ly_lz = libgrpp_spherical_to_cartesian_coef(L, m, lx, ly);
94 13550418 : int index = (m + L) * ncart + icomb;
95 13550418 : coef_table->coeffs[index] = u_lm_lx_ly_lz;
96 : }
97 : }
98 :
99 738 : return coef_table;
100 : }
101 :
102 : /**
103 : * Access to the table for the angular momentum value L.
104 : */
105 55473634 : rsh_coef_table_t *libgrpp_get_real_spherical_harmonic_table(int L) {
106 55473634 : if (L > rsh_tables_lmax) {
107 0 : printf("get_real_spherical_harmonic_table(): %d > Lmax\n", L);
108 0 : return NULL;
109 : }
110 :
111 55473634 : return rsh_coef_tables[L];
112 : }
113 :
114 : /**
115 : * For the given real spherical harmonic (RSH) S_lm calculates the coefficient
116 : * C_{l,m}^{lx,ly,lz} before the unitary spherical polynomial (USP) in its
117 : * expansion.
118 : *
119 : * The formula is taken from:
120 : * R. Flores-Moreno et al, J. Comput. Chem. 27, 1009 (2006)
121 : * doi: 10.1002/jcc.20410
122 : * (formula 32)
123 : */
124 13550418 : double libgrpp_spherical_to_cartesian_coef(int l, int m, int lx, int ly) {
125 13550418 : int j = lx + ly - abs(m);
126 13550418 : if (j % 2 != 0) {
127 : return 0.0;
128 : }
129 6779178 : j /= 2;
130 :
131 6779178 : if (!((m > 0 && (abs(m) - lx) % 2 == 0) || (m == 0 && lx % 2 == 0) ||
132 3333960 : (m < 0 && (abs(m) - lx) % 2 != 0))) {
133 : return 0.0;
134 : }
135 :
136 3393558 : double prefactor =
137 3393558 : sqrt((2 * l + 1) / (2 * M_PI) * libgrpp_factorial(l - abs(m)) /
138 3393558 : libgrpp_factorial(l + abs(m)));
139 3393558 : prefactor /= pow(2, l) * libgrpp_factorial(l);
140 :
141 3393558 : double u_lm_lx_ly_lz = 0.0;
142 24165540 : for (int i = j; i <= (l - abs(m)) / 2; i++) {
143 : // any term that implies the factorial of a negative number is neglected
144 20771982 : if (2 * l - 2 * i < 0) {
145 : u_lm_lx_ly_lz = 0.0;
146 : break;
147 : }
148 20771982 : if (l - abs(m) - 2 * i < 0) {
149 : u_lm_lx_ly_lz = 0.0;
150 : break;
151 : }
152 :
153 41543964 : double factor_1 =
154 20771982 : libgrpp_binomial(l, i) * libgrpp_binomial(i, j) * pow(-1, i) *
155 20771982 : libgrpp_factorial_ratio(2 * l - 2 * i, l - abs(m) - 2 * i);
156 :
157 20771982 : double sum = 0.0;
158 85548078 : for (int k = 0; k <= j; k++) {
159 64776096 : sum += libgrpp_binomial(j, k) * libgrpp_binomial(abs(m), lx - 2 * k) *
160 64776096 : pow(-1, (abs(m) - lx + 2 * k) / 2);
161 : }
162 :
163 20771982 : u_lm_lx_ly_lz += factor_1 * sum;
164 : }
165 :
166 3393558 : u_lm_lx_ly_lz *= prefactor;
167 3393558 : if (m == 0 && (lx % 2 == 0)) {
168 59598 : u_lm_lx_ly_lz *= M_SQRT1_2; // x 1/sqrt(2)
169 : }
170 :
171 : return u_lm_lx_ly_lz;
172 : }
173 :
174 : /**
175 : * Calculates value of the real spherical harmonic S_lm at the point k/|k| of
176 : * the unit sphere.
177 : */
178 0 : double libgrpp_evaluate_real_spherical_harmonic(const int l, const int m,
179 : const double *k) {
180 0 : double unitary_kx;
181 0 : double unitary_ky;
182 0 : double unitary_kz;
183 0 : double kx_powers[200];
184 0 : double ky_powers[200];
185 0 : double kz_powers[200];
186 :
187 0 : rsh_coef_table_t *rsh_coef_l = libgrpp_get_real_spherical_harmonic_table(l);
188 :
189 0 : double length_k = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
190 :
191 0 : if (length_k > LIBGRPP_ZERO_THRESH) {
192 0 : unitary_kx = k[0] / length_k;
193 0 : unitary_ky = k[1] / length_k;
194 0 : unitary_kz = k[2] / length_k;
195 : } else {
196 : unitary_kx = 0.0;
197 : unitary_ky = 0.0;
198 : unitary_kz = 0.0;
199 : }
200 :
201 0 : kx_powers[0] = 1.0;
202 0 : ky_powers[0] = 1.0;
203 0 : kz_powers[0] = 1.0;
204 :
205 0 : for (int i = 1; i <= l; i++) {
206 0 : kx_powers[i] = kx_powers[i - 1] * unitary_kx;
207 0 : ky_powers[i] = ky_powers[i - 1] * unitary_ky;
208 0 : kz_powers[i] = kz_powers[i - 1] * unitary_kz;
209 : }
210 :
211 0 : double value = 0.0;
212 :
213 0 : int ncart = rsh_coef_l->n_cart_comb;
214 0 : for (int icomb = 0; icomb < ncart; icomb++) {
215 0 : int r = rsh_coef_l->cartesian_comb[3 * icomb];
216 0 : int s = rsh_coef_l->cartesian_comb[3 * icomb + 1];
217 0 : int t = rsh_coef_l->cartesian_comb[3 * icomb + 2];
218 0 : double y_lm_rst = rsh_coef_l->coeffs[(m + l) * ncart + icomb];
219 :
220 0 : value += y_lm_rst * kx_powers[r] * ky_powers[s] * kz_powers[t];
221 : }
222 :
223 0 : return value;
224 : }
225 :
226 : /**
227 : * Calculates values of the real spherical harmonic S_lm at the point k/|k| of
228 : * the unit sphere for all m = -l, ..., +l
229 : */
230 1583800 : void libgrpp_evaluate_real_spherical_harmonics_array(const int l,
231 : const double *k,
232 : double *rsh_array) {
233 1583800 : double unitary_kx;
234 1583800 : double unitary_ky;
235 1583800 : double unitary_kz;
236 1583800 : double kx_powers[200];
237 1583800 : double ky_powers[200];
238 1583800 : double kz_powers[200];
239 :
240 1583800 : rsh_coef_table_t *rsh_coef_l = libgrpp_get_real_spherical_harmonic_table(l);
241 :
242 1583800 : double length_k = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
243 :
244 1583800 : if (length_k > LIBGRPP_ZERO_THRESH) {
245 1103878 : double inv_length = 1.0 / length_k;
246 1103878 : unitary_kx = k[0] * inv_length;
247 1103878 : unitary_ky = k[1] * inv_length;
248 1103878 : unitary_kz = k[2] * inv_length;
249 : } else {
250 : unitary_kx = 0.0;
251 : unitary_ky = 0.0;
252 : unitary_kz = 0.0;
253 : }
254 :
255 1583800 : kx_powers[0] = 1.0;
256 1583800 : ky_powers[0] = 1.0;
257 1583800 : kz_powers[0] = 1.0;
258 :
259 4507686 : for (int i = 1; i <= l; i++) {
260 2923886 : kx_powers[i] = kx_powers[i - 1] * unitary_kx;
261 2923886 : ky_powers[i] = ky_powers[i - 1] * unitary_ky;
262 2923886 : kz_powers[i] = kz_powers[i - 1] * unitary_kz;
263 : }
264 :
265 1583800 : memset(rsh_array, 0, (2 * l + 1) * sizeof(double));
266 :
267 1583800 : int ncart = rsh_coef_l->n_cart_comb;
268 1583800 : int *rst_array = rsh_coef_l->cartesian_comb;
269 :
270 12248382 : for (int icomb = 0; icomb < ncart; icomb++) {
271 10664582 : int r = rst_array[3 * icomb];
272 10664582 : int s = rst_array[3 * icomb + 1];
273 10664582 : int t = rst_array[3 * icomb + 2];
274 :
275 10664582 : double k_xyz = kx_powers[r] * ky_powers[s] * kz_powers[t];
276 :
277 91634704 : for (int m = -l; m <= l; m++) {
278 80970122 : double y_lm_rst = rsh_coef_l->coeffs[(m + l) * ncart + icomb];
279 80970122 : rsh_array[m + l] += y_lm_rst * k_xyz;
280 : }
281 : }
282 1583800 : }
283 :
284 738 : static int *generate_cartesian_combinations(int L, int *num) {
285 738 : *num = (L + 1) * (L + 2) / 2;
286 :
287 738 : int *combinations = (int *)calloc(*num, 3 * sizeof(int));
288 :
289 738 : int n = 0;
290 16236 : for (int i = 0; i <= L; i++) {
291 444276 : for (int j = 0; j <= L; j++) {
292 13772556 : for (int k = 0; k <= L; k++) {
293 13343778 : if (i + j + k == L) {
294 222138 : combinations[3 * n + 0] = i;
295 222138 : combinations[3 * n + 1] = j;
296 222138 : combinations[3 * n + 2] = k;
297 222138 : n++;
298 : }
299 : }
300 : }
301 : }
302 :
303 738 : return combinations;
304 : }
|
