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 kinetic-energy 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 :
27 : #ifndef M_PI
28 : #define M_PI 3.14159265358979323846
29 : #endif
30 :
31 : #include "grpp_kinetic.h"
32 :
33 : #include "grpp_norm_gaussian.h"
34 : #include "grpp_utils.h"
35 : #include "libgrpp.h"
36 :
37 : static void kinetic_energy_integrals_shell_pair_obara_saika(
38 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
39 : double alpha_B, double *kinetic_matrix);
40 :
41 0 : void libgrpp_kinetic_energy_integrals(libgrpp_shell_t *shell_A,
42 : libgrpp_shell_t *shell_B,
43 : double *kinetic_matrix) {
44 0 : int size_A = libgrpp_get_shell_size(shell_A);
45 0 : int size_B = libgrpp_get_shell_size(shell_B);
46 :
47 0 : double *buf = calloc(size_A * size_B, sizeof(double));
48 :
49 0 : memset(kinetic_matrix, 0, size_A * size_B * sizeof(double));
50 :
51 : // loop over primitives in contractions
52 0 : for (int i = 0; i < shell_A->num_primitives; i++) {
53 0 : for (int j = 0; j < shell_B->num_primitives; j++) {
54 0 : double alpha_i = shell_A->alpha[i];
55 0 : double alpha_j = shell_B->alpha[j];
56 0 : double coef_A_i = shell_A->coeffs[i];
57 0 : double coef_B_j = shell_B->coeffs[j];
58 :
59 0 : kinetic_energy_integrals_shell_pair_obara_saika(shell_A, shell_B, alpha_i,
60 : alpha_j, buf);
61 :
62 0 : libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf, kinetic_matrix);
63 : }
64 : }
65 :
66 0 : free(buf);
67 0 : }
68 :
69 0 : static void kinetic_energy_integrals_shell_pair_obara_saika(
70 : libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
71 : double alpha_B, double *kinetic_matrix) {
72 0 : int size_A = libgrpp_get_shell_size(shell_A);
73 0 : int size_B = libgrpp_get_shell_size(shell_B);
74 0 : int L_A = shell_A->L;
75 0 : int L_B = shell_B->L;
76 0 : double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
77 0 : double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
78 :
79 0 : double p = alpha_A + alpha_B;
80 0 : double mu = alpha_A * alpha_B / (alpha_A + alpha_B);
81 0 : double *A = shell_A->origin;
82 0 : double *B = shell_B->origin;
83 :
84 : // calculate S_ij
85 0 : double S[3][LIBGRPP_MAX_BASIS_L + 2][LIBGRPP_MAX_BASIS_L + 2];
86 :
87 0 : for (int coord = 0; coord < 3; coord++) {
88 0 : double P = (alpha_A * A[coord] + alpha_B * B[coord]) / p;
89 :
90 0 : double X_AB = A[coord] - B[coord];
91 0 : double X_PA = P - A[coord];
92 0 : double X_PB = P - B[coord];
93 0 : double pfac = 1.0 / (2.0 * p);
94 :
95 0 : for (int i = 0; i <= L_A + 2; i++) {
96 0 : for (int j = 0; j <= L_B + 2; j++) {
97 0 : double S_ij = 0.0;
98 :
99 0 : if (i + j == 0) {
100 0 : S[coord][0][0] = sqrt(M_PI / p) * exp(-mu * X_AB * X_AB);
101 0 : continue;
102 : }
103 :
104 0 : if (i == 0) { // upward by j
105 0 : S_ij += X_PB * S[coord][i][j - 1];
106 0 : if (j - 1 > 0) {
107 0 : S_ij += (j - 1) * pfac * S[coord][i][j - 2];
108 : }
109 : } else { // upward by i
110 0 : S_ij += X_PA * S[coord][i - 1][j];
111 0 : if (i - 1 > 0) {
112 0 : S_ij += (i - 1) * pfac * S[coord][i - 2][j];
113 : }
114 0 : if (j > 0) {
115 0 : S_ij += j * pfac * S[coord][i - 1][j - 1];
116 : }
117 : }
118 :
119 0 : S[coord][i][j] = S_ij;
120 : }
121 : }
122 : }
123 :
124 : // calculate D^2_ij
125 :
126 : double D2[3][LIBGRPP_MAX_BASIS_L][LIBGRPP_MAX_BASIS_L];
127 :
128 0 : for (int coord = 0; coord < 3; coord++) {
129 0 : for (int i = 0; i <= L_A; i++) {
130 0 : for (int j = 0; j <= L_B; j++) {
131 :
132 0 : double D2_ij = 0.0;
133 0 : D2_ij += 4.0 * alpha_A * alpha_A * S[coord][i + 2][j];
134 0 : D2_ij -= 2.0 * alpha_A * (2 * i + 1) * S[coord][i][j];
135 0 : if (i >= 2) {
136 0 : D2_ij += i * (i - 1) * S[coord][i - 2][j];
137 : }
138 :
139 0 : D2[coord][i][j] = D2_ij;
140 : }
141 : }
142 : }
143 :
144 : // loop over cartesian functions inside the shells
145 0 : for (int m = 0; m < size_A; m++) {
146 0 : for (int n = 0; n < size_B; n++) {
147 0 : int n_A = shell_A->cart_list[3 * m + 0];
148 0 : int l_A = shell_A->cart_list[3 * m + 1];
149 0 : int m_A = shell_A->cart_list[3 * m + 2];
150 0 : int n_B = shell_B->cart_list[3 * n + 0];
151 0 : int l_B = shell_B->cart_list[3 * n + 1];
152 0 : int m_B = shell_B->cart_list[3 * n + 2];
153 :
154 0 : kinetic_matrix[m * size_B + n] =
155 0 : -0.5 * N_A * N_B *
156 0 : (D2[0][n_A][n_B] * S[1][l_A][l_B] * S[2][m_A][m_B] +
157 0 : S[0][n_A][n_B] * D2[1][l_A][l_B] * S[2][m_A][m_B] +
158 0 : S[0][n_A][n_B] * S[1][l_A][l_B] * D2[2][m_A][m_B]);
159 : }
160 : }
161 0 : }
|
