orbitalspace.h

//
// orbitalspace.h
//
// Copyright (C) 2004 Edward Valeev
//
// Author: Edward Valeev <evaleev@vt.edu>
// Maintainer: EV
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit 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 Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
 
#ifndef _chemistry_qc_mbptr12_orbitalspace_h
#define _chemistry_qc_mbptr12_orbitalspace_h
 
#include <bitset>
#include <vector>
#include <stdexcept>
#include <algorithm>
#include <util/ref/ref.h>
#include <util/state/statein.h>
#include <util/state/stateout.h>
#include <util/group/thread.h>
#include <math/scmat/abstract.h>
#include <util/state/statein.h>
#include <chemistry/qc/basis/basis.h>
#include <chemistry/qc/basis/operator.h>
#include <chemistry/qc/wfn/spin.h>
#include <util/misc/registry.h>
 
namespace sc {
 
 
  template <typename Attributes>
  class DecoratedOrbital {
    public:
      DecoratedOrbital() {}
      DecoratedOrbital(size_t index, const Attributes& attr) : index_(index), attr_(attr) {}
 
      size_t index() const { return index_; }
      const Attributes& attr() const { return attr_; }
 
    private:
      size_t index_;
      Attributes attr_;
  };
 
  typedef DecoratedOrbital< unsigned int > BlockedOrbital;
 
    struct MolecularOrbitalAttributes {
      public:
        MolecularOrbitalAttributes(unsigned int irrep, double energy,
                                   double occnum) :
          irrep_(irrep), energy_(energy), occnum_(occnum) {
        }
 
        unsigned int irrep() const {
          return irrep_;
        }
        double energy() const {
          return energy_;
        }
        double occnum() const {
          return occnum_;
        }
 
      private:
        unsigned int irrep_;
        double energy_;
        double occnum_;
    };
    struct MolecularSpinOrbitalAttributes : public MolecularOrbitalAttributes {
      public:
        MolecularSpinOrbitalAttributes(unsigned int irrep,
                                       double energy,
                                       double occnum,
                                       SpinCase1 spin) :
          MolecularOrbitalAttributes(irrep, energy, occnum), spin_(spin)
          {}
 
        SpinCase1 spin() const { return spin_; }
 
      private:
        SpinCase1 spin_;
    };
 
 
  struct ParticleHoleOrbitalAttributes : public std::bitset<2> {
    ParticleHoleOrbitalAttributes() {}
    ParticleHoleOrbitalAttributes(unsigned long val) : std::bitset<2>(val) {}
    ParticleHoleOrbitalAttributes(const std::bitset<2>& bs) : std::bitset<2>(bs) {}
 
    static ParticleHoleOrbitalAttributes Hole;          
    static ParticleHoleOrbitalAttributes Particle;      
    static ParticleHoleOrbitalAttributes Any;           
    static ParticleHoleOrbitalAttributes None;          
  };
 
  typedef DecoratedOrbital< MolecularOrbitalAttributes > MolecularOrbital;
  typedef DecoratedOrbital< MolecularSpinOrbitalAttributes > MolecularSpinOrbital;
 
  struct SymmetryMOOrder {
    public:
      SymmetryMOOrder(unsigned int nirreps) : nirreps_(nirreps) {}
 
      bool operator()(const MolecularOrbital& o1, const MolecularOrbital& o2) const {
        if (o1.attr().irrep() < o2.attr().irrep())
          return true;
        else if (o1.attr().irrep() == o2.attr().irrep()) {
          if (o1.attr().energy() < o2.attr().energy())
            return true;
          else if (o1.attr().energy() == o2.attr().energy()) {
            if (o1.attr().occnum() < o2.attr().occnum())
              return true;
          }
        }
        return false;
      }
      unsigned int block(const MolecularOrbital& o) const {
        return o.attr().irrep();
      }
      unsigned int nblocks() const {
        return nirreps_;
      }
 
    private:
      unsigned int nirreps_;
  };
 
  template <typename EnergyCompare = std::less<double> > struct EnergyMOOrder {
    public:
      bool operator()(const MolecularOrbital& o1, const MolecularOrbital& o2) const {
        if ( ecompare(o1.attr().energy(), o2.attr().energy()) )
          return true;
        else if (o1.attr().energy() == o2.attr().energy()) {
          if (o1.attr().irrep() < o2.attr().irrep())
            return true;
        }
        return false;
      }
      unsigned int block(const MolecularOrbital& o) const {
        return 0;
      }
      unsigned int nblocks() const {
        return 1;
      }
    private:
      EnergyCompare ecompare;
  };
 
