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 <algorithm>
21 #include <stdlib.h>
22 #include <assert.h>
23 #include <math.h>
24
25 #include "SyBookkeeper.h"
26 #include "Irreps.h"
27 #include "Options.h"
28
SyBookkeeper(const Problem * Prob,const int D)29 CheMPS2::SyBookkeeper::SyBookkeeper( const Problem * Prob, const int D ){
30
31 this->Prob = Prob;
32 Irreps temp( Prob->gSy() );
33 this->num_irreps = temp.getNumberOfIrreps();
34
35 // Allocate the arrays
36 allocate_arrays();
37
38 // Fill FCIdim
39 fillFCIdim();
40
41 // Copy FCIdim to CURdim
42 CopyDim( FCIdim, CURdim );
43
44 // Scale the CURdim
45 ScaleCURdim( D, 1, gL() - 1 );
46
47 assert( IsPossible() );
48
49 }
50
SyBookkeeper(const SyBookkeeper & tocopy)51 CheMPS2::SyBookkeeper::SyBookkeeper( const SyBookkeeper & tocopy ){
52
53 this->Prob = tocopy.gProb();
54 Irreps temp( Prob->gSy() );
55 this->num_irreps = temp.getNumberOfIrreps();
56
57 // Allocate the arrays
58 allocate_arrays();
59
60 // Fill FCIdim
61 fillFCIdim();
62
63 // Copy the CURdim
64 for ( int boundary = 0; boundary <= gL(); boundary++ ){
65 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
66 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
67 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
68 CURdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ] = tocopy.gCurrentDim( boundary, N, TwoS, irrep );
69 }
70 }
71 }
72 }
73
74 }
75
allocate_arrays()76 void CheMPS2::SyBookkeeper::allocate_arrays(){
77
78 // Set the min and max particle number and spin
79 Nmin = new int[ gL() + 1 ];
80 Nmax = new int[ gL() + 1 ];
81 TwoSmin = new int*[ gL() + 1 ];
82 TwoSmax = new int*[ gL() + 1 ];
83 for ( int boundary = 0; boundary <= gL(); boundary++ ){
84 Nmin[ boundary ] = std::max( std::max( 0, gN() + 2 * ( boundary - gL() ) ), boundary - gL() + ( gN() + gTwoS() ) / 2 );
85 Nmax[ boundary ] = std::min( std::min( 2 * boundary, gN() ), boundary + ( gN() - gTwoS() ) / 2 );
86 TwoSmin[ boundary ] = new int[ Nmax[ boundary ] - Nmin[ boundary ] + 1 ];
87 TwoSmax[ boundary ] = new int[ Nmax[ boundary ] - Nmin[ boundary ] + 1 ];
88 for ( int N = Nmin[ boundary ]; N <= Nmax[ boundary ]; N++ ){
89 const int temporary = gL() - boundary - abs( gN() - N - gL() + boundary );
90 TwoSmin[ boundary ][ N - Nmin[ boundary ] ] = std::max( N % 2, gTwoS() - temporary );
91 TwoSmax[ boundary ][ N - Nmin[ boundary ] ] = std::min( boundary - abs( boundary - N ), gTwoS() + temporary );
92 }
93 }
94
95 // FCIdim & CURdim memory allocation
96 FCIdim = new int***[ gL() + 1 ];
97 CURdim = new int***[ gL() + 1 ];
98 for ( int boundary = 0; boundary <= gL(); boundary++ ){
99 FCIdim[ boundary ] = new int**[ gNmax( boundary ) - gNmin( boundary ) + 1 ];
100 CURdim[ boundary ] = new int**[ gNmax( boundary ) - gNmin( boundary ) + 1 ];
101 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
102 FCIdim[ boundary ][ N - gNmin( boundary ) ] = new int*[ ( gTwoSmax( boundary, N ) - gTwoSmin( boundary, N ) ) / 2 + 1 ];
103 CURdim[ boundary ][ N - gNmin( boundary ) ] = new int*[ ( gTwoSmax( boundary, N ) - gTwoSmin( boundary, N ) ) / 2 + 1 ];
104 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
105 FCIdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ] = new int[ num_irreps ];
106 CURdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ] = new int[ num_irreps ];
107 }
108 }
109 }
110
111 }
112
~SyBookkeeper()113 CheMPS2::SyBookkeeper::~SyBookkeeper(){
114
115 for ( int boundary = 0; boundary <= gL(); boundary++ ){
116 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
117 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
118 delete [] FCIdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ];
119 delete [] CURdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ];
120 }
121 delete [] FCIdim[ boundary ][ N - gNmin( boundary ) ];
122 delete [] CURdim[ boundary ][ N - gNmin( boundary ) ];
123 }
124 delete [] FCIdim[ boundary ];
125 delete [] CURdim[ boundary ];
126 }
127 delete [] FCIdim;
128 delete [] CURdim;
129
130 for ( int boundary = 0; boundary <= gL(); boundary++ ){
131 delete [] TwoSmin[ boundary ];
132 delete [] TwoSmax[ boundary ];
133 }
134 delete [] TwoSmin;
135 delete [] TwoSmax;
136 delete [] Nmin;
137 delete [] Nmax;
138
139 }
140
gProb() const141 const CheMPS2::Problem * CheMPS2::SyBookkeeper::gProb() const{ return Prob; }
142
gL() const143 int CheMPS2::SyBookkeeper::gL() const{ return Prob->gL(); }
144
gIrrep(const int orbital) const145 int CheMPS2::SyBookkeeper::gIrrep( const int orbital ) const{ return Prob->gIrrep( orbital ); }
146
gTwoS() const147 int CheMPS2::SyBookkeeper::gTwoS() const{ return Prob->gTwoS(); }
148
gN() const149 int CheMPS2::SyBookkeeper::gN() const{ return Prob->gN(); }
150
gIrrep() const151 int CheMPS2::SyBookkeeper::gIrrep() const{ return Prob->gIrrep(); }
152
getNumberOfIrreps() const153 int CheMPS2::SyBookkeeper::getNumberOfIrreps() const{ return num_irreps; }
154
gNmin(const int boundary) const155 int CheMPS2::SyBookkeeper::gNmin( const int boundary ) const{ return Nmin[ boundary ]; }
156
gNmax(const int boundary) const157 int CheMPS2::SyBookkeeper::gNmax( const int boundary ) const{ return Nmax[ boundary ]; }
158
gTwoSmin(const int boundary,const int N) const159 int CheMPS2::SyBookkeeper::gTwoSmin( const int boundary, const int N) const{ return TwoSmin[ boundary ][ N - Nmin[ boundary ] ]; }
160
gTwoSmax(const int boundary,const int N) const161 int CheMPS2::SyBookkeeper::gTwoSmax( const int boundary, const int N) const{ return TwoSmax[ boundary ][ N - Nmin[ boundary ] ]; }
162
gFCIdim(const int boundary,const int N,const int TwoS,const int irrep) const163 int CheMPS2::SyBookkeeper::gFCIdim( const int boundary, const int N, const int TwoS, const int irrep ) const{ return gDimPrivate( FCIdim, boundary, N, TwoS, irrep ); }
164
gCurrentDim(const int boundary,const int N,const int TwoS,const int irrep) const165 int CheMPS2::SyBookkeeper::gCurrentDim( const int boundary, const int N, const int TwoS, const int irrep ) const{ return gDimPrivate( CURdim, boundary, N, TwoS, irrep ); }
166
SetDim(const int boundary,const int N,const int TwoS,const int irrep,const int value)167 void CheMPS2::SyBookkeeper::SetDim( const int boundary, const int N, const int TwoS, const int irrep, const int value ){
168
169 if ( gFCIdim( boundary, N, TwoS, irrep ) != 0 ){
170 CURdim[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ] = value;
171 }
172
173 }
174
fillFCIdim()175 void CheMPS2::SyBookkeeper::fillFCIdim(){
176
177 // On the left-hand side only the trivial symmetry sector is allowed
178 for ( int irrep = 0; irrep < num_irreps; irrep++ ){ FCIdim[ 0 ][ 0 ][ 0 ][ irrep ] = 0; }
179 FCIdim[ 0 ][ 0 ][ 0 ][ 0 ] = 1;
180
181 // Fill boundaries 1 to L from left to right
182 fill_fci_dim_right( FCIdim, 1, gL() );
183
184 // Remember the FCI virtual dimension at the RHS
185 const int rhs = FCIdim[ gL() ][ 0 ][ 0 ][ gIrrep() ];
186
187 // On the right-hand side only the targeted symmetry sector is allowed
188 for ( int irrep = 0; irrep < num_irreps; irrep++ ){ FCIdim[ gL() ][ 0 ][ 0 ][ irrep ] = 0; }
189 FCIdim[ gL() ][ 0 ][ 0 ][ gIrrep() ] = std::min( 1, rhs );
190
191 // Fill boundarties 0 to L - 1 from right to left
192 fill_fci_dim_left( FCIdim, 0, gL() - 1 );
193
194 }
195
fill_fci_dim_right(int **** storage,const int start,const int stop)196 void CheMPS2::SyBookkeeper::fill_fci_dim_right( int **** storage, const int start, const int stop ){
197
198 for ( int boundary = start; boundary <= stop; boundary++ ){
199 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
200 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
201 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
202 const int value = std::min( CheMPS2::SYBK_dimensionCutoff,
203 gDimPrivate( storage, boundary - 1, N, TwoS , irrep )
204 + gDimPrivate( storage, boundary - 1, N - 2, TwoS , irrep )
205 + gDimPrivate( storage, boundary - 1, N - 1, TwoS + 1, Irreps::directProd( irrep, gIrrep( boundary - 1 ) ) )
206 + gDimPrivate( storage, boundary - 1, N - 1, TwoS - 1, Irreps::directProd( irrep, gIrrep( boundary - 1 ) ) ) );
207 storage[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ] = value;
208 }
209 }
210 }
211 }
212
213 }
214
fill_fci_dim_left(int **** storage,const int start,const int stop)215 void CheMPS2::SyBookkeeper::fill_fci_dim_left( int **** storage, const int start, const int stop ){
216
217 for ( int boundary = stop; boundary >= start; boundary-- ){
218 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
219 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
220 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
221 const int value = std::min( gDimPrivate( storage, boundary, N, TwoS, irrep ),
222 std::min( CheMPS2::SYBK_dimensionCutoff,
223 gDimPrivate( storage, boundary + 1, N, TwoS , irrep )
224 + gDimPrivate( storage, boundary + 1, N + 2, TwoS , irrep )
225 + gDimPrivate( storage, boundary + 1, N + 1, TwoS + 1, Irreps::directProd( irrep, gIrrep( boundary ) ) )
226 + gDimPrivate( storage, boundary + 1, N + 1, TwoS - 1, Irreps::directProd( irrep, gIrrep( boundary ) ) ) ) );
227 storage[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ] = value;
228 }
229 }
230 }
231 }
232
233 }
234
CopyDim(int **** origin,int **** target)235 void CheMPS2::SyBookkeeper::CopyDim( int **** origin, int **** target ){
236
237 for ( int boundary = 0; boundary <= gL(); boundary++ ){
238 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
239 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
240 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
241 target[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ]
242 = origin[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ];
243 }
244 }
245 }
246 }
247
248 }
249
ScaleCURdim(const int virtual_dim,const int start,const int stop)250 void CheMPS2::SyBookkeeper::ScaleCURdim( const int virtual_dim, const int start, const int stop ){
251
252 for ( int boundary = start; boundary <= stop; boundary++ ){
253
254 const int totaldim = gTotDimAtBound( boundary );
255
256 if ( totaldim > virtual_dim ){
257 double factor = ( 1.0 * virtual_dim ) / totaldim;
258 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
259 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
260 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
261 const int value = ( int )( ceil( factor * gCurrentDim( boundary, N, TwoS, irrep ) ) + 0.1 );
262 SetDim( boundary, N, TwoS, irrep, value );
263 }
264 }
265 }
266 }
267 }
268
269 }
270
gDimPrivate(int **** storage,const int boundary,const int N,const int TwoS,const int irrep) const271 int CheMPS2::SyBookkeeper::gDimPrivate( int **** storage, const int boundary, const int N, const int TwoS, const int irrep ) const{
272
273 if (( boundary < 0 ) || ( boundary > gL() )){ return 0; }
274 if (( N > gNmax( boundary ) ) || ( N < gNmin( boundary ) )){ return 0; }
275 if (( TwoS % 2 ) != ( gTwoSmin( boundary, N ) % 2 )){ return 0; }
276 if (( TwoS < gTwoSmin( boundary, N ) ) || ( TwoS > gTwoSmax( boundary, N ) )){ return 0; }
277 if (( irrep < 0 ) || ( irrep >= num_irreps )){ return 0; }
278 return storage[ boundary ][ N - gNmin( boundary ) ][ ( TwoS - gTwoSmin( boundary, N ) ) / 2 ][ irrep ];
279
280 }
281
gMaxDimAtBound(const int boundary) const282 int CheMPS2::SyBookkeeper::gMaxDimAtBound( const int boundary ) const{
283
284 int max_dim = 0;
285 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
286 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
287 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
288 const int dim = gCurrentDim( boundary, N, TwoS, irrep );
289 if ( dim > max_dim ){ max_dim = dim; }
290 }
291 }
292 }
293 return max_dim;
294
295 }
296
gTotDimAtBound(const int boundary) const297 int CheMPS2::SyBookkeeper::gTotDimAtBound( const int boundary ) const{
298
299 int tot_dim = 0;
300 for ( int N = gNmin( boundary ); N <= gNmax( boundary ); N++ ){
301 for ( int TwoS = gTwoSmin( boundary, N ); TwoS <= gTwoSmax( boundary, N ); TwoS += 2 ){
302 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
303 tot_dim += gCurrentDim( boundary, N, TwoS, irrep );
304 }
305 }
306 }
307 return tot_dim;
308
309 }
310
restart(const int start,const int stop,const int virtual_dim)311 void CheMPS2::SyBookkeeper::restart( const int start, const int stop, const int virtual_dim ){
312
313 fill_fci_dim_right( CURdim, start, stop );
314 fill_fci_dim_left( CURdim, start, stop );
315 ScaleCURdim( virtual_dim, start, stop );
316
317 }
318
IsPossible() const319 bool CheMPS2::SyBookkeeper::IsPossible() const{
320
321 return ( gCurrentDim( gL(), gN(), gTwoS(), gIrrep() ) == 1 );
322
323 }
324
325
326