DMRGSCFunitary.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 DMRGSCFUNITARY_CHEMPS2_H
#define DMRGSCFUNITARY_CHEMPS2_H

#include "Options.h"
#include "DMRGSCFindices.h"
#include "DMRGSCFmatrix.h"
#include "DIIS.h"

namespace CheMPS2{
   class DMRGSCFunitary : public DMRGSCFmatrix{

      public:
      

         DMRGSCFunitary( const DMRGSCFindices * iHandler );
         
         virtual ~DMRGSCFunitary();


         int getNumVariablesX() const;


         void updateUnitary( double * workmem1, double * workmem2, double * vector, const bool multiply, const bool compact );
         

         void rotateActiveSpaceVectors( double * eigenvecs, double * work );
         

         void CheckDeviationFromUnitary( double * work ) const;
         

         void getLog( double * vector, double * temp1, double * temp2 ) const;


         void makeSureAllBlocksDetOne( double * temp1, double * temp2 );


         void saveU( const string filename=DMRGSCF_unitary_storage_name ) const;


         void loadU( const string filename=DMRGSCF_unitary_storage_name );

      private:

         //Number of variables in the x-matrix
         int x_linearlength;

         //Helper arrays to jump from linear x-matrix index to orbital indices and back
         int ** jumper;

         // Build in result the skew symmetric matrix X for irrep block irrep based on the elements in Xelem. If compact==true, they are stored in gradient form.
         void buildSkewSymmX( const int irrep, double * result, double * Xelem, const bool compact ) const;

         // Get the determinant, and in the process fill work1 with eigvec( U[ irrep ] + U^T[ irrep ] ) and work2 with work1^T U work1 (TRIDIAGONAL matrix).
         double get_determinant( const int irrep, double * work1, double * work2, double * work_eig, int lwork_eig ) const;

   };
}

#endif