  struct CorrelatedMOOrder {
    public:
      CorrelatedMOOrder(unsigned int nirreps) : nirreps_(nirreps) {}
 
      bool operator()(const MolecularOrbital& o1, const MolecularOrbital& o2) const {
        // occupieds come before virtuals
        if (o1.attr().occnum() > o2.attr().occnum())
          return true;
        else if (o1.attr().occnum() == o2.attr().occnum()) {
          if (o1.attr().irrep() < o2.attr().irrep())
            return true;
          else if (o1.attr().irrep() == o2.attr().irrep()) {
            if (o1.attr().energy() < o2.attr().energy())
              return true;
          }
        }
        return false;
      }
 
      unsigned int block(const MolecularOrbital& o) const {
        const unsigned int irrep = o.attr().irrep();
        const double occnum = o.attr().occnum();
        // occupieds come before virtuals
        const int occblock = (occnum == 1.0) ? 0 : 1;
        return occblock * nirreps_ + irrep;
      }
      unsigned int nblocks() const {
        return nirreps_ * 2;
      }
 
    private:
      unsigned int nirreps_;
  };
 
  struct CorrelatedSpinMOOrder {
    public:
      CorrelatedSpinMOOrder(unsigned int nirreps) : nirreps_(nirreps) {}
 
      bool operator()(const MolecularSpinOrbital& o1, const MolecularSpinOrbital& o2) const {
        // occupieds come before virtuals
        if (o1.attr().occnum() > o2.attr().occnum())
          return true;
        else if (o1.attr().occnum() == o2.attr().occnum()) {
          if (o1.attr().spin() < o2.attr().spin())
            return true;
          else if (o1.attr().spin() == o2.attr().spin()) {
            if (o1.attr().irrep() < o2.attr().irrep())
              return true;
            else if (o1.attr().irrep() == o2.attr().irrep()) {
              if (o1.attr().energy() < o2.attr().energy())
                return true;
            }
          }
        }
        return false;
      }
 
      unsigned int block(const MolecularSpinOrbital& o) const {
        const unsigned int irrep = o.attr().irrep();
        const SpinCase1 spin = o.attr().spin();
        const unsigned int spincase = (spin == Alpha) ? 0 : 1;
        const double occnum = o.attr().occnum();
        const unsigned int occblock = (occnum == 1.0) ? 0 : 1;
        // occupieds come before virtuals, Alpha before Beta
        const unsigned int result = (occblock * 2 + spincase) * nirreps_ + irrep;
        return result;
      }
      unsigned int nblocks() const {
        return nirreps_ * 4;
      }
 
    private:
      unsigned int nirreps_;
  };
 
  namespace detail {
 
    template <typename Container> struct ContainerAdaptor {
      public:
        typedef typename Container::value_type value_type;
        ContainerAdaptor(const Container& cont) : cont_(cont) {}
        size_t size() const { return cont_.size(); }
        value_type elem(size_t i) const { return cont_[i]; }
 
      private:
        const Container& cont_;
    };
 
    template<> struct ContainerAdaptor<RefDiagSCMatrix> {
      public:
        typedef RefDiagSCMatrix Container;
        typedef double value_type;
        ContainerAdaptor(const Container& cont) : cont_(cont) {}
        size_t size() const { return cont_.dim().n(); }
        value_type elem(size_t i) const { return cont_(i); }
 
      private:
        const Container& cont_;
    };
 
  } // namespace detail
 
  template <typename Attribute,
            typename AttributeContainer,
            typename Compare = std::less<Attribute> >
  struct MolecularOrbitalMask {
    private:
      typedef DecoratedOrbital<Attribute> MO;
      struct _compare {
        bool operator()(const MO& mo1,
            const MO& mo2) const {
          Compare comp;
          return comp(mo1.attr(), mo2.attr());
        }
      };
 
