r12int_eval.h

//
// r12int_eval.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_r12inteval_h
#define _chemistry_qc_mbptr12_r12inteval_h
 
#include <util/ref/ref.h>
#include <chemistry/qc/mbptr12/r12wfnworld.h>
#include <chemistry/qc/mbptr12/r12technology.h>
#include <chemistry/qc/mbptr12/twoparticlecontraction.h>
#include <chemistry/qc/wfn/spin.h>
#include <chemistry/qc/mbptr12/fixedcoefficient.h>
#include <chemistry/qc/mbptr12/twobodytensorinfo.h>
 
namespace sc {
 
  class TwoBodyMOIntsTransform;
  class DistArray4;
  class R12Amplitudes;
  class R12EnergyIntermediates;
 
class R12IntEval : virtual public SavableState {
  private:
  // change to false to use the old fock builder
  static const bool USE_FOCKBUILD = true;
 
 
  bool evaluated_;
 
  // Calculation information (number of basis functions, R12 approximation, etc.)
  Ref<R12WavefunctionWorld> r12world_;
 
  RefSCMatrix V_[NSpinCases2];
  RefSCMatrix X_[NSpinCases2];
  // Note that intermediate B is symmetric but is stored as a full matrix to simplify the code
  // that computes asymmetric form of B
  RefSCMatrix B_[NSpinCases2];
  RefSCMatrix BB_[NSpinCases2];  // The difference between B intermediate of approximation B and A'
  RefSCMatrix A_[NSpinCases2];
  Ref<CuspConsistentGeminalCoefficient> cuspconsistentgeminalcoefficient_[NSpinCases2];
 
  RefSCVector emp2pair_[NSpinCases2];
  RefSCDimension dim_oo_[NSpinCases2];
  RefSCDimension dim_vv_[NSpinCases2];
  RefSCDimension dim_aa_[NSpinCases2];
  RefSCDimension dim_f12_[NSpinCases2];
  RefSCDimension dim_GG_[NSpinCases2];
  RefSCDimension dim_gg_[NSpinCases2];
 
  Ref<R12Amplitudes> Amps_;  // First-order amplitudes of various contributions to the pair functions
  double emp2_obs_singles_;
  double emp2_cabs_singles_;
  RefSCMatrix T1_cabs_[NSpinCases1];
  int debug_;
 
  mutable RefSymmSCMatrix ordm_[NSpinCases1];  
 
  void spinadapt_mospace_labels(SpinCase1 spin, std::string& id, std::string& name) const;
 
