r12technology.h

//
// r12technology.h
//
// Copyright (C) 2007 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_r12technology_h
#define _chemistry_qc_mbptr12_r12technology_h
 
#include <util/state/state.h>
#include <chemistry/qc/basis/intdescr.h>
#include <math/optimize/gaussianfit.h>
 
namespace sc {
 
// //////////////////////////////////////////////////////////////////////////
 
class R12Technology: virtual public SavableState {
  public:
 
  enum Projector {Projector_0 = 0,
    Projector_1 = 1,
    Projector_2 = 2,
    Projector_3 = 3};
  enum StandardApproximation {
    //StdApprox_A = 0, // is now obsolete
    StdApprox_Ap = 1,
    StdApprox_App = 2,
    StdApprox_B = 3,
    StdApprox_C = 4,
    StdApprox_Cp = 5
    };
  enum ABSMethod {ABS_CABS = 2,
    ABS_CABSPlus = 3};
 
  enum OrbitalProduct_GG {
    OrbProdGG_ij = 0,
    OrbProdGG_pq = 1,
  };
 
  enum OrbitalProduct_gg {
    OrbProdgg_ij = 0,
    OrbProdgg_pq = 1,
  };
 
  enum PositiveDefiniteB {
    PositiveDefiniteB_no = 0,
    PositiveDefiniteB_yes = 1,
    PositiveDefiniteB_weak = 2
  };
 
  enum GeminalAmplitudeAnsatz {
    GeminalAmplitudeAnsatz_fullopt = 0,
    GeminalAmplitudeAnsatz_fixed = 1,
    GeminalAmplitudeAnsatz_scaledfixed = 2
  };
 
  enum H0_dk_approx_pauli {
    H0_dk_approx_pauli_true = 0,
    H0_dk_approx_pauli_fHf = 1,
    H0_dk_approx_pauli_fHf_Q = 2,
    H0_dk_approx_pauli_false = 3
  };
 
  class R12Ansatz : virtual public SavableState {
    public:
    R12Ansatz(const Ref<KeyVal>&);
    R12Ansatz(StateIn&);
    R12Ansatz();
    ~R12Ansatz();
 
    void save_data_state(StateOut&);
    void print(std::ostream& o =ExEnv::out0()) const;
 
    R12Technology::Projector projector() const;
    bool diag() const;
    R12Technology::GeminalAmplitudeAnsatz amplitudes() const;
    bool wof() const;
    R12Technology::OrbitalProduct_GG orbital_product_GG() const;
    R12Technology::OrbitalProduct_gg orbital_product_gg() const;
 
    private:
    R12Technology::Projector projector_;
    bool diag_;
    R12Technology::GeminalAmplitudeAnsatz amplitudes_;
    bool scaled_;     //<
    bool wof_;
    R12Technology::OrbitalProduct_GG orbital_product_GG_;
    R12Technology::OrbitalProduct_gg orbital_product_gg_;
  }; // end of R12Ansatz declaration
 
 
  class GeminalDescriptor : public RefCount {
    private:
      std::string type_;
      std::vector<std::string> params_;
    public:
      GeminalDescriptor();
      ~GeminalDescriptor(){}
      GeminalDescriptor(const std::string& type, const std::vector<std::string> &params);
      GeminalDescriptor(const GeminalDescriptor& source);
      std::string type() const;
      std::vector<std::string> params() const;
      void print(std::ostream &o=ExEnv::out0());
  };
 
  static bool invalid(const Ref<GeminalDescriptor>& gdesc);
  static bool R12(const Ref<GeminalDescriptor>& gdesc);
  static bool STG(const Ref<GeminalDescriptor>& gdesc);
  static bool G12(const Ref<GeminalDescriptor>& gdesc);
  static double single_slater_exponent(const Ref<GeminalDescriptor>& gdesc);
 
