1 /*
2 * Copyright (C) 2004-2021 Edward F. Valeev
3 *
4 * This file is part of Libint.
5 *
6 * Libint is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU 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 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 General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21 #include <boost/preprocessor/list/for_each.hpp>
22
23 #include <rr.h>
24 #include <comp_xyz.h>
25 #include <integral_1_1.h>
26 #include <prefactors.h>
27 #include <strategy.h>
28 #include <dg.h>
29 #include <rr.h>
30 #include <graph_registry.h>
31 #include <intset_to_ints.h>
32 #include <uncontract.h>
33 #include <singl_stack.h>
34
35 using namespace std;
36
37 namespace libint2 {
38
39 template<>
compute(const CGShell & a,const CGShell & b,const OperType &)40 void CR_XYZ_1_1<CGShell,OverlapOper,EmptySet>::compute(const CGShell& a, const CGShell& b, const OperType&) {
41 using namespace libint2::algebra;
42 using namespace libint2::prefactor;
43 using namespace libint2::braket;
44
45 // why do I need this? when this function is in header clang-602.0.49 on OS X works without this
46 using libint2::algebra::operator*;
47
48 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_X>,OverlapOper,EmptySet>> factory_X(this);
49 auto ab_X = factory_X.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm());
50
51 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Y>,OverlapOper,EmptySet>> factory_Y(this);
52 auto ab_Y = factory_Y.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm());
53
54 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Z>,OverlapOper,EmptySet>> factory_Z(this);
55 auto ab_Z = factory_Z.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm());
56 }
57
58 template<>
compute(const CGShell1d<CartesianAxis_X> & a,const CGShell1d<CartesianAxis_X> & b,const OperType &)59 void CR_XYZ_1_1<CGShell1d<CartesianAxis_X>,OverlapOper,EmptySet>::compute(const CGShell1d<CartesianAxis_X>& a,
60 const CGShell1d<CartesianAxis_X>& b,
61 const OperType&) {
62 this->add_child(prefactor::Scalar("_0_Overlap_0_x"));
63 }
64 template<>
compute(const CGShell1d<CartesianAxis_Y> & a,const CGShell1d<CartesianAxis_Y> & b,const OperType &)65 void CR_XYZ_1_1<CGShell1d<CartesianAxis_Y>,OverlapOper,EmptySet>::compute(const CGShell1d<CartesianAxis_Y>& a,
66 const CGShell1d<CartesianAxis_Y>& b,
67 const OperType&) {
68 this->add_child(prefactor::Scalar("_0_Overlap_0_y"));
69 }
70 template<>
compute(const CGShell1d<CartesianAxis_Z> & a,const CGShell1d<CartesianAxis_Z> & b,const OperType &)71 void CR_XYZ_1_1<CGShell1d<CartesianAxis_Z>,OverlapOper,EmptySet>::compute(const CGShell1d<CartesianAxis_Z>& a,
72 const CGShell1d<CartesianAxis_Z>& b,
73 const OperType&) {
74 this->add_child(prefactor::Scalar("_0_Overlap_0_z"));
75 }
76
77 template<>
compute(const CGF & a,const CGF & b,const OperType &)78 void CR_XYZ_1_1<CGF,OverlapOper,EmptySet>::compute(const CGF& a, const CGF& b, const OperType&) {
79
80 using namespace libint2::algebra;
81 using namespace libint2::prefactor;
82 using namespace libint2::braket;
83
84 // why do I need this? when this function is in header clang-602.0.49 on OS X works without this
85 using libint2::algebra::operator*;
86
87 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_X>,OverlapOper,EmptySet>> factory_X(this);
88 auto ab_X = factory_X.make_child(a[CartesianAxis_X],b[CartesianAxis_X]);
89
90 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Y>,OverlapOper,EmptySet>> factory_Y(this);
91 auto ab_Y = factory_Y.make_child(a[CartesianAxis_Y],b[CartesianAxis_Y]);
92
93 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Z>,OverlapOper,EmptySet>> factory_Z(this);
94 auto ab_Z = factory_Z.make_child(a[CartesianAxis_Z],b[CartesianAxis_Z]);
95
96 if (is_simple()) {
97 expr_ = ab_X * ab_Y * ab_Z;
98 nflops_ += 2;
99 }
100
101 } // end of Overlap compute
102
103 template<>
compute(const CGShell & a,const CGShell & b,const OperType &)104 void CR_XYZ_1_1<CGShell,KineticOper,EmptySet>::compute(const CGShell& a, const CGShell& b, const OperType&) {
105 using namespace libint2::algebra;
106 using namespace libint2::prefactor;
107 using namespace libint2::braket;
108
109 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_X>,OverlapOper,EmptySet>> factory_X(this);
110 auto ab_X = factory_X.make_child(a.norm() + a.deriv().norm() + 1,b.norm() + b.deriv().norm() + 1);
111
112 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Y>,OverlapOper,EmptySet>> factory_Y(this);
113 auto ab_Y = factory_Y.make_child(a.norm() + a.deriv().norm() + 1,b.norm() + b.deriv().norm() + 1);
114
115 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Z>,OverlapOper,EmptySet>> factory_Z(this);
116 auto ab_Z = factory_Z.make_child(a.norm() + a.deriv().norm() + 1,b.norm() + b.deriv().norm() + 1);
117 }
118
119
120
121 template<>
compute(const CGF & a,const CGF & b,const OperType &)122 void CR_XYZ_1_1<CGF,KineticOper,EmptySet>::compute(const CGF& a, const CGF& b, const OperType&) {
123
124 using namespace libint2::algebra;
125 // why do I need this? when this function is in header clang-602.0.49 on OS X works without this
126 using libint2::algebra::operator*;
127
128 // this is what I want this code to look like
129 //this->wedge(dot(Nabla(a), Nabla(b)), Scalar(-0.25), EmptySet(), OverlapOper());
130
131 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_X>,OverlapOper,EmptySet>> factory_X(this);
132 auto s_X = factory_X.make_child(a[CartesianAxis_X],b[CartesianAxis_X]);
133
134 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Y>,OverlapOper,EmptySet>> factory_Y(this);
135 auto s_Y = factory_Y.make_child(a[CartesianAxis_Y],b[CartesianAxis_Y]);
136
137 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Z>,OverlapOper,EmptySet>> factory_Z(this);
138 auto s_Z = factory_Z.make_child(a[CartesianAxis_Z],b[CartesianAxis_Z]);
139
140 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_X>,KineticOper,EmptySet>> tfactory_X(this);
141 auto t_X = tfactory_X.make_child(a[CartesianAxis_X],b[CartesianAxis_X]);
142
143 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Y>,KineticOper,EmptySet>> tfactory_Y(this);
144 auto t_Y = tfactory_Y.make_child(a[CartesianAxis_Y],b[CartesianAxis_Y]);
145
146 ChildFactory<ThisType,GenIntegralSet_1_1<CGF1d<CartesianAxis_Z>,KineticOper,EmptySet>> tfactory_Z(this);
147 auto t_Z = tfactory_Z.make_child(a[CartesianAxis_Z],b[CartesianAxis_Z]);
148
149 expr_ = t_X * s_Y * s_Z;
150 expr_ += s_X * t_Y * s_Z;
151 expr_ += s_X * s_Y * t_Z;
152 nflops_ += 8;
153 } // end of Kinetic 3-d compute
154
155 #define BOOST_PP_CR_XYZ_1_1_1D_KINETIC_COMPUTE(r,data,elem) \
156 template<> \
157 void CR_XYZ_1_1<CGF1d<CartesianAxis_ ## elem>,KineticOper,EmptySet>::compute(const BasisFunctionType& a, \
158 const BasisFunctionType& b, \
159 const OperType&) { \
160 using namespace libint2::algebra; \
161 using namespace libint2::prefactor; \
162 using libint2::algebra::operator*; \
163 \
164 ChildFactory<ThisType, \
165 GenIntegralSet_1_1<BasisFunctionType,OverlapOper,EmptySet>> factory(this); \
166 \
167 const BasisFunctionType& _1 = unit<BasisFunctionType>(0); \
168 auto ap1 = a + _1; \
169 auto bp1 = b + _1; \
170 auto ap1bp1 = factory.make_child(ap1,bp1); \
171 expr_ = Scalar(0.5) * Scalar("two_alpha0_bra") * Scalar("two_alpha0_ket") * ap1bp1; \
172 nflops_ += 3; \
173 auto am1 = a - _1; \
174 auto bm1 = b - _1; \
175 if (exists(am1)) { \
176 auto am1bp1 = factory.make_child(am1,bp1); \
177 expr_ -= Scalar(0.5 * a[0]) * Scalar("two_alpha0_ket") * am1bp1; \
178 nflops_ += 3; \
179 } \
180 if (exists(bm1)) { \
181 auto ap1bm1 = factory.make_child(ap1,bm1); \
182 expr_ -= Scalar(0.5 * b[0]) * Scalar("two_alpha0_bra") * ap1bm1; \
183 nflops_ += 3; \
184 } \
185 if (exists(am1) && exists(bm1)) { \
186 auto am1bm1 = factory.make_child(am1,bm1); \
187 expr_ += Scalar(0.5 * a[0] * b[0]) * am1bm1; \
188 nflops_ += 2; \
189 } \
190 }
191
192 #define BOOST_PP_XYZ_LIST (X, (Y, (Z, BOOST_PP_NIL)))
BOOST_PP_LIST_FOR_EACH(BOOST_PP_CR_XYZ_1_1_1D_KINETIC_COMPUTE,_,BOOST_PP_XYZ_LIST)193 BOOST_PP_LIST_FOR_EACH ( BOOST_PP_CR_XYZ_1_1_1D_KINETIC_COMPUTE, _, BOOST_PP_XYZ_LIST)
194 #undef BOOST_PP_XYZ_LIST
195
196
197 template<>
198 void CR_XYZ_1_1<CGShell,CartesianMultipoleOper<3u>,EmptySet>::compute(const CGShell& a,
199 const CGShell& b,
200 const OperType& oper) {
201 using namespace libint2::algebra;
202 using namespace libint2::prefactor;
203 using namespace libint2::braket;
204
205 // why do I need this? when this function is in header clang-602.0.49 on OS X works without this
206 using libint2::algebra::operator*;
207
208 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_X>,OverlapOper,EmptySet>> factory_X(this);
209 auto ab_X = factory_X.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm() + oper.descr().norm());
210
211 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Y>,OverlapOper,EmptySet>> factory_Y(this);
212 auto ab_Y = factory_Y.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm() + oper.descr().norm());
213
214 ChildFactory<ThisType,GenIntegralSet_1_1<CGShell1d<CartesianAxis_Z>,OverlapOper,EmptySet>> factory_Z(this);
215 auto ab_Z = factory_Z.make_child(a.norm() + a.deriv().norm(),b.norm() + b.deriv().norm() + oper.descr().norm());
216 }
217
218 template<>
compute(const CGF & a,const CGF & b,const OperType & oper)219 void CR_XYZ_1_1<CGF,CartesianMultipoleOper<3u>,EmptySet>::compute(const CGF& a,
220 const CGF& b,
221 const OperType& oper) {
222
223 using namespace libint2::algebra;
224 using namespace libint2::prefactor;
225 using namespace libint2::braket;
226
227 // why do I need this? when this function is in header clang-602.0.49 on OS X works without this
228 using libint2::algebra::operator*;
229
230 using Descr1d = CartesianMultipoleOper<1u>::Descriptor;
231
232 ChildFactory<ThisType,
233 GenIntegralSet_1_1<CGF1d<CartesianAxis_X>,
234 CartesianMultipoleOper<1u>,
235 EmptySet>
236 > factory_X(this);
237 auto ab_X = factory_X.make_child(a[CartesianAxis_X],b[CartesianAxis_X],
238 EmptySet(),
239 Descr1d(oper.descr()[CartesianAxis_X]));
240
241 ChildFactory<ThisType,
242 GenIntegralSet_1_1<CGF1d<CartesianAxis_Y>,
243 CartesianMultipoleOper<1u>,
244 EmptySet>
245 > factory_Y(this);
246 auto ab_Y = factory_Y.make_child(a[CartesianAxis_Y],b[CartesianAxis_Y],
247 EmptySet(),
248 Descr1d(oper.descr()[CartesianAxis_Y]));
249
250 ChildFactory<ThisType,
251 GenIntegralSet_1_1<CGF1d<CartesianAxis_Z>,
252 CartesianMultipoleOper<1u>,
253 EmptySet>
254 > factory_Z(this);
255 auto ab_Z = factory_Z.make_child(a[CartesianAxis_Z],b[CartesianAxis_Z],
256 EmptySet(),
257 Descr1d(oper.descr()[CartesianAxis_Z]));
258
259 expr_ = ab_X * ab_Y * ab_Z;
260 nflops_ += 2;
261
262 }
263
264 #define BOOST_PP_CR_XYZ_1_1_1D_CMULTIPOLE_COMPUTE(r,data,elem) \
265 template<> \
266 void CR_XYZ_1_1<CGF1d<CartesianAxis_ ## elem>,CartesianMultipoleOper<1u>,EmptySet>::compute( \
267 const BasisFunctionType& a, \
268 const BasisFunctionType& b, \
269 const OperType& oper) { \
270 using namespace libint2::algebra; \
271 using namespace libint2::prefactor; \
272 using libint2::algebra::operator*; \
273 \
274 const BasisFunctionType& _1 = unit<BasisFunctionType>(0); \
275 auto bp1 = b + _1; \
276 auto descr_m1 = oper.descr(); \
277 if (descr_m1[0] > 0) { \
278 ChildFactory<ThisType, \
279 GenIntegralSet_1_1<BasisFunctionType,OperType,EmptySet>> factory(this); \
280 descr_m1.dec(0); auto oper_m1 = OperType(descr_m1); \
281 expr_ = factory.make_child(a,bp1,EmptySet(),oper_m1) + \
282 Vector("BO")[CartesianAxis_ ## elem] * factory.make_child(a,b,EmptySet(),oper_m1); \
283 } else { \
284 ChildFactory<ThisType, \
285 GenIntegralSet_1_1<BasisFunctionType,OverlapOper,EmptySet>> sfactory(this); \
286 expr_ = Scalar(0u) + sfactory.make_child(a,b); \
287 } \
288 nflops_ += 2; \
289 }
290
291 #define BOOST_PP_XYZ_LIST (X, (Y, (Z, BOOST_PP_NIL)))
292 BOOST_PP_LIST_FOR_EACH ( BOOST_PP_CR_XYZ_1_1_1D_CMULTIPOLE_COMPUTE, _, BOOST_PP_XYZ_LIST)
293 #undef BOOST_PP_XYZ_LIST
294
295 } // namespace libint2
296