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 26944917 : 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 26944917 : double sum = 0.0;
78 :
79 26944917 : rsh_coef_table_t *rsh_coef_lambda =
80 26944917 : libgrpp_get_real_spherical_harmonic_table(lambda);
81 26944917 : rsh_coef_table_t *rsh_coef_L = libgrpp_get_real_spherical_harmonic_table(L);
82 :
83 26944917 : int ncomb_rst = rsh_coef_lambda->n_cart_comb;
84 26944917 : int ncomb_uvw = rsh_coef_L->n_cart_comb;
85 :
86 170896556 : for (int mu = -lambda; mu <= lambda; mu++) {
87 143951639 : double sum2 = 0.0;
88 143951639 : double rsh_value_k = rsh_values[mu + lambda];
89 143951639 : if (fabs(rsh_value_k) < LIBGRPP_ZERO_THRESH) {
90 90255984 : continue;
91 : }
92 :
93 466290920 : for (int icomb_uvw = 0; icomb_uvw < ncomb_uvw; icomb_uvw++) {
94 :
95 412595265 : int u = rsh_coef_L->cartesian_comb[3 * icomb_uvw];
96 412595265 : int v = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 1];
97 412595265 : int w = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 2];
98 412595265 : double y_lm_uvw = rsh_coef_L->coeffs[(m + L) * ncomb_uvw + icomb_uvw];
99 412595265 : if (fabs(y_lm_uvw) < LIBGRPP_ZERO_THRESH) {
100 309223129 : continue;
101 : }
102 :
103 1230534956 : for (int icomb_rst = 0; icomb_rst < ncomb_rst; icomb_rst++) {
104 :
105 1127162820 : int r = rsh_coef_lambda->cartesian_comb[3 * icomb_rst];
106 1127162820 : int s = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 1];
107 1127162820 : int t = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 2];
108 1127162820 : double y_lam_mu_rst =
109 1127162820 : rsh_coef_lambda->coeffs[(mu + lambda) * ncomb_rst + icomb_rst];
110 1127162820 : if (fabs(y_lam_mu_rst) < LIBGRPP_ZERO_THRESH) {
111 811701207 : continue;
112 : }
113 :
114 630923226 : double usp_int = integrate_unitary_sphere_polynomial(
115 315461613 : a + r + u, b + s + v, c + t + w);
116 :
117 315461613 : sum2 += y_lam_mu_rst * y_lm_uvw * usp_int;
118 : }
119 : }
120 :
121 53695655 : sum += sum2 * rsh_value_k;
122 : }
123 :
124 26944917 : return sum;
125 : }
126 :
127 : /**
128 : * Integral of the unitary sphere polynomial over full solid angle.
129 : * (see MrMurchie & Davidson, formula (30))
130 : */
131 315461613 : static double integrate_unitary_sphere_polynomial(int i, int j, int k) {
132 315461613 : if ((i % 2 == 0) && (j % 2 == 0) && (k % 2 == 0)) {
133 102591450 : double dfac_i = (double)libgrpp_double_factorial(i - 1);
134 102591450 : double dfac_j = (double)libgrpp_double_factorial(j - 1);
135 102591450 : double dfac_k = (double)libgrpp_double_factorial(k - 1);
136 102591450 : double dfac_ijk = (double)libgrpp_double_factorial(i + j + k + 1);
137 102591450 : return 4 * M_PI * dfac_i * dfac_j * dfac_k / dfac_ijk;
138 : } else {
139 : return 0.0;
140 : }
141 : }
|
