1 /* 2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry 3 Copyright (C) 2013-2021 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 DMRG_CHEMPS2_H 21 #define DMRG_CHEMPS2_H 22 23 #include <string> 24 25 #include "Options.h" 26 #include "Problem.h" 27 #include "TensorT.h" 28 #include "TensorX.h" 29 #include "TensorL.h" 30 #include "SyBookkeeper.h" 31 #include "TensorF1.h" 32 #include "TensorF0.h" 33 #include "TensorS0.h" 34 #include "TensorS1.h" 35 #include "TensorOperator.h" 36 #include "TensorQ.h" 37 #include "TensorO.h" 38 #include "Tensor3RDM.h" 39 #include "TensorGYZ.h" 40 #include "TensorKM.h" 41 #include "TwoDM.h" 42 #include "ThreeDM.h" 43 #include "Correlations.h" 44 #include "Heff.h" 45 #include "Sobject.h" 46 #include "ConvergenceScheme.h" 47 48 #include <hdf5.h> 49 50 //For the timings of the different parts of DMRG 51 #define CHEMPS2_TIME_S_JOIN 0 52 #define CHEMPS2_TIME_S_SOLVE 1 53 #define CHEMPS2_TIME_S_SPLIT 2 54 #define CHEMPS2_TIME_TENS_TOTAL 3 55 #define CHEMPS2_TIME_TENS_ALLOC 4 56 #define CHEMPS2_TIME_TENS_FREE 5 57 #define CHEMPS2_TIME_DISK_WRITE 6 58 #define CHEMPS2_TIME_DISK_READ 7 59 #define CHEMPS2_TIME_TENS_CALC 8 60 #define CHEMPS2_TIME_VECLENGTH 9 61 62 namespace CheMPS2{ 63 /** DMRG class. 64 \author Sebastian Wouters <sebastianwouters@gmail.com> 65 \date July 31, 2013 66 67 The DMRG class solves the Problem with its given parameters. A fully SU(2) symmetric MPS wavefunction is variationally optimized in a two-site sweep algorithm. When the solution has been reached, the converged energy and spin contracted 2DMs can be accessed. For more information, please take a look at 68 69 S. Wouters, W. Poelmans, P.W. Ayers and D. Van Neck, \n 70 CheMPS2: a free open-source spin-adapted implementation of the density matrix renormalization group for ab initio quantum chemistry, \n 71 Computer Physics Communications 185, 1501-1514 (2014) \n 72 http://dx.doi.org/10.1016/j.cpc.2014.01.019 \n 73 http://arxiv.org/abs/1312.2415 74 75 S. Wouters and D. Van Neck, \n 76 The density matrix renormalization group for ab initio quantum chemistry, \n 77 European Physical Journal D 68, 272 (2014) \n 78 http://dx.doi.org/10.1140/epjd/e2014-50500-1 \n 79 http://arxiv.org/abs/1407.2040 80 81 The user manual: http://sebwouters.github.io/CheMPS2/index.html 82 */ 83 class DMRG{ 84 85 public: 86 87 //! Constructor 88 /** \param Probin The problem to be solved 89 \param OptSchemeIn The optimization scheme for the DMRG sweeps 90 \param makechkpt Whether or not to save MPS checkpoints in the working directory 91 \param tmpfolder Temporary folder on a large partition to store the renormalized operators on disk (by default "/tmp") 92 \param occupancies ROHF occupancies of a determinant to enlarge in the initial guess in HAM order */ 93 DMRG(Problem * Probin, ConvergenceScheme * OptSchemeIn, const bool makechkpt=CheMPS2::DMRG_storeMpsOnDisk, const string tmpfolder=CheMPS2::defaultTMPpath, int * occupancies=NULL); 94 95 //! Destructor 96 virtual ~DMRG(); 97 98 //! Solver 99 /** \return The min. energy encountered so far during the sweeps. */ 100 double Solve(); 101 102 //! Reconstruct the renormalized operators when you overwrite the matrix elements with Prob->setMxElement() 103 void PreSolve(); 104 105 //! Calculate the 2-RDM and correlations. Afterwards the MPS is again in LLLLLLLC gauge. calc2DMandCorrelations()106 void calc2DMandCorrelations(){ calc_rdms_and_correlations(false); } 107 108 //! Calculate the reduced density matrices and correlations. Afterwards the MPS is again in LLLLLLLC gauge. 109 /** \param do_3rdm Whether or not to calculate the 3-RDM 110 \param disk_3rdm Whether or not to use disk in order to avoid storing the full 3-RDM of size L^6 */ 111 void calc_rdms_and_correlations( const bool do_3rdm, const bool disk_3rdm = false ); 112 113 //! Get the pointer to the 2-RDM 114 /** \return The 2-RDM. Returns a NULL pointer if not yet calculated. */ get2DM()115 TwoDM * get2DM(){ return the2DM; } 116 117 //! Get the pointer to the 3-RDM 118 /** \return The 3-RDM. Returns a NULL pointer if not yet calculated. */ get3DM()119 ThreeDM * get3DM(){ return the3DM; } 120 121 //! Obtain the symmetrized 4-RDM terms 0.5 * ( Gamma4_ijkl,pqrt + Gamma4_ijkt,pqrl ) with l and t fixed, after the 3-RDM has been calculated. 122 /** \param output Array to store the symmetrized 4-RDM terms in Hamiltonian index notation: output[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] = 0.5 * ( Gamma4_ijkl,pqrt + Gamma4_ijkt,pqrl ). 123 \param ham_orb1 The Hamiltonian index of the first fixed orbital. 124 \param ham_orb2 The Hamiltonian index of the second fixed orbital. 125 \param last_case If true, everything will be set up to allow to continue sweeping. */ 126 void Symm4RDM( double * output, const int ham_orb1, const int ham_orb2, const bool last_case ); 127 128 //! Get the pointer to the Correlations 129 /** \return The Correlations. Returns a NULL pointer if not yet calculated. */ getCorrelations()130 Correlations * getCorrelations(){ return theCorr; } 131 132 //! Get a specific FCI coefficient. The array coeff contains the occupation numbers of the L Hamiltonian orbitals. It is assumed that the unpaired electrons are all alpha electrons, and that this number equals twice the total targeted spin. 133 /** \param coeff Array containing the occupation numbers of the L Hamiltonian orbitals (occupations can be 0, 1, or 2). 134 \return The desired FCI coefficient */ 135 double getSpecificCoefficient(int * coeff) const; 136 137 //! Get a specific FCI coefficient. The arrays alpha and beta contain the alpha and beta occupation numbers of the L Hamiltonian orbitals. 138 /** \param alpha Array containing the alpha electron occupation numbers of the L Hamiltonian orbitals (occupations can be 0 or 1). 139 \param beta Array containing the beta electron occupation numbers of the L Hamiltonian orbitals (occupations can be 0 or 1). 140 \param mpi_chemps2_master_only When running with MPI, whether only the master process should calculate the FCI coefficient. If false, any process can calculate the FCI coefficient, and it won't be broadcasted (allows for parallel FCI coefficient calculation). 141 \return The desired FCI coefficient */ 142 double getFCIcoefficient(int * alpha, int * beta, const bool mpi_chemps2_master_only=true) const; 143 144 //! Call "rm " + CheMPS2::DMRG_MPS_storage_prefix + "*.h5" 145 void deleteStoredMPS(); 146 147 //! Call "rm " + tempfolder + "/" + CheMPS2::DMRG_OPERATOR_storage_prefix + string(thePID) + "*.h5"; 148 void deleteStoredOperators(); 149 150 //! Activate the necessary storage and machinery to handle excitations 151 /** \param maxExcIn The max. number of excitations desired */ 152 void activateExcitations(const int maxExcIn); 153 154 //! Push back current calculation and set everything up to calculate a (new) excitation 155 /** \param EshiftIn To the Hamiltonian, a level shift is introduced to exclude the previously calculated MPS: Hnew = Hold + EshiftIn * | prev> <prev| */ 156 void newExcitation(const double EshiftIn); 157 158 //! Print the license 159 static void PrintLicense(); 160 161 //! Get the number of MPS variables 162 int get_num_mps_var() const; 163 164 private: 165 166 //Setup the DMRG SyBK and MPS (in separate function to allow pushbacks and recreations for excited states) 167 void setupBookkeeperAndMPS( int * occupancies=NULL ); 168 169 //! DMRG MPS + virt. dim. storage filename 170 string MPSstoragename; 171 172 //The optimization scheme for the DMRG sweeps (externally allocated, filled and deleted) 173 ConvergenceScheme * OptScheme; 174 175 //Whether or not the MPS was loaded from disk to memory at the start 176 bool loadedMPS; 177 178 //Integer to distinguish storage between different calculations 179 int thePID; 180 181 //Pointer to the Problem --> constructed and destructed outside of this class 182 Problem * Prob; 183 184 //The number of orbitals: copied here so that the DMRG destructor doesn't depend on whether Prob still exists 185 int L; 186 187 //Minimum energy encountered during all performed micro-iterations (as opposed to the 2DM/edge energy) 188 double TotalMinEnergy; 189 190 //Minimum energy encountered during the micro-iterations of the last performed sweep (as opposed to the 2DM/edge energy) 191 double LastMinEnergy; 192 193 //Max. discarded weight of last sweep 194 double MaxDiscWeightLastSweep; 195 196 //Symmetry information object 197 SyBookkeeper * denBK; 198 199 //The MPS 200 TensorT ** MPS; 201 202 //The TwoDM 203 TwoDM * the2DM; 204 205 //The ThreeDM 206 ThreeDM * the3DM; 207 208 //The Correlations 209 Correlations * theCorr; 210 211 //Whether or not allocated 212 int * isAllocated; 213 214 //TensorL's 215 TensorL *** Ltensors; 216 217 //TensorX's 218 TensorX ** Xtensors; 219 220 //TensorF0's 221 TensorF0 **** F0tensors; 222 223 //TensorF1's 224 TensorF1 **** F1tensors; 225 226 //TensorS0's 227 TensorS0 **** S0tensors; 228 229 //TensorS1's 230 TensorS1 **** S1tensors; 231 232 //ABCD-tensors 233 TensorOperator **** Atensors; 234 TensorOperator **** Btensors; 235 TensorOperator **** Ctensors; 236 TensorOperator **** Dtensors; 237 238 //TensorQ's 239 TensorQ *** Qtensors; 240 241 //Tensors required for the 3-RDM calculation 242 Tensor3RDM ***** tensor_3rdm_a_J0_doublet; 243 Tensor3RDM ***** tensor_3rdm_a_J1_doublet; 244 Tensor3RDM ***** tensor_3rdm_a_J1_quartet; 245 Tensor3RDM ***** tensor_3rdm_b_J0_doublet; 246 Tensor3RDM ***** tensor_3rdm_b_J1_doublet; 247 Tensor3RDM ***** tensor_3rdm_b_J1_quartet; 248 Tensor3RDM ***** tensor_3rdm_c_J0_doublet; 249 Tensor3RDM ***** tensor_3rdm_c_J1_doublet; 250 Tensor3RDM ***** tensor_3rdm_c_J1_quartet; 251 Tensor3RDM ***** tensor_3rdm_d_J0_doublet; 252 Tensor3RDM ***** tensor_3rdm_d_J1_doublet; 253 Tensor3RDM ***** tensor_3rdm_d_J1_quartet; 254 255 //Tensors required for the Correlations calculation 256 TensorGYZ ** Gtensors; 257 TensorGYZ ** Ytensors; 258 TensorGYZ ** Ztensors; 259 TensorKM ** Ktensors; 260 TensorKM ** Mtensors; 261 262 // Sweeps 263 double sweepleft( const bool change, const int instruction, const bool am_i_master ); 264 double sweepright( const bool change, const int instruction, const bool am_i_master ); 265 double solve_site( const int index, const double dvdson_rtol, const double noise_level, const int virtual_dimension, const bool am_i_master, const bool moving_right, const bool change ); 266 267 //Load and save functions 268 void MY_HDF5_WRITE_BATCH(const hid_t file_id, const int number, Tensor ** batch, const long long totalsize, const std::string tag); 269 void MY_HDF5_READ_BATCH( const hid_t file_id, const int number, Tensor ** batch, const long long totalsize, const std::string tag); 270 void OperatorsOnDisk(const int index, const bool movingRight, const bool store); 271 string tempfolder; 272 273 void saveMPS(const std::string name, TensorT ** MPSlocation, SyBookkeeper * BKlocation, bool isConverged) const; 274 void loadDIM(const std::string name, SyBookkeeper * BKlocation); 275 void loadMPS(const std::string name, TensorT ** MPSlocation, bool * isConverged); 276 bool makecheckpoints; 277 278 //Helper functions for making the boundary operators 279 void updateMovingRight(const int index); 280 void updateMovingLeft(const int index); 281 void deleteTensors(const int index, const bool movingRight); 282 void allocateTensors(const int index, const bool movingRight); 283 void updateMovingRightSafe(const int cnt); 284 void updateMovingRightSafeFirstTime(const int cnt); 285 void updateMovingRightSafe2DM(const int cnt); 286 void updateMovingLeftSafe(const int cnt); 287 void updateMovingLeftSafeFirstTime(const int cnt); 288 void updateMovingLeftSafe2DM(const int cnt); 289 void deleteAllBoundaryOperators(); 290 291 // Helper functions for making the 3-RDM boundary operators 292 void update_safe_3rdm_operators( const int boundary ); 293 void allocate_3rdm_operators( const int boundary ); 294 void update_3rdm_operators( const int boundary ); 295 void delete_3rdm_operators( const int boundary ); 296 297 // Helper functions for making the symmetrized 4-RDM 298 void solve_fock( const int dmrg_orb1, const int dmrg_orb2, const double alpha, const double beta ); 299 static void solve_fock_update_helper( const int index, const int dmrg_orb1, const int dmrg_orb2, const bool moving_right, TensorT ** new_mps, TensorT ** old_mps, SyBookkeeper * new_bk, SyBookkeeper * old_bk, TensorO ** overlaps, TensorL ** regular, TensorL ** trans ); 300 static void left_normalize( TensorT * left_mps, TensorT * right_mps ); 301 static void right_normalize( TensorT * left_mps, TensorT * right_mps ); 302 void symm_4rdm_helper( double * output, const int ham_orb1, const int ham_orb2, const double alpha, const double beta, const bool add, const double factor ); 303 304 //Helper functions for making the Correlations boundary operators 305 void update_correlations_tensors(const int siteindex); 306 307 //The storage and functions to handle excited states 308 int nStates; 309 bool Exc_activated; 310 int maxExc; 311 double * Exc_Eshifts; 312 TensorT *** Exc_MPSs; 313 SyBookkeeper ** Exc_BKs; 314 TensorO *** Exc_Overlaps; 315 double ** prepare_excitations(Sobject * denS); 316 void cleanup_excitations(double ** VeffTilde) const; 317 void calcVeffTilde(double * result, Sobject * currentS, int state_number); 318 void calc_overlaps( const bool moving_right ); 319 320 // Performance counters 321 double timings[ CHEMPS2_TIME_VECLENGTH ]; 322 long long num_double_write_disk; 323 long long num_double_read_disk; 324 void print_tensor_update_performance() const; 325 326 }; 327 } 328 329 #endif 330