  void f_bra_ket(SpinCase1 spin,
         bool make_F,
         bool make_hJ,
         bool make_K,
         Ref<OrbitalSpace>& F,
         Ref<OrbitalSpace>& hJ,
         Ref<OrbitalSpace>& K,
         const Ref<OrbitalSpace>& extspace,
         const Ref<OrbitalSpace>& intspace
      );
  Ref<OrbitalSpace> hj_i_p_[NSpinCases1];
  Ref<OrbitalSpace> hj_i_A_[NSpinCases1];
  Ref<OrbitalSpace> hj_i_P_[NSpinCases1];
  Ref<OrbitalSpace> hj_i_m_[NSpinCases1];
  Ref<OrbitalSpace> hj_i_a_[NSpinCases1];
  Ref<OrbitalSpace> hj_m_m_[NSpinCases1];
  Ref<OrbitalSpace> hj_m_p_[NSpinCases1];
  Ref<OrbitalSpace> hj_a_A_[NSpinCases1];
  Ref<OrbitalSpace> hj_p_p_[NSpinCases1];
  Ref<OrbitalSpace> hj_p_A_[NSpinCases1];
  Ref<OrbitalSpace> hj_p_P_[NSpinCases1];
  Ref<OrbitalSpace> hj_P_P_[NSpinCases1];
  Ref<OrbitalSpace> hj_p_m_[NSpinCases1];
  Ref<OrbitalSpace> hj_p_a_[NSpinCases1];
  Ref<OrbitalSpace> K_i_p_[NSpinCases1];
  Ref<OrbitalSpace> K_i_m_[NSpinCases1];
  Ref<OrbitalSpace> K_i_a_[NSpinCases1];
  Ref<OrbitalSpace> K_i_A_[NSpinCases1];
  Ref<OrbitalSpace> K_i_P_[NSpinCases1];
  Ref<OrbitalSpace> K_m_a_[NSpinCases1];
  Ref<OrbitalSpace> K_a_a_[NSpinCases1];
  Ref<OrbitalSpace> K_a_p_[NSpinCases1];
  Ref<OrbitalSpace> K_a_P_[NSpinCases1];
  Ref<OrbitalSpace> K_p_p_[NSpinCases1];
  Ref<OrbitalSpace> K_p_m_[NSpinCases1];
  Ref<OrbitalSpace> K_p_a_[NSpinCases1];
  Ref<OrbitalSpace> K_p_A_[NSpinCases1];
  Ref<OrbitalSpace> K_p_P_[NSpinCases1];
  Ref<OrbitalSpace> K_A_P_[NSpinCases1];
  Ref<OrbitalSpace> K_P_P_[NSpinCases1];
  Ref<OrbitalSpace> F_P_P_[NSpinCases1];
  Ref<OrbitalSpace> F_p_A_[NSpinCases1];
  Ref<OrbitalSpace> F_p_p_[NSpinCases1];
  Ref<OrbitalSpace> F_p_m_[NSpinCases1];
  Ref<OrbitalSpace> F_p_a_[NSpinCases1];
  Ref<OrbitalSpace> F_m_m_[NSpinCases1];
  Ref<OrbitalSpace> F_m_a_[NSpinCases1];
  Ref<OrbitalSpace> F_m_p_[NSpinCases1];
  Ref<OrbitalSpace> F_m_P_[NSpinCases1];
  Ref<OrbitalSpace> F_gg_P_[NSpinCases1];
  Ref<OrbitalSpace> F_m_A_[NSpinCases1];
  Ref<OrbitalSpace> F_i_A_[NSpinCases1];
  Ref<OrbitalSpace> F_i_m_[NSpinCases1];
  Ref<OrbitalSpace> F_i_a_[NSpinCases1];
  Ref<OrbitalSpace> F_i_p_[NSpinCases1];
  Ref<OrbitalSpace> F_i_P_[NSpinCases1];
  Ref<OrbitalSpace> F_a_a_[NSpinCases1];
  Ref<OrbitalSpace> F_a_A_[NSpinCases1];
  Ref<OrbitalSpace> h_P_P_[NSpinCases1];
  Ref<OrbitalSpace> J_i_p_[NSpinCases1];
  Ref<OrbitalSpace> J_i_P_[NSpinCases1];
  Ref<OrbitalSpace> J_P_P_[NSpinCases1];
  Ref<OrbitalSpace> F_A_A_[NSpinCases1];
  Ref<OrbitalSpace> F_p_P_[NSpinCases1];
  Ref<OrbitalSpace> gamma_p_p_[NSpinCases1];
  Ref<OrbitalSpace> gamma_m_m_[NSpinCases1];
  Ref<OrbitalSpace> gammaFgamma_p_p_[NSpinCases1];
  Ref<OrbitalSpace> Fgamma_p_P_[NSpinCases1];
 
  Ref<OrbitalSpace> cabs_space_fockcanonical(SpinCase1 s,
                                             double scale_H,
                                             double scale_J,
                                             double scale_K);
 
  void init_intermeds_();
  void init_intermeds_r12_();
  void init_intermeds_g12_();
  void r2_contrib_to_X_new_();
  RefSCMatrix compute_I_(const Ref<OrbitalSpace>& space1,
             const Ref<OrbitalSpace>& space2,
             const Ref<OrbitalSpace>& space3,
             const Ref<OrbitalSpace>& space4);
  RefSCMatrix compute_r2_(const Ref<OrbitalSpace>& space1,
                          const Ref<OrbitalSpace>& space2,
                          const Ref<OrbitalSpace>& space3,
                          const Ref<OrbitalSpace>& space4);
 
#if 0
 
