CheMPS2

#include <CASSCF.h>
Public Member Functions  
CASSCF (Hamiltonian *ham_in, int *docc, int *socc, int *nocc, int *ndmrg, int *nvirt, const string tmp_folder=CheMPS2::defaultTMPpath)  
Constructor. More...  
virtual  ~CASSCF () 
Destructor.  
int  get_num_irreps () 
Get the number of irreps. More...  
double  solve (const int Nelectrons, const int TwoS, const int Irrep, ConvergenceScheme *OptScheme, const int rootNum, DMRGSCFoptions *scf_options) 
Do the CASSCF cycles with the augmented Hessian NewtonRaphson method. More...  
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) 
Calculate the caspt2 correction energy for a converged casscf wavefunction. More...  
Static Public Member Functions  
static void  deleteStoredUnitary (const string filename=CheMPS2::DMRGSCF_unitary_storage_name) 
CASSCF unitary rotation remove call.  
static void  deleteStoredDIIS (const string filename=CheMPS2::DMRGSCF_diis_storage_name) 
CASSCF DIIS vectors remove call.  
static void  buildFmat (DMRGSCFmatrix *localFmat, const DMRGSCFmatrix *localTmat, const DMRGSCFmatrix *localJKocc, const DMRGSCFmatrix *localJKact, const DMRGSCFindices *localIdx, const DMRGSCFintegrals *theInts, double *local2DM, double *local1DM) 
Build the Fmatrix (Eq. (11) in the Siegbahn paper [CAS3]) More...  
static void  buildWtilde (DMRGSCFwtilde *localwtilde, const DMRGSCFmatrix *localTmat, const DMRGSCFmatrix *localJKocc, const DMRGSCFmatrix *localJKact, const DMRGSCFindices *localIdx, const DMRGSCFintegrals *theInts, double *local2DM, double *local1DM) 
Build the Wtildematrix (Eq. (20b) in the Siegbahn paper [CAS3]) More...  
static void  augmentedHessianNR (DMRGSCFmatrix *localFmat, DMRGSCFwtilde *localwtilde, const DMRGSCFindices *localIdx, const DMRGSCFunitary *localUmat, double *theupdate, double *updateNorm, double *gradNorm) 
Calculate the augmented Hessian NewtonRaphson update for the orthogonal orbital rotation matrix. More...  
static void  copy2DMover (TwoDM *theDMRG2DM, const int LAS, double *two_dm) 
Copy over the DMRG 2RDM. More...  
static void  setDMRG1DM (const int num_elec, const int LAS, double *one_dm, double *two_dm) 
Construct the 1RDM from the 2RDM. More...  
static void  copy_active (double *origin, DMRGSCFmatrix *result, const DMRGSCFindices *idx, const bool one_rdm) 
Copy a oneorbital quantity from array format to DMRGSCFmatrix format. More...  
static void  copy_active (const DMRGSCFmatrix *origin, double *result, const DMRGSCFindices *idx) 
Copy a oneorbital quantity from DMRGSCFmatrix format to array format. More...  
static void  fillLocalizedOrbitalRotations (DMRGSCFunitary *umat, DMRGSCFindices *idx, double *eigenvecs) 
From an EdmistonRuedenberg active space rotation, fetch the eigenvectors and store them in eigenvecs. More...  
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) 
Blockdiagonalize Mat. More...  
static void  construct_fock (DMRGSCFmatrix *Fock, const DMRGSCFmatrix *Tmat, const DMRGSCFmatrix *Qocc, const DMRGSCFmatrix *Qact, const DMRGSCFindices *idx) 
Construct the Fock matrix. More...  
static double  deviation_from_blockdiag (DMRGSCFmatrix *matrix, const DMRGSCFindices *idx) 
Return the RMS deviation from blockdiagonal. More...  
static void  write_f4rdm_checkpoint (const string f4rdm_file, int *hamorb1, int *hamorb2, const int tot_dmrg_power6, double *contract) 
Write the checkpoint file for the contraction of the generalized Fock operator with the 4RDM to disk. More...  
static bool  read_f4rdm_checkpoint (const string f4rdm_file, int *hamorb1, int *hamorb2, const int tot_dmrg_power6, double *contract) 
Read the checkpoint file for the contraction of the generalized Fock operator with the 4RDM from disk. More...  
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) 
Build the contraction of the fock matrix with the 4RDM. More...  
CASSCF class.
In methods which use a FCI solver, this solver can be replaced by DMRG. DMRG allows for an efficient extraction of the 2RDM [CAS1, CAS2]. The 2RDM of the active space is required in the complete active space selfconsistent field (CASSCF) method to compute the gradient and the Hessian with respect to orbital rotations [CAS3]. It is therefore natural to introduce a CASSCF variant with DMRG as active space solver, called DMRGSCF [CAS2, CAS4, CAS5], which allows to treat static correlation in large active spaces. In CheMPS2, I have implemented the augmented Hessian NewtonRaphson DMRGSCF method, with exact Hessian [CAS5, CAS6]. It can be called with the function CheMPS2::CASSCF::doCASSCFnewtonraphson.
The calculation of the orbital gradient and Hessian for DMRGSCF is based on [CAS3]. The basic idea is to express the energy with the unitary group generators:
The variables only connect orbitals with the same irrep ( ). Assuming that DMRG is exact, in addition only connects orbitals when they belong to different occupation blocks: occupied, active, virtual. With some algebra, the derivatives can be rewritten. Realvalued symmetric oneelectron integrals and realvalued eightfold permutation symmetric twoelectron integrals are assumed (chemical notation for the latter).
In the calculation of , the indices can only be occupied or active due to their appearance in the density matrices, and the only index which can be virtual is hence . Moreover, due to the irrep symmetry of the integrals and density matrices, is diagonal in the irreps: . Alternatively, this can be understood by the fact that only connects orbitals with the same irrep.
In the calculation of , the indices can only be occupied or active due to their appearance in the density matrices, and the only indices which can be virtual are hence . Together with the remark for , this can save time for the twoelectron integral rotation. Moreover, as only connects orbitals with the same irrep, and in .
By rewriting the density matrices, the calculation of and can be simplified. In the following, occ and act denote the doubly occupied and active orbital spaces, respectively.
Define the following symmetric charge (Coulomb + exchange) matrices:
They can be calculated efficiently by (1) rotating the occupied and active density matrices from the current basis to the original basis, (2) contracting the rotated density matrices with the twoelectron integrals in the original basis, and (3) rotating these contractions to the current basis. The constant part and the oneelectron integrals of the active space Hamiltonian are:
The calculation of boils down to:
And the calculation of (remember that and ):
The CASSCF energy is a function of . Up to second order, the energy is given by
The vector is the gradient and the matrix the Hessian for orbital rotations [CAS3]. They have been described in the previous section. The minimum of is found at . The variables parametrize an additional orbital rotation , with a realvalued skewsymmetric matrix. The additional orbital rotation transforms the current orbitals to the new orbitals
This updating scheme is called the NewtonRaphson method [CAS3]. If the Hessian is positive definite, these updates are stable. For saddle points in the energy landscape, the Hessian has negative eigenvalues, and these updates can be unstable. It is therefore better to use the augmented Hessian NewtonRaphson method [CAS7]:
The eigenvector with smallest algebraic eigenvalue determines a stable update , as is well explained in Ref. [CAS7].
As a final remark in this section, I would like to say that orbitals have gauge freedom. One can always multiply them with a phase factor. It is therefore possible to choose the orbital gauges so that all are always special orthogonal: .
When the update norm is small enough, the convergence can be accelerated by the direct inversion of the iterative subspace (DIIS) [CAS5, CAS8, CAS9, CAS10]. For a given set of orbitals , the update is calculated with the augmented Hessian NewtonRaphson method. This update defines the next set of orbitals:
In DIIS, the error vector and the state vector are added to a list. The error
is minimized under the constraint . is the size of the list memory, i.e. the number of retained vectors. The minimization of the error can be performed with Lagrangian calculus:
where is the Lagrangian multiplier and
The new state vector is then defined as
The new state vector is calculated by the function CheMPS2::DIIS::calculateParam. The current orbitals are then set to .
[CAS1] D. Zgid and M. Nooijen, Journal of Chemical Physics 128, 144115 (2008). http://dx.doi.org/10.1063/1.2883980
[CAS2] D. Ghosh, J. Hachmann, T. Yanai and G.K.L. Chan, Journal of Chemical Physics 128, 144117 (2008). http://dx.doi.org/10.1063/1.2883976
[CAS3] P.E.M. Siegbahn, J. Almlof, A. Heiberg and B.O. Roos, Journal of Chemical Physics 74, 23842396 (1981). http://dx.doi.org/10.1063/1.441359
[CAS4] D. Zgid and M. Nooijen, Journal of Chemical Physics 128, 144116 (2008). http://dx.doi.org/10.1063/1.2883981
[CAS5] T. Yanai, Y. Kurashige, D. Ghosh and G.K.L. Chan, International Journal of Quantum Chemistry 109, 21782190 (2009). http://dx.doi.org/10.1002/qua.22099
[CAS6] S. Wouters, W. Poelmans, P.W. Ayers and D. Van Neck, Computer Physics Communications 185, 15011514 (2014). http://dx.doi.org/10.1016/j.cpc.2014.01.019
[CAS7] A. Banerjee, N. Adams, J. Simons and R. Shepard, Journal of Physical Chemistry 89, 5257 (1985). http://dx.doi.org/10.1021/j100247a015
[CAS8] P. Pulay, Chemical Physics Letters 73, 393398 (1980). http://dx.doi.org/10.1016/00092614(80)803964
[CAS9] C.D. Sherrill, Programming DIIS, http://vergil.chemistry.gatech.edu/notes/diis/node3.html (2000).
[CAS10] T. Rohwedder and R. Schneider, Journal of Mathematical Chemistry 49, 18891914 (2011). http://dx.doi.org/10.1007/s109100119863y
CheMPS2::CASSCF::CASSCF  (  Hamiltonian *  ham_in, 
int *  docc,  
int *  socc,  
int *  nocc,  
int *  ndmrg,  
int *  nvirt,  
const string  tmp_folder = CheMPS2::defaultTMPpath 

) 
Constructor.
ham_in  Hamiltonian containing the matrix elements of the Hamiltonian for which a CASSCF calculation is desired 
docc  Array containing the number of doubly occupied HF orbitals per irrep 
socc  Array containing the number of singly occupied HF orbitals per irrep 
nocc  Array containing the number of doubly occupied (inactive) orbitals per irrep 
ndmrg  Array containing the number of active orbitals per irrep 
nvirt  Array containing the number of virtual (secondary) orbitals per irrep 
tmp_folder  Temporary work folder for the DMRG renormalized operators and the ERI rotations 
Definition at line 39 of file CASSCF.cpp.

