LIBINT  2.6.0
vrr_11_twoprep_11.h
1 /*
2  * Copyright (C) 2004-2019 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 #ifndef _libint2_src_bin_libint_vrr11twoprep11_h_
22 #define _libint2_src_bin_libint_vrr11twoprep11_h_
23 
24 #include <generic_rr.h>
25 #include <twoprep_11_11.h>
26 
27 using namespace std;
28 
29 namespace libint2 {
30 
34  template <class BFSet, int part, FunctionPosition where>
35  class VRR_11_TwoPRep_11 : public GenericRecurrenceRelation< VRR_11_TwoPRep_11<BFSet,part,where>,
36  BFSet,
37  GenIntegralSet_11_11<BFSet,TwoPRep,mType> >
38  {
39  public:
41  typedef BFSet BasisFunctionType;
45  static const unsigned int max_nchildren = 26;
46 
47  using ParentType::Instance;
48 
50  static bool directional() { return ParentType::default_directional(); }
51 
52  private:
53  using ParentType::RecurrenceRelation::expr_;
54  using ParentType::RecurrenceRelation::nflops_;
55  using ParentType::target_;
56  using ParentType::is_simple;
57 
59  VRR_11_TwoPRep_11(const SafePtr<TargetType>&, unsigned int dir);
60 
61  static std::string descr() { return "OSVRR"; }
65  std::string generate_label() const
66  {
67  typedef typename TargetType::AuxIndexType mType;
68  static SafePtr<mType> aux0(new mType(0u));
69  ostringstream os;
70  os << descr() << "P" << part << to_string(where)
71  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
72  return os.str();
73  }
74 
75  #if LIBINT_ENABLE_GENERIC_CODE
76  bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
79  std::string generic_header() const;
81  std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
82  #endif
83  };
84 
85  template <class F, int part, FunctionPosition where>
86  VRR_11_TwoPRep_11<F,part,where>::VRR_11_TwoPRep_11(const SafePtr< TargetType >& Tint,
87  unsigned int dir) :
88  ParentType(Tint,dir)
89  {
90  using namespace libint2::algebra;
91  using namespace libint2::prefactor;
92  using namespace libint2::braket;
93  const unsigned int m = Tint->aux()->elem(0);
94  const F& _1 = unit<F>(dir);
95 
96  { // can't apply to contracted basis functions
97  F a(Tint->bra(0,0));
98  F b(Tint->ket(0,0));
99  F c(Tint->bra(1,0));
100  F d(Tint->ket(1,0));
101  if (a.contracted() ||
102  b.contracted() ||
103  c.contracted() ||
104  d.contracted())
105  return;
106  }
107 
108  // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
109  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
110  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
111  const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
112  const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
113  const bool deriv = dA.zero() == false ||
114  dB.zero() == false ||
115  dC.zero() == false ||
116  dD.zero() == false;
117 
118  // This is a hack to avoid creating recurrence relations for which generic code has not been yet implemented
119 #if LIBINT_ENABLE_GENERIC_CODE
120  {
121  F sh_a(target_->bra(0,0));
122  F sh_b(target_->ket(0,0));
123  F sh_c(target_->bra(1,0));
124  F sh_d(target_->ket(1,0));
125  // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
126  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
127  if (sh_b.zero() && sh_d.zero() &&
128  (sh_a.norm() > 1u && sh_c.norm() > 1u)
129  ) { // have a generic implemented ...
130  if (part != 0) // ... but only implemented build on A in this case
131  return;
132  }
133  if (sh_a.zero() && sh_c.zero() &&
134  (sh_b.norm() > 1u && sh_d.norm() > 1u)
135  ) {
136  if (part != 0) // but only implemented build on B in this case
137  return;
138  }
139  }
140 #endif
141 
142  typedef TargetType ChildType;
143  ChildFactory<ThisType,ChildType> factory(this);
144 
145  bool part0_has_unit=false, part1_has_unit=false;
146 
147  // Build on A
148  if (part == 0 && where == InBra) {
149  F a(Tint->bra(0,0) - _1);
150  if (!exists(a)) return;
151  F b(Tint->ket(0,0)); const bool unit_b = (b == F::unit()); part0_has_unit |= unit_b;
152  F c(Tint->bra(1,0));
153  F d(Tint->ket(1,0));
154 
155  SafePtr<DGVertex> ABCD_m; if (not unit_b) ABCD_m = factory.make_child(a,b,c,d,m);
156  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
157  if (is_simple()) {
158  if (not unit_b) {
159  expr_ = Vector("PA")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1; nflops_+=3;
160  }
161  else {
162  expr_ = Vector("WP")[dir] * ABCD_mp1; nflops_+=1;
163  }
164  }
165 
166  // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
167  const bool ahlrichs_simplification = a.pure_sh() && unit_b;
168  auto am1 = a - _1;
169  if (exists(am1) && not ahlrichs_simplification) {
170  auto Am1BCD_m = factory.make_child(am1,b,c,d,m);
171  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
172 #if LIBINT_GENERATE_FMA
173  // this form is amenable to generation of fmsub
174  if (is_simple()) { expr_ -= Scalar(a[dir]) * Scalar("oo2z") * (Scalar("roz") * Am1BCD_mp1 - Am1BCD_m); nflops_+=5; }
175 #else
176  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1); nflops_+=5; }
177 #endif
178  }
179  const F& bm1 = b - _1;
180  if (exists(bm1)) {
181  auto ABm1CD_m = factory.make_child(a,bm1,c,d,m);
182  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
183 #if LIBINT_GENERATE_FMA
184  // this form is amenable to generation of fmsub
185  if (is_simple()) { expr_ -= Scalar(b[dir]) * Scalar("oo2z") * (Scalar("roz") * ABm1CD_mp1 - ABm1CD_m); nflops_+=5; }
186 #else
187  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1); nflops_+=5; }
188 #endif
189  }
190  const F& cm1 = c - _1;
191  if (exists(cm1)) {
192  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
193  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1; nflops_+=3; }
194  }
195  const F& dm1 = d - _1;
196  if (exists(dm1)) {
197  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
198  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1; nflops_+=3; }
199  }
200  }
201  // Build on B
202  if (part == 0 && where == InKet) {
203  F a(Tint->bra(0,0)); const bool unit_a = (a == F::unit()); part0_has_unit |= unit_a;
204  F b(Tint->ket(0,0) - _1);
205  if (!exists(b)) return;
206  F c(Tint->bra(1,0));
207  F d(Tint->ket(1,0));
208 
209  SafePtr<DGVertex> ABCD_m; if (not unit_a) ABCD_m = factory.make_child(a,b,c,d,m);
210  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
211  if (is_simple()) {
212  if (not unit_a) {
213  expr_ = Vector("PB")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1; nflops_+=3;
214  }
215  else {
216  expr_ = Vector("WP")[dir] * ABCD_mp1; nflops_+=1;
217  }
218  }
219 
220  const F& am1 = a - _1;
221  if (exists(am1)) {
222  auto Am1BCD_m = factory.make_child(am1,b,c,d,m);
223  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
224 #if LIBINT_GENERATE_FMA
225  // this form is amenable to generation of fmsub
226  if (is_simple()) { expr_ -= Scalar(a[dir]) * Scalar("oo2z") * (Scalar("roz") * Am1BCD_mp1 - Am1BCD_m); nflops_+=5; }
227 #else
228  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1); nflops_+=5; }
229 #endif
230  }
231  // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
232  const bool ahlrichs_simplification = b.pure_sh() && unit_a;
233  const F& bm1 = b - _1;
234  if (exists(bm1) && not ahlrichs_simplification) {
235  auto ABm1CD_m = factory.make_child(a,bm1,c,d,m);
236  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
237 #if LIBINT_GENERATE_FMA
238  // this form is amenable to generation of fmsub
239  if (is_simple()) { expr_ -= Scalar(b[dir]) * Scalar("oo2z") * (Scalar("roz") * ABm1CD_mp1 - ABm1CD_m); nflops_+=5; }
240 #else
241  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1); nflops_+=5; }
242 #endif
243  }
244  const F& cm1 = c - _1;
245  if (exists(cm1)) {
246  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
247  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1; nflops_+=3; }
248  }
249  const F& dm1 = d - _1;
250  if (exists(dm1)) {
251  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
252  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1; nflops_+=3; }
253  }
254  }
255  // Build on C
256  if (part == 1 && where == InBra) {
257  F a(Tint->bra(0,0));
258  F b(Tint->ket(0,0));
259  F c(Tint->bra(1,0) - _1);
260  if (!exists(c)) return;
261  F d(Tint->ket(1,0)); const bool unit_d = (d == F::unit()); part1_has_unit |= unit_d;
262 
263  SafePtr<DGVertex> ABCD_m; if (not unit_d) ABCD_m = factory.make_child(a,b,c,d,m);
264  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
265  if (is_simple()) {
266  if (not unit_d) {
267  expr_ = Vector("QC")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1; nflops_+=3;
268  }
269  else {
270  expr_ = Vector("WQ")[dir] * ABCD_mp1; nflops_+=1;
271  }
272  }
273 
274  // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
275  const bool ahlrichs_simplification = c.pure_sh() && unit_d;
276  const F& cm1 = c - _1;
277  if (exists(cm1) && not ahlrichs_simplification) {
278  auto ABCm1D_m = factory.make_child(a,b,cm1,d,m);
279  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
280 #if LIBINT_GENERATE_FMA
281  // this form is amenable to generation of fmsub
282  if (is_simple()) { expr_ -= Scalar(c[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCm1D_mp1 - ABCm1D_m); nflops_+=5; }
283 #else
284  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1); nflops_+=5; }
285 #endif
286  }
287  const F& dm1 = d - _1;
288  if (exists(dm1)) {
289  auto ABCDm1_m = factory.make_child(a,b,c,dm1,m);
290  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
291 #if LIBINT_GENERATE_FMA
292  // this form is amenable to generation of fmsub
293  if (is_simple()) { expr_ -= Scalar(d[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCDm1_mp1 - ABCDm1_m); nflops_+=5; }
294 #else
295  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1); nflops_+=5; }
296 #endif
297  }
298  const F& am1 = a - _1;
299  if (exists(am1)) {
300  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
301  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1; nflops_+=3; }
302  }
303  const F& bm1 = b - _1;
304  if (exists(bm1)) {
305  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
306  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1; nflops_+=3; }
307  }
308  }
309  // Build on D
310  if (part == 1 && where == InKet) {
311  F a(Tint->bra(0,0));
312  F b(Tint->ket(0,0));
313  F c(Tint->bra(1,0)); const bool unit_c = (c == F::unit()); part1_has_unit |= unit_c;
314  F d(Tint->ket(1,0) - _1);
315  if (!exists(d)) return;
316 
317  SafePtr<DGVertex> ABCD_m; if (not unit_c) ABCD_m = factory.make_child(a,b,c,d,m);
318  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
319  if (is_simple()) {
320  if (not unit_c) {
321  expr_ = Vector("QD")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1; nflops_+=3;
322  }
323  else {
324  expr_ = Vector("WQ")[dir] * ABCD_mp1; nflops_+=1;
325  }
326  }
327 
328  const F& cm1 = c - _1;
329  if (exists(cm1)) {
330  auto ABCm1D_m = factory.make_child(a,b,cm1,d,m);
331  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
332 #if LIBINT_GENERATE_FMA
333  // this form is amenable to generation of fmsub
334  if (is_simple()) { expr_ -= Scalar(c[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCm1D_mp1 - ABCm1D_m); nflops_+=5; }
335 #else
336  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1); nflops_+=5; }
337 #endif
338  }
339  // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
340  const bool ahlrichs_simplification = d.pure_sh() && unit_c;
341  const F& dm1 = d - _1;
342  if (exists(dm1) && not ahlrichs_simplification) {
343  auto ABCDm1_m = factory.make_child(a,b,c,dm1,m);
344  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
345 #if LIBINT_GENERATE_FMA
346  // this form is amenable to generation of fmsub
347  if (is_simple()) { expr_ -= Scalar(d[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCDm1_mp1 - ABCDm1_m); nflops_+=5; }
348 #else
349  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1); nflops_+=5; }
350 #endif
351  }
352  const F& am1 = a - _1;
353  if (exists(am1)) {
354  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
355  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1; nflops_+=3; }
356  }
357  const F& bm1 = b - _1;
358  if (exists(bm1)) {
359  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
360  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1; nflops_+=3; }
361  }
362  }
363 
364  // if got here, can decrement by at least 1 quantum
365  // add additional derivative terms
366  if (deriv) {
367  // bf quantum on the build center subtracted by 1
368  F a( part == 0 && where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
369  F b( part == 0 && where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
370  F c( part == 1 && where == InBra ? Tint->bra(1,0) - _1 : Tint->bra(1,0) );
371  F d( part == 1 && where == InKet ? Tint->ket(1,0) - _1 : Tint->ket(1,0) );
372 
373  // treatment of derivative terms differs for shell sets and integrals
374  // since in computing shell sets transfer/build will occur in all 3 directions
375  // change in up to all three derivative indices will occur
376  for(unsigned int dxyz=0; dxyz<3; ++dxyz) {
377 
378  if (is_simple() && dxyz != dir) // for integrals only consider derivatives in THE build direction
379  continue;
380 
381  OriginDerivative<3u> _d1; _d1.inc(dxyz);
382 
383  SafePtr<DGVertex> _nullptr;
384 
385  // dA - _1?
386  {
387  const OriginDerivative<3u> dAm1(dA - _d1);
388  if (exists(dAm1)) { // yes
389  a.deriv() = dAm1;
390  auto ABCD_m = (part == 0 && not part0_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
391  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
392  if (is_simple()) {
393  if (part == 0 && where == InBra) { // building on A
394  if (not part0_has_unit) {
395  expr_ -= Vector(dA)[dxyz] * (Scalar("rho12_over_alpha1") * ABCD_m + Scalar("alpha1_rho_over_zeta2") * ABCD_mp1); nflops_ += 5;
396  }
397  else {
398  expr_ -= Vector(dA)[dxyz] * (Scalar("alpha1_rho_over_zeta2") * ABCD_mp1); nflops_ += 3;
399  }
400  }
401  if (part == 0 && where == InKet) { // building on B
402  if (not part0_has_unit) {
403  expr_ += Vector(dA)[dxyz] * (Scalar("rho12_over_alpha2") * ABCD_m - Scalar("alpha1_rho_over_zeta2") * ABCD_mp1); nflops_ += 5;
404  }
405  else {
406  expr_ -= Vector(dA)[dxyz] * (Scalar("alpha1_rho_over_zeta2") * ABCD_mp1); nflops_ += 3;
407  }
408  }
409  if (part == 1) { // building on C or D
410  expr_ += Vector(dA)[dxyz] * Scalar("alpha1_over_zetapluseta") * ABCD_mp1; nflops_ += 3; }
411  }
412  a.deriv() = dA;
413  }
414  }
415 
416  // dB - _1?
417  {
418  const OriginDerivative<3u> dBm1(dB - _d1);
419  if (exists(dBm1)) { // yes
420  b.deriv() = dBm1;
421  auto ABCD_m = (part == 0 && not part0_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
422  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
423  if (is_simple()) {
424  if (part == 0 && where == InBra) { // building on A
425  if (not part0_has_unit) {
426  expr_ += Vector(dB)[dxyz] * (Scalar("rho12_over_alpha1") * ABCD_m - Scalar("alpha2_rho_over_zeta2") * ABCD_mp1); nflops_ += 5;
427  }
428  else {
429  expr_ -= Vector(dB)[dxyz] * (Scalar("alpha2_rho_over_zeta2") * ABCD_mp1); nflops_ += 3;
430  }
431  }
432  if (part == 0 && where == InKet) { // building on B
433  if (not part0_has_unit) {
434  expr_ -= Vector(dB)[dxyz] * (Scalar("rho12_over_alpha2") * ABCD_m + Scalar("alpha2_rho_over_zeta2") * ABCD_mp1); nflops_ += 5;
435  }
436  else {
437  expr_ -= Vector(dB)[dxyz] * (Scalar("alpha2_rho_over_zeta2") * ABCD_mp1); nflops_ += 3;
438  }
439  }
440  if (part == 1) { // building on C or D
441  expr_ += Vector(dB)[dxyz] * Scalar("alpha2_over_zetapluseta") * ABCD_mp1; nflops_ += 3; }
442  }
443  b.deriv() = dB;
444  }
445  }
446 
447  // dC - _1?
448  {
449  const OriginDerivative<3u> dCm1(dC - _d1);
450  if (exists(dCm1)) { // yes
451  c.deriv() = dCm1;
452  auto ABCD_m = (part == 1 && not part1_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
453  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
454  if (is_simple()) {
455  if (part == 0) { // building on A or B
456  expr_ += Vector(dC)[dxyz] * Scalar("alpha3_over_zetapluseta") * ABCD_mp1; nflops_ += 3; }
457  if (part == 1 && where == InBra) { // building on C
458  if (not part1_has_unit) {
459  expr_ -= Vector(dC)[dxyz] * (Scalar("rho34_over_alpha3") * ABCD_m + Scalar("alpha3_rho_over_eta2") * ABCD_mp1); nflops_ += 5;
460  }
461  else {
462  expr_ -= Vector(dC)[dxyz] * (Scalar("alpha3_rho_over_eta2") * ABCD_mp1); nflops_ += 3;
463  }
464  }
465  if (part == 1 && where == InKet) { // building on D
466  if (not part1_has_unit) {
467  expr_ += Vector(dC)[dxyz] * (Scalar("rho34_over_alpha4") * ABCD_m - Scalar("alpha3_rho_over_eta2") * ABCD_mp1); nflops_ += 5;
468  }
469  else {
470  expr_ -= Vector(dC)[dxyz] * (Scalar("alpha3_rho_over_eta2") * ABCD_mp1); nflops_ += 3;
471  }
472  }
473  }
474  c.deriv() = dC;
475  }
476  }
477 
478  // dD - _1?
479  {
480  const OriginDerivative<3u> dDm1(dD - _d1);
481  if (exists(dDm1)) { // yes
482  d.deriv() = dDm1;
483  auto ABCD_m = (part == 1 && not part1_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
484  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
485  if (is_simple()) {
486  if (part == 0) { // building on A or B
487  expr_ += Vector(dD)[dxyz] * Scalar("alpha4_over_zetapluseta") * ABCD_mp1; nflops_ += 3; }
488  if (part == 1 && where == InBra) { // building on C
489  if (not part1_has_unit) {
490  expr_ += Vector(dD)[dxyz] * (Scalar("rho34_over_alpha3") * ABCD_m - Scalar("alpha4_rho_over_eta2") * ABCD_mp1); nflops_ += 5;
491  }
492  else {
493  expr_ -= Vector(dD)[dxyz] * (Scalar("alpha4_rho_over_eta2") * ABCD_mp1); nflops_ += 3;
494  }
495  }
496  if (part == 1 && where == InKet) { // building on D
497  if (not part1_has_unit) {
498  expr_ -= Vector(dD)[dxyz] * (Scalar("rho34_over_alpha4") * ABCD_m + Scalar("alpha4_rho_over_eta2") * ABCD_mp1); nflops_ += 5;
499  }
500  else {
501  expr_ -= Vector(dD)[dxyz] * (Scalar("alpha4_rho_over_eta2") * ABCD_mp1); nflops_ += 3;
502  }
503  }
504  }
505  d.deriv() = dD;
506  }
507  }
508  }
509  } // end of deriv
510 
511  return;
512  }
513 
514 #if LIBINT_ENABLE_GENERIC_CODE
515  template <class F, int part, FunctionPosition where>
516  bool
517  VRR_11_TwoPRep_11<F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
518  {
519  if (TrivialBFSet<F>::result)
520  return false;
521 
522  F sh_a(target_->bra(0,0));
523  F sh_b(target_->ket(0,0));
524  F sh_c(target_->bra(1,0));
525  F sh_d(target_->ket(1,0));
526  const unsigned int max_opt_am = cparams->max_am_opt();
527  // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
528  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
529  if (sh_b.zero() && sh_d.zero() &&
530  (sh_a.norm() > std::max(2*max_opt_am,1u) ||
531  sh_c.norm() > std::max(2*max_opt_am,1u)
532  ) &&
533  (sh_a.norm() > 1u && sh_c.norm() > 1u)
534  ) {
535  assert(part == 0); // has only implemented build on A in this case
536  return true;
537  }
538  if (sh_a.zero() && sh_c.zero() &&
539  (sh_b.norm() > std::max(2*max_opt_am,1u) ||
540  sh_d.norm() > std::max(2*max_opt_am,1u)
541  ) &&
542  (sh_b.norm() > 1u && sh_d.norm() > 1u)
543  ) {
544  assert(part == 0); // has only implemented build on B in this case
545  return true;
546  }
547  return false;
548  }
549 
550  template <class F, int part, FunctionPosition where>
551  std::string
552  VRR_11_TwoPRep_11<F,part,where>::generic_header() const
553  {
554  F sh_a(target_->bra(0,0));
555  F sh_b(target_->ket(0,0));
556  F sh_c(target_->bra(1,0));
557  F sh_d(target_->ket(1,0));
558  const bool xsxs = sh_b.zero() && sh_d.zero();
559  const bool sxsx = sh_a.zero() && sh_c.zero();
560 
561  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
562  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
563  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
564  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
565  const bool deriv = dA.zero() == false ||
566  dB.zero() == false ||
567  dC.zero() == false ||
568  dD.zero() == false;
569 
570  if (deriv == false) {
571  if (xsxs) return std::string("OSVRR_xs_xs.h");
572  if (sxsx) return std::string("OSVRR_sx_sx.h");
573  }
574  else {
575  if (xsxs) return std::string("OSVRR_xs_xs_deriv.h");
576  if (sxsx) return std::string("OSVRR_sx_sx_deriv.h");
577  }
578  abort(); // unreachable
579  }
580 
581  template <class F, int part, FunctionPosition where>
582  std::string
583  VRR_11_TwoPRep_11<F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
584  {
585  std::ostringstream oss;
586  F sh_a(target_->bra(0,0));
587  F sh_b(target_->ket(0,0));
588  F sh_c(target_->bra(1,0));
589  F sh_d(target_->ket(1,0));
590  const bool xsxs = sh_b.zero() && sh_d.zero();
591  const bool sxsx = sh_a.zero() && sh_c.zero();
592  bool ahlrichs_simplification = false;
593  bool unit_s = false;
594  bool part0_has_unit = false;
595  bool part1_has_unit = false;
596  if (xsxs) {
597  ahlrichs_simplification = (sh_a.pure_sh() && sh_b.is_unit()) ||
598  (sh_c.pure_sh() && sh_d.is_unit());
599  unit_s = sh_b.is_unit();
600  part0_has_unit = sh_b.is_unit();
601  part1_has_unit = sh_d.is_unit();
602  }
603  if (sxsx) {
604  ahlrichs_simplification = (sh_b.pure_sh() && sh_a.is_unit()) ||
605  (sh_d.pure_sh() && sh_c.is_unit());
606  unit_s = sh_a.is_unit();
607  part0_has_unit = sh_a.is_unit();
608  part1_has_unit = sh_c.is_unit();
609  }
610 
611  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
612  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
613  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
614  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
615  const bool deriv = dA.zero() == false ||
616  dB.zero() == false ||
617  dC.zero() == false ||
618  dD.zero() == false;
619 
620  oss << "using namespace libint2;" << endl;
621 
622  if (deriv == false) { // for regular integrals I know exactly how many prerequisites I need
623  if(xsxs) {
624  oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
625  << "VRR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
626  }
627  if (sxsx) {
628  oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
629  << "VRR_sx_sx<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
630  }
631  if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
632  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
633  oss << ">::compute(inteval";
634 
635  oss << "," << args->symbol(0); // target
636  if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
637  oss << ",0"; // src0-> 0x0
638  const unsigned int nargs = args->n();
639  for(unsigned int a=1; a<nargs; a++) {
640  oss << "," << args->symbol(a);
641  }
642  oss << ");";
643  }
644  else { // deriv == true -> only some arguments are needed
645  if(xsxs) {
646  oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
647  << "VRR_xs_xs_deriv<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
648  }
649  if(sxsx) {
650  oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
651  << "VRR_sx_sx_deriv<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
652  }
653 
654  for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_a.deriv().d(xyz) << ",";
655  for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_b.deriv().d(xyz) << ",";
656  for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_c.deriv().d(xyz) << ",";
657  for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_d.deriv().d(xyz) << ",";
658  if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
659  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
660  oss << ">::compute(inteval";
661  // out of all 22 possible prerequisites first few have same derivative degree as the target
662  // 5 if standard 4-center integral
663  // 1+4 if the center opposite the build center carries a unit function
664  // will pass 0 instead of src0
665  // 2 if Ahlrichs
666  oss << "," << args->symbol(0); // target
667  if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
668  oss << ",0"; // src0-> 0x0
669  //const unsigned int nargs = args->n();
670  unsigned int arg = 1;
671  for(;
672  arg<(ahlrichs_simplification? 3 : (unit_s ? 5 : 6)); // nargs + 1 target
673  arg++) {
674  oss << "," << args->symbol(arg);
675  }
676  for(unsigned int xyz=0; xyz<3; ++xyz) {
677  if (sh_a.deriv().d(xyz) > 0) {
678  // see the dA-1 clause: (ab|cd)^m is skipped
679  oss << "," << (part0_has_unit ? std::to_string(0) : args->symbol(arg++));
680  oss << "," << args->symbol(arg++);
681  }
682  else
683  oss << ",0,0";
684  if (sh_b.deriv().d(xyz) > 0) {
685  // see the dB-1 clause: (ab|cd)^m is skipped
686  oss << "," << (part0_has_unit ? std::to_string(0) : args->symbol(arg++));
687  oss << "," << args->symbol(arg++);
688  }
689  else
690  oss << ",0,0";
691  if (sh_c.deriv().d(xyz) > 0) {
692  oss << "," << args->symbol(arg++);
693  }
694  else
695  oss << ",0";
696  if (sh_d.deriv().d(xyz) > 0) {
697  oss << "," << args->symbol(arg++);
698  }
699  else
700  oss << ",0";
701  }
702  oss << ");";
703  }
704 
705  return oss.str();
706  }
707 #endif // #if !LIBINT_ENABLE_GENERIC_CODE
708 
709 };
710 
711 #endif
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: generic_rr.h:81
Generic integral over a two-body operator with one bfs for each particle in bra and ket.
Definition: integral_11_11.h:33
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:49
VRR Recurrence Relation for 2-e ERI.
Definition: vrr_11_twoprep_11.h:35
these objects help to construct BraketPairs
Definition: braket.h:270
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
Set of basis functions.
Definition: bfset.h:42
static bool directional()
Default directionality.
Definition: vrr_11_twoprep_11.h:50