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 #include <math.h>
21 #include <stdlib.h>
22 #include <iostream>
23 #include <string.h>
24 #include <sstream>
25 #include <sys/stat.h>
26 #include <sys/time.h>
27 #include <assert.h>
28 #include <unistd.h>
29
30 #include "DMRG.h"
31 #include "MPIchemps2.h"
32
33 using std::cout;
34 using std::endl;
35
DMRG(Problem * ProbIn,ConvergenceScheme * OptSchemeIn,const bool makechkpt,const string tmpfolder,int * occupancies)36 CheMPS2::DMRG::DMRG( Problem * ProbIn, ConvergenceScheme * OptSchemeIn, const bool makechkpt, const string tmpfolder, int * occupancies ){
37
38 #ifdef CHEMPS2_MPI_COMPILATION
39 if ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER ){ PrintLicense(); }
40 #else
41 PrintLicense();
42 #endif
43
44 assert( ProbIn->checkConsistency() );
45 Prob = ProbIn;
46 L = Prob->gL();
47 Prob->construct_mxelem();
48 OptScheme = OptSchemeIn;
49 thePID = getpid(); //PID is unique for each MPI process
50 nStates = 1;
51
52 Ltensors = new TensorL ** [ L - 1 ];
53 F0tensors = new TensorF0 *** [ L - 1 ];
54 F1tensors = new TensorF1 *** [ L - 1 ];
55 S0tensors = new TensorS0 *** [ L - 1 ];
56 S1tensors = new TensorS1 *** [ L - 1 ];
57 Atensors = new TensorOperator *** [ L - 1 ];
58 Btensors = new TensorOperator *** [ L - 1 ];
59 Ctensors = new TensorOperator *** [ L - 1 ];
60 Dtensors = new TensorOperator *** [ L - 1 ];
61 Qtensors = new TensorQ ** [ L - 1 ];
62 Xtensors = new TensorX * [ L - 1 ];
63 isAllocated = new int[ L - 1 ]; // 0 not allocated; 1 moving right; 2 moving left
64
65 tensor_3rdm_a_J0_doublet = NULL;
66 tensor_3rdm_a_J1_doublet = NULL;
67 tensor_3rdm_a_J1_quartet = NULL;
68 tensor_3rdm_b_J0_doublet = NULL;
69 tensor_3rdm_b_J1_doublet = NULL;
70 tensor_3rdm_b_J1_quartet = NULL;
71 tensor_3rdm_c_J0_doublet = NULL;
72 tensor_3rdm_c_J1_doublet = NULL;
73 tensor_3rdm_c_J1_quartet = NULL;
74 tensor_3rdm_d_J0_doublet = NULL;
75 tensor_3rdm_d_J1_doublet = NULL;
76 tensor_3rdm_d_J1_quartet = NULL;
77
78 Gtensors = NULL;
79 Ytensors = NULL;
80 Ztensors = NULL;
81 Ktensors = NULL;
82 Mtensors = NULL;
83
84 for ( int cnt = 0; cnt < L - 1; cnt++ ){ isAllocated[ cnt ] = 0; }
85 for ( int timecnt = 0; timecnt < CHEMPS2_TIME_VECLENGTH; timecnt++ ){ timings[ timecnt ] = 0.0; }
86 num_double_write_disk = 0;
87 num_double_read_disk = 0;
88
89 the2DM = NULL;
90 the3DM = NULL;
91 theCorr = NULL;
92 Exc_activated = false;
93 makecheckpoints = makechkpt;
94 tempfolder = tmpfolder;
95
96 setupBookkeeperAndMPS( occupancies );
97 PreSolve();
98
99 }
100
setupBookkeeperAndMPS(int * occupancies)101 void CheMPS2::DMRG::setupBookkeeperAndMPS( int * occupancies ){
102
103 denBK = new SyBookkeeper( Prob, OptScheme->get_D( 0 ) );
104 assert( denBK->IsPossible() );
105
106 std::stringstream sstream;
107 sstream << CheMPS2::DMRG_MPS_storage_prefix << nStates-1 << ".h5";
108 MPSstoragename.assign( sstream.str() );
109 struct stat stFileInfo;
110 int intStat = stat( MPSstoragename.c_str(), &stFileInfo );
111 loadedMPS = (( makecheckpoints ) && ( intStat==0 )) ? true : false;
112 #ifdef CHEMPS2_MPI_COMPILATION
113 assert( MPIchemps2::all_booleans_equal( loadedMPS ) );
114 #endif
115
116 if ( loadedMPS ){ loadDIM( MPSstoragename, denBK ); }
117
118 // Convert occupancies from HAM to DMRG orbitals
119 if (( occupancies != NULL ) && ( Prob->gReorder() )){
120 int * tmp_cpy_occ = new int[ L ];
121 for ( int cnt = 0; cnt < L; cnt++ ){ tmp_cpy_occ[ cnt ] = occupancies[ cnt ]; }
122 for ( int cnt = 0; cnt < L; cnt++ ){ occupancies[ cnt ] = tmp_cpy_occ[ Prob->gf2( cnt ) ]; }
123 delete [] tmp_cpy_occ;
124 }
125
126 // Set to ROHF dimensions
127 /*if (( !loadedMPS ) && ( occupancies != NULL )){
128 int left_n = 0;
129 int left_i = 0;
130 int left_2s = 0;
131 for ( int site = 0; site < denBK->gL(); site++ ){
132 for ( int N = denBK->gNmin( site ); N <= denBK->gNmax( site ); N++ ){
133 for ( int TwoS = denBK->gTwoSmin( site, N ); TwoS <= denBK->gTwoSmax( site, N ); TwoS+=2 ){
134 for ( int Irrep = 0; Irrep < denBK->getNumberOfIrreps(); Irrep++ ){
135 denBK->SetDim( site, N, TwoS, Irrep, 0 );
136 }
137 }
138 }
139 denBK->SetDim( site, left_n, left_2s, left_i, 1 );
140 left_n = left_n + occupancies[ site ];
141 left_i = (( occupancies[ site ] == 1 ) ? Irreps::directProd( left_i, Prob->gIrrep( site ) ) : left_i );
142 left_2s = (( occupancies[ site ] == 1 ) ? ( left_2s + 1 ) : left_2s );
143 }
144 assert( left_n == Prob->gN() );
145 assert( left_i == Prob->gIrrep() );
146 assert( left_2s == Prob->gTwoS() );
147 }*/
148
149 MPS = new TensorT * [ L ];
150 for ( int cnt = 0; cnt < L; cnt++ ){ MPS[ cnt ] = new TensorT( cnt, denBK ); }
151
152 if ( loadedMPS ){
153 bool isConverged;
154 loadMPS( MPSstoragename, MPS, &isConverged );
155 #ifdef CHEMPS2_MPI_COMPILATION
156 if ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER )
157 #endif
158 { cout << "Loaded MPS " << MPSstoragename << " converged y/n? : " << isConverged << endl; }
159 } else {
160 #ifdef CHEMPS2_MPI_COMPILATION
161 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
162 #else
163 const bool am_i_master = true;
164 #endif
165 if ( occupancies == NULL ){
166 for ( int site = 0; site < L; site++ ){
167 if ( am_i_master ){ MPS[ site ]->random(); }
168 left_normalize( MPS[ site ], NULL );
169 }
170 } else {
171 assert( Prob->check_rohf_occ( occupancies ) ); // Check compatibility
172 int left_n = 0;
173 int left_i = 0;
174 int left_2s = 0;
175 for ( int site = 0; site < L; site++ ){
176 const int right_n = left_n + occupancies[ site ];
177 const int right_i = (( occupancies[ site ] == 1 ) ? Irreps::directProd( left_i, Prob->gIrrep( site ) ) : left_i );
178 const int right_2s = (( occupancies[ site ] == 1 ) ? ( left_2s + 1 ) : left_2s );
179 const int dimL = denBK->gCurrentDim( site, left_n, left_2s, left_i );
180 const int dimR = denBK->gCurrentDim( site + 1, right_n, right_2s, right_i );
181 assert( dimL > 0 );
182 assert( dimR > 0 );
183 if ( am_i_master ){
184 MPS[ site ]->random();
185 for ( int NL = right_n - 2; NL <= right_n; NL++ ){
186 const int DS = (( right_n == NL + 1 ) ? 1 : 0 );
187 const int IL = (( right_n == NL + 1 ) ? Irreps::directProd( right_i, Prob->gIrrep( site ) ) : right_i );
188 for ( int TwoSL = right_2s - DS; TwoSL <= right_2s + DS; TwoSL+=2 ){
189 const int dimL2 = denBK->gCurrentDim( site, NL, TwoSL, IL );
190 if ( dimL2 > 0 ){
191 double * space = MPS[ site ]->gStorage( NL, TwoSL, IL, right_n, right_2s, right_i );
192 for ( int row = 0; row < dimL2; row++ ){ space[ row + dimL2 * 0 ] = 0.0; }
193 if (( NL == left_n ) && ( TwoSL == left_2s ) && ( IL == left_i )){ space[ 0 + dimL * 0 ] = 42; }
194 }
195 }
196 }
197 }
198 left_normalize( MPS[ site ], NULL );
199 left_n = right_n;
200 left_i = right_i;
201 left_2s = right_2s;
202 }
203 assert( left_n == Prob->gN() );
204 assert( left_i == Prob->gIrrep() );
205 assert( left_2s == Prob->gTwoS() );
206 }
207 }
208
209 }
210
~DMRG()211 CheMPS2::DMRG::~DMRG(){
212
213 if ( the2DM != NULL ){ delete the2DM; }
214 if ( the3DM != NULL ){ delete the3DM; }
215 if ( theCorr != NULL ){ delete theCorr; }
216
217 deleteAllBoundaryOperators();
218
219 delete [] Ltensors;
220 delete [] F0tensors;
221 delete [] F1tensors;
222 delete [] S0tensors;
223 delete [] S1tensors;
224 delete [] Atensors;
225 delete [] Btensors;
226 delete [] Ctensors;
227 delete [] Dtensors;
228 delete [] Qtensors;
229 delete [] Xtensors;
230 delete [] isAllocated;
231
232 for ( int site = 0; site < L; site++ ){ delete MPS[ site ]; }
233 delete [] MPS;
234
235 if ( Exc_activated ){
236 delete [] Exc_Eshifts;
237 for ( int state = 0; state < nStates - 1; state++ ){
238 #ifdef CHEMPS2_MPI_COMPILATION
239 if ( MPIchemps2::owner_specific_excitation( L, state ) == MPIchemps2::mpi_rank() )
240 #endif
241 {
242 for ( int orb = 0; orb < L; orb++ ){ delete Exc_MPSs[ state ][ orb ]; }
243 delete [] Exc_MPSs[ state ];
244 delete Exc_BKs[ state ];
245 delete [] Exc_Overlaps[ state ]; // The rest is allocated and deleted at DMRGoperators.cpp
246 }
247 }
248 delete [] Exc_MPSs;
249 delete [] Exc_BKs;
250 delete [] Exc_Overlaps;
251 }
252
253 delete denBK;
254
255 }
256
PreSolve()257 void CheMPS2::DMRG::PreSolve(){
258
259 deleteAllBoundaryOperators();
260
261 for ( int cnt = 0; cnt < L - 2; cnt++ ){ updateMovingRightSafeFirstTime( cnt ); }
262
263 TotalMinEnergy = 1e8;
264 MaxDiscWeightLastSweep = 0.0;
265
266 }
267
Solve()268 double CheMPS2::DMRG::Solve(){
269
270 bool change = ( TotalMinEnergy < 1e8 ) ? true : false; // 1 sweep from right to left: fixed virtual dimensions
271
272 double Energy = 0.0;
273
274 #ifdef CHEMPS2_MPI_COMPILATION
275 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
276 #else
277 const bool am_i_master = true;
278 #endif
279
280 for ( int instruction = 0; instruction < OptScheme->get_number(); instruction++ ){
281
282 int nIterations = 0;
283 double EnergyPrevious = Energy + 10 * OptScheme->get_energy_conv( instruction ); // Guarantees that there's always at least 1 left-right sweep
284
285 while (( fabs( Energy - EnergyPrevious ) > OptScheme->get_energy_conv( instruction ) ) && ( nIterations < OptScheme->get_max_sweeps( instruction ) )){
286
287 for ( int timecnt = 0; timecnt < CHEMPS2_TIME_VECLENGTH; timecnt++ ){ timings[ timecnt ] = 0.0; }
288 num_double_write_disk = 0;
289 num_double_read_disk = 0;
290 struct timeval start, end;
291 EnergyPrevious = Energy;
292 gettimeofday( &start, NULL );
293 Energy = sweepleft( change, instruction, am_i_master ); // Only relevant call in this block of code
294 gettimeofday( &end, NULL );
295 double elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
296 if ( am_i_master ){
297 cout << "******************************************************************" << endl;
298 cout << "*** Information on left sweep " << nIterations << " of instruction " << instruction << ":" << endl;
299 cout << "*** Elapsed wall time = " << elapsed << " seconds" << endl;
300 cout << "*** |--> S.join = " << timings[ CHEMPS2_TIME_S_JOIN ] << " seconds" << endl;
301 cout << "*** |--> S.solve = " << timings[ CHEMPS2_TIME_S_SOLVE ] << " seconds" << endl;
302 cout << "*** |--> S.split = " << timings[ CHEMPS2_TIME_S_SPLIT ] << " seconds" << endl;
303 print_tensor_update_performance();
304 cout << "*** Minimum energy = " << LastMinEnergy << endl;
305 cout << "*** Maximum discarded weight = " << MaxDiscWeightLastSweep << endl;
306 }
307 if ( Exc_activated ){ calc_overlaps( false ); }
308 if ( am_i_master ){
309 cout << "******************************************************************" << endl;
310 }
311 change = true; //rest of sweeps: variable virtual dimensions
312 for ( int timecnt = 0; timecnt < CHEMPS2_TIME_VECLENGTH; timecnt++ ){ timings[ timecnt ] = 0.0; }
313 num_double_write_disk = 0;
314 num_double_read_disk = 0;
315 gettimeofday( &start, NULL );
316 Energy = sweepright( change, instruction, am_i_master ); // Only relevant call in this block of code
317 gettimeofday( &end, NULL );
318 elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
319 if ( am_i_master ){
320 cout << "******************************************************************" << endl;
321 cout << "*** Information on right sweep " << nIterations << " of instruction " << instruction << ":" << endl;
322 cout << "*** Elapsed wall time = " << elapsed << " seconds" << endl;
323 cout << "*** |--> S.join = " << timings[ CHEMPS2_TIME_S_JOIN ] << " seconds" << endl;
324 cout << "*** |--> S.solve = " << timings[ CHEMPS2_TIME_S_SOLVE ] << " seconds" << endl;
325 cout << "*** |--> S.split = " << timings[ CHEMPS2_TIME_S_SPLIT ] << " seconds" << endl;
326 print_tensor_update_performance();
327 cout << "*** Minimum energy = " << LastMinEnergy << endl;
328 cout << "*** Maximum discarded weight = " << MaxDiscWeightLastSweep << endl;
329 cout << "*** Energy difference with respect to previous leftright sweep = " << fabs(Energy-EnergyPrevious) << endl;
330 }
331 if ( Exc_activated ){ calc_overlaps( true ); }
332 if ( am_i_master ){
333 cout << "******************************************************************" << endl;
334 if ( makecheckpoints ){ saveMPS( MPSstoragename, MPS, denBK, false ); } // Only the master proc makes MPS checkpoints !!
335 }
336
337 nIterations++;
338
339 }
340
341 if ( am_i_master ){
342 cout << "*** Information on completed instruction " << instruction << ":" << endl;
343 cout << "*** The reduced virtual dimension DSU(2) = " << OptScheme->get_D(instruction) << endl;
344 cout << "*** The total number of reduced MPS variables = " << get_num_mps_var() << endl;
345 cout << "*** Minimum energy encountered during all instructions = " << TotalMinEnergy << endl;
346 cout << "*** Minimum energy encountered during the last sweep = " << LastMinEnergy << endl;
347 cout << "*** Maximum discarded weight during the last sweep = " << MaxDiscWeightLastSweep << endl;
348 cout << "******************************************************************" << endl;
349 }
350
351 }
352
353 return TotalMinEnergy;
354
355 }
356
sweepleft(const bool change,const int instruction,const bool am_i_master)357 double CheMPS2::DMRG::sweepleft( const bool change, const int instruction, const bool am_i_master ){
358
359 double Energy = 0.0;
360 const double noise_level = fabs( OptScheme->get_noise_prefactor( instruction ) ) * MaxDiscWeightLastSweep;
361 const double dvdson_rtol = OptScheme->get_dvdson_rtol( instruction );
362 const int vir_dimension = OptScheme->get_D( instruction );
363 MaxDiscWeightLastSweep = 0.0;
364 LastMinEnergy = 1e8;
365
366 for ( int index = L - 2; index > 0; index-- ){
367
368 Energy = solve_site( index, dvdson_rtol, noise_level, vir_dimension, am_i_master, false, change );
369 if ( Energy < TotalMinEnergy ){ TotalMinEnergy = Energy; }
370 if ( Energy < LastMinEnergy ){ LastMinEnergy = Energy; }
371 if ( am_i_master ){
372 cout << "Energy at sites (" << index << ", " << index + 1 << ") is " << Energy << endl;
373 }
374
375 // Prepare for next step
376 struct timeval start, end;
377 gettimeofday( &start, NULL );
378 updateMovingLeftSafe( index );
379 gettimeofday( &end, NULL );
380 timings[ CHEMPS2_TIME_TENS_TOTAL ] += ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
381
382 }
383
384 return Energy;
385
386 }
387
sweepright(const bool change,const int instruction,const bool am_i_master)388 double CheMPS2::DMRG::sweepright( const bool change, const int instruction, const bool am_i_master ){
389
390 double Energy = 0.0;
391 const double noise_level = fabs( OptScheme->get_noise_prefactor( instruction ) ) * MaxDiscWeightLastSweep;
392 const double dvdson_rtol = OptScheme->get_dvdson_rtol( instruction );
393 const int vir_dimension = OptScheme->get_D( instruction );
394 MaxDiscWeightLastSweep = 0.0;
395 LastMinEnergy = 1e8;
396
397 for ( int index = 0; index < L - 2; index++ ){
398
399 Energy = solve_site( index, dvdson_rtol, noise_level, vir_dimension, am_i_master, true, change );
400 if ( Energy < TotalMinEnergy ){ TotalMinEnergy = Energy; }
401 if ( Energy < LastMinEnergy ){ LastMinEnergy = Energy; }
402 if ( am_i_master ){
403 cout << "Energy at sites (" << index << ", " << index + 1 << ") is " << Energy << endl;
404 }
405
406 // Prepare for next step
407 struct timeval start, end;
408 gettimeofday( &start, NULL );
409 updateMovingRightSafe( index );
410 gettimeofday( &end, NULL );
411 timings[ CHEMPS2_TIME_TENS_TOTAL ] += ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
412
413 }
414
415 return Energy;
416
417 }
418
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)419 double CheMPS2::DMRG::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 ){
420
421 struct timeval start, end;
422
423 // Construct two-site object S. Each MPI process joins the MPS tensors. Before a matrix-vector multiplication the vector is broadcasted anyway.
424 gettimeofday( &start, NULL );
425 Sobject * denS = new Sobject( index, denBK );
426 denS->Join( MPS[ index ], MPS[ index + 1 ] );
427 gettimeofday( &end, NULL );
428 timings[ CHEMPS2_TIME_S_JOIN ] += ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
429
430 // Feed everything to the solver. Each MPI process returns the correct energy. Only MPI_CHEMPS2_MASTER has the correct denS solution.
431 gettimeofday( &start, NULL );
432 Heff Solver( denBK, Prob, dvdson_rtol );
433 double ** VeffTilde = NULL;
434 if ( Exc_activated ){ VeffTilde = prepare_excitations( denS ); }
435 double Energy = Solver.SolveDAVIDSON( denS, Ltensors, Atensors, Btensors, Ctensors, Dtensors, S0tensors, S1tensors, F0tensors, F1tensors, Qtensors, Xtensors, nStates - 1, VeffTilde );
436 Energy += Prob->gEconst();
437 if ( Exc_activated ){ cleanup_excitations( VeffTilde ); }
438 gettimeofday( &end, NULL );
439 timings[ CHEMPS2_TIME_S_SOLVE ] += ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
440
441 // Decompose the S-object. MPI_CHEMPS2_MASTER decomposes denS. Each MPI process returns the correct discWeight. Each MPI process has the new MPS tensors set.
442 gettimeofday( &start, NULL );
443 if (( noise_level > 0.0 ) && ( am_i_master )){ denS->addNoise( noise_level ); }
444 const double discWeight = denS->Split( MPS[ index ], MPS[ index + 1 ], virtual_dimension, moving_right, change );
445 delete denS;
446 if ( discWeight > MaxDiscWeightLastSweep ){ MaxDiscWeightLastSweep = discWeight; }
447 gettimeofday( &end, NULL );
448 timings[ CHEMPS2_TIME_S_SPLIT ] += ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
449
450 return Energy;
451
452 }
453
get_num_mps_var() const454 int CheMPS2::DMRG::get_num_mps_var() const{
455
456 int num_var = 0;
457 for ( int site = 0; site < L; site++ ){
458 num_var += MPS[ site ]->gKappa2index( MPS[ site ]->gNKappa() );
459 }
460 return num_var;
461
462 }
463
activateExcitations(const int maxExcIn)464 void CheMPS2::DMRG::activateExcitations( const int maxExcIn ){
465
466 Exc_activated = true;
467 maxExc = maxExcIn;
468 Exc_Eshifts = new double[ maxExc ];
469 Exc_MPSs = new TensorT ** [ maxExc ];
470 Exc_BKs = new SyBookkeeper * [ maxExc ];
471 Exc_Overlaps = new TensorO ** [ maxExc ];
472
473 }
474
newExcitation(const double EshiftIn)475 void CheMPS2::DMRG::newExcitation( const double EshiftIn ){
476
477 assert( Exc_activated );
478 assert( nStates - 1 < maxExc );
479
480 if ( the2DM != NULL ){ delete the2DM; the2DM = NULL; }
481 if ( the3DM != NULL ){ delete the3DM; the3DM = NULL; }
482 if ( theCorr != NULL ){ delete theCorr; theCorr = NULL; }
483 deleteAllBoundaryOperators();
484
485 Exc_Eshifts[ nStates - 1 ] = EshiftIn;
486 #ifdef CHEMPS2_MPI_COMPILATION
487 if ( MPIchemps2::owner_specific_excitation( L, nStates - 1 ) == MPIchemps2::mpi_rank() ){
488 #endif
489 Exc_MPSs[ nStates - 1 ] = MPS;
490 Exc_BKs[ nStates - 1 ] = denBK;
491 Exc_Overlaps[ nStates - 1 ] = new TensorO*[ L - 1 ];
492 #ifdef CHEMPS2_MPI_COMPILATION
493 } else {
494 for ( int site = 0; site < L; site++ ){ delete MPS[ site ]; }
495 delete [] MPS;
496 delete denBK;
497 }
498 #endif
499
500 nStates++;
501
502 setupBookkeeperAndMPS();
503 PreSolve();
504
505 }
506
507