ref.h

//
// ref.h
//
// Copyright (C) 2005 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_nbody_ref_h
#define _chemistry_qc_nbody_ref_h
 
#include <chemistry/qc/wfn/wfn.h>
#include <chemistry/qc/scf/scf.h>
#include <chemistry/qc/lcao/wfnworld.h>
#include <chemistry/qc/wfn/orbitalspace.h>
#include <chemistry/qc/wfn/spin.h>
 
namespace sc {
 
  class Wavefunction;
 
 
  class PopulatedOrbitalSpace : virtual public SavableState {
    public:
      PopulatedOrbitalSpace(const Ref<OrbitalSpaceRegistry>& oreg,
                            SpinCase1 spin, const Ref<GaussianBasisSet>& bs,
                            const Ref<Integral>& integral,
                            const RefSCMatrix& coefs,
                            const std::vector<double>& occs,
                            const std::vector<ParticleHoleOrbitalAttributes>& active,
                            const RefDiagSCMatrix& energies,
                            bool eorder_increasing = true,
                            Ref<OrbitalSpace> vbs = 0,
                            Ref<FockBuildRuntime> fbrun = 0
                           );
      PopulatedOrbitalSpace(StateIn& si);
      ~PopulatedOrbitalSpace();
      void save_data_state(StateOut& so);
 
      static double zero_occupancy() {
        return 1e-8; // \approx sqrt(DBL_EPSILON)
      }
 
      const Ref<OrbitalSpaceRegistry>& orbital_registry() const { return oreg_; }
 
 
 
      const Ref<OrbitalSpace>& orbs_sb() const { return orbs_sb_; }
      const Ref<OrbitalSpace>& orbs() const { return orbs_; }
      const Ref<OrbitalSpace>& occ_sb() const { return occ_sb_; }
      const Ref<OrbitalSpace>& occ() const { return occ_; }
      const Ref<OrbitalSpace>& occ_act_sb() const { return occ_act_sb_; }
      const Ref<OrbitalSpace>& occ_act() const { return occ_act_; }
      const Ref<OrbitalSpace>& occ_frz_sb() const { return occ_frz_sb_; }
      const Ref<OrbitalSpace>& occ_frz() const { return occ_frz_; }
      const Ref<OrbitalSpace>& uocc_sb() const { return uocc_sb_; }
      const Ref<OrbitalSpace>& uocc() const { return uocc_; }
      const Ref<OrbitalSpace>& uocc_act_sb() const { return uocc_act_sb_; }
      const Ref<OrbitalSpace>& uocc_act() const { return uocc_act_; }
 
 
    private:
      Ref<OrbitalSpaceRegistry> oreg_;
      Ref<OrbitalSpace> orbs_sb_;
      Ref<OrbitalSpace> orbs_;
      Ref<OrbitalSpace> occ_sb_;
      Ref<OrbitalSpace> occ_act_sb_;
      Ref<OrbitalSpace> occ_frz_sb_;
      Ref<OrbitalSpace> occ_;
      Ref<OrbitalSpace> occ_act_;
      Ref<OrbitalSpace> occ_frz_;
      Ref<OrbitalSpace> uocc_sb_;
      Ref<OrbitalSpace> uocc_act_sb_;
      Ref<OrbitalSpace> uocc_;
      Ref<OrbitalSpace> uocc_act_;
 
      void purge();
  };
 
  class RefWavefunction : virtual public SavableState {
    protected:
 
      bool use_world_dfinfo() const { return use_world_dfinfo_; }
 
      RefWavefunction(const Ref<KeyVal>& kv);
 
      RefWavefunction(StateIn&);
      RefWavefunction(const Ref<WavefunctionWorld>& world,
                      const Ref<GaussianBasisSet>& basis,
                      const Ref<Integral>& integral,
                      bool use_world_df = false);
 
    public:
    ~RefWavefunction();
    void save_data_state(StateOut&);
 
    virtual void obsolete();
    virtual void purge();
 
    const Ref<WavefunctionWorld>& world() const { return world_; }
    const Ref<GaussianBasisSet>& basis() const { return basis_; }
    const Ref<Integral>& integral() const { return integral_; }
    virtual const Ref<GaussianBasisSet>& uocc_basis() const { return basis(); }
    const Ref<OrbitalSpace>& valence_orbs() const;
 
