bfset.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_bfset_h_
#define _libint2_src_bin_libint_bfset_h_
 
#include <iostream>
#include <string>
#include <cassert>
#include <numeric>
#include <sstream>
#include <stdarray.h>
#include <smart_ptr.h>
#include <polyconstr.h>
#include <hashable.h>
#include <contractable.h>
#include <global_macros.h>
#include <util_types.h>
 
namespace libint2 {
 
  class BFSet : public ConstructablePolymorphically {
 
  public:
    virtual ~BFSet() {}
    virtual unsigned int num_bf() const =0;
    virtual std::string label() const =0;
 
  protected:
    BFSet() {}
 
  };
 
  class IncableBFSet : public BFSet {
 
  public:
    virtual ~IncableBFSet() {}
 
    virtual void inc(unsigned int xyz, unsigned int c = 1u) =0;
    virtual void dec(unsigned int xyz, unsigned int c = 1u) =0;
    virtual unsigned int norm() const =0;
    bool zero() const { return norm() == 0; }
    bool valid() const { return valid_; }
 
  protected:
    IncableBFSet() : valid_(true) {}
 
    void invalidate() { valid_ = false; }
 
  private:
    bool valid_;
  };
 
  template <typename F>
  F unit(unsigned int X) { F tmp; tmp.inc(X,1); return tmp; }
  inline bool exists(const IncableBFSet& A) { return A.valid(); }
 
  template <unsigned NDIM = 3>
  class OriginDerivative : public Hashable<LIBINT2_UINT_LEAST64,ReferToKey> {
 
    static_assert(NDIM == 1 || NDIM == 3, "OriginDerivative<NDIM>: only NDIM=1 or NDIM=3 are supported");
  public:
    OriginDerivative() : valid_(true) {
      for(auto d=0u; d!=NDIM; ++d) d_[d] = 0u;
    }
    OriginDerivative(const OriginDerivative& other) : valid_(true) {
      std::copy(other.d_, other.d_ + NDIM, d_);
    }
    OriginDerivative& operator=(const OriginDerivative& other) {
      valid_ = other.valid_;
      std::copy(other.d_, other.d_ + NDIM, d_);
      return *this;
    }
    OriginDerivative& operator+=(const OriginDerivative& other) {
      assert(valid_);
      for(auto d=0u; d!=NDIM; ++d) d_[d] += other.d_[d];
      return *this;
    }
    OriginDerivative& operator-=(const OriginDerivative& other) {
      assert(valid_);
      for(auto d=0u; d!=NDIM; ++d) d_[d] -= other.d_[d];
      return *this;
    }
    ~OriginDerivative() {
      static_assert(NDIM == 3u || NDIM == 1u, "OriginDerivative with NDIM=1,3 are implemented");
    }
 
    unsigned int d(unsigned int xyz) const {
      assert(xyz < NDIM);
      return d_[xyz];
    }
    unsigned int operator[](unsigned int xyz) const {
      return this->d(xyz);
    }
    void inc(unsigned int xyz, unsigned int c = 1u) {
      assert(xyz < NDIM);
      assert(valid_);
      d_[xyz] += c;
    }
    void dec(unsigned int xyz, unsigned int c = 1u) {
      assert(xyz < NDIM);
      //assert(valid_);
      if (d_[xyz] >= c)
        d_[xyz] -= c;
      else
        valid_ = false;
    }
    unsigned int norm() const {
      return std::accumulate(d_, d_+NDIM, 0u);
    }
    bool zero() const { return norm() == 0; }
    bool valid() const { return valid_; }
    LIBINT2_UINT_LEAST64 key() const override {
      if (NDIM == 3u) {
        unsigned nxy = d_[1] + d_[2];
        unsigned l = nxy + d_[0];
        LIBINT2_UINT_LEAST64 key = nxy*(nxy+1)/2 + d_[2];
        const auto result = key + key_l_offset.at(l);
        assert(result < max_key);
        return result;
      }
      if (NDIM == 1u) {
        const auto result = d_[0];
        assert(result < max_key);
        return result;
      }
      assert(false);
    }
    std::string label() const {
      char result[NDIM+1];
      for(auto xyz=0u; xyz<NDIM; ++xyz) result[xyz] = '0' + d_[xyz];
      result[NDIM] = '\0';
      return std::string(result);
    }
 
