src/bin/libint/braket.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_braket_h_
#define _libint2_src_bin_libint_braket_h_
 
#include <polyconstr.h>
#include <bfset.h>
#include <algebra.h>
#include <libint2/util/intrinsic_types.h>
#include <global_macros.h>
 
namespace libint2 {
 
 
  template <class BFS, unsigned int NP> class ArrayBraket {
  public:
    typedef ArrayBraket<BFS,NP> this_type;
    typedef struct{} parent_type;
    typedef ArrayBraket<typename BFS::iter_type,NP> iter_type;
    static constexpr auto num_particles = NP;
    static constexpr auto num_bf = NP;
 
    typedef BFS bfs_type;
    typedef bfs_type& bfs_ref;
    typedef const BFS& bfs_cref;
 
    ArrayBraket();
    ArrayBraket(const BFS* braket);
    ArrayBraket(const std::vector<std::vector<BFS> >& braket);
    ~ArrayBraket();
 
    bool operator==(const this_type&) const;
    bfs_cref member(unsigned int p, unsigned int i=0) const;
    bfs_ref member(unsigned int p, unsigned int i=0);
    SubIterator* member_subiter(unsigned int p, unsigned int i=0) const;
    void set_member(const BFS&, unsigned int p, unsigned int i=0);
    void set_member(const ConstructablePolymorphically&, unsigned int p, unsigned int i=0);
    unsigned int num_members(unsigned int p) const { assert(p<NP); return 1; }
    unsigned int num_part() const { return NP; }
    inline LIBINT2_UINT_LEAST64 key() const;
    LIBINT2_UINT_LEAST64 max_key() const;
#if COMPUTE_SIZE_DIRECTLY
    unsigned int size() const;
#endif
 
  private:
    BFS bfs_[NP];
 
    // this function resets all cached values
    void reset_cache();
#if COMPUTE_SIZE_DIRECTLY
    mutable unsigned int size_;
#endif
  };
 
  template <class BFS, unsigned int NP>
  ArrayBraket<BFS,NP>::ArrayBraket()
  {
#if COMPUTE_SIZE_DIRECTLY
    size_ = 0;
#endif
  }
 
  template <class BFS, unsigned int NP>
  ArrayBraket<BFS,NP>::ArrayBraket(const BFS* braket) {
    for(int i=0; i<NP; i++)
      bfs_[i] = braket[i];
#if COMPUTE_SIZE_DIRECTLY
    size_ = 0;
#endif
  }
 
  template <class BFS, unsigned int NP>
  ArrayBraket<BFS,NP>::ArrayBraket(const std::vector<std::vector<BFS> >& braket) {
    assert(braket.size()==NP);
    for(unsigned int i=0; i<NP; i++) {
      assert(braket[i].size()==1);
      bfs_[i] = braket[i][0];
    }
#if COMPUTE_SIZE_DIRECTLY
    size_ = 0;
#endif
  }
 
  template <class BFS, unsigned int NP>
  ArrayBraket<BFS,NP>::~ArrayBraket() {}
 
  template <class BFS, unsigned int NP>
  bool
  ArrayBraket<BFS,NP>::operator==(const ArrayBraket<BFS,NP>& a) const {
    for(int i=0; i<NP; i++)
      if (bfs_[i] != a[i])
        return false;
    return true;
  }
 
  template <class BFS, unsigned int NP>
  typename ArrayBraket<BFS,NP>::bfs_cref
  ArrayBraket<BFS,NP>::member(unsigned int p, unsigned int i) const {
    assert(i==0 && p<NP);
    return bfs_[p];
  }
 
  template <class BFS, unsigned int NP>
  typename ArrayBraket<BFS,NP>::bfs_ref
  ArrayBraket<BFS,NP>::member(unsigned int p, unsigned int i) {
    assert(i==0 && p<NP);
    reset_cache();
    return bfs_[p];
  }
 
  template <class BFS, unsigned int NP>
  SubIterator*
  ArrayBraket<BFS,NP>::member_subiter(unsigned int p, unsigned int i) const {
    assert(i==0 && p<NP);
    return static_cast<SubIterator*>(new SubIteratorBase<BFS>( member(p,i) ) );
  }
 