  class CorrelationFactor : public RefCount {
    public:
      //typedef IntParamsG12::PrimitiveGeminal PrimitiveGeminal;
      //typedef IntParamsG12::ContractedGeminal ContractedGeminal;
      //typedef std::vector<ContractedGeminal> CorrelationParameters;
 
      CorrelationFactor(const std::string& label, const Ref<GeminalDescriptor> &geminaldescriptor);
      CorrelationFactor();
      virtual ~CorrelationFactor();
 
      // return true if this is equivalent to cf
      virtual bool equiv(const Ref<CorrelationFactor>& cf) const =0;
 
      const std::string& label() const;
      virtual unsigned int nfunctions() const;
      virtual unsigned int nprimitives(unsigned int c) const;
 
      virtual double value(unsigned int c, double r12) const =0;
 
      virtual Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f) const;
      virtual Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int fbra, unsigned int fket) const;
      virtual unsigned int max_num_tbint_types() const =0;
 
      //
      // These functions are used to map the logical type of integrals ([T1,F12], etc.) to concrete types as produced by TwoBodyInt
      //
 
      virtual TwoBodyOper::type tbint_type_eri() const;
      virtual TwoBodyOper::type tbint_type_f12() const;
      virtual TwoBodyOper::type tbint_type_t1f12() const;
      virtual TwoBodyOper::type tbint_type_t2f12() const;
      virtual TwoBodyOper::type tbint_type_f12eri() const;
      virtual TwoBodyOper::type tbint_type_f12f12() const;
      virtual TwoBodyOper::type tbint_type_f12t1f12() const;
      virtual TwoBodyOper::type tbint_type_f12f12_anti() const;
 
      void print(std::ostream& os = ExEnv::out0()) const;
      Ref<GeminalDescriptor> geminaldescriptor();
 
    protected:
      std::string label_;
      Ref<GeminalDescriptor> geminaldescriptor_;
 