    /* ---------------
     * key computation
     * --------------- */
    const static unsigned max_deriv = 4;
 
    // nkeys_l[L] is the number of possible OriginDerivative's with \c norm() == L
    // \note for NDIM=1 nkeys_l[L] = 1
    // \note for NDIM=2 nkeys_l[L] = (L+1)
    // \note for NDIM=3 nkeys_l[L] = (L+1)(L+2)/2
    //static std::array<LIBINT2_UINT_LEAST64, OriginDerivative::max_deriv+1> nkeys;
 
    const static unsigned max_key = NDIM == 3 ? (1 + max_deriv)*(2 + max_deriv)*(3 + max_deriv)/6 : (1+max_deriv);
 
    void print(std::ostream& os = std::cout) const;
 
  protected:
 
    void invalidate() { valid_ = false; }
 
  private:
    unsigned int d_[NDIM];
    bool valid_;  // indicates valid/invalid state
    static std::array<LIBINT2_UINT_LEAST64, OriginDerivative::max_deriv+1> key_l_offset;
 
  };
 
  // explicit instantiation declaration, definitions are gauss.cc
  template<>
  std::array<LIBINT2_UINT_LEAST64, OriginDerivative<1u>::max_deriv+1> OriginDerivative<1u>::key_l_offset;
  template<>
  std::array<LIBINT2_UINT_LEAST64, OriginDerivative<3u>::max_deriv+1> OriginDerivative<3u>::key_l_offset;
 
  template <unsigned NDIM>
  OriginDerivative<NDIM> operator-(const OriginDerivative<NDIM>& A, const OriginDerivative<NDIM>& B) {
    OriginDerivative<NDIM> Diff(A);
    for(unsigned int xyz=0; xyz<3; ++xyz)
      Diff.dec(xyz,B.d(xyz));
    return Diff;
  }
 
  template <unsigned NDIM>
  bool operator==(const OriginDerivative<NDIM>& A, const OriginDerivative<NDIM>& B) {
    for(unsigned d=0; d!=NDIM; ++d)
      if (A.d(d) != B.d(d))
        return false;
    return true;
  }
 
  template <unsigned NDIM> inline bool exists(const OriginDerivative<NDIM>& A) { return A.valid(); }
 
  class CGF;   // forward declaration of CGF
 
  class CGShell : public IncableBFSet, public Hashable<LIBINT2_UINT_LEAST64,ReferToKey>,
                  public Contractable<CGShell> {
 
    unsigned int qn_[1];
    OriginDerivative<3> deriv_;
    bool pure_sh_;  //< if true, assumed to contain solid harmonics with quantum number qn_[0] only
    bool unit_;
 
    friend CGShell operator+(const CGShell& A, const CGShell& B);
    friend CGShell operator-(const CGShell& A, const CGShell& B);
 
  public:
    typedef CGF iter_type;
    typedef IncableBFSet parent_type;
 
    CGShell();
    CGShell(unsigned int qn, bool pure_sh = false);
    CGShell(const CGShell&);
    virtual ~CGShell();
    CGShell& operator=(const CGShell&);
 
    const OriginDerivative<3u>& deriv() const { return deriv_; }
    OriginDerivative<3u>& deriv() { return deriv_; }
 
    std::string label() const override;
    unsigned int num_bf() const override { return (qn_[0]+1)*(qn_[0]+2)/2; }
    unsigned int qn(unsigned int xyz=0) const { return qn_[0]; }
    unsigned int operator[](unsigned int xyz) const {
      return this->qn(xyz);
    }
 
    bool operator==(const CGShell&) const;
 
    bool pure_sh() const { return pure_sh_; }
    void pure_sh(bool p) { pure_sh_ = p; }
 
    void inc(unsigned int xyz, unsigned int c = 1u) override;
    void dec(unsigned int xyz, unsigned int c = 1u) override;
    unsigned int norm() const override;
    LIBINT2_UINT_LEAST64 key() const override {
      if (is_unit()) return max_key-1;
      const LIBINT2_UINT_LEAST64 result =
             ((deriv().key() * 2 +
               (contracted() ? 1 : 0)
              ) * (max_qn+1) +
              qn_[0]
             ) * 2 +
             (pure_sh() ? 1 : 0);
      assert(result < max_key-1);
      return result;
    }
    const static LIBINT2_UINT_LEAST64 max_qn = LIBINT_CARTGAUSS_MAX_AM;
    const static LIBINT2_UINT_LEAST64 max_key = OriginDerivative<3u>::max_key * 2 * (max_qn + 1) * 2 + 1;
 
    void print(std::ostream& os = std::cout) const;
 
    static CGShell unit();
    bool is_unit() const { return unit_; }
 
  };
 
