psiwfn.h

//
// psiwfn.h
//
// Copyright (C) 2002 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_psi_psiwfn_h
#define _chemistry_qc_psi_psiwfn_h
 
#include <chemistry/qc/wfn/wfn.h>
#include <chemistry/qc/wfn/obwfn.h>
#include <chemistry/qc/psi/psiexenv.h>
#include <chemistry/qc/wfn/spin.h>
#include <chemistry/qc/wfn/orbitalspace.h>
#include <chemistry/qc/mbptr12/mbptr12.h>
 
namespace sc {
 
 
  class PsiWavefunction : public Wavefunction {
 
      Ref<PsiExEnv> exenv_;
      int value_implemented() const {
        return 1;
      }
 
    protected:
      Ref<PsiWavefunction> prerequisite_;
      int nirrep_;
      size_t memory_;
      char *memory_str_;
      bool compute_1rdm_;
      char *dertype_;
      virtual void write_input(int conv) =0;
 
      std::vector<unsigned int> read_occ(const Ref<KeyVal> &keyval, const char *name,
                                         size_t nirrep);
 
      int debug() const;
 
    public:
      PsiWavefunction(const Ref<KeyVal>&);
      PsiWavefunction(StateIn&);
      ~PsiWavefunction();
 
      void save_data_state(StateOut&);
 
      static Integral::CartesianOrdering cartesian_ordering();
 
      virtual void write_basic_input(int conv);
      void compute();
      void print(std::ostream&o=ExEnv::out0()) const;
      int nirrep() const { return nirrep_; }
 
      virtual RefSymmSCMatrix density();
      virtual RefSymmSCMatrix mo_density(SpinCase1 spin = AnySpinCase1) =0;
 
      Ref<PsiExEnv> exenv() const {
        return exenv_;
      }
 
      Ref<PsiInput> get_psi_input() const {
        return exenv_->get_psi_input();
      }
 
#if 0
      std::vector< std::pair<unsigned int,unsigned int> > shell_map();
      std::vector<unsigned int> ao_map();
#endif
 
      double nuclear_repulsion_energy();
 
      void obsolete();
      void symmetry_changed();
  };
 
  class PsiSCF;
 
 
  class PsiCorrWavefunction : public PsiWavefunction {
    protected:
      Ref<PsiSCF> reference_;
      unsigned int nfzc_;
      unsigned int nfzv_;
      mutable std::vector<unsigned int> frozen_docc_;
      mutable std::vector<unsigned int> frozen_uocc_;
      RefSymmSCMatrix mo_density_[NSpinCases1];
 
      void write_input(int conv);
      // set fz2restr to use RESTRICTED instead of FROZEN keywords
      void write_input_frozen2restricted(int conv, bool fz2restr);
 
      double valacc_to_refacc() const { return 100.0; }
 
      virtual std::vector<unsigned int> map_density_to_sb();
 
    public:
      PsiCorrWavefunction(const Ref<KeyVal>&);
      PsiCorrWavefunction(StateIn&);
      ~PsiCorrWavefunction();
      void save_data_state(StateOut&);
 
      void print(std::ostream& os) const;
      double magnetic_moment() const;
      void set_desired_value_accuracy(double acc);
 
      void compute(); // compute is overloaded because reference object needs to be marked computed
      void obsolete();
      void symmetry_changed();
      RefSymmSCMatrix density();
 
      const Ref<PsiSCF>& reference() const;
      int nelectron();
 
#if 0   // these should be implemented in PsiCC classes
      const Ref<OrbitalSpace>& occ_act_sb(SpinCase1);
      const Ref<OrbitalSpace>& vir_act_sb(SpinCase1);
#endif
      unsigned int nfzc() const;
      unsigned int nfzv() const;
      virtual const Ref<OrbitalSpace>&  orbs_sb(SpinCase1 spin);
 
      const std::vector<unsigned int>& frozen_docc() const;
      const std::vector<unsigned int>& frozen_uocc() const;
      const std::vector<unsigned int> docc_act();
      const std::vector<unsigned int> socc();
      const std::vector<unsigned int> uocc_act();
      const std::vector<unsigned int> docc();
      const std::vector<unsigned int> uocc();
 
      virtual double reference_energy();
 
      RefSymmSCMatrix mo_density(SpinCase1 spin);
#if 0
      RefSymmSCMatrix onepdm(const SpinCase1 &spin);
      RefSymmSCMatrix onepdm();
      RefSymmSCMatrix twopdm();
      RefSymmSCMatrix twopdm_dirac();
      RefSymmSCMatrix twopdm_dirac_from_components();
#endif
      RefSymmSCMatrix twopdm_dirac(const SpinCase2 &pairspin);
      RefSymmSCMatrix twopdm_dirac();
      void print_onepdm_vec(FILE *output,const RefSCVector &opdm,double TOL);
      void print_onepdm_mat(FILE *output,const RefSymmSCMatrix &opdm,double TOL);
      void print_twopdm_mat(FILE *output,const RefSymmSCMatrix &tpdm, double TOL);
      void print_twopdm_arr(FILE *output,double *tpdm,double TOL);
  };
 
 
 
