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 : * Calculation of momentum integrals.
17 : *
18 : * For details, see:
19 : * T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-Structure Theory,
20 : * John Wiley & Sons Ltd, 2000.
21 : * Chapter 9.3.4, "Momentum and kinetic-energy integrals"
22 : */
23 : #include <math.h>
24 : #include <stdlib.h>
25 : #include <string.h>
26 : #ifndef M_PI
27 : #define M_PI 3.14159265358979323846
28 : #endif
29 : #include "grpp_momentum.h"
30 :
31 : #include "grpp_norm_gaussian.h"
32 : #include "grpp_utils.h"
33 : #include "libgrpp.h"
34 :
35 : static void momentum_integrals_shell_pair_obara_saika(
36 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
37 : double alpha_B, double *momentum_x_matrix, double *momentum_y_matrix,
38 : double *momentum_z_matrix);
39 :
40 : /**
41 : * returns imaginary(!) part of integrals over the momentum operator p = -i
42 : * \hbar \nabla. The "minus" sign is included.
43 : */
44 0 : void libgrpp_momentum_integrals(libgrpp_shell_t *shell_A,
45 : libgrpp_shell_t *shell_B,
46 : double *momentum_x_matrix,
47 : double *momentum_y_matrix,
48 : double *momentum_z_matrix) {
49 0 : int size_A = libgrpp_get_shell_size(shell_A);
50 0 : int size_B = libgrpp_get_shell_size(shell_B);
51 :
52 0 : double *buf_x = calloc(size_A * size_B, sizeof(double));
53 0 : double *buf_y = calloc(size_A * size_B, sizeof(double));
54 0 : double *buf_z = calloc(size_A * size_B, sizeof(double));
55 :
56 0 : memset(momentum_x_matrix, 0, size_A * size_B * sizeof(double));
57 0 : memset(momentum_y_matrix, 0, size_A * size_B * sizeof(double));
58 0 : memset(momentum_z_matrix, 0, size_A * size_B * sizeof(double));
59 :
60 : // loop over primitives in contractions
61 0 : for (int i = 0; i < shell_A->num_primitives; i++) {
62 0 : for (int j = 0; j < shell_B->num_primitives; j++) {
63 0 : double alpha_i = shell_A->alpha[i];
64 0 : double alpha_j = shell_B->alpha[j];
65 0 : double coef_A_i = shell_A->coeffs[i];
66 0 : double coef_B_j = shell_B->coeffs[j];
67 :
68 0 : momentum_integrals_shell_pair_obara_saika(shell_A, shell_B, alpha_i,
69 : alpha_j, buf_x, buf_y, buf_z);
70 :
71 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_x,
72 : momentum_x_matrix);
73 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_y,
74 : momentum_y_matrix);
75 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_z,
76 : momentum_z_matrix);
77 : }
78 : }
79 :
80 0 : free(buf_x);
81 0 : free(buf_y);
82 0 : free(buf_z);
83 0 : }
84 :
85 0 : static void momentum_integrals_shell_pair_obara_saika(
86 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
87 : double alpha_B, double *momentum_x_matrix, double *momentum_y_matrix,
88 : double *momentum_z_matrix) {
89 0 : int size_A = libgrpp_get_shell_size(shell_A);
90 0 : int size_B = libgrpp_get_shell_size(shell_B);
91 0 : int L_A = shell_A->L;
92 0 : int L_B = shell_B->L;
93 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
94 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
95 :
96 0 : double p = alpha_A + alpha_B;
97 0 : double mu = alpha_A * alpha_B / (alpha_A + alpha_B);
98 0 : double *A = shell_A->origin;
99 0 : double *B = shell_B->origin;
100 :
101 : // calculate S_ij
102 0 : double S[3][LIBGRPP_MAX_BASIS_L + 1][LIBGRPP_MAX_BASIS_L + 1];
103 :
104 0 : for (int coord = 0; coord < 3; coord++) {
105 0 : double P = (alpha_A * A[coord] + alpha_B * B[coord]) / p;
106 :
107 0 : double X_AB = A[coord] - B[coord];
108 0 : double X_PA = P - A[coord];
109 0 : double X_PB = P - B[coord];
110 0 : double pfac = 1.0 / (2.0 * p);
111 :
112 0 : for (int i = 0; i <= L_A + 1; i++) {
113 0 : for (int j = 0; j <= L_B + 1; j++) {
114 0 : double S_ij = 0.0;
115 :
116 0 : if (i + j == 0) {
117 0 : S[coord][0][0] = sqrt(M_PI / p) * exp(-mu * X_AB * X_AB);
118 0 : continue;
119 : }
120 :
121 0 : if (i == 0) { // upward by j
122 0 : S_ij += X_PB * S[coord][i][j - 1];
123 0 : if (j - 1 > 0) {
124 0 : S_ij += (j - 1) * pfac * S[coord][i][j - 2];
125 : }
126 : } else { // upward by i
127 0 : S_ij += X_PA * S[coord][i - 1][j];
128 0 : if (i - 1 > 0) {
129 0 : S_ij += (i - 1) * pfac * S[coord][i - 2][j];
130 : }
131 0 : if (j > 0) {
132 0 : S_ij += j * pfac * S[coord][i - 1][j - 1];
133 : }
134 : }
135 :
136 0 : S[coord][i][j] = S_ij;
137 : }
138 : }
139 : }
140 :
141 : // calculate D^1_ij
142 :
143 : double D1[3][LIBGRPP_MAX_BASIS_L][LIBGRPP_MAX_BASIS_L];
144 :
145 0 : for (int coord = 0; coord < 3; coord++) {
146 0 : for (int i = 0; i <= L_A; i++) {
147 0 : for (int j = 0; j <= L_B; j++) {
148 :
149 0 : double D1_ij = 0.0;
150 0 : D1_ij += 2.0 * alpha_A * S[coord][i + 1][j];
151 0 : if (i >= 1) {
152 0 : D1_ij -= i * S[coord][i - 1][j];
153 : }
154 :
155 0 : D1[coord][i][j] = D1_ij;
156 : }
157 : }
158 : }
159 :
160 : // loop over cartesian functions inside the shells
161 0 : for (int m = 0; m < size_A; m++) {
162 0 : for (int n = 0; n < size_B; n++) {
163 0 : int n_A = shell_A->cart_list[3 * m + 0];
164 0 : int l_A = shell_A->cart_list[3 * m + 1];
165 0 : int m_A = shell_A->cart_list[3 * m + 2];
166 0 : int n_B = shell_B->cart_list[3 * n + 0];
167 0 : int l_B = shell_B->cart_list[3 * n + 1];
168 0 : int m_B = shell_B->cart_list[3 * n + 2];
169 :
170 0 : momentum_x_matrix[m * size_B + n] =
171 0 : -N_A * N_B * D1[0][n_A][n_B] * S[1][l_A][l_B] * S[2][m_A][m_B];
172 0 : momentum_y_matrix[m * size_B + n] =
173 0 : -N_A * N_B * S[0][n_A][n_B] * D1[1][l_A][l_B] * S[2][m_A][m_B];
174 0 : momentum_z_matrix[m * size_B + n] =
175 0 : -N_A * N_B * S[0][n_A][n_B] * S[1][l_A][l_B] * D1[2][m_A][m_B];
176 : }
177 : }
178 0 : }
|
