1 /* 2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry 3 Copyright (C) 2013-2018 Sebastian Wouters 4 5 This program is free software; you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation; either version 2 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License along 16 with this program; if not, write to the Free Software Foundation, Inc., 17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 18 */ 19 20 #ifndef CASPT2_CHEMPS2_H 21 #define CASPT2_CHEMPS2_H 22 23 #include "DMRGSCFindices.h" 24 #include "DMRGSCFintegrals.h" 25 #include "DMRGSCFmatrix.h" 26 27 #define CHEMPS2_CASPT2_A 0 28 #define CHEMPS2_CASPT2_B_SINGLET 1 29 #define CHEMPS2_CASPT2_B_TRIPLET 2 30 #define CHEMPS2_CASPT2_C 3 31 #define CHEMPS2_CASPT2_D 4 32 #define CHEMPS2_CASPT2_E_SINGLET 5 33 #define CHEMPS2_CASPT2_E_TRIPLET 6 34 #define CHEMPS2_CASPT2_F_SINGLET 7 35 #define CHEMPS2_CASPT2_F_TRIPLET 8 36 #define CHEMPS2_CASPT2_G_SINGLET 9 37 #define CHEMPS2_CASPT2_G_TRIPLET 10 38 #define CHEMPS2_CASPT2_H_SINGLET 11 39 #define CHEMPS2_CASPT2_H_TRIPLET 12 40 #define CHEMPS2_CASPT2_NUM_CASES 13 41 42 namespace CheMPS2{ 43 /** CASPT2 class. 44 \author Sebastian Wouters <sebastianwouters@gmail.com> 45 \date December 11, 2015 46 47 \section theo_caspt2 Information 48 49 The CASPT2 class contains the functions to perform second order multireference perturbation theory on top of a CASSCF wavefuntion [CASPT1, CASPT2]. CASPT2 has recently also been used with DMRG as active space solver: the active space 4-RDM contracted with the Fock operator, together with the 1-, 2- and 3-RDM are required thereto [CASPT3]. Alternatively, cumulant approximations can be used as well [CASPT4]. To mitigate problems with CASPT2 several modifications of the zeroth order Hamiltonian have been introduced: IPEA corrections [CASPT5], the g1 term [CASPT6], real level shifts [CASPT7] and imaginary level shifts [CASPT8]. 50 51 \section biblio_caspt2 References 52 53 [CASPT1] K. Andersson, P.-A. Malmqvist, B.O. Roos, A.J. Sadlej and K. Wolinski, Journal of Physical Chemistry 94, 5483-5488 (1990). http://dx.doi.org/10.1021/j100377a012 \n 54 [CASPT2] K. Andersson, P.‐A. Malmqvist and B.O. Roos, Journal of Chemical Physics 96, 1218-1226 (1992). http://dx.doi.org/10.1063/1.462209 \n 55 [CASPT3] Y. Kurashige and T. Yanai, Journal of Chemical Physics 135, 094104 (2011). http://dx.doi.org/10.1063/1.3629454 \n 56 [CASPT4] Y. Kurashige, J. Chalupsky, T.N. Lan and T. Yanai, Journal of Chemical Physics 141, 174111 (2014). http://dx.doi.org/10.1063/1.4900878 \n 57 [CASPT5] G. Ghigo, B.O. Roos and P.-A. Malmqvist, Chemical Physics Letters 396, 142-149 (2004). http://dx.doi.org/10.1016/j.cplett.2004.08.032 \n 58 [CASPT6] K. Andersson, Theoretica Chimica Acta 91, 31-46 (1995). http://dx.doi.org/10.1007/BF01113860 \n 59 [CASPT7] B.O. Roos and K. Andersson, Chemical Physics Letters 245, 215-223 (1995). http://dx.doi.org/10.1016/0009-2614(95)01010-7 \n 60 [CASPT8] N. Forsberg and P.-A. Malmqvist, Chemical Physics Letters 274, 196-204 (1997). http://dx.doi.org/10.1016/S0009-2614(97)00669-6 \n 61 */ 62 class CASPT2{ 63 64 public: 65 66 //! Constructor 67 /** \param idx DMRGSCFindices which contain the partitioning into occupied, active, and virtual orbitals per irrep 68 \param ints The two-electron integrals needed for CASSCF and CASPT2, in pseudocanonical orbitals 69 \param oei All one-electron integrals, in pseudocanonical orbitals 70 \param fock The fock matrix of CASPT2, in pseudocanonical orbitals 71 \param one_dm The spin-summed one-particle density matrix one_dm[i+L*j] = sum_sigma < a^+_i,sigma a_j,sigma > (with L the number DMRG orbitals), in pseudocanonical orbitals 72 \param two_dm The spin-summed two-particle density matrix two_dm[i+L*(j+L*(k+L*l))] = sum_sigma,tau < a^+_i,sigma a^+_j,tau a_l,tau a_k,sigma > (with L the number DMRG orbitals), in pseudocanonical orbitals 73 \param three_dm The spin-summed three-particle density matrix three_dm[i+L*(j+L*(k+L*(l+L*(m+L*n))))] = sum_z,tau,s < a^+_{i,z} a^+_{j,tau} a^+_{k,s} a_{n,s} a_{m,tau} a_{l,z} > (with L the number DMRG orbitals), in pseudocanonical orbitals 74 \param contract The spin-summed four-particle density matrix contracted with the fock operator contract[i+L*(j+L*(k+L*(p+L*(q+L*r))))] = sum_{t,sigma,tau,s} fock(t,t) < a^+_{i,sigma} a^+_{j,tau} a^+_{k,s} E_{tt} a_{r,s} a_{q,tau} a_{p,sigma} > (with L the number DMRG orbitals), in pseudocanonical orbitals 75 \param IPEA The CASPT2 IPEA shift from Ghigo, Roos and Malmqvist, Chemical Physics Letters 396, 142-149 (2004) */ 76 CASPT2(DMRGSCFindices * idx, DMRGSCFintegrals * ints, DMRGSCFmatrix * oei, DMRGSCFmatrix * fock, double * one_dm, double * two_dm, double * three_dm, double * contract, const double IPEA); 77 78 //! Destructor 79 virtual ~CASPT2(); 80 81 //! Solve for the CASPT2 energy (note that the IPEA shift has been set in the constructor) 82 /** \param imag_shift The CASPT2 imaginary shift from Forsberg and Malmqvist, Chemical Physics Letters 274, 196-204 (1997) 83 \param CONJUGATE_GRADIENT If true (false), the conjugate gradient (Davidson) algorithm is used to solve the CASPT2 equation 84 \return The CASPT2 variational correction energy */ 85 double solve( const double imag_shift, const bool CONJUGATE_GRADIENT = false ) const; 86 87 //! Return the vector length for the CASPT2 first order wavefunction (before diagonalization of the overlap matrix) 88 /** \param idx The number of core, active, and virtual orbitals per irrep 89 \return The vector length for the CASPT2 first order wavefunction (before diagonalization of the overlap matrix) */ 90 static long long vector_length( const DMRGSCFindices * idx ); 91 92 private: 93 94 // The number of occupied, active, and virtual orbitals per irrep (externally allocated and deleted) 95 const DMRGSCFindices * indices; 96 97 // The Fock matrix in pseudocanonical orbitals (externally allocated and deleted) 98 const DMRGSCFmatrix * fock; 99 100 // The active space 1-RDM (externally allocated and deleted) 101 double * one_rdm; 102 103 // The active space 2-RDM (externally allocated and deleted) 104 double * two_rdm; 105 106 // The active space 3-RDM (externally allocated and deleted) 107 double * three_rdm; 108 109 // The active space 4-RDM contracted with the Fock operator (externally allocated and deleted) 110 double * f_dot_4dm; 111 112 // The active space 3-RDM contracted with the Fock operator (allocated and deleted in this class) 113 double * f_dot_3dm; 114 115 // The active space 2-RDM contracted with the Fock operator (allocated and deleted in this class) 116 double * f_dot_2dm; 117 118 // The active space 1-RDM contracted with the Fock operator 119 double f_dot_1dm; 120 121 // The number of irreps 122 int num_irreps; 123 124 // Calculate the expectation value of the Fock operator 125 void create_f_dots(); 126 127 // Calculate the total vector length and the partitioning of the vector in blocks 128 int vector_helper(); 129 130 // Once make_S**() has been calles, these overlap matrices can be used to contruct the RHS of the linear problem 131 void construct_rhs( const DMRGSCFmatrix * oei, const DMRGSCFintegrals * integrals ); 132 133 // Fill result with the diagonal elements of the Fock operator 134 void diagonal( double * result ) const; 135 136 // Fill result with Fock operator times vector 137 void matvec( double * vector, double * result, double * diag_fock ) const; 138 static void matmat( char totrans, int rowdim, int coldim, int sumdim, double alpha, double * matrix, int ldaM, double * origin, int ldaO, double * target, int ldaT ); 139 140 // Helper functions for solve 141 void add_shift( double * vector, double * result, double * diag_fock, const double shift, const int * normalizations ) const; 142 double inproduct_vectors( double * first, double * second, const int * normalizations ) const; 143 void energy_per_sector( double * solution ) const; 144 145 // Variables for the partitioning of the vector in blocks 146 int * jump; 147 int * size_A; 148 int * size_C; 149 int * size_D; 150 int * size_E; 151 int * size_G; 152 int * size_B_singlet; 153 int * size_B_triplet; 154 int * size_F_singlet; 155 int * size_F_triplet; 156 157 // Functions for the partitioning of the vector in blocks 158 int get_maxsize() const; 159 static int jump_AC_active( const DMRGSCFindices * idx, const int irrep_t, const int irrep_u, const int irrep_v ); 160 static int jump_BF_active( const DMRGSCFindices * idx, const int irrep_t, const int irrep_u, const int ST ); 161 static int shift_D_nonactive( const DMRGSCFindices * idx, const int irrep_i, const int irrep_a ); 162 static int shift_B_nonactive( const DMRGSCFindices * idx, const int irrep_i, const int irrep_j, const int ST ); 163 static int shift_F_nonactive( const DMRGSCFindices * idx, const int irrep_a, const int irrep_b, const int ST ); 164 static int shift_E_nonactive( const DMRGSCFindices * idx, const int irrep_a, const int irrep_i, const int irrep_j, const int ST ); 165 static int shift_G_nonactive( const DMRGSCFindices * idx, const int irrep_i, const int irrep_a, const int irrep_b, const int ST ); 166 static int shift_H_nonactive( const DMRGSCFindices * idx, const int irrep_i, const int irrep_j, const int irrep_a, const int irrep_b, const int ST ); 167 168 // The RHS of the linear problem 169 double * vector_rhs; 170 171 // Variables for the overlap (only allocated during creation of the CASPT2 object) 172 double ** SAA; 173 double ** SCC; 174 double ** SDD; 175 double ** SEE; 176 double ** SGG; 177 double ** SBB_singlet; 178 double ** SBB_triplet; 179 double ** SFF_singlet; 180 double ** SFF_triplet; 181 182 // Variables for the diagonal part of the Fock operator 183 double ** FAA; 184 double ** FCC; 185 double ** FDD; 186 double ** FEE; 187 double ** FGG; 188 double ** FBB_singlet; 189 double ** FBB_triplet; 190 double ** FFF_singlet; 191 double ** FFF_triplet; 192 193 // Variables for the off-diagonal part of the Fock operator: Operator[ IL ][ IR ][ w ][ left + SIZE * right ] 194 double **** FAD; 195 double **** FCD; 196 double *** FEH; 197 double *** FGH; 198 double **** FAB_singlet; 199 double **** FAB_triplet; 200 double **** FCF_singlet; 201 double **** FCF_triplet; 202 double **** FBE_singlet; 203 double **** FBE_triplet; 204 double **** FFG_singlet; 205 double **** FFG_triplet; 206 double **** FDE_singlet; 207 double **** FDE_triplet; 208 double **** FDG_singlet; 209 double **** FDG_triplet; 210 211 // Fill overlap and Fock matrices 212 void make_AA_CC( const bool OVLP, const double IPEA ); 213 void make_DD( const bool OVLP, const double IPEA ); 214 void make_EE_GG( const bool OVLP, const double IPEA ); 215 void make_BB_FF_singlet( const bool OVLP, const double IPEA ); 216 void make_BB_FF_triplet( const bool OVLP, const double IPEA ); 217 218 void make_FAD_FCD(); 219 void make_FEH_FGH(); 220 void make_FAB_FCF_singlet(); 221 void make_FAB_FCF_triplet(); 222 void make_FBE_FFG_singlet(); 223 void make_FBE_FFG_triplet(); 224 void make_FDE_FDG(); 225 226 // Diagonalize the overlap matrices and adjust jump, vector_rhs, and FXX accordingly 227 void recreate(); 228 static int recreatehelper1( double * FOCK, double * OVLP, int SIZE, double * work, double * eigs, int lwork ); 229 static void recreatehelper2( double * LEFT, double * RIGHT, double ** matrix, double * work, int OLD_LEFT, int NEW_LEFT, int OLD_RIGHT, int NEW_RIGHT, const int number ); 230 static void recreatehelper3( double * OVLP, int OLDSIZE, int NEWSIZE, double * rhs_old, double * rhs_new, const int num_rhs ); 231 232 }; 233 } 234 235 #endif 236