  CGShell operator+(const CGShell& A, const CGShell& B);
  CGShell operator-(const CGShell& A, const CGShell& B);
 
  class CGF : public IncableBFSet, public Hashable<LIBINT2_UINT_LEAST64,ComputeKey>,
              public Contractable<CGF> {
 
    unsigned int qn_[3];
    OriginDerivative<3u> deriv_;
    bool pure_sh_;  //< if true, assumed to contain solid harmonics with quantum number qn_[0] only
    bool unit_; //< if true, this is a unit Gaussian (exponent = 0)
 
    friend CGF operator+(const CGF& A, const CGF& B);
    friend CGF operator-(const CGF& A, const CGF& B);
 
  public:
    typedef CGF iter_type;
    typedef IncableBFSet parent_type;
    //typedef typename Hashable<unsigned,ComputeKey>::KeyReturnType KeyReturnType;
 
    CGF();
    CGF(unsigned int qn[3], bool pure_sh = false);
    CGF(const CGF&);
    explicit CGF(const ConstructablePolymorphically&);
    virtual ~CGF();
    CGF& operator=(const CGF&);
 
    const OriginDerivative<3u>& deriv() const { return deriv_; }
    OriginDerivative<3u>& deriv() { return deriv_; }
 
    std::string label() const override;
    unsigned int num_bf() const override { return 1; }
    unsigned int qn(unsigned int axis) const;
    unsigned int operator[](unsigned int axis) const {
      return qn(axis);
    }
 
    bool pure_sh() const { return pure_sh_; }
    void pure_sh(bool p) { pure_sh_ = p; }
 
    bool operator==(const CGF&) const;
 
    void inc(unsigned int xyz, unsigned int c = 1u) override;
    void dec(unsigned int xyz, unsigned int c = 1u) override;
    unsigned int norm() const override;
    LIBINT2_UINT_LEAST64 key() const override {
      if (is_unit()) return max_key-1;
      unsigned nxy = qn_[1] + qn_[2];
      unsigned l = nxy + qn_[0];
      LIBINT2_UINT_LEAST64 key = nxy*(nxy+1)/2 + qn_[2];
      const LIBINT2_UINT_LEAST64 result =
          ( ( deriv().key() * 2 +
              (contracted() ? 1 : 0)
            ) * max_num_qn +
            key + key_l_offset.at(l)
          ) * 2
          + (pure_sh() ? 1 : 0);
      if (result >= max_key-1) {
        this->print(std::cout);
        std::cout << "result,max_key-1 = " << result << "," << max_key-1 << std::endl;
        assert(result < max_key-1);
      }
      return result;
    }
    const static LIBINT2_UINT_LEAST64 max_num_qn = ((1 + (CGShell::max_qn+1)) * (2 + (CGShell::max_qn+1)) * (3 + (CGShell::max_qn+1)) /6);
    // deriv_key_range = OriginDerivative<3u>::max_key
    // contracted = 2 (yes or no)
    // qn_range = max_num_qn
    // puresh_key_range = 2
    // +1 to account for unit function
    const static LIBINT2_UINT_LEAST64 max_key = OriginDerivative<3u>::max_key * 2ul * max_num_qn * 2ul + 1;
 
    void print(std::ostream& os = std::cout) const;
 
    static CGF unit();
    bool is_unit() const { return unit_; }
 
  private:
    static std::array<LIBINT2_UINT_LEAST64, CGShell::max_qn+1> key_l_offset;
  };
 