      virtual void print_params(std::ostream& os, unsigned int f) const;
 
  };
 
  class NullCorrelationFactor : public CorrelationFactor {
    public:
    NullCorrelationFactor();
 
    bool equiv(const Ref<CorrelationFactor>& cf) const;
    unsigned int max_num_tbint_types() const { return 1; }
    double value(unsigned int c, double r12) const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f) const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int fbra, unsigned int fket) const;
  };
 
  class R12CorrelationFactor : public CorrelationFactor {
    public:
    R12CorrelationFactor();
 
    bool equiv(const Ref<CorrelationFactor>& cf) const;
    unsigned int max_num_tbint_types() const { return 4; }
    TwoBodyOper::type tbint_type_f12() const;
    TwoBodyOper::type tbint_type_t1f12() const;
    TwoBodyOper::type tbint_type_t2f12() const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f) const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f, unsigned int g) const;
    double value(unsigned int c, double r12) const;
  };
 
  template<class IntParam>
  struct CorrParamCompare {
      typedef typename IntParam::PrimitiveGeminal PrimitiveGeminal;
      typedef typename IntParam::ContractedGeminal ContractedGeminal;
      typedef std::vector<ContractedGeminal> ContractedGeminals;
 
      // 2 parameters are equivalent if their values differ by less than epsilon
      static double epsilon;
      static bool equiv(const ContractedGeminals& A,
                        const ContractedGeminals& B);
 
    private:
      static bool equiv(const PrimitiveGeminal& A,
                        const PrimitiveGeminal& B);
  };
 
  class G12CorrelationFactor : public CorrelationFactor {
    public:
    typedef IntParamsG12::PrimitiveGeminal PrimitiveGeminal;
    typedef IntParamsG12::ContractedGeminal ContractedGeminal;
    typedef std::vector<ContractedGeminal> CorrelationParameters;
 
    G12CorrelationFactor(const CorrelationParameters& params, const Ref<GeminalDescriptor> &geminaldescriptor = 0);
 
    bool equiv(const Ref<CorrelationFactor>& cf) const;
    unsigned int nfunctions() const;
    const ContractedGeminal& function(unsigned int c) const;
    unsigned int nprimitives(unsigned int c) const;
    const PrimitiveGeminal& primitive(unsigned int c, unsigned int p) const;
    unsigned int max_num_tbint_types() const { return 6; }
    TwoBodyOper::type tbint_type_f12() const;
    TwoBodyOper::type tbint_type_f12eri() const;
    TwoBodyOper::type tbint_type_t1f12() const;
    TwoBodyOper::type tbint_type_t2f12() const;
    TwoBodyOper::type tbint_type_f12f12() const;
    TwoBodyOper::type tbint_type_f12t1f12() const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f) const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int fbra, unsigned int fket) const;
    double value(unsigned int c, double r12) const;
 
  private:
    CorrelationParameters params_;
 
  void print_params(std::ostream& os, unsigned int f) const;
 
  };
 
  class G12NCCorrelationFactor : public CorrelationFactor {
    public:
    typedef IntParamsG12::PrimitiveGeminal PrimitiveGeminal;
    typedef IntParamsG12::ContractedGeminal ContractedGeminal;
    typedef std::vector<ContractedGeminal> CorrelationParameters;
 
    G12NCCorrelationFactor(const CorrelationParameters& params, const Ref<GeminalDescriptor> &geminaldescriptor = 0);
 
    bool equiv(const Ref<CorrelationFactor>& cf) const;
    unsigned int nfunctions() const;
    const ContractedGeminal& function(unsigned int c) const;
    unsigned int nprimitives(unsigned int c) const;
    const PrimitiveGeminal& primitive(unsigned int c, unsigned int p) const;
    unsigned int max_num_tbint_types() const { return 6; }
    TwoBodyOper::type tbint_type_f12() const;
    TwoBodyOper::type tbint_type_f12eri() const;
    TwoBodyOper::type tbint_type_f12f12() const;
    TwoBodyOper::type tbint_type_f12t1f12() const;
    TwoBodyOper::type tbint_type_f12f12_anti() const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int f) const;
    Ref<TwoBodyIntDescr> tbintdescr(const Ref<Integral>& IF, unsigned int fbra, unsigned int fket) const;
    double value(unsigned int c, double r12) const;
 
    static ContractedGeminal product(const ContractedGeminal& A,
                                     const ContractedGeminal& B);
 
  private:
    CorrelationParameters params_;
 
    void print_params(std::ostream& os, unsigned int f) const;
 
  };
 
  class GeminalDescriptorFactory : public RefCount {
    private:
      const char* invalid_id_;  // = "invalid"
      const char* r12_id_;  // = "R12"
      const char* stg_id_;  // = "STG"
      const char* g12_id_;  // = "G12"
    public:
      GeminalDescriptorFactory();
      Ref<GeminalDescriptor> null_geminal();
      Ref<GeminalDescriptor> r12_geminal();
      Ref<GeminalDescriptor> slater_geminal(double gamma);
      Ref<GeminalDescriptor> slater_geminal(const std::vector<double> &gamma);
      Ref<GeminalDescriptor> gaussian_geminal(double gamma);
      Ref<GeminalDescriptor> contracted_gaussian_geminal(const std::vector<double> &coeff,
                                                         const std::vector<double> &gamma);
      Ref<GeminalDescriptor> gaussian_geminal(const G12CorrelationFactor::CorrelationParameters &corrparams);
  };
 
  template<class CF>
      static Ref<CF> direct_product(const Ref<CF>& A, const Ref<CF>& B) {
        const unsigned int nf_A = A->nfunctions();
        const unsigned int nf_B = B->nfunctions();
        typedef typename CF::CorrelationParameters CorrParams;
        CorrParams corrparams;
        for (int f = 0; f < nf_A; ++f) {
          for (int g = 0; g < nf_B; ++g) {
            corrparams.push_back(CF::product(A->function(f), B->function(g)));
          }
        }
        return new CF(corrparams);
      }
 
  private:
 
    // implements KeyVal constructor
    void init_from_kv(const Ref<KeyVal>& keyval,
                      bool abs_eq_obs = true,
                      bool vbs_eq_obs = true);
 
    bool abs_eq_obs_;
    bool vbs_eq_obs_;
 
    Ref<CorrelationFactor> corrfactor_;
    StandardApproximation stdapprox_;
    Ref<R12Ansatz> ansatz_;
    ABSMethod abs_method_;
    double abs_lindep_tol_;
    int abs_nlindep_;
    unsigned int maxnabs_;
    bool gbc_;
    bool ebc_;
    bool coupling_;
    bool compute_1rdm_;  // compute mp2r12 1e density
    bool coupling_1rdm_f12b_; // compute coupling contri. to ccsd-f12b 1e density
    bool omit_P_;
    H0_dk_approx_pauli H0_dk_approx_pauli_;
    bool H0_dk_keep_;
    bool safety_check_;
    PositiveDefiniteB posdef_B_;
 
    // for debugging purposes only
    bool omit_B_;
 
    // no need to store this guy
