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 <algorithm>
23 #include <assert.h>
24
25 #include "Sobject.h"
26 #include "TensorT.h"
27 #include "SyBookkeeper.h"
28 #include "Lapack.h"
29 #include "MPIchemps2.h"
30 #include "Wigner.h"
31 #include "Special.h"
32
33 using std::min;
34 using std::max;
35
Sobject(const int index,SyBookkeeper * denBK)36 CheMPS2::Sobject::Sobject( const int index, SyBookkeeper * denBK ){
37
38 this->index = index;
39 this->denBK = denBK;
40 Ilocal1 = denBK->gIrrep( index );
41 Ilocal2 = denBK->gIrrep( index + 1 );
42
43 nKappa = 0;
44
45 for ( int NL = denBK->gNmin( index ); NL <= denBK->gNmax( index ); NL++ ){
46 for ( int TwoSL = denBK->gTwoSmin( index, NL ); TwoSL <= denBK->gTwoSmax( index, NL ); TwoSL += 2 ){
47 for ( int IL = 0; IL < denBK->getNumberOfIrreps(); IL++ ){
48 const int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL );
49 if ( dimL > 0 ){
50 for ( int N1 = 0; N1 <= 2; N1++ ){
51 for ( int N2 = 0; N2 <= 2; N2++ ){
52 const int NR = NL + N1 + N2;
53 const int IM = (( N1 == 1 ) ? Irreps::directProd( IL, Ilocal1 ) : IL );
54 const int IR = (( N2 == 1 ) ? Irreps::directProd( IM, Ilocal2 ) : IM );
55 const int TwoJmin = ( N1 + N2 ) % 2;
56 const int TwoJmax = ((( N1 == 1 ) && ( N2 == 1 )) ? 2 : TwoJmin );
57 for ( int TwoJ = TwoJmin; TwoJ <= TwoJmax; TwoJ += 2 ){
58 for ( int TwoSR = TwoSL - TwoJ; TwoSR <= TwoSL + TwoJ; TwoSR += 2 ){
59 if ( TwoSR >= 0 ){
60 const int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR );
61 if ( dimR > 0 ){
62 nKappa++;
63 }
64 }
65 }
66 }
67 }
68 }
69 }
70 }
71 }
72 }
73
74 sectorNL = new int[ nKappa ];
75 sectorTwoSL = new int[ nKappa ];
76 sectorIL = new int[ nKappa ];
77 sectorN1 = new int[ nKappa ];
78 sectorN2 = new int[ nKappa ];
79 sectorTwoJ = new int[ nKappa ];
80 sectorNR = new int[ nKappa ];
81 sectorTwoSR = new int[ nKappa ];
82 sectorIR = new int[ nKappa ];
83 kappa2index = new int[ nKappa + 1 ];
84 kappa2index[ 0 ] = 0;
85
86 nKappa = 0;
87
88 for ( int NL = denBK->gNmin( index ); NL <= denBK->gNmax( index ); NL++ ){
89 for ( int TwoSL = denBK->gTwoSmin( index, NL ); TwoSL <= denBK->gTwoSmax( index, NL ); TwoSL += 2 ){
90 for ( int IL = 0; IL < denBK->getNumberOfIrreps(); IL++ ){
91 const int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL );
92 if ( dimL > 0 ){
93 for ( int N1 = 0; N1 <= 2; N1++ ){
94 for ( int N2 = 0; N2 <= 2; N2++ ){
95 const int NR = NL + N1 + N2;
96 const int IM = (( N1 == 1 ) ? Irreps::directProd( IL, Ilocal1 ) : IL );
97 const int IR = (( N2 == 1 ) ? Irreps::directProd( IM, Ilocal2 ) : IM );
98 const int TwoJmin = ( N1 + N2 ) % 2;
99 const int TwoJmax = ((( N1 == 1 ) && ( N2 == 1 )) ? 2 : TwoJmin );
100 for ( int TwoJ = TwoJmin; TwoJ <= TwoJmax; TwoJ += 2 ){
101 for ( int TwoSR = TwoSL - TwoJ; TwoSR <= TwoSL + TwoJ; TwoSR += 2 ){
102 if ( TwoSR >= 0 ){
103 const int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR );
104 if ( dimR > 0 ){
105 sectorNL [ nKappa ] = NL;
106 sectorTwoSL[ nKappa ] = TwoSL;
107 sectorIL [ nKappa ] = IL;
108 sectorN1 [ nKappa ] = N1;
109 sectorN2 [ nKappa ] = N2;
110 sectorTwoJ [ nKappa ] = TwoJ;
111 sectorNR [ nKappa ] = NR;
112 sectorTwoSR[ nKappa ] = TwoSR;
113 sectorIR [ nKappa ] = IR;
114 nKappa++;
115 kappa2index[ nKappa ] = kappa2index[ nKappa - 1 ] + dimL * dimR;
116 }
117 }
118 }
119 }
120 }
121 }
122 }
123 }
124 }
125 }
126
127 storage = new double[ kappa2index[ nKappa ] ];
128
129 reorder = new int[ nKappa ];
130 for ( int cnt = 0; cnt < nKappa; cnt++ ){ reorder[ cnt ] = cnt; }
131 bool sorted = false;
132 while ( sorted == false ){ //Bubble sort so that blocksize(reorder[i]) >= blocksize(reorder[i+1]), with blocksize(k) = kappa2index[k+1]-kappa2index[k]
133 sorted = true;
134 for ( int cnt = 0; cnt < nKappa - 1; cnt++ ){
135 const int index1 = reorder[ cnt ];
136 const int index2 = reorder[ cnt + 1 ];
137 const int size1 = kappa2index[ index1 + 1 ] - kappa2index[ index1 ];
138 const int size2 = kappa2index[ index2 + 1 ] - kappa2index[ index2 ];
139 if ( size1 < size2 ){
140 sorted = false;
141 reorder[ cnt ] = index2;
142 reorder[ cnt + 1 ] = index1;
143 }
144 }
145 }
146
147 }
148
~Sobject()149 CheMPS2::Sobject::~Sobject(){
150
151 delete [] sectorNL;
152 delete [] sectorTwoSL;
153 delete [] sectorIL;
154 delete [] sectorN1;
155 delete [] sectorN2;
156 delete [] sectorTwoJ;
157 delete [] sectorNR;
158 delete [] sectorTwoSR;
159 delete [] sectorIR;
160 delete [] kappa2index;
161 delete [] storage;
162 delete [] reorder;
163
164 }
165
gNKappa() const166 int CheMPS2::Sobject::gNKappa() const { return nKappa; }
167
gStorage()168 double * CheMPS2::Sobject::gStorage() { return storage; }
169
gReorder(const int ikappa) const170 int CheMPS2::Sobject::gReorder( const int ikappa ) const{ return reorder[ ikappa ]; }
171
gKappa(const int NL,const int TwoSL,const int IL,const int N1,const int N2,const int TwoJ,const int NR,const int TwoSR,const int IR) const172 int CheMPS2::Sobject::gKappa( const int NL, const int TwoSL, const int IL, const int N1, const int N2, const int TwoJ, const int NR, const int TwoSR, const int IR ) const{
173
174 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
175 if (( sectorNL [ ikappa ] == NL ) &&
176 ( sectorTwoSL[ ikappa ] == TwoSL ) &&
177 ( sectorIL [ ikappa ] == IL ) &&
178 ( sectorN1 [ ikappa ] == N1 ) &&
179 ( sectorN2 [ ikappa ] == N2 ) &&
180 ( sectorTwoJ [ ikappa ] == TwoJ ) &&
181 ( sectorNR [ ikappa ] == NR ) &&
182 ( sectorTwoSR[ ikappa ] == TwoSR ) &&
183 ( sectorIR [ ikappa ] == IR )){ return ikappa; }
184 }
185
186 return -1;
187
188 }
189
gKappa2index(const int kappa) const190 int CheMPS2::Sobject::gKappa2index( const int kappa ) const{ return kappa2index[ kappa ]; }
191
gStorage(const int NL,const int TwoSL,const int IL,const int N1,const int N2,const int TwoJ,const int NR,const int TwoSR,const int IR)192 double * CheMPS2::Sobject::gStorage( const int NL, const int TwoSL, const int IL, const int N1, const int N2, const int TwoJ, const int NR, const int TwoSR, const int IR ){
193
194 const int kappa = gKappa( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
195 if ( kappa == -1 ){ return NULL; }
196 return storage + kappa2index[ kappa ];
197
198 }
199
gIndex() const200 int CheMPS2::Sobject::gIndex() const { return index; }
201
gNL(const int ikappa) const202 int CheMPS2::Sobject::gNL ( const int ikappa ) const{ return sectorNL [ ikappa ]; }
gTwoSL(const int ikappa) const203 int CheMPS2::Sobject::gTwoSL( const int ikappa ) const{ return sectorTwoSL[ ikappa ]; }
gIL(const int ikappa) const204 int CheMPS2::Sobject::gIL ( const int ikappa ) const{ return sectorIL [ ikappa ]; }
gN1(const int ikappa) const205 int CheMPS2::Sobject::gN1 ( const int ikappa ) const{ return sectorN1 [ ikappa ]; }
gN2(const int ikappa) const206 int CheMPS2::Sobject::gN2 ( const int ikappa ) const{ return sectorN2 [ ikappa ]; }
gTwoJ(const int ikappa) const207 int CheMPS2::Sobject::gTwoJ ( const int ikappa ) const{ return sectorTwoJ [ ikappa ]; }
gNR(const int ikappa) const208 int CheMPS2::Sobject::gNR ( const int ikappa ) const{ return sectorNR [ ikappa ]; }
gTwoSR(const int ikappa) const209 int CheMPS2::Sobject::gTwoSR( const int ikappa ) const{ return sectorTwoSR[ ikappa ]; }
gIR(const int ikappa) const210 int CheMPS2::Sobject::gIR ( const int ikappa ) const{ return sectorIR [ ikappa ]; }
211
Join(TensorT * Tleft,TensorT * Tright)212 void CheMPS2::Sobject::Join( TensorT * Tleft, TensorT * Tright ){
213
214 #pragma omp parallel for schedule(dynamic)
215 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
216
217 const int NL = sectorNL [ ikappa ];
218 const int TwoSL = sectorTwoSL[ ikappa ];
219 const int IL = sectorIL [ ikappa ];
220
221 const int NR = sectorNR [ ikappa ];
222 const int TwoSR = sectorTwoSR[ ikappa ];
223 const int IR = sectorIR [ ikappa ];
224
225 const int TwoJ = sectorTwoJ [ ikappa ];
226 const int N1 = sectorN1 [ ikappa ];
227 const int N2 = sectorN2 [ ikappa ];
228 const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
229 const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
230 const int fase = Special::phase( TwoSL + TwoSR + TwoS1 + TwoS2 );
231
232 // Clear block
233 int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL ); // dimL > 0, checked at creation
234 int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR ); // dimR > 0, checked at creation
235 double * block_s = storage + kappa2index[ ikappa ];
236 for ( int cnt = 0; cnt < dimL * dimR; cnt++ ){ block_s[ cnt ] = 0.0; }
237
238 // Central symmetry sectors
239 const int NM = NL + N1;
240 const int IM = (( TwoS1 == 1 ) ? Irreps::directProd( IL, Ilocal1 ) : IL );
241 const int TwoJMlower = max( abs( TwoSL - TwoS1 ), abs( TwoSR - TwoS2 ) );
242 const int TwoJMupper = min( ( TwoSL + TwoS1 ), ( TwoSR + TwoS2 ) );
243 for ( int TwoJM = TwoJMlower; TwoJM <= TwoJMupper; TwoJM += 2 ){
244 int dimM = denBK->gCurrentDim( index + 1, NM, TwoJM, IM );
245 if ( dimM > 0 ){
246 double * block_left = Tleft ->gStorage( NL, TwoSL, IL, NM, TwoJM, IM );
247 double * block_right = Tright->gStorage( NM, TwoJM, IM, NR, TwoSR, IR );
248 double prefactor = fase
249 * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoJM + 1 ) )
250 * Wigner::wigner6j( TwoSL, TwoSR, TwoJ, TwoS2, TwoS1, TwoJM );
251 char notrans = 'N';
252 double add = 1.0;
253 dgemm_( ¬rans, ¬rans, &dimL, &dimR, &dimM, &prefactor, block_left, &dimL, block_right, &dimM, &add, block_s, &dimL );
254 }
255 }
256 }
257
258 }
259
Split(TensorT * Tleft,TensorT * Tright,const int virtualdimensionD,const bool movingright,const bool change)260 double CheMPS2::Sobject::Split( TensorT * Tleft, TensorT * Tright, const int virtualdimensionD, const bool movingright, const bool change ){
261
262 #ifdef CHEMPS2_MPI_COMPILATION
263 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
264 #endif
265
266 // Get the number of central sectors
267 int nCenterSectors = 0;
268 for ( int NM = denBK->gNmin( index + 1 ); NM <= denBK->gNmax( index + 1 ); NM++ ){
269 for ( int TwoSM = denBK->gTwoSmin( index + 1, NM ); TwoSM <= denBK->gTwoSmax( index + 1, NM ); TwoSM += 2 ){
270 for ( int IM = 0; IM < denBK->getNumberOfIrreps(); IM++ ){
271 const int dimM = denBK->gFCIdim( index + 1, NM, TwoSM, IM ); //FCIdim !! Whether possible hence.
272 if ( dimM > 0 ){
273 nCenterSectors++;
274 }
275 }
276 }
277 }
278
279 // Get the labels of the central sectors
280 int * SplitSectNM = new int[ nCenterSectors ];
281 int * SplitSectTwoJM = new int[ nCenterSectors ];
282 int * SplitSectIM = new int[ nCenterSectors ];
283 nCenterSectors = 0;
284 for ( int NM = denBK->gNmin( index + 1 ); NM <= denBK->gNmax( index + 1 ); NM++ ){
285 for ( int TwoSM = denBK->gTwoSmin( index + 1, NM ); TwoSM <= denBK->gTwoSmax( index + 1, NM ); TwoSM += 2 ){
286 for ( int IM = 0; IM < denBK->getNumberOfIrreps(); IM++ ){
287 const int dimM = denBK->gFCIdim( index + 1, NM, TwoSM, IM ); //FCIdim !! Whether possible hence.
288 if ( dimM > 0 ){
289 SplitSectNM [ nCenterSectors ] = NM;
290 SplitSectTwoJM[ nCenterSectors ] = TwoSM;
291 SplitSectIM [ nCenterSectors ] = IM;
292 nCenterSectors++;
293 }
294 }
295 }
296 }
297
298 // Only MPI_CHEMPS2_MASTER performs SVD --> Allocate memory
299 double ** Lambdas = NULL;
300 double ** Us = NULL;
301 double ** VTs = NULL;
302 int * CenterDims = NULL;
303 int * DimLtotal = NULL;
304 int * DimRtotal = NULL;
305
306 #ifdef CHEMPS2_MPI_COMPILATION
307 if ( am_i_master ){
308 #endif
309
310 Lambdas = new double*[ nCenterSectors ];
311 Us = new double*[ nCenterSectors ];
312 VTs = new double*[ nCenterSectors ];
313 CenterDims = new int[ nCenterSectors ];
314 DimLtotal = new int[ nCenterSectors ];
315 DimRtotal = new int[ nCenterSectors ];
316
317 //PARALLEL
318 #pragma omp parallel for schedule(dynamic)
319 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
320
321 //Determine left and right dimensions contributing to the center block iCenter
322 DimLtotal[ iCenter ] = 0;
323 for ( int NL = SplitSectNM[ iCenter ] - 2; NL <= SplitSectNM[ iCenter ]; NL++ ){
324 const int TwoS1 = (( NL + 1 == SplitSectNM[ iCenter ] ) ? 1 : 0 );
325 for ( int TwoSL = SplitSectTwoJM[ iCenter ] - TwoS1; TwoSL <= SplitSectTwoJM[ iCenter ] + TwoS1; TwoSL += 2 ){
326 if ( TwoSL >= 0 ){
327 const int IL = (( TwoS1 == 1 ) ? Irreps::directProd( Ilocal1, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
328 const int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL );
329 if ( dimL > 0 ){
330 DimLtotal[ iCenter ] += dimL;
331 }
332 }
333 }
334 }
335 DimRtotal[ iCenter ] = 0;
336 for ( int NR = SplitSectNM[ iCenter ]; NR <= SplitSectNM[ iCenter ] + 2; NR++ ){
337 const int TwoS2 = (( NR == SplitSectNM[ iCenter ] + 1 ) ? 1 : 0 );
338 for ( int TwoSR = SplitSectTwoJM[ iCenter ] - TwoS2; TwoSR <= SplitSectTwoJM[ iCenter ] + TwoS2; TwoSR += 2 ){
339 if ( TwoSR >= 0 ){
340 const int IR = (( TwoS2 == 1 ) ? Irreps::directProd( Ilocal2, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
341 const int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR );
342 if ( dimR > 0 ){
343 DimRtotal[ iCenter ] += dimR;
344 }
345 }
346 }
347 }
348 CenterDims[ iCenter ] = min( DimLtotal[ iCenter ], DimRtotal[ iCenter ] ); // CenterDims contains the min. amount
349
350 //Allocate memory to copy the different parts of the S-object. Use prefactor sqrt((2jR+1)/(2jM+1) * (2jM+1) * (2j+1)) W6J (-1)^(jL+jR+s1+s2) and sum over j.
351 if ( CenterDims[ iCenter ] > 0 ){
352
353 // Only if CenterDims[ iCenter ] exists should you allocate the following three arrays
354 Lambdas[ iCenter ] = new double[ CenterDims[ iCenter ] ];
355 Us[ iCenter ] = new double[ CenterDims[ iCenter ] * DimLtotal[ iCenter ] ];
356 VTs[ iCenter ] = new double[ CenterDims[ iCenter ] * DimRtotal[ iCenter ] ];
357
358 const int memsize = DimLtotal[ iCenter ] * DimRtotal[ iCenter ];
359 double * mem = new double[ memsize ];
360 for ( int cnt = 0; cnt < memsize; cnt++ ){ mem[ cnt ] = 0.0; }
361
362 int dimLtotal2 = 0;
363 for ( int NL = SplitSectNM[ iCenter ] - 2; NL <= SplitSectNM[ iCenter ]; NL++ ){
364 const int TwoS1 = (( NL + 1 == SplitSectNM[ iCenter ] ) ? 1 : 0 );
365 for ( int TwoSL = SplitSectTwoJM[ iCenter ] - TwoS1; TwoSL <= SplitSectTwoJM[ iCenter ] + TwoS1; TwoSL += 2 ){
366 if ( TwoSL >= 0 ){
367 const int IL = (( TwoS1 == 1 ) ? Irreps::directProd( Ilocal1, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
368 const int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL );
369 if ( dimL > 0 ){
370 int dimRtotal2 = 0;
371 for ( int NR = SplitSectNM[ iCenter ]; NR <= SplitSectNM[ iCenter ] + 2; NR++ ){
372 const int TwoS2 = (( NR == SplitSectNM[ iCenter ] + 1 ) ? 1 : 0 );
373 for ( int TwoSR = SplitSectTwoJM[ iCenter ] - TwoS2; TwoSR <= SplitSectTwoJM[ iCenter ] + TwoS2; TwoSR += 2 ){
374 if ( TwoSR >= 0 ){
375 const int IR = (( TwoS2 == 1 ) ? Irreps::directProd( Ilocal2, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
376 const int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR );
377 if ( dimR > 0 ){
378 // Loop over contributing TwoJ's
379 const int fase = Special::phase( TwoSL + TwoSR + TwoS1 + TwoS2 );
380 const int TwoJmin = max( abs( TwoSR - TwoSL ), abs( TwoS2 - TwoS1 ) );
381 const int TwoJmax = min( TwoS1 + TwoS2, TwoSL + TwoSR );
382 for ( int TwoJ = TwoJmin; TwoJ <= TwoJmax; TwoJ += 2 ){
383 // Calc prefactor
384 const double prefactor = fase
385 * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoSR + 1 ) )
386 * Wigner::wigner6j( TwoSL, TwoSR, TwoJ, TwoS2, TwoS1, SplitSectTwoJM[ iCenter ] );
387
388 // Add them to mem --> += because several TwoJ
389 double * Block = gStorage( NL, TwoSL, IL, SplitSectNM[ iCenter ] - NL, NR - SplitSectNM[ iCenter ], TwoJ, NR, TwoSR, IR );
390 for ( int l = 0; l < dimL; l++ ){
391 for ( int r = 0; r < dimR; r++ ){
392 mem[ dimLtotal2 + l + DimLtotal[ iCenter ] * ( dimRtotal2 + r ) ] += prefactor * Block[ l + dimL * r ];
393 }
394 }
395 }
396 dimRtotal2 += dimR;
397 }
398 }
399 }
400 }
401 dimLtotal2 += dimL;
402 }
403 }
404 }
405 }
406
407 // Now mem contains sqrt((2jR+1)/(2jM+1)) * (TT)^{jM nM IM) --> SVD per central symmetry
408 char jobz = 'S'; // M x min(M,N) in U and min(M,N) x N in VT
409 int lwork = 3 * CenterDims[ iCenter ] + max( max( DimLtotal[ iCenter ], DimRtotal[ iCenter ] ), 4 * CenterDims[ iCenter ] * ( CenterDims[ iCenter ] + 1 ) );
410 double * work = new double[ lwork ];
411 int * iwork = new int[ 8 * CenterDims[ iCenter ] ];
412 int info;
413
414 // dgesdd is not thread-safe in every implementation ( intel MKL is safe, Atlas is not safe )
415 #ifndef CHEMPS2_MKL
416 #pragma omp critical
417 #endif
418 dgesdd_( &jobz, DimLtotal + iCenter, DimRtotal + iCenter, mem, DimLtotal + iCenter,
419 Lambdas[ iCenter ], Us[ iCenter ], DimLtotal + iCenter, VTs[ iCenter ], CenterDims + iCenter, work, &lwork, iwork, &info );
420
421 delete [] work;
422 delete [] iwork;
423 delete [] mem;
424 }
425 }
426
427 #ifdef CHEMPS2_MPI_COMPILATION
428 }
429 #endif
430
431 double discardedWeight = 0.0; // Only if change==true; will the discardedWeight be meaningful and different from zero.
432 int updateSectors = 0;
433 int * NewDims = NULL;
434
435 // If change: determine new virtual dimensions.
436 #ifdef CHEMPS2_MPI_COMPILATION
437 if (( change ) && ( am_i_master )){
438 #else
439 if ( change ){
440 #endif
441
442 NewDims = new int[ nCenterSectors ];
443 // First determine the total number of singular values
444 int totalDimSVD = 0;
445 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
446 NewDims[ iCenter ] = CenterDims[ iCenter ];
447 totalDimSVD += NewDims[ iCenter ];
448 }
449
450 // If larger then the required virtualdimensionD, new virtual dimensions will be set in NewDims.
451 if ( totalDimSVD > virtualdimensionD ){
452 // Copy them all in 1 array
453 double * values = new double[ totalDimSVD ];
454 totalDimSVD = 0;
455 int inc = 1;
456 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
457 if ( NewDims[ iCenter ] > 0 ){
458 dcopy_( NewDims + iCenter, Lambdas[ iCenter ], &inc, values + totalDimSVD, &inc );
459 totalDimSVD += NewDims[ iCenter ];
460 }
461 }
462
463 // Sort them in decreasing order
464 char ID = 'D';
465 int info;
466 dlasrt_( &ID, &totalDimSVD, values, &info ); // Quicksort
467
468 // The D+1'th value becomes the lower bound Schmidt value. Every value smaller than or equal to the D+1'th value is thrown out (hence Dactual <= Ddesired).
469 const double lowerBound = values[ virtualdimensionD ];
470 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
471 for ( int cnt = 0; cnt < NewDims[ iCenter ]; cnt++ ){
472 if ( Lambdas[ iCenter ][ cnt ] <= lowerBound ){ NewDims[ iCenter ] = cnt; }
473 }
474 }
475
476 // Discarded weight
477 double totalSum = 0.0;
478 double discardedSum = 0.0;
479 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
480 for ( int iLocal = 0; iLocal < CenterDims[ iCenter ]; iLocal++ ){
481 double temp = ( SplitSectTwoJM[ iCenter ] + 1 ) * Lambdas[ iCenter ][ iLocal ] * Lambdas[ iCenter ][ iLocal ];
482 totalSum += temp;
483 if ( Lambdas[ iCenter ][ iLocal ] <= lowerBound ){ discardedSum += temp; }
484 }
485 }
486 discardedWeight = discardedSum / totalSum;
487
488 // Clean-up
489 delete [] values;
490 }
491
492 // Check if there is a sector which differs
493 updateSectors = 0;
494 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
495 const int MPSdim = denBK->gCurrentDim( index + 1, SplitSectNM[ iCenter ], SplitSectTwoJM[ iCenter ], SplitSectIM[ iCenter ] );
496 if ( NewDims[ iCenter ] != MPSdim ){ updateSectors = 1; }
497 }
498
499 }
500
501 #ifdef CHEMPS2_MPI_COMPILATION
502 MPIchemps2::broadcast_array_int( &updateSectors, 1, MPI_CHEMPS2_MASTER );
503 MPIchemps2::broadcast_array_double( &discardedWeight, 1, MPI_CHEMPS2_MASTER );
504 #endif
505
506 if ( updateSectors == 1 ){
507
508 #ifdef CHEMPS2_MPI_COMPILATION
509 if ( NewDims == NULL ){ NewDims = new int[ nCenterSectors ]; }
510 MPIchemps2::broadcast_array_int( NewDims, nCenterSectors, MPI_CHEMPS2_MASTER );
511 #endif
512
513 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
514 denBK->SetDim( index + 1, SplitSectNM[ iCenter ], SplitSectTwoJM[ iCenter ], SplitSectIM[ iCenter ], NewDims[ iCenter ] );
515 }
516 Tleft ->Reset();
517 Tright->Reset();
518
519 }
520
521 if ( NewDims != NULL ){ delete [] NewDims; }
522
523 #ifdef CHEMPS2_MPI_COMPILATION
524 if ( am_i_master ){
525 #endif
526
527 // Copy first dimM per central symmetry sector to the relevant parts
528 #pragma omp parallel for schedule(dynamic)
529 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
530 const int dimM = denBK->gCurrentDim( index + 1, SplitSectNM[ iCenter ], SplitSectTwoJM[ iCenter ], SplitSectIM[ iCenter ] );
531 if ( dimM > 0 ){
532 // U-part: copy
533 int dimLtotal2 = 0;
534 for ( int NL = SplitSectNM[ iCenter ] - 2; NL <= SplitSectNM[ iCenter ]; NL++ ){
535 const int TwoS1 = (( NL + 1 == SplitSectNM[ iCenter ] ) ? 1 : 0 );
536 for ( int TwoSL = SplitSectTwoJM[ iCenter ] - TwoS1; TwoSL <= SplitSectTwoJM[ iCenter ] + TwoS1; TwoSL += 2 ){
537 if ( TwoSL >= 0 ){
538 const int IL = (( TwoS1 == 1 ) ? Irreps::directProd( Ilocal1, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
539 const int dimL = denBK->gCurrentDim( index, NL, TwoSL, IL );
540 if ( dimL > 0 ){
541 double * TleftBlock = Tleft->gStorage( NL, TwoSL, IL, SplitSectNM[ iCenter ], SplitSectTwoJM[ iCenter ], SplitSectIM[ iCenter ] );
542 const int dimension_limit_right = min( dimM, CenterDims[ iCenter ] );
543 for ( int r = 0; r < dimension_limit_right; r++ ){
544 const double factor = (( movingright ) ? 1.0 : Lambdas[ iCenter ][ r ] );
545 for ( int l = 0; l < dimL; l++ ){
546 TleftBlock[ l + dimL * r ] = factor * Us[ iCenter ][ dimLtotal2 + l + DimLtotal[ iCenter ] * r ];
547 }
548 }
549 for ( int r = dimension_limit_right; r < dimM; r++ ){
550 for ( int l = 0; l < dimL; l++ ){
551 TleftBlock[ l + dimL * r ] = 0.0;
552 }
553 }
554 dimLtotal2 += dimL;
555 }
556 }
557 }
558 }
559
560 // VT-part: copy
561 int dimRtotal2 = 0;
562 for ( int NR = SplitSectNM[ iCenter ]; NR <= SplitSectNM[ iCenter ] + 2; NR++ ){
563 const int TwoS2 = (( NR == SplitSectNM[ iCenter ] + 1 ) ? 1 : 0 );
564 for ( int TwoSR = SplitSectTwoJM[ iCenter ] - TwoS2; TwoSR <= SplitSectTwoJM[ iCenter ] + TwoS2; TwoSR += 2 ){
565 if ( TwoSR >= 0 ){
566 const int IR = (( TwoS2 == 1 ) ? Irreps::directProd( Ilocal2, SplitSectIM[ iCenter ] ) : SplitSectIM[ iCenter ] );
567 const int dimR = denBK->gCurrentDim( index + 2, NR, TwoSR, IR );
568 if ( dimR > 0 ){
569 double * TrightBlock = Tright->gStorage( SplitSectNM[ iCenter ], SplitSectTwoJM[ iCenter ], SplitSectIM[ iCenter ], NR, TwoSR, IR );
570 const int dimension_limit_left = min( dimM, CenterDims[ iCenter ] );
571 const double factor_base = sqrt( ( SplitSectTwoJM[ iCenter ] + 1.0 ) / ( TwoSR + 1 ) );
572 for ( int l = 0; l < dimension_limit_left; l++ ){
573 const double factor = factor_base * (( movingright ) ? Lambdas[ iCenter ][ l ] : 1.0 );
574 for ( int r = 0; r < dimR; r++ ){
575 TrightBlock[ l + dimM * r ] = factor * VTs[ iCenter ][ l + CenterDims[ iCenter ] * ( dimRtotal2 + r ) ];
576 }
577 }
578 for ( int r = 0; r < dimR; r++ ){
579 for ( int l = dimension_limit_left; l < dimM; l++ ){
580 TrightBlock[ l + dimM * r ] = 0.0;
581 }
582 }
583 dimRtotal2 += dimR;
584 }
585 }
586 }
587 }
588 }
589 }
590
591 #ifdef CHEMPS2_MPI_COMPILATION
592 }
593 MPIchemps2::broadcast_tensor( Tleft, MPI_CHEMPS2_MASTER );
594 MPIchemps2::broadcast_tensor( Tright, MPI_CHEMPS2_MASTER );
595 #endif
596
597 // Clean up
598 delete [] SplitSectNM;
599 delete [] SplitSectTwoJM;
600 delete [] SplitSectIM;
601 #ifdef CHEMPS2_MPI_COMPILATION
602 if ( am_i_master )
603 #endif
604 {
605 for ( int iCenter = 0; iCenter < nCenterSectors; iCenter++ ){
606 if ( CenterDims[ iCenter ] > 0 ){
607 delete [] Us[ iCenter ];
608 delete [] Lambdas[ iCenter ];
609 delete [] VTs[ iCenter ];
610 }
611 }
612 delete [] Us;
613 delete [] Lambdas;
614 delete [] VTs;
615 delete [] CenterDims;
616 delete [] DimLtotal;
617 delete [] DimRtotal;
618 }
619
620 return discardedWeight;
621
622 }
623
624 void CheMPS2::Sobject::prog2symm(){
625
626 #pragma omp parallel for schedule(dynamic)
627 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
628
629 int dim = kappa2index[ ikappa + 1 ] - kappa2index[ ikappa ];
630 double alpha = sqrt( sectorTwoSR[ ikappa ] + 1.0 );
631 int inc = 1;
632 dscal_( &dim, &alpha, storage + kappa2index[ ikappa ], &inc );
633
634 }
635
636 }
637
638 void CheMPS2::Sobject::symm2prog(){
639
640 #pragma omp parallel for schedule(dynamic)
641 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
642
643 int dim = kappa2index[ ikappa + 1 ] - kappa2index[ ikappa ];
644 double alpha = 1.0 / sqrt( sectorTwoSR[ ikappa ] + 1.0 );
645 int inc = 1;
646 dscal_( &dim, &alpha, storage + kappa2index[ ikappa ], &inc );
647
648 }
649
650 }
651
652 void CheMPS2::Sobject::addNoise( const double NoiseLevel ){
653
654 for ( int cnt = 0; cnt < gKappa2index( gNKappa() ); cnt++ ){
655 const double RN = ( ( double ) rand() ) / RAND_MAX - 0.5;
656 gStorage()[ cnt ] += RN * NoiseLevel;
657 }
658
659 }
660
661
662