  CGF operator+(const CGF& A, const CGF& B);
  CGF operator-(const CGF& A, const CGF& B);
 
#if 1
  template <CartesianAxis Axis>
  class CGF1d : public IncableBFSet,
                public Hashable<LIBINT2_UINT_LEAST64, ComputeKey>,
                public Contractable<CGF1d<Axis> > {
    unsigned int qn_[1];
    OriginDerivative<1u> deriv_;
    bool unit_; //< if true, this is a unit Gaussian (exponent = 0)
 
  public:
 
    static constexpr auto axis = Axis;
 
    typedef CGF1d iter_type;
    typedef IncableBFSet parent_type;
 
    CGF1d() : unit_(false) { qn_[0] = 0; }
    CGF1d(unsigned int qn) : unit_(false) { qn_[0] = qn; }
    CGF1d(unsigned int qn[1]) : unit_(false) { qn_[0] = qn[0]; }
    CGF1d(const CGF1d& source) : Contractable<CGF1d>(source),
        deriv_(source.deriv_), unit_(source.unit_)
    {
      qn_[0] = source.qn_[0];
    }
    explicit CGF1d(const ConstructablePolymorphically& sptr) :
      Contractable<CGF1d>(dynamic_cast<const CGF1d&>(sptr))
    {
      const CGF1d& sptr_cast = dynamic_cast<const CGF1d&>(sptr);
      qn_[0] = sptr_cast.qn_[0];
      deriv_ = sptr_cast.deriv_;
      unit_ = sptr_cast.unit_;
    }
    virtual ~CGF1d() {
    }
 
    CGF1d& operator=(const CGF1d& source)
    {
      qn_[0] = source.qn_[0];
      deriv_ = source.deriv_;
      unit_ = source.unit_;
      Contractable<CGF1d>::operator=(source);
      if (!source.valid()) invalidate();
      return *this;
    }
 
    CGF1d operator+(const CGF1d& B) const {
      //assert(this->is_unit() == false && B.is_unit() == false);
      CGF1d<Axis> Sum(*this);
      Sum.inc(0,B.qn(0));
      Sum.deriv_ += B.deriv_;
      return Sum;
    }
    CGF1d operator-(const CGF1d& B) const {
      //assert(A.is_unit() == false && B.is_unit() == false);
      CGF1d Diff(*this);
      Diff.dec(0,B.qn(0));
      Diff.deriv_ -= B.deriv_;
      return Diff;
    }
 
    const OriginDerivative<1u>& deriv() const { return deriv_; }
    OriginDerivative<1u>& deriv() { return deriv_; }
 
    std::string label() const override {
      // unit *functions* are treated as regular qn-0 functions so that (00|00)^(m) = (unit 0|00)^(m)
      std::ostringstream oss;
      oss << to_string(Axis) << qn_[0];
      if (!deriv_.zero()) oss << "_" << deriv_.label();
 
      // I don't handle labels of contracted CGF1d because I don't think I need them
      // make sure just in case
      assert(!this->contracted());
 
      return oss.str();
    }
 
 
    unsigned int num_bf() const override { return 1; }
    unsigned int qn(unsigned int dir = 0) const {
      assert(dir == 0);
      return qn_[0];
    }
    unsigned int operator[](unsigned int dir) const {
      return this->qn(dir);
    }
 
    bool operator==(const CGF1d& a) const {
      return ( qn_[0] == a.qn_[0] &&
               this->contracted() == a.contracted() &&
               deriv_ == a.deriv_ &&
               unit_ == a.unit_);
    }
 
    void inc(unsigned int dir, unsigned int c = 1u) override {
      assert(!is_unit());
      assert(dir==0);
      if (valid())
        qn_[0] += c;
    }
    void dec(unsigned int dir, unsigned int c = 1u) override {
      if (is_unit()) { invalidate(); return; }
      assert(dir==0);
      if (valid()) {
        if (qn_[0] < c) {
          invalidate();
          return;
        }
        qn_[0] -= c;
      }
    }
    unsigned int norm() const override { return qn_[0]; }
    LIBINT2_UINT_LEAST64 key() const override {
      if (is_unit()) return max_key-1;
      const LIBINT2_UINT_LEAST64 result =
            ( deriv().key() * 2ul +
              (this->contracted() ? 1ul : 0ul)
            ) * max_num_qn +
            qn_[0];
      if (result >= max_key-1) {
        this->print(std::cout);
        std::cout << "result,max_key-1 = " << result << "," << max_key-1 << std::endl;
        assert(result < max_key-1);
      }
      return result;
    }
    const static LIBINT2_UINT_LEAST64 max_num_qn = CGShell::max_qn+1;
    // deriv_key_range = OriginDerivative<1u>::max_key
    // contracted = 2 (yes or no)
    // qn_range = max_num_qn
    // +1 to account for unit function
    const static LIBINT2_UINT_LEAST64 max_key = OriginDerivative<1u>::max_key * OriginDerivative<1u>::max_key * max_num_qn + 1;
 
    void print(std::ostream& os = std::cout) const {
      os << "CGF1d<" << to_string(Axis) << ">: " << label() << std::endl;
    }
 
    static CGF1d unit() {
      CGF1d result;
      result.unit_ = true;
      result.uncontract();
      return result;
    }
    bool is_unit() const { return unit_; }
 
  private:
    static std::array<LIBINT2_UINT_LEAST64, CGShell::max_qn+1> key_l_offset;
  };
 
//  template <CartesianAxis Axis>
//  inline CGF1d<Axis> operator+(const CGF1d<Axis>& A, const CGF1d<Axis>& B) {
//    assert(A.is_unit() == false && B.is_unit() == false);
//    CGF1d<Axis> Sum(A);
//    Sum.inc(0,B.qn(0));
//    Sum.deriv_ += B.deriv_;
//    return Sum;
//  }
//  template <CartesianAxis Axis>
//  inline CGF1d<Axis> operator-(const CGF1d<Axis>& A, const CGF1d<Axis>& B) {
//    //assert(A.is_unit() == false && B.is_unit() == false);
//    CGF1d<Axis> Diff(A);
//    Diff.dec(0,B.qn(0));
//    Diff.deriv_ -= B.deriv_;
//
//    return Diff;
//  }
 
 
  template <CartesianAxis Axis>
  class CGShell1d : public IncableBFSet, public Hashable<LIBINT2_UINT_LEAST64,ComputeKey>,
                    public Contractable< CGShell1d<Axis> > {
 
    unsigned int qn_[1];
    OriginDerivative<1u> deriv_;
    bool unit_; //< if true, this is a unit Gaussian (exponent = 0)
 
  public:
 
    static constexpr auto axis = Axis;
 
    typedef CGF1d<Axis> iter_type;
    typedef IncableBFSet parent_type;
 
    CGShell1d() : unit_(false) { qn_[0] = 0; }
    CGShell1d(unsigned int qn) : unit_(false) { qn_[0] = qn; }
    CGShell1d(unsigned int qn[1]) : unit_(false) { qn_[0] = qn[0]; }
    CGShell1d(const CGShell1d& source) : Contractable<CGShell1d>(source),
        deriv_(source.deriv_), unit_(source.unit_)
    {
      qn_[0] = source.qn_[0];
    }
    virtual ~CGShell1d() {
    }
 
    CGShell1d& operator=(const CGShell1d& source)
    {
      qn_[0] = source.qn_[0];
      deriv_ = source.deriv_;
      unit_ = source.unit_;
      Contractable<CGShell1d>::operator=(source);
      if (!source.valid()) invalidate();
      return *this;
    }
 
    const OriginDerivative<1u>& deriv() const { return deriv_; }
    OriginDerivative<1u>& deriv() { return deriv_; }
 
    std::string label() const override {
      // unit *functions* are treated as regular qn-0 functions so that (00|00)^(m) = (unit 0|00)^(m)
      std::ostringstream oss;
      auto axis_label = to_string(Axis);
      axis_label[0] = std::toupper(axis_label[0]);
      oss << axis_label << qn_[0];
      if (!deriv_.zero()) oss << "_" << deriv_.label();
 
      // I don't handle labels of contracted CGF1d because I don't think I need them
      // make sure just in case
      assert(!this->contracted());
 
      return oss.str();
    }
 
 
    unsigned int num_bf() const override { return qn_[0]+1; }
    unsigned int qn(unsigned int dir=0) const {
      assert(dir == 0);
      return qn_[0];
    }
 
    bool operator==(const CGShell1d& a) const {
      return ( qn_[0] == a.qn_[0] &&
               this->contracted() == a.contracted() &&
               deriv_ == a.deriv_ &&
               unit_ == a.unit_);
    }
 
    void inc(unsigned int dir, unsigned int c = 1u) override {
      assert(false);
    }
    void dec(unsigned int dir, unsigned int c = 1u) override {
      assert(false);
    }
    unsigned int norm() const override { return qn_[0]; }
    LIBINT2_UINT_LEAST64 key() const override {
      if (is_unit()) return max_key-1;
      const LIBINT2_UINT_LEAST64 result =
            ( deriv().key() * 2ul +
              (this->contracted() ? 1ul : 0ul)
            ) * max_num_qn +
            qn_[0];
      if (result >= max_key-1) {
        this->print(std::cout);
        std::cout << "result,max_key-1 = " << result << "," << max_key-1 << std::endl;
        assert(result < max_key-1);
      }
      return result;
    }
    const static LIBINT2_UINT_LEAST64 max_num_qn = CGShell::max_qn+1;
    // deriv_key_range = OriginDerivative<1u>::max_key
    // contracted = 2 (yes or no)
    // qn_range = max_num_qn
    // +1 to account for unit function
    const static LIBINT2_UINT_LEAST64 max_key = OriginDerivative<1u>::max_key * OriginDerivative<1u>::max_key * max_num_qn + 1;
 
    void print(std::ostream& os = std::cout) const {
      os << "CGShell1d<" << to_string(Axis) << ">: " << label() << std::endl;
    }
 
    static CGShell1d unit() {
      CGShell1d result;
      result.unit_ = true;
      result.uncontract();
      return result;
    }
    bool is_unit() const { return unit_; }
 
  private:
    static std::array<LIBINT2_UINT_LEAST64, CGShell::max_qn+1> key_l_offset;
  };
 
#endif
 
#if 0
  class SHGF; // forward declaration
 