  class PsiSCF : public PsiWavefunction {
      Ref<OrbitalSpace> orbs_sb_[NSpinCases1];
      RefDiagSCMatrix evals_[NSpinCases1];
      RefSCMatrix coefs_[NSpinCases1];
      RefSymmSCMatrix mo_density_[NSpinCases1];
      std::vector<unsigned int> occpi_[NSpinCases1];
      std::vector<unsigned int> uoccpi_[NSpinCases1];
      std::vector<unsigned int> mopi_;
      std::vector<double> occupation_[NSpinCases1];
      void compute_occupations(SpinCase1 spin);
 
      // PsiCorrWavefunction needs to be able to set value of PsiSCF
      friend void PsiCorrWavefunction::compute();
    protected:
      std::vector<unsigned int> docc_;
      std::vector<unsigned int> socc_;
      int multp_;
      int charge_;
      int maxiter_;
      double diisdamp_;
      double levelshift_;
      bool diis_;
      static const int default_maxiter = 200;
      static double guess_acc_ratio() { return 1e4; }
 
      Ref<OneBodyWavefunction> guess_wfn_;
 
    public:
      PsiSCF(const Ref<KeyVal>&);
      PsiSCF(StateIn&);
      ~PsiSCF();
      void save_data_state(StateOut&);
 
      double magnetic_moment() const;
 
      void import_occupations(const Ref<OneBodyWavefunction>& obwfn);
      enum RefType {rhf, hsoshf, uhf};
      virtual PsiSCF::RefType reftype() const =0;
      virtual const RefDiagSCMatrix& evals(SpinCase1 spin = AnySpinCase1);
      virtual const RefSCMatrix& coefs(SpinCase1 spin = AnySpinCase1);
      const std::vector<unsigned int>& occpi(SpinCase1 S);
      const std::vector<unsigned int>& uoccpi(SpinCase1 S);
      const std::vector<unsigned int>& mopi();
      const Ref<OrbitalSpace>&  orbs_sb(SpinCase1 spin);
      int nelectron();
      double occupation(int mo);
      double alpha_occupation(int mo);
      double beta_occupation(int mo);
 
      RefSymmSCMatrix density();
      RefSymmSCMatrix alpha_density();
      RefSymmSCMatrix beta_density();
 
      RefSymmSCMatrix ao_density(SpinCase1);
      RefSymmSCMatrix ao_density();
      RefSymmSCMatrix mo_density(SpinCase1 spin = AnySpinCase1);
 
      unsigned int nmo();
      unsigned int nocc(SpinCase1 spin);
      unsigned int multiplicity() const { return multp_; }
 
      void obsolete();
      void symmetry_changed();
  };
 
 
  class PsiCLHF : public PsiSCF {
    protected:
      void write_input(int conv);
    public:
      PsiCLHF(const Ref<KeyVal>&);
      PsiCLHF(StateIn&);
      ~PsiCLHF();
      void print(std::ostream& os = ExEnv::out0()) const;
 
      void write_basic_input(int conv);
      bool analytic_gradient_implemented() const {
        return true;
      }
      PsiSCF::RefType reftype() const {
        return rhf;
      }
  };
 
 
  class PsiHSOSHF : public PsiSCF {
    protected:
      void write_input(int conv);
    public:
      PsiHSOSHF(const Ref<KeyVal>&);
      PsiHSOSHF(StateIn&);
      ~PsiHSOSHF();
      void print(std::ostream& os = ExEnv::out0()) const;
 
      void write_basic_input(int conv);
      bool analytic_gradient_implemented() const {
        return true;
      }
      PsiSCF::RefType reftype() const {
        return hsoshf;
      }
 
      const RefSCMatrix& coefs_semicanonical(SpinCase1 s);
      const RefDiagSCMatrix& evals_semicanonical(SpinCase1 s);
 
    private:
      RefSCMatrix coefs_sc_[NSpinCases1];
      RefDiagSCMatrix evals_sc_[NSpinCases1];
 
      void semicanonical();
 
      void obsolete();
  };
 
 
  class PsiUHF : public PsiSCF {
    protected:
      void write_input(int conv);
    public:
      PsiUHF(const Ref<KeyVal>&);
      PsiUHF(StateIn&);
      ~PsiUHF();
      void print(std::ostream& os = ExEnv::out0()) const;
 
      void write_basic_input(int conv);
      bool analytic_gradient_implemented() const {
        return true;
      }
      PsiSCF::RefType reftype() const {
        return uhf;
      }
  };
 
  namespace detail {
    RefSymmSCMatrix rdopdm(SpinCase1 spin,
                           const std::vector<unsigned int>& mopi,
                           const std::vector<unsigned int>& dmap,
                           Ref<SCMatrixKit> kit = 0);
    RefSymmSCMatrix rdtpdm(SpinCase2 pairspin,
                           const std::vector<unsigned int>& dmap,
                           bool spinfree = false,
                           Ref<SCMatrixKit> kit = 0);
  }
}
#endif