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 FCI_CHEMPS2_H 21 #define FCI_CHEMPS2_H 22 23 #include "Hamiltonian.h" 24 25 namespace CheMPS2{ 26 /** FCI class. 27 \author Sebastian Wouters <sebastianwouters@gmail.com> 28 \date November 6, 2014 29 30 The FCI class performs the full configuration interaction ground state calculation in a given particle number, irrep, and spin projection sector of a given Hamiltonian. It also contains the functionality to calculate Green's functions. 31 */ 32 class FCI{ 33 34 public: 35 36 //! Constructor 37 /** \param Ham The Hamiltonian matrix elements 38 \param Nel_up The number of up (alpha) electrons 39 \param Nel_down The number of down (beta) electrons 40 \param TargetIrrep The targeted point group irrep 41 \param maxMemWorkMB Maximum workspace size in MB to be used for matrix vector product (this does not include the FCI vectors as stored for example in GSDavidson!!) 42 \param FCIverbose The FCI verbose level: 0 print nothing, 1 print start and solution, 2 print everything */ 43 FCI(CheMPS2::Hamiltonian * Ham, const unsigned int Nel_up, const unsigned int Nel_down, const int TargetIrrep, const double maxMemWorkMB=100.0, const int FCIverbose=2); 44 45 //! Destructor 46 virtual ~FCI(); 47 48 //==========> Getters of basic information: all these variables are set by the constructor 49 50 //! Getter for the number of orbitals 51 /** \return The number of orbitals */ getL()52 unsigned int getL() const{ return L; } 53 54 //! Getter for the number of up or alpha electrons 55 /** \return The number of up or alpha electrons */ getNel_up()56 unsigned int getNel_up() const{ return Nel_up; } 57 58 //! Getter for the number of down or beta electrons 59 /** \return The number of down or beta electrons */ getNel_down()60 unsigned int getNel_down() const{ return Nel_down; } 61 62 //! Getter for the number of variables in the vector " E_ij | FCI vector > " ; where irrep_center = I_i x I_j 63 /** \param irrep_center The single electron excitation irrep I_i x I_j 64 \return The number of variables in the corresponding vector */ getVecLength(const int irrep_center)65 unsigned int getVecLength(const int irrep_center) const{ return irrep_center_jumps[ irrep_center ][ num_irreps ]; } 66 67 //! Get the target irrep 68 /** \return The target irrep */ getTargetIrrep()69 int getTargetIrrep() const{ return TargetIrrep; } 70 71 //! Get an orbital irrep 72 /** \param orb The orbital index 73 \return The irrep of orbital orb */ getOrb2Irrep(const int orb)74 int getOrb2Irrep(const int orb) const{ return orb2irrep[ orb ]; } 75 76 //! Function which returns the electron repulsion integral (in chemists notation: integral dr1 dr2 orb1(r1) * orb2(r1) * orb3(r2) * orb4(r2) / |r1-r2|) 77 /** \param orb1 First orbital index (electron at position r1) 78 \param orb2 Second orbital index (electron at position r1) 79 \param orb3 Third orbital index (electron at position r2) 80 \param orb4 Fourth orbital index (electron at position r2) 81 \return The desired electron repulsion integral */ getERI(const int orb1,const int orb2,const int orb3,const int orb4)82 double getERI(const int orb1, const int orb2, const int orb3, const int orb4) const{ return ERI[ orb1 + L * ( orb2 + L * ( orb3 + L * orb4 ) ) ]; } 83 84 //! Function which returns the AUGMENTED one-body matrix elements ( G_{ij} = T_{ij} - 0.5 * sum_k ERI_{ikkj} ; see Chem. Phys. Lett. 111 (4-5), 315-321 (1984) ) 85 /** \param orb1 First orbital index 86 \param orb2 Second orbital index 87 \return The desired AUGMENTED one-body matrix element */ getGmat(const int orb1,const int orb2)88 double getGmat(const int orb1, const int orb2) const{ return Gmat[ orb1 + L * orb2 ]; } 89 90 //! Function which returns the nuclear repulsion energy 91 /** \return The nuclear repulsion energy */ getEconst()92 double getEconst() const{ return Econstant; } 93 94 //==========> The core routines for users 95 96 //! Calculates the FCI ground state with Davidson's algorithm 97 /** \param inoutput If inoutput!=NULL, vector with getVecLength(0) variables which contains the initial guess at the start, and on exit the solution of the FCI calculation 98 \param DVDSN_NUM_VEC The maximum number of vectors to use in Davidson's algorithm; adjustable in case memory becomes an issue 99 \return The ground state energy */ 100 double GSDavidson(double * inoutput=NULL, const int DVDSN_NUM_VEC=CheMPS2::DAVIDSON_NUM_VEC) const; 101 102 //! Return the global counter of the Slater determinant with the lowest energy 103 /** \return The global counter of the Slater determinant with the lowest energy */ 104 unsigned int LowestEnergyDeterminant() const; 105 106 //! Construct the (spin-summed) 2-RDM of a FCI vector: Gamma^2(i,j,k,l) = sum_sigma,tau < a^+_i,sigma a^+_j,tau a_l,tau a_k,sigma > = TwoRDM[ i + L * ( j + L * ( k + L * l ) ) ] 107 /** \param vector The FCI vector of length getVecLength(0) 108 \param TwoRDM To store the 2-RDM; needs to be of size getL()^4; point group symmetry shows in 2-RDM elements being zero 109 \return The energy of the given FCI vector, calculated by contraction of the 2-RDM with Gmat and ERI */ 110 double Fill2RDM(double * vector, double * TwoRDM) const; 111 112 //! Construct the (spin-summed) 3-RDM of a FCI vector: Gamma^3(i,j,k,l,m,n) = sum_sigma,tau,s < a^+_{i,sigma} a^+_{j,tau} a^+_{k,s} a_{n,s} a_{m,tau} a_{l,sigma} > = ThreeRDM[ i + L * ( j + L * ( k + L * ( l + L * ( m + L * n ) ) ) ) ] 113 /** \param vector The FCI vector of length getVecLength(0) 114 \param ThreeRDM To store the 3-RDM; needs to be of size getL()^6; point group symmetry shows in 3-RDM elements being zero */ 115 void Fill3RDM(double * vector, double * ThreeRDM) const; 116 117 //! Construct the (spin-summed) 4-RDM of a FCI vector: Gamma^4(i,j,k,l,p,q,r,t) = sum_sigma,tau,s,z < a^+_{i,sigma} a^+_{j,tau} a^+_{k,s} a^+_{l,z} a_{t,z} a_{r,s} a_{q,tau} a_{p,sigma} > = FourRDM[ i + L * ( j + L * ( k + L * ( l + L * ( p + L * ( q + L * ( r + L * t ) ) ) ) ) ) ] 118 /** \param vector The FCI vector of length getVecLength(0) 119 \param FourRDM To store the 4-RDM; needs to be of size getL()^8; point group symmetry shows in 4-RDM elements being zero */ 120 void Fill4RDM(double * vector, double * FourRDM) const; 121 122 //! Construct the (spin-summed) contraction of the 4-RDM with the Fock operator: output(i,j,k,p,q,r) = sum_{l,t} Fock(l,t) * Gamma^4(i,j,k,l,p,q,r,t) 123 /** \param vector The FCI vector of length getVecLength(0) 124 \param ThreeRDM The spin-summed 3-RDM as calculated by Fill3RDM 125 \param Fock The symmetric Fock operator Fock(i,j) = Fock[ i + L * j ] = Fock[ j + L * i ] 126 \param output To store the contraction output(i,j,k,p,q,r) = output[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r ) ) ) ) ]; needs to be of size getL()^6; point group symmetry shows in elements being zero; has 12-fold permutation symmetry just like 3-RDM */ 127 void Fock4RDM(double * vector, double * ThreeRDM, double * Fock, double * output) const; 128 129 //! Construct part of the 4-RDM: output(i,j,k,p,q,r) = Gamma^4(i,j,k,z,p,q,r,z) 130 /** \param vector The FCI vector of length getVecLength(0) 131 \param three_rdm The spin-summed 3-RDM as calculated by Fill3RDM 132 \param orbz The orbital z which is fixed in Gamma^4(i,j,k,z,p,q,r,z) 133 \param output To store part of the 4-RDM output(i,j,k,p,q,r) = output[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r ) ) ) ) ]; needs to be of size getL()^6; point group symmetry shows in elements being zero; has 12-fold permutation symmetry just like 3-RDM */ 134 void Diag4RDM( double * vector, double * three_rdm, const unsigned int orbz, double * output ) const; 135 136 //! Measure S(S+1) (spin squared) 137 /** \param vector The FCI vector of length getVecLength(0) 138 \return Measured value of S(S+1) */ 139 double CalcSpinSquared(double * vector) const; 140 141 //! Fill a vector with random numbers in the interval [-1,1[; used when output for GSDavidson is desired but no specific input can be given 142 /** \param vecLength The length of the vector; when used for GSDavidson it should be getVecLength(0) 143 \param vec The vector to fill with random numbers */ 144 static void FillRandom(const unsigned int vecLength, double * vec); 145 146 //! Set the entries of a vector to zero 147 /** \param vecLength The vector length 148 \param vec The vector which has to be set to zero */ 149 static void ClearVector(const unsigned int vecLength, double * vec); 150 151 //==========> Green's functions functionality 152 153 //! Calculate the retarded Green's function (= addition + removal amplitude) 154 /** \param omega The frequency value 155 \param eta The regularization parameter (... + I*eta in the denominator) 156 \param orb_alpha The first orbital index 157 \param orb_beta The second orbital index 158 \param isUp If true, the spin projection value of the second quantized operators is up, otherwise it will be down 159 \param GSenergy The ground state energy returned by GSDavidson 160 \param GSvector The ground state vector as calculated by GSDavidson 161 \param Ham The Hamiltonian, which contains the matrix elements 162 \param RePartGF On exit RePartGF[0] contains the real part of the retarded Green's function 163 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the retarded Green's function */ 164 void RetardedGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF) const; 165 166 //! Calculate the addition part of the retarded Green's function: <GSvector| a_{alpha, spin(isUp)} [ omega - Ham + GSenergy + I*eta ]^{-1} a^+_{beta, spin(isUp)} |GSvector> 167 /** \param omega The frequency value 168 \param eta The regularization parameter 169 \param orb_alpha The first orbital index 170 \param orb_beta The second orbital index 171 \param isUp If true, the spin projection value of the second quantized operators is up, otherwise it will be down 172 \param GSenergy The ground state energy returned by GSDavidson 173 \param GSvector The ground state vector as calculated by GSDavidson 174 \param Ham The Hamiltonian, which contains the matrix elements 175 \param RePartGF On exit RePartGF[0] contains the real part of the addition part of the retarded Green's function 176 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the addition part of the retarded Green's function 177 \param TwoRDMreal If not NULL, on exit the 2-RDM of Re{ [ omega - Ham + GSenergy + I*eta ]^{-1} a^+_{beta, spin(isUp)} | GSvector > } 178 \param TwoRDMimag If not NULL, on exit the 2-RDM of Im{ [ omega - Ham + GSenergy + I*eta ]^{-1} a^+_{beta, spin(isUp)} | GSvector > } 179 \param TwoRDMadd If not NULL, on exit the 2-RDM of a^+_{beta, spin(isUp)} | GSvector > */ 180 void RetardedGF_addition(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMadd=NULL) const; 181 182 //! Calculate the addition Green's function: GF[i+numLeft*j] = <GSvector| a_{orbsLeft[i], spin} [ alpha + beta * Ham + I*eta ]^{-1} a^+_{orbsRight[j], spin} |GSvector> 183 /** \param alpha Constant parameter in the resolvent 184 \param beta Prefector of the Hamiltonian in the resolvent 185 \param eta The regularization parameter 186 \param orbsLeft The left orbital indices 187 \param numLeft The number of left orbital indices 188 \param orbsRight The right orbital indices 189 \param numRight The number of right orbital indices 190 \param isUp If true, the spin projection value of the second quantized operators is up, otherwise it will be down 191 \param GSvector The ground state vector as calculated by GSDavidson 192 \param Ham The Hamiltonian, which contains the matrix elements 193 \param RePartsGF On exit RePartsGF[i+numLeft*j] contains the real part of the addition Green's function 194 \param ImPartsGF On exit ImPartsGF[i+numLeft*j] contains the imaginary part of the addition Green's function 195 \param TwoRDMreal If not NULL, TwoRDMreal[j] contains on exit the 2-RDM of Re{ [ alpha + beta * Ham + I*eta ]^{-1} a^+_{orbsRight[j], spin} | GSvector > } 196 \param TwoRDMimag If not NULL, TwoRDMimag[j] contains on exit the 2-RDM of Im{ [ alpha + beta * Ham + I*eta ]^{-1} a^+_{orbsRight[j], spin} | GSvector > } 197 \param TwoRDMadd If not NULL, TwoRDMadd[j] contains on exit the 2-RDM of a^+_{orbsRight[j], spin} | GSvector > */ 198 void GFmatrix_addition(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal=NULL, double ** TwoRDMimag=NULL, double ** TwoRDMadd=NULL) const; 199 200 //! Calculate the removal part of the retarded Green's function: <GSvector| a^+_{beta, spin(isUp)} [ omega + Ham - GSenergy + I*eta ]^{-1} a_{alpha, spin(isUp)} |GSvector> 201 /** \param omega The frequency value 202 \param eta The regularization parameter 203 \param orb_alpha The first orbital index 204 \param orb_beta The second orbital index 205 \param isUp If true, the spin projection value of the second quantized operators is up, otherwise it will be down 206 \param GSenergy The ground state energy returned by GSDavidson 207 \param GSvector The ground state vector as calculated by GSDavidson 208 \param Ham The Hamiltonian, which contains the matrix elements 209 \param RePartGF On exit RePartGF[0] contains the real part of the removal part of the retarded Green's function 210 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the removal part of the retarded Green's function 211 \param TwoRDMreal If not NULL, on exit the 2-RDM of Re{ [ omega + Ham - GSenergy + I*eta ]^{-1} a_{alpha, spin(isUp)} | GSvector > } 212 \param TwoRDMimag If not NULL, on exit the 2-RDM of Im{ [ omega + Ham - GSenergy + I*eta ]^{-1} a_{alpha, spin(isUp)} | GSvector > } 213 \param TwoRDMrem If not NULL, on exit the 2-RDM of a_{alpha, spin(isUp)} | GSvector > */ 214 void RetardedGF_removal(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMrem=NULL) const; 215 216 //! Calculate the removal Green's function: GF[i+numLeft*j] = <GSvector| a^+_{orbsLeft[i], spin} [ alpha + beta * Ham + I*eta ]^{-1} a_{orbsRight[j], spin} |GSvector> 217 /** \param alpha Constant parameter in the resolvent 218 \param beta Prefector of the Hamiltonian in the resolvent 219 \param eta The regularization parameter 220 \param orbsLeft The left orbital indices 221 \param numLeft The number of left orbital indices 222 \param orbsRight The right orbital indices 223 \param numRight The number of right orbital indices 224 \param isUp If true, the spin projection value of the second quantized operators is up, otherwise it will be down 225 \param GSvector The ground state vector as calculated by GSDavidson 226 \param Ham The Hamiltonian, which contains the matrix elements 227 \param RePartsGF On exit RePartsGF[i+numLeft*j] contains the real part of the removal Green's function GF[i+numLeft*j] 228 \param ImPartsGF On exit ImPartsGF[i+numLeft*j] contains the imaginary part of the removal Green's function GF[i+numLeft*j] 229 \param TwoRDMreal If not NULL, TwoRDMreal[j] contains on exit the 2-RDM of Re{ [ alpha + beta * Ham + I*eta ]^{-1} a_{orbsRight[j], spin} | GSvector > } 230 \param TwoRDMimag If not NULL, TwoRDMimag[j] contains on exit the 2-RDM of Im{ [ alpha + beta * Ham + I*eta ]^{-1} a_{orbsRight[j], spin} | GSvector > } 231 \param TwoRDMrem If not NULL, TwoRDMrem[j] contains on exit the 2-RDM of a_{orbsRight[j], spin} | GSvector > */ 232 void GFmatrix_removal(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal=NULL, double ** TwoRDMimag=NULL, double ** TwoRDMrem=NULL) const; 233 234 //! Calculate the density response Green's function (= forward - backward propagating part) 235 /** \param omega The frequency value 236 \param eta The regularization parameter (... + I*eta in the denominator) 237 \param orb_alpha The first orbital index 238 \param orb_beta The second orbital index 239 \param GSenergy The ground state energy returned by GSDavidson 240 \param GSvector The ground state vector as calculated by GSDavidson 241 \param RePartGF On exit RePartGF[0] contains the real part of the density response Green's function 242 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the density response Green's function */ 243 void DensityResponseGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF) const; 244 245 //! Calculate the forward propagating part of the density response Green's function: <GSvector| ( n_alpha - <GSvector| n_alpha |GSvector> ) [ omega - Ham + GSenergy + I*eta ]^{-1} ( n_beta - <GSvector| n_beta |GSvector> ) |GSvector> 246 /** \param omega The frequency value 247 \param eta The regularization parameter 248 \param orb_alpha The first orbital index 249 \param orb_beta The second orbital index 250 \param GSenergy The ground state energy returned by GSDavidson 251 \param GSvector The ground state vector as calculated by GSDavidson 252 \param RePartGF On exit RePartGF[0] contains the real part of the forward propagating part of the density response Green's function 253 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the forward propagating part of the density response Green's function 254 \param TwoRDMreal If not NULL, on exit the 2-RDM of Re{ [ omega - Ham + GSenergy + I*eta ]^{-1} ( n_beta - <GSvector| n_beta |GSvector> ) |GSvector> } 255 \param TwoRDMimag If not NULL, on exit the 2-RDM of Im{ [ omega - Ham + GSenergy + I*eta ]^{-1} ( n_beta - <GSvector| n_beta |GSvector> ) |GSvector> } 256 \param TwoRDMdens If not NULL, on exit the 2-RDM of ( n_beta - <GSvector| n_beta |GSvector> ) |GSvector> */ 257 void DensityResponseGF_forward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMdens=NULL) const; 258 259 //! Calculate the backward propagating part of the density response Green's function: <GSvector| ( n_beta - <GSvector| n_beta |GSvector> ) [ omega + Ham - GSenergy + I*eta ]^{-1} ( n_alpha - <GSvector| n_alpha |GSvector> ) |GSvector> 260 /** \param omega The frequency value 261 \param eta The regularization parameter 262 \param orb_alpha The first orbital index 263 \param orb_beta The second orbital index 264 \param GSenergy The ground state energy returned by GSDavidson 265 \param GSvector The ground state vector as calculated by GSDavidson 266 \param RePartGF On exit RePartGF[0] contains the real part of the backward propagating part of the density response Green's function 267 \param ImPartGF On exit ImPartGF[0] contains the imaginary part of the backward propagating part of the density response Green's function 268 \param TwoRDMreal If not NULL, on exit the 2-RDM of Re{ [ omega + Ham - GSenergy + I*eta ]^{-1} ( n_alpha - <GSvector| n_alpha |GSvector> ) |GSvector> } 269 \param TwoRDMimag If not NULL, on exit the 2-RDM of Im{ [ omega + Ham - GSenergy + I*eta ]^{-1} ( n_alpha - <GSvector| n_alpha |GSvector> ) |GSvector> } 270 \param TwoRDMdens If not NULL, on exit the 2-RDM of ( n_alpha - <GSvector| n_alpha |GSvector> ) |GSvector> */ 271 void DensityResponseGF_backward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal=NULL, double * TwoRDMimag=NULL, double * TwoRDMdens=NULL) const; 272 273 //! Calculate the solution of the equation ( alpha + beta * Hamiltonian + I * eta ) Solution = RHS with conjugate gradient 274 /** \param alpha The real part of the scalar in the operator 275 \param beta The real-valued prefactor of the Hamiltonian in the operator 276 \param eta The imaginary part of the scalar in the operator 277 \param RHS The real-valued right-hand side of the equation with length getVecLength(0) 278 \param RealSol On exit this array of length getVecLength(0) contains the real part of the solution 279 \param ImagSol On exit this array of length getVecLength(0) contains the imaginary part of the solution 280 \param checkError If true, the RMS error without preconditioner will be calculated and printed after convergence */ 281 void CGSolveSystem(const double alpha, const double beta, const double eta, double * RHS, double * RealSol, double * ImagSol, const bool checkError=true) const; 282 283 //void CheckHamDEBUG() const; 284 285 //! Function which returns a FCI coefficient 286 /** \param bits_up The bit string representation of the up or alpha electron Slater determinant 287 \param bits_down The bit string representation of the down or beta electron Slater determinant 288 \param vector The FCI vector with getVecLength(0) variables from which a coefficient is desired 289 \return The corresponding FCI coefficient; 0.0 if bits_up and bits_down do not form a valid FCI determinant */ 290 double getFCIcoeff(int * bits_up, int * bits_down, double * vector) const; 291 292 protected: 293 294 //==========> Functions involving Hamiltonian matrix elements 295 296 //! Function which returns the diagonal elements of the FCI Hamiltonian (without Econstant!!) = Slater determinant energies 297 /** \param diag Vector with getVecLength(0) variables which contains on exit the diagonal elements of the FCI Hamiltonian */ 298 void DiagHam(double * diag) const; 299 300 //! Function which returns the diagonal elements of the FCI Hamiltonian squared (without Econstant!!) 301 /** \param output Vector with getVecLength(0) variables which contains on exit the diagonal elements of the FCI Hamiltonian squared */ 302 void DiagHamSquared(double * output) const; 303 304 //! Function which performs the Hamiltonian times Vector product (without Econstant!!) making use of the (ij|kl) = (ji|kl) = (ij|lk) = (ji|lk) symmetry of the electron repulsion integrals 305 /** \param input The vector of length getVecLength(0) on which the Hamiltonian should act 306 \param output Vector of length getVecLength(0) which contains on exit the Hamiltonian times input */ 307 void matvec( double * input, double * output ) const; 308 309 //! Sandwich the Hamiltonian between two Slater determinants (return a specific element) (without Econstant!!) 310 /** \param bits_bra_up Bit representation of the <bra| Slater determinant of the up (alpha) electrons (length L) 311 \param bits_bra_down Bit representation of the <bra| Slater determinant of the down (beta) electrons (length L) 312 \param bits_ket_up Bit representation of the |ket> Slater determinant of the up (alpha) electrons (length L) 313 \param bits_ket_down Bit representation of the |ket> Slater determinant of the down (beta) electrons (length L) 314 \param work Work array of length 8 315 \return The FCI Hamiltonian element which connects the given two Slater determinants */ 316 double GetMatrixElement(int * bits_bra_up, int * bits_bra_down, int * bits_ket_up, int * bits_ket_down, int * work) const; 317 318 //==========> Basic conversions between bit string representations 319 320 //! Find the bit representation of a global counter corresponding to " E_ij | FCI vector > " ; where irrep_center = I_i x I_j 321 /** \param irrep_center The single electron excitation irrep I_i x I_j 322 \param counter The given global counter corresponding to " E_ij | FCI vector > " 323 \param bits_up Array of length L to store the bit representation of the up (alpha) electrons in 324 \param bits_down Array of length L to store the bit representation of the down (beta) electrons in */ 325 void getBitsOfCounter(const int irrep_center, const unsigned int counter, int * bits_up, int * bits_down) const; 326 327 //! Convertor between two representations of a same spin-projection Slater determinant 328 /** \param Lvalue The number of orbitals 329 \param bitstring The input integer, whos bits are the occupation numbers of the orbitals 330 \param bits Contains on exit the Lvalue bits of bitstring */ 331 static void str2bits(const unsigned int Lvalue, const unsigned int bitstring, int * bits); 332 333 //! Convertor between two representations of a same spin-projection Slater determinant 334 /** \param Lvalue The number of orbitals 335 \param bits Array with the Lvalue bits which should be combined to a single integer 336 \return The single integer which represents the bits in bits */ 337 static unsigned int bits2str(const unsigned int Lvalue, int * bits); 338 339 //! Find the irrep of the up Slater determinant of a global counter corresponding to " E_ij | FCI vector > " ; where irrep_center = I_i x I_j 340 /** \param irrep_center The single electron excitation irrep I_i x I_j 341 \param counter The given global counter corresponding to " E_ij | FCI vector > " 342 \return The corresponding irrep of the up Slater determinant */ 343 int getUpIrrepOfCounter(const int irrep_center, const unsigned int counter) const; 344 345 //==========> Some lapack like routines 346 347 //! Take the inproduct of two vectors 348 /** \param vecLength The vector length 349 \param vec1 The first vector 350 \param vec2 The second vector 351 \return The inproduct < vec1 | vec2 > */ 352 static double FCIddot(const unsigned int vecLength, double * vec1, double * vec2); 353 354 //! Copy a vector 355 /** \param vecLength The vector length 356 \param origin Vector to be copied 357 \param target Where to copy the vector to */ 358 static void FCIdcopy(const unsigned int vecLength, double * origin, double * target); 359 360 //! Calculate the 2-norm of a vector 361 /** \param vecLength The vector length 362 \param vec The vector 363 \return The 2-norm of vec */ 364 static double FCIfrobeniusnorm(const unsigned int vecLength, double * vec); 365 366 //! Do lapack's daxpy vec_y += alpha * vec_x 367 /** \param vecLength The vector length 368 \param alpha The scalar factor 369 \param vec_x The vector which has to be added to vec_y in rescaled form 370 \param vec_y The target vector */ 371 static void FCIdaxpy(const unsigned int vecLength, const double alpha, double * vec_x, double * vec_y); 372 373 //! Do lapack's dscal vec *= alpha 374 /** \param vecLength The vector length 375 \param alpha The scalar factor 376 \param vec The vector which has to be rescaled */ 377 static void FCIdscal(const unsigned int vecLength, const double alpha, double * vec); 378 379 //==========> Protected functions regarding the Green's functions 380 381 //! Set thisVector to a creator/annihilator acting on otherVector 382 /** \param whichOperator With which operator should be acted on the other FCI state: C means creator and A means annihilator 383 \param isUp Boolean which denotes if the operator corresponds to an up (alpha) or down (beta) electron 384 \param orbIndex Orbital index on which the operator acts 385 \param thisVector Vector with length getVecLength(0) where the result of the operation should be stored 386 \param otherFCI FCI instance which corresponds to the FCI vector otherVector on which is acted 387 \param otherVector Vector with length otherFCI->getVecLength(0) which contains the FCI vector on which is acted */ 388 void ActWithSecondQuantizedOperator(const char whichOperator, const bool isUp, const unsigned int orbIndex, double * thisVector, const FCI * otherFCI, double * otherVector) const; 389 390 //! Set resultVector to the number operator of a specific site acting on sourceVector 391 /** \param orbIndex Orbital index of the number operator 392 \param resultVector Vector with length getVecLength(0) where the result of the operation should be stored 393 \param sourceVector Vector with length getVecLength(0) on which the number operator acts */ 394 void ActWithNumberOperator(const unsigned int orbIndex, double * resultVector, double * sourceVector) const; 395 396 //! Calculate the solution of the system Operator |Sol> = |RESID> with Operator = precon * [ ( alpha + beta * H )^2 + eta^2 ] * precon 397 /** \param alpha The parameter alpha of the operator 398 \param beta The parameter beta of the operator 399 \param eta The parameter eta of the operator 400 \param precon The diagonal preconditioner 401 \param Sol On entry this array of size getVecLength(0) contains the initial guess; on exit the solution is stored here 402 \param RESID On entry the RHS of the equation is stored in this array of size getVecLength(0), on exit it is overwritten 403 \param PVEC Workspace of size getVecLength(0) 404 \param OxPVEC Workspace of size getVecLength(0) 405 \param temp Workspace of size getVecLength(0) 406 \param temp2 Workspace of size getVecLength(0) */ 407 void CGCoreSolver(const double alpha, const double beta, const double eta, double * precon, double * Sol, double * RESID, double * PVEC, double * OxPVEC, double * temp, double * temp2) const; 408 409 //! Calculate out = (alpha + beta * Hamiltonian) * in (Econstant is taken into account!!) 410 /** \param alpha The parameter alpha of the operator 411 \param beta The parameter beta of the operator 412 \param in On entry the vector of size getVecLength(0) on which the operator should be applied; unchanged on exit 413 \param out Array of size getVecLength(0), which contains on exit (alpha + beta * Hamiltonian) * in */ 414 void CGAlphaPlusBetaHAM(const double alpha, const double beta, double * in, double * out) const; 415 416 //! Calculate out = [(alpha + beta * Hamiltonian)^2 + eta^2] * in (Econstant is taken into account!!) 417 /** \param alpha The parameter alpha of the operator 418 \param beta The parameter beta of the operator 419 \param eta The parameter eta of the operator 420 \param in On entry the vector of size getVecLength(0) on which the operator should be applied; unchanged on exit 421 \param out Array of size getVecLength(0), which contains on exit [(alpha + beta * Hamiltonian)^2 + eta^2] * in 422 \param temp Workspace of size getVecLength(0) */ 423 void CGoperator(const double alpha, const double beta, const double eta, double * in, double * temp, double * out) const; 424 425 //! Calculate (without approximation) diagonal = diag[ (alpha + beta * Hamiltonian)^2 + eta^2 ] (Econstant is taken into account!!) 426 /** \param alpha The parameter alpha of the operator 427 \param beta The parameter beta of the operator 428 \param eta The parameter eta of the operator 429 \param diagonal Array of size getVecLength(0), which contains on exit diag[ (alpha + beta * Hamiltonian)^2 + eta^2 ] 430 \param workspace Workspace of size getVecLength(0) */ 431 void CGdiagonal(const double alpha, const double beta, const double eta, double * diagonal, double * workspace) const; 432 433 //==========> Protected functions regarding the higher order RDMs 434 435 //! Apply E_{crea,anni} to |orig_vector> and store the result in |result_vector> 436 /** \param orig_vector The original vector 437 \param result_vector The result vector 438 \param crea The orbital index of the creator 439 \param anni The orbital index of the annihilator 440 \param orig_target_irrep The irrep of the orig_vector */ 441 void apply_excitation( double * orig_vector, double * result_vector, const int crea, const int anni, const int orig_target_irrep ) const; 442 443 private: 444 445 //! The FCI verbose level: 0 print nothing, 1 print start and solution, 2 print everything 446 int FCIverbose; 447 448 //! The maximum number of MB which can be used to store both HXVworkbig1 and HXVworkbig2 449 double maxMemWorkMB; 450 451 //! The constant term of the Hamiltonian 452 double Econstant; 453 454 //! The AUGMENTED one-body matrix elements Gmat[ i + L * j ] = T[ i + L * j ] - 0.5 * sum_k getERI(i, k, k, j) 455 double * Gmat; 456 457 //! The electron repulsion integrals (in chemists notation: \int dr1 dr2 orb1(r1) * orb2(r1) * orb3(r2) * orb4(r2) / |r1-r2| = ERI[ orb1 + L * ( orb2 + L * ( orb3 + L * orb4 ) ) ]) 458 double * ERI; 459 460 //! The number of irreps of the molecule's Abelian point group with real-valued character table 461 unsigned int num_irreps; 462 463 //! The targeted irrep of the total FCI wavefunction ( 0 <= TargetIrrep < num_irreps ) 464 int TargetIrrep; 465 466 //! Array of length L which contains for each orbital the irrep ( 0 <= orb < L : 0 <= orb2irrep[ orb ] < num_irreps ) 467 int * orb2irrep; 468 469 //! The number of orbitals ( Hubbard type orbitals with 4 possible occupations each ) 470 unsigned int L; 471 472 //! The number of up (alpha) electrons 473 unsigned int Nel_up; 474 475 //! The number of down (beta) electrons 476 unsigned int Nel_down; 477 478 //! The number of up (alpha) Slater determinants with irrep "irrep" and Nel_up electrons is given by numPerIrrep_up[ irrep ] 479 unsigned int * numPerIrrep_up; 480 481 //! The number of down (beta) Slater determinants with irrep "irrep" and Nel_down electrons is given by numPerIrrep_down[ irrep ] 482 unsigned int * numPerIrrep_down; 483 484 //! For irrep "irrep" and bit string "bitstring", str2cnt_up[ irrep ][ bitstring ] returns the counter (0 <= counter < numPerIrrep_up[ irrep ]) of the up (alpha) Slater determinant within the irrep block (or -1 if invalid) 485 int ** str2cnt_up; 486 487 //! For irrep "irrep" and bit string "bitstring", str2cnt_down[ irrep ][ bitstring ] returns the counter (0 <= counter < numPerIrrep_down[ irrep ]) of the down (beta) Slater determinant within the irrep block (or -1 if invalid) 488 int ** str2cnt_down; 489 490 //! For irrep "irrep" and counter of the up (alpha) Slater determinant "counter" (0 <= counter < numPerIrrep_up[ irrep ]) cnt2str_up[ irrep ][ counter ] returns the bitstring representation of the corresponding up (alpha) Slater determinant 491 unsigned int ** cnt2str_up; 492 493 //! For irrep "irrep" and counter of the down (beta) Slater determinant "counter" (0 <= counter < numPerIrrep_down[ irrep ]) cnt2str_down[ irrep ][ counter ] returns the bitstring representation of the corresponding down (beta) Slater determinant 494 unsigned int ** cnt2str_down; 495 496 //! For irrep "irrep_result" and up (alpha) Slater determinant counter "result" lookup_cnt_alpha[ irrep_result ][ i + L * ( j + L * result ) ] returns the counter "origin" which corresponds to | result > = +/- E^{alpha}_ij | origin > or | origin > = +/- E^{alpha}_ji | result > 497 int *** lookup_cnt_alpha; 498 499 //! For irrep "irrep_result" and down (beta) Slater determinant counter "result" lookup_cnt_beta[ irrep_result ][ i + L * ( j + L * result ) ] returns the counter "origin" which corresponds to | result > = +/- E^{beta}_ij | origin > or | origin > = +/- E^{beta}_ji | result > 500 int *** lookup_cnt_beta; 501 502 //! For irrep "irrep_result" and up (alpha) Slater determinant counter "result" lookup_sign_alpha[ irrep_result ][ i + L * ( j + L * result ) ] returns the sign s which corresponds to | result > = s * E^{alpha}_ij | origin > 503 int *** lookup_sign_alpha; 504 505 //! For irrep "irrep_result" and down (beta) Slater determinant counter "result" lookup_sign_beta[ irrep_result ][ i + L * ( j + L * result ) ] returns the sign s which corresponds to | result > = s * E^{beta}_ij | origin > 506 int *** lookup_sign_beta; 507 508 //! For irrep_center = irrep_creator x irrep_annihilator the number of corresponding excitation pairs E_{creator <= annihilator} is given by irrep_center_num[ irrep_center ] 509 unsigned int * irrep_center_num; 510 511 //! For irrep_center = irrep_creator x irrep_annihilator the creator orbital index of the nth (0 <= n < irrep_center_num[ irrep_center ]) excitation pair is given by irrep_center_crea_orb[ irrep_center ][ n ] 512 unsigned int ** irrep_center_crea_orb; 513 514 //! For irrep_center = irrep_creator x irrep_annihilator the annihilator orbital index of the nth (0 <= n < irrep_center_num[ irrep_center ]) excitation pair is given by irrep_center_anni_orb[ irrep_center ][ n ] 515 unsigned int ** irrep_center_anni_orb; 516 517 //! The global index corresponding to a vector " E_{ij} | FCI vector > " with irrep_center = irrep_i x irrep_j is given by irrep_center_jumps[ irrep_center ][ irrep_alpha ] + count_alpha + numPerIrrep_up[ irrep_alpha ] * cnt_beta where count_alpha and count_beta are the up (alpha) and down (beta) Slater determinants with resp. irreps irrep_alpha and irrep_beta = TargetIrrep x irrep_alpha x irrep_center 518 unsigned int ** irrep_center_jumps; 519 520 //! Number of doubles in each of the HVXworkbig arrays 521 unsigned long long HXVsizeWorkspace; 522 523 //! Work space of size L*L*L*L 524 double * HXVworksmall; 525 526 //! Work space of size HXVsizeWorkspace 527 double * HXVworkbig1; 528 529 //! Work space of size HXVsizeWorkspace 530 double * HXVworkbig2; 531 532 //! Initialize a part of the private variables 533 void StartupCountersVsBitstrings(); 534 535 //! Initialize a part of the private variables 536 void StartupLookupTables(); 537 538 //! Initialize a part of the private variables 539 void StartupIrrepCenter(); 540 541 //! Actual routine used by Fill3RDM, Fock4RDM, Diag4RDM 542 double Driver3RDM(double * vector, double * output, double * three_rdm, double * fock, const unsigned int orbz) const; 543 544 //! Alpha excitation kernels 545 static void excite_alpha_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int dim_down, double * origin, double * result, int * signmap, int * countmap ); 546 static void excite_alpha_first( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ); 547 static void excite_alpha_second_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ); 548 549 //! Beta excitation kernels 550 static void excite_beta_omp( const unsigned int dim_up, const unsigned int dim_new_down, double * origin, double * result, int * signmap, int * countmap ); 551 static void excite_beta_first( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ); 552 static void excite_beta_second_omp( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ); 553 554 }; 555 556 } 557 558 #endif 559