FCI.h

/*
   CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
   Copyright (C) 2013-2016 Sebastian Wouters

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program 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 General Public License for more details.

   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/

#ifndef FCI_CHEMPS2_H
#define FCI_CHEMPS2_H

#include "Hamiltonian.h"

namespace CheMPS2{
   class FCI{

      public:
         

         FCI(CheMPS2::Hamiltonian * Ham, const unsigned int Nel_up, const unsigned int Nel_down, const int TargetIrrep, const double maxMemWorkMB=100.0, const int FCIverbose=2);
         
         virtual ~FCI();
         
//==========> Getters of basic information: all these variables are set by the constructor
         

         unsigned int getL() const{ return L; }


         unsigned int getNel_up() const{ return Nel_up; }


         unsigned int getNel_down() const{ return Nel_down; }
         

         unsigned int getVecLength(const int irrep_center) const{ return irrep_center_jumps[ irrep_center ][ num_irreps ]; }
         

         int getTargetIrrep() const{ return TargetIrrep; }
         

         int getOrb2Irrep(const int orb) const{ return orb2irrep[ orb ]; }
         

         double getERI(const int orb1, const int orb2, const int orb3, const int orb4) const{ return ERI[ orb1 + L * ( orb2 + L * ( orb3 + L * orb4 ) ) ]; }
         

         double getGmat(const int orb1, const int orb2) const{ return Gmat[ orb1 + L * orb2 ]; }
         

         double getEconst() const{ return Econstant; }
         
//==========> The core routines for users
         

         double GSDavidson(double * inoutput=NULL, const int DVDSN_NUM_VEC=CheMPS2::DAVIDSON_NUM_VEC) const;
         

         unsigned int LowestEnergyDeterminant() const;
         

         double Fill2RDM(double * vector, double * TwoRDM) const;
         

         void Fill3RDM(double * vector, double * ThreeRDM) const;
         

         void Fill4RDM(double * vector, double * FourRDM) const;
         

         void Fock4RDM(double * vector, double * ThreeRDM, double * Fock, double * output) const;
         

         void Diag4RDM( double * vector, double * three_rdm, const unsigned int orbz, double * output ) const;
         

         double CalcSpinSquared(double * vector) const;


         static void FillRandom(const unsigned int vecLength, double * vec);
         

         static void ClearVector(const unsigned int vecLength, double * vec);
         
//==========> Green's functions functionality
         

         void RetardedGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF) const;
         

         void RetardedGF_addition(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMadd=NULL) const;
         

         void GFmatrix_addition(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal=NULL, double ** TwoRDMimag=NULL, double ** TwoRDMadd=NULL) const;
         

         void RetardedGF_removal(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMrem=NULL) const;
         

         void GFmatrix_removal(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal=NULL, double ** TwoRDMimag=NULL, double ** TwoRDMrem=NULL) const;
         

         void DensityResponseGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF) const;
         

         void DensityResponseGF_forward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMdens=NULL) const;
         

         void DensityResponseGF_backward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMdens=NULL) const;
         

         void CGSolveSystem(const double alpha, const double beta, const double eta, double * RHS, double * RealSol, double * ImagSol, const bool checkError=true) const;
         
         //void CheckHamDEBUG() const;
         

         double getFCIcoeff(int * bits_up, int * bits_down, double * vector) const;
         
      protected:
      
//==========> Functions involving Hamiltonian matrix elements
         

         void DiagHam(double * diag) const;
         

         void DiagHamSquared(double * output) const;
      

         void matvec( double * input, double * output ) const;
         

         double GetMatrixElement(int * bits_bra_up, int * bits_bra_down, int * bits_ket_up, int * bits_ket_down, int * work) const;
         
//==========> Basic conversions between bit string representations
      

         void getBitsOfCounter(const int irrep_center, const unsigned int counter, int * bits_up, int * bits_down) const;
         

         static void str2bits(const unsigned int Lvalue, const unsigned int bitstring, int * bits);
         

         static unsigned int bits2str(const unsigned int Lvalue, int * bits);
         

         int getUpIrrepOfCounter(const int irrep_center, const unsigned int counter) const;
         
//==========> Some lapack like routines


         static double FCIddot(const unsigned int vecLength, double * vec1, double * vec2);
         

         static void FCIdcopy(const unsigned int vecLength, double * origin, double * target);
         

         static double FCIfrobeniusnorm(const unsigned int vecLength, double * vec);
         

         static void FCIdaxpy(const unsigned int vecLength, const double alpha, double * vec_x, double * vec_y);
         

         static void FCIdscal(const unsigned int vecLength, const double alpha, double * vec);
         
//==========> Protected functions regarding the Green's functions
         

         void ActWithSecondQuantizedOperator(const char whichOperator, const bool isUp, const unsigned int orbIndex, double * thisVector, const FCI * otherFCI, double * otherVector) const;
         

         void ActWithNumberOperator(const unsigned int orbIndex, double * resultVector, double * sourceVector) const;


         void CGCoreSolver(const double alpha, const double beta, const double eta, double * precon, double * Sol, double * RESID, double * PVEC, double * OxPVEC, double * temp, double * temp2) const;


         void CGAlphaPlusBetaHAM(const double alpha, const double beta, double * in, double * out) const;
         

         void CGoperator(const double alpha, const double beta, const double eta, double * in, double * temp, double * out) const;
         

         void CGdiagonal(const double alpha, const double beta, const double eta, double * diagonal, double * workspace) const;
         
//==========> Protected functions regarding the higher order RDMs


         void apply_excitation( double * orig_vector, double * result_vector, const int crea, const int anni, const int orig_target_irrep ) const;
      
      private:
      
         int FCIverbose;
         
         double maxMemWorkMB;
         
         double Econstant;
         
         double * Gmat;
         
         double * ERI;
         
         unsigned int num_irreps;
         
         int TargetIrrep;
         
         int * orb2irrep;
         
         unsigned int L;
         
         unsigned int Nel_up;
         
         unsigned int Nel_down;
         
         unsigned int * numPerIrrep_up;
         
         unsigned int * numPerIrrep_down;
         
         int ** str2cnt_up;
         
         int ** str2cnt_down;
         
         unsigned int ** cnt2str_up;
         
         unsigned int ** cnt2str_down;
         
         int *** lookup_cnt_alpha;
         
         int *** lookup_cnt_beta;
         
         int *** lookup_sign_alpha;
         
         int *** lookup_sign_beta;
         
         unsigned int * irrep_center_num;
         
         unsigned int ** irrep_center_crea_orb;
         
         unsigned int ** irrep_center_anni_orb;
         
         unsigned int ** irrep_center_jumps;
         
         unsigned long long HXVsizeWorkspace;
         
         double * HXVworksmall;
         
         double * HXVworkbig1;
         
         double * HXVworkbig2;
         
         void StartupCountersVsBitstrings();
         
         void StartupLookupTables();
         
         void StartupIrrepCenter();
         
         double Driver3RDM(double * vector, double * output, double * three_rdm, double * fock, const unsigned int orbz) const;

         static void excite_alpha_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int dim_down, double * origin, double * result, int * signmap, int * countmap );
         static void excite_alpha_first( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap );
         static void excite_alpha_second_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap );

         static void excite_beta_omp( const unsigned int dim_up, const unsigned int dim_new_down, double * origin, double * result, int * signmap, int * countmap );
         static void excite_beta_first( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap );
         static void excite_beta_second_omp( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap );

   };

}

#endif