  class SHGShell : public IncableBFSet, public Hashable<unsigned,ReferToKey>,
                   public Contractable<SHGShell> {
 
    unsigned int qn_[1];
    OriginDerivative deriv_;
 
    friend SHGShell operator+(const SHGShell& A, const SHGShell& B);
    friend SHGShell operator-(const SHGShell& A, const SHGShell& B);
 
  public:
    typedef SHGF iter_type;
    typedef IncableBFSet parent_type;
 
    SHGShell();
    SHGShell(unsigned int qn);
    SHGShell(const SHGShell&);
    virtual ~SHGShell();
    SHGShell& operator=(const SHGShell&);
 
    const OriginDerivative& deriv() const { return deriv_; }
    OriginDerivative& deriv() { return deriv_; }
 
    std::string label() const;
    unsigned int num_bf() const { return 2*qn_[0]+1; };
    unsigned int qn(unsigned int m=0) const { return qn_[0]; }
 
    bool operator==(const SHGShell&) const;
 
    void inc(unsigned int xyz, unsigned int c = 1u);
    void dec(unsigned int xyz, unsigned int c = 1u);
    unsigned int norm() const;
    unsigned key() const { return (deriv().key() * 2 + (contracted() ? 1 : 0)) * (max_qn+1) + qn_[0]; }
    const static unsigned max_qn = LIBINT_CARTGAUSS_MAX_AM;
    // The range of keys is [0,max_key]
    const static unsigned max_key = 2 * (max_qn + 1) * OriginDerivative::max_key * 2; // deriv_key_range = 2
                                                                                      // qn_range = max_qn + 1
                                                                                      // puresh_key_range = 2
 
    void print(std::ostream& os = std::cout) const;
 
  };
 
