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          14 : void libgrpp_create_real_spherical_harmonic_coeffs_tables(int Lmax) {
      61          14 :   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          14 :     rsh_coef_tables = (rsh_coef_table_t **)realloc(
      66          14 :         rsh_coef_tables, (Lmax + 1) * sizeof(rsh_coef_table_t *));
      67             : 
      68         588 :     for (int L = rsh_tables_lmax + 1; L <= Lmax; L++) {
      69         574 :       rsh_coef_tables[L] = libgrpp_tabulate_real_spherical_harmonic_coeffs(L);
      70             :     }
      71          14 :     rsh_tables_lmax = Lmax;
      72             :   }
      73          14 : }
      74             : 
      75             : /**
      76             :  * Calculates all C_{l,m}^{lx,ly,lz} coefficients for the given angular momentum
      77             :  * L.
      78             :  */
      79         574 : rsh_coef_table_t *libgrpp_tabulate_real_spherical_harmonic_coeffs(int L) {
      80         574 :   int ncart = (L + 1) * (L + 2) / 2;
      81             : 
      82         574 :   rsh_coef_table_t *coef_table =
      83         574 :       (rsh_coef_table_t *)calloc(1, sizeof(rsh_coef_table_t));
      84         574 :   coef_table->n_cart_comb = ncart;
      85         574 :   coef_table->cartesian_comb = generate_cartesian_combinations(L, &ncart);
      86         574 :   coef_table->coeffs = (double *)calloc((2 * L + 1) * ncart, sizeof(double));
      87             : 
      88       24108 :   for (int m = -L; m <= L; m++) {
      89    10562748 :     for (int icomb = 0; icomb < ncart; icomb++) {
      90    10539214 :       int lx = coef_table->cartesian_comb[3 * icomb];
      91    10539214 :       int ly = coef_table->cartesian_comb[3 * icomb + 1];
      92             :       // int lz = coef_table->cartesian_comb[3 * icomb + 2];
      93    10539214 :       double u_lm_lx_ly_lz = libgrpp_spherical_to_cartesian_coef(L, m, lx, ly);
      94    10539214 :       int index = (m + L) * ncart + icomb;
      95    10539214 :       coef_table->coeffs[index] = u_lm_lx_ly_lz;
      96             :     }
      97             :   }
      98             : 
      99         574 :   return coef_table;
     100             : }
     101             : 
     102             : /**
     103             :  * Access to the table for the angular momentum value L.
     104             :  */
     105    50585262 : rsh_coef_table_t *libgrpp_get_real_spherical_harmonic_table(int L) {
     106    50585262 :   if (L > rsh_tables_lmax) {
     107           0 :     printf("get_real_spherical_harmonic_table(): %d > Lmax\n", L);
     108           0 :     return NULL;
     109             :   }
     110             : 
     111    50585262 :   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    10539214 : double libgrpp_spherical_to_cartesian_coef(int l, int m, int lx, int ly) {
     125    10539214 :   int j = lx + ly - abs(m);
     126    10539214 :   if (j % 2 != 0) {
     127             :     return 0.0;
     128             :   }
     129     5272694 :   j /= 2;
     130             : 
     131     5272694 :   if (!((m > 0 && (abs(m) - lx) % 2 == 0) || (m == 0 && lx % 2 == 0) ||
     132     2593080 :         (m < 0 && (abs(m) - lx) % 2 != 0))) {
     133             :     return 0.0;
     134             :   }
     135             : 
     136     2639434 :   double prefactor =
     137     2639434 :       sqrt((2 * l + 1) / (2 * M_PI) * libgrpp_factorial(l - abs(m)) /
     138     2639434 :            libgrpp_factorial(l + abs(m)));
     139     2639434 :   prefactor /= pow(2, l) * libgrpp_factorial(l);
     140             : 
     141     2639434 :   double u_lm_lx_ly_lz = 0.0;
     142    18795420 :   for (int i = j; i <= (l - abs(m)) / 2; i++) {
     143             :     // any term that implies the factorial of a negative number is neglected
     144    16155986 :     if (2 * l - 2 * i < 0) {
     145             :       u_lm_lx_ly_lz = 0.0;
     146             :       break;
     147             :     }
     148    16155986 :     if (l - abs(m) - 2 * i < 0) {
     149             :       u_lm_lx_ly_lz = 0.0;
     150             :       break;
     151             :     }
     152             : 
     153    32311972 :     double factor_1 =
     154    16155986 :         libgrpp_binomial(l, i) * libgrpp_binomial(i, j) * pow(-1, i) *
     155    16155986 :         libgrpp_factorial_ratio(2 * l - 2 * i, l - abs(m) - 2 * i);
     156             : 
     157    16155986 :     double sum = 0.0;
     158    66537394 :     for (int k = 0; k <= j; k++) {
     159    50381408 :       sum += libgrpp_binomial(j, k) * libgrpp_binomial(abs(m), lx - 2 * k) *
     160    50381408 :              pow(-1, (abs(m) - lx + 2 * k) / 2);
     161             :     }
     162             : 
     163    16155986 :     u_lm_lx_ly_lz += factor_1 * sum;
     164             :   }
     165             : 
     166     2639434 :   u_lm_lx_ly_lz *= prefactor;
     167     2639434 :   if (m == 0 && (lx % 2 == 0)) {
     168       46354 :     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     1281430 : void libgrpp_evaluate_real_spherical_harmonics_array(const int l,
     231             :                                                      const double *k,
     232             :                                                      double *rsh_array) {
     233     1281430 :   double unitary_kx;
     234     1281430 :   double unitary_ky;
     235     1281430 :   double unitary_kz;
     236     1281430 :   double kx_powers[200];
     237     1281430 :   double ky_powers[200];
     238     1281430 :   double kz_powers[200];
     239             : 
     240     1281430 :   rsh_coef_table_t *rsh_coef_l = libgrpp_get_real_spherical_harmonic_table(l);
     241             : 
     242     1281430 :   double length_k = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
     243             : 
     244     1281430 :   if (length_k > LIBGRPP_ZERO_THRESH) {
     245     1008601 :     double inv_length = 1.0 / length_k;
     246     1008601 :     unitary_kx = k[0] * inv_length;
     247     1008601 :     unitary_ky = k[1] * inv_length;
     248     1008601 :     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     1281430 :   kx_powers[0] = 1.0;
     256     1281430 :   ky_powers[0] = 1.0;
     257     1281430 :   kz_powers[0] = 1.0;
     258             : 
     259     3779668 :   for (int i = 1; i <= l; i++) {
     260     2498238 :     kx_powers[i] = kx_powers[i - 1] * unitary_kx;
     261     2498238 :     ky_powers[i] = ky_powers[i - 1] * unitary_ky;
     262     2498238 :     kz_powers[i] = kz_powers[i - 1] * unitary_kz;
     263             :   }
     264             : 
     265     1281430 :   memset(rsh_array, 0, (2 * l + 1) * sizeof(double));
     266             : 
     267     1281430 :   int ncart = rsh_coef_l->n_cart_comb;
     268     1281430 :   int *rst_array = rsh_coef_l->cartesian_comb;
     269             : 
     270    10469930 :   for (int icomb = 0; icomb < ncart; icomb++) {
     271     9188500 :     int r = rst_array[3 * icomb];
     272     9188500 :     int s = rst_array[3 * icomb + 1];
     273     9188500 :     int t = rst_array[3 * icomb + 2];
     274             : 
     275     9188500 :     double k_xyz = kx_powers[r] * ky_powers[s] * kz_powers[t];
     276             : 
     277    81421352 :     for (int m = -l; m <= l; m++) {
     278    72232852 :       double y_lm_rst = rsh_coef_l->coeffs[(m + l) * ncart + icomb];
     279    72232852 :       rsh_array[m + l] += y_lm_rst * k_xyz;
     280             :     }
     281             :   }
     282     1281430 : }
     283             : 
     284         574 : static int *generate_cartesian_combinations(int L, int *num) {
     285         574 :   *num = (L + 1) * (L + 2) / 2;
     286             : 
     287         574 :   int *combinations = (int *)calloc(*num, 3 * sizeof(int));
     288             : 
     289         574 :   int n = 0;
     290       12628 :   for (int i = 0; i <= L; i++) {
     291      345548 :     for (int j = 0; j <= L; j++) {
     292    10711988 :       for (int k = 0; k <= L; k++) {
     293    10378494 :         if (i + j + k == L) {
     294      172774 :           combinations[3 * n + 0] = i;
     295      172774 :           combinations[3 * n + 1] = j;
     296      172774 :           combinations[3 * n + 2] = k;
     297      172774 :           n++;
     298             :         }
     299             :       }
     300             :     }
     301             :   }
     302             : 
     303         574 :   return combinations;
     304             : }