DMRG.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 DMRG_CHEMPS2_H
#define DMRG_CHEMPS2_H

#include <string>

#include "Options.h"
#include "Problem.h"
#include "TensorT.h"
#include "TensorX.h"
#include "TensorL.h"
#include "SyBookkeeper.h"
#include "TensorF1.h"
#include "TensorF0.h"
#include "TensorS0.h"
#include "TensorS1.h"
#include "TensorOperator.h"
#include "TensorQ.h"
#include "TensorO.h"
#include "Tensor3RDM.h"
#include "TensorGYZ.h"
#include "TensorKM.h"
#include "TwoDM.h"
#include "ThreeDM.h"
#include "Correlations.h"
#include "Heff.h"
#include "Sobject.h"
#include "ConvergenceScheme.h"
#include "MyHDF5.h"

//For the timings of the different parts of DMRG
#define CHEMPS2_TIME_S_JOIN      0
#define CHEMPS2_TIME_S_SOLVE     1
#define CHEMPS2_TIME_S_SPLIT     2
#define CHEMPS2_TIME_TENS_TOTAL  3
#define CHEMPS2_TIME_TENS_ALLOC  4
#define CHEMPS2_TIME_TENS_FREE   5
#define CHEMPS2_TIME_DISK_WRITE  6
#define CHEMPS2_TIME_DISK_READ   7
#define CHEMPS2_TIME_TENS_CALC   8
#define CHEMPS2_TIME_VECLENGTH   9

namespace CheMPS2{
   class DMRG{

      public:
      

         DMRG(Problem * Probin, ConvergenceScheme * OptSchemeIn, const bool makechkpt=CheMPS2::DMRG_storeMpsOnDisk, const string tmpfolder=CheMPS2::defaultTMPpath);
         
         virtual ~DMRG();
         

         double Solve();
         
         void PreSolve();
         
         void calc2DMandCorrelations(){ calc_rdms_and_correlations(false); }
         

         void calc_rdms_and_correlations( const bool do_3rdm, const bool disk_3rdm = false );
         

         TwoDM * get2DM(){ return the2DM; }
         

         ThreeDM * get3DM(){ return the3DM; }


         void Symm4RDM( double * output, const int ham_orb1, const int ham_orb2, const bool last_case );


         Correlations * getCorrelations(){ return theCorr; }
         

         double getSpecificCoefficient(int * coeff) const;
         

         double getFCIcoefficient(int * alpha, int * beta, const bool mpi_chemps2_master_only=true) const;
         
         void deleteStoredMPS();
         
         void deleteStoredOperators();
         

         void activateExcitations(const int maxExcIn);
         

         void newExcitation(const double EshiftIn);
         
         static void PrintLicense();
         
      private:
      
         //Setup the DMRG SyBK and MPS (in separate function to allow pushbacks and recreations for excited states)
         void setupBookkeeperAndMPS();
      
         string MPSstoragename;
         
         //The optimization scheme for the DMRG sweeps (externally allocated, filled and deleted)
         ConvergenceScheme * OptScheme;
         
         //Whether or not the MPS was loaded from disk to memory at the start
         bool loadedMPS;
         
         //Integer to distinguish storage between different calculations
         int thePID;
      
         //Pointer to the Problem --> constructed and destructed outside of this class
         Problem * Prob;
         
         //The number of orbitals: copied here so that the DMRG destructor doesn't depend on whether Prob still exists
         int L;
         
         //Minimum energy encountered during all performed micro-iterations (as opposed to the 2DM/edge energy)
         double TotalMinEnergy;
         
         //Minimum energy encountered during the micro-iterations of the last performed sweep (as opposed to the 2DM/edge energy)
         double LastMinEnergy;
         
         //Max. discarded weight of last sweep
         double MaxDiscWeightLastSweep;
         
         //Symmetry information object
         SyBookkeeper * denBK;
         
         //The MPS
         TensorT ** MPS;
         
         //The TwoDM
         TwoDM * the2DM;
         
         //The ThreeDM
         ThreeDM * the3DM;
         
         //The Correlations
         Correlations * theCorr;
         
         //Whether or not allocated
         int * isAllocated;
         
         //TensorL's
         TensorL *** Ltensors;
         
         //TensorX's
         TensorX ** Xtensors;
         
         //TensorF0's
         TensorF0 **** F0tensors;
         
         //TensorF1's
         TensorF1 **** F1tensors;
         
         //TensorS0's
         TensorS0 **** S0tensors;
         
         //TensorS1's
         TensorS1 **** S1tensors;
         
         //ABCD-tensors
         TensorOperator **** Atensors;
         TensorOperator **** Btensors;
         TensorOperator **** Ctensors;
         TensorOperator **** Dtensors;
         
         //TensorQ's
         TensorQ *** Qtensors;
         