  RefSCMatrix Delta_DKH_(const Ref<OrbitalSpace>& bra_space,
                         const Ref<OrbitalSpace>& ket_space,
                         SpinCase1 S = Alpha);
  // Computes T, V and the mass-velocity term in the momentum basis.
  // It's a modified version of Wavefunction::core_hamiltonian_dk
  RefSymmSCMatrix hcore_plus_massvelocity_(const Ref<GaussianBasisSet> &bas,
                                           const Ref<GaussianBasisSet> &p_bas,
                                           bool include_T = true,
                                           bool include_V = true,
                                           bool include_MV = true);
#endif
 
  // Computes the skalar Pauli-Hamiltonian (T + V + mass_velocity + Darwin),
  // with the mass-velocity term evaluated in the momentum basis.
  RefSymmSCMatrix pauli(const Ref<GaussianBasisSet> &bas,
                        const Ref<GaussianBasisSet> &pbas = 0,
                        const bool momentum=false);
  RefSymmSCMatrix pauli_realspace(const Ref<GaussianBasisSet> &bas);
  RefSymmSCMatrix pauli_momentumspace(const Ref<GaussianBasisSet> &bas,
                        const Ref<GaussianBasisSet> &p_bas);
 
  RefSCMatrix coulomb_(const SpinCase1 &spin, const Ref<OrbitalSpace>& bra_space,
                       const Ref<OrbitalSpace>& ket_space);
  RefSCMatrix coulomb_(const Ref<OrbitalSpace>& occ_space, const Ref<OrbitalSpace>& bra_space,
                       const Ref<OrbitalSpace>& ket_space);
  RefSCMatrix exchange_(const SpinCase1 &spin, const Ref<OrbitalSpace>& bra_space,
                        const Ref<OrbitalSpace>& ket_space);
  RefSCMatrix exchange_(const Ref<OrbitalSpace>& occ_space, const Ref<OrbitalSpace>& bra_space,
                        const Ref<OrbitalSpace>& ket_space);
 
  void checkpoint_() const;
 
  void contrib_to_VXB_a_();
  void contrib_to_VXB_c_ansatz1_();
  void contrib_to_VX_GenRefansatz2_();
  void contrib_to_VX_GenRefansatz2_spinfree_();
  void contrib_to_VXB_a_vbsneqobs_();
  void contrib_to_B_DKH_a_();
 
  void compute_mp2_pair_energies_(RefSCVector& emp2pair,
                                  SpinCase2 S,
                                  const Ref<OrbitalSpace>& space1,
                                  const Ref<OrbitalSpace>& space2,
                                  const Ref<OrbitalSpace>& space3,
                                  const Ref<OrbitalSpace>& space4,
                                  const std::string& tform_key);
 
  void compute_A(); //< computes A ("coupling") matrix for all spincases
 
  void compute_A_direct_(RefSCMatrix& A,
                         const Ref<OrbitalSpace>& space1,
                         const Ref<OrbitalSpace>& space2,
                         const Ref<OrbitalSpace>& space3,
                         const Ref<OrbitalSpace>& space4,
                         const Ref<OrbitalSpace>& rispace2,
                         const Ref<OrbitalSpace>& rispace4,
                         bool antisymmetrize);
 
  // make compute_tbint_tensor_ public for now
  public:
  template <typename DataProcess, bool CorrFactorInBra, bool CorrFactorInKet>
  DEPRECATED void compute_tbint_tensor(RefSCMatrix& T,
                              TwoBodyOper::type tbint_type,
                              const Ref<OrbitalSpace>& space1,
                              const Ref<OrbitalSpace>& space2,
                              const Ref<OrbitalSpace>& space3,
                              const Ref<OrbitalSpace>& space4,
                              bool antisymmetrize,
                              const std::vector<std::string>& tform_keys);
 
