1 /* 2 * Copyright (C) 2004-2021 Edward F. Valeev 3 * 4 * This file is part of Libint. 5 * 6 * Libint is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * Libint is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with Libint. If not, see <http://www.gnu.org/licenses/>. 18 * 19 */ 20 21 #ifndef _libint2_src_bin_libint_cr11r2dotr2g1211_h_ 22 #define _libint2_src_bin_libint_cr11r2dotr2g1211_h_ 23 24 #include <iostream> 25 #include <sstream> 26 #include <string> 27 #include <vector> 28 #include <stdexcept> 29 #include <cassert> 30 #include <dgvertex.h> 31 #include <rr.h> 32 #include <rr.templ.h> 33 #include <integral.h> 34 #include <r2dotr2g12_11_11.h> 35 #include <algebra.h> 36 #include <flop.h> 37 #include <prefactors.h> 38 #include <context.h> 39 #include <default_params.h> 40 41 namespace libint2 { 42 43 /** Compute relation for 2-e integrals of the r2.r2 x G12 operators. 44 I<BFSet> is the integral set specialization that describes the 45 integrals of the R2dotR2_G12 operator. 46 */ 47 template <template <class> class I, class BFSet> 48 class CR_11_R2dotR2G12_11 : public RecurrenceRelation 49 { 50 51 public: 52 typedef RecurrenceRelation ParentType; 53 typedef BFSet BasisFunctionType; 54 typedef CR_11_R2dotR2G12_11<I,BFSet> ThisType; 55 typedef I<BFSet> TargetType; 56 typedef R12kG12_11_11<BFSet,0> ChildType; 57 /// The type of expressions in which RecurrenceRelations result. 58 typedef RecurrenceRelation::ExprType ExprType; 59 60 /** Use Instance() to obtain an instance of RR. This function is provided to avoid 61 issues with getting a SafePtr from constructor (as needed for registry to work). 62 */ 63 static SafePtr<ThisType> Instance(const SafePtr<TargetType>&); ~CR_11_R2dotR2G12_11()64 virtual ~CR_11_R2dotR2G12_11() {} 65 66 /// Implementation of RecurrenceRelation::num_children() num_children()67 unsigned int num_children() const override { return nchildren_; } 68 /// target() returns pointer to the i-th child target()69 SafePtr<TargetType> target() const { return target_; }; 70 /// child(i) returns pointer to the i-th child 71 SafePtr<ChildType> child(unsigned int i) const; 72 /// Implementation of RecurrenceRelation::rr_target() rr_target()73 SafePtr<DGVertex> rr_target() const override { return static_pointer_cast<DGVertex,TargetType>(target()); } 74 /// Implementation of RecurrenceRelation::rr_child() rr_child(unsigned int i)75 SafePtr<DGVertex> rr_child(unsigned int i) const override { return dynamic_pointer_cast<DGVertex,ChildType>(child(i)); } 76 /// Implementation of RecurrenceRelation::is_simple() is_simple()77 bool is_simple() const override { 78 return TrivialBFSet<BFSet>::result; 79 } 80 81 private: 82 /** 83 dir specifies which quantum number is incremented. 84 For example, dir can be 0 (x), 1(y), or 2(z) if BFSet is 85 a Cartesian Gaussian. 86 */ 87 CR_11_R2dotR2G12_11(const SafePtr<TargetType>&); 88 89 #if 0 90 /// registers this RR with the stack, if needed 91 bool register_with_rrstack() const; 92 #endif 93 94 static const unsigned int max_nchildren_ = 3; 95 SafePtr<TargetType> target_; 96 SafePtr<ChildType> children_[max_nchildren_]; 97 unsigned int nchildren_; 98 generate_label()99 std::string generate_label() const override 100 { 101 std::ostringstream os; 102 os << "RR ( " << rr_target()->label() << " )"; 103 return os.str(); 104 } 105 106 107 }; 108 109 template <template <class> class I, class F> 110 SafePtr< CR_11_R2dotR2G12_11<I,F> > Instance(const SafePtr<TargetType> & Tint)111 CR_11_R2dotR2G12_11<I,F>::Instance(const SafePtr<TargetType>& Tint) 112 { 113 SafePtr<ThisType> this_ptr(new ThisType(Tint)); 114 // Do post-construction duties 115 if (this_ptr->num_children() != 0) { 116 this_ptr->register_with_rrstack<ThisType>(); 117 return this_ptr; 118 } 119 return SafePtr<ThisType>(); 120 } 121 122 template <template <class> class I, class F> CR_11_R2dotR2G12_11(const SafePtr<I<F>> & Tint)123 CR_11_R2dotR2G12_11<I,F>::CR_11_R2dotR2G12_11(const SafePtr<I<F> >& Tint) : 124 ParentType(), target_(Tint), nchildren_(0) 125 { 126 F sh_a(Tint->bra(0,0)); 127 F sh_b(Tint->ket(0,0)); 128 F sh_c(Tint->bra(1,0)); 129 F sh_d(Tint->ket(1,0)); 130 131 std::vector<F> bra; 132 std::vector<F> ket; 133 bra.push_back(sh_a); 134 bra.push_back(sh_c); 135 ket.push_back(sh_b); 136 ket.push_back(sh_d); 137 138 auto* bra_ref = &bra; 139 auto* ket_ref = &ket; 140 141 const unsigned int ndirs = is_simple() ? 3 : 1; 142 for(int xyz=0; xyz<ndirs; xyz++) { 143 144 // c+1_i d+1_i 145 bra_ref->operator[](1).inc(xyz); 146 ket_ref->operator[](1).inc(xyz); 147 int next_child = nchildren_; 148 children_[next_child] = ChildType::Instance(bra[0],ket[0],bra[1],ket[1],0); 149 ++nchildren_; 150 151 if (is_simple()) { 152 SafePtr<ExprType> expr0_ptr(new ExprType(ExprType::OperatorTypes::Times,Scalar(1.0),rr_child(next_child))); 153 add_expr(expr0_ptr); 154 nflops_ += 1; 155 } 156 bra_ref->operator[](1).dec(xyz); 157 ket_ref->operator[](1).dec(xyz); 158 } 159 } 160 161 template <template <class> class I, class F> 162 SafePtr< typename CR_11_R2dotR2G12_11<I,F>::ChildType > child(unsigned int i)163 CR_11_R2dotR2G12_11<I,F>::child(unsigned int i) const 164 { 165 assert(i>=0 && i<nchildren_); 166 unsigned int nc=0; 167 for(int c=0; c<max_nchildren_; c++) { 168 if (children_[c] != 0) { 169 if (nc == i) 170 return children_[c]; 171 nc++; 172 } 173 } 174 }; 175 176 /// Useful typedefs 177 typedef CR_11_R2dotR2G12_11<R2dotR2G12_11_11,CGShell> CR_11_R2dotR2G12_11_sq; 178 typedef CR_11_R2dotR2G12_11<R2dotR2G12_11_11,CGF> CR_11_R2dotR2G12_11_int; 179 180 }; 181 182 #endif 183