    virtual double energy() =0;
    void set_desired_value_accuracy(double);
    virtual double actual_value_accuracy () const =0;
    virtual double desired_value_accuracy() const =0;
    virtual bool desired_value_accuracy_set_to_default() const;
    virtual int nelectron() const =0;
    virtual double magnetic_moment() const =0;
    bool spin_polarized() const {
      return magnetic_moment() != 0.0;
    }
    virtual int dk() const =0;
    virtual Ref<GaussianBasisSet> momentum_basis() const =0;
    virtual RefSymmSCMatrix core_hamiltonian_for_basis(const Ref<GaussianBasisSet> &basis,
                                                       const Ref<GaussianBasisSet> &p_basis) =0;
    virtual RefSymmSCMatrix ordm(SpinCase1 spin) const =0;
    virtual RefSymmSCMatrix ordm_orbs_sb(SpinCase1 spin) const;
    virtual RefSymmSCMatrix ordm_occ_sb(SpinCase1 spin) const;
 
    virtual bool sdref() const =0;
    virtual Ref<DensityFittingInfo> dfinfo() const =0;
 
    const Ref<OrbitalSpace>& oso_space() const;
    const Ref<OrbitalSpace>& orbs_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& orbs(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ_act_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ_frz_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ_act(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& occ_frz(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& uocc_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& uocc_act_sb(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& uocc(SpinCase1 spin = AnySpinCase1) const;
    const Ref<OrbitalSpace>& uocc_act(SpinCase1 spin = AnySpinCase1) const;
 
    virtual void print(std::ostream& os = ExEnv::out0()) const =0;
 
    protected:
 
    mutable Ref<PopulatedOrbitalSpace> spinspaces_[NSpinCases1];
 
    bool force_average_AB_rdm1_;
 
    void init() const;
    virtual void reset();
    virtual void init_spaces() = 0;
 
    private:
    Ref<WavefunctionWorld> world_;   // who owns this?
    Ref<GaussianBasisSet> basis_;
    Ref<Integral> integral_;
    bool use_world_dfinfo_; 
    Ref<OrbitalSpace> valence_orbs_;
 
    virtual void _set_desired_value_accuracy(double eps) =0;
 
  };
 
  class SD_RefWavefunction : public RefWavefunction {
    public:
 
      SD_RefWavefunction(const Ref<KeyVal>& kv);
 
      SD_RefWavefunction(const Ref<WavefunctionWorld>& world,
                            const Ref<OneBodyWavefunction>& obwfn,
                            bool spin_restricted = true,
                            unsigned int nfzc = 0,
                            unsigned int nfzv = 0,
                            Ref<OrbitalSpace> vir_space = 0,
                            std::string occ_orbs = std::string("canonical"));
      SD_RefWavefunction(StateIn&);
      ~SD_RefWavefunction();
      void save_data_state(StateOut&);
 
      void print(std::ostream& os = ExEnv::out0()) const;
 
      void purge();
 
      bool sdref() const { return true; }
      const Ref<OneBodyWavefunction>& obwfn() const { return obwfn_; }
      const Ref<OrbitalSpace>& vir_space() const { return vir_space_; }
      const Ref<GaussianBasisSet>& uocc_basis() const {
        if (vir_space_) return vir_space_->basis();
        else return this->basis();
      }
 
      double energy() { return obwfn()->energy(); }
      double actual_value_accuracy () const { return obwfn()->actual_value_accuracy(); }
      double desired_value_accuracy() const { return obwfn()->desired_value_accuracy(); }
      bool desired_value_accuracy_set_to_default() const { return obwfn()->desired_value_accuracy_set_to_default(); }
      int nelectron() const { return obwfn()->nelectron(); }
      double magnetic_moment() const { return obwfn()->magnetic_moment(); }
      bool spin_restricted() const { return spin_restricted_; }
      int dk() const { return obwfn()->dk(); }
      Ref<GaussianBasisSet> momentum_basis() const { return obwfn()->momentum_basis(); }
      RefSymmSCMatrix core_hamiltonian_for_basis(const Ref<GaussianBasisSet> &basis,
                                                 const Ref<GaussianBasisSet> &p_basis);
      unsigned int nfzc() const { return nfzc_; }
      unsigned int nfzv() const { return nfzv_; }
      RefSymmSCMatrix ordm(SpinCase1 spin) const;
      Ref<DensityFittingInfo> dfinfo() const;
      const std::string& occ_orbitals() const { return occ_orbitals_; }
    private:
      Ref<OneBodyWavefunction> obwfn_;
      Ref<OrbitalSpace> vir_space_;
      bool spin_restricted_;
      std::string occ_orbitals_;
      unsigned int nfzc_;
      unsigned int nfzv_;
      void init_spaces();
      void init_spaces_restricted();
      void init_spaces_unrestricted();
      void _set_desired_value_accuracy(double eps) { obwfn_->set_desired_value_accuracy(eps); }
  };
 
  class Extern_RefWavefunction : public RefWavefunction {
    public:
      Extern_RefWavefunction(const Ref<WavefunctionWorld>& world,
                  const Ref<GaussianBasisSet>& basis,
                  const Ref<Integral>& integral,
                  const RefSCMatrix& orbs,
                  const std::vector<unsigned int>& orbsymm,
                  const RefSymmSCMatrix& alpha_1rdm,
                  const RefSymmSCMatrix& beta_1rdm,
                  unsigned int nocc,
                  unsigned int nfzc = 0,
                  unsigned int nfzv = 0,
                  bool omit_uocc = false);
      Extern_RefWavefunction(const Ref<WavefunctionWorld>& world,
                  const Ref<GaussianBasisSet>& basis,
                  const Ref<Integral>& integral,
                  const RefSCMatrix& orbs,
                  const std::vector<unsigned int>& orbsymm,
                  const RefSymmSCMatrix& alpha_1rdm,
                  const RefSymmSCMatrix& beta_1rdm,
                  std::vector<unsigned int> occpi,
                  std::vector<unsigned int> fzcpi,
                  std::vector<unsigned int> fzvpi,
                  std::vector<unsigned int> holepi = std::vector<unsigned int>(),
                  std::vector<unsigned int> partpi = std::vector<unsigned int>(),
                  bool omit_uocc = false);
      Extern_RefWavefunction(StateIn&);
      virtual ~Extern_RefWavefunction();
      void save_data_state(StateOut&);
      RefSymmSCMatrix ordm(SpinCase1 spin) const { return rdm_[spin]; }
 
      void print(std::ostream& os = ExEnv::out0()) const;
 
      void obsolete();
 
      bool sdref() const;
      double energy() { return 0.0; }
      double actual_value_accuracy () const { return DBL_EPSILON; }
      double desired_value_accuracy() const { return DBL_EPSILON; }
      int nelectron() const { return nelectron_; }
      bool spin_polarized() const { return rdm_[Alpha] != rdm_[Beta]; }
      double magnetic_moment() const;
      bool spin_restricted() const { return true; }
      int dk() const { return 0; }
      Ref<GaussianBasisSet> momentum_basis() const { return 0; }
      RefSymmSCMatrix core_hamiltonian_for_basis(const Ref<GaussianBasisSet> &basis,
                                                 const Ref<GaussianBasisSet> &p_basis);
      unsigned int nfzc() const { return nfzc_; }
      unsigned int nfzv() const { return nfzv_; }
      bool omit_uocc() const { return omit_uocc_; }
      bool ordm_idempotent() const { return ordm_idempotent_; }
    private:
      RefSymmSCMatrix rdm_[NSpinCases1];
      unsigned int nfzc_;
      unsigned int nfzv_;
      unsigned int nelectron_;
      unsigned int magmom_;
      bool omit_uocc_;
      mutable bool ordm_idempotent_;
 
      void init(const RefSCMatrix& orbs,
                const std::vector<unsigned int>& orbsym,
                std::vector<unsigned int> occpi,
                std::vector<unsigned int> fzcpi,
                std::vector<unsigned int> fzvpi,
                std::vector<unsigned int> holepi,
                std::vector<unsigned int> partpi);
 
      void init_spaces() {  throw sc::ProgrammingError("For Extern_RefWavefunction, spaces must be init-ed in constructor");}
      //void init_spaces(unsigned int nocc, const RefSCMatrix& orbs,
      //                 const std::vector<unsigned int>& orbsym);
      void init_spaces(const RefSCMatrix& orbs,
                       const std::vector<unsigned int>& orbsym,
                       std::vector<unsigned int> occpi,
                       std::vector<unsigned int> fzcpi,
                       std::vector<unsigned int> fzvpi,
                       std::vector<unsigned int> holepi,
                       std::vector<unsigned int> partpi);
      void _set_desired_value_accuracy(double eps) {
        // do nothing
      }
 
      Ref<DensityFittingInfo> dfinfo() const;
  };
 
  struct RefWavefunctionFactory {
      static Ref<RefWavefunction>
        make(const Ref<WavefunctionWorld>& world,
             const Ref<Wavefunction>& ref,
             bool spin_restricted = true,
             unsigned int nfzc = 0,
             unsigned int nfzv = 0,
             Ref<OrbitalSpace> vir_space = 0);
  };
 
  Ref<OrbitalSpace>
  compute_canonvir_space(const Ref<FockBuildRuntime>& fb_rtime,
                         const Ref<OrbitalSpace>& A,
                         SpinCase1 spin);
 
  void add_ao_space(const Ref<GaussianBasisSet>& bs,
                    const Ref<Integral>& ints,
                    const Ref<AOSpaceRegistry>& aoreg,
                    const Ref<OrbitalSpaceRegistry> oreg);
  void remove_ao_space(const Ref<GaussianBasisSet>& bs,
                       const Ref<AOSpaceRegistry>& aoreg,
                       const Ref<OrbitalSpaceRegistry> oreg);
 
  // end of addtogroup ChemistryElectronicStructureOneBody
 
};
 
#endif