CASSCF.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 CASSCF_CHEMPS2_H
#define CASSCF_CHEMPS2_H

#include <string>

#include "Hamiltonian.h"
#include "Irreps.h"
#include "TwoDM.h"
#include "ThreeDM.h"
#include "DMRG.h"
#include "Problem.h"
#include "Options.h"
#include "ConvergenceScheme.h"
#include "DMRGSCFindices.h"
#include "DMRGSCFunitary.h"
#include "DMRGSCFoptions.h"
#include "DMRGSCFwtilde.h"
#include "DMRGSCFmatrix.h"
#include "DMRGSCFintegrals.h"

namespace CheMPS2{
   class CASSCF{

      public:
      

         CASSCF( Hamiltonian * ham_in, int * docc, int * socc, int * nocc, int * ndmrg, int * nvirt, const string tmp_folder=CheMPS2::defaultTMPpath );
         
         virtual ~CASSCF();
         

         int get_num_irreps();
         

         double solve( const int Nelectrons, const int TwoS, const int Irrep, ConvergenceScheme * OptScheme, const int rootNum, DMRGSCFoptions * scf_options );


         double caspt2( const int Nelectrons, const int TwoS, const int Irrep, ConvergenceScheme * OptScheme, const int rootNum, DMRGSCFoptions * scf_options, const double IPEA, const double IMAG, const bool PSEUDOCANONICAL, const bool CHECKPOINT = false, const bool CUMULANT = false );

         /* \param filename File to delete */
         static void deleteStoredUnitary( const string filename=CheMPS2::DMRGSCF_unitary_storage_name ){ delete_file( filename ); }

         /* \param filename File to delete */
         static void deleteStoredDIIS( const string filename=CheMPS2::DMRGSCF_diis_storage_name ){ delete_file( filename ); }


         static void buildFmat( DMRGSCFmatrix * localFmat, const DMRGSCFmatrix * localTmat, const DMRGSCFmatrix * localJKocc, const DMRGSCFmatrix * localJKact, const DMRGSCFindices * localIdx, const DMRGSCFintegrals * theInts, double * local2DM, double * local1DM );


         static void buildWtilde( DMRGSCFwtilde * localwtilde, const DMRGSCFmatrix * localTmat, const DMRGSCFmatrix * localJKocc, const DMRGSCFmatrix * localJKact, const DMRGSCFindices * localIdx, const DMRGSCFintegrals * theInts, double * local2DM, double * local1DM );


         static void augmentedHessianNR( DMRGSCFmatrix * localFmat, DMRGSCFwtilde * localwtilde, const DMRGSCFindices * localIdx, const DMRGSCFunitary * localUmat, double * theupdate, double * updateNorm, double * gradNorm );


         static void copy2DMover( TwoDM * theDMRG2DM, const int LAS, double * two_dm );


         static void setDMRG1DM( const int num_elec, const int LAS, double * one_dm, double * two_dm );


         static void copy_active( double * origin, DMRGSCFmatrix * result, const DMRGSCFindices * idx, const bool one_rdm );


         static void copy_active( const DMRGSCFmatrix * origin, double * result, const DMRGSCFindices * idx );


         static void fillLocalizedOrbitalRotations( DMRGSCFunitary * umat, DMRGSCFindices * idx, double * eigenvecs );


         static void block_diagonalize( const char space, const DMRGSCFmatrix * Mat, DMRGSCFunitary * Umat, double * work1, double * work2, const DMRGSCFindices * idx, const bool invert, double * two_dm, double * three_dm, double * contract );


         static void construct_fock( DMRGSCFmatrix * Fock, const DMRGSCFmatrix * Tmat, const DMRGSCFmatrix * Qocc, const DMRGSCFmatrix * Qact, const DMRGSCFindices * idx );


         static double deviation_from_blockdiag( DMRGSCFmatrix * matrix, const DMRGSCFindices * idx );