static 
Calculate the augmented Hessian NewtonRaphson update for the orthogonal orbital rotation matrix.
localFmat  Matrix which contains the Fock operator (Eq. (11) in the Siegbahn paper [CAS3]) 
localwtilde  Object which contains the second order derivative of the energy with respect to the unitary (Eq. (20b) in the Siegbahn paper [CAS3]) 
localIdx  Orbital index bookkeeper for the CASSCF calculations 
localUmat  The unitary matrix for CASSCF calculations (in this function it is used to fetch the orbital ordering convention of the skewsymmetric parametrization) 
theupdate  Where the augmented Hessian NewtonRaphson update will be stored 
updateNorm  Pointer to one double to store the update norm 
gradNorm  Pointer to one double to store the gradient norm 
Definition at line 316 of file CASSCFnewtonraphson.cpp.

static 
Blockdiagonalize Mat.
space  Can be 'O', 'A', or 'V' and denotes which block of Mat should be considered 
Mat  Matrix to blockdiagonalize 
Umat  The unitary rotation will be updated so that Mat is blockdiagonal in the orbitals 'space' 
work1  Workspace 
work2  Workspace 
idx  Object which handles the index conventions for CASSCF 
invert  If true, the eigenvectors are sorted from large to small instead of the other way around 
two_dm  If not NULL, this 4index array will be rotated to the new eigenvecs if space == 'A' 
three_dm  If not NULL, this 6index array will be rotated to the new eigenvecs if space == 'A' 
contract  If not NULL, this 6index array will be rotated to the new eigenvecs if space == 'A' 
Definition at line 433 of file CASSCF.cpp.