#if 0
    // determines the weight function used to fit the correlation factor
    struct GTGFitWeight {
      typedef enum {TewKlopper, Cusp} Type;
    };
    GTGFitWeight::Type gtg_fit_weight_;
#endif
 
  public:
    R12Technology(StateIn&);
    R12Technology(const Ref<KeyVal>&);
    R12Technology(const Ref<KeyVal>&,
          const Ref<GaussianBasisSet>& bs,
          const Ref<GaussianBasisSet>& vbs,
          const Ref<GaussianBasisSet>& abs
    );
    ~R12Technology();
 
    void save_data_state(StateOut&);
 
    const Ref<CorrelationFactor>& corrfactor() const;
    void corrfactor(const Ref<CorrelationFactor>&);
    unsigned int maxnabs() const;
    bool gbc() const;
    bool ebc() const;
    bool coupling() const;
    bool compute_1rdm() const;
    bool coupling_1rdm_f12b() const;
    ABSMethod abs_method() const;
    int abs_nlindep() const;
    double abs_lindep_tol() const;
    StandardApproximation stdapprox() const;
    const Ref<R12Ansatz>& ansatz() const;
    bool spinadapted() const;
    bool omit_P() const;
    H0_dk_approx_pauli H0_dk_approx() const;
    bool H0_dk_keep() const;
    bool safety_check() const;
    PositiveDefiniteB posdef_B() const;
 
    //
    // these are for debugging only
    //
    // omit expensive parts of B
    bool omit_B() const;
 
    // This checks if ints is suitable for R12 calculations. Throws, if false.
    void check_integral_factory(const Ref<Integral>& ints);
 
    void print(std::ostream&o=ExEnv::out0()) const;
 
    static Ref<GaussianBasisSet> make_auto_cabs(const Ref<GaussianBasisSet>& bs);
    static std::string default_cabs_name(const std::string& obs_name);
    static double default_stg_exponent(const std::string& obs_name);
};
 
template <class IntParam> double R12Technology::CorrParamCompare<IntParam>::epsilon(1e-6);
 
template<class IntParam>
  bool R12Technology::CorrParamCompare<IntParam>::equiv(
                                                        const ContractedGeminals& A,
                                                        const ContractedGeminals& B) {
    unsigned int nf = A.size();
    if (nf != B.size())
      return false;
 
    for (unsigned int f = 0; f < nf; ++f) {
      const ContractedGeminal& Af = A[f];
      const ContractedGeminal& Bf = B[f];
      unsigned int np = Af.size();
      if (np != Bf.size())
        return false;
      for (unsigned int p = 0; p < np; ++p) {
        if (!equiv(Af[p], Bf[p]))
          return false;
      }
    }
 
    return true;
  }
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: