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      149632 : 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      149632 :   assert(libgrpp_is_initialized());
      43             : 
      44      149632 :   int size_A = libgrpp_get_shell_size(shell_A);
      45      149632 :   int size_B = libgrpp_get_shell_size(shell_B);
      46             : 
      47      149632 :   memset(matrix, 0, size_A * size_B * sizeof(double));
      48             : 
      49      149632 :   int L = potential->L;
      50      149632 :   int L_A =
      51      149632 :       shell_A->cart_list[0] + shell_A->cart_list[1] + shell_A->cart_list[2];
      52      149632 :   int L_B =
      53      149632 :       shell_B->cart_list[0] + shell_B->cart_list[1] + shell_B->cart_list[2];
      54             : 
      55      149632 :   double *A = shell_A->origin;
      56      149632 :   double *B = shell_B->origin;
      57      149632 :   double *C = rpp_origin;
      58             : 
      59      149632 :   double CA_x = C[0] - A[0];
      60      149632 :   double CA_y = C[1] - A[1];
      61      149632 :   double CA_z = C[2] - A[2];
      62      149632 :   double CB_x = C[0] - B[0];
      63      149632 :   double CB_y = C[1] - B[1];
      64      149632 :   double CB_z = C[2] - B[2];
      65      149632 :   double CA_2 = CA_x * CA_x + CA_y * CA_y + CA_z * CA_z;
      66      149632 :   double CB_2 = CB_x * CB_x + CB_y * CB_y + CB_z * CB_z;
      67             : 
      68      149632 :   double alpha_A = shell_A->alpha[0];
      69      149632 :   double alpha_B = shell_B->alpha[0];
      70      149632 :   double kA_x = -2.0 * (alpha_A * CA_x);
      71      149632 :   double kA_y = -2.0 * (alpha_A * CA_y);
      72      149632 :   double kA_z = -2.0 * (alpha_A * CA_z);
      73      149632 :   double kB_x = -2.0 * (alpha_B * CB_x);
      74      149632 :   double kB_y = -2.0 * (alpha_B * CB_y);
      75      149632 :   double kB_z = -2.0 * (alpha_B * CB_z);
      76      149632 :   double kA_vec[3];
      77      149632 :   kA_vec[0] = kA_x;
      78      149632 :   kA_vec[1] = kA_y;
      79      149632 :   kA_vec[2] = kA_z;
      80      149632 :   double kB_vec[3];
      81      149632 :   kB_vec[0] = kB_x;
      82      149632 :   kB_vec[1] = kB_y;
      83      149632 :   kB_vec[2] = kB_z;
      84             : 
      85      149632 :   int lambda1_max = L + L_A;
      86      149632 :   int lambda2_max = L + L_B;
      87      149632 :   int N_max = L_A + L_B;
      88             : 
      89             :   /*
      90             :    * for further evaluation of angular integrals
      91             :    */
      92      149632 :   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      149632 :   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      149632 :   double rsh_values_kA[3 * LMAX][6 * LMAX];
     105      149632 :   double rsh_values_kB[3 * LMAX][6 * LMAX];
     106             : 
     107      790347 :   for (int lambda = 0; lambda <= lmax; lambda++) {
     108      640715 :     libgrpp_evaluate_real_spherical_harmonics_array(lambda, kA_vec,
     109      640715 :                                                     rsh_values_kA[lambda]);
     110      640715 :     libgrpp_evaluate_real_spherical_harmonics_array(lambda, kB_vec,
     111      640715 :                                                     rsh_values_kB[lambda]);
     112             :   }
     113             : 
     114             :   /*
     115             :    * main loop
     116             :    * over shell pairs
     117             :    */
     118      624684 :   for (int icart = 0; icart < size_A; icart++) {
     119     2002687 :     for (int jcart = 0; jcart < size_B; jcart++) {
     120             : 
     121     1527635 :       double gamma_AB = 0.0;
     122             : 
     123     1527635 :       int n_A = shell_A->cart_list[3 * icart + 0];
     124     1527635 :       int l_A = shell_A->cart_list[3 * icart + 1];
     125     1527635 :       int m_A = shell_A->cart_list[3 * icart + 2];
     126     1527635 :       int n_B = shell_B->cart_list[3 * jcart + 0];
     127     1527635 :       int l_B = shell_B->cart_list[3 * jcart + 1];
     128     1527635 :       int m_B = shell_B->cart_list[3 * jcart + 2];
     129             : 
     130     3823736 :       for (int a = 0; a <= n_A; a++) {
     131             : 
     132     2296101 :         double C_nA_a = libgrpp_binomial(n_A, a);
     133     2296101 :         double pow_CA_x = pow(CA_x, n_A - a);
     134             : 
     135     5543335 :         for (int b = 0; b <= l_A; b++) {
     136             : 
     137     3247234 :           double C_lA_b = libgrpp_binomial(l_A, b);
     138     3247234 :           double pow_CA_y = pow(CA_y, l_A - b);
     139             : 
     140     7643851 :           for (int c = 0; c <= m_A; c++) {
     141             : 
     142     4396617 :             double C_mA_c = libgrpp_binomial(m_A, c);
     143     4396617 :             double pow_CA_z = pow(CA_z, m_A - c);
     144             : 
     145    11034695 :             for (int d = 0; d <= n_B; d++) {
     146             : 
     147     6638078 :               double C_nB_d = libgrpp_binomial(n_B, d);
     148     6638078 :               double pow_CB_x = pow(CB_x, n_B - d);
     149             : 
     150    16061097 :               for (int e = 0; e <= l_B; e++) {
     151             : 
     152     9423019 :                 double C_lB_e = libgrpp_binomial(l_B, e);
     153     9423019 :                 double pow_CB_y = pow(CB_y, l_B - e);
     154             : 
     155    22222091 :                 for (int f = 0; f <= m_B; f++) {
     156             : 
     157    12799072 :                   double C_mB_f = libgrpp_binomial(m_B, f);
     158    12799072 :                   double pow_CB_z = pow(CB_z, m_B - f);
     159             : 
     160    12799072 :                   double factor = C_nA_a * C_lA_b * C_mA_c * C_nB_d * C_lB_e *
     161    12799072 :                                   C_mB_f * pow_CA_x * pow_CA_y * pow_CA_z *
     162    12799072 :                                   pow_CB_x * pow_CB_y * pow_CB_z;
     163             : 
     164    12799072 :                   if (fabs(factor) < 1e-13) {
     165     9037990 :                     continue;
     166             :                   }
     167             : 
     168     3761082 :                   int N = a + b + c + d + e + f;
     169     3761082 :                   double sum_omega_Q = 0.0;
     170             : 
     171     3761082 :                   int lambda1_lower = int_max2(L - a - b - c, 0);
     172     3761082 :                   int lambda2_lower = int_max2(L - d - e - f, 0);
     173     3761082 :                   int lambda1_upper = L + a + b + c;
     174     3761082 :                   int lambda2_upper = L + d + e + f;
     175             : 
     176    15108023 :                   for (int lambda_1 = lambda1_lower; lambda_1 <= lambda1_upper;
     177    11346941 :                        lambda_1++) {
     178    11346941 :                     if ((L + a + b + c - lambda_1) % 2 != 0) {
     179     4112712 :                       continue;
     180             :                     }
     181             : 
     182             :                     for (int lambda_2 = lambda2_lower;
     183    29756695 :                          lambda_2 <= lambda2_upper; lambda_2++) {
     184    22522466 :                       if ((L + d + e + f - lambda_2) % 2 != 0) {
     185     8151692 :                         continue;
     186             :                       }
     187             : 
     188    14370774 :                       double QN = libgrpp_get_radial_type2_integral(
     189             :                           radial_table, lambda_1, lambda_2, N);
     190    14370774 :                       if (fabs(QN) < 1e-16) {
     191    10249936 :                         continue;
     192             :                       }
     193             : 
     194     8241676 :                       double sum_angular = type2_angular_sum(
     195             :                           L, lambda_1, a, b, c, lambda_2, d, e, f,
     196     4120838 :                           rsh_values_kA[lambda_1], rsh_values_kB[lambda_2]);
     197             : 
     198     4120838 :                       sum_omega_Q += QN * sum_angular;
     199             :                     } // loop over lambda_2
     200             :                   } // loop over lambda_1
     201             : 
     202     3761082 :                   gamma_AB += factor * sum_omega_Q;
     203             :                 }
     204             :               }
     205             :             }
     206             :           }
     207             :         }
     208             :       }
     209             : 
     210     1527635 :       gamma_AB *= 16 * M_PI * M_PI;
     211             : 
     212     1527635 :       matrix[icart * size_B + jcart] = gamma_AB;
     213             :     }
     214             :   }
     215             : 
     216      149632 :   libgrpp_delete_radial_type2_integrals(radial_table);
     217      149632 : }
     218             : 
     219             : /*
     220             :  * Sum of products of type 2 angular integrals
     221             :  * (McMurchie, Davidson, 1981, formulas (23) and (24))
     222             :  */
     223     4120838 : 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     4120838 :   double sum_angular = 0.0;
     227             : 
     228             :   /*
     229             :    * contract tensors with angular integrals
     230             :    */
     231    22874458 :   for (int m = -L; m <= L; m++) {
     232    18753620 :     double omega_1 =
     233    18753620 :         libgrpp_angular_type2_integral(lambda_1, L, m, a, b, c, rsh_values_kA);
     234    18753620 :     if (fabs(omega_1) < 1e-16) {
     235    12855324 :       continue;
     236             :     }
     237             : 
     238     5898296 :     double omega_2 =
     239     5898296 :         libgrpp_angular_type2_integral(lambda_2, L, m, d, e, f, rsh_values_kB);
     240             : 
     241     5898296 :     sum_angular += omega_1 * omega_2;
     242             :   }
     243             : 
     244     4120838 :   return sum_angular;
     245             : }