static 
Build the Fmatrix (Eq. (11) in the Siegbahn paper [CAS3])
localFmat  Matrix where the result should be stored 
localTmat  Matrix which contains the oneelectron integrals 
localJKocc  Matrix which contains the Coulomb and exchange interaction due to the frozen core orbitals 
localJKact  Matrix which contains the Coulomb and exchange interaction due to the active space 
localIdx  Orbital index bookkeeper for the CASSCF calculations 
theInts  The rotated twoelectron integrals (at most 2 virtual indices) 
local2DM  The DMRG 2RDM 
local1DM  The DMRG 1RDM 
Definition at line 823 of file CASSCFnewtonraphson.cpp.

static 
Build the Wtildematrix (Eq. (20b) in the Siegbahn paper [CAS3])
localwtilde  Where the result should be stored 
localTmat  Matrix which contains the oneelectron integrals 
localJKocc  Matrix which contains the Coulomb and exchange interaction due to the frozen core orbitals 
localJKact  Matrix which contains the Coulomb and exchange interaction due to the active space 
localIdx  Orbital index bookkeeper for the CASSCF calculations 
theInts  The rotated twoelectron integrals (at most 2 virtual indices) 
local2DM  The DMRG 2RDM 
local1DM  The DMRG 1RDM 
Definition at line 611 of file CASSCFnewtonraphson.cpp.
double CheMPS2::CASSCF::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 

) 
Calculate the caspt2 correction energy for a converged casscf wavefunction.
Nelectrons  Total number of electrons in the system: occupied HF orbitals + active space 
TwoS  Twice the targeted spin 
Irrep  Desired wavefunction irrep 
OptScheme  The optimization scheme to run the inner DMRG loop. If NULL: use FCI instead of DMRG. 
rootNum  Denotes the targeted state in statespecific CASSCF; 1 means ground state, 2 first excited state etc. 
scf_options  Contains the DMRGSCF options 
IPEA  The CASPT2 IPEA shift from Ghigo, Roos and Malmqvist, Chemical Physics Letters 396, 142149 (2004) 
IMAG  The CASPT2 imaginary shift from Forsberg and Malmqvist, Chemical Physics Letters 274, 196204 (1997) 
PSEUDOCANONICAL  If true, use the exact DMRG 4RDM in the pseudocanonical basis. If false, use the cumulant approximated DMRG 4RDM in the unrotated basis. 
CHECKPOINT  If true, write checkpoints to disk and read them back in again in order to perform the contraction of the generalized Fock operator with the 4RDM in multiple runs. 
CUMULANT  If true, a cumulant approximation is used for the 4RDM and CHECKPOINT is overwritten to false. If false, the full 4RDM is used. 
Definition at line 189 of file CASSCFpt2.cpp.