    public:
      MolecularOrbitalMask(unsigned int n,
                           const AttributeContainer& attributes) :
        mask_(detail::ContainerAdaptor<AttributeContainer>(attributes).size(),true)
        {
          // validate input
          if (n == 0) return;
          const size_t nmos = mask_.size();
          MPQC_ASSERT(n < nmos);
 
          // copy attributes to vector of MOs
          std::vector<MO> mos;
          detail::ContainerAdaptor<AttributeContainer> contadaptor(attributes);
          for(size_t mo=0; mo<nmos; ++mo) {
            mos.push_back(MO(mo,contadaptor.elem(mo)));
          }
 
          // sort
          _compare comp;
          std::stable_sort(mos.begin(), mos.end(), comp);
 
          // copy to mask
          for(unsigned int t=0; t<n; ++t) {
            mask_[ mos[t].index() ] = false;
          }
        }
 
      const std::vector<bool>& mask() const { return mask_; }
 
      bool operator[](size_t o) const { return mask_[o]; }
 
    private:
      std::vector<bool> mask_;
  };
 
 
 
  class OrbitalSpace: virtual public SavableState {
 
    public:
 
      enum IndexOrder {
        symmetry = 0, 
        energy = 1, 
        correlated = 2, 
        general = 3 
      };
 
    private:
 
      static const unsigned int max_id_length = 30; // max length of string returned by id()
      std::string id_; // see documentation for id()
      std::string name_; // String identifier for the orbital space
 
      Ref<GaussianBasisSet> basis_; // The AO basis
      Ref<Integral> integral_;
      RefSCMatrix coefs_; // AO->MO transformation coefficients (nao by nmo matrix)
      RefDiagSCMatrix evals_; // "eigenvalues" associated with the MOs
      RefSCDimension dim_; // only here to allow dim() return const &
      std::vector<unsigned int> orbsym_; // irrep of each orbital
      std::vector<unsigned int> block_offsets_; // Start of each block
      std::vector<unsigned int> block_sizes_; // Number of orbitals in each block
 
      // checks orbsym_ for irrep indices outside of the allowed range
      void check_orbsym() const;
 
      // determines block_sizes_ from nfzc, nfzv, and evals and returns offsets for each block
      std::vector<unsigned int>
          frozen_to_blockinfo(unsigned int nfzc, unsigned int nfzv,
                              const RefDiagSCMatrix& evals);
 
      // computes coefficient matrix from the full coefficient matrix. If moorder == energy
      // then the MO vectors will be sorted by their eigenvalues
      void full_coefs_to_coefs(const RefSCMatrix& full_coefs,
                               const RefDiagSCMatrix& evals,
                               const std::vector<unsigned int>& offsets,
                               IndexOrder moorder);
      void init();
 
      // sorting functions borrowed from mbpt.cc
      static void dquicksort(double *item, unsigned int *index, unsigned int n);
 
    protected:
      OrbitalSpace();
      void init(const std::string& id, const std::string& name,
                const Ref<GaussianBasisSet>& basis, const Ref<Integral>& integral,
                const RefSCMatrix& coefs,
                const RefDiagSCMatrix& evals,
                const std::vector<unsigned int>& orbsyms,
                unsigned int nblocks,
                const std::vector<BlockedOrbital>& indexmap);
 
    public:
      OrbitalSpace(StateIn&);
      OrbitalSpace(const OrbitalSpace&);
      OrbitalSpace(const std::string& id, const std::string& name,
                   const RefSCMatrix& full_coefs, const RefDiagSCMatrix& evals,
                   const Ref<GaussianBasisSet>& basis,
                   const Ref<Integral>& integral,
                   const std::vector<unsigned int>& block_offsets,
                   const std::vector<unsigned int>& block_sizes);
 
      OrbitalSpace(const std::string& id, const std::string& name,
                   const RefSCMatrix& full_coefs,
                   const Ref<GaussianBasisSet>& basis,
                   const Ref<Integral>& integral,
                   const std::vector<unsigned int>& block_offsets,
                   const std::vector<unsigned int>& block_sizes,
                   const IndexOrder& moorder = symmetry,
                   const RefDiagSCMatrix& evals = 0);
      OrbitalSpace(const std::string& id, const std::string& name,
                   const RefSCMatrix& full_coefs,
                   const Ref<GaussianBasisSet>& basis,
                   const Ref<Integral>& integral, const RefDiagSCMatrix& evals,
                   unsigned int nfzc, unsigned int nfzv,
                   const IndexOrder& moorder = energy);
      OrbitalSpace(const std::string& id, const std::string& name,
                   const RefSCMatrix& full_coefs,
                   const Ref<GaussianBasisSet>& basis,
                   const Ref<Integral>& integral);
 