         static void write_f4rdm_checkpoint( const string f4rdm_file, int * hamorb1, int * hamorb2, const int tot_dmrg_power6, double * contract );


         static bool read_f4rdm_checkpoint( const string f4rdm_file, int * hamorb1, int * hamorb2, const int tot_dmrg_power6, double * contract );


         static void fock_dot_4rdm( double * fockmx, CheMPS2::DMRG * dmrgsolver, CheMPS2::Hamiltonian * ham, int next_orb1, int next_orb2, double * work, double * result, const bool CHECKPOINT, const bool PSEUDOCANONICAL );

      private:

         // CASSCF tmp folder
         string tmp_folder;

         // Index convention handler
         DMRGSCFindices * iHandler;

         // Unitary matrix storage and manipulator
         DMRGSCFunitary * unitary;

         // Whether CheMPS2::CASSCF::solve has been successfully terminated
         bool successful_solve;

         // The original Hamiltonian
         double NUCL_ORIG;
         const TwoIndex  * TMAT_ORIG;
         const FourIndex * VMAT_ORIG;

         // Irreps controller
         Irreps SymmInfo;

         // The number of orbitals
         int L;

         // The number of irreps
         int num_irreps;

         // Number of DMRG orbitals
         int nOrbDMRG;

         // Space for the DMRG 1DM
         double * DMRG1DM;

         // Space for the DMRG 2DM
         double * DMRG2DM;

         // Helper function to check HF
         void checkHF( int * docc, int * socc );

         // Fill Econst and Tmat of HamDMRG
         void fillConstAndTmatDMRG( Hamiltonian * HamDMRG ) const;

         // Calculate the gradient, return function is the gradient 2-norm
         static double construct_gradient( DMRGSCFmatrix * Fmatrix, const DMRGSCFindices * idx, double * gradient );

         // Add hessian * origin to target
         static void add_hessian( DMRGSCFmatrix * Fmatrix, DMRGSCFwtilde * Wtilde, const DMRGSCFindices * idx, double * origin, double * target );

         // Construct the diagonal of the hessian
         static void diag_hessian( DMRGSCFmatrix * Fmatrix, const DMRGSCFwtilde * Wtilde, const DMRGSCFindices * idx, double * diagonal );

         // Apply augmented hessian
         static void DGEMM_WRAP( double prefactor, char transA, char transB, double * A, double * B, double * C, int m, int n, int k, int lda, int ldb, int ldc );
         static void DGEMV_WRAP( double prefactor, double * matrix, double * result, double * vector, int rowdim, int coldim, int ldmat, int incres, int incvec );
         static void augmented_hessian( DMRGSCFmatrix * Fmatrix, DMRGSCFwtilde * Wtilde, const DMRGSCFindices * idx, double * origin, double * target, double * gradient, const int linsize );

         // Rotate an active space object
         static void rotate_active_space_object( const int num_indices, double * object, double * work, double * rotation, const int LAS, const int NJUMP, const int NROTATE );

         // Fmat function as defined by Eq. (11) in the Siegbahn paper.
         DMRGSCFmatrix * theFmatrix;

         // The Coulomb and exchange interaction with the occupied and active electrons respectively
         DMRGSCFmatrix * theQmatOCC;
         DMRGSCFmatrix * theQmatACT;
         DMRGSCFmatrix * theQmatWORK;
         DMRGSCFmatrix * theTmatrix;
         void rotateOldToNew(DMRGSCFmatrix * myMatrix);
         void buildTmatrix();
         void constructCoulombAndExchangeMatrixInOrigIndices( DMRGSCFmatrix * density, DMRGSCFmatrix * result );
         void buildQmatOCC();
         void buildQmatACT();

         // Function to get coefficients of certain Slater determinants for Fe2. Important to figure out diatomic D(inf)h symmetries when calculating them in D2h symmetry.
         void coeff_fe2( DMRG * theDMRG );

         static void delete_file( const string filename );

   };
}

#endif