mp2r12_energy_util.h

//
// mp2r12_energy_util.h
//
// Copyright (C) 2006 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_mp2r12energyutil_h
#define _chemistry_qc_mbptr12_mp2r12energyutil_h
 
#include <util/ref/ref.h>
#include <math/scmat/matrix.h>
#include <util/misc/scexception.h>
#include <chemistry/qc/wfn/spin.h>
#include <math/scmat/svd.h>
#include <chemistry/qc/wfn/orbitalspace.h>
#include <chemistry/qc/mbptr12/r12technology.h>
 
namespace sc {
 
  class MP2R12EnergyUtil_base : virtual public RefCount {
    protected:
      unsigned int nocc_act_;
      RefSCDimension oodim_;
      RefSCDimension xydim_;
      RefSCDimension f12dim_;
      RefSCDimension gdim_;
      unsigned int nf12_;
 
    public:
    MP2R12EnergyUtil_base();
    MP2R12EnergyUtil_base(const RefSCDimension& oodim,
                          const RefSCDimension& xydim,
                          const RefSCDimension& f12dim,
                          const unsigned int nocc_act);
    virtual ~MP2R12EnergyUtil_base();
 
    void check_dims(const RefSCMatrix& A) const;
    void check_dims(const RefSymmSCMatrix& A) const;
 
