TwoDM.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 TWODM_CHEMPS2_H
#define TWODM_CHEMPS2_H

#include "TensorT.h"
#include "TensorL.h"
#include "TensorF0.h"
#include "TensorF1.h"
#include "TensorS0.h"
#include "TensorS1.h"
#include "SyBookkeeper.h"

namespace CheMPS2{
   class TwoDM{

      public:
      

         TwoDM(const SyBookkeeper * denBKIn, const Problem * ProbIn);
         
         virtual ~TwoDM();
         

         double getTwoDMA_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const;
         

         double getTwoDMB_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const;
         

         double get1RDM_DMRG(const int cnt1, const int cnt2) const;
         

         double spin_density_dmrg( const int cnt1, const int cnt2 ) const;
         

         double getTwoDMA_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const;
         

         double getTwoDMB_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const;
         

         double get1RDM_HAM(const int cnt1, const int cnt2) const;
         

         double spin_density_ham( const int cnt1, const int cnt2 ) const;
         

         void FillSite(TensorT * denT, TensorL *** Ltens, TensorF0 **** F0tens, TensorF1 **** F1tens, TensorS0 **** S0tens, TensorS1 **** S1tens);
         
         void correct_higher_multiplicities();
             

         double trace() const;
         

         double energy() const;
         
         void print_noon() const;
         
         void save() const;
         
         void read();


         void save_HAM( const string filename ) const;


         void write2DMAfile(const string filename) const;
         
         void mpi_allreduce();
         
      private:
      
         //The BK containing all the irrep information
         const SyBookkeeper * denBK;
         
         //The problem containing orbital reshuffling and symmetry information
         const Problem * Prob;
         
         //The chain length
         int L;
         
         //Two 2DM^{A,B} objects
         double * two_rdm_A;
         double * two_rdm_B;
         
         // Set 2DM terms, using the DMRG indices
         void set_2rdm_A_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value);
         void set_2rdm_B_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value);
         
         //Helper functions
         double doD1(TensorT * denT);
         double doD2(TensorT * denT, TensorL * Lright, double * workmem);
         double doD3(TensorT * denT, TensorS0 * S0right, double * workmem);
         double doD4(TensorT * denT, TensorF0 * F0right, double * workmem);
         double doD5(TensorT * denT, TensorF0 * F0right, double * workmem);
         double doD6(TensorT * denT, TensorF1 * F1right, double * workmem);
         double doD7(TensorT * denT, TensorL * Lleft, double * workmem);
         double doD8(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g);
         void doD9andD10andD11(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, double * d9, double * d10, double * d11, int Irrep_g);
         double doD12(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g);
         double doD13(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g);
         double doD14(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g);
         double doD15(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g);
         double doD16(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g);
         double doD17orD21(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD17);
         double doD18orD22(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD18);
         double doD19orD23(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD19);
         double doD20orD24(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD20);
         
   };
}

#endif