  SHGShell operator+(const SHGShell& A, const SHGShell& B);
  SHGShell operator-(const SHGShell& A, const SHGShell& B);
 
  class SHGF : public IncableBFSet, public Hashable<unsigned,ComputeKey>,
               public Contractable<SHGF> {
 
    unsigned int qn_[3];
    OriginDerivative deriv_;
 
    friend SHGF operator+(const SHGF& A, const SHGF& B);
    friend SHGF operator-(const SHGF& A, const SHGF& B);
 
  public:
    typedef SHGF iter_type;
    typedef IncableBFSet parent_type;
    //typedef typename Hashable<unsigned,ComputeKey>::KeyReturnType KeyReturnType;
 
    SHGF();
    SHGF(unsigned int qn[3]);
    SHGF(const SHGF&);
    SHGF(const ConstructablePolymorphically&);
    virtual ~SHGF();
    SHGF& operator=(const SHGF&);
 
    const OriginDerivative& deriv() const { return deriv_; }
    OriginDerivative& deriv() { return deriv_; }
 
    std::string label() const;
    unsigned int num_bf() const { return 1; };
    unsigned int qn(unsigned int xyz) const;
 
    bool operator==(const SHGF&) const;
 
    void inc(unsigned int xyz, unsigned int c = 1u);
    void dec(unsigned int xyz, unsigned int c = 1u);
    unsigned int norm() const;
    unsigned key() const {
      unsigned nxy = qn_[1] + qn_[2];
      unsigned l = nxy + qn_[0];
      unsigned key = nxy*(nxy+1)/2 + qn_[2];
      return ( deriv().key() * 2 + (contracted() ? 1 : 0)) * max_num_qn + key + key_l_offset.at(l);
    }
    const static unsigned max_num_qn = ((1 + (SHGShell::max_qn+1)) * (2 + (SHGShell::max_qn+1)) * (3 + (SHGShell::max_qn+1)) /6);
    const static unsigned max_key = 2 * OriginDerivative::max_key * max_num_qn;
 
    void print(std::ostream& os = std::cout) const;
 
  private:
    static unsigned key_l_offset[SHGShell::max_key+2];
  };
 
  SHGF operator+(const SHGF& A, const SHGF& B);
  SHGF operator-(const SHGF& A, const SHGF& B);
#endif
 
  template <class T>
    struct TrivialBFSet;
  template <>
    struct TrivialBFSet<CGShell> {
      static const bool result = false;
    };
  template <>
    struct TrivialBFSet<CGF> {
      static const bool result = true;
    };
  template <CartesianAxis Axis>
    struct TrivialBFSet< CGShell1d<Axis> > {
      static const bool result = false;
    };
  template <CartesianAxis Axis>
    struct TrivialBFSet< CGF1d<Axis> > {
      static const bool result = true;
    };
#if 0
  template <>
    struct TrivialBFSet<SHGShell> {
      static const bool result = false;
    };
  template <>
    struct TrivialBFSet<SHGF> {
      static const bool result = true;
    };
#endif
 
};
 
#endif