static 
Construct the Fock matrix.
Fock  Matrix to store the Fock operator in 
Tmat  Matrix with the oneelectron integrals 
Qocc  Matrix with the Coulomb and exchange contributions of the occupied (inactive) orbitals 
Qact  Matrix with the Coulomb and exchange contributions of the active space orbitals 
idx  Object which handles the index conventions for CASSCF 
Definition at line 395 of file CASSCFpt2.cpp.

static 

static 
Copy a oneorbital quantity from array format to DMRGSCFmatrix format.
origin  Array to copy 
result  DMRGSCFmatrix to store the copy 
idx  Object which handles the index conventions for CASSCF 
one_rdm  If true, the occupied orbitals get occupation 2 
Definition at line 375 of file CASSCF.cpp.

static 
Copy a oneorbital quantity from DMRGSCFmatrix format to array format.
origin  DMRGSCFmatrix to copy 
result  Array to store the copy 
idx  Object which handles the index conventions for CASSCF 
Definition at line 407 of file CASSCF.cpp.

static 
Return the RMS deviation from blockdiagonal.
matrix  Matrix to be assessed 
idx  Object which handles the index conventions for CASSCF 
Definition at line 411 of file CASSCFpt2.cpp.

static 
From an EdmistonRuedenberg active space rotation, fetch the eigenvectors and store them in eigenvecs.
umat  The EdmistonRuedenberg active space rotation 
idx  Object which handles the index conventions for CASSCF 
eigenvecs  Where the eigenvectors are stored 
Definition at line 160 of file CASSCF.cpp.