  private:
  template <typename DataProcessBra,
            typename DataProcessKet,
            typename DataProcessBraKet,
            bool CorrFactorInBra,
            bool CorrFactorInKet,
            bool CorrFactorInInt>
  DEPRECATED void contract_tbint_tensor(
           RefSCMatrix& T,
           TwoBodyOper::type tbint_type_bra,
           TwoBodyOper::type tbint_type_ket,
           const Ref<OrbitalSpace>& space1_bra,
           const Ref<OrbitalSpace>& space2_bra,
           const Ref<OrbitalSpace>& space1_intb,
           const Ref<OrbitalSpace>& space2_intb,
           const Ref<OrbitalSpace>& space1_ket,
           const Ref<OrbitalSpace>& space2_ket,
           const Ref<OrbitalSpace>& space1_intk,
           const Ref<OrbitalSpace>& space2_intk,
           const Ref<mbptr12::TwoParticleContraction>& tpcontract,
           bool antisymmetrize,
           const std::vector<std::string>& tformkeys_bra,
           const std::vector<std::string>& tformkeys_ket);
 
  template <bool CorrFactorInBra,
            bool CorrFactorInKet>
    void DEPRECATED contract_tbint_tensor(
           RefSCMatrix& T,
           TwoBodyOper::type tbint_type_bra,
           TwoBodyOper::type tbint_type_ket,
           double scale,
           const Ref<OrbitalSpace>& space1_bra,
           const Ref<OrbitalSpace>& space2_bra,
           const Ref<OrbitalSpace>& space1_intb,
           const Ref<OrbitalSpace>& space2_intb,
           const Ref<OrbitalSpace>& space1_ket,
           const Ref<OrbitalSpace>& space2_ket,
           const Ref<OrbitalSpace>& space1_intk,
           const Ref<OrbitalSpace>& space2_intk,
           bool antisymmetrize,
           const std::vector<std::string>& tformkeys_bra,
           const std::vector<std::string>& tformkeys_ket);
 
  template <bool CorrFactorInBra,
            bool CorrFactorInKet>
  void
  contract_tbint_tensor(
           std::vector< Ref<DistArray4> >& results,
           TwoBodyOper::type tbint_type_bra,
           TwoBodyOper::type tbint_type_ket,
           double scale,
           const Ref<OrbitalSpace>& space1_bra,
           const Ref<OrbitalSpace>& space2_bra,
           const Ref<OrbitalSpace>& space1_intb,
           const Ref<OrbitalSpace>& space2_intb,
           const Ref<OrbitalSpace>& space1_ket,
           const Ref<OrbitalSpace>& space2_ket,
           const Ref<OrbitalSpace>& space1_intk,
           const Ref<OrbitalSpace>& space2_intk,
           bool antisymmetrize,
           const std::vector<std::string>& tformkeys_bra,
           const std::vector<std::string>& tformkeys_ket);
 
  // for now make it public
  public:
  template <typename DataProcessBra,
            typename DataProcessKet,
            bool CorrFactorInBra,
            bool CorrFactorInKet,
            bool CorrFactorInInt>
    void contract_tbint_tensors_to_obtensor(RefSCMatrix& T,
                                            SpinCase2 pairspin,
                                            TwoBodyTensorInfo tbtensor_type_bra,
                                            TwoBodyTensorInfo tbtensor_type_ket,
                                            const Ref<OrbitalSpace>& space1_bra,
                                            const Ref<OrbitalSpace>& space1_intb,
                                            const Ref<OrbitalSpace>& space2_intb,
                                            const Ref<OrbitalSpace>& space3_intb,
                                            const Ref<OrbitalSpace>& space1_ket,
                                            const Ref<OrbitalSpace>& space1_intk,
                                            const Ref<OrbitalSpace>& space2_intk,
                                            const Ref<OrbitalSpace>& space3_intk,
                                            const Ref<mbptr12::TwoParticleContraction>& tpcontract,
                                            const std::vector<std::string>& tformkeys_bra,
                                            const std::vector<std::string>& tformkeys_ket);
 