    virtual void print(const char* label, const RefSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
    virtual void print(const char* label, const RefSymmSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
    virtual void print(const char* label, const RefDiagSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
    virtual void invert(RefSymmSCMatrix& A) const = 0;
    virtual RefDiagSCMatrix eigenvalues(const RefSymmSCMatrix& A) const = 0;
    virtual void diagonalize(const RefSymmSCMatrix& A, RefDiagSCMatrix& evals, RefSCMatrix& evecs) const = 0;
    virtual void transform(const RefSymmSCMatrix& B, const RefDiagSCMatrix& A, const RefSCMatrix& U) const = 0;
    virtual void solve_linear_system(const RefSymmSCMatrix& A,
                                     RefSCMatrix& X,
                                     const RefSCMatrix& B) const = 0;
    virtual void solve_linear_system(unsigned int ij,
                                     const RefSymmSCMatrix& A,
                                     RefSCMatrix& X,
                                     const RefSCMatrix& B) const = 0;
    virtual void times(const RefSymmSCMatrix& A,
                       const RefSCMatrix& x,
                       RefSCMatrix& y) const = 0;
    virtual void times(unsigned int ij,
                       const RefSymmSCMatrix& A,
                       const RefSCMatrix& x,
                       RefSCMatrix& y) const = 0;
    virtual RefSCVector dot_product(const RefSCMatrix& A,
                                    const RefSCMatrix& B) const = 0;
  };
 
  class MP2R12EnergyUtil_Diag : virtual public RefCount {
    protected:
      unsigned int nocc_act_;
      RefSCDimension oodim_;
      RefSCDimension xydim_;
      RefSCDimension f12dim_;
      RefSCDimension gdim_;
      unsigned int nf12_;
 
    public:
    MP2R12EnergyUtil_Diag();
    MP2R12EnergyUtil_Diag(const RefSCDimension& oodim,
                          const RefSCDimension& xydim,
                          const RefSCDimension& f12dim,
                          const unsigned int nocc_act);
    virtual ~MP2R12EnergyUtil_Diag();
 
    void check_dims(const RefSCMatrix& A) const;
    void check_dims(const RefSymmSCMatrix& A) const;
 
    virtual void print(const char* label, const RefSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
    virtual void print(const char* label, const RefSymmSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
    virtual void print(const char* label, const RefDiagSCMatrix& A, std::ostream& os = ExEnv::out0()) const = 0;
 
    virtual void get(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const = 0;
    virtual void get(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const =0;
    virtual void get(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const = 0;
    virtual void get(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const =0;
 
    virtual void put(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const =0;
    virtual void put(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const =0;
    virtual void put(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const = 0;
    virtual void put(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const = 0;
 
    virtual void invert(RefSymmSCMatrix& A) const = 0;
    virtual RefDiagSCMatrix eigenvalues(const RefSymmSCMatrix& A) const = 0;
    virtual void diagonalize(const RefSymmSCMatrix& A, RefDiagSCMatrix& evals, RefSCMatrix& evecs) const = 0;
    virtual void transform(const RefSymmSCMatrix& B, const RefDiagSCMatrix& A, const RefSCMatrix& U) const = 0;
    virtual void solve_linear_system(const RefSymmSCMatrix& A,
                                     RefSCMatrix& X,
                                     const RefSCMatrix& B) const = 0;
    virtual void solve_linear_system(unsigned int ij,
                                     const RefSymmSCMatrix& A,
                                     RefSCMatrix& X,
                                     const RefSCMatrix& B) const = 0;
    virtual void times(const RefSymmSCMatrix& A,
                       const RefSCMatrix& x,
                       RefSCMatrix& y) const = 0;
    virtual void times(unsigned int ij,
                       const RefSymmSCMatrix& A,
                       const RefSCMatrix& x,
                       RefSCMatrix& y) const = 0;
    virtual RefSCVector dot_product(const RefSCMatrix& A,
                                    const RefSCMatrix& B) const = 0;
  };
 
  class MP2R12EnergyUtil_Diag_DifferentSpin : public MP2R12EnergyUtil_Diag {
    private:
      SpinCase2 spincase2_;
    public:
      MP2R12EnergyUtil_Diag_DifferentSpin(const RefSCDimension& oodim,
                                           const RefSCDimension& xydim,
                                           const RefSCDimension& f12dim,
                                           const unsigned int nocc_act);
      ~MP2R12EnergyUtil_Diag_DifferentSpin() {}
 
      unsigned int nrowblks(const RefSCMatrix& A) const;
      unsigned int ncolblks(const RefSCMatrix& A) const;
      unsigned int nblks(const RefSymmSCMatrix& A) const;
 
      void get(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const;
      void get(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const;
      void get(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const;
      void get(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const;
 
      void put(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const;
      void put(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const;
      void put(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const;
      void put(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const;
 
      void invert(RefSymmSCMatrix& A) const;
      RefDiagSCMatrix eigenvalues(const RefSymmSCMatrix& A) const;
      void diagonalize(const RefSymmSCMatrix& A, RefDiagSCMatrix& evals, RefSCMatrix& evecs) const;
      void transform(const RefSymmSCMatrix& B, const RefDiagSCMatrix& A, const RefSCMatrix& U) const;
      void solve_linear_system(const RefSymmSCMatrix& A,
                               RefSCMatrix& X,
                               const RefSCMatrix& B) const;
      void solve_linear_system(unsigned int ij,
                               const RefSymmSCMatrix& A,
                               RefSCMatrix& X,
                               const RefSCMatrix& B) const;
      void times(const RefSymmSCMatrix& A,
                 const RefSCMatrix& x,
                 RefSCMatrix& y) const;
      void times(unsigned int ij,
                 const RefSymmSCMatrix& A,
                 const RefSCMatrix& x,
                 RefSCMatrix& y) const;
      RefSCVector dot_product(const RefSCMatrix& A,
                              const RefSCMatrix& B) const;
 
      void print(const char* label, const RefSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
      void print(const char* label, const RefSymmSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
      void print(const char* label, const RefDiagSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
  };
 
  class MP2R12EnergyUtil_Diag_SameSpin : public MP2R12EnergyUtil_Diag {
    private:
      SpinCase2 spincase2_;
    public:
      MP2R12EnergyUtil_Diag_SameSpin(const RefSCDimension& oodim,
                                      const RefSCDimension& xydim,
                                      const RefSCDimension& f12dim,
                                      const unsigned int nocc_act);
      ~MP2R12EnergyUtil_Diag_SameSpin() {}
 
      unsigned int nrowblks(const RefSCMatrix& A) const;
      unsigned int ncolblks(const RefSCMatrix& A) const;
      unsigned int nblks(const RefSymmSCMatrix& A) const;
 
      void get(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const;
      void get(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const;
      void get(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const;
      void get(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const;
 
      void put(unsigned int ij, const RefSCMatrix& A, const RefSCVector& Aij) const;
      void put(unsigned int ij, const RefSCMatrix& A, const RefSCMatrix& Aij) const;
      void put(unsigned int ij, const RefSymmSCMatrix& A, const RefSymmSCMatrix& Aij) const;
      void put(unsigned int ij, const RefDiagSCMatrix& A, const RefDiagSCMatrix& Aij) const;
 
      void invert(RefSymmSCMatrix& A) const;
      RefDiagSCMatrix eigenvalues(const RefSymmSCMatrix& A) const;
      void diagonalize(const RefSymmSCMatrix& A, RefDiagSCMatrix& evals, RefSCMatrix& evecs) const;
      void transform(const RefSymmSCMatrix& B, const RefDiagSCMatrix& A, const RefSCMatrix& U) const;
      void solve_linear_system(const RefSymmSCMatrix& A,
                               RefSCMatrix& X,
                               const RefSCMatrix& B) const;
      void solve_linear_system(unsigned int ij,
                               const RefSymmSCMatrix& A,
                               RefSCMatrix& X,
                               const RefSCMatrix& B) const;
      void times(const RefSymmSCMatrix& A,
                 const RefSCMatrix& x,
                 RefSCMatrix& y) const;
      void times(unsigned int ij,
                 const RefSymmSCMatrix& A,
                 const RefSCMatrix& x,
                 RefSCMatrix& y) const;
      RefSCVector dot_product(const RefSCMatrix& A,
                              const RefSCMatrix& B) const;
 
      void print(const char* label, const RefSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
      void print(const char* label, const RefSymmSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
      void print(const char* label, const RefDiagSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
  };
 
    class MP2R12EnergyUtil_Nondiag : public MP2R12EnergyUtil_base {
    public:
    MP2R12EnergyUtil_Nondiag(const RefSCDimension& oodim,
                             const RefSCDimension& xydim,
                             const RefSCDimension& f12dim,
                             const unsigned int nocc_act);
    ~MP2R12EnergyUtil_Nondiag() {}
 
    void print(const char* label, const RefSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
    void print(const char* label, const RefSymmSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
    void print(const char* label, const RefDiagSCMatrix& A, std::ostream& os = ExEnv::out0()) const;
    void invert(RefSymmSCMatrix& A) const;
    RefDiagSCMatrix eigenvalues(const RefSymmSCMatrix& A) const;
    void diagonalize(const RefSymmSCMatrix& A, RefDiagSCMatrix& evals, RefSCMatrix& evecs) const;
    void transform(const RefSymmSCMatrix& B, const RefDiagSCMatrix& A, const RefSCMatrix& U) const;
    void solve_linear_system(const RefSymmSCMatrix& A,
                             RefSCMatrix& X,
                             const RefSCMatrix& B) const;
    void solve_linear_system(unsigned int ij,
                             const RefSymmSCMatrix& A,
                             RefSCMatrix& X,
                             const RefSCMatrix& B) const;
    void times(const RefSymmSCMatrix& A,
               const RefSCMatrix& x,
               RefSCMatrix& y) const;
    void times(unsigned int ij,
               const RefSymmSCMatrix& A,
               const RefSCMatrix& x,
               RefSCMatrix& y) const;
    RefSCVector dot_product(const RefSCMatrix& A,
                            const RefSCMatrix& B) const;
 
    unsigned int nrowblks(const RefSCMatrix& A) const;
    unsigned int ncolblks(const RefSCMatrix& A) const;
    unsigned int nblks(const RefSymmSCMatrix& A) const;
 
  };
 
  Ref<MP2R12EnergyUtil_Diag> generate_MP2R12EnergyUtil_Diag(SpinCase2 spincase2,
                                                            const RefSCDimension& oodim,
                                                            const RefSCDimension& xydim,
                                                            const RefSCDimension& f12dim,
                                                            const unsigned int nocc_act);
 
  void firstorder_cusp_coefficients(const SpinCase2 &spincase2,
                                      RefSCMatrix& C,
                                      const Ref<OrbitalSpace>& i1,
                                      const Ref<OrbitalSpace>& i2,
                                      const Ref<R12Technology::CorrelationFactor>& corrfactor);
}
 
#endif // include guard