static 
Build the contraction of the fock matrix with the 4RDM.
fockmx  Array of size ham>getL() x ham>getL() containing the Fock matrix elements 
dmrgsolver  DMRG object which is solved, and for which the 2RDM and 3RDM have been calculated as well 
ham  Active space Hamiltonian, which is needed for the size of the active space and the orbital irreps 
next_orb1  The next first orbital for which the contraction should be continued 
next_orb2  The next second orbital for which the contraction should be continued 
work  Work array of size ham>getL()**6 
result  On entry, contains the partial contraction corresponding to (next_orb1, next_orb2). On exit, contains the full contraction. 
CHECKPOINT  Whether or not the standard CheMPS2 F.4RDM checkpoint should be created/updated to continue the contraction at later times 
PSEUDOCANONICAL  Whether or not pseudocanonical orbitals are used in the active space 
Definition at line 140 of file CASSCFpt2.cpp.
int CheMPS2::CASSCF::get_num_irreps  (  ) 

static 
Read the checkpoint file for the contraction of the generalized Fock operator with the 4RDM from disk.
f4rdm_file  The filename 
hamorb1  The next hamiltonian orbital 1 
hamorb2  The next hamiltonian orbital 2 
tot_dmrg_power6  The size of the array contract 
contract  The current partial contraction 
Definition at line 87 of file CASSCFpt2.cpp.

static 
Construct the 1RDM from the 2RDM.
num_elec  The number of DMRG active space electrons 
LAS  The total number of DMRG orbitals 
one_dm  The CASSCF 1RDM 
two_dm  The CASSCF 2RDM 
Definition at line 134 of file CASSCF.cpp.
double CheMPS2::CASSCF::solve  (  const int  Nelectrons, 
const int  TwoS,  
const int  Irrep,  
ConvergenceScheme *  OptScheme,  
const int  rootNum,  
DMRGSCFoptions *  scf_options  
) 
Do the CASSCF cycles with the augmented Hessian NewtonRaphson method.
Nelectrons  Total number of electrons in the system: occupied HF orbitals + active space 
TwoS  Twice the targeted spin 
Irrep  Desired wavefunction irrep 
OptScheme  The optimization scheme to run the inner DMRG loop. If NULL: use FCI instead of DMRG. 
rootNum  Denotes the targeted state in statespecific CASSCF; 1 means ground state, 2 first excited state etc. 
scf_options  Contains the DMRGSCF options 
Definition at line 56 of file CASSCFnewtonraphson.cpp.

static 
Write the checkpoint file for the contraction of the generalized Fock operator with the 4RDM to disk.
f4rdm_file  The filename 
hamorb1  The next hamiltonian orbital 1 
hamorb2  The next hamiltonian orbital 2 
tot_dmrg_power6  The size of the array contract 
contract  The current partial contraction 
Definition at line 44 of file CASSCFpt2.cpp.