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             :  * This file contains subroutines for evaluation of angular integrals of the 1st
      17             :  * and 2nd type. It also contains functions for construction of matrices of the
      18             :  * angular momentum operator in the bases of either real or complex spherical
      19             :  * harmonics.
      20             :  *
      21             :  * For more details on the angular integrals used in RPP integration, see:
      22             :  *
      23             :  * L. E. McMurchie, E. R. Davidson. Calculation of integrals over ab initio
      24             :  * pseudopotentials. J. Comput. Phys. 44(2), 289 (1981).
      25             :  * doi: 10.1016/0021-9991(81)90053-x
      26             :  */
      27             : 
      28             : #include <math.h>
      29             : 
      30             : #ifndef M_PI
      31             : #define M_PI 3.14159265358979323846
      32             : #endif
      33             : 
      34             : #include "grpp_angular_integrals.h"
      35             : #include "grpp_factorial.h"
      36             : #include "grpp_spherical_harmonics.h"
      37             : #include "libgrpp.h"
      38             : 
      39             : static double integrate_unitary_sphere_polynomial(int i, int j, int k);
      40             : 
      41             : /**
      42             :  * Type 1 angular integral.
      43             :  * (see MrMurchie & Davidson, formula (28))
      44             :  */
      45           0 : double libgrpp_angular_type1_integral(int lambda, int II, int JJ, int KK,
      46             :                                       double *k) {
      47           0 :   double sum = 0.0;
      48             : 
      49           0 :   for (int mu = -lambda; mu <= lambda; mu++) {
      50             :     double sum2 = 0.0;
      51           0 :     for (int r = 0; r <= lambda; r++) {
      52           0 :       for (int s = 0; s <= lambda; s++) {
      53           0 :         for (int t = 0; t <= lambda; t++) {
      54           0 :           if (r + s + t == lambda) {
      55           0 :             double y_lm_rst =
      56           0 :                 libgrpp_spherical_to_cartesian_coef(lambda, mu, r, s);
      57           0 :             double usp_int =
      58           0 :                 integrate_unitary_sphere_polynomial(II + r, JJ + s, KK + t);
      59           0 :             sum2 += y_lm_rst * usp_int;
      60             :           }
      61             :         }
      62             :       }
      63             :     }
      64           0 :     sum += sum2 * libgrpp_evaluate_real_spherical_harmonic(lambda, mu, k);
      65             :   }
      66             : 
      67           0 :   return sum;
      68             : }
      69             : 
      70             : /**
      71             :  * Type 2 angular integral.
      72             :  * (see MrMurchie & Davidson, formula (29))
      73             :  */
      74    24651916 : double libgrpp_angular_type2_integral(const int lambda, const int L,
      75             :                                       const int m, const int a, const int b,
      76             :                                       const int c, const double *rsh_values) {
      77    24651916 :   double sum = 0.0;
      78             : 
      79    24651916 :   rsh_coef_table_t *rsh_coef_lambda =
      80    24651916 :       libgrpp_get_real_spherical_harmonic_table(lambda);
      81    24651916 :   rsh_coef_table_t *rsh_coef_L = libgrpp_get_real_spherical_harmonic_table(L);
      82             : 
      83    24651916 :   int ncomb_rst = rsh_coef_lambda->n_cart_comb;
      84    24651916 :   int ncomb_uvw = rsh_coef_L->n_cart_comb;
      85             : 
      86   161031732 :   for (int mu = -lambda; mu <= lambda; mu++) {
      87   136379816 :     double sum2 = 0.0;
      88   136379816 :     double rsh_value_k = rsh_values[mu + lambda];
      89   136379816 :     if (fabs(rsh_value_k) < LIBGRPP_ZERO_THRESH) {
      90    90245036 :       continue;
      91             :     }
      92             : 
      93   421180984 :     for (int icomb_uvw = 0; icomb_uvw < ncomb_uvw; icomb_uvw++) {
      94             : 
      95   375046204 :       int u = rsh_coef_L->cartesian_comb[3 * icomb_uvw];
      96   375046204 :       int v = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 1];
      97   375046204 :       int w = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 2];
      98   375046204 :       double y_lm_uvw = rsh_coef_L->coeffs[(m + L) * ncomb_uvw + icomb_uvw];
      99   375046204 :       if (fabs(y_lm_uvw) < LIBGRPP_ZERO_THRESH) {
     100   282446707 :         continue;
     101             :       }
     102             : 
     103  1145537575 :       for (int icomb_rst = 0; icomb_rst < ncomb_rst; icomb_rst++) {
     104             : 
     105  1052938078 :         int r = rsh_coef_lambda->cartesian_comb[3 * icomb_rst];
     106  1052938078 :         int s = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 1];
     107  1052938078 :         int t = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 2];
     108  1052938078 :         double y_lam_mu_rst =
     109  1052938078 :             rsh_coef_lambda->coeffs[(mu + lambda) * ncomb_rst + icomb_rst];
     110  1052938078 :         if (fabs(y_lam_mu_rst) < LIBGRPP_ZERO_THRESH) {
     111   756465567 :           continue;
     112             :         }
     113             : 
     114   592945022 :         double usp_int = integrate_unitary_sphere_polynomial(
     115   296472511 :             a + r + u, b + s + v, c + t + w);
     116             : 
     117   296472511 :         sum2 += y_lam_mu_rst * y_lm_uvw * usp_int;
     118             :       }
     119             :     }
     120             : 
     121    46134780 :     sum += sum2 * rsh_value_k;
     122             :   }
     123             : 
     124    24651916 :   return sum;
     125             : }
     126             : 
     127             : /**
     128             :  * Integral of the unitary sphere polynomial over full solid angle.
     129             :  * (see MrMurchie & Davidson, formula (30))
     130             :  */
     131   296472511 : static double integrate_unitary_sphere_polynomial(int i, int j, int k) {
     132   296472511 :   if ((i % 2 == 0) && (j % 2 == 0) && (k % 2 == 0)) {
     133    96429012 :     double dfac_i = (double)libgrpp_double_factorial(i - 1);
     134    96429012 :     double dfac_j = (double)libgrpp_double_factorial(j - 1);
     135    96429012 :     double dfac_k = (double)libgrpp_double_factorial(k - 1);
     136    96429012 :     double dfac_ijk = (double)libgrpp_double_factorial(i + j + k + 1);
     137    96429012 :     return 4 * M_PI * dfac_i * dfac_j * dfac_k / dfac_ijk;
     138             :   } else {
     139             :     return 0.0;
     140             :   }
     141             : }