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