LIBINT 2.9.0
VRR_GTG_1d_xx_xx.h
1/*
2 * Copyright (C) 2004-2024 Edward F. Valeev
3 *
4 * This file is part of Libint library.
5 *
6 * Libint library is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * Libint library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with Libint library. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef _libint2_src_lib_libint_vrrgtg1dxxxx_h_
22#define _libint2_src_lib_libint_vrrgtg1dxxxx_h_
23
24#include <libint2.h>
25#include <util_types.h>
26
27#include <cassert>
28#include <cstdlib>
29
30namespace libint2 {
31
42template <unsigned int CartesianAxis, int La, int Lb, int Lc, int Ld,
43 bool vectorize>
44struct VRR_GTG_1d_xx_xx {
45 static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
46 const LIBINT2_REALTYPE* src0) {
47 enum XYZ { x = 0, y = 1, z = 2 };
48 assert(CartesianAxis == x || CartesianAxis == y || CartesianAxis == z);
49 assert(!vectorize);
50
51 const unsigned int veclen = vectorize ? inteval->veclen : 1;
52
53 // corner case: (00|00)
54 if (La == 0 && Lb == 0 && Lc == 0 && Ld == 0) {
55 for (unsigned int v = 0; v != veclen; ++v) target[v] = src0[v];
56 return;
57 }
58
59 //---------------------------------------------
60 // Part (1): build (a+b 0|c+d 0)
61 //---------------------------------------------
62
63 LIBINT2_REALTYPE apb_0_GTG_cpd_0[La + Lb + 1][Lc + Ld + 1];
64 apb_0_GTG_cpd_0[0][0] = src0[0];
65
66 const LIBINT2_REALTYPE *pfac0_0, *pfac0_1;
67 const LIBINT2_REALTYPE* pfac1_0 = inteval->R12kG12_pfac1_0;
68 const LIBINT2_REALTYPE* pfac1_1 = inteval->R12kG12_pfac1_1;
69 const LIBINT2_REALTYPE* pfac2 = inteval->R12kG12_pfac2;
70 switch (CartesianAxis) {
71 case x:
72 pfac0_0 = inteval->R12kG12_pfac0_0_x;
73 pfac0_1 = inteval->R12kG12_pfac0_1_x;
74 break;
75 case y:
76 pfac0_0 = inteval->R12kG12_pfac0_0_y;
77 pfac0_1 = inteval->R12kG12_pfac0_1_y;
78 break;
79 case z:
80 pfac0_0 = inteval->R12kG12_pfac0_0_z;
81 pfac0_1 = inteval->R12kG12_pfac0_1_z;
82 break;
83 default:
84 assert(false);
85 }
86
87 // build (0 0|1 0)
88 if (Lc + Ld > 0) {
89 apb_0_GTG_cpd_0[0][1] = pfac0_1[0] * apb_0_GTG_cpd_0[0][0];
90#if LIBINT2_FLOP_COUNT
91 inteval->nflops[0] += 1;
92#endif
93 }
94
95 // build (0 0|c+d 0)
96 if (Lc + Ld > 1) {
97 for (int c_plus_d = 1; c_plus_d != Lc + Ld; ++c_plus_d) {
98 apb_0_GTG_cpd_0[0][c_plus_d + 1] =
99 pfac0_1[0] * apb_0_GTG_cpd_0[0][c_plus_d] +
100 c_plus_d * pfac1_1[0] * apb_0_GTG_cpd_0[0][c_plus_d - 1];
101 }
102#if LIBINT2_FLOP_COUNT
103 inteval->nflops[0] += 4 * (Lc + Ld - 1);
104#endif
105 }
106
107 // build (1 0|0 0)
108 if (La + Lb > 0) {
109 apb_0_GTG_cpd_0[1][0] = pfac0_0[0] * apb_0_GTG_cpd_0[0][0];
110#if LIBINT2_FLOP_COUNT
111 inteval->nflops[0] += 1;
112#endif
113 }
114
115 // build (a+b 0|0 0)
116 if (La + Lb > 1) {
117 for (int a_plus_b = 1; a_plus_b != La + Lb; ++a_plus_b) {
118 apb_0_GTG_cpd_0[a_plus_b + 1][0] =
119 pfac0_0[0] * apb_0_GTG_cpd_0[a_plus_b][0] +
120 a_plus_b * pfac1_0[0] * apb_0_GTG_cpd_0[a_plus_b - 1][0];
121 }
122#if LIBINT2_FLOP_COUNT
123 inteval->nflops[0] += 4 * (La + Lb - 1);
124#endif
125 }
126
127 // build (1 0|c+d 0)
128 if (La + Lb > 0 && Lc + Ld > 0) {
129 for (int c_plus_d = 1; c_plus_d <= Lc + Ld; ++c_plus_d) {
130 apb_0_GTG_cpd_0[1][c_plus_d] =
131 pfac0_0[0] * apb_0_GTG_cpd_0[0][c_plus_d] +
132 c_plus_d * pfac2[0] * apb_0_GTG_cpd_0[0][c_plus_d - 1];
133 }
134#if LIBINT2_FLOP_COUNT
135 inteval->nflops[0] += 4 * (Lc + Ld - 1);
136#endif
137 }
138
139 // build (a+b 0|c+d 0)
140 if (La + Lb > 1 && Lc + Ld > 0) {
141 for (int a_plus_b = 1; a_plus_b != La + Lb; ++a_plus_b) {
142 for (int c_plus_d = 1; c_plus_d <= Lc + Ld; ++c_plus_d) {
143 apb_0_GTG_cpd_0[a_plus_b + 1][c_plus_d] =
144 pfac0_0[0] * apb_0_GTG_cpd_0[a_plus_b][c_plus_d] +
145 a_plus_b * pfac1_0[0] * apb_0_GTG_cpd_0[a_plus_b - 1][c_plus_d] +
146 c_plus_d * pfac2[0] * apb_0_GTG_cpd_0[a_plus_b][c_plus_d - 1];
147 }
148 }
149#if LIBINT2_FLOP_COUNT
150 inteval->nflops[0] += 7 * (La + Lb - 1) * (Lc + Ld - 1);
151#endif
152 }
153
154 //---------------------------------------------
155 // Part (2): build (a b|c+d 0)
156 //---------------------------------------------
157
158 const LIBINT2_REALTYPE* AB;
159 switch (CartesianAxis) {
160 case x:
161 AB = inteval->AB_x;
162 break;
163 case y:
164 AB = inteval->AB_y;
165 break;
166 case z:
167 AB = inteval->AB_z;
168 break;
169 default:
170 assert(false);
171 }
172
173 LIBINT2_REALTYPE a_b_GTG_cpd_0[La + 1][Lb + 1][Lc + Ld + 1];
174 for (int c_plus_d = 0; c_plus_d <= Lc + Ld; ++c_plus_d) {
175 // copy (a+b 0| to a local 0,a+b buffer
176 LIBINT2_REALTYPE b_a_GTG[La + Lb + 1][La + Lb + 1];
177 for (int a_plus_b = 0; a_plus_b <= La + Lb; ++a_plus_b) {
178 b_a_GTG[0][a_plus_b] = apb_0_GTG_cpd_0[a_plus_b][c_plus_d];
179 }
180 // use HRR to compute b,a
181 for (int b = 1; b <= Lb; ++b) {
182 for (int a = 0; a <= La + Lb - b; ++a) {
183 b_a_GTG[b][a] = b_a_GTG[b - 1][a + 1] + AB[0] * b_a_GTG[b - 1][a];
184 }
185#if LIBINT2_FLOP_COUNT
186 inteval->nflops[0] += 2 * (La + Lb - b + 1);
187#endif
188 }
189 // copy b,a to (a b|
190 for (int b = 0; b <= Lb; ++b) {
191 for (int a = 0; a <= La; ++a) {
192 a_b_GTG_cpd_0[a][b][c_plus_d] = b_a_GTG[b][a];
193 }
194 }
195 }
196
197 //---------------------------------------------
198 // Part (3): build (a b|c d)
199 //---------------------------------------------
200
201 const LIBINT2_REALTYPE* CD;
202 switch (CartesianAxis) {
203 case x:
204 CD = inteval->CD_x;
205 break;
206 case y:
207 CD = inteval->CD_y;
208 break;
209 case z:
210 CD = inteval->CD_z;
211 break;
212 default:
213 assert(false);
214 }
215
216 LIBINT2_REALTYPE* target_a_b_blk_ptr = target;
217 const int Nd = (Ld + 1);
218 const int Ncd = (Lc + 1) * Nd;
219 for (int a = 0; a <= La; ++a) {
220 for (int b = 0; b <= Lb; ++b, target_a_b_blk_ptr += Ncd) {
221 // copy |c+d 0) to a local 0,c+d buffer
222 LIBINT2_REALTYPE d_c_GTG[Lc + Ld + 1][Lc + Ld + 1];
223 for (int c_plus_d = 0; c_plus_d <= Lc + Ld; ++c_plus_d) {
224 d_c_GTG[0][c_plus_d] = a_b_GTG_cpd_0[a][b][c_plus_d];
225 }
226 // use HRR to compute d,c
227 for (int d = 1; d <= Ld; ++d) {
228 for (int c = 0; c <= Lc + Ld - d; ++c) {
229 d_c_GTG[d][c] = d_c_GTG[d - 1][c + 1] + CD[0] * d_c_GTG[d - 1][c];
230 }
231#if LIBINT2_FLOP_COUNT
232 inteval->nflops[0] += 2 * (Lc + Ld - d + 1);
233#endif
234 }
235 // copy d,c to |c d)
236 for (int d = 0; d <= Ld; ++d) {
237 for (int c = 0, cd = d; c <= Lc; ++c, cd += Nd) {
238 target_a_b_blk_ptr[cd] = d_c_GTG[d][c];
239 }
240 }
241 }
242 }
243
244 // done
245 }
246};
247
248}; // namespace libint2
249
250#endif // header guard
libint2
Defaults definitions for various parameters assumed by Libint.
Definition algebra.cc:24
libint2::VRR_GTG_1d_xx_xx
builds (ab| GTG_1d |cd), the shell set of 2-dimensional integrals needed for Rys quadrature evaluatio...
Definition VRR_GTG_1d_xx_xx.h:44
Generated by doxygen 1.10.0