      OrbitalSpace(const std::string& id, const std::string& name,
                   const Ref<OrbitalSpace>& orig_space,
                   const RefSCMatrix& new_coefs,
                   const Ref<GaussianBasisSet>& new_basis);
      ~OrbitalSpace();
 
      void save_data_state(StateOut&);
 
      OrbitalSpace& operator=(const OrbitalSpace& other);
 
      const std::string& name() const;
 
      const std::string& id() const;
      const RefSCDimension& dim() const;
      const Ref<GaussianBasisSet>& basis() const;
      const Ref<Integral>& integral() const;
      const RefSCMatrix& coefs() const;
      RefSCMatrix coefs_nb() const;
      const RefDiagSCMatrix& evals() const;
      const std::vector<unsigned int>& orbsym() const;
      unsigned int rank() const;
      unsigned int nblocks() const;
      const std::vector<unsigned int>& block_sizes() const;
 
      size_t memory_in_use() const;
 
      void print(std::ostream&o = ExEnv::out0()) const;
      void print_summary(std::ostream& os) const;
      void print_detail(std::ostream&o = ExEnv::out0()) const;
 
  };
 
  bool operator==(const OrbitalSpace& space1, const OrbitalSpace& space2);
 
  class CannotConstructMap: public std::logic_error {
    public:
      CannotConstructMap() :
        std::logic_error("Cannot map given OrbitalSpaces") {
      }
  };
 
  typedef std::vector<unsigned int> MOIndexMap;
  MOIndexMap operator<<(const OrbitalSpace& space2, const OrbitalSpace& space1);
 
  std::vector<int> map(const OrbitalSpace& space2, const OrbitalSpace& space1, bool expect_same_bases = true);
 
  typedef std::vector<std::pair<unsigned int, double> > SparseMOIndexMap;
  SparseMOIndexMap sparsemap(const OrbitalSpace& space2,
                             const OrbitalSpace& space1, double hardzero =
                                 1e-12);
 
  RefSCMatrix transform(const OrbitalSpace& space2, const OrbitalSpace& space1,
                        const Ref<SCMatrixKit>& kit =
                            SCMatrixKit::default_matrixkit());
 
  RefSCMatrix overlap(const OrbitalSpace& space2, const OrbitalSpace& space1,
                      const Ref<SCMatrixKit>& kit =
                          SCMatrixKit::default_matrixkit());
 
  bool in(const OrbitalSpace& s1, const OrbitalSpace& s2);
 
  class ParsedOrbitalSpaceKey {
    public:
      struct exception {};
 
      ParsedOrbitalSpaceKey(const std::string& key);
 
      const std::string& key() const {
        return key_;
      }
      const std::string& label() const {
        return label_;
      }
      SpinCase1 spin() const {
        return spin_;
      }
 
      static std::string key(const std::string& label, SpinCase1 spin);
 
    private:
      std::string key_;
      std::string label_;
      SpinCase1 spin_;
  };
 
  class ParsedTransformedOrbitalSpaceKey {
    public:
      struct exception { };
 
      ParsedTransformedOrbitalSpaceKey(const std::string& key);
 
      const std::string& key() const {
        return key_;
      }
      const std::string& label() const {
        return label_;
      }
      SpinCase1 spin() const {
        return spin_;
      }
      const std::string& support_label() const {
        return support_label_;
      }
      SpinCase1 support_spin() const {
        return spin_;
      }
      OneBodyOper::type transform_operator() const {
        return transform_operator_;
      }
 
      static std::string key(const std::string& label, SpinCase1 spin,
                             const std::string& original_label,
                             SpinCase1 original_spin, OneBodyOper::type oper);
 
      static bool valid_key(const std::string& key);
 
    private:
      std::string key_;
      std::string label_;
      SpinCase1 spin_;
      std::string support_label_;
      SpinCase1 support_spin_;
      OneBodyOper::type transform_operator_;
  };
 
