src/lib/chemistry/qc/mbptr12/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.
//

#ifdef __GNUG__
#pragma interface
#endif

#ifndef _chemistry_qc_mbptr12_orbitalspace_h
#define _chemistry_qc_mbptr12_orbitalspace_h

#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/mbptr12/spin.h>
#include <chemistry/qc/mbptr12/registry.h>

using namespace std;

namespace sc {

  template <typename Attributes>
  class DecoratedOrbital {
    public:
      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;

  namespace detail {
    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_;
    };
  } // namespace detail

  typedef DecoratedOrbital< detail::MolecularOrbitalAttributes > MolecularOrbital;
  typedef DecoratedOrbital< detail::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_;
  };

  struct EnergyMOOrder {
    public:
      bool operator()(const MolecularOrbital& o1, const MolecularOrbital& o2) const {
        if (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 1;
      }
      unsigned int nblocks() const {
        return 1;
      }
  };

  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();
          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_; }

    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;
      const RefDiagSCMatrix& evals() const;
      const vector<unsigned int>& orbsym() const;
      unsigned int rank() const;
      unsigned int nblocks() const;
      const 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);

  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:
      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:
      typedef enum {
        Fock = 0,
        Coulomb = 1,
        Exchange = 2,
        KineticEnergy = 3,
        CoreHamiltonian = 4,
        CorePlusCoulomb = 5,
        InvalidOneBodyOperator = 6
      } OneBodyOperator;
#define ParsedTransformedOrbitalSpaceKey_OneBodyOperatorLabelInitializer {"F", "J", "K", "T", "h", "h+J"}
      static const char* OneBodyOperatorLabel[]; // = ParsedTransformedOrbitalSpaceKey_OneBodyOperatorLabelInitializer

      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& original_label() const {
        return original_label_;
      }
      SpinCase1 original_spin() const {
        return spin_;
      }
      OneBodyOperator transform_operator() const {
        return transform_operator_;
      }

      static std::string key(const std::string& label, SpinCase1 spin,
                             const std::string& original_label,
                             SpinCase1 original_spin, OneBodyOperator oper);

    private:
      std::string key_;
      std::string label_;
      SpinCase1 spin_;
      std::string original_label_;
      SpinCase1 original_spin_;
      OneBodyOperator transform_operator_;
  };

  typedef Registry<std::string, Ref<OrbitalSpace> ,
      detail::SingletonCreationPolicy> OrbitalSpaceRegistry;
  std::pair<std::string, Ref<OrbitalSpace> > make_keyspace_pair(const Ref<
      OrbitalSpace>& space, SpinCase1 spin = AnySpinCase1);

  typedef Registry<Ref<GaussianBasisSet> , Ref<OrbitalSpace> ,
      detail::SingletonCreationPolicy> 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&);

  };

}

#include <chemistry/qc/mbptr12/orbitalspace.timpl.h>

#endif

// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End:

---