  private:
  void compute_X_(RefSCMatrix& X,
                  SpinCase2 sc2,
                  const Ref<OrbitalSpace>& bra1,
                  const Ref<OrbitalSpace>& bra2,
                  const Ref<OrbitalSpace>& ket1,
                  const Ref<OrbitalSpace>& ket2,
                  bool F2_only = false);
  void compute_FxF_(RefSCMatrix& FxF,
                    SpinCase2 sc2,
                    const Ref<OrbitalSpace>& bra1,
                    const Ref<OrbitalSpace>& bra2,
                    const Ref<OrbitalSpace>& ket1,
                    const Ref<OrbitalSpace>& ket2,
                    const Ref<OrbitalSpace>& int1,
                    const Ref<OrbitalSpace>& int2,
                    const Ref<OrbitalSpace>& intk1,
                    const Ref<OrbitalSpace>& intk2,
                    const Ref<OrbitalSpace>& intkx1,
                    const Ref<OrbitalSpace>& intkx2
                   );
 
  void AT2_contrib_to_V_();
 
  void AF12_contrib_to_B_();
 
  void compute_B_gbc_();
 
  void compute_B_fX_();
 
  void compute_BB_();
 
  void compute_BC_();
  void compute_BC_ansatz1_();
  void compute_BC_GenRefansatz2_();
  void compute_BC_GenRefansatz2_spinfree();
 
 
  void compute_BCp_();
 
  void compute_BApp_();
 
  void compute_B_DKH_();
 
  const RefSCVector& compute_emp2(SpinCase2 s);
  double compute_emp2_obs_singles(bool obs_singles);
#if 0
 
  double compute_emp2_cabs_singles();
#endif
 
  double compute_emp2_cabs_singles_noncanonical(bool vir_cabs_coupling);
  double compute_emp2_cabs_singles_noncanonical_ccsd(const RefSCMatrix& T1_ia_alpha,
                                                     const RefSCMatrix& T1_ia_beta);
 
  void compute_R_vbsneqobs_(const Ref<TwoBodyMOIntsTransform>& ipjq_tform,
                            RefSCMatrix& Raa, RefSCMatrix& Rab);
 
  void compute_amps_();
 
  void globally_sum_scmatrix_(RefSCMatrix& A, bool to_all_tasks = false, bool to_average = false);
 
  void globally_sum_scvector_(RefSCVector& A, bool to_all_tasks = false, bool to_average = false);
 
  void globally_sum_intermeds_(bool to_all_tasks = false);
 
public:
  R12IntEval(StateIn&);
  R12IntEval(const Ref<R12WavefunctionWorld>& r12w);
  ~R12IntEval();
 
  void save_data_state(StateOut&);
  virtual void obsolete();
 
  void debug(int d) { debug_ = d; }
 
  int dk() const { return r12world()->refwfn()->dk(); }
 
#if 1
  const Ref<R12Technology::CorrelationFactor>& corrfactor() const { return r12world()->r12tech()->corrfactor(); }
  R12Technology::ABSMethod abs_method() const { return r12world()->r12tech()->abs_method(); }
  const Ref<R12Technology::R12Ansatz>& ansatz() const { return r12world()->r12tech()->ansatz(); }
  bool spin_polarized() const { return r12world()->refwfn()->spin_polarized(); }
  bool gbc() const { return r12world()->r12tech()->gbc(); }
  bool ebc() const { return r12world()->r12tech()->ebc(); }
  bool coupling() const { return r12world()->r12tech()->coupling(); }
  bool compute_1rdm() const { return r12world()->r12tech()->compute_1rdm(); }
  bool coupling_1rdm_f12b() const { return r12world()->r12tech()->coupling_1rdm_f12b(); }
  R12Technology::StandardApproximation stdapprox() const { return r12world()->r12tech()->stdapprox(); }
  bool omit_P() const { return r12world()->r12tech()->omit_P(); }
  const Ref<MOIntsTransformFactory>& tfactory() const { return r12world()->world()->tfactory(); };
  const Ref<MOIntsRuntime>& moints_runtime() const { return r12world()->world()->moints_runtime(); };
  const Ref<TwoBodyFourCenterMOIntsRuntime>& moints_runtime4() const { return r12world()->world()->moints_runtime()->runtime_4c(); };
  const Ref<FockBuildRuntime>& fockbuild_runtime() const { return r12world()->world()->fockbuild_runtime(); };
 
#endif
 