         //Tensors required for the 3-RDM calculation
         Tensor3RDM ***** tensor_3rdm_a_J0_doublet;
         Tensor3RDM ***** tensor_3rdm_a_J1_doublet;
         Tensor3RDM ***** tensor_3rdm_a_J1_quartet;
         Tensor3RDM ***** tensor_3rdm_b_J0_doublet;
         Tensor3RDM ***** tensor_3rdm_b_J1_doublet;
         Tensor3RDM ***** tensor_3rdm_b_J1_quartet;
         Tensor3RDM ***** tensor_3rdm_c_J0_doublet;
         Tensor3RDM ***** tensor_3rdm_c_J1_doublet;
         Tensor3RDM ***** tensor_3rdm_c_J1_quartet;
         Tensor3RDM ***** tensor_3rdm_d_J0_doublet;
         Tensor3RDM ***** tensor_3rdm_d_J1_doublet;
         Tensor3RDM ***** tensor_3rdm_d_J1_quartet;
         
         //Tensors required for the Correlations calculation
         TensorGYZ ** Gtensors;
         TensorGYZ ** Ytensors;
         TensorGYZ ** Ztensors;
         TensorKM  ** Ktensors;
         TensorKM  ** Mtensors;

         // Sweeps
         double sweepleft(  const bool change, const int instruction, const bool am_i_master );
         double sweepright( const bool change, const int instruction, const bool am_i_master );
         double solve_site( const int index, const double dvdson_rtol, const double noise_level, const int virtual_dimension, const bool am_i_master, const bool moving_right, const bool change );

         //Load and save functions
         void MY_HDF5_WRITE_BATCH(const hid_t file_id, const int number, Tensor ** batch, const long long totalsize, const std::string tag);
         void MY_HDF5_READ_BATCH( const hid_t file_id, const int number, Tensor ** batch, const long long totalsize, const std::string tag);
         void OperatorsOnDisk(const int index, const bool movingRight, const bool store);
         string tempfolder;
         
         void saveMPS(const std::string name, TensorT ** MPSlocation, SyBookkeeper * BKlocation, bool isConverged) const;
         void loadDIM(const std::string name, SyBookkeeper * BKlocation);
         void loadMPS(const std::string name, TensorT ** MPSlocation, bool * isConverged);
         bool makecheckpoints;
         
         //Helper functions for making the boundary operators
         void updateMovingRight(const int index);
         void updateMovingLeft(const int index);
         void deleteTensors(const int index, const bool movingRight);
         void allocateTensors(const int index, const bool movingRight);
         void updateMovingRightSafe(const int cnt);
         void updateMovingRightSafeFirstTime(const int cnt);
         void updateMovingRightSafe2DM(const int cnt);
         void updateMovingLeftSafe(const int cnt);
         void updateMovingLeftSafeFirstTime(const int cnt);
         void updateMovingLeftSafe2DM(const int cnt);
         void deleteAllBoundaryOperators();

         // Helper functions for making the 3-RDM boundary operators
         void update_safe_3rdm_operators( const int boundary );
         void allocate_3rdm_operators( const int boundary );
         void update_3rdm_operators( const int boundary );
         void delete_3rdm_operators( const int boundary );
         
         // Helper functions for making the symmetrized 4-RDM
         void solve_fock( const int dmrg_orb1, const int dmrg_orb2, const double alpha, const double beta );
         static void solve_fock_update_helper( const int index, const int dmrg_orb1, const int dmrg_orb2, const bool moving_right, TensorT ** new_mps, TensorT ** old_mps, SyBookkeeper * new_bk, SyBookkeeper * old_bk, TensorO ** overlaps, TensorL ** regular, TensorL ** trans );
         static void  left_normalize( TensorT * left_mps, TensorT * right_mps );
         static void right_normalize( TensorT * left_mps, TensorT * right_mps );
         void symm_4rdm_helper( double * output, const int ham_orb1, const int ham_orb2, const double alpha, const double beta, const bool add, const double factor );

         //Helper functions for making the Correlations boundary operators
         void update_correlations_tensors(const int siteindex);

         //The storage and functions to handle excited states
         int nStates;
         bool Exc_activated;
         int maxExc;
         double * Exc_Eshifts;
         TensorT *** Exc_MPSs;
         SyBookkeeper ** Exc_BKs;
         TensorO *** Exc_Overlaps;
         double ** prepare_excitations(Sobject * denS);
         void cleanup_excitations(double ** VeffTilde) const;
         void calcVeffTilde(double * result, Sobject * currentS, int state_number);
         void calc_overlaps( const bool moving_right );
         
         // Performance counters
         double timings[ CHEMPS2_TIME_VECLENGTH ];
         long long num_double_write_disk;
         long long num_double_read_disk;
         void print_tensor_update_performance() const;
         
   };
}

#endif