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 TWODM_CHEMPS2_H 21 #define TWODM_CHEMPS2_H 22 23 #include "TensorT.h" 24 #include "TensorL.h" 25 #include "TensorF0.h" 26 #include "TensorF1.h" 27 #include "TensorS0.h" 28 #include "TensorS1.h" 29 #include "SyBookkeeper.h" 30 31 namespace CheMPS2{ 32 /** TwoDM class. 33 \author Sebastian Wouters <sebastianwouters@gmail.com> 34 \date June 13, 2013 35 36 The TwoDM class stores the result of a converged DMRG calculation. With the 2DM \n 37 \f$ \Gamma_{(\alpha \sigma) (\beta \tau) ; (\gamma \sigma) (\delta \tau)} = \braket{ a^{\dagger}_{\alpha \sigma} a^{\dagger}_{\beta \tau} a_{\delta \tau} a_{\gamma \sigma} } \f$\n 38 we can define two spin-reduced versions of interest:\n 39 \f$ \Gamma^{2A}_{\alpha \beta ; \gamma \delta} = \sum_{\sigma \tau} \Gamma_{(\alpha \sigma) (\beta \tau) ; (\gamma \sigma) (\delta \tau)} \f$ \n 40 \f$ \Gamma^{2B}_{\alpha \beta ; \gamma \delta} = \sum_{\sigma} \left( \Gamma_{(\alpha \sigma) (\beta \sigma) ; (\gamma \sigma) (\delta \sigma)} - \Gamma_{(\alpha \sigma) (\beta -\sigma) ; (\gamma \sigma) (\delta -\sigma)} \right) \f$. \n 41 Because the wave-function belongs to a certain Abelian irrep, \f$ I_{\alpha} \otimes I_{\beta} \otimes I_{\gamma} \otimes I_{\delta} = I_{trivial} \f$ must be valid before the corresponding element \f$ \Gamma^{A,B}_{\alpha \beta ; \gamma \delta} \f$ is non-zero. \n 42 \n 43 We can also define spin-densities in the spin-ensemble as:\n 44 \f{eqnarray*}{ 45 \Gamma^{spin}_{ij} & = & \frac{3(1 - \delta_{S,0})}{(S+1)(2S+1)} \sum_{S^z} S^z \braket{ S S^z \mid a^{\dagger}_{i \uparrow} a_{j \uparrow} - a^{\dagger}_{i \downarrow} a_{j \downarrow} \mid S S^z } \\ 46 & = & \frac{3(1 - \delta_{S,0})}{2(S+1)} \left[ ( 2 - N_{elec} ) \Gamma^1_{ij} - \sum_r \Gamma^{2A}_{ir;rj} - \sum_r \Gamma^{2B}_{ir;rj} \right] 47 \f} \n 48 The normalization factor is chosen so that \f$ trace( \Gamma^{spin} ) = 2 S \f$. 49 */ 50 class TwoDM{ 51 52 public: 53 54 //! Constructor 55 /** \param denBKIn Symmetry sector bookkeeper 56 \param ProbIn The problem to be solved */ 57 TwoDM(const SyBookkeeper * denBKIn, const Problem * ProbIn); 58 59 //! Destructor 60 virtual ~TwoDM(); 61 62 //! Get a 2DM_A term, using the DMRG indices 63 /** \param cnt1 the first index 64 \param cnt2 the second index 65 \param cnt3 the third index 66 \param cnt4 the fourth index 67 \return the desired value */ 68 double getTwoDMA_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const; 69 70 //! Get a 2DM_B term, using the DMRG indices 71 /** \param cnt1 the first index 72 \param cnt2 the second index 73 \param cnt3 the third index 74 \param cnt4 the fourth index 75 \return the desired value */ 76 double getTwoDMB_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const; 77 78 //! Get a 1-RDM term, using the DMRG indices 79 /** \param cnt1 the first index 80 \param cnt2 the second index 81 \return the desired value */ 82 double get1RDM_DMRG(const int cnt1, const int cnt2) const; 83 84 //! Get a spin-density term, using the DMRG indices 85 /** \param cnt1 the first index 86 \param cnt2 the second index 87 \return the desired value */ 88 double spin_density_dmrg( const int cnt1, const int cnt2 ) const; 89 90 //! Get a 2DM_A term, using the HAM indices 91 /** \param cnt1 the first index 92 \param cnt2 the second index 93 \param cnt3 the third index 94 \param cnt4 the fourth index 95 \return the desired value */ 96 double getTwoDMA_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const; 97 98 //! Get a 2DM_B term, using the HAM indices 99 /** \param cnt1 the first index 100 \param cnt2 the second index 101 \param cnt3 the third index 102 \param cnt4 the fourth index 103 \return the desired value */ 104 double getTwoDMB_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const; 105 106 //! Get a 1-RDM term, using the HAM indices 107 /** \param cnt1 the first index 108 \param cnt2 the second index 109 \return the desired value */ 110 double get1RDM_HAM(const int cnt1, const int cnt2) const; 111 112 //! Get a spin-density term, using the HAM indices 113 /** \param cnt1 the first index 114 \param cnt2 the second index 115 \return the desired value */ 116 double spin_density_ham( const int cnt1, const int cnt2 ) const; 117 118 //! Fill the 2DM terms with as second site index denT->gIndex() 119 /** \param denT DMRG site-matrices 120 \param Ltens Ltensors 121 \param F0tens F0tensors 122 \param F1tens F1tensors 123 \param S0tens S0tensors 124 \param S1tens S1tensors*/ 125 void FillSite(TensorT * denT, TensorL *** Ltens, TensorF0 **** F0tens, TensorF1 **** F1tens, TensorS0 **** S0tens, TensorS1 **** S1tens); 126 127 //! After the whole 2-RDM is filled, a prefactor for higher multiplicities should be applied 128 void correct_higher_multiplicities(); 129 130 //! Return the double trace of 2DM-A (should be N(N-1)) 131 /** \return Double trace of 2DM-A */ 132 double trace() const; 133 134 //! Calculate the energy based on the 2DM-A 135 /** \return The energy calculated as 0.5*trace(2DM-A * Ham) */ 136 double energy() const; 137 138 //! Print the natural orbital occupation numbers 139 void print_noon() const; 140 141 //! Save the TwoDMs to disk 142 void save() const; 143 144 //! Load the TwoDMs from disk 145 void read(); 146 147 //! Save the 2-RDM-A to disk in Hamiltonian indices 148 /** \param filename The filename to store the 2-RDM at */ 149 void save_HAM( const string filename ) const; 150 151 //! Write the 2-RDM-A to a file 152 /** param filename where to write the 2-RDM-A elements to */ 153 void write2DMAfile(const string filename) const; 154 155 //! Add the 2-RDM elements of all MPI processes 156 void mpi_allreduce(); 157 158 private: 159 160 //The BK containing all the irrep information 161 const SyBookkeeper * denBK; 162 163 //The problem containing orbital reshuffling and symmetry information 164 const Problem * Prob; 165 166 //The chain length 167 int L; 168 169 //Two 2DM^{A,B} objects 170 double * two_rdm_A; 171 double * two_rdm_B; 172 173 // Set 2DM terms, using the DMRG indices 174 void set_2rdm_A_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value); 175 void set_2rdm_B_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value); 176 177 //Helper functions 178 double doD1(TensorT * denT); 179 double doD2(TensorT * denT, TensorL * Lright, double * workmem); 180 double doD3(TensorT * denT, TensorS0 * S0right, double * workmem); 181 double doD4(TensorT * denT, TensorF0 * F0right, double * workmem); 182 double doD5(TensorT * denT, TensorF0 * F0right, double * workmem); 183 double doD6(TensorT * denT, TensorF1 * F1right, double * workmem); 184 double doD7(TensorT * denT, TensorL * Lleft, double * workmem); 185 double doD8(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g); 186 void doD9andD10andD11(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, double * d9, double * d10, double * d11, int Irrep_g); 187 double doD12(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g); 188 double doD13(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g); 189 double doD14(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g); 190 double doD15(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g); 191 double doD16(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g); 192 double doD17orD21(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD17); 193 double doD18orD22(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD18); 194 double doD19orD23(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD19); 195 double doD20orD24(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD20); 196 197 }; 198 } 199 200 #endif 201