  const Ref<R12WavefunctionWorld>& r12world() const { return r12world_; }
 
  RefSCDimension dim_oo(SpinCase2 S) const { return dim_oo_[S]; }
  RefSCDimension dim_vv(SpinCase2 S) const { return dim_vv_[S]; }
  RefSCDimension dim_aa(SpinCase2 S) const { return dim_aa_[S]; }
  RefSCDimension dim_f12(SpinCase2 S) const { return dim_f12_[S]; }
  RefSCDimension dim_GG(SpinCase2 S) const { return dim_GG_[S]; }
  RefSCDimension dim_gg(SpinCase2 S) const { return dim_gg_[S]; }
 
  int nspincases1() const {
    if (r12world()->spinadapted()) return 1;
    return ::sc::nspincases1(r12world()->refwfn()->spin_polarized());
  }
  int nspincases2() const {
    if (r12world()->spinadapted()) return 1;
    return ::sc::nspincases2(r12world()->refwfn()->spin_polarized());
  }
 
  virtual void compute();
 
  bool bc() const;
 
  Ref<R12Amplitudes> amps();
 
  const RefSCMatrix& V(SpinCase2 S);
  const RefSCMatrix& V();
  RefSymmSCMatrix X(SpinCase2 S);
  RefSymmSCMatrix X();
  RefSymmSCMatrix B(SpinCase2 S);
  RefSymmSCMatrix B();
  RefSymmSCMatrix BB(SpinCase2 S);
  const RefSCMatrix& A(SpinCase2 S);
  const RefSCMatrix& T2(SpinCase2 S);
  const RefSCMatrix& F12(SpinCase2 S);
 
  RefSCMatrix V(SpinCase2 spincase2,
                const Ref<OrbitalSpace>& p,
                const Ref<OrbitalSpace>& q);
  std::vector<Ref<DistArray4> > V_distarray4(SpinCase2 spincase2,
                               const Ref<OrbitalSpace>& p,
                               const Ref<OrbitalSpace>& q);
  std::vector<Ref<DistArray4> > U_distarray4(SpinCase2 spincase2,
                                             const Ref<OrbitalSpace>& p,
                                             const Ref<OrbitalSpace>& q);
  RefSymmSCMatrix P(SpinCase2 S);
  double C_CuspConsistent(int i,int j,int k,int l,SpinCase2 pairspin);
 
  double emp2_obs_singles();
  double emp2_cabs_singles(bool vir_cabs_coupling = true);
  double emp2_cabs_singles(const RefSCMatrix& T1_ia_alpha,
                           const RefSCMatrix& T1_ia_beta);
  const RefSCVector& emp2(SpinCase2 S);
 
  const RefSCMatrix& T1_cabs(SpinCase1 spin) const;
 