  template <class BFS, unsigned int NP>
  void
  ArrayBraket<BFS,NP>::set_member(const BFS& f, unsigned int p, unsigned int i) {
    assert(i==0 && p<NP);
    reset_cache();
    bfs_[p] = f;
  }
 
  template <class BFS, unsigned int NP>
  void
  ArrayBraket<BFS,NP>::set_member(const ConstructablePolymorphically& f, unsigned int p, unsigned int i) {
    assert(i==0 && p<NP);
    reset_cache();
    bfs_[p] = BFS(f);
  }
 
  template <class BFS, unsigned int NP>
  LIBINT2_UINT_LEAST64
  ArrayBraket<BFS,NP>::key() const {
    LIBINT2_UINT_LEAST64 pfac = 1;
    LIBINT2_UINT_LEAST64 key = 0;
    for(int p=NP-1; p>=0; p--) {
      key += pfac*bfs_[p].key();
      pfac *= BFS::max_key;
    }
    assert(key < this->max_key());
    return key;
  }
 
  template <class BFS, unsigned int NP>
  LIBINT2_UINT_LEAST64
  ArrayBraket<BFS,NP>::max_key() const {
    LIBINT2_UINT_LEAST64 max_key = 1;
    for(int p=NP-1; p>=0; p--) {
      max_key *= BFS::max_key;
    }
    return max_key;
  }
 
#if COMPUTE_SIZE_DIRECTLY
  template <class BFS, unsigned int NP>
  unsigned int
  ArrayBraket<BFS,NP>::size() const {
    if (size_ > 0) return size_;
    size_ = 1;
    if (!TrivialBFSet<BFS>::result) {
      for(int p=NP-1; p>=0; p--) {
        size_ *= bfs_[p].num_bf();
      }
    }
    return size_;
  }
#endif
 
  template <class BFS, unsigned int NP>
  void
  ArrayBraket<BFS,NP>::reset_cache() {
    size_ = 0;
  }
 
  // really need to have typedef template!
  template <typename BFS>
    struct DefaultOnePBraket {
      typedef ArrayBraket<BFS,1> Result;
    };
 
  // really need to have typedef template!
  template <typename BFS>
    struct DefaultTwoPBraket {
      typedef ArrayBraket<BFS,2> Result;
    };
 
 
  enum BraketType {CBra, CKet, PBra, PKet};
  template <class BFS, BraketType BKType> class BraketPair {
  public:
    typedef BraketPair<BFS,BKType> this_type;
    typedef BFS bfs_type;
 
    BraketPair(const BFS& f1, const BFS& f2) : bfs_(f1,f2) {}
    BraketPair(const BraketPair& source) : bfs_(source.bfs_) {}
    BraketPair& operator=(const BraketPair& rhs) { bfs_ = rhs.bfs_; return *this; }
    const BFS& operator[](unsigned int i) const {
      if (i == 0) return bfs_.first;
      if (i == 1) return bfs_.second;
      throw std::logic_error("BraketPair::operator[] -- argument out of range");
    }
    bool operator==(const this_type& rhs) const {
      return rhs.bfs_ == bfs_;
    }
 
  private:
    std::pair<BFS,BFS> bfs_;
  };
 
  namespace braket {
    template <class F> BraketPair<F,PBra> _pbra(const F& f1, const F& f2) {
      return BraketPair<F,PBra>(f1,f2);
    }
    template <class F> BraketPair<F,PKet> _pket(const F& f1, const F& f2) {
      return BraketPair<F,PKet>(f1,f2);
    }
    template <class F> BraketPair<F,CBra> _cbra(const F& f1, const F& f2) {
      return BraketPair<F,CBra>(f1,f2);
    }
    template <class F> BraketPair<F,CKet> _cket(const F& f1, const F& f2) {
      return BraketPair<F,CKet>(f1,f2);
    }
  };
 
  template <class BFS, BraketType BKTypeL,  BraketType BKTypeR>
  algebra::Wedge< BraketPair<BFS,BKTypeL>, BraketPair<BFS,BKTypeR> >
  operator^(const BraketPair<BFS,BKTypeL>& L, const BraketPair<BFS,BKTypeR>& R) {
    return algebra::make_wedge(L,R);
  }
 
};
 
#endif