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 <stdio.h>
18 : #include <stdlib.h>
19 : #include <string.h>
20 :
21 : #ifndef M_PI
22 : #define M_PI 3.14159265358979323846
23 : #endif
24 :
25 : #include "grpp_angular_integrals.h"
26 : #include "grpp_binomial.h"
27 : #include "grpp_norm_gaussian.h"
28 : #include "grpp_radial_type1_integral.h"
29 : #include "grpp_type1_mcmurchie_davidson.h"
30 : #include "grpp_utils.h"
31 : #include "libgrpp.h"
32 : #include "libgrpp_types.h"
33 : /* for the old (numerical) version:
34 : void evaluate_type1_integral_primitive_gaussians(double *A, int n_cart_A, int
35 : *cart_list_A, double alpha_A, double *B, int n_cart_B, int *cart_list_B, double
36 : alpha_B, double *C, libgrpp_potential_t *potential, double *matrix);
37 : */
38 :
39 : extern void libgrpp_delete_radial_type1_integrals(radial_type1_table_t *table);
40 :
41 : void libgrpp_evaluate_radially_local_potential_integral_primitive_gaussians(
42 : double *A, int n_cart_A, int *cart_list_A, double alpha_A, double *B,
43 : int n_cart_B, int *cart_list_B, double alpha_B, double *C,
44 : double (*potential)(double r, void *params), void *potential_params,
45 : double *matrix);
46 :
47 : static double evaluate_pseudopotential(double r, void *params);
48 :
49 : /**
50 : * Evaluation of type 1 RPP integrals (scalar-relativistic radially local RPP).
51 : */
52 0 : void libgrpp_type1_integrals(libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B,
53 : double *rpp_origin, libgrpp_potential_t *potential,
54 : double *matrix) {
55 0 : assert(libgrpp_is_initialized());
56 :
57 0 : int size_A = libgrpp_get_shell_size(shell_A);
58 0 : int size_B = libgrpp_get_shell_size(shell_B);
59 0 : memset(matrix, 0, size_A * size_B * sizeof(double));
60 :
61 0 : if (potential == NULL) {
62 : return;
63 : }
64 :
65 : /*
66 : * RPP terms with n = 1, 2 are evaluated in a completely analytic manner
67 : * using the Obara-Saika-like recurrence relations
68 : */
69 :
70 0 : double *buf = calloc(size_A * size_B, sizeof(double));
71 :
72 0 : for (int k = 0; k < potential->num_primitives; k++) {
73 0 : double pot_coef = potential->coeffs[k];
74 0 : double pot_alpha = potential->alpha[k];
75 0 : int pot_n = potential->powers[k];
76 :
77 0 : libgrpp_type1_integrals_mcmurchie_davidson_1978(
78 : shell_A, shell_B, rpp_origin, pot_alpha, pot_n, buf);
79 :
80 0 : libgrpp_daxpy(size_A * size_B, pot_coef, buf, matrix);
81 : }
82 :
83 0 : free(buf);
84 :
85 : /*
86 : * old (numerical) version
87 : */
88 :
89 : /*for (int i = 0; i < shell_A->num_primitives; i++) {
90 : for (int j = 0; j < shell_B->num_primitives; j++) {
91 : double coef_A_i = shell_A->coeffs[i];
92 : double coef_B_j = shell_B->coeffs[j];
93 :
94 : if (fabs(coef_A_i * coef_B_j) < 1e-15) {
95 : continue;
96 : }
97 :
98 : evaluate_type1_integral_primitive_gaussians(
99 : shell_A->origin, size_A, shell_A->cart_list,
100 : shell_A->alpha[i], shell_B->origin, size_B, shell_B->cart_list,
101 : shell_B->alpha[j], rpp_origin, potential, buf
102 : );
103 :
104 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf, matrix);
105 : }
106 : }*/
107 : }
108 :
109 : /**
110 : * Evaluation of type 1 RPP integrals (scalar-relativistic radially local RPP)
111 : * for the pair of shells constructed from primitive Gaussians.
112 : */
113 0 : void evaluate_type1_integral_primitive_gaussians(
114 : double *A, int n_cart_A, int *cart_list_A, double alpha_A, double *B,
115 : int n_cart_B, int *cart_list_B, double alpha_B, double *C,
116 : libgrpp_potential_t *potential, double *matrix) {
117 0 : libgrpp_potential_t *potential_shrinked = libgrpp_shrink_potential(potential);
118 :
119 0 : libgrpp_evaluate_radially_local_potential_integral_primitive_gaussians(
120 : A, n_cart_A, cart_list_A, alpha_A, B, n_cart_B, cart_list_B, alpha_B, C,
121 : evaluate_pseudopotential, potential_shrinked, matrix);
122 :
123 0 : libgrpp_delete_potential(potential_shrinked);
124 0 : }
125 :
126 0 : static double evaluate_pseudopotential(double r, void *params) {
127 0 : libgrpp_potential_t *potential = (libgrpp_potential_t *)params;
128 :
129 0 : double u = libgrpp_potential_value(potential, r);
130 :
131 0 : return u;
132 : }
133 :
134 : /**
135 : * Evaluation of AO integrals for an arbitrary radially-local operator
136 : * for the pair of shells constructed from primitive Gaussians.
137 : */
138 0 : void libgrpp_evaluate_radially_local_potential_integral_primitive_gaussians(
139 : double *A, int n_cart_A, int *cart_list_A, double alpha_A, double *B,
140 : int n_cart_B, int *cart_list_B, double alpha_B, double *C,
141 : double (*potential)(double r, void *params), void *potential_params,
142 : double *matrix) {
143 0 : assert(n_cart_A > 0);
144 0 : assert(n_cart_B > 0);
145 :
146 0 : memset(matrix, 0, n_cart_A * n_cart_B * sizeof(double));
147 :
148 0 : double CA_x = C[0] - A[0];
149 0 : double CA_y = C[1] - A[1];
150 0 : double CA_z = C[2] - A[2];
151 0 : double CB_x = C[0] - B[0];
152 0 : double CB_y = C[1] - B[1];
153 0 : double CB_z = C[2] - B[2];
154 0 : double CA_2 = CA_x * CA_x + CA_y * CA_y + CA_z * CA_z;
155 0 : double CB_2 = CB_x * CB_x + CB_y * CB_y + CB_z * CB_z;
156 :
157 0 : double kx = -2.0 * (alpha_A * CA_x + alpha_B * CB_x);
158 0 : double ky = -2.0 * (alpha_A * CA_y + alpha_B * CB_y);
159 0 : double kz = -2.0 * (alpha_A * CA_z + alpha_B * CB_z);
160 0 : double k = sqrt(kx * kx + ky * ky + kz * kz);
161 0 : double kvec[3];
162 0 : kvec[0] = kx;
163 0 : kvec[1] = ky;
164 0 : kvec[2] = kz;
165 :
166 0 : int L_A = cart_list_A[0] + cart_list_A[1] + cart_list_A[2];
167 0 : int L_B = cart_list_B[0] + cart_list_B[1] + cart_list_B[2];
168 :
169 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
170 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
171 0 : double D_ABC = 4 * M_PI * N_A * N_B;
172 :
173 0 : int lambda_max = L_A + L_B;
174 0 : int n_max = lambda_max;
175 : // create_real_spherical_harmonic_coeffs_tables(lambda_max);
176 :
177 : /*
178 : * pre-compute type 1 radial integrals
179 : */
180 0 : radial_type1_table_t *radial_table = libgrpp_tabulate_radial_type1_integrals(
181 : lambda_max, n_max, CA_2, CB_2, alpha_A, alpha_B, k, D_ABC, potential,
182 : potential_params);
183 :
184 : /*
185 : * main loop
186 : * over shell pairs
187 : */
188 0 : for (int icart = 0; icart < n_cart_A; icart++) {
189 0 : for (int jcart = 0; jcart < n_cart_B; jcart++) {
190 :
191 0 : double chi_AB = 0.0;
192 :
193 0 : int n_A = cart_list_A[3 * icart + 0];
194 0 : int l_A = cart_list_A[3 * icart + 1];
195 0 : int m_A = cart_list_A[3 * icart + 2];
196 0 : int n_B = cart_list_B[3 * jcart + 0];
197 0 : int l_B = cart_list_B[3 * jcart + 1];
198 0 : int m_B = cart_list_B[3 * jcart + 2];
199 :
200 0 : for (int a = 0; a <= n_A; a++) {
201 0 : double C_nA_a = libgrpp_binomial(n_A, a);
202 0 : double pow_CA_x = pow(CA_x, n_A - a);
203 :
204 0 : for (int b = 0; b <= l_A; b++) {
205 0 : double C_lA_b = libgrpp_binomial(l_A, b);
206 0 : double pow_CA_y = pow(CA_y, l_A - b);
207 :
208 0 : for (int c = 0; c <= m_A; c++) {
209 0 : double C_mA_c = libgrpp_binomial(m_A, c);
210 0 : double pow_CA_z = pow(CA_z, m_A - c);
211 :
212 0 : for (int d = 0; d <= n_B; d++) {
213 0 : double C_nB_d = libgrpp_binomial(n_B, d);
214 0 : double pow_CB_x = pow(CB_x, n_B - d);
215 :
216 0 : for (int e = 0; e <= l_B; e++) {
217 0 : double C_lB_e = libgrpp_binomial(l_B, e);
218 0 : double pow_CB_y = pow(CB_y, l_B - e);
219 :
220 0 : for (int f = 0; f <= m_B; f++) {
221 0 : double C_mB_f = libgrpp_binomial(m_B, f);
222 0 : double pow_CB_z = pow(CB_z, m_B - f);
223 :
224 0 : double factor = C_nA_a * C_lA_b * C_mA_c * C_nB_d * C_lB_e *
225 0 : C_mB_f * pow_CA_x * pow_CA_y * pow_CA_z *
226 0 : pow_CB_x * pow_CB_y * pow_CB_z;
227 :
228 0 : if (fabs(factor) < 1e-13) {
229 0 : continue;
230 : }
231 :
232 0 : int N = a + b + c + d + e + f;
233 0 : double sum_omega_Q = 0.0;
234 0 : for (int lambda = 0; lambda <= lambda_max; lambda++) {
235 :
236 0 : double Q = libgrpp_get_radial_type1_integral(radial_table,
237 : lambda, N);
238 0 : if (fabs(Q) < 1e-16) {
239 0 : continue;
240 : }
241 :
242 0 : double omega = libgrpp_angular_type1_integral(
243 : lambda, a + d, b + e, c + f, kvec);
244 :
245 0 : sum_omega_Q += omega * Q;
246 : }
247 :
248 0 : chi_AB += factor * sum_omega_Q;
249 : }
250 : }
251 : }
252 : }
253 : }
254 : }
255 :
256 0 : matrix[icart * n_cart_B + jcart] = chi_AB;
257 : }
258 : }
259 :
260 0 : libgrpp_delete_radial_type1_integrals(radial_table);
261 0 : }
|