  const Ref<OrbitalSpace>& occ_act(SpinCase1 S) const;
  const Ref<OrbitalSpace>& occ(SpinCase1 S) const;
  const Ref<OrbitalSpace>& vir_act(SpinCase1 S) const;
  const Ref<OrbitalSpace>& vir(SpinCase1 S) const;
  const Ref<OrbitalSpace>& orbs(SpinCase1 S) const;
  const Ref<OrbitalSpace>& GGspace(SpinCase1 S) const;
  const Ref<OrbitalSpace>& ggspace(SpinCase1 S) const;
  const Ref<OrbitalSpace>& cabs_space_canonical(SpinCase1 s);
  const Ref<OrbitalSpace>& cabs_space_hcanonical(SpinCase1 s);
  RefSymmSCMatrix ordm(SpinCase1 S) const;
  RefSymmSCMatrix ordm() const;
  RefSymmSCMatrix ordm_av() const;
  RefSymmSCMatrix ordm_occ_av() const;
  const Ref<OrbitalSpace>& hj_x_P(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_x_A(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_x_p(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_x_m(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_x_a(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_i_P(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_i_A(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_i_p(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_i_m(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_i_a(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_m_m(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_m_p(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_a_A(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_p_P(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_p_A(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_p_p(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_p_m(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_p_a(SpinCase1 S);
  const Ref<OrbitalSpace>& hj_P_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_x_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_x_A(SpinCase1 S);
  const Ref<OrbitalSpace>& K_x_p(SpinCase1 S);
  const Ref<OrbitalSpace>& K_x_m(SpinCase1 S);
  const Ref<OrbitalSpace>& K_x_a(SpinCase1 S);
  const Ref<OrbitalSpace>& K_i_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_i_A(SpinCase1 S);
  const Ref<OrbitalSpace>& K_i_p(SpinCase1 S);
  const Ref<OrbitalSpace>& K_i_m(SpinCase1 S);
  const Ref<OrbitalSpace>& K_i_a(SpinCase1 S);
  const Ref<OrbitalSpace>& K_m_a(SpinCase1 S);
  const Ref<OrbitalSpace>& K_a_a(SpinCase1 S);
  const Ref<OrbitalSpace>& K_a_p(SpinCase1 S);
  const Ref<OrbitalSpace>& K_a_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_p_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_p_A(SpinCase1 S);
  const Ref<OrbitalSpace>& K_p_p(SpinCase1 S);
  const Ref<OrbitalSpace>& K_p_m(SpinCase1 S);
  const Ref<OrbitalSpace>& K_p_a(SpinCase1 S);
  const Ref<OrbitalSpace>& K_A_P(SpinCase1 S);
  const Ref<OrbitalSpace>& K_P_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_x_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_x_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_x_p(SpinCase1 S);
  const Ref<OrbitalSpace>& F_x_m(SpinCase1 S);
  const Ref<OrbitalSpace>& F_x_a(SpinCase1 S);
  const Ref<OrbitalSpace>& F_i_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_i_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_i_p(SpinCase1 S);
  const Ref<OrbitalSpace>& F_i_m(SpinCase1 S);
  const Ref<OrbitalSpace>& F_i_a(SpinCase1 S);
  const Ref<OrbitalSpace>& F_m_m(SpinCase1 S);
  const Ref<OrbitalSpace>& F_m_a(SpinCase1 S);
  const Ref<OrbitalSpace>& F_m_p(SpinCase1 S);
  const Ref<OrbitalSpace>& F_m_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_gg_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_m_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_a_a(SpinCase1 S);
  const Ref<OrbitalSpace>& F_a_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_A_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_p_P(SpinCase1 S);
  const Ref<OrbitalSpace>& F_p_A(SpinCase1 S);
  const Ref<OrbitalSpace>& F_p_p(SpinCase1 S);
  const Ref<OrbitalSpace>& F_p_m(SpinCase1 S);
  const Ref<OrbitalSpace>& F_p_a(SpinCase1 S);
  const Ref<OrbitalSpace>& F_P_P(SpinCase1 S);
  const Ref<OrbitalSpace>& h_P_P(SpinCase1 S);
  const Ref<OrbitalSpace>& J_x_p(SpinCase1 S);
  const Ref<OrbitalSpace>& J_i_p(SpinCase1 S);
  const Ref<OrbitalSpace>& J_i_P(SpinCase1 S);
  const Ref<OrbitalSpace>& J_P_P(SpinCase1 S);
 
  const Ref<OrbitalSpace>& gamma_p_p(SpinCase1 S);
    const Ref<OrbitalSpace>& gamma_p_p_av();
    const Ref<OrbitalSpace>& gamma_m_m_av();
  const Ref<OrbitalSpace>& gammaFgamma_p_p(SpinCase1 S);
  const Ref<OrbitalSpace>& gammaFgamma_p_p();
  const Ref<OrbitalSpace>& Fgamma_p_P(SpinCase1 S);
  const Ref<OrbitalSpace>& Fgamma_p_P();
//  const Ref<OrbitalSpace>& Fgamma_P_p(SpinCase1 S,const RefSymmSCMatrix &gamma,Ref<OrbitalSpace>& FGamma);
//  const Ref<OrbitalSpace>& gammaF_p_P(SpinCase1 S,const RefSymmSCMatrix &gamma,Ref<OrbitalSpace>& GammaF);
  Ref<OrbitalSpace> obtensor_p_A(const RefSCMatrix &obtensor,SpinCase1 S);
 
  std::string transform_label(const Ref<OrbitalSpace>& space1,
                              const Ref<OrbitalSpace>& space2,
                              const Ref<OrbitalSpace>& space3,
                              const Ref<OrbitalSpace>& space4,
                              const std::string& operator_label = std::string()) const;
  std::string transform_label(const Ref<OrbitalSpace>& space1,
                              const Ref<OrbitalSpace>& space2,
                              const Ref<OrbitalSpace>& space3,
                              const Ref<OrbitalSpace>& space4,
                              unsigned int f12,
                              const std::string& operator_label = std::string()) const;
  std::string transform_label(const Ref<OrbitalSpace>& space1,
                              const Ref<OrbitalSpace>& space2,
                              const Ref<OrbitalSpace>& space3,
                              const Ref<OrbitalSpace>& space4,
                              unsigned int f12_left,
                              unsigned int f12_right,
                              const std::string& operator_label = std::string()) const;
 
  void compute_T2_(RefSCMatrix& T2,
                   const Ref<OrbitalSpace>& space1,
                   const Ref<OrbitalSpace>& space2,
                   const Ref<OrbitalSpace>& space3,
                   const Ref<OrbitalSpace>& space4,
                   bool antisymmetrize,
                   const std::string& tform_key);
  void compute_F12_(RefSCMatrix& F12,
                   const Ref<OrbitalSpace>& space1,
                   const Ref<OrbitalSpace>& space2,
                   const Ref<OrbitalSpace>& space3,
                   const Ref<OrbitalSpace>& space4,
                   bool antisymmetrize,
                   const std::vector<std::string>& transform_keys);
 
  RefSCMatrix fock(const Ref<OrbitalSpace>& bra_space,
                   const Ref<OrbitalSpace>& ket_space, SpinCase1 S = Alpha,
                   double scale_J = 1.0, double scale_K = 1.0, double scale_H = 1.0,
                   // leave at default to use R12Technology's pauli flag, else set to 0 or 1
                   int override_pauli = -1);
 
 
  RefSCMatrix V_cabs(SpinCase2 spincase2,
                     const Ref<OrbitalSpace>& p,
                     const Ref<OrbitalSpace>& q);
 
  RefSCMatrix V_genref_spinfree(const Ref<OrbitalSpace>& p,
                     const Ref<OrbitalSpace>& q);
 
#if defined(MPQC_NEW_FEATURES)
 
  void V_diag_ta();
 
  // \param orbital index of the orbital for self-energy calculation, -1 = HOMO (default), +1 = LUMO, 0 = do nothing
  void gf2_r12(int orbital = -1);
 
  // TiledArray-based MP2-F12 one-electron properties
  void compute_TA_mp2f12_1rdm();
#endif
 
  void compute_ccr12_1rdm(const RefSCMatrix& T1, const Ref<DistArray4> (&T2)[NSpinCases2]);
 
  const Ref<OrbitalSpaceRegistry>& orbital_registry() const {
    return this->r12world()->world()->tfactory()->orbital_registry();
  }
  const Ref<AOSpaceRegistry>& ao_registry() const {
    return this->r12world()->world()->tfactory()->ao_registry();
  }
 
};
 
std::vector< Ref<DistArray4> >
A_distarray4(SpinCase2 spincase2, const Ref<R12IntEval>& r12eval);
 
// compute orbital Z-vector from F12 contribution
RefSCMatrix Onerdm_X_F12(SpinCase1 spin, const Ref<R12IntEval>& r12eval, int debug);
 
// compute orbital relaxation contributions from CABS Singles
RefSCMatrix Onerdm_X_CABS_Singles(SpinCase1 spin,
                                  const Ref<R12IntEval>& r12eval,
                                  const Ref<R12EnergyIntermediates>& r12intermediates,
                                  int debug);
 
} // end of namespace sc
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: