Loading [MathJax]/extensions/tex2jax.js
 (git:aabdcc8)
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
grpp_spherical_harmonics.c
Go to the documentation of this file.
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"
40#include "libgrpp.h"
41/*
42 * Tables with pretabulated expansion coefficients
43 */
44
46static int rsh_tables_lmax = -1;
47
48/*
49 * Function pre-definitions
50 */
52
53static 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 */
61 if (Lmax <= rsh_tables_lmax) {
62 // nothing to do
63 } else {
64 // expand tables: realloc memory and add tables for the highest L values
66 rsh_coef_tables, (Lmax + 1) * sizeof(rsh_coef_table_t *));
67
68 for (int L = rsh_tables_lmax + 1; L <= Lmax; L++) {
70 }
71 rsh_tables_lmax = Lmax;
72 }
73}
74
75/**
76 * Calculates all C_{l,m}^{lx,ly,lz} coefficients for the given angular momentum
77 * L.
78 */
80 int ncart = (L + 1) * (L + 2) / 2;
81
82 rsh_coef_table_t *coef_table =
83 (rsh_coef_table_t *)calloc(1, sizeof(rsh_coef_table_t));
84 coef_table->n_cart_comb = ncart;
85 coef_table->cartesian_comb = generate_cartesian_combinations(L, &ncart);
86 coef_table->coeffs = (double *)calloc((2 * L + 1) * ncart, sizeof(double));
87
88 for (int m = -L; m <= L; m++) {
89 for (int icomb = 0; icomb < ncart; icomb++) {
90 int lx = coef_table->cartesian_comb[3 * icomb];
91 int ly = coef_table->cartesian_comb[3 * icomb + 1];
92 // int lz = coef_table->cartesian_comb[3 * icomb + 2];
93 double u_lm_lx_ly_lz = libgrpp_spherical_to_cartesian_coef(L, m, lx, ly);
94 int index = (m + L) * ncart + icomb;
95 coef_table->coeffs[index] = u_lm_lx_ly_lz;
96 }
97 }
98
99 return coef_table;
100}
101
102/**
103 * Access to the table for the angular momentum value L.
104 */
106 if (L > rsh_tables_lmax) {
107 printf("get_real_spherical_harmonic_table(): %d > Lmax\n", L);
108 return NULL;
109 }
110
111 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 */
124double libgrpp_spherical_to_cartesian_coef(int l, int m, int lx, int ly) {
125 int j = lx + ly - abs(m);
126 if (j % 2 != 0) {
127 return 0.0;
128 }
129 j /= 2;
130
131 if (!((m > 0 && (abs(m) - lx) % 2 == 0) || (m == 0 && lx % 2 == 0) ||
132 (m < 0 && (abs(m) - lx) % 2 != 0))) {
133 return 0.0;
134 }
135
136 double prefactor =
137 sqrt((2 * l + 1) / (2 * M_PI) * libgrpp_factorial(l - abs(m)) /
138 libgrpp_factorial(l + abs(m)));
139 prefactor /= pow(2, l) * libgrpp_factorial(l);
140
141 double u_lm_lx_ly_lz = 0.0;
142 for (int i = j; i <= (l - abs(m)) / 2; i++) {
143 // any term that implies the factorial of a negative number is neglected
144 if (2 * l - 2 * i < 0) {
145 u_lm_lx_ly_lz = 0.0;
146 break;
147 }
148 if (l - abs(m) - 2 * i < 0) {
149 u_lm_lx_ly_lz = 0.0;
150 break;
151 }
152
153 double factor_1 =
154 libgrpp_binomial(l, i) * libgrpp_binomial(i, j) * pow(-1, i) *
155 libgrpp_factorial_ratio(2 * l - 2 * i, l - abs(m) - 2 * i);
156
157 double sum = 0.0;
158 for (int k = 0; k <= j; k++) {
159 sum += libgrpp_binomial(j, k) * libgrpp_binomial(abs(m), lx - 2 * k) *
160 pow(-1, (abs(m) - lx + 2 * k) / 2);
161 }
162
163 u_lm_lx_ly_lz += factor_1 * sum;
164 }
165
166 u_lm_lx_ly_lz *= prefactor;
167 if (m == 0 && (lx % 2 == 0)) {
168 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 */
178double libgrpp_evaluate_real_spherical_harmonic(const int l, const int m,
179 const double *k) {
180 double unitary_kx;
181 double unitary_ky;
182 double unitary_kz;
183 double kx_powers[200];
184 double ky_powers[200];
185 double kz_powers[200];
186
188
189 double length_k = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
190
191 if (length_k > LIBGRPP_ZERO_THRESH) {
192 unitary_kx = k[0] / length_k;
193 unitary_ky = k[1] / length_k;
194 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 kx_powers[0] = 1.0;
202 ky_powers[0] = 1.0;
203 kz_powers[0] = 1.0;
204
205 for (int i = 1; i <= l; i++) {
206 kx_powers[i] = kx_powers[i - 1] * unitary_kx;
207 ky_powers[i] = ky_powers[i - 1] * unitary_ky;
208 kz_powers[i] = kz_powers[i - 1] * unitary_kz;
209 }
210
211 double value = 0.0;
212
213 int ncart = rsh_coef_l->n_cart_comb;
214 for (int icomb = 0; icomb < ncart; icomb++) {
215 int r = rsh_coef_l->cartesian_comb[3 * icomb];
216 int s = rsh_coef_l->cartesian_comb[3 * icomb + 1];
217 int t = rsh_coef_l->cartesian_comb[3 * icomb + 2];
218 double y_lm_rst = rsh_coef_l->coeffs[(m + l) * ncart + icomb];
219
220 value += y_lm_rst * kx_powers[r] * ky_powers[s] * kz_powers[t];
221 }
222
223 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 */
231 const double *k,
232 double *rsh_array) {
233 double unitary_kx;
234 double unitary_ky;
235 double unitary_kz;
236 double kx_powers[200];
237 double ky_powers[200];
238 double kz_powers[200];
239
241
242 double length_k = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
243
244 if (length_k > LIBGRPP_ZERO_THRESH) {
245 double inv_length = 1.0 / length_k;
246 unitary_kx = k[0] * inv_length;
247 unitary_ky = k[1] * inv_length;
248 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 kx_powers[0] = 1.0;
256 ky_powers[0] = 1.0;
257 kz_powers[0] = 1.0;
258
259 for (int i = 1; i <= l; i++) {
260 kx_powers[i] = kx_powers[i - 1] * unitary_kx;
261 ky_powers[i] = ky_powers[i - 1] * unitary_ky;
262 kz_powers[i] = kz_powers[i - 1] * unitary_kz;
263 }
264
265 memset(rsh_array, 0, (2 * l + 1) * sizeof(double));
266
267 int ncart = rsh_coef_l->n_cart_comb;
268 int *rst_array = rsh_coef_l->cartesian_comb;
269
270 for (int icomb = 0; icomb < ncart; icomb++) {
271 int r = rst_array[3 * icomb];
272 int s = rst_array[3 * icomb + 1];
273 int t = rst_array[3 * icomb + 2];
274
275 double k_xyz = kx_powers[r] * ky_powers[s] * kz_powers[t];
276
277 for (int m = -l; m <= l; m++) {
278 double y_lm_rst = rsh_coef_l->coeffs[(m + l) * ncart + icomb];
279 rsh_array[m + l] += y_lm_rst * k_xyz;
280 }
281 }
282}
283
284static int *generate_cartesian_combinations(int L, int *num) {
285 *num = (L + 1) * (L + 2) / 2;
286
287 int *combinations = (int *)calloc(*num, 3 * sizeof(int));
288
289 int n = 0;
290 for (int i = 0; i <= L; i++) {
291 for (int j = 0; j <= L; j++) {
292 for (int k = 0; k <= L; k++) {
293 if (i + j + k == L) {
294 combinations[3 * n + 0] = i;
295 combinations[3 * n + 1] = j;
296 combinations[3 * n + 2] = k;
297 n++;
298 }
299 }
300 }
301 }
302
303 return combinations;
304}
static void const int const int i
uint64_t libgrpp_binomial(uint64_t n, uint64_t k)
double libgrpp_factorial(int n)
double libgrpp_factorial_ratio(int n, int m)
rsh_coef_table_t * libgrpp_get_real_spherical_harmonic_table(int L)
void libgrpp_evaluate_real_spherical_harmonics_array(const int l, const double *k, double *rsh_array)
rsh_coef_table_t * libgrpp_tabulate_real_spherical_harmonic_coeffs(int L)
static int * generate_cartesian_combinations(int L, int *num)
static rsh_coef_table_t ** rsh_coef_tables
#define M_SQRT1_2
void libgrpp_create_real_spherical_harmonic_coeffs_tables(int Lmax)
#define M_PI
double libgrpp_spherical_to_cartesian_coef(int l, int m, int lx, int ly)
static int rsh_tables_lmax
double libgrpp_evaluate_real_spherical_harmonic(const int l, const int m, const double *k)
#define LIBGRPP_ZERO_THRESH