21 #ifndef _libint2_src_bin_libint_vrr1onep1_h_ 22 #define _libint2_src_bin_libint_vrr1onep1_h_ 25 #include <generic_rr.h> 35 template <
class BFSet, FunctionPosition where>
38 GenIntegralSet_1_1<BFSet,OverlapOper,EmptySet> >
46 static const unsigned int max_nchildren = 9;
48 using ParentType::Instance;
51 static bool directional() {
return ParentType::default_directional(); }
54 using ParentType::RecurrenceRelation::expr_;
55 using ParentType::RecurrenceRelation::nflops_;
56 using ParentType::target_;
57 using ParentType::is_simple;
62 static std::string descr() {
return "OSVRROverlap"; }
65 template <
class F, FunctionPosition where>
66 VRR_1_Overlap_1<F,where>::VRR_1_Overlap_1(
const SafePtr< TargetType >& Tint,
70 using namespace libint2::algebra;
71 using namespace libint2::prefactor;
73 const F& _1 = unit<F>(dir);
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 ||
89 typedef TargetType ChildType;
90 ChildFactory<ThisType,ChildType> factory(
this);
95 auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
97 auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
100 auto AB = factory.make_child(a,b);
101 if (
is_simple()) { expr_ = Vector(where == InBra ?
"PA" :
"PB")[dir] * AB; nflops_+=1; }
103 auto am1 = a - _1;
auto am1_exists =
exists(am1);
104 auto bm1 = b - _1;
auto bm1_exists =
exists(bm1);
107 auto Am1B = factory.make_child(am1,b);
108 if (
is_simple()) { expr_ += (Scalar(a[dir]) * Scalar(
"oo2z")) * Am1B; nflops_+=3; }
111 auto ABm1 = factory.make_child(a,bm1);
112 if (
is_simple()) { expr_ += (Scalar(b[dir]) * Scalar(
"oo2z")) * ABm1; nflops_+=3; }
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) );
126 for(
unsigned int dxyz=0; dxyz<3; ++dxyz) {
131 OriginDerivative<3u> _d1; _d1.inc(dxyz);
133 SafePtr<DGVertex> _nullptr;
137 const OriginDerivative<3u> dAm1(dA - _d1);
140 auto AB = factory.make_child(a,b);
142 if (where == InBra) {
143 expr_ -= Vector(dA)[dxyz] * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
144 if (where == InKet) {
145 expr_ += Vector(dA)[dxyz] * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
153 const OriginDerivative<3u> dBm1(dB - _d1);
156 auto AB = factory.make_child(a,b);
158 if (where == InBra) {
159 expr_ += Vector(dB)[dxyz] * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
160 if (where == InKet) {
161 expr_ -= Vector(dB)[dxyz] * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
176 template <CartesianAxis Axis, FunctionPosition where>
179 GenIntegralSet_1_1<CGF1d<Axis>,OverlapOper,EmptySet> >
187 static const unsigned int max_nchildren = 9;
188 static constexpr
auto axis = Axis;
196 using ParentType::RecurrenceRelation::expr_;
197 using ParentType::RecurrenceRelation::nflops_;
198 using ParentType::target_;
204 static std::string descr() {
return std::string(
"OSVRROverlap") +
to_string(axis); }
207 template <CartesianAxis Axis, FunctionPosition where>
208 VRR_1_Overlap_1_1d<Axis,where>::VRR_1_Overlap_1_1d(
const SafePtr< TargetType >& Tint,
214 using namespace libint2::algebra;
215 using namespace libint2::prefactor;
217 typedef CGF1d<Axis> F;
218 const F& _1 = unit<F>(dir);
223 if (a.contracted() ||
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 ||
234 typedef TargetType ChildType;
235 ChildFactory<ThisType,ChildType> factory(
this);
240 auto zero = Tint->bra(0,0).zero() and Tint->ket(0,0).zero() and not deriv;
242 SafePtr<DGVertex> int00 = Vector(
"_0_Overlap_0")[Axis];
243 expr_ = Scalar(0u) + int00;
251 auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
253 auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
256 auto AB = factory.make_child(a,b);
257 if (
is_simple()) { expr_ = Vector(where == InBra ?
"PA" :
"PB")[Axis] * AB; nflops_+=1; }
259 auto am1 = a - _1;
auto am1_exists =
exists(am1);
260 auto bm1 = b - _1;
auto bm1_exists =
exists(bm1);
263 auto Am1B = factory.make_child(am1,b);
264 if (
is_simple()) { expr_ += (Scalar(a.qn()) * Scalar(
"oo2z")) * Am1B; nflops_+=3; }
267 auto ABm1 = factory.make_child(a,bm1);
268 if (
is_simple()) { expr_ += (Scalar(b.qn()) * Scalar(
"oo2z")) * ABm1; nflops_+=3; }
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) );
280 OriginDerivative<1u> _d1; _d1.inc(0);
282 SafePtr<DGVertex> _nullptr;
286 const OriginDerivative<1u> dAm1(dA - _d1);
289 auto AB = factory.make_child(a,b);
291 if (where == InBra) {
292 expr_ -= Scalar(dA[0]) * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
293 if (where == InKet) {
294 expr_ += Scalar(dA[0]) * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
302 const OriginDerivative<1u> dBm1(dB - _d1);
305 auto AB = factory.make_child(a,b);
307 if (where == InBra) {
308 expr_ += Scalar(dB[0]) * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
309 if (where == InKet) {
310 expr_ -= Scalar(dB[0]) * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
326 template <
class BFSet, FunctionPosition where>
329 GenIntegralSet_1_1<BFSet,KineticOper,EmptySet> >
337 static const unsigned int max_nchildren = 9;
345 using ParentType::RecurrenceRelation::expr_;
346 using ParentType::RecurrenceRelation::nflops_;
347 using ParentType::target_;
353 static std::string descr() {
return "OSVRRKinetic"; }
356 template <
class F, FunctionPosition where>
357 VRR_1_Kinetic_1<F,where>::VRR_1_Kinetic_1(
const SafePtr< TargetType >& Tint,
361 using namespace libint2::algebra;
362 using namespace libint2::prefactor;
364 const F& _1 = unit<F>(dir);
369 if (a.contracted() ||
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 ||
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);
388 auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
390 auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
393 auto AB = factory.make_child(a,b);
394 if (
is_simple()) { expr_ = Vector(where == InBra ?
"PA" :
"PB")[dir] * AB; nflops_+=1; }
398 auto Am1B = factory.make_child(am1,b);
399 auto S_Am1B = (where == InBra) ? overlap_factory.make_child(am1,b) : SafePtr<DGVertex>();
401 if (where == InBra) {
402 expr_ += Scalar(a[dir]) * ( Scalar(
"oo2z") * Am1B -
403 Scalar(
"rho12_over_alpha1") * S_Am1B );
407 expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * Am1B;
414 auto ABm1 = factory.make_child(a,bm1);
415 auto S_ABm1 = (where == InKet) ? overlap_factory.make_child(a,bm1) : SafePtr<DGVertex>();
417 if (where == InKet) {
418 expr_ += Scalar(b[dir]) * ( Scalar(
"oo2z") * ABm1 -
419 Scalar(
"rho12_over_alpha2") * S_ABm1 );
423 expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * ABm1;
430 auto S_AB_target = where == InBra ? overlap_factory.make_child(a + _1,b) : overlap_factory.make_child(a,b + _1);
432 expr_ += Scalar(
"two_rho12") * S_AB_target;
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) );
451 for(
unsigned int dxyz=0; dxyz<3; ++dxyz) {
456 OriginDerivative<3u> _d1; _d1.inc(dxyz);
458 SafePtr<DGVertex> _nullptr;
462 const OriginDerivative<3u> dAm1(dA - _d1);
465 auto AB = factory.make_child(a,b);
467 if (where == InBra) {
468 expr_ -= Vector(dA)[dxyz] * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
469 if (where == InKet) {
470 expr_ += Vector(dA)[dxyz] * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
478 const OriginDerivative<3u> dBm1(dB - _d1);
481 auto AB = factory.make_child(a,b);
483 if (where == InBra) {
484 expr_ += Vector(dB)[dxyz] * Scalar(
"rho12_over_alpha1") * AB; nflops_ += 3; }
485 if (where == InKet) {
486 expr_ -= Vector(dB)[dxyz] * Scalar(
"rho12_over_alpha2") * AB; nflops_ += 3; }
501 template <
class BFSet, FunctionPosition where>
504 GenIntegralSet_1_1<BFSet,ElecPotOper,mType> >
512 static const unsigned int max_nchildren = 9;
520 using ParentType::RecurrenceRelation::expr_;
521 using ParentType::RecurrenceRelation::nflops_;
522 using ParentType::target_;
528 static std::string descr() {
return "OSVRRElecPot"; }
534 typedef typename TargetType::AuxIndexType
mType;
535 static SafePtr<mType> aux0(
new mType(0u));
538 << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
543 template <
class F, FunctionPosition where>
544 VRR_1_ElecPot_1<F,where>::VRR_1_ElecPot_1(
const SafePtr< TargetType >& Tint,
548 using namespace libint2::algebra;
549 using namespace libint2::prefactor;
551 const unsigned int m = Tint->aux()->elem(0);
552 const F& _1 = unit<F>(dir);
557 if (a.contracted() ||
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 ||
568 typedef TargetType ChildType;
569 ChildFactory<ThisType,ChildType> factory(
this);
574 auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
576 auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
579 auto AB_m = factory.make_child(a,b,m);
580 auto AB_mp1 = factory.make_child(a,b,m+1);
582 expr_ = Vector(where == InBra ?
"PA" :
"PB")[dir] * AB_m -
583 Vector(
"PC")[dir] * AB_mp1;
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; }
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; }
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) );
611 for(
unsigned int dxyz=0; dxyz<3; ++dxyz) {
616 OriginDerivative<3u> _d1; _d1.inc(dxyz);
618 SafePtr<DGVertex> _nullptr;
622 const OriginDerivative<3u> dAm1(dA - _d1);
625 auto AB_m = factory.make_child(a,b,m);
626 auto AB_mp1 = factory.make_child(a,b,m+1);
628 if (where == InBra) {
629 expr_ -= Vector(dA)[dxyz] * (Scalar(
"rho12_over_alpha1") * AB_m + Scalar(
"rho12_over_alpha2") * AB_mp1); nflops_ += 5; }
630 if (where == InKet) {
631 expr_ += Vector(dA)[dxyz] * Scalar(
"rho12_over_alpha2") * (AB_m - AB_mp1); nflops_ += 4; }
639 const OriginDerivative<3u> dBm1(dB - _d1);
642 auto AB_m = factory.make_child(a,b,m);
643 auto AB_mp1 = factory.make_child(a,b,m+1);
645 if (where == InBra) {
646 expr_ += Vector(dB)[dxyz] * Scalar(
"rho12_over_alpha1") * (AB_m - AB_mp1); nflops_ += 4; }
647 if (where == InKet) {
648 expr_ -= Vector(dB)[dxyz] * (Scalar(
"rho12_over_alpha2") * AB_m + Scalar(
"rho12_over_alpha1") * AB_mp1); nflops_ += 5; }
663 template <
class BFSet, FunctionPosition where>
666 GenIntegralSet_1_1<BFSet,SphericalMultipoleOper,EmptySet> >
675 static const unsigned int max_nchildren = 8;
683 using ParentType::RecurrenceRelation::expr_;
684 using ParentType::RecurrenceRelation::nflops_;
685 using ParentType::target_;
687 using typename ParentType::RecurrenceRelation::ExprType;
689 static std::string descr() {
return "OSVRRSMultipole"; }
692 VRR_1_SMultipole_1(
const SafePtr<TargetType>& Tint,
unsigned int dir) :
695 using namespace libint2::algebra;
696 using namespace libint2::prefactor;
698 using Sign = SphericalMultipoleQuanta::Sign;
700 const F& _1 = unit<F>(dir);
705 if (a.contracted() ||
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 ||
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();
725 typedef TargetType ChildType;
726 ChildFactory<ThisType,ChildType> factory(
this);
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);
733 if (Tint->bra(0,0).norm() != 0 || Tint->ket(0,0).norm() != 0) {
737 auto a = ( where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
739 auto b = ( where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
746 if (
is_simple()) { expr_ = Vector(where == InBra ?
"PA" :
"PB")[dir] * AB; nflops_+=1; }
748 auto am1 = a - _1;
auto am1_exists =
exists(am1);
749 auto bm1 = b - _1;
auto bm1_exists =
exists(bm1);
751 SafePtr<ExprType> subexpr;
755 if (
is_simple()) { subexpr += Scalar(a[dir]) * Am1B; nflops_+=1; }
759 if (
is_simple()) { subexpr += Scalar(b[dir]) * ABm1; nflops_+=1; }
764 auto O_lm1_mp1 = SphericalMultipole_Descr(l-1, m+1, sign);
766 auto AB_O_lm1_mp1 =
make_child(a, b, O_lm1_mp1);
768 subexpr += Scalar(O_lm1_mp1.phase() > 0 ? 0.5 : -0.5) * AB_O_lm1_mp1;
773 auto O_lm1_mm1 = SphericalMultipole_Descr(l-1, m-1, sign);
775 auto AB_O_lm1_mm1 =
make_child(a, b, O_lm1_mm1);
777 subexpr -= Scalar(O_lm1_mm1.phase() > 0 ? 0.5 : -0.5) * AB_O_lm1_mm1;
784 const auto m_pfac = sign == Sign::plus ? 1 : -1;
785 auto Om_lm1_mp1 = SphericalMultipole_Descr(l-1, m+1,
flip(sign));
787 auto AB_Om_lm1_mp1 =
make_child(a, b, Om_lm1_mp1);
789 subexpr += Scalar(m_pfac * (Om_lm1_mp1.phase() > 0 ? 0.5 : -0.5)) * AB_Om_lm1_mp1;
794 auto Om_lm1_mm1 = SphericalMultipole_Descr(l-1, m-1,
flip(sign));
796 auto AB_Om_lm1_mm1 =
make_child(a, b, Om_lm1_mm1);
798 subexpr += Scalar(m_pfac * (Om_lm1_mm1.phase() > 0 ? 0.5 : -0.5)) * AB_Om_lm1_mm1;
804 auto O_lm1_m = SphericalMultipole_Descr(l-1, m, sign);
808 subexpr += AB_O_lm1_m;
815 expr_ += Scalar(
"oo2z") * subexpr;
820 if (dir != 0)
return;
821 auto a = Tint->bra(0,0);
822 auto b = Tint->ket(0,0);
826 assert(sign == Sign::plus);
829 typedef GenIntegralSet_1_1<BFSet, CartesianMultipoleOper<3u>,
EmptySet>
831 ChildFactory<ThisType, OverlapType> overlap_factory(
this);
832 auto S00 = overlap_factory.make_child(a, b);
834 expr_ = Scalar(1) * S00;
837 SafePtr<ExprType> subexpr;
843 SphericalMultipole_Descr(l - 1, m - 1, Sign::plus);
844 auto AB_Op_lm1_mm1 =
make_child(a, b, Op_lm1_mm1);
846 subexpr = Scalar(sign == Sign::plus ?
"PO_x" :
"PO_y") *
850 SphericalMultipole_Descr(l - 1, m - 1, Sign::minus);
852 auto AB_Om_lm1_mm1 =
make_child(a, b, Om_lm1_mm1);
854 if (sign == Sign::plus)
855 subexpr -= Scalar(
"PO_y") * AB_Om_lm1_mm1;
857 subexpr += Scalar(
"PO_x") * AB_Om_lm1_mm1;
862 expr_ = Scalar(-0.5 / m) * subexpr;
866 SafePtr<ExprType> subexpr;
868 auto O_lm1_m = SphericalMultipole_Descr(l - 1, m, sign);
872 subexpr = Scalar((2 * l - 1)) * Scalar(
"PO_z") * AB_O_lm1_m;
875 auto O_lm2_m = SphericalMultipole_Descr(l - 2, m, sign);
879 subexpr -= Scalar(
"PO2") * AB_O_lm2_m;
884 expr_ = Scalar(1.0 / ((l + m) * (l - m))) * subexpr;
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