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

#ifdef __GNUG__
#pragma interface
#endif

#include <util/ref/ref.h>
#include <chemistry/qc/mbptr12/vxb_eval_info.h>
#include <chemistry/qc/mbptr12/linearr12.h>
#include <chemistry/qc/mbptr12/twoparticlecontraction.h>
#include <chemistry/qc/mbptr12/spin.h>

#ifndef _chemistry_qc_mbptr12_r12inteval_h
#define _chemistry_qc_mbptr12_r12inteval_h

namespace sc {

  class TwoBodyMOIntsTransform;
  class R12IntsAcc;
  class F12Amplitudes;

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<R12IntEvalInfo> r12info_;

  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];

  RefSCVector emp2pair_[NSpinCases2];
  RefSCDimension dim_oo_[NSpinCases2];
  RefSCDimension dim_vv_[NSpinCases2];
  RefSCDimension dim_aa_[NSpinCases2];
  RefSCDimension dim_f12_[NSpinCases2];
  RefSCDimension dim_xy_[NSpinCases2];

  Ref<F12Amplitudes> Amps_;  // First-order amplitudes of various contributions to the pair functions
  RefSCDimension dim_ij_s_, dim_ij_t_;
  double emp2_singles_;
  int debug_;

  // Map to TwoBodyMOIntsTransform objects that have been computed previously
  typedef std::map<std::string, Ref<TwoBodyMOIntsTransform> > TformMap;
  TformMap tform_map_;

  void spinadapt_mospace_labels(SpinCase1 spin, std::string& id, std::string& name) const;

  void form_canonvir_space_();

  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_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_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_a_a_[NSpinCases1];
  Ref<OrbitalSpace> F_a_A_[NSpinCases1];

  void init_tforms_();
  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);
  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);
  // 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 Ref<OrbitalSpace>& occ_space, 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_a_vbsneqobs_();

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

  template <typename DataProcess, bool CorrFactorInBra, bool CorrFactorInKet>
    void compute_tbint_tensor(RefSCMatrix& T,
                              TwoBodyInt::tbint_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);

  template <typename DataProcessBra,
            typename DataProcessKet,
            typename DataProcessBraKet,
            bool CorrFactorInBra,
            bool CorrFactorInKet,
            bool CorrFactorInInt>
    void contract_tbint_tensor(
           RefSCMatrix& T,
           TwoBodyInt::tbint_type tbint_type_bra,
           TwoBodyInt::tbint_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<LinearR12::TwoParticleContraction>& tpcontract,
           bool antisymmetrize,
           const std::vector<std::string>& tformkeys_bra,
           const std::vector<std::string>& tformkeys_ket);

  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_bc_();

  void compute_BB_();

  void compute_BC_();

  void compute_BApp_();

  void compute_B_DKH_();

  void compute_singles_emp2_();

  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<R12IntEvalInfo>& info);
  /*R12IntEval(const Ref<R12IntEvalInfo>& info, bool gbc = true, bool ebc = true,
             LinearR12::ABSMethod abs_method = LinearR12::ABS_CABSPlus,
             LinearR12::StandardApproximation stdapprox = LinearR12::StdApprox_Ap);*/
  ~R12IntEval();

  void save_data_state(StateOut&);
  virtual void obsolete();

  void set_debug(int debug);
  void set_dynamic(bool dynamic);
  void set_print_percent(double print_percent);
  void set_memory(size_t nbytes);

  int dk() const;

  const Ref<LinearR12::CorrelationFactor>& corrfactor() const { return r12info()->corrfactor(); }
  LinearR12::ABSMethod abs_method() const { return r12info()->abs_method(); }
  const Ref<LinearR12Ansatz>& ansatz() const { return r12info()->ansatz(); }
  bool spin_polarized() const { return r12info()->refinfo()->spin_polarized(); }
  bool gbc() const { return r12info()->gbc(); }
  bool ebc() const { return r12info()->ebc(); }
  LinearR12::StandardApproximation stdapprox() const { return r12info()->stdapprox(); }
  bool omit_P() const { return r12info()->omit_P(); }

  const Ref<R12IntEvalInfo>& r12info() const;

  RefSCDimension dim_oo_s() const;
  RefSCDimension dim_oo_t() const;
  RefSCDimension dim_oo(SpinCase2 S) const;
  RefSCDimension dim_vv(SpinCase2 S) const;
  RefSCDimension dim_aa(SpinCase2 S) const;
  RefSCDimension dim_f12(SpinCase2 S) const;
  RefSCDimension dim_xy(SpinCase2 S) const;
  int nspincases1() const { return ::sc::nspincases1(spin_polarized()); }
  int nspincases2() const { return ::sc::nspincases2(spin_polarized()); }

  virtual void compute();

  Ref<F12Amplitudes> amps();

  const RefSCMatrix& V(SpinCase2 S);
  RefSymmSCMatrix X(SpinCase2 S);
  RefSymmSCMatrix B(SpinCase2 S);
  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);
  RefSymmSCMatrix P(SpinCase2 S);

  double emp2_singles();
  const RefSCVector& emp2(SpinCase2 S);
  const RefDiagSCMatrix& evals(SpinCase1 S) const;

  RefDiagSCMatrix evals() const;
  RefDiagSCMatrix evals_a() const;
  RefDiagSCMatrix evals_b() 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>& xspace(SpinCase1 S) 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_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>& 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_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);

  Ref<TwoBodyMOIntsTransform> get_tform_(const std::string&) const;
  void add_tform(const std::string& label,
                 const Ref<TwoBodyMOIntsTransform>& T);
  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);

};

class TransformNotFound: public ProgrammingError {
  public:
  TransformNotFound(const char *description=0, const char *file=0, int line=0) : ProgrammingError(description,file,line) {}
};

}

#endif

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

---