  typedef Registry<std::string, Ref<OrbitalSpace> ,
                   detail::NonsingletonCreationPolicy,
                   std::equal_to<std::string>,
                   RefObjectEqual<OrbitalSpace> > OrbitalSpaceRegistry;
  std::pair<std::string, Ref<OrbitalSpace> > make_keyspace_pair(const Ref<
      OrbitalSpace>& space, SpinCase1 spin = AnySpinCase1);
  std::string new_unique_key(const Ref<OrbitalSpaceRegistry>& oreg);
 
  typedef Registry<Ref<GaussianBasisSet>, Ref<OrbitalSpace>,
                   detail::NonsingletonCreationPolicy,
                   std::equal_to< Ref<GaussianBasisSet> >,
                   RefObjectEqual<OrbitalSpace> > AOSpaceRegistry;
 
 
  template <typename Order>
  class OrderedOrbitalSpace : public OrbitalSpace {
    public:
      OrderedOrbitalSpace(const std::string& id, const std::string& name,
                          const Ref<GaussianBasisSet>& basis,
                          const Ref<Integral>& integral,
                          const RefSCMatrix& coefs, const RefDiagSCMatrix& evals,
                          const RefDiagSCMatrix& occnums,
                          const std::vector<unsigned int>& orbsyms,
                          const Order& order);
 
      OrderedOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
      ~OrderedOrbitalSpace();
 
    private:
 
      typedef OrderedOrbitalSpace this_type;
      // ClassDesc object
      static ClassDesc class_desc_;
  };
 
  template <typename Order>
  class OrderedSpinOrbitalSpace : public OrbitalSpace {
    public:
      OrderedSpinOrbitalSpace(const std::string& id, const std::string& name,
                          const Ref<GaussianBasisSet>& basis,
                          const Ref<Integral>& integral,
                          const RefSCMatrix& coefs_alpha,
                          const RefSCMatrix& coefs_beta,
                          const RefDiagSCMatrix& evals_alpha,
                          const RefDiagSCMatrix& evals_beta,
                          const RefDiagSCMatrix& occnums_alpha,
                          const RefDiagSCMatrix& occnums_beta,
                          const std::vector<unsigned int>& orbsyms_alpha,
                          const std::vector<unsigned int>& orbsyms_beta,
                          const Order& order);
 
      OrderedSpinOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
      ~OrderedSpinOrbitalSpace();
 
    private:
 
      typedef OrderedSpinOrbitalSpace this_type;
      // ClassDesc object
      static ClassDesc class_desc_;
  };
 
 
  class MaskedOrbitalSpace : public OrbitalSpace {
    public:
      MaskedOrbitalSpace(const std::string& id, const std::string& name,
                         const Ref<OrbitalSpace>& orig_space,
                         const std::vector<bool>& mask);
 
      MaskedOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
 
  };
 
 
  class NonblockedOrbitalSpace : public OrbitalSpace {
    public:
      NonblockedOrbitalSpace(const std::string& id, const std::string& name,
                       const Ref<OrbitalSpace>& orig_space);
 
      NonblockedOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
 
  };
 
 
  class AtomicOrbitalSpace : public OrbitalSpace {
    public:
      AtomicOrbitalSpace(const std::string& id, const std::string& name,
                         const Ref<GaussianBasisSet>& basis,
                         const Ref<Integral>& integral);
 
      AtomicOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
 
  };
 
 
  class OrbitalSpaceUnion : public OrbitalSpace {
    public:
      OrbitalSpaceUnion(const std::string& id, const std::string& name,
                        const OrbitalSpace& s1, const OrbitalSpace& s2,
                        bool merge_blocks = true);
 
      OrbitalSpaceUnion(StateIn&);
      void save_data_state(StateOut&);
  };
 
 
  class EmptyOrbitalSpace : public OrbitalSpace {
    public:
      EmptyOrbitalSpace(const std::string& id, const std::string& name,
                        const Ref<GaussianBasisSet>& basis,
                        const Ref<Integral>& integral,
                        const IndexOrder& moorder = symmetry
                        );
 
      EmptyOrbitalSpace(StateIn&);
      void save_data_state(StateOut&);
    private:
      static ClassDesc class_desc_;
  };
 
  // end of addtogroup ChemistryElectronicStructureOneBody
 
}
 
#include <chemistry/qc/wfn/orbitalspace.timpl.h>
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: