LIBINT  2.6.0
vrr_1_onep_1.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_vrr1onep1_h_
22 #define _libint2_src_bin_libint_vrr1onep1_h_
23 
24 #include <cassert>
25 #include <generic_rr.h>
26 #include <onep_1_1.h>
27 
28 using namespace std;
29 
30 namespace libint2 {
31 
35  template <class BFSet, FunctionPosition where>
36  class VRR_1_Overlap_1 : public GenericRecurrenceRelation< VRR_1_Overlap_1<BFSet,where>,
37  BFSet,
38  GenIntegralSet_1_1<BFSet,OverlapOper,EmptySet> >
39  {
40  public:
41  typedef VRR_1_Overlap_1 ThisType;
42  typedef BFSet BasisFunctionType;
46  static const unsigned int max_nchildren = 9;
47 
48  using ParentType::Instance;
49 
51  static bool directional() { return ParentType::default_directional(); }
52 
53  private:
54  using ParentType::RecurrenceRelation::expr_;
55  using ParentType::RecurrenceRelation::nflops_;
56  using ParentType::target_;
57  using ParentType::is_simple;
58 
60  VRR_1_Overlap_1(const SafePtr<TargetType>&, unsigned int dir);
61 
62  static std::string descr() { return "OSVRROverlap"; }
63  };
64 
65  template <class F, FunctionPosition where>
66  VRR_1_Overlap_1<F,where>::VRR_1_Overlap_1(const SafePtr< TargetType >& Tint,
67  unsigned int dir) :
68  ParentType(Tint,dir)
69  {
70  using namespace libint2::algebra;
71  using namespace libint2::prefactor;
72  using namespace libint2::braket;
73  const F& _1 = unit<F>(dir);
74 
75  { // can't apply to contracted basis functions
76  F a(Tint->bra(0,0));
77  F b(Tint->ket(0,0));
78  if (a.contracted() ||
79  b.contracted())
80  return;
81  }
82 
83  // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
84  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
85  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
86  const bool deriv = dA.zero() == false ||
87  dB.zero() == false;
88 
89  typedef TargetType ChildType;
90  ChildFactory<ThisType,ChildType> factory(this);
91 
92  // Build on A or B
93  {
94  // bf quantum on the build center subtracted by 1
95  auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
96  if (!exists(a)) return;
97  auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
98  if (!exists(b)) return;
99 
100  auto AB = factory.make_child(a,b);
101  if (is_simple()) { expr_ = Vector(where == InBra ? "PA" : "PB")[dir] * AB; nflops_+=1; }
102 
103  auto am1 = a - _1; auto am1_exists = exists(am1);
104  auto bm1 = b - _1; auto bm1_exists = exists(bm1);
105 
106  if (am1_exists) {
107  auto Am1B = factory.make_child(am1,b);
108  if (is_simple()) { expr_ += (Scalar(a[dir]) * Scalar("oo2z")) * Am1B; nflops_+=3; }
109  }
110  if (bm1_exists) {
111  auto ABm1 = factory.make_child(a,bm1);
112  if (is_simple()) { expr_ += (Scalar(b[dir]) * Scalar("oo2z")) * ABm1; nflops_+=3; }
113  }
114  }
115 
116  // if got here, can decrement by at least 1 quantum
117  // add additional derivative terms
118  if (deriv) {
119  // bf quantum on the build center subtracted by 1
120  F a( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
121  F b( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
122 
123  // treatment of derivative terms differs for shell sets and integrals
124  // since in computing shell sets transfer/build will occur in all 3 directions
125  // change in up to all three derivative indices will occur
126  for(unsigned int dxyz=0; dxyz<3; ++dxyz) {
127 
128  if (is_simple() && dxyz != dir) // for integrals only consider derivatives in THE build direction
129  continue;
130 
131  OriginDerivative<3u> _d1; _d1.inc(dxyz);
132 
133  SafePtr<DGVertex> _nullptr;
134 
135  // dA - _1?
136  {
137  const OriginDerivative<3u> dAm1(dA - _d1);
138  if (exists(dAm1)) { // yes
139  a.deriv() = dAm1;
140  auto AB = factory.make_child(a,b);
141  if (is_simple()) {
142  if (where == InBra) { // building on A
143  expr_ -= Vector(dA)[dxyz] * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
144  if (where == InKet) { // building on B
145  expr_ += Vector(dA)[dxyz] * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
146  }
147  a.deriv() = dA;
148  }
149  }
150 
151  // dB - _1?
152  {
153  const OriginDerivative<3u> dBm1(dB - _d1);
154  if (exists(dBm1)) { // yes
155  b.deriv() = dBm1;
156  auto AB = factory.make_child(a,b);
157  if (is_simple()) {
158  if (where == InBra) { // building on A
159  expr_ += Vector(dB)[dxyz] * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
160  if (where == InKet) { // building on B
161  expr_ -= Vector(dB)[dxyz] * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
162  }
163  b.deriv() = dB;
164  }
165  }
166 
167  }
168  } // end of deriv
169 
170  return;
171  }
172 
176  template <CartesianAxis Axis, FunctionPosition where>
177  class VRR_1_Overlap_1_1d : public GenericRecurrenceRelation< VRR_1_Overlap_1_1d<Axis,where>,
178  CGF1d<Axis>,
179  GenIntegralSet_1_1<CGF1d<Axis>,OverlapOper,EmptySet> >
180  {
181  public:
187  static const unsigned int max_nchildren = 9;
188  static constexpr auto axis = Axis;
189 
190  using ParentType::Instance;
191 
193  static bool directional() { return ParentType::default_directional(); }
194 
195  private:
196  using ParentType::RecurrenceRelation::expr_;
197  using ParentType::RecurrenceRelation::nflops_;
198  using ParentType::target_;
199  using ParentType::is_simple;
200 
202  VRR_1_Overlap_1_1d(const SafePtr<TargetType>&, unsigned int dir);
203 
204  static std::string descr() { return std::string("OSVRROverlap") + to_string(axis); }
205  };
206 
207  template <CartesianAxis Axis, FunctionPosition where>
208  VRR_1_Overlap_1_1d<Axis,where>::VRR_1_Overlap_1_1d(const SafePtr< TargetType >& Tint,
209  unsigned int dir) :
210  ParentType(Tint,dir)
211  {
212  assert(dir == 0); // this integral is along 1 axis only
213 
214  using namespace libint2::algebra;
215  using namespace libint2::prefactor;
216  using namespace libint2::braket;
217  typedef CGF1d<Axis> F;
218  const F& _1 = unit<F>(dir);
219 
220  { // can't apply to contracted basis functions
221  F a(Tint->bra(0,0));
222  F b(Tint->ket(0,0));
223  if (a.contracted() ||
224  b.contracted())
225  return;
226  }
227 
228  // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
229  const OriginDerivative<1u> dA = Tint->bra(0,0).deriv();
230  const OriginDerivative<1u> dB = Tint->ket(0,0).deriv();
231  const bool deriv = dA.zero() == false ||
232  dB.zero() == false;
233 
234  typedef TargetType ChildType;
235  ChildFactory<ThisType,ChildType> factory(this);
236 
237  // handle the special case of (0|0) integral
238  // to avoid complications with non-precomputed shell blocks it's "computed
239  // by copying from inteval
240  auto zero = Tint->bra(0,0).zero() and Tint->ket(0,0).zero() and not deriv;
241  if (zero) {
242  SafePtr<DGVertex> int00 = Vector("_0_Overlap_0")[Axis];
243  expr_ = Scalar(0u) + int00;
244  this->add_child(int00);
245  return;
246  }
247 
248  // Build on A or B
249  {
250  // bf quantum on the build center subtracted by 1
251  auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
252  if (!exists(a)) return;
253  auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
254  if (!exists(b)) return;
255 
256  auto AB = factory.make_child(a,b);
257  if (is_simple()) { expr_ = Vector(where == InBra ? "PA" : "PB")[Axis] * AB; nflops_+=1; }
258 
259  auto am1 = a - _1; auto am1_exists = exists(am1);
260  auto bm1 = b - _1; auto bm1_exists = exists(bm1);
261 
262  if (am1_exists) {
263  auto Am1B = factory.make_child(am1,b);
264  if (is_simple()) { expr_ += (Scalar(a.qn()) * Scalar("oo2z")) * Am1B; nflops_+=3; }
265  }
266  if (bm1_exists) {
267  auto ABm1 = factory.make_child(a,bm1);
268  if (is_simple()) { expr_ += (Scalar(b.qn()) * Scalar("oo2z")) * ABm1; nflops_+=3; }
269  }
270  }
271 
272  // if got here, can decrement by at least 1 quantum
273  // add additional derivative terms
274  if (deriv) {
275  // bf quantum on the build center subtracted by 1
276  F a( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
277  F b( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
278 
279  {
280  OriginDerivative<1u> _d1; _d1.inc(0);
281 
282  SafePtr<DGVertex> _nullptr;
283 
284  // dA - _1?
285  {
286  const OriginDerivative<1u> dAm1(dA - _d1);
287  if (exists(dAm1)) { // yes
288  a.deriv() = dAm1;
289  auto AB = factory.make_child(a,b);
290  if (is_simple()) {
291  if (where == InBra) { // building on A
292  expr_ -= Scalar(dA[0]) * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
293  if (where == InKet) { // building on B
294  expr_ += Scalar(dA[0]) * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
295  }
296  a.deriv() = dA;
297  }
298  }
299 
300  // dB - _1?
301  {
302  const OriginDerivative<1u> dBm1(dB - _d1);
303  if (exists(dBm1)) { // yes
304  b.deriv() = dBm1;
305  auto AB = factory.make_child(a,b);
306  if (is_simple()) {
307  if (where == InBra) { // building on A
308  expr_ += Scalar(dB[0]) * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
309  if (where == InKet) { // building on B
310  expr_ -= Scalar(dB[0]) * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
311  }
312  b.deriv() = dB;
313  }
314  }
315 
316  }
317  } // end of deriv
318 
319  return;
320  }
321 
322 
326  template <class BFSet, FunctionPosition where>
327  class VRR_1_Kinetic_1 : public GenericRecurrenceRelation< VRR_1_Kinetic_1<BFSet,where>,
328  BFSet,
329  GenIntegralSet_1_1<BFSet,KineticOper,EmptySet> >
330  {
331  public:
332  typedef VRR_1_Kinetic_1 ThisType;
333  typedef BFSet BasisFunctionType;
337  static const unsigned int max_nchildren = 9;
338 
339  using ParentType::Instance;
340 
342  static bool directional() { return ParentType::default_directional(); }
343 
344  private:
345  using ParentType::RecurrenceRelation::expr_;
346  using ParentType::RecurrenceRelation::nflops_;
347  using ParentType::target_;
348  using ParentType::is_simple;
349 
351  VRR_1_Kinetic_1(const SafePtr<TargetType>&, unsigned int dir);
352 
353  static std::string descr() { return "OSVRRKinetic"; }
354  };
355 
356  template <class F, FunctionPosition where>
357  VRR_1_Kinetic_1<F,where>::VRR_1_Kinetic_1(const SafePtr< TargetType >& Tint,
358  unsigned int dir) :
359  ParentType(Tint,dir)
360  {
361  using namespace libint2::algebra;
362  using namespace libint2::prefactor;
363  using namespace libint2::braket;
364  const F& _1 = unit<F>(dir);
365 
366  { // can't apply to contracted basis functions
367  F a(Tint->bra(0,0));
368  F b(Tint->ket(0,0));
369  if (a.contracted() ||
370  b.contracted())
371  return;
372  }
373 
374  // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
375  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
376  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
377  const bool deriv = dA.zero() == false ||
378  dB.zero() == false;
379 
380  typedef TargetType Child1Type;
381  ChildFactory<ThisType,Child1Type> factory(this);
382  typedef GenIntegralSet_1_1<F,OverlapOper,EmptySet> Child2Type;
383  ChildFactory<ThisType,Child2Type> overlap_factory(this);
384 
385  // Build on A or B
386  {
387  // bf quantum on the build center subtracted by 1
388  auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
389  if (!exists(a)) return;
390  auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
391  if (!exists(b)) return;
392 
393  auto AB = factory.make_child(a,b);
394  if (is_simple()) { expr_ = Vector(where == InBra ? "PA" : "PB")[dir] * AB; nflops_+=1; }
395 
396  auto am1 = a - _1;
397  if (exists(am1)) {
398  auto Am1B = factory.make_child(am1,b);
399  auto S_Am1B = (where == InBra) ? overlap_factory.make_child(am1,b) : SafePtr<DGVertex>();
400  if (is_simple()) {
401  if (where == InBra) {
402  expr_ += Scalar(a[dir]) * ( Scalar("oo2z") * Am1B -
403  Scalar("rho12_over_alpha1") * S_Am1B );
404  nflops_+=5;
405  }
406  else {
407  expr_ += Scalar(a[dir]) * Scalar("oo2z") * Am1B;
408  nflops_+=3;
409  }
410  }
411  }
412  auto bm1 = b - _1;
413  if (exists(bm1)) {
414  auto ABm1 = factory.make_child(a,bm1);
415  auto S_ABm1 = (where == InKet) ? overlap_factory.make_child(a,bm1) : SafePtr<DGVertex>();
416  if (is_simple()) {
417  if (where == InKet) {
418  expr_ += Scalar(b[dir]) * ( Scalar("oo2z") * ABm1 -
419  Scalar("rho12_over_alpha2") * S_ABm1 );
420  nflops_+=5;
421  }
422  else {
423  expr_ += Scalar(b[dir]) * Scalar("oo2z") * ABm1;
424  nflops_+=3;
425  }
426  }
427  }
428 
429  {
430  auto S_AB_target = where == InBra ? overlap_factory.make_child(a + _1,b) : overlap_factory.make_child(a,b + _1);
431  if (is_simple()) {
432  expr_ += Scalar("two_rho12") * S_AB_target;
433  nflops_+=2;
434  }
435  }
436 
437  }
438 
439  // if got here, can decrement by at least 1 quantum
440  // add additional derivative terms
441  if (deriv) {
442  assert(false); // not yet implemented
443 
444  // bf quantum on the build center subtracted by 1
445  F a( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
446  F b( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
447 
448  // treatment of derivative terms differs for shell sets and integrals
449  // since in computing shell sets transfer/build will occur in all 3 directions
450  // change in up to all three derivative indices will occur
451  for(unsigned int dxyz=0; dxyz<3; ++dxyz) {
452 
453  if (is_simple() && dxyz != dir) // for integrals only consider derivatives in THE build direction
454  continue;
455 
456  OriginDerivative<3u> _d1; _d1.inc(dxyz);
457 
458  SafePtr<DGVertex> _nullptr;
459 
460  // dA - _1?
461  {
462  const OriginDerivative<3u> dAm1(dA - _d1);
463  if (exists(dAm1)) { // yes
464  a.deriv() = dAm1;
465  auto AB = factory.make_child(a,b);
466  if (is_simple()) {
467  if (where == InBra) { // building on A
468  expr_ -= Vector(dA)[dxyz] * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
469  if (where == InKet) { // building on B
470  expr_ += Vector(dA)[dxyz] * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
471  }
472  a.deriv() = dA;
473  }
474  }
475 
476  // dB - _1?
477  {
478  const OriginDerivative<3u> dBm1(dB - _d1);
479  if (exists(dBm1)) { // yes
480  b.deriv() = dBm1;
481  auto AB = factory.make_child(a,b);
482  if (is_simple()) {
483  if (where == InBra) { // building on A
484  expr_ += Vector(dB)[dxyz] * Scalar("rho12_over_alpha1") * AB; nflops_ += 3; }
485  if (where == InKet) { // building on B
486  expr_ -= Vector(dB)[dxyz] * Scalar("rho12_over_alpha2") * AB; nflops_ += 3; }
487  }
488  b.deriv() = dB;
489  }
490  }
491 
492  }
493  } // end of deriv
494 
495  return;
496  }
497 
501  template <class BFSet, FunctionPosition where>
502  class VRR_1_ElecPot_1 : public GenericRecurrenceRelation< VRR_1_ElecPot_1<BFSet,where>,
503  BFSet,
504  GenIntegralSet_1_1<BFSet,ElecPotOper,mType> >
505  {
506  public:
507  typedef VRR_1_ElecPot_1 ThisType;
508  typedef BFSet BasisFunctionType;
512  static const unsigned int max_nchildren = 9;
513 
514  using ParentType::Instance;
515 
517  static bool directional() { return ParentType::default_directional(); }
518 
519  private:
520  using ParentType::RecurrenceRelation::expr_;
521  using ParentType::RecurrenceRelation::nflops_;
522  using ParentType::target_;
523  using ParentType::is_simple;
524 
526  VRR_1_ElecPot_1(const SafePtr<TargetType>&, unsigned int dir);
527 
528  static std::string descr() { return "OSVRRElecPot"; }
532  std::string generate_label() const
533  {
534  typedef typename TargetType::AuxIndexType mType;
535  static SafePtr<mType> aux0(new mType(0u));
536  ostringstream os;
537  os << descr() << to_string(where)
538  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
539  return os.str();
540  }
541  };
542 
543  template <class F, FunctionPosition where>
544  VRR_1_ElecPot_1<F,where>::VRR_1_ElecPot_1(const SafePtr< TargetType >& Tint,
545  unsigned int dir) :
546  ParentType(Tint,dir)
547  {
548  using namespace libint2::algebra;
549  using namespace libint2::prefactor;
550  using namespace libint2::braket;
551  const unsigned int m = Tint->aux()->elem(0);
552  const F& _1 = unit<F>(dir);
553 
554  { // can't apply to contracted basis functions
555  F a(Tint->bra(0,0));
556  F b(Tint->ket(0,0));
557  if (a.contracted() ||
558  b.contracted())
559  return;
560  }
561 
562  // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
563  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
564  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
565  const bool deriv = dA.zero() == false ||
566  dB.zero() == false;
567 
568  typedef TargetType ChildType;
569  ChildFactory<ThisType,ChildType> factory(this);
570 
571  // Build on A or B
572  {
573  // bf quantum on the build center subtracted by 1
574  auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
575  if (!exists(a)) return;
576  auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
577  if (!exists(b)) return;
578 
579  auto AB_m = factory.make_child(a,b,m);
580  auto AB_mp1 = factory.make_child(a,b,m+1);
581  if (is_simple()) {
582  expr_ = Vector(where == InBra ? "PA" : "PB")[dir] * AB_m -
583  Vector("PC")[dir] * AB_mp1;
584  nflops_+=3;
585  }
586 
587  auto am1 = a - _1;
588  if (exists(am1)) {
589  auto Am1B_m = factory.make_child(am1,b,m);
590  auto Am1B_mp1 = factory.make_child(am1,b,m+1);
591  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1B_m - Am1B_mp1); nflops_+=4; }
592  }
593  auto bm1 = b - _1;
594  if (exists(bm1)) {
595  auto ABm1_m = factory.make_child(a,bm1,m);
596  auto ABm1_mp1 = factory.make_child(a,bm1,m+1);
597  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1_m - ABm1_mp1); nflops_+=4; }
598  }
599  }
600 
601  // if got here, can decrement by at least 1 quantum
602  // add additional derivative terms
603  if (deriv) {
604  // bf quantum on the build center subtracted by 1
605  F a( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
606  F b( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
607 
608  // treatment of derivative terms differs for shell sets and integrals
609  // since in computing shell sets transfer/build will occur in all 3 directions
610  // change in up to all three derivative indices will occur
611  for(unsigned int dxyz=0; dxyz<3; ++dxyz) {
612 
613  if (is_simple() && dxyz != dir) // for integrals only consider derivatives in THE build direction
614  continue;
615 
616  OriginDerivative<3u> _d1; _d1.inc(dxyz);
617 
618  SafePtr<DGVertex> _nullptr;
619 
620  // dA - _1?
621  {
622  const OriginDerivative<3u> dAm1(dA - _d1);
623  if (exists(dAm1)) { // yes
624  a.deriv() = dAm1;
625  auto AB_m = factory.make_child(a,b,m);
626  auto AB_mp1 = factory.make_child(a,b,m+1);
627  if (is_simple()) {
628  if (where == InBra) { // building on A -> derivative of (PA)_i and (PC)_i prefactors w.r.t A_i
629  expr_ -= Vector(dA)[dxyz] * (Scalar("rho12_over_alpha1") * AB_m + Scalar("rho12_over_alpha2") * AB_mp1); nflops_ += 5; }
630  if (where == InKet) { // building on B -> derivative of (PB)_i and (PC)_i prefactors w.r.t A_i
631  expr_ += Vector(dA)[dxyz] * Scalar("rho12_over_alpha2") * (AB_m - AB_mp1); nflops_ += 4; }
632  }
633  a.deriv() = dA;
634  }
635  }
636 
637  // dB - _1?
638  {
639  const OriginDerivative<3u> dBm1(dB - _d1);
640  if (exists(dBm1)) { // yes
641  b.deriv() = dBm1;
642  auto AB_m = factory.make_child(a,b,m);
643  auto AB_mp1 = factory.make_child(a,b,m+1);
644  if (is_simple()) {
645  if (where == InBra) { // building on A -> derivative of (PA)_i and (PC)_i prefactors w.r.t B_i
646  expr_ += Vector(dB)[dxyz] * Scalar("rho12_over_alpha1") * (AB_m - AB_mp1); nflops_ += 4; }
647  if (where == InKet) { // building on B -> derivative of (PB)_i and (PC)_i prefactors w.r.t B_i
648  expr_ -= Vector(dB)[dxyz] * (Scalar("rho12_over_alpha2") * AB_m + Scalar("rho12_over_alpha1") * AB_mp1); nflops_ += 5; }
649  }
650  b.deriv() = dB;
651  }
652  }
653 
654  }
655  } // end of deriv
656 
657  return;
658  }
659 
663  template <class BFSet, FunctionPosition where>
664  class VRR_1_SMultipole_1 : public GenericRecurrenceRelation< VRR_1_SMultipole_1<BFSet,where>,
665  BFSet,
666  GenIntegralSet_1_1<BFSet,SphericalMultipoleOper,EmptySet> >
667  {
668  public:
670  typedef BFSet BasisFunctionType;
675  static const unsigned int max_nchildren = 8;
676 
677  using ParentType::Instance;
678 
680  static bool directional() { return ParentType::default_directional(); }
681 
682  private:
683  using ParentType::RecurrenceRelation::expr_;
684  using ParentType::RecurrenceRelation::nflops_;
685  using ParentType::target_;
686  using ParentType::is_simple;
687  using typename ParentType::RecurrenceRelation::ExprType;
688 
689  static std::string descr() { return "OSVRRSMultipole"; }
690 
692  VRR_1_SMultipole_1(const SafePtr<TargetType>& Tint, unsigned int dir) :
693  ParentType(Tint,dir)
694  {
695  using namespace libint2::algebra;
696  using namespace libint2::prefactor;
697  using namespace libint2::braket;
698  using Sign = SphericalMultipoleQuanta::Sign;
699  using F = BFSet;
700  const F& _1 = unit<F>(dir);
701 
702  { // can't apply to contracted basis functions
703  F a(Tint->bra(0,0));
704  F b(Tint->ket(0,0));
705  if (a.contracted() ||
706  b.contracted())
707  return;
708  }
709 
710  // implementation follows J.M. Pérez-Jordá and W. Yang, J Chem Phys 107, 1218 (1997).
711  // Eqs. (58-61) for gaussian quanta, and
712  // J.M. Pérez-Jordá and W. Yang, J Chem Phys 104, 8003 (1996), Eqs. (23-27) for multipole quanta
713  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
714  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
715  const bool deriv = dA.zero() == false ||
716  dB.zero() == false;
717  if (deriv) // derivative relations not yet implemented
718  return;
719 
720  auto O_l_m = Tint->oper()->descr();
721  const auto l = O_l_m.l();
722  const auto m = O_l_m.m();
723  const auto sign = O_l_m.sign();
724 
725  typedef TargetType ChildType;
726  ChildFactory<ThisType,ChildType> factory(this);
727 
728  auto make_child = [&](const F& bra, const F& ket, const OperType::Descriptor& descr) -> SafePtr<DGVertex> {
729  return factory.make_child(bra, ket, EmptySet(), descr);
730  };
731 
732  // Build on A or B if either bra or ket has quanta, otherwise build multipole quanta
733  if (Tint->bra(0,0).norm() != 0 || Tint->ket(0,0).norm() != 0) {
734  // build A or B
735 
736  // bf quantum on the build center subtracted by 1
737  auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
738  if (!exists(a)) return;
739  auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
740  if (!exists(b)) return;
741 
742  //
743  // first 2 terms are common to Eqs. (58-60)
744  //
745  auto AB = make_child(a, b, O_l_m);
746  if (is_simple()) { expr_ = Vector(where == InBra ? "PA" : "PB")[dir] * AB; nflops_+=1; }
747 
748  auto am1 = a - _1; auto am1_exists = exists(am1);
749  auto bm1 = b - _1; auto bm1_exists = exists(bm1);
750 
751  SafePtr<ExprType> subexpr; // to be multiplied by 1/(2 zeta)
752 
753  if (am1_exists) {
754  auto Am1B = make_child(am1, b, O_l_m);
755  if (is_simple()) { subexpr += Scalar(a[dir]) * Am1B; nflops_+=1; }
756  }
757  if (bm1_exists) {
758  auto ABm1 = make_child(a, bm1, O_l_m);
759  if (is_simple()) { subexpr += Scalar(b[dir]) * ABm1; nflops_+=1; }
760  }
761 
762  // Eq. (58)
763  // (a|N_{l,m}|b) += 0.5 * ((a|N_{l-1,m+1}|b) - (a|N_{l-1,m-1}|b))
764  auto O_lm1_mp1 = SphericalMultipole_Descr(l-1, m+1, sign);
765  if (exists(O_lm1_mp1)) {
766  auto AB_O_lm1_mp1 = make_child(a, b, O_lm1_mp1);
767  if (is_simple() && dir == 0) {
768  subexpr += Scalar(O_lm1_mp1.phase() > 0 ? 0.5 : -0.5) * AB_O_lm1_mp1;
769  nflops_ += 1;
770  }
771  }
772 
773  auto O_lm1_mm1 = SphericalMultipole_Descr(l-1, m-1, sign);
774  if (exists(O_lm1_mm1)) {
775  auto AB_O_lm1_mm1 = make_child(a, b, O_lm1_mm1);
776  if (is_simple() && dir == 0) {
777  subexpr -= Scalar(O_lm1_mm1.phase() > 0 ? 0.5 : -0.5) * AB_O_lm1_mm1;
778  nflops_ += 1;
779  }
780  }
781 
782  // Eq. (59)
783  // (a|N^{+-}_{l,m}|b) {+-}= 0.5 * ((a|N^{-+}_{l-1,m+1}|b) + (a|N^{-+}_{l-1,m-1}|b))
784  const auto m_pfac = sign == Sign::plus ? 1 : -1;
785  auto Om_lm1_mp1 = SphericalMultipole_Descr(l-1, m+1, flip(sign));
786  if (exists(Om_lm1_mp1)) {
787  auto AB_Om_lm1_mp1 = make_child(a, b, Om_lm1_mp1);
788  if (is_simple() && dir == 1) {
789  subexpr += Scalar(m_pfac * (Om_lm1_mp1.phase() > 0 ? 0.5 : -0.5)) * AB_Om_lm1_mp1;
790  nflops_ += 1;
791  }
792  }
793 
794  auto Om_lm1_mm1 = SphericalMultipole_Descr(l-1, m-1, flip(sign));
795  if (exists(Om_lm1_mm1)) {
796  auto AB_Om_lm1_mm1 = make_child(a, b, Om_lm1_mm1);
797  if (is_simple() && dir == 1) {
798  subexpr += Scalar(m_pfac * (Om_lm1_mm1.phase() > 0 ? 0.5 : -0.5)) * AB_Om_lm1_mm1;
799  nflops_ += 1;
800  }
801  }
802 
803  // Eq. (60)
804  auto O_lm1_m = SphericalMultipole_Descr(l-1, m, sign);
805  if (exists(O_lm1_m)) {
806  auto AB_O_lm1_m = make_child(a, b, O_lm1_m);
807  if (is_simple() && dir == 2) {
808  subexpr += AB_O_lm1_m;
809  nflops_ += 1;
810  }
811  }
812 
813  if (is_simple()) {
814  if(subexpr)
815  expr_ += Scalar("oo2z") * subexpr;
816  }
817  }
818  else {
819  // build multipole quanta only if dir == 0
820  if (dir != 0) return;
821  auto a = Tint->bra(0,0);
822  auto b = Tint->ket(0,0);
823 
824  if (l == m) {
825  if (l == 0) { // Eq.
826  assert(sign == Sign::plus); // Eq. (24) should not be needed
827  // since N^-_{0,0} should just be
828  // omitted above
829  typedef GenIntegralSet_1_1<BFSet, CartesianMultipoleOper<3u>, EmptySet>
830  OverlapType;
831  ChildFactory<ThisType, OverlapType> overlap_factory(this);
832  auto S00 = overlap_factory.make_child(a, b);
833  if (is_simple()) {
834  expr_ = Scalar(1) * S00; // Eq. (25) and Eq. (61)
835  }
836  } else {
837  SafePtr<ExprType> subexpr; // to be multiplied by - 1/(2 m)
838 
839  // Eqs. (25-26)
840  // (0|N^+_{m,m}|0) = -(1/(2m)) ( x (0|N^+_{m-1,m-1}|0) - y
841  // (0|N^-_{m-1,m-1}|0) )
842  auto Op_lm1_mm1 =
843  SphericalMultipole_Descr(l - 1, m - 1, Sign::plus);
844  auto AB_Op_lm1_mm1 = make_child(a, b, Op_lm1_mm1);
845  if (is_simple()) {
846  subexpr = Scalar(sign == Sign::plus ? "PO_x" : "PO_y") *
847  AB_Op_lm1_mm1;
848  }
849  auto Om_lm1_mm1 =
850  SphericalMultipole_Descr(l - 1, m - 1, Sign::minus);
851  if (exists(Om_lm1_mm1)) {
852  auto AB_Om_lm1_mm1 = make_child(a, b, Om_lm1_mm1);
853  if (is_simple()) {
854  if (sign == Sign::plus)
855  subexpr -= Scalar("PO_y") * AB_Om_lm1_mm1;
856  else // sign == Sign::minus
857  subexpr += Scalar("PO_x") * AB_Om_lm1_mm1;
858  }
859  }
860  if (is_simple()) {
861  assert(subexpr);
862  expr_ = Scalar(-0.5 / m) * subexpr;
863  }
864  }
865  } else { // l != m, use Eq. 27
866  SafePtr<ExprType> subexpr; // to be multiplied by 1/((l+m)(l-m))
867 
868  auto O_lm1_m = SphericalMultipole_Descr(l - 1, m, sign);
869  if (exists(O_lm1_m)) {
870  auto AB_O_lm1_m = make_child(a, b, O_lm1_m);
871  if (is_simple())
872  subexpr = Scalar((2 * l - 1)) * Scalar("PO_z") * AB_O_lm1_m;
873  }
874 
875  auto O_lm2_m = SphericalMultipole_Descr(l - 2, m, sign);
876  if (exists(O_lm2_m)) {
877  auto AB_O_lm2_m = make_child(a, b, O_lm2_m);
878  if (is_simple())
879  subexpr -= Scalar("PO2") * AB_O_lm2_m;
880  }
881 
882  if (is_simple()) {
883  assert(subexpr);
884  expr_ = Scalar(1.0 / ((l + m) * (l - m))) * subexpr;
885  }
886  }
887  }
888 
889  return;
890  }
891 
892  }; // VRR_1_SMultipole_1
893 
894 }; // namespace libint2
895 
896 #endif
static bool directional()
Default directionality.
Definition: vrr_1_onep_1.h:193
static bool default_directional()
is this recurrence relation parameterized by a direction (x, y, or z).
Definition: generic_rr.h:101
Cartesian components of 3D CGF = 1D CGF.
Definition: bfset.h:438
static SafePtr< RRImpl > Instance(const SafePtr< TargetType > &Tint, unsigned int dir)
Return an instance if applicable, or a null pointer otherwise.
Definition: generic_rr.h:57
static bool directional()
Default directionality.
Definition: vrr_1_onep_1.h:517
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
GenOper is a single operator described by descriptor Descr.
Definition: oper.h:162
VRR Recurrence Relation for 1-e overlap integrals.
Definition: vrr_1_onep_1.h:36
VRR Recurrence Relation for 1-e electrostatic potential integrals.
Definition: vrr_1_onep_1.h:502
DefaultQuantumNumbers< unsigned int, 1 >::Result mType
mType is the type that describes the auxiliary index of standard 2-body repulsion integrals
Definition: quanta.h:411
VRR Recurrence Relation for 1-e spherical multipole moment aka regular solid harmonics integrals.
Definition: vrr_1_onep_1.h:664
static bool directional()
Default directionality.
Definition: vrr_1_onep_1.h:51
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:49
VRR Recurrence Relation for 1-d overlap integrals.
Definition: vrr_1_onep_1.h:177
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:74
these objects help to construct BraketPairs
Definition: braket.h:270
static bool directional()
Default directionality.
Definition: vrr_1_onep_1.h:680
Generic integral over a one-body operator with one bfs for each particle in bra and ket.
Definition: integral_1_1.h:33
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
const SafePtr< DGVertex > & add_child(const SafePtr< DGVertex > &child)
add child
Definition: generic_rr.h:110
Set of basis functions.
Definition: bfset.h:42
static bool directional()
Default directionality.
Definition: vrr_1_onep_1.h:342
SphericalMultipoleQuanta::Sign flip(const SphericalMultipoleQuanta::Sign &s)
plus <-> minus
Definition: multipole.h:240
std::string generate_label() const
Implementation of RecurrenceRelation::generate_label()
Definition: generic_rr.h:86
const SafePtr< DGVertex > & make_child(const typename RealChildType::BasisFunctionType &A, const typename RealChildType::BasisFunctionType &B, const typename RealChildType::BasisFunctionType &C, const typename RealChildType::BasisFunctionType &D, const typename RealChildType::AuxIndexType &aux=typename RealChildType::AuxIndexType(), const typename RealChildType::OperType &oper=typename RealChildType::OperType())
make_child should really looks something like this, but gcc 4.3.0 craps out TODO test is this works
Definition: generic_rr.h:128
VRR Recurrence Relation for 1-e kinetic energy integrals.
Definition: vrr_1_onep_1.h:327
DefaultQuantumNumbers< int, 0 >::Result EmptySet
EmptySet is the type that describes null set of auxiliary indices.
Definition: quanta.h:407