vrr_11_twoprep_11.h

/*
 *  Copyright (C) 2004-2021 Edward F. Valeev
 *
 *  This file is part of Libint.
 *
 *  Libint is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_bin_libint_vrr11twoprep11_h_
#define _libint2_src_bin_libint_vrr11twoprep11_h_
 
#include <generic_rr.h>
#include <twoprep_11_11.h>
 
namespace libint2 {
 
    template <class BFSet, int part, FunctionPosition where>
      class VRR_11_TwoPRep_11 : public GenericRecurrenceRelation< VRR_11_TwoPRep_11<BFSet,part,where>,
                                                                  BFSet,
                                                                  GenIntegralSet_11_11<BFSet,TwoPRep,mType> >
    {
    public:
      typedef VRR_11_TwoPRep_11 ThisType;
      typedef BFSet BasisFunctionType;
      typedef GenIntegralSet_11_11<BFSet,TwoPRep,mType> TargetType;
      typedef GenericRecurrenceRelation<ThisType,BFSet,TargetType> ParentType;
      friend class GenericRecurrenceRelation<ThisType,BFSet,TargetType>;
      static const unsigned int max_nchildren = 26;
 
      using ParentType::Instance;
 
      static bool directional() { return ParentType::default_directional(); }
 
    private:
      using ParentType::RecurrenceRelation::expr_;
      using ParentType::RecurrenceRelation::nflops_;
      using ParentType::target_;
      using ParentType::is_simple;
 
      VRR_11_TwoPRep_11(const SafePtr<TargetType>&, unsigned int dir);
 
      static std::string descr() { return "OSVRR"; }
      std::string generate_label() const override
      {
        typedef typename TargetType::AuxIndexType mType;
        static SafePtr<mType> aux0(new mType(0u));
        std::ostringstream os;
        os << descr() << "P" << part << to_string(where)
           << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
        return os.str();
      }
 
  #if LIBINT_ENABLE_GENERIC_CODE
      bool has_generic(const SafePtr<CompilationParameters>& cparams) const override;
      std::string generic_header() const override;
      std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const override;
  #endif
    };
 
  template <class F, int part, FunctionPosition where>
    VRR_11_TwoPRep_11<F,part,where>::VRR_11_TwoPRep_11(const SafePtr< TargetType >& Tint,
                                                       unsigned int dir) :
    ParentType(Tint,dir)
    {
      using namespace libint2::algebra;
      using namespace libint2::prefactor;
      using namespace libint2::braket;
      const unsigned int m = Tint->aux()->elem(0);
      const F& _1 = unit<F>(dir);
 
      { // can't apply to contracted basis functions
        F a(Tint->bra(0,0));
        F b(Tint->ket(0,0));
        F c(Tint->bra(1,0));
        F d(Tint->ket(1,0));
        if (a.contracted() ||
          b.contracted() ||
          c.contracted() ||
          d.contracted())
          return;
      }
 
      // if derivative integrals, there will be extra terms (Eq. (143) in Obara & Saika JCP 89)
      const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
      const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
      const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
      const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
      const bool deriv = !dA.zero() ||
          !dB.zero() ||
          !dC.zero() ||
          !dD.zero();
 
      // This is a hack to avoid creating recurrence relations for which generic code has not been yet implemented
#if LIBINT_ENABLE_GENERIC_CODE
      {
        F sh_a(target_->bra(0,0));
        F sh_b(target_->ket(0,0));
        F sh_c(target_->bra(1,0));
        F sh_d(target_->ket(1,0));
        // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
        // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
        if (sh_b.zero() && sh_d.zero() &&
           (sh_a.norm() > 1u && sh_c.norm() > 1u)
           ) { // have a generic implemented ...
          if (part != 0) // ... but only implemented build on A in this case
            return;
        }
        if (sh_a.zero() && sh_c.zero() &&
           (sh_b.norm() > 1u && sh_d.norm() > 1u)
           ) {
          if (part != 0) // but only implemented build on B in this case
            return;
        }
      }
#endif
 
      typedef TargetType ChildType;
      ChildFactory<ThisType,ChildType> factory(this);
 
      bool part0_has_unit=false, part1_has_unit=false;
 
      // Build on A
      if (part == 0 && where == InBra) {
        F a(Tint->bra(0,0) - _1);
        if (!exists(a)) return;
        F b(Tint->ket(0,0)); const bool unit_b = (b == F::unit()); part0_has_unit |= unit_b;
        F c(Tint->bra(1,0));
        F d(Tint->ket(1,0));
 
        SafePtr<DGVertex> ABCD_m; if (not unit_b) ABCD_m = factory.make_child(a,b,c,d,m);
        auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
        if (is_simple()) {
          if (not unit_b) {
            expr_ = Vector("PA")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1;  nflops_+=3;
          }
          else {
            expr_ = Vector("WP")[dir] * ABCD_mp1;  nflops_+=1;
          }
        }
 
        // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
        const bool ahlrichs_simplification = a.pure_sh() && unit_b;
        auto am1 = a - _1;
        if (exists(am1) && not ahlrichs_simplification) {
          auto Am1BCD_m = factory.make_child(am1,b,c,d,m);
          auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(a[dir]) * Scalar("oo2z") * (Scalar("roz") * Am1BCD_mp1 - Am1BCD_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1);  nflops_+=5; }
#endif
        }
        const F& bm1 = b - _1;
        if (exists(bm1)) {
          auto ABm1CD_m = factory.make_child(a,bm1,c,d,m);
          auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(b[dir]) * Scalar("oo2z") * (Scalar("roz") * ABm1CD_mp1 - ABm1CD_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1);  nflops_+=5; }
#endif
        }
        const F& cm1 = c - _1;
        if (exists(cm1)) {
          auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
          if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1;  nflops_+=3; }
        }
        const F& dm1 = d - _1;
        if (exists(dm1)) {
          auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
          if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1;  nflops_+=3; }
        }
      }
      // Build on B
      if (part == 0 && where == InKet) {
        F a(Tint->bra(0,0));  const bool unit_a = (a == F::unit()); part0_has_unit |= unit_a;
        F b(Tint->ket(0,0) - _1);
        if (!exists(b)) return;
        F c(Tint->bra(1,0));
        F d(Tint->ket(1,0));
 
        SafePtr<DGVertex> ABCD_m; if (not unit_a) ABCD_m = factory.make_child(a,b,c,d,m);
        auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
        if (is_simple()) {
          if (not unit_a) {
            expr_ = Vector("PB")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1;  nflops_+=3;
          }
          else {
            expr_ = Vector("WP")[dir] * ABCD_mp1;  nflops_+=1;
          }
        }
 
        const F& am1 = a - _1;
        if (exists(am1)) {
          auto Am1BCD_m = factory.make_child(am1,b,c,d,m);
          auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(a[dir]) * Scalar("oo2z") * (Scalar("roz") * Am1BCD_mp1 - Am1BCD_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1);  nflops_+=5; }
#endif
        }
        // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
        const bool ahlrichs_simplification = b.pure_sh() && unit_a;
        const F& bm1 = b - _1;
        if (exists(bm1) && not ahlrichs_simplification) {
          auto ABm1CD_m = factory.make_child(a,bm1,c,d,m);
          auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(b[dir]) * Scalar("oo2z") * (Scalar("roz") * ABm1CD_mp1 - ABm1CD_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1);  nflops_+=5; }
#endif
        }
        const F& cm1 = c - _1;
        if (exists(cm1)) {
          auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
          if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1;  nflops_+=3; }
        }
        const F& dm1 = d - _1;
        if (exists(dm1)) {
          auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
          if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1;  nflops_+=3; }
        }
      }
      // Build on C
      if (part == 1 && where == InBra) {
        F a(Tint->bra(0,0));
        F b(Tint->ket(0,0));
        F c(Tint->bra(1,0) - _1);
        if (!exists(c)) return;
        F d(Tint->ket(1,0)); const bool unit_d = (d == F::unit()); part1_has_unit |= unit_d;
 
        SafePtr<DGVertex> ABCD_m; if (not unit_d) ABCD_m = factory.make_child(a,b,c,d,m);
        auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
        if (is_simple()) {
          if (not unit_d) {
            expr_ = Vector("QC")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1;  nflops_+=3;
          }
          else {
            expr_ = Vector("WQ")[dir] * ABCD_mp1;  nflops_+=1;
          }
        }
 
        // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
        const bool ahlrichs_simplification = c.pure_sh() && unit_d;
        const F& cm1 = c - _1;
        if (exists(cm1) && not ahlrichs_simplification) {
          auto ABCm1D_m = factory.make_child(a,b,cm1,d,m);
          auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(c[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCm1D_mp1 - ABCm1D_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1);  nflops_+=5; }
#endif
        }
        const F& dm1 = d - _1;
        if (exists(dm1)) {
          auto ABCDm1_m = factory.make_child(a,b,c,dm1,m);
          auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(d[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCDm1_mp1 - ABCDm1_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1);  nflops_+=5; }
#endif
        }
        const F& am1 = a - _1;
        if (exists(am1)) {
          auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
          if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1;  nflops_+=3; }
        }
        const F& bm1 = b - _1;
        if (exists(bm1)) {
          auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
          if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1;  nflops_+=3; }
        }
      }
      // Build on D
      if (part == 1 && where == InKet) {
        F a(Tint->bra(0,0));
        F b(Tint->ket(0,0));
        F c(Tint->bra(1,0));  const bool unit_c = (c == F::unit()); part1_has_unit |= unit_c;
        F d(Tint->ket(1,0) - _1);
        if (!exists(d)) return;
 
        SafePtr<DGVertex> ABCD_m; if (not unit_c) ABCD_m = factory.make_child(a,b,c,d,m);
        auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
        if (is_simple()) {
          if (not unit_c) {
            expr_ = Vector("QD")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1;  nflops_+=3;
          }
          else {
            expr_ = Vector("WQ")[dir] * ABCD_mp1;  nflops_+=1;
          }
        }
 
        const F& cm1 = c - _1;
        if (exists(cm1)) {
          auto ABCm1D_m = factory.make_child(a,b,cm1,d,m);
          auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(c[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCm1D_mp1 - ABCm1D_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1);  nflops_+=5; }
#endif
        }
        // simplified 3-center VRR due to Ahlrichs (PCCP 6, 5119 (2004))
        const bool ahlrichs_simplification = d.pure_sh() && unit_c;
        const F& dm1 = d - _1;
        if (exists(dm1) && not ahlrichs_simplification) {
          auto ABCDm1_m = factory.make_child(a,b,c,dm1,m);
          auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1);
#if LIBINT_GENERATE_FMA
          // this form is amenable to generation of fmsub
          if (is_simple()) { expr_ -= Scalar(d[dir]) * Scalar("oo2e") * (Scalar("roe") * ABCDm1_mp1 - ABCDm1_m);  nflops_+=5; }
#else
          if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1);  nflops_+=5; }
#endif
        }
        const F& am1 = a - _1;
        if (exists(am1)) {
          auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1);
          if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1;  nflops_+=3; }
        }
        const F& bm1 = b - _1;
        if (exists(bm1)) {
          auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1);
          if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1;  nflops_+=3; }
        }
      }
 
      // if got here, can decrement by at least 1 quantum
      // add additional derivative terms
      if (deriv) {
        // bf quantum on the build center subtracted by 1
        F a( part == 0 && where == InBra ? Tint->bra(0,0) - _1 : Tint->bra(0,0) );
        F b( part == 0 && where == InKet ? Tint->ket(0,0) - _1 : Tint->ket(0,0) );
        F c( part == 1 && where == InBra ? Tint->bra(1,0) - _1 : Tint->bra(1,0) );
        F d( part == 1 && where == InKet ? Tint->ket(1,0) - _1 : Tint->ket(1,0) );
 
        // treatment of derivative terms differs for shell sets and integrals
        // since in computing shell sets transfer/build will occur in all 3 directions
        // change in up to all three derivative indices will occur
        for(unsigned int dxyz=0; dxyz<3; ++dxyz) {
 
          if (is_simple() && dxyz != dir) // for integrals only consider derivatives in THE build direction
            continue;
 
          OriginDerivative<3u> _d1; _d1.inc(dxyz);
 
          SafePtr<DGVertex> _nullptr;
 
          // dA - _1?
          {
            const OriginDerivative<3u> dAm1(dA - _d1);
            if (exists(dAm1)) { // yes
              a.deriv() = dAm1;
              auto ABCD_m = (part == 0 && not part0_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
              auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
              if (is_simple()) {
                if (part == 0 && where == InBra) { // building on A
                  if (not part0_has_unit) {
                    expr_ -= Vector(dA)[dxyz] * (Scalar("rho12_over_alpha1") * ABCD_m + Scalar("alpha1_rho_over_zeta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dA)[dxyz] * (Scalar("alpha1_rho_over_zeta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 0 && where == InKet) { // building on B
                  if (not part0_has_unit) {
                    expr_ += Vector(dA)[dxyz] * (Scalar("rho12_over_alpha2") * ABCD_m - Scalar("alpha1_rho_over_zeta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dA)[dxyz] * (Scalar("alpha1_rho_over_zeta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 1) { // building on C or D
                  expr_ += Vector(dA)[dxyz] * Scalar("alpha1_over_zetapluseta") * ABCD_mp1;  nflops_ += 3; }
              }
              a.deriv() = dA;
            }
          }
 
          // dB - _1?
          {
            const OriginDerivative<3u> dBm1(dB - _d1);
            if (exists(dBm1)) { // yes
              b.deriv() = dBm1;
              auto ABCD_m = (part == 0 && not part0_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
              auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
              if (is_simple()) {
                if (part == 0 && where == InBra) { // building on A
                  if (not part0_has_unit) {
                    expr_ += Vector(dB)[dxyz] * (Scalar("rho12_over_alpha1") * ABCD_m - Scalar("alpha2_rho_over_zeta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dB)[dxyz] * (Scalar("alpha2_rho_over_zeta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 0 && where == InKet) { // building on B
                  if (not part0_has_unit) {
                    expr_ -= Vector(dB)[dxyz] * (Scalar("rho12_over_alpha2") * ABCD_m + Scalar("alpha2_rho_over_zeta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dB)[dxyz] * (Scalar("alpha2_rho_over_zeta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 1) { // building on C or D
                  expr_ += Vector(dB)[dxyz] * Scalar("alpha2_over_zetapluseta") * ABCD_mp1;  nflops_ += 3; }
              }
              b.deriv() = dB;
            }
          }
 
          // dC - _1?
          {
            const OriginDerivative<3u> dCm1(dC - _d1);
            if (exists(dCm1)) { // yes
              c.deriv() = dCm1;
              auto ABCD_m = (part == 1 && not part1_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
              auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
              if (is_simple()) {
                if (part == 0) { // building on A or B
                  expr_ += Vector(dC)[dxyz] * Scalar("alpha3_over_zetapluseta") * ABCD_mp1;  nflops_ += 3; }
                if (part == 1 && where == InBra) { // building on C
                  if (not part1_has_unit) {
                    expr_ -= Vector(dC)[dxyz] * (Scalar("rho34_over_alpha3") * ABCD_m + Scalar("alpha3_rho_over_eta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dC)[dxyz] * (Scalar("alpha3_rho_over_eta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 1 && where == InKet) { // building on D
                  if (not part1_has_unit) {
                    expr_ += Vector(dC)[dxyz] * (Scalar("rho34_over_alpha4") * ABCD_m - Scalar("alpha3_rho_over_eta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dC)[dxyz] * (Scalar("alpha3_rho_over_eta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
              }
              c.deriv() = dC;
            }
          }
 
          // dD - _1?
          {
            const OriginDerivative<3u> dDm1(dD - _d1);
            if (exists(dDm1)) { // yes
              d.deriv() = dDm1;
              auto ABCD_m = (part == 1 && not part1_has_unit) ? factory.make_child(a,b,c,d,m) : _nullptr;
              auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1);
              if (is_simple()) {
                if (part == 0) { // building on A or B
                  expr_ += Vector(dD)[dxyz] * Scalar("alpha4_over_zetapluseta") * ABCD_mp1;  nflops_ += 3; }
                if (part == 1 && where == InBra) { // building on C
                  if (not part1_has_unit) {
                    expr_ += Vector(dD)[dxyz] * (Scalar("rho34_over_alpha3") * ABCD_m - Scalar("alpha4_rho_over_eta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dD)[dxyz] * (Scalar("alpha4_rho_over_eta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
                if (part == 1 && where == InKet) { // building on D
                  if (not part1_has_unit) {
                    expr_ -= Vector(dD)[dxyz] * (Scalar("rho34_over_alpha4") * ABCD_m + Scalar("alpha4_rho_over_eta2") * ABCD_mp1);  nflops_ += 5;
                  }
                  else {
                    expr_ -= Vector(dD)[dxyz] * (Scalar("alpha4_rho_over_eta2") * ABCD_mp1);  nflops_ += 3;
                  }
                }
              }
              d.deriv() = dD;
            }
          }
        }
      } // end of deriv
 
      return;
    }
 
#if LIBINT_ENABLE_GENERIC_CODE
  template <class F, int part, FunctionPosition where>
    bool
    VRR_11_TwoPRep_11<F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
    {
      if (TrivialBFSet<F>::result)
        return false;
 
      F sh_a(target_->bra(0,0));
      F sh_b(target_->ket(0,0));
      F sh_c(target_->bra(1,0));
      F sh_d(target_->ket(1,0));
      const unsigned int max_opt_am = cparams->max_am_opt();
      // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
      // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
      if (sh_b.zero() && sh_d.zero() &&
          (sh_a.norm() > std::max(2*max_opt_am,1u) ||
           sh_c.norm() > std::max(2*max_opt_am,1u)
          ) &&
          (sh_a.norm() > 1u && sh_c.norm() > 1u)
         ) {
        assert(part == 0); // has only implemented build on A in this case
        return true;
      }
      if (sh_a.zero() && sh_c.zero() &&
          (sh_b.norm() > std::max(2*max_opt_am,1u) ||
           sh_d.norm() > std::max(2*max_opt_am,1u)
          ) &&
          (sh_b.norm() > 1u && sh_d.norm() > 1u)
         ) {
        assert(part == 0); // has only implemented build on B in this case
        return true;
      }
      return false;
    }
 
  template <class F, int part, FunctionPosition where>
    std::string
    VRR_11_TwoPRep_11<F,part,where>::generic_header() const
    {
      F sh_a(target_->bra(0,0));
      F sh_b(target_->ket(0,0));
      F sh_c(target_->bra(1,0));
      F sh_d(target_->ket(1,0));
      const bool xsxs = sh_b.zero() && sh_d.zero();
      const bool sxsx = sh_a.zero() && sh_c.zero();
 
      const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
      const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
      const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
      const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
      const bool deriv = !dA.zero() ||
          !dB.zero() ||
          !dC.zero() ||
          !dD.zero();
 
      if (!deriv) {
        if (xsxs) return std::string("OSVRR_xs_xs.h");
        if (sxsx) return std::string("OSVRR_sx_sx.h");
      }
      else {
        if (xsxs) return std::string("OSVRR_xs_xs_deriv.h");
        if (sxsx) return std::string("OSVRR_sx_sx_deriv.h");
      }
      abort(); // unreachable
    }
 
  template <class F, int part, FunctionPosition where>
    std::string
    VRR_11_TwoPRep_11<F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
    {
      std::ostringstream oss;
      F sh_a(target_->bra(0,0));
      F sh_b(target_->ket(0,0));
      F sh_c(target_->bra(1,0));
      F sh_d(target_->ket(1,0));
      const bool xsxs = sh_b.zero() && sh_d.zero();
      const bool sxsx = sh_a.zero() && sh_c.zero();
      bool ahlrichs_simplification = false;
      bool unit_s = false;
      bool part0_has_unit = false;
      bool part1_has_unit = false;
      if (xsxs) {
        ahlrichs_simplification = (sh_a.pure_sh() && sh_b.is_unit()) ||
                                  (sh_c.pure_sh() && sh_d.is_unit());
        unit_s = sh_b.is_unit();
        part0_has_unit = sh_b.is_unit();
        part1_has_unit = sh_d.is_unit();
      }
      if (sxsx) {
        ahlrichs_simplification = (sh_b.pure_sh() && sh_a.is_unit()) ||
                                  (sh_d.pure_sh() && sh_c.is_unit());
        unit_s = sh_a.is_unit();
        part0_has_unit = sh_a.is_unit();
        part1_has_unit = sh_c.is_unit();
      }
 
      const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
      const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
      const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
      const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
      const bool deriv = !dA.zero() ||
          !dB.zero() ||
          !dC.zero() ||
          !dD.zero();
 
      oss << "using namespace libint2;" << std::endl;
 
      if (!deriv) { // for regular integrals I know exactly how many prerequisites I need
        if(xsxs) {
          oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
              << "VRR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
        }
        if (sxsx) {
          oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
              << "VRR_sx_sx<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
        }
        if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
        oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
        oss << ">::compute(inteval";
 
        oss << "," << args->symbol(0); // target
        if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
          oss << ",0"; // src0-> 0x0
        const unsigned int nargs = args->n();
        for(unsigned int a=1; a<nargs; a++) {
          oss << "," << args->symbol(a);
        }
        oss << ");";
      }
      else { // deriv == true -> only some arguments are needed
        if(xsxs) {
          oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
              << "VRR_xs_xs_deriv<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
        }
        if(sxsx) {
          oss << "libint2::OS" << (ahlrichs_simplification ? "A" : "")
              << "VRR_sx_sx_deriv<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
        }
 
        for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_a.deriv().d(xyz) << ",";
        for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_b.deriv().d(xyz) << ",";
        for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_c.deriv().d(xyz) << ",";
        for(unsigned int xyz=0; xyz<3; ++xyz) oss << sh_d.deriv().d(xyz) << ",";
        if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
        oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
        oss << ">::compute(inteval";
        // out of all 22 possible prerequisites first few have same derivative degree as the target
        // 5 if standard 4-center integral
        // 1+4 if the center opposite the build center carries a unit function
        //     will pass 0 instead of src0
        // 2 if Ahlrichs
        oss << "," << args->symbol(0); // target
        if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
          oss << ",0"; // src0-> 0x0
        //const unsigned int nargs = args->n();
        unsigned int arg = 1;
        for(;
            arg<(ahlrichs_simplification? 3 : (unit_s ? 5 : 6));  // nargs + 1 target
            arg++) {
          oss << "," << args->symbol(arg);
        }
        for(unsigned int xyz=0; xyz<3; ++xyz) {
          if (sh_a.deriv().d(xyz) > 0) {
            // see the dA-1 clause: (ab|cd)^m is skipped
            oss << "," << (part0_has_unit ? std::to_string(0) : args->symbol(arg++));
            oss << "," << args->symbol(arg++);
          }
          else
            oss << ",0,0";
          if (sh_b.deriv().d(xyz) > 0) {
            // see the dB-1 clause: (ab|cd)^m is skipped
            oss << "," << (part0_has_unit ? std::to_string(0) : args->symbol(arg++));
            oss << "," << args->symbol(arg++);
          }
          else
            oss << ",0,0";
          if (sh_c.deriv().d(xyz) > 0) {
            oss << "," << args->symbol(arg++);
          }
          else
            oss << ",0";
          if (sh_d.deriv().d(xyz) > 0) {
            oss << "," << args->symbol(arg++);
          }
          else
            oss << ",0";
        }
        oss << ");";
      }
 
      return oss.str();
    }
#endif // #if !LIBINT_ENABLE_GENERIC_CODE
 
};
 
#endif