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 <climits>
21 #include <assert.h>
22 #include <iostream>
23 #include <math.h>
24 #include <sys/stat.h>
25 #include <sys/time.h>
26 #include <algorithm>
27
28 using std::cout;
29 using std::endl;
30
31 #include "FCI.h"
32 #include "Irreps.h"
33 #include "Lapack.h"
34 #include "Davidson.h"
35 #include "ConjugateGradient.h"
36
FCI(Hamiltonian * Ham,const unsigned int theNel_up,const unsigned int theNel_down,const int TargetIrrep_in,const double maxMemWorkMB_in,const int FCIverbose_in)37 CheMPS2::FCI::FCI(Hamiltonian * Ham, const unsigned int theNel_up, const unsigned int theNel_down, const int TargetIrrep_in, const double maxMemWorkMB_in, const int FCIverbose_in){
38
39 // Copy the basic information
40 FCIverbose = FCIverbose_in;
41 maxMemWorkMB = maxMemWorkMB_in;
42 L = Ham->getL();
43 assert( theNel_up <= L );
44 assert( theNel_down <= L );
45 assert( maxMemWorkMB > 0.0 );
46 Nel_up = theNel_up;
47 Nel_down = theNel_down;
48
49 // Construct the irrep product table and the list with the orbitals irreps
50 num_irreps = Irreps::getNumberOfIrreps( Ham->getNGroup() );
51 TargetIrrep = TargetIrrep_in;
52 orb2irrep = new int[ L ];
53 for (unsigned int orb = 0; orb < L; orb++){ orb2irrep[ orb ] = Ham->getOrbitalIrrep( orb ); }
54
55 /* Copy the Hamiltonian over:
56 G_ij = T_ij - 0.5 \sum_k <ik|kj> and ERI_{ijkl} = <ij|kl>
57 <ij|kl> is the electron repulsion integral, int dr1 dr2 i(r1) j(r1) k(r2) l(r2) / |r1-r2| */
58 Econstant = Ham->getEconst();
59 Gmat = new double[ L * L ];
60 ERI = new double[ L * L * L * L ];
61 for (unsigned int orb1 = 0; orb1 < L; orb1++){
62 for (unsigned int orb2 = 0; orb2 < L; orb2++){
63 double tempvar = 0.0;
64 for (unsigned int orb3 = 0; orb3 < L; orb3++){
65 tempvar += Ham->getVmat( orb1, orb3, orb3, orb2 );
66 for (unsigned int orb4 = 0; orb4 < L; orb4++){
67 // CheMPS2::Hamiltonian uses physics notation ; ERI chemists notation.
68 ERI[ orb1 + L * ( orb2 + L * ( orb3 + L * orb4 ) ) ] = Ham->getVmat( orb1 , orb3 , orb2 , orb4 );
69 }
70 }
71 Gmat[ orb1 + L * orb2 ] = Ham->getTmat( orb1 , orb2 ) - 0.5 * tempvar;
72 }
73 }
74
75 // Set all other internal variables
76 StartupCountersVsBitstrings();
77 StartupLookupTables();
78 StartupIrrepCenter();
79
80 }
81
~FCI()82 CheMPS2::FCI::~FCI(){
83
84 // FCI::FCI
85 delete [] orb2irrep;
86 delete [] Gmat;
87 delete [] ERI;
88
89 // FCI::StartupCountersVsBitstrings
90 for ( unsigned int irrep=0; irrep<num_irreps; irrep++ ){
91 delete [] str2cnt_up[irrep];
92 delete [] str2cnt_down[irrep];
93 delete [] cnt2str_up[irrep];
94 delete [] cnt2str_down[irrep];
95 }
96 delete [] str2cnt_up;
97 delete [] str2cnt_down;
98 delete [] cnt2str_up;
99 delete [] cnt2str_down;
100 delete [] numPerIrrep_up;
101 delete [] numPerIrrep_down;
102
103 // FCI::StartupLookupTables
104 for ( unsigned int irrep = 0; irrep < num_irreps; irrep++ ){
105 for ( unsigned int ij = 0; ij < L * L; ij++ ){
106 delete [] lookup_cnt_alpha[irrep][ij];
107 delete [] lookup_cnt_beta[irrep][ij];
108 delete [] lookup_sign_alpha[irrep][ij];
109 delete [] lookup_sign_beta[irrep][ij];
110 }
111 delete [] lookup_cnt_alpha[irrep];
112 delete [] lookup_cnt_beta[irrep];
113 delete [] lookup_sign_alpha[irrep];
114 delete [] lookup_sign_beta[irrep];
115 }
116 delete [] lookup_cnt_alpha;
117 delete [] lookup_cnt_beta;
118 delete [] lookup_sign_alpha;
119 delete [] lookup_sign_beta;
120
121 // FCI::StartupIrrepCenter
122 for ( unsigned int irrep=0; irrep<num_irreps; irrep++ ){
123 delete [] irrep_center_crea_orb[irrep];
124 delete [] irrep_center_anni_orb[irrep];
125 delete [] irrep_center_jumps[irrep];
126 }
127 delete [] irrep_center_crea_orb;
128 delete [] irrep_center_anni_orb;
129 delete [] irrep_center_jumps;
130 delete [] irrep_center_num;
131 delete [] HXVworksmall;
132 delete [] HXVworkbig1;
133 delete [] HXVworkbig2;
134
135 }
136
StartupCountersVsBitstrings()137 void CheMPS2::FCI::StartupCountersVsBitstrings(){
138
139 // Can you represent the alpha and beta Slater determinants as unsigned integers?
140 assert( L <= CHAR_BIT * sizeof(unsigned int) );
141
142 // Variable which is only needed here: 2^L
143 unsigned int TwoPowL = 1;
144 for (unsigned int orb = 0; orb < L; orb++){ TwoPowL *= 2; }
145
146 // Create the required arrays to perform the conversions between counters and bitstrings
147 numPerIrrep_up = new unsigned int[ num_irreps ];
148 numPerIrrep_down = new unsigned int[ num_irreps ];
149 str2cnt_up = new int*[ num_irreps ];
150 str2cnt_down = new int*[ num_irreps ];
151 cnt2str_up = new unsigned int*[ num_irreps ];
152 cnt2str_down = new unsigned int*[ num_irreps ];
153
154 for (unsigned int irrep = 0; irrep < num_irreps; irrep++){
155 numPerIrrep_up [ irrep ] = 0;
156 numPerIrrep_down[ irrep ] = 0;
157 str2cnt_up [ irrep ] = new int[ TwoPowL ];
158 str2cnt_down[ irrep ] = new int[ TwoPowL ];
159 }
160
161 int * bits = new int[ L ]; // Temporary helper array
162
163 // Loop over all allowed bit strings in the spinless fermion Fock space
164 for (unsigned int bitstring = 0; bitstring < TwoPowL; bitstring++){
165
166 // Find the number of particles and the irrep which correspond to each basis vector
167 str2bits( L , bitstring , bits );
168 unsigned int Nparticles = 0;
169 int irrep = 0;
170 for (unsigned int orb=0; orb<L; orb++){
171 if ( bits[orb] ){
172 Nparticles++;
173 irrep = Irreps::directProd( irrep, getOrb2Irrep( orb ) );
174 }
175 }
176
177 // If allowed: set the corresponding str2cnt to the correct counter and keep track of the number of allowed vectors
178 for ( unsigned int irr = 0; irr < num_irreps; irr++ ){
179 str2cnt_up [ irr ][ bitstring ] = -1;
180 str2cnt_down[ irr ][ bitstring ] = -1;
181 }
182 if ( Nparticles == Nel_up ){
183 str2cnt_up[ irrep ][ bitstring ] = numPerIrrep_up[ irrep ];
184 numPerIrrep_up[ irrep ]++;
185 }
186 if ( Nparticles == Nel_down ){
187 str2cnt_down[ irrep ][ bitstring ] = numPerIrrep_down[ irrep ];
188 numPerIrrep_down[ irrep ]++;
189 }
190
191 }
192
193 // Fill the reverse info array: cnt2str
194 for ( unsigned int irrep = 0; irrep < num_irreps; irrep++ ){
195
196 if ( FCIverbose>1 ){
197 cout << "FCI::Startup : For irrep " << irrep << " there are " << numPerIrrep_up [ irrep ] << " alpha Slater determinants and "
198 << numPerIrrep_down[ irrep ] << " beta Slater determinants." << endl;
199 }
200
201 cnt2str_up [ irrep ] = new unsigned int[ numPerIrrep_up [ irrep ] ];
202 cnt2str_down[ irrep ] = new unsigned int[ numPerIrrep_down[ irrep ] ];
203 for (unsigned int bitstring = 0; bitstring < TwoPowL; bitstring++){
204 if ( str2cnt_up [ irrep ][ bitstring ] != -1 ){ cnt2str_up [ irrep ][ str2cnt_up [ irrep ][ bitstring ] ] = bitstring; }
205 if ( str2cnt_down[ irrep ][ bitstring ] != -1 ){ cnt2str_down[ irrep ][ str2cnt_down[ irrep ][ bitstring ] ] = bitstring; }
206 }
207
208 }
209
210 delete [] bits; // Delete temporary helper array
211
212 }
213
StartupLookupTables()214 void CheMPS2::FCI::StartupLookupTables(){
215
216 // Create a bunch of stuff
217 lookup_cnt_alpha = new int**[ num_irreps ];
218 lookup_cnt_beta = new int**[ num_irreps ];
219 lookup_sign_alpha = new int**[ num_irreps ];
220 lookup_sign_beta = new int**[ num_irreps ];
221
222 int * bits = new int[ L ]; // Temporary helper array
223
224 // Quick lookup tables for " sign | new > = E^spinproj_{ij} | old >
225 for ( unsigned int irrep = 0; irrep < num_irreps; irrep++ ){
226
227 const unsigned int num_up = numPerIrrep_up [ irrep ];
228 const unsigned int num_down = numPerIrrep_down[ irrep ];
229
230 lookup_cnt_alpha [ irrep ] = new int*[ L * L ];
231 lookup_cnt_beta [ irrep ] = new int*[ L * L ];
232 lookup_sign_alpha[ irrep ] = new int*[ L * L ];
233 lookup_sign_beta [ irrep ] = new int*[ L * L ];
234
235 for ( unsigned int ij = 0; ij < L * L; ij++ ){
236
237 lookup_cnt_alpha [ irrep ][ ij ] = new int[ num_up ];
238 lookup_cnt_beta [ irrep ][ ij ] = new int[ num_down ];
239 lookup_sign_alpha[ irrep ][ ij ] = new int[ num_up ];
240 lookup_sign_beta [ irrep ][ ij ] = new int[ num_down ];
241
242 for ( unsigned int cnt_new_alpha = 0; cnt_new_alpha < num_up; cnt_new_alpha++ ){
243 // Check for the sign. If no check for sign, you multiply with sign 0 and everything should be OK...
244 lookup_cnt_alpha [ irrep ][ ij ][ cnt_new_alpha ] = 0;
245 lookup_sign_alpha[ irrep ][ ij ][ cnt_new_alpha ] = 0;
246 }
247 for ( unsigned int cnt_new_beta = 0; cnt_new_beta < num_down; cnt_new_beta++ ){
248 // Check for the sign. If no check for sign, you multiply with sign and everything should be OK...
249 lookup_cnt_beta [ irrep ][ ij ][ cnt_new_beta ] = 0;
250 lookup_sign_beta[ irrep ][ ij ][ cnt_new_beta ] = 0;
251 }
252 }
253
254 for ( unsigned int cnt_new_alpha = 0; cnt_new_alpha < num_up; cnt_new_alpha++ ){
255
256 str2bits( L , cnt2str_up[ irrep ][ cnt_new_alpha ] , bits );
257
258 int phase_crea = 1;
259 for ( unsigned int crea = 0; crea < L; crea++ ){
260 if ( bits[ crea ] ){
261 bits[ crea ] = 0;
262
263 int phase_anni = 1;
264 for ( unsigned int anni = 0; anni < L; anni++ ){
265 if ( !(bits[ anni ]) ){
266 bits[ anni ] = 1;
267
268 const int irrep_old = Irreps::directProd( irrep , Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) ) );
269 const int cnt_old = str2cnt_up[ irrep_old ][ bits2str( L , bits ) ];
270 const int phase = phase_crea * phase_anni;
271
272 lookup_cnt_alpha [ irrep ][ crea + L * anni ][ cnt_new_alpha ] = cnt_old;
273 lookup_sign_alpha[ irrep ][ crea + L * anni ][ cnt_new_alpha ] = phase;
274
275 bits[ anni ] = 0;
276 } else {
277 phase_anni *= -1;
278 }
279 }
280
281 bits[ crea ] = 1;
282 phase_crea *= -1;
283 }
284 }
285 }
286
287 for ( unsigned int cnt_new_beta = 0; cnt_new_beta < num_down; cnt_new_beta++ ){
288
289 str2bits( L , cnt2str_down[ irrep ][ cnt_new_beta ] , bits );
290
291 int phase_crea = 1;
292 for ( unsigned int crea = 0; crea < L; crea++ ){
293 if ( bits[ crea ] ){
294 bits[ crea ] = 0;
295
296 int phase_anni = 1;
297 for ( unsigned int anni = 0; anni < L; anni++ ){
298 if ( !(bits[ anni ]) ){
299 bits[ anni ] = 1;
300
301 const int irrep_old = Irreps::directProd( irrep , Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) ) );
302 const int cnt_old = str2cnt_down[ irrep_old ][ bits2str( L , bits ) ];
303 const int phase = phase_crea * phase_anni;
304
305 lookup_cnt_beta [ irrep ][ crea + L * anni ][ cnt_new_beta ] = cnt_old;
306 lookup_sign_beta [ irrep ][ crea + L * anni ][ cnt_new_beta ] = phase;
307
308 bits[ anni ] = 0;
309 } else {
310 phase_anni *= -1;
311 }
312 }
313
314 bits[ crea ] = 1;
315 phase_crea *= -1;
316 }
317 }
318 }
319
320 }
321
322 delete [] bits; // Delete temporary helper array
323
324 }
325
StartupIrrepCenter()326 void CheMPS2::FCI::StartupIrrepCenter(){
327
328 // Find the orbital combinations which can form a center irrep
329 irrep_center_num = new unsigned int [ num_irreps ];
330 irrep_center_crea_orb = new unsigned int*[ num_irreps ];
331 irrep_center_anni_orb = new unsigned int*[ num_irreps ];
332
333 for ( unsigned int irrep_center = 0; irrep_center < num_irreps; irrep_center++ ){
334 const int irrep_center_const_signed = irrep_center;
335
336 irrep_center_num[ irrep_center ] = 0;
337 for ( unsigned int crea = 0; crea < L; crea++ ){
338 for ( unsigned int anni = crea; anni < L; anni++ ){
339 if ( Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) ) == irrep_center_const_signed ){
340 irrep_center_num[ irrep_center ] += 1;
341 }
342 }
343 }
344 irrep_center_crea_orb[ irrep_center ] = new unsigned int[ irrep_center_num[ irrep_center ] ];
345 irrep_center_anni_orb[ irrep_center ] = new unsigned int[ irrep_center_num[ irrep_center ] ];
346 irrep_center_num[ irrep_center ] = 0;
347 for ( unsigned int creator = 0; creator < L; creator++ ){
348 for ( unsigned int annihilator = creator; annihilator < L; annihilator++){
349 if ( Irreps::directProd( getOrb2Irrep( creator ) , getOrb2Irrep( annihilator ) ) == irrep_center_const_signed ){
350 irrep_center_crea_orb[ irrep_center ][ irrep_center_num[ irrep_center ] ] = creator;
351 irrep_center_anni_orb[ irrep_center ][ irrep_center_num[ irrep_center ] ] = annihilator;
352 irrep_center_num[ irrep_center ] += 1;
353 }
354 }
355 }
356
357 }
358
359 irrep_center_jumps = new unsigned int*[ num_irreps ];
360 HXVsizeWorkspace = 0;
361 for ( unsigned int irrep_center = 0; irrep_center < num_irreps; irrep_center++ ){
362 unsigned long long check = 0;
363 irrep_center_jumps[ irrep_center ] = new unsigned int[ num_irreps + 1 ];
364 const int localTargetIrrep = Irreps::directProd( irrep_center, getTargetIrrep() );
365 irrep_center_jumps[ irrep_center ][ 0 ] = 0;
366 for ( unsigned int irrep_up = 0; irrep_up < num_irreps; irrep_up++ ){
367 const int irrep_down = Irreps::directProd( irrep_up, localTargetIrrep );
368 check += ((unsigned long long) numPerIrrep_up[ irrep_up ] ) * ((unsigned long long) numPerIrrep_down[ irrep_down ] );
369 unsigned int temp = numPerIrrep_up[ irrep_up ] * numPerIrrep_down[ irrep_down ];
370 irrep_center_jumps[ irrep_center ][ irrep_up + 1 ] = irrep_center_jumps[ irrep_center ][ irrep_up ] + temp;
371 HXVsizeWorkspace = std::max( HXVsizeWorkspace, ((unsigned long long) irrep_center_num[ irrep_center ] ) * ((unsigned long long) temp ) );
372 }
373 assert( check <= ((unsigned int) INT_MAX ) ); // Length of FCI vectors should be less then the SIGNED integer size (to be able to call lapack)
374 }
375 if ( FCIverbose > 0 ){
376 cout << "FCI::Startup : Number of variables in the FCI vector = " << getVecLength( 0 ) << endl;
377 double num_megabytes = ( 2.0 * sizeof(double) * HXVsizeWorkspace ) / 1048576;
378 cout << "FCI::Startup : Without additional loops the FCI matrix-vector product requires a workspace of " << num_megabytes << " MB memory." << endl;
379 if ( maxMemWorkMB < num_megabytes ){
380 HXVsizeWorkspace = (unsigned int) ceil( ( maxMemWorkMB * 1048576 ) / ( 2 * sizeof(double) ) );
381 num_megabytes = ( 2.0 * sizeof(double) * HXVsizeWorkspace ) / 1048576;
382 cout << " For practical purposes, the workspace is constrained to " << num_megabytes << " MB memory." << endl;
383 }
384 }
385 HXVworksmall = new double[ L * L * L * L ];
386 HXVworkbig1 = new double[ HXVsizeWorkspace ];
387 HXVworkbig2 = new double[ HXVsizeWorkspace ];
388
389 }
390
str2bits(const unsigned int Lval,const unsigned int bitstring,int * bits)391 void CheMPS2::FCI::str2bits(const unsigned int Lval, const unsigned int bitstring, int * bits){
392
393 for (unsigned int bit = 0; bit < Lval; bit++){ bits[ bit ] = ( bitstring & ( 1 << bit ) ) >> bit; }
394
395 }
396
bits2str(const unsigned int Lval,int * bits)397 unsigned int CheMPS2::FCI::bits2str(const unsigned int Lval, int * bits){
398
399 unsigned int factor = 1;
400 unsigned int result = 0;
401 for (unsigned int bit = 0; bit < Lval; bit++){
402 result += bits[ bit ] * factor;
403 factor *= 2;
404 }
405 return result;
406
407 }
408
getUpIrrepOfCounter(const int irrep_center,const unsigned int counter) const409 int CheMPS2::FCI::getUpIrrepOfCounter(const int irrep_center, const unsigned int counter) const{
410
411 int irrep_up = num_irreps;
412 while ( counter < irrep_center_jumps[ irrep_center ][ irrep_up-1 ] ){ irrep_up--; }
413 return irrep_up-1;
414
415 }
416
getBitsOfCounter(const int irrep_center,const unsigned int counter,int * bits_up,int * bits_down) const417 void CheMPS2::FCI::getBitsOfCounter(const int irrep_center, const unsigned int counter, int * bits_up, int * bits_down) const{
418
419 const int localTargetIrrep = Irreps::directProd( irrep_center, TargetIrrep );
420
421 const int irrep_up = getUpIrrepOfCounter( irrep_center, counter );
422 const int irrep_down = Irreps::directProd( irrep_up, localTargetIrrep );
423
424 const unsigned int count_up = ( counter - irrep_center_jumps[ irrep_center ][ irrep_up ] ) % numPerIrrep_up[ irrep_up ];
425 const unsigned int count_down = ( counter - irrep_center_jumps[ irrep_center ][ irrep_up ] ) / numPerIrrep_up[ irrep_up ];
426
427 const unsigned int string_up = cnt2str_up [ irrep_up ][ count_up ];
428 const unsigned int string_down = cnt2str_down[ irrep_down ][ count_down ];
429
430 str2bits( L , string_up , bits_up );
431 str2bits( L , string_down , bits_down );
432
433 }
434
getFCIcoeff(int * bits_up,int * bits_down,double * vector) const435 double CheMPS2::FCI::getFCIcoeff(int * bits_up, int * bits_down, double * vector) const{
436
437 const unsigned string_up = bits2str(L, bits_up );
438 const unsigned string_down = bits2str(L, bits_down);
439
440 int irrep_up = 0;
441 int irrep_down = 0;
442 for ( unsigned int orb = 0; orb < L; orb++ ){
443 if ( bits_up [ orb ] ){ irrep_up = Irreps::directProd( irrep_up , getOrb2Irrep( orb ) ); }
444 if ( bits_down[ orb ] ){ irrep_down = Irreps::directProd( irrep_down , getOrb2Irrep( orb ) ); }
445 }
446
447 const int counter_up = str2cnt_up [ irrep_up ][ string_up ];
448 const int counter_down = str2cnt_down[ irrep_down ][ string_down ];
449
450 if (( counter_up == -1 ) || ( counter_down == -1 )){ return 0.0; }
451
452 return vector[ irrep_center_jumps[ 0 ][ irrep_up ] + counter_up + numPerIrrep_up[ irrep_up ] * counter_down ];
453
454 }
455
456 /*void CheMPS2::FCI::CheckHamDEBUG() const{
457
458 const unsigned int vecLength = getVecLength( 0 );
459
460 // Building Ham by matvec
461 double * HamHXV = new double[ vecLength * vecLength ];
462 double * workspace = new double[ vecLength ];
463 for (unsigned int count = 0; count < vecLength; count++){
464
465 ClearVector( vecLength , workspace );
466 workspace[ count ] = 1.0;
467 matvec( workspace , HamHXV + count*vecLength );
468
469 }
470
471 // Building Diag by HamDiag
472 DiagHam( workspace );
473 double RMSdiagdifference = 0.0;
474 for (unsigned int row = 0; row < vecLength; row++){
475 double diff = workspace[ row ] - HamHXV[ row + vecLength * row ];
476 RMSdiagdifference += diff * diff;
477 }
478 RMSdiagdifference = sqrt( RMSdiagdifference );
479 cout << "The RMS difference of DiagHam() and diag(HamHXV) = " << RMSdiagdifference << endl;
480
481 // Building Ham by getMatrixElement
482 int * work = new int[ 8 ];
483 int * ket_up = new int[ L ];
484 int * ket_down = new int[ L ];
485 int * bra_up = new int[ L ];
486 int * bra_down = new int[ L ];
487 double RMSconstructiondiff = 0.0;
488 for (unsigned int row = 0; row < vecLength; row++){
489 for (unsigned int col = 0; col < vecLength; col++){
490 getBitsOfCounter( 0 , row , bra_up , bra_down );
491 getBitsOfCounter( 0 , col , ket_up , ket_down );
492 double tempvar = HamHXV[ row + vecLength * col ] - GetMatrixElement( bra_up , bra_down , ket_up , ket_down , work );
493 RMSconstructiondiff += tempvar * tempvar;
494 }
495 }
496 cout << "The RMS difference of HamHXV - HamMXELEM = " << RMSconstructiondiff << endl;
497 delete [] work;
498 delete [] ket_up;
499 delete [] ket_down;
500 delete [] bra_up;
501 delete [] bra_down;
502
503 // Building Ham^2 by matvec
504 double * workspace2 = new double[ vecLength ];
505 for (unsigned int count = 0; count < vecLength; count++){
506
507 ClearVector( vecLength , workspace );
508 workspace[ count ] = 1.0;
509 matvec( workspace , workspace2 );
510 matvec( workspace2 , HamHXV + count*vecLength );
511
512 }
513
514 // Building diag( Ham^2 ) by DiagHamSquared
515 DiagHamSquared( workspace );
516 double RMSdiagdifference2 = 0.0;
517 for (unsigned int row = 0; row < vecLength; row++){
518 double diff = workspace[ row ] - HamHXV[ row + vecLength * row ];
519 RMSdiagdifference2 += diff * diff;
520 }
521 RMSdiagdifference2 = sqrt( RMSdiagdifference2 );
522 cout << "The RMS difference of DiagHamSquared() and diag(HamSquared by HXV) = " << RMSdiagdifference2 << endl;
523
524 delete [] workspace2;
525 delete [] workspace;
526 delete [] HamHXV;
527
528 }*/
529
excite_alpha_omp(const unsigned int dim_new_up,const unsigned int dim_old_up,const unsigned int dim_down,double * origin,double * result,int * signmap,int * countmap)530 void CheMPS2::FCI::excite_alpha_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int dim_down, double * origin, double * result, int * signmap, int * countmap ){
531
532 #pragma omp parallel for schedule(static)
533 for ( unsigned int cnt_new_up = 0; cnt_new_up < dim_new_up; cnt_new_up++ ){
534 const int sign_up = signmap[ cnt_new_up ];
535 if ( sign_up != 0 ){
536 const int cnt_old_up = countmap[ cnt_new_up ];
537 for ( unsigned int cnt_down = 0; cnt_down < dim_down; cnt_down++ ){
538 result[ cnt_new_up + dim_new_up * cnt_down ] += sign_up * origin[ cnt_old_up + dim_old_up * cnt_down ];
539 }
540 }
541 }
542
543 }
544
excite_beta_omp(const unsigned int dim_up,const unsigned int dim_new_down,double * origin,double * result,int * signmap,int * countmap)545 void CheMPS2::FCI::excite_beta_omp( const unsigned int dim_up, const unsigned int dim_new_down, double * origin, double * result, int * signmap, int * countmap ){
546
547 #pragma omp parallel for schedule(static)
548 for ( unsigned int cnt_new_down = 0; cnt_new_down < dim_new_down; cnt_new_down++ ){
549 const int sign_down = signmap[ cnt_new_down ];
550 if ( sign_down != 0 ){
551 const int cnt_old_down = countmap[ cnt_new_down ];
552 for ( unsigned int cnt_up = 0; cnt_up < dim_up; cnt_up++ ){
553 result[ cnt_up + dim_up * cnt_new_down ] += sign_down * origin[ cnt_up + dim_up * cnt_old_down ];
554 }
555 }
556 }
557
558 }
559
excite_alpha_first(const unsigned int dim_new_up,const unsigned int dim_old_up,const unsigned int start_down,const unsigned int stop_down,double * origin,double * result,int * signmap,int * countmap)560 void CheMPS2::FCI::excite_alpha_first( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ){
561
562 for ( unsigned int cnt_new_up = 0; cnt_new_up < dim_new_up; cnt_new_up++ ){
563 const int sign_up = signmap[ cnt_new_up ];
564 if ( sign_up != 0 ){
565 const int cnt_old_up = countmap[ cnt_new_up ];
566 for ( unsigned int cnt_down = start_down; cnt_down < stop_down; cnt_down++ ){
567 result[ cnt_new_up + dim_new_up * ( cnt_down - start_down ) ] += sign_up * origin[ cnt_old_up + dim_old_up * cnt_down ];
568 }
569 }
570 }
571
572 }
573
excite_beta_first(const unsigned int dim_up,const unsigned int start_down,const unsigned int stop_down,double * origin,double * result,int * signmap,int * countmap)574 void CheMPS2::FCI::excite_beta_first( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ){
575
576 for ( unsigned int cnt_new_down = start_down; cnt_new_down < stop_down; cnt_new_down++ ){
577 const int sign_down = signmap[ cnt_new_down ];
578 if ( sign_down != 0 ){
579 const int cnt_old_down = countmap[ cnt_new_down ];
580 for ( unsigned int cnt_up = 0; cnt_up < dim_up; cnt_up++ ){
581 result[ cnt_up + dim_up * ( cnt_new_down - start_down ) ] += sign_down * origin[ cnt_up + dim_up * cnt_old_down ];
582 }
583 }
584 }
585
586 }
587
excite_alpha_second_omp(const unsigned int dim_new_up,const unsigned int dim_old_up,const unsigned int start_down,const unsigned int stop_down,double * origin,double * result,int * signmap,int * countmap)588 void CheMPS2::FCI::excite_alpha_second_omp( const unsigned int dim_new_up, const unsigned int dim_old_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ){
589
590 #pragma omp parallel for schedule(static)
591 for ( unsigned int cnt_old_up = 0; cnt_old_up < dim_old_up; cnt_old_up++ ){
592 const int sign_up = signmap[ cnt_old_up ];
593 if ( sign_up != 0 ){ // Required for thread safety
594 const int cnt_new_up = countmap[ cnt_old_up ];
595 for ( unsigned int cnt_down = start_down; cnt_down < stop_down; cnt_down++ ){
596 result[ cnt_new_up + dim_new_up * cnt_down ] += sign_up * origin[ cnt_old_up + dim_old_up * ( cnt_down - start_down ) ];
597 }
598 }
599 }
600
601 }
602
excite_beta_second_omp(const unsigned int dim_up,const unsigned int start_down,const unsigned int stop_down,double * origin,double * result,int * signmap,int * countmap)603 void CheMPS2::FCI::excite_beta_second_omp( const unsigned int dim_up, const unsigned int start_down, const unsigned int stop_down, double * origin, double * result, int * signmap, int * countmap ){
604
605 #pragma omp parallel for schedule(static)
606 for ( unsigned int cnt_old_down = start_down; cnt_old_down < stop_down; cnt_old_down++ ){
607 const int sign_down = signmap[ cnt_old_down ];
608 if ( sign_down != 0 ){ // Required for thread safety
609 const int cnt_new_down = countmap[ cnt_old_down ];
610 for ( unsigned int cnt_up = 0; cnt_up < dim_up; cnt_up++ ){
611 result[ cnt_up + dim_up * cnt_new_down ] += sign_down * origin[ cnt_up + dim_up * ( cnt_old_down - start_down ) ];
612 }
613 }
614 }
615
616 }
617
matvec(double * input,double * output) const618 void CheMPS2::FCI::matvec( double * input, double * output ) const{
619
620 struct timeval start, end;
621 gettimeofday( &start, NULL );
622
623 ClearVector( getVecLength( 0 ), output );
624
625 // P.J. Knowles and N.C. Handy, A new determinant-based full configuration interaction method, Chemical Physics Letters 111 (4-5), 315-321 (1984)
626
627 // irrep_center is the center irrep of the ERI : (ij|kl) --> irrep_center = I_i x I_j = I_k x I_l
628 for ( unsigned int irrep_center = 0; irrep_center < num_irreps; irrep_center++ ){
629
630 const int irrep_target_center = Irreps::directProd( TargetIrrep, irrep_center );
631 const unsigned int num_pairs = irrep_center_num[ irrep_center ];
632 const unsigned int * center_crea_orb = irrep_center_crea_orb[ irrep_center ];
633 const unsigned int * center_anni_orb = irrep_center_anni_orb[ irrep_center ];
634 const unsigned int * zero_jumps = irrep_center_jumps[ 0 ];
635
636 for ( unsigned int irrep_center_up = 0; irrep_center_up < num_irreps; irrep_center_up++ ){
637 const int irrep_center_down = Irreps::directProd( irrep_target_center, irrep_center_up );
638 const unsigned int dim_center_up = numPerIrrep_up [ irrep_center_up ];
639 const unsigned int dim_center_down = numPerIrrep_down[ irrep_center_down ];
640 if ( dim_center_up * dim_center_down > 0 ){
641 const unsigned int blocksize_beta = HXVsizeWorkspace / std::max( (unsigned int) 1, dim_center_up * num_pairs );
642 assert( blocksize_beta > 0 ); // At least one full column should fit in the workspaces...
643 unsigned int num_block_beta = dim_center_down / blocksize_beta;
644 while ( blocksize_beta * num_block_beta < dim_center_down ){ num_block_beta++; }
645 for ( unsigned int block = 0; block < num_block_beta; block++ ){
646 const unsigned int start_center_down = block * blocksize_beta;
647 const unsigned int stop_center_down = std::min( ( block + 1 ) * blocksize_beta, dim_center_down );
648 const unsigned int size_center = dim_center_up * ( stop_center_down - start_center_down );
649 if ( size_center > 0 ){
650
651 // First build workbig1[ veccounter + size_center * pair ] = E_{i<=j} + ( 1 - delta_i==j ) E_{j>i} (irrep_center) | input > */
652 #pragma omp parallel for schedule(static)
653 for ( unsigned int pair = 0; pair < num_pairs; pair++ ){
654 double * target_space = HXVworkbig1 + size_center * pair;
655 const unsigned int crea = center_crea_orb[ pair ];
656 const unsigned int anni = center_anni_orb[ pair ];
657 const int irrep_excited = Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) );
658 const int irrep_zero_up = Irreps::directProd( irrep_excited, irrep_center_up );
659 const unsigned int dim_zero_up = numPerIrrep_up[ irrep_zero_up ];
660 for ( unsigned int count = 0; count < size_center; count++ ){ target_space[ count ] = 0.0; }
661
662 excite_alpha_first( dim_center_up, dim_zero_up, start_center_down, stop_center_down,
663 input + zero_jumps[ irrep_zero_up ],
664 target_space,
665 lookup_sign_alpha[ irrep_center_up ][ crea + L * anni ],
666 lookup_cnt_alpha [ irrep_center_up ][ crea + L * anni ] );
667
668 excite_beta_first( dim_center_up, start_center_down, stop_center_down,
669 input + zero_jumps[ irrep_center_up ],
670 target_space,
671 lookup_sign_beta[ irrep_center_down ][ crea + L * anni ],
672 lookup_cnt_beta [ irrep_center_down ][ crea + L * anni ] );
673
674 if ( anni > crea ){
675
676 excite_alpha_first( dim_center_up, dim_zero_up, start_center_down, stop_center_down,
677 input + zero_jumps[ irrep_zero_up ],
678 target_space,
679 lookup_sign_alpha[ irrep_center_up ][ anni + L * crea ],
680 lookup_cnt_alpha [ irrep_center_up ][ anni + L * crea ] );
681
682 excite_beta_first( dim_center_up, start_center_down, stop_center_down,
683 input + zero_jumps[ irrep_center_up ],
684 target_space,
685 lookup_sign_beta[ irrep_center_down ][ anni + L * crea ],
686 lookup_cnt_beta [ irrep_center_down ][ anni + L * crea ] );
687
688 }
689 }
690
691 // If irrep_center == 0, do the one-body terms
692 if ( irrep_center == 0 ){
693 for ( unsigned int pair = 0; pair < num_pairs; pair++ ){
694 HXVworksmall[ pair ] = getGmat( center_crea_orb[ pair ], center_anni_orb[ pair ] );
695 }
696 char notrans = 'N';
697 double one = 1.0;
698 int mdim = size_center;
699 int kdim = num_pairs;
700 int ndim = 1;
701 double * target = output + zero_jumps[ irrep_center_up ] + dim_center_up * start_center_down;
702 dgemm_( ¬rans, ¬rans, &mdim, &ndim, &kdim, &one, HXVworkbig1, &mdim, HXVworksmall, &kdim, &one, target, &mdim );
703 }
704
705 // Now build workbig2[ veccounter + size_center * new_pair ] = 0.5 * ( new_pair | old_pair ) * workbig1[ veccounter + size_center * old_pair ]
706 {
707 for ( unsigned int pair1 = 0; pair1 < num_pairs; pair1++ ){
708 for ( unsigned int pair2 = 0; pair2 < num_pairs; pair2++ ){
709 HXVworksmall[ pair1 + num_pairs * pair2 ]
710 = 0.5 * getERI( center_crea_orb[ pair1 ], center_anni_orb[ pair1 ] ,
711 center_crea_orb[ pair2 ], center_anni_orb[ pair2 ] );
712 }
713 }
714 char notrans = 'N';
715 double one = 1.0;
716 double set = 0.0;
717 int mdim = size_center;
718 int kdim = num_pairs;
719 int ndim = num_pairs;
720 dgemm_( ¬rans, ¬rans, &mdim, &ndim, &kdim, &one, HXVworkbig1, &mdim, HXVworksmall, &kdim, &set, HXVworkbig2, &mdim );
721 }
722
723 // Finally do output <-- E_{i<=j} + (1 - delta_{i==j}) E_{j>i} workbig2[ veccounter + size_center * pair ]
724 for ( unsigned int pair = 0; pair < num_pairs; pair++ ){
725 double * origin_space = HXVworkbig2 + size_center * pair;
726 const unsigned int crea = center_crea_orb[ pair ];
727 const unsigned int anni = center_anni_orb[ pair ];
728 const int irrep_excited = Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) );
729 const int irrep_zero_up = Irreps::directProd( irrep_excited, irrep_center_up );
730 const unsigned int dim_zero_up = numPerIrrep_up[ irrep_zero_up ];
731
732 excite_alpha_second_omp( dim_zero_up, dim_center_up, start_center_down, stop_center_down,
733 origin_space,
734 output + zero_jumps[ irrep_zero_up ],
735 lookup_sign_alpha[ irrep_center_up ][ anni + L * crea ],
736 lookup_cnt_alpha [ irrep_center_up ][ anni + L * crea ] );
737
738 excite_beta_second_omp( dim_center_up, start_center_down, stop_center_down,
739 origin_space,
740 output + zero_jumps[ irrep_center_up ],
741 lookup_sign_beta[ irrep_center_down ][ anni + L * crea ],
742 lookup_cnt_beta [ irrep_center_down ][ anni + L * crea ] );
743
744 if ( anni > crea ){
745
746 excite_alpha_second_omp( dim_zero_up, dim_center_up, start_center_down, stop_center_down,
747 origin_space,
748 output + zero_jumps[ irrep_zero_up ],
749 lookup_sign_alpha[ irrep_center_up ][ crea + L * anni ],
750 lookup_cnt_alpha [ irrep_center_up ][ crea + L * anni ] );
751
752 excite_beta_second_omp( dim_center_up, start_center_down, stop_center_down,
753 origin_space,
754 output + zero_jumps[ irrep_center_up ],
755 lookup_sign_beta[ irrep_center_down ][ crea + L * anni ],
756 lookup_cnt_beta [ irrep_center_down ][ crea + L * anni ] );
757
758 }
759 }
760 }
761 }
762 }
763 }
764 }
765
766 gettimeofday( &end, NULL );
767 const double elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
768 if ( FCIverbose >= 1 ){ cout << "FCI::matvec : Wall time = " << elapsed << " seconds" << endl; }
769
770 }
771
apply_excitation(double * orig_vector,double * result_vector,const int crea,const int anni,const int orig_target_irrep) const772 void CheMPS2::FCI::apply_excitation( double * orig_vector, double * result_vector, const int crea, const int anni, const int orig_target_irrep ) const{
773
774 const int excitation_irrep = Irreps::directProd( getOrb2Irrep( crea ), getOrb2Irrep( anni ) );
775 const int result_target_irrep = Irreps::directProd( excitation_irrep, orig_target_irrep );
776 const int orig_irrep_center = Irreps::directProd( TargetIrrep, orig_target_irrep );
777 const int result_irrep_center = Irreps::directProd( TargetIrrep, result_target_irrep );
778
779 ClearVector( getVecLength( result_irrep_center ) , result_vector );
780
781 for ( unsigned int result_irrep_up = 0; result_irrep_up < num_irreps; result_irrep_up++ ){
782
783 const int result_irrep_down = Irreps::directProd( result_irrep_up, result_target_irrep );
784 const int orig_irrep_up = Irreps::directProd( excitation_irrep, result_irrep_up );
785
786 excite_alpha_omp( numPerIrrep_up [ result_irrep_up ], // dim_new_up
787 numPerIrrep_up [ orig_irrep_up ], // dim_old_up
788 numPerIrrep_down[ result_irrep_down ], // dim_down
789 orig_vector + irrep_center_jumps[ orig_irrep_center ][ orig_irrep_up ], // origin
790 result_vector + irrep_center_jumps[ result_irrep_center ][ result_irrep_up ], // result
791 lookup_sign_alpha[ result_irrep_up ][ crea + L * anni ], // signmap
792 lookup_cnt_alpha [ result_irrep_up ][ crea + L * anni ] ); // countmap
793
794 excite_beta_omp( numPerIrrep_up [ result_irrep_up ], // dim_up
795 numPerIrrep_down[ result_irrep_down ], // dim_new_down
796 orig_vector + irrep_center_jumps[ orig_irrep_center ][ result_irrep_up ], // origin
797 result_vector + irrep_center_jumps[ result_irrep_center ][ result_irrep_up ], // result
798 lookup_sign_beta[ result_irrep_down ][ crea + L * anni ], // signmap
799 lookup_cnt_beta [ result_irrep_down ][ crea + L * anni ] ); // countmap
800
801 }
802
803 }
804
Fill2RDM(double * vector,double * two_rdm) const805 double CheMPS2::FCI::Fill2RDM(double * vector, double * two_rdm) const{
806
807 assert( Nel_up + Nel_down >= 2 );
808
809 struct timeval start, end;
810 gettimeofday(&start, NULL);
811
812 ClearVector( L*L*L*L, two_rdm );
813 const unsigned int orig_length = getVecLength( 0 );
814 unsigned int max_length = 0;
815 for ( unsigned int irrep = 0; irrep < num_irreps; irrep++ ){
816 if ( getVecLength( irrep ) > max_length ){ max_length = getVecLength( irrep ); }
817 }
818 double * workspace1 = new double[ max_length ];
819 double * workspace2 = new double[ orig_length ];
820
821 // Gamma_{ijkl} = < E_ik E_jl > - delta_jk < E_il >
822 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){ // anni1 = l
823 for ( unsigned int crea1 = anni1; crea1 < L; crea1++ ){ // crea1 = j >= l
824
825 const int irrep_center1 = Irreps::directProd( getOrb2Irrep( crea1 ), getOrb2Irrep( anni1 ) );
826 const int target_irrep1 = Irreps::directProd( TargetIrrep, irrep_center1 );
827 apply_excitation( vector, workspace1, crea1, anni1, TargetIrrep );
828
829 if ( irrep_center1 == 0 ){
830 const double value = FCIddot( orig_length, workspace1, vector ); // < E_{crea1,anni1} >
831 for ( unsigned int jk = anni1; jk < L; jk++ ){
832 two_rdm[ crea1 + L * ( jk + L * ( jk + L * anni1 ) ) ] -= value;
833 }
834 }
835
836 for ( unsigned int crea2 = anni1; crea2 < L; crea2++ ){ // crea2 = i >= l
837 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){ // anni2 = k >= l
838
839 const int irrep_center2 = Irreps::directProd( getOrb2Irrep( crea2 ), getOrb2Irrep( anni2 ) );
840 if ( irrep_center2 == irrep_center1 ){
841
842 apply_excitation( workspace1, workspace2, crea2, anni2, target_irrep1 );
843 const double value = FCIddot( orig_length, workspace2, vector ); // < E_{crea2,anni2} E_{crea1,anni1} >
844 two_rdm[ crea2 + L * ( crea1 + L * ( anni2 + L * anni1 ) ) ] += value;
845
846 }
847 }
848 }
849 }
850 }
851 delete [] workspace1;
852 delete [] workspace2;
853
854 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){
855 for ( unsigned int crea1 = anni1; crea1 < L; crea1++ ){
856 const int irrep_center1 = Irreps::directProd( getOrb2Irrep( crea1 ) , getOrb2Irrep( anni1 ) );
857 for ( unsigned int crea2 = anni1; crea2 < L; crea2++ ){
858 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){
859 const int irrep_center2 = Irreps::directProd( getOrb2Irrep( crea2 ) , getOrb2Irrep( anni2 ) );
860 if ( irrep_center2 == irrep_center1 ){
861 const double value = two_rdm[ crea2 + L * ( crea1 + L * ( anni2 + L * anni1 ) ) ];
862 two_rdm[ crea1 + L * ( crea2 + L * ( anni1 + L * anni2 ) ) ] = value;
863 two_rdm[ anni2 + L * ( anni1 + L * ( crea2 + L * crea1 ) ) ] = value;
864 two_rdm[ anni1 + L * ( anni2 + L * ( crea1 + L * crea2 ) ) ] = value;
865 }
866 }
867 }
868 }
869 }
870
871 // Calculate the FCI energy
872 double FCIenergy = getEconst();
873 for ( unsigned int orb1 = 0; orb1 < L; orb1++ ){
874 for ( unsigned int orb2 = 0; orb2 < L; orb2++ ){
875 double tempvar = 0.0;
876 double tempvar2 = 0.0;
877 for ( unsigned int orb3 = 0; orb3 < L; orb3++ ){
878 tempvar += getERI( orb1 , orb3 , orb3 , orb2 );
879 tempvar2 += two_rdm[ orb1 + L * ( orb3 + L * ( orb2 + L * orb3 ) ) ];
880 for ( unsigned int orb4 = 0; orb4 < L; orb4++ ){
881 FCIenergy += 0.5 * two_rdm[ orb1 + L * ( orb2 + L * ( orb3 + L * orb4 ) ) ] * getERI( orb1 , orb3 , orb2 , orb4 );
882 }
883 }
884 FCIenergy += ( getGmat( orb1 , orb2 ) + 0.5 * tempvar ) * tempvar2 / ( Nel_up + Nel_down - 1.0);
885 }
886 }
887
888 gettimeofday(&end, NULL);
889 const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
890 if ( FCIverbose > 0 ){ cout << "FCI::Fill2RDM : Wall time = " << elapsed << " seconds" << endl; }
891 if ( FCIverbose > 0 ){ cout << "FCI::Fill2RDM : Energy (Ham * 2-RDM) = " << FCIenergy << endl; }
892 return FCIenergy;
893
894 }
895
Fill4RDM(double * vector,double * four_rdm) const896 void CheMPS2::FCI::Fill4RDM(double * vector, double * four_rdm) const{
897
898 assert( Nel_up + Nel_down >= 4 );
899
900 struct timeval start, end;
901 gettimeofday(&start, NULL);
902
903 /*
904 Gamma_{ijkl,pqrt} = < E_ip E_jq E_kr E_lt >
905 - delta_lr < E_ip E_jq E_kt >
906 - delta_lq < E_ip E_kr E_jt >
907 - delta_lp < E_jq E_kr E_it >
908 - delta_kq < E_ip E_jr E_lt >
909 - delta_kp < E_ir E_jq E_lt >
910 - delta_jp < E_iq E_kr E_lt >
911 + delta_lr delta_kq < E_ip E_jt >
912 + delta_lr delta_kp < E_jq E_it >
913 + delta_lr delta_jp < E_iq E_kt >
914 + delta_lq delta_jr < E_ip E_kt >
915 + delta_lq delta_kp < E_ir E_jt >
916 + delta_lq delta_jp < E_kr E_it >
917 + delta_lp delta_kq < E_jr E_it >
918 + delta_lp delta_iq < E_kr E_jt >
919 + delta_lp delta_ir < E_jq E_kt >
920 + delta_kq delta_jp < E_ir E_lt >
921 + delta_kp delta_jr < E_iq E_lt >
922 - delta_jp delta_kq delta_lr < E_it >
923 - delta_kp delta_lq delta_jr < E_it >
924 - delta_kp delta_iq delta_lr < E_jt >
925 - delta_lp delta_kq delta_ir < E_jt >
926 - delta_jp delta_lq delta_ir < E_kt >
927 - delta_lp delta_iq delta_jr < E_kt >
928 */
929
930 ClearVector( L*L*L*L*L*L*L*L, four_rdm );
931 const unsigned int orig_length = getVecLength( 0 );
932 unsigned int max_length = getVecLength( 0 );
933 for ( unsigned int irrep = 1; irrep < num_irreps; irrep++ ){
934 if ( getVecLength( irrep ) > max_length ){ max_length = getVecLength( irrep ); }
935 }
936 double * workspace1 = new double[ max_length ];
937 double * workspace2 = new double[ max_length ];
938 double * workspace3 = new double[ max_length ];
939 double * workspace4 = new double[ orig_length ];
940
941 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){ // anni1 = t
942 for ( unsigned int crea1 = anni1; crea1 < L; crea1++ ){ // crea1 = l >= t
943
944 const int irrep_center1 = Irreps::directProd( getOrb2Irrep( crea1 ), getOrb2Irrep( anni1 ) );
945 const int target_irrep1 = Irreps::directProd( TargetIrrep, irrep_center1 );
946 apply_excitation( vector, workspace1, crea1, anni1, TargetIrrep );
947
948 if ( irrep_center1 == 0 ){
949
950 // value = < E_{crea1,anni1} >
951 const double value = FCIddot( orig_length, workspace1, vector );
952 for ( unsigned int p = anni1; p < L; p++ ){
953 for ( unsigned int q = anni1; q < L; q++ ){
954 for ( unsigned int r = anni1; r < L; r++ ){
955 // i j k l p q r t
956 four_rdm[ crea1 + L*( p + L*( q + L*( r + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_jp delta_kq delta_lr < E_it >
957 four_rdm[ crea1 + L*( r + L*( p + L*( q + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_kp delta_lq delta_jr < E_it >
958 four_rdm[ q + L*( crea1 + L*( p + L*( r + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_kp delta_iq delta_lr < E_jt >
959 four_rdm[ r + L*( crea1 + L*( q + L*( p + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_lp delta_kq delta_ir < E_jt >
960 four_rdm[ r + L*( p + L*( crea1 + L*( q + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_jp delta_lq delta_ir < E_kt >
961 four_rdm[ q + L*( r + L*( crea1 + L*( p + L*( p + L*( q + L*( r + L * anni1 ))))))] -= value; // - delta_lp delta_iq delta_jr < E_kt >
962 }
963 }
964 }
965
966 }
967
968 for ( unsigned int crea2 = anni1; crea2 < L; crea2++ ){ // crea2 = k >= t
969 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){ // anni2 = r >= t
970
971 const int irrep_center2 = Irreps::directProd( getOrb2Irrep( crea2 ), getOrb2Irrep( anni2 ) );
972 const int target_irrep2 = Irreps::directProd( target_irrep1, irrep_center2 );
973 apply_excitation( workspace1, workspace2, crea2, anni2, target_irrep1 );
974
975 if ( irrep_center1 == irrep_center2 ){
976
977 // value = < E_{crea2,anni2} E_{crea1,anni1} >
978 const double value = FCIddot( orig_length, workspace2, vector );
979 for ( unsigned int orb = anni1; orb < L; orb++ ){
980 for ( unsigned int xyz = anni1; xyz < L; xyz++ ){
981 // i j k l p q r t
982 four_rdm[ crea2 + L*( crea1 + L*( xyz + L*( orb + L*( anni2 + L*( xyz + L*( orb + L * anni1 ))))))] += value; // + delta_lr delta_kq < E_ip E_jt >
983 four_rdm[ crea1 + L*( crea2 + L*( xyz + L*( orb + L*( xyz + L*( anni2 + L*( orb + L * anni1 ))))))] += value; // + delta_lr delta_kp < E_jq E_it >
984 four_rdm[ crea2 + L*( xyz + L*( crea1 + L*( orb + L*( xyz + L*( anni2 + L*( orb + L * anni1 ))))))] += value; // + delta_lr delta_jp < E_iq E_kt >
985 four_rdm[ crea2 + L*( xyz + L*( crea1 + L*( orb + L*( anni2 + L*( orb + L*( xyz + L * anni1 ))))))] += value; // + delta_lq delta_jr < E_ip E_kt >
986 four_rdm[ crea2 + L*( crea1 + L*( xyz + L*( orb + L*( xyz + L*( orb + L*( anni2 + L * anni1 ))))))] += value; // + delta_lq delta_kp < E_ir E_jt >
987 four_rdm[ crea1 + L*( xyz + L*( crea2 + L*( orb + L*( xyz + L*( orb + L*( anni2 + L * anni1 ))))))] += value; // + delta_lq delta_jp < E_kr E_it >
988 four_rdm[ crea1 + L*( crea2 + L*( xyz + L*( orb + L*( orb + L*( xyz + L*( anni2 + L * anni1 ))))))] += value; // + delta_lp delta_kq < E_jr E_it >
989 four_rdm[ xyz + L*( crea1 + L*( crea2 + L*( orb + L*( orb + L*( xyz + L*( anni2 + L * anni1 ))))))] += value; // + delta_lp delta_iq < E_kr E_jt >
990 four_rdm[ xyz + L*( crea2 + L*( crea1 + L*( orb + L*( orb + L*( anni2 + L*( xyz + L * anni1 ))))))] += value; // + delta_lp delta_ir < E_jq E_kt >
991 four_rdm[ crea2 + L*( xyz + L*( orb + L*( crea1 + L*( xyz + L*( orb + L*( anni2 + L * anni1 ))))))] += value; // + delta_kq delta_jp < E_ir E_lt >
992 four_rdm[ crea2 + L*( xyz + L*( orb + L*( crea1 + L*( orb + L*( anni2 + L*( xyz + L * anni1 ))))))] += value; // + delta_kp delta_jr < E_iq E_lt >
993 }
994 }
995
996 }
997
998 for ( unsigned int crea3 = anni1; crea3 < L; crea3++ ){ // crea3 = j >= t = anni1
999 for ( unsigned int anni3 = (( anni1 == anni2 ) ? anni1 + 1 : anni1 ); anni3 < L; anni3++ ){ // anni3 = q >= t
1000
1001 const int irrep_center3 = Irreps::directProd( getOrb2Irrep( crea3 ), getOrb2Irrep( anni3 ) );
1002 const int target_irrep3 = Irreps::directProd( target_irrep2, irrep_center3 );
1003 const int irrep_center4 = Irreps::directProd( Irreps::directProd( irrep_center1 , irrep_center2 ), irrep_center3 );
1004 apply_excitation( workspace2, workspace3, crea3, anni3, target_irrep2 );
1005
1006 if ( irrep_center4 == 0 ){
1007
1008 // value = < E_{crea3,anni3} E_{crea2,anni2} E_{crea1,anni1} >
1009 const double value = FCIddot( orig_length, workspace3, vector );
1010 for ( unsigned int orb = anni1; orb < L; orb++ ){
1011 // i j k l p q r t
1012 four_rdm[ crea3 + L*( crea2 + L*( crea1 + L*( orb + L*( anni3 + L*( anni2 + L*( orb + L * anni1 ))))))] -= value; // - delta_lr < E_ip E_jq E_kt >
1013 four_rdm[ crea3 + L*( crea1 + L*( crea2 + L*( orb + L*( anni3 + L*( orb + L*( anni2 + L * anni1 ))))))] -= value; // - delta_lq < E_ip E_kr E_jt >
1014 four_rdm[ crea1 + L*( crea3 + L*( crea2 + L*( orb + L*( orb + L*( anni3 + L*( anni2 + L * anni1 ))))))] -= value; // - delta_lp < E_jq E_kr E_it >
1015 four_rdm[ crea3 + L*( crea2 + L*( orb + L*( crea1 + L*( anni3 + L*( orb + L*( anni2 + L * anni1 ))))))] -= value; // - delta_kq < E_ip E_jr E_lt >
1016 four_rdm[ crea3 + L*( crea2 + L*( orb + L*( crea1 + L*( orb + L*( anni2 + L*( anni3 + L * anni1 ))))))] -= value; // - delta_kp < E_ir E_jq E_lt >
1017 four_rdm[ crea3 + L*( orb + L*( crea2 + L*( crea1 + L*( orb + L*( anni3 + L*( anni2 + L * anni1 ))))))] -= value; // - delta_jp < E_iq E_kr E_lt >
1018 }
1019
1020 }
1021
1022 if ( crea3 >= (( crea1 == crea2 ) ? crea2 + 1 : crea2 ) ){ // crea3 = j >= k = crea2 >= t = anni1
1023
1024 for ( unsigned int crea4 = (( crea2 == crea3 ) ? crea3 + 1 : crea3 ); crea4 < L; crea4++ ){ // crea4 = i >= j = crea3 >= k = crea2 >= t = anni1
1025 for ( unsigned int anni4 = (( anni1 == anni2 ) ? anni1 + 1 : anni1 ); anni4 < L; anni4++ ){ // anni4 = p >= t
1026
1027 if ( (( anni2 == anni3 ) && ( anni3 == anni4 )) == false ){
1028
1029 const int irrep_product4 = Irreps::directProd( getOrb2Irrep( crea4 ), getOrb2Irrep( anni4 ) );
1030 if ( irrep_product4 == irrep_center4 ){
1031
1032 apply_excitation( workspace3, workspace4, crea4, anni4, target_irrep3 );
1033 // value = < E_{crea4,anni4} E_{crea3,anni3} E_{crea2,anni2} E_{crea1,anni1} >
1034 const double value = FCIddot( orig_length, workspace4, vector );
1035 four_rdm[ crea4 + L*( crea3 + L*( crea2 + L*( crea1 + L*( anni4 + L*( anni3 + L*( anni2 + L * anni1 ))))))] += value;
1036
1037 }
1038 }
1039 }
1040 }
1041 }
1042 }
1043 }
1044 }
1045 }
1046 }
1047 }
1048 delete [] workspace1;
1049 delete [] workspace2;
1050 delete [] workspace3;
1051 delete [] workspace4;
1052
1053 // Make 48-fold permutation symmetric
1054 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){ // anni1 = t
1055 for ( unsigned int crea1 = anni1; crea1 < L; crea1++ ){ // crea1 = l >= t = anni1
1056 const int irrep_center1 = Irreps::directProd( getOrb2Irrep( crea1 ), getOrb2Irrep( anni1 ) );
1057 for ( unsigned int crea2 = anni1; crea2 < L; crea2++ ){ // crea2 = k >= t = anni1
1058 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){ // anni2 = r >= t = anni1
1059 const int irrep_center2 = Irreps::directProd( getOrb2Irrep( crea2 ), getOrb2Irrep( anni2 ) );
1060 const int irrep_12 = Irreps::directProd( irrep_center1, irrep_center2 );
1061 for ( unsigned int crea3 = crea2; crea3 < L; crea3++ ){ // crea3 = j >= k = crea2 >= t = anni1
1062 for ( unsigned int anni3 = anni1; anni3 < L; anni3++ ){ // anni3 = q >= t = anni1
1063 const int irrep_center3 = Irreps::directProd( getOrb2Irrep( crea3 ), getOrb2Irrep( anni3 ) );
1064 const int irrep_123 = Irreps::directProd( irrep_12, irrep_center3 );
1065 for ( unsigned int crea4 = crea3; crea4 < L; crea4++ ){ // crea4 = i >= j = crea3 >= k = crea2 >= t = anni1
1066 for ( unsigned int anni4 = anni1; anni4 < L; anni4++ ){ // anni4 = p >= t = anni1
1067 const int irrep_center4 = Irreps::directProd( getOrb2Irrep( crea4 ), getOrb2Irrep( anni4 ) );
1068 if ( irrep_123 == irrep_center4 ){
1069
1070 /* crea4 >= crea3 >= crea2 >= anni1
1071 crea1, anni2, anni3, anni4 >= anni1 */
1072
1073 const double value = four_rdm[ crea4 + L*( crea3 + L*( crea2 + L*( crea1 + L*( anni4 + L*( anni3 + L*( anni2 + L * anni1 )))))) ];
1074 four_rdm[ crea4 + L*( crea2 + L*( crea1 + L*( crea3 + L*( anni4 + L*( anni2 + L*( anni1 + L * anni3 )))))) ] = value;
1075 four_rdm[ crea4 + L*( crea1 + L*( crea3 + L*( crea2 + L*( anni4 + L*( anni1 + L*( anni3 + L * anni2 )))))) ] = value;
1076 four_rdm[ crea4 + L*( crea1 + L*( crea2 + L*( crea3 + L*( anni4 + L*( anni1 + L*( anni2 + L * anni3 )))))) ] = value;
1077 four_rdm[ crea4 + L*( crea2 + L*( crea3 + L*( crea1 + L*( anni4 + L*( anni2 + L*( anni3 + L * anni1 )))))) ] = value;
1078 four_rdm[ crea4 + L*( crea3 + L*( crea1 + L*( crea2 + L*( anni4 + L*( anni3 + L*( anni1 + L * anni2 )))))) ] = value;
1079
1080 four_rdm[ crea3 + L*( crea4 + L*( crea2 + L*( crea1 + L*( anni3 + L*( anni4 + L*( anni2 + L * anni1 )))))) ] = value;
1081 four_rdm[ crea3 + L*( crea2 + L*( crea1 + L*( crea4 + L*( anni3 + L*( anni2 + L*( anni1 + L * anni4 )))))) ] = value;
1082 four_rdm[ crea3 + L*( crea1 + L*( crea4 + L*( crea2 + L*( anni3 + L*( anni1 + L*( anni4 + L * anni2 )))))) ] = value;
1083 four_rdm[ crea3 + L*( crea1 + L*( crea2 + L*( crea4 + L*( anni3 + L*( anni1 + L*( anni2 + L * anni4 )))))) ] = value;
1084 four_rdm[ crea3 + L*( crea2 + L*( crea4 + L*( crea1 + L*( anni3 + L*( anni2 + L*( anni4 + L * anni1 )))))) ] = value;
1085 four_rdm[ crea3 + L*( crea4 + L*( crea1 + L*( crea2 + L*( anni3 + L*( anni4 + L*( anni1 + L * anni2 )))))) ] = value;
1086
1087 four_rdm[ crea2 + L*( crea3 + L*( crea4 + L*( crea1 + L*( anni2 + L*( anni3 + L*( anni4 + L * anni1 )))))) ] = value;
1088 four_rdm[ crea2 + L*( crea4 + L*( crea1 + L*( crea3 + L*( anni2 + L*( anni4 + L*( anni1 + L * anni3 )))))) ] = value;
1089 four_rdm[ crea2 + L*( crea1 + L*( crea3 + L*( crea4 + L*( anni2 + L*( anni1 + L*( anni3 + L * anni4 )))))) ] = value;
1090 four_rdm[ crea2 + L*( crea1 + L*( crea4 + L*( crea3 + L*( anni2 + L*( anni1 + L*( anni4 + L * anni3 )))))) ] = value;
1091 four_rdm[ crea2 + L*( crea4 + L*( crea3 + L*( crea1 + L*( anni2 + L*( anni4 + L*( anni3 + L * anni1 )))))) ] = value;
1092 four_rdm[ crea2 + L*( crea3 + L*( crea1 + L*( crea4 + L*( anni2 + L*( anni3 + L*( anni1 + L * anni4 )))))) ] = value;
1093
1094 four_rdm[ crea1 + L*( crea3 + L*( crea2 + L*( crea4 + L*( anni1 + L*( anni3 + L*( anni2 + L * anni4 )))))) ] = value;
1095 four_rdm[ crea1 + L*( crea2 + L*( crea4 + L*( crea3 + L*( anni1 + L*( anni2 + L*( anni4 + L * anni3 )))))) ] = value;
1096 four_rdm[ crea1 + L*( crea4 + L*( crea3 + L*( crea2 + L*( anni1 + L*( anni4 + L*( anni3 + L * anni2 )))))) ] = value;
1097 four_rdm[ crea1 + L*( crea4 + L*( crea2 + L*( crea3 + L*( anni1 + L*( anni4 + L*( anni2 + L * anni3 )))))) ] = value;
1098 four_rdm[ crea1 + L*( crea2 + L*( crea3 + L*( crea4 + L*( anni1 + L*( anni2 + L*( anni3 + L * anni4 )))))) ] = value;
1099 four_rdm[ crea1 + L*( crea3 + L*( crea4 + L*( crea2 + L*( anni1 + L*( anni3 + L*( anni4 + L * anni2 )))))) ] = value;
1100
1101 four_rdm[ anni4 + L*( anni3 + L*( anni2 + L*( anni1 + L*( crea4 + L*( crea3 + L*( crea2 + L * crea1 )))))) ] = value;
1102 four_rdm[ anni4 + L*( anni2 + L*( anni1 + L*( anni3 + L*( crea4 + L*( crea2 + L*( crea1 + L * crea3 )))))) ] = value;
1103 four_rdm[ anni4 + L*( anni1 + L*( anni3 + L*( anni2 + L*( crea4 + L*( crea1 + L*( crea3 + L * crea2 )))))) ] = value;
1104 four_rdm[ anni4 + L*( anni1 + L*( anni2 + L*( anni3 + L*( crea4 + L*( crea1 + L*( crea2 + L * crea3 )))))) ] = value;
1105 four_rdm[ anni4 + L*( anni2 + L*( anni3 + L*( anni1 + L*( crea4 + L*( crea2 + L*( crea3 + L * crea1 )))))) ] = value;
1106 four_rdm[ anni4 + L*( anni3 + L*( anni1 + L*( anni2 + L*( crea4 + L*( crea3 + L*( crea1 + L * crea2 )))))) ] = value;
1107
1108 four_rdm[ anni3 + L*( anni4 + L*( anni2 + L*( anni1 + L*( crea3 + L*( crea4 + L*( crea2 + L * crea1 )))))) ] = value;
1109 four_rdm[ anni3 + L*( anni2 + L*( anni1 + L*( anni4 + L*( crea3 + L*( crea2 + L*( crea1 + L * crea4 )))))) ] = value;
1110 four_rdm[ anni3 + L*( anni1 + L*( anni4 + L*( anni2 + L*( crea3 + L*( crea1 + L*( crea4 + L * crea2 )))))) ] = value;
1111 four_rdm[ anni3 + L*( anni1 + L*( anni2 + L*( anni4 + L*( crea3 + L*( crea1 + L*( crea2 + L * crea4 )))))) ] = value;
1112 four_rdm[ anni3 + L*( anni2 + L*( anni4 + L*( anni1 + L*( crea3 + L*( crea2 + L*( crea4 + L * crea1 )))))) ] = value;
1113 four_rdm[ anni3 + L*( anni4 + L*( anni1 + L*( anni2 + L*( crea3 + L*( crea4 + L*( crea1 + L * crea2 )))))) ] = value;
1114
1115 four_rdm[ anni2 + L*( anni3 + L*( anni4 + L*( anni1 + L*( crea2 + L*( crea3 + L*( crea4 + L * crea1 )))))) ] = value;
1116 four_rdm[ anni2 + L*( anni4 + L*( anni1 + L*( anni3 + L*( crea2 + L*( crea4 + L*( crea1 + L * crea3 )))))) ] = value;
1117 four_rdm[ anni2 + L*( anni1 + L*( anni3 + L*( anni4 + L*( crea2 + L*( crea1 + L*( crea3 + L * crea4 )))))) ] = value;
1118 four_rdm[ anni2 + L*( anni1 + L*( anni4 + L*( anni3 + L*( crea2 + L*( crea1 + L*( crea4 + L * crea3 )))))) ] = value;
1119 four_rdm[ anni2 + L*( anni4 + L*( anni3 + L*( anni1 + L*( crea2 + L*( crea4 + L*( crea3 + L * crea1 )))))) ] = value;
1120 four_rdm[ anni2 + L*( anni3 + L*( anni1 + L*( anni4 + L*( crea2 + L*( crea3 + L*( crea1 + L * crea4 )))))) ] = value;
1121
1122 four_rdm[ anni1 + L*( anni3 + L*( anni2 + L*( anni4 + L*( crea1 + L*( crea3 + L*( crea2 + L * crea4 )))))) ] = value;
1123 four_rdm[ anni1 + L*( anni2 + L*( anni4 + L*( anni3 + L*( crea1 + L*( crea2 + L*( crea4 + L * crea3 )))))) ] = value;
1124 four_rdm[ anni1 + L*( anni4 + L*( anni3 + L*( anni2 + L*( crea1 + L*( crea4 + L*( crea3 + L * crea2 )))))) ] = value;
1125 four_rdm[ anni1 + L*( anni4 + L*( anni2 + L*( anni3 + L*( crea1 + L*( crea4 + L*( crea2 + L * crea3 )))))) ] = value;
1126 four_rdm[ anni1 + L*( anni2 + L*( anni3 + L*( anni4 + L*( crea1 + L*( crea2 + L*( crea3 + L * crea4 )))))) ] = value;
1127 four_rdm[ anni1 + L*( anni3 + L*( anni4 + L*( anni2 + L*( crea1 + L*( crea3 + L*( crea4 + L * crea2 )))))) ] = value;
1128
1129 }
1130 }
1131 }
1132 }
1133 }
1134 }
1135 }
1136 }
1137 }
1138
1139 for ( unsigned int ann = 0; ann < L; ann++ ){
1140 for ( unsigned int orb = 0; orb < L; orb++ ){
1141 for ( unsigned int combo = 0; combo < L*L*L*L; combo++ ){
1142 four_rdm[ combo + L*L*L*L*( ann + L * ( ann + L * ( ann + L * orb ))) ] = 0.0;
1143 four_rdm[ combo + L*L*L*L*( ann + L * ( ann + L * ( orb + L * ann ))) ] = 0.0;
1144 four_rdm[ combo + L*L*L*L*( ann + L * ( orb + L * ( ann + L * ann ))) ] = 0.0;
1145 four_rdm[ combo + L*L*L*L*( orb + L * ( ann + L * ( ann + L * ann ))) ] = 0.0;
1146 four_rdm[ L*L*L*L * combo + ann + L * ( ann + L * ( ann + L * orb )) ] = 0.0;
1147 four_rdm[ L*L*L*L * combo + ann + L * ( ann + L * ( orb + L * ann )) ] = 0.0;
1148 four_rdm[ L*L*L*L * combo + ann + L * ( orb + L * ( ann + L * ann )) ] = 0.0;
1149 four_rdm[ L*L*L*L * combo + orb + L * ( ann + L * ( ann + L * ann )) ] = 0.0;
1150 }
1151 }
1152 }
1153
1154 gettimeofday(&end, NULL);
1155 const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
1156 if ( FCIverbose > 0 ){ cout << "FCI::Fill4RDM : Wall time = " << elapsed << " seconds" << endl; }
1157
1158 }
1159
Fock4RDM(double * vector,double * three_rdm,double * fock,double * output) const1160 void CheMPS2::FCI::Fock4RDM( double * vector, double * three_rdm, double * fock, double * output ) const{
1161
1162 assert( Nel_up + Nel_down >= 4 );
1163 const double elapsed = Driver3RDM( vector, output, three_rdm, fock, L + 1 );
1164 if ( FCIverbose > 0 ){ cout << "FCI::Fock4RDM : Wall time = " << elapsed << " seconds" << endl; }
1165
1166 }
1167
Fill3RDM(double * vector,double * output) const1168 void CheMPS2::FCI::Fill3RDM( double * vector, double * output ) const{
1169
1170 assert( Nel_up + Nel_down >= 3 );
1171 const double elapsed = Driver3RDM( vector, output, NULL, NULL, L + 1 );
1172 if ( FCIverbose > 0 ){ cout << "FCI::Fill3RDM : Wall time = " << elapsed << " seconds" << endl; }
1173
1174 }
1175
Diag4RDM(double * vector,double * three_rdm,const unsigned int orbz,double * output) const1176 void CheMPS2::FCI::Diag4RDM( double * vector, double * three_rdm, const unsigned int orbz, double * output ) const{
1177
1178 assert( Nel_up + Nel_down >= 4 );
1179 const double elapsed = Driver3RDM( vector, output, three_rdm, NULL, orbz );
1180 if ( FCIverbose > 0 ){ cout << "FCI::Diag4RDM : Wall time = " << elapsed << " seconds" << endl; }
1181
1182 }
1183
Driver3RDM(double * vector,double * output,double * three_rdm,double * fock,const unsigned int orbz) const1184 double CheMPS2::FCI::Driver3RDM( double * vector, double * output, double * three_rdm, double * fock, const unsigned int orbz ) const{
1185
1186 struct timeval start, end;
1187 gettimeofday(&start, NULL);
1188
1189 ClearVector( L*L*L*L*L*L, output );
1190 const unsigned int orig_length = getVecLength( 0 );
1191 unsigned int max_length = getVecLength( 0 );
1192 for ( unsigned int irrep = 1; irrep < num_irreps; irrep++ ){
1193 if ( getVecLength( irrep ) > max_length ){ max_length = getVecLength( irrep ); }
1194 }
1195 double * workspace1 = new double[ max_length ];
1196 double * workspace2 = new double[ max_length ];
1197 double * workspace3 = new double[ orig_length ];
1198
1199 double * chi = NULL;
1200 const bool task_fock = ( fock != NULL );
1201 const bool task_E_zz = (( fock == NULL ) && ( three_rdm != NULL ));
1202 const bool task_3rdm = (( fock == NULL ) && ( three_rdm == NULL ));
1203
1204 if ( task_3rdm ){
1205
1206 assert( orbz == L + 1 );
1207 chi = vector; // | Chi > = | 0 >
1208
1209 /* Calculate the 3-RDM:
1210 Gamma_{ijk,pqr} = < 0 | E_ip E_jq E_kr | Chi >
1211 - delta_kq < 0 | E_ip E_jr | Chi >
1212 - delta_kp < 0 | E_ir E_jq | Chi >
1213 - delta_jp < 0 | E_iq E_kr | Chi >
1214 + delta_kq delta_jp < 0 | E_ir | Chi >
1215 + delta_kp delta_jr < 0 | E_iq | Chi >
1216 */
1217 }
1218
1219 if ( task_E_zz ){
1220
1221 assert( orbz < L );
1222 chi = new double[ orig_length ];
1223 apply_excitation( vector, chi, orbz, orbz, TargetIrrep ); // | Chi > = E_zz | 0 >
1224
1225 /* Calculate the 4-RDM elements with orbital z fixed:
1226 Gamma_{ijkz,pqrz} = < 0 | E_ip E_jq E_kr | Chi >
1227 - delta_kq < 0 | E_ip E_jr | Chi >
1228 - delta_kp < 0 | E_ir E_jq | Chi >
1229 - delta_jp < 0 | E_iq E_kr | Chi >
1230 + delta_kq delta_jp < 0 | E_ir | Chi >
1231 + delta_kp delta_jr < 0 | E_iq | Chi >
1232 - ( delta_pz + delta_qz + delta_rz ) Gamma_{ijk,pqr}
1233 */
1234 }
1235
1236 if ( task_fock ){
1237
1238 assert( orbz == L + 1 );
1239 chi = new double[ orig_length ];
1240 ClearVector( orig_length, chi );
1241 for ( unsigned int anni = 0; anni < L; anni++ ){
1242 for ( unsigned int crea = 0; crea < L; crea++ ){
1243 if ( getOrb2Irrep( crea ) == getOrb2Irrep( anni ) ){
1244 apply_excitation( vector, workspace1, crea, anni, TargetIrrep );
1245 FCIdaxpy( orig_length, fock[ crea + L * anni ], workspace1, chi ); // | Chi > = sum_{l,t} fock[l,t] E_lt | 0 >
1246 }
1247 }
1248 }
1249
1250 /* Calculate the 4-RDM contracted with the Fock operator:
1251 sum_{l,t} fock[l,t] Gamma_{ijkl,pqrt} = < 0 | E_ip E_jq E_kr | Chi >
1252 - delta_kq < 0 | E_ip E_jr | Chi >
1253 - delta_kp < 0 | E_ir E_jq | Chi >
1254 - delta_jp < 0 | E_iq E_kr | Chi >
1255 + delta_kq delta_jp < 0 | E_ir | Chi >
1256 + delta_kp delta_jr < 0 | E_iq | Chi >
1257 - sum_{t} ( fock[r,t] Gamma_{ijk,pqt}
1258 + fock[q,t] Gamma_{ijk,ptr}
1259 + fock[p,t] Gamma_{ijk,tqr} )
1260 */
1261 }
1262
1263 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){ // anni1 = i ( works in on the bra ) ( smaller than j, k )
1264 for ( unsigned int crea1 = 0; crea1 < L; crea1++ ){ // crea1 = p ( can be anything )
1265
1266 const int irrep_center1 = Irreps::directProd( getOrb2Irrep( crea1 ), getOrb2Irrep( anni1 ) );
1267 const int target_irrep1 = Irreps::directProd( TargetIrrep, irrep_center1 );
1268 apply_excitation( vector, workspace1, crea1, anni1, TargetIrrep );
1269
1270 if ( irrep_center1 == 0 ){
1271
1272 // value = < Chi | E_{crea1,anni1} | 0 >
1273 const double value = FCIddot( orig_length, workspace1, chi );
1274 for ( unsigned int j = anni1; j < L; j++ ){
1275 for ( unsigned int k = anni1; k < L; k++ ){
1276 // i j k p q r
1277 output[ anni1 + L*( j + L*( k + L*( j + L*( k + L * crea1 )))) ] += value; // + delta_kq delta_jp < 0 | E_ir | Chi >
1278 output[ anni1 + L*( j + L*( k + L*( k + L*( crea1 + L * j )))) ] += value; // + delta_kp delta_jr < 0 | E_iq | Chi >
1279 }
1280 }
1281
1282 }
1283
1284 for ( unsigned int crea2 = 0; crea2 < L; crea2++ ){ // crea2 = q
1285 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){ // anni2 = j >= ( i = anni1 )
1286
1287 const int irrep_center2 = Irreps::directProd( getOrb2Irrep( crea2 ), getOrb2Irrep( anni2 ) );
1288 const int target_irrep2 = Irreps::directProd( target_irrep1, irrep_center2 );
1289 const int irrep_center3 = Irreps::directProd( irrep_center1, irrep_center2 );
1290 apply_excitation( workspace1, workspace2, crea2, anni2, target_irrep1 );
1291
1292 if ( irrep_center1 == irrep_center2 ){
1293
1294 // value = < Chi | E_{crea2,anni2} E_{crea1,anni1} | 0 >
1295 const double value = FCIddot( orig_length, workspace2, chi );
1296 for ( unsigned int orb = anni1; orb < L; orb++ ){
1297 // i j k p q r
1298 output[ anni1 + L*( anni2 + L*( orb + L*( crea1 + L*( orb + L * crea2 )))) ] -= value; // - delta_kq < 0 | E_ip E_jr | Chi >
1299 output[ anni1 + L*( anni2 + L*( orb + L*( orb + L*( crea2 + L * crea1 )))) ] -= value; // - delta_kp < 0 | E_ir E_jq | Chi >
1300 output[ anni1 + L*( orb + L*( anni2 + L*( orb + L*( crea1 + L * crea2 )))) ] -= value; // - delta_jp < 0 | E_iq E_kr | Chi >
1301 }
1302
1303 }
1304
1305 if (( crea1 >= anni1 ) && ( crea2 >= anni1 )){
1306
1307 for ( unsigned int crea3 = (( crea1 == crea2 ) ? crea2 + 1 : crea2 ); crea3 < L; crea3++ ){ // crea3 = r >= ( q = crea2 ) >= ( i = anni1 )
1308 for ( unsigned int anni3 = (( anni1 == anni2 ) ? anni1 + 1 : anni1 ); anni3 < L; anni3++ ){ // anni3 = k >= ( i = anni1 )
1309
1310 const int irrep_product3 = Irreps::directProd( getOrb2Irrep( crea3 ), getOrb2Irrep( anni3 ) );
1311
1312 if ( irrep_center3 == irrep_product3 ){ // I1 x I2 x I3 = Itrivial
1313
1314 apply_excitation( workspace2, workspace3, crea3, anni3, target_irrep2 );
1315
1316 // value = < Chi | E_{crea3,anni3} E_{crea2,anni2} E_{crea1,anni1} | 0 >
1317 double value = FCIddot( orig_length, workspace3, chi );
1318
1319 if ( task_fock ){
1320 for ( unsigned int t = 0; t < L; t++ ){
1321 // Irrep diagonality of fock is checked by three_rdm values being zero
1322 value -= ( fock[ crea3 + L * t ] * three_rdm[ anni1 + L*( anni2 + L*( anni3 + L*( crea1 + L*( crea2 + L * t )))) ]
1323 + fock[ crea2 + L * t ] * three_rdm[ anni1 + L*( anni2 + L*( anni3 + L*( crea1 + L*( t + L * crea3 )))) ]
1324 + fock[ crea1 + L * t ] * three_rdm[ anni1 + L*( anni2 + L*( anni3 + L*( t + L*( crea2 + L * crea3 )))) ] );
1325 }
1326 }
1327
1328 if ( task_E_zz ){
1329 const int number = (( orbz == crea1 ) ? 1 : 0 ) + (( orbz == crea2 ) ? 1 : 0 ) + (( orbz == crea3 ) ? 1 : 0 );
1330 if ( number > 0 ){
1331 value -= number * three_rdm[ anni1 + L*( anni2 + L*( anni3 + L*( crea1 + L*( crea2 + L * crea3 )))) ];
1332 }
1333 }
1334
1335 output[ anni1 + L*( anni2 + L*( anni3 + L*( crea1 + L*( crea2 + L * crea3 )))) ] += value;
1336
1337 }
1338 }
1339 }
1340 }
1341 }
1342 }
1343 }
1344 }
1345 delete [] workspace1;
1346 delete [] workspace2;
1347 delete [] workspace3;
1348 if (( task_fock ) || ( task_E_zz )){ delete [] chi; }
1349
1350 // Make 12-fold permutation symmetric
1351 for ( unsigned int anni1 = 0; anni1 < L; anni1++ ){
1352 for ( unsigned int crea1 = anni1; crea1 < L; crea1++ ){
1353 const int irrep_prod1 = Irreps::directProd( getOrb2Irrep( crea1 ) , getOrb2Irrep( anni1 ) ); // Ic1 x Ia1
1354 for ( unsigned int crea2 = anni1; crea2 < L; crea2++ ){
1355 const int irrep_prod2 = Irreps::directProd( irrep_prod1 , getOrb2Irrep( crea2 ) ); // Ic1 x Ia1 x Ic2
1356 for ( unsigned int anni2 = anni1; anni2 < L; anni2++ ){
1357 const int irrep_prod3 = Irreps::directProd( irrep_prod2 , getOrb2Irrep( anni2 ) ); // Ic1 x Ia1 x Ic2 x Ia2
1358 for ( unsigned int crea3 = crea2; crea3 < L; crea3++ ){
1359 const int irrep_prod4 = Irreps::directProd( irrep_prod3 , getOrb2Irrep( crea3 ) ); // Ic1 x Ia1 x Ic2 x Ia2 x Ic3
1360 for ( unsigned int anni3 = anni1; anni3 < L; anni3++ ){
1361 if ( irrep_prod4 == getOrb2Irrep( anni3 )){ // Ic1 x Ia1 x Ic2 x Ia2 x Ic3 == Ia3
1362
1363 /* crea3 >= crea2 >= anni1
1364 crea1, anni3, anni2 >= anni1 */
1365
1366 const double value = output[ anni1 + L * ( anni2 + L * ( anni3 + L * ( crea1 + L * ( crea2 + L * crea3 ) ) ) ) ];
1367 output[ anni1 + L * ( anni3 + L * ( anni2 + L * ( crea1 + L * ( crea3 + L * crea2 ) ) ) ) ] = value;
1368
1369 output[ anni3 + L * ( anni2 + L * ( anni1 + L * ( crea3 + L * ( crea2 + L * crea1 ) ) ) ) ] = value;
1370 output[ anni2 + L * ( anni3 + L * ( anni1 + L * ( crea2 + L * ( crea3 + L * crea1 ) ) ) ) ] = value;
1371
1372 output[ anni2 + L * ( anni1 + L * ( anni3 + L * ( crea2 + L * ( crea1 + L * crea3 ) ) ) ) ] = value;
1373 output[ anni3 + L * ( anni1 + L * ( anni2 + L * ( crea3 + L * ( crea1 + L * crea2 ) ) ) ) ] = value;
1374
1375 output[ crea3 + L * ( crea2 + L * ( crea1 + L * ( anni3 + L * ( anni2 + L * anni1 ) ) ) ) ] = value;
1376 output[ crea2 + L * ( crea3 + L * ( crea1 + L * ( anni2 + L * ( anni3 + L * anni1 ) ) ) ) ] = value;
1377
1378 output[ crea2 + L * ( crea1 + L * ( crea3 + L * ( anni2 + L * ( anni1 + L * anni3 ) ) ) ) ] = value;
1379 output[ crea3 + L * ( crea1 + L * ( crea2 + L * ( anni3 + L * ( anni1 + L * anni2 ) ) ) ) ] = value;
1380
1381 output[ crea1 + L * ( crea3 + L * ( crea2 + L * ( anni1 + L * ( anni3 + L * anni2 ) ) ) ) ] = value;
1382 output[ crea1 + L * ( crea2 + L * ( crea3 + L * ( anni1 + L * ( anni2 + L * anni3 ) ) ) ) ] = value;
1383
1384 }
1385 }
1386 }
1387 }
1388 }
1389 }
1390 }
1391
1392 for ( unsigned int anni = 0; anni < L; anni++ ){
1393 for ( unsigned int combo = 0; combo < L*L*L; combo++ ){
1394 output[ combo + L * L * L * anni * ( 1 + L + L * L ) ] = 0.0;
1395 output[ anni * ( 1 + L + L * L ) + L * L * L * combo ] = 0.0;
1396 }
1397 }
1398
1399 gettimeofday(&end, NULL);
1400 const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
1401 return elapsed;
1402
1403 }
1404
CalcSpinSquared(double * vector) const1405 double CheMPS2::FCI::CalcSpinSquared(double * vector) const{
1406
1407 const unsigned int vecLength = getVecLength( 0 );
1408 double result = 0.0;
1409
1410 #pragma omp parallel for schedule(static) reduction(+:result)
1411 for ( unsigned int counter = 0; counter < vecLength; counter++ ){
1412 for ( unsigned int orbi = 0; orbi < L; orbi++ ){
1413
1414 const int irrep_up = getUpIrrepOfCounter( 0 , counter );
1415 const int irrep_down = Irreps::directProd( irrep_up , TargetIrrep );
1416 const int count_up = ( counter - irrep_center_jumps[ 0 ][ irrep_up ] ) % numPerIrrep_up[ irrep_up ];
1417 const int count_down = ( counter - irrep_center_jumps[ 0 ][ irrep_up ] ) / numPerIrrep_up[ irrep_up ];
1418
1419 // Diagonal terms
1420 const int diff_ii = lookup_sign_alpha[ irrep_up ][ orbi + L * orbi ][ count_up ]
1421 - lookup_sign_beta [ irrep_down ][ orbi + L * orbi ][ count_down ]; //Signed integers so subtracting is OK
1422 const double vector_at_counter_squared = vector[ counter ] * vector[ counter ];
1423 result += 0.75 * diff_ii * diff_ii * vector_at_counter_squared;
1424
1425 for ( unsigned int orbj = orbi+1; orbj < L; orbj++ ){
1426
1427 // Sz Sz
1428 const int diff_jj = lookup_sign_alpha[ irrep_up ][ orbj + L * orbj ][ count_up ]
1429 - lookup_sign_beta [ irrep_down ][ orbj + L * orbj ][ count_down ]; //Signed integers so subtracting is OK
1430 result += 0.5 * diff_ii * diff_jj * vector_at_counter_squared;
1431
1432 const int irrep_up_bis = Irreps::directProd( irrep_up , Irreps::directProd( getOrb2Irrep( orbi ) , getOrb2Irrep( orbj ) ) );
1433
1434 // - ( a_i,up^+ a_j,up )( a_j,down^+ a_i,down )
1435 const int sign_down_ji = lookup_sign_beta [ irrep_down ][ orbj + L * orbi ][ count_down ];
1436 const int sign_up_ij = lookup_sign_alpha[ irrep_up ][ orbi + L * orbj ][ count_up ];
1437 const int sign_product1 = sign_up_ij * sign_down_ji;
1438 if ( sign_product1 != 0 ){
1439 const int cnt_down_ji = lookup_cnt_beta [ irrep_down ][ orbj + L * orbi ][ count_down ];
1440 const int cnt_up_ij = lookup_cnt_alpha[ irrep_up ][ orbi + L * orbj ][ count_up ];
1441 result -= sign_product1 * vector[ irrep_center_jumps[ 0 ][ irrep_up_bis ] + cnt_up_ij + numPerIrrep_up[ irrep_up_bis ] * cnt_down_ji ] * vector[ counter ];
1442 }
1443
1444 // - ( a_j,up^+ a_i,up )( a_i,down^+ a_j,down )
1445 const int sign_down_ij = lookup_sign_beta [ irrep_down ][ orbi + L * orbj ][ count_down ];
1446 const int sign_up_ji = lookup_sign_alpha[ irrep_up ][ orbj + L * orbi ][ count_up ];
1447 const int sign_product2 = sign_up_ji * sign_down_ij;
1448 if ( sign_product2 != 0 ){
1449 const int cnt_down_ij = lookup_cnt_beta [ irrep_down ][ orbi + L * orbj ][ count_down ];
1450 const int cnt_up_ji = lookup_cnt_alpha[ irrep_up ][ orbj + L * orbi ][ count_up ];
1451 result -= sign_product2 * vector[ irrep_center_jumps[ 0 ][ irrep_up_bis ] + cnt_up_ji + numPerIrrep_up[ irrep_up_bis ] * cnt_down_ij ] * vector[ counter ];
1452 }
1453
1454 }
1455 }
1456 }
1457
1458 if ( FCIverbose > 0 ){
1459 const double intendedS = fabs( 0.5 * Nel_up - 0.5 * Nel_down ); // Be careful with subtracting unsigned integers...
1460 cout << "FCI::CalcSpinSquared : For intended spin " << intendedS
1461 << " the measured S(S+1) = " << result << " and intended S(S+1) = " << intendedS * (intendedS + 1.0) << endl;
1462 }
1463 return result;
1464
1465 }
1466
DiagHam(double * diag) const1467 void CheMPS2::FCI::DiagHam(double * diag) const{
1468
1469 const unsigned int vecLength = getVecLength( 0 );
1470
1471 #pragma omp parallel
1472 {
1473
1474 int * bits_up = new int[ L ];
1475 int * bits_down = new int[ L ];
1476
1477 #pragma omp for schedule(static)
1478 for ( unsigned int counter = 0; counter < vecLength; counter++ ){
1479
1480 double myResult = 0.0;
1481 getBitsOfCounter( 0 , counter , bits_up , bits_down ); // Fetch the corresponding bits
1482
1483 for ( unsigned int orb1 = 0; orb1 < L; orb1++ ){
1484 const int n_tot_orb1 = bits_up[ orb1 ] + bits_down[ orb1 ];
1485 myResult += n_tot_orb1 * getGmat( orb1 , orb1 );
1486 for ( unsigned int orb2 = 0; orb2 < L; orb2++ ){
1487 myResult += 0.5 * n_tot_orb1 * ( bits_up[ orb2 ] + bits_down[ orb2 ] ) * getERI( orb1 , orb1 , orb2 , orb2 );
1488 myResult += 0.5 * ( n_tot_orb1 - bits_up[ orb1 ] * bits_up[ orb2 ] - bits_down[ orb1 ] * bits_down[ orb2 ] ) * getERI( orb1 , orb2 , orb2 , orb1 );
1489 }
1490 }
1491
1492 diag[ counter ] = myResult;
1493
1494 }
1495
1496 delete [] bits_up;
1497 delete [] bits_down;
1498
1499 }
1500
1501 }
1502
1503
DiagHamSquared(double * output) const1504 void CheMPS2::FCI::DiagHamSquared(double * output) const{
1505
1506 struct timeval start, end;
1507 gettimeofday(&start, NULL);
1508
1509 const unsigned int vecLength = getVecLength( 0 );
1510
1511 //
1512 // Wick's theorem to evaluate the Hamiltonian squared:
1513 //
1514 // H = g_ij E_ij + 0.5 * (ij|kl) E_ij E_kl
1515 //
1516 // H^2 = g_ij g_kl E_ij E_kl
1517 // + 0.5 * [ g_ab (ij|kl) + (ab|ij) g_kl ] E_ab E_ij E_kl
1518 // + 0.25 * (ab|cd) * (ij|kl) E_ab E_cd E_ij E_kl
1519 //
1520 // Short illustration of what is being done:
1521 //
1522 // E_ij E_kl = (i,s1)^+ (j,s1) (k,s2)^+ (l,s2)
1523 //
1524 // = (i,s1)^+ (j,s1) (k,s2)^+ (l,s2)
1525 // | | | |
1526 // ---------- ----------
1527 // + (i,s1)^+ (j,s1) (k,s2)^+ (l,s2)
1528 // | | | |
1529 // | -------- |
1530 // --------------------------
1531 // = num(i,s1) delta(i,j) num(k,s2) delta(k,l)
1532 // + num(i,s1) delta(i,l) delta(s1,s2) [1-num(k,s2)] delta(k,j)
1533 //
1534 // g_ij g_kl E_ij E_kl = g_ii g_kk num(i,s1) num(k,s2) + g_ik g_ki num(i,s1) [1-num(k,s2)] delta(s1,s2)
1535 //
1536
1537 #pragma omp parallel
1538 {
1539
1540 int * bits_up = new int[ L ];
1541 int * bits_down = new int[ L ];
1542
1543 double * Jmat = new double[ L * L ]; // (ij|kk)( n_k,up + n_k,down )
1544 double * K_reg_up = new double[ L * L ]; // (ik|kj)( n_k,up )
1545 double * K_reg_down = new double[ L * L ]; // (ik|kj)( n_k,down )
1546 double * K_bar_up = new double[ L * L ]; // (ik|kj)( 1 - n_k,up )
1547 double * K_bar_down = new double[ L * L ]; // (ik|kj)( 1 - n_k,down )
1548
1549 int * specific_orbs_irrep = new int[ num_irreps * ( L + 1 ) ];
1550 for ( unsigned int irrep = 0; irrep < num_irreps; irrep++ ){
1551 int count = 0;
1552 for ( unsigned int orb = 0; orb < L; orb++){
1553 specific_orbs_irrep[ orb + ( L + 1 ) * irrep ] = 0;
1554 if ( getOrb2Irrep(orb) == irrep ){
1555 specific_orbs_irrep[ count + ( L + 1 ) * irrep ] = orb;
1556 count++;
1557 }
1558 }
1559 specific_orbs_irrep[ L + ( L + 1 ) * irrep ] = count;
1560 }
1561
1562 #pragma omp for schedule(static)
1563 for (unsigned int counter = 0; counter < vecLength; counter++){
1564
1565 getBitsOfCounter( 0 , counter , bits_up , bits_down ); // Fetch the corresponding bits
1566
1567 // Construct the J and K matrices properly
1568 for ( unsigned int i = 0; i < L; i++ ){
1569 for ( unsigned int j = i; j < L; j++ ){
1570
1571 double val_J = 0.0;
1572 double val_KregUP = 0.0;
1573 double val_KregDOWN = 0.0;
1574 double val_KbarUP = 0.0;
1575 double val_KbarDOWN = 0.0;
1576
1577 if ( getOrb2Irrep(i) == getOrb2Irrep(j) ){
1578 for ( unsigned int k = 0; k < L; k++ ){
1579 const double temp = getERI(i,k,k,j);
1580 val_J += getERI(i,j,k,k) * ( bits_up[k] + bits_down[k] );
1581 val_KregUP += temp * bits_up[k];
1582 val_KregDOWN += temp * bits_down[k];
1583 val_KbarUP += temp * ( 1 - bits_up[k] );
1584 val_KbarDOWN += temp * ( 1 - bits_down[k] );
1585 }
1586 }
1587
1588 Jmat[ i + L * j ] = val_J;
1589 Jmat[ j + L * i ] = val_J;
1590 K_reg_up[ i + L * j ] = val_KregUP;
1591 K_reg_up[ j + L * i ] = val_KregUP;
1592 K_reg_down[ i + L * j ] = val_KregDOWN;
1593 K_reg_down[ j + L * i ] = val_KregDOWN;
1594 K_bar_up[ i + L * j ] = val_KbarUP;
1595 K_bar_up[ j + L * i ] = val_KbarUP;
1596 K_bar_down[ i + L * j ] = val_KbarDOWN;
1597 K_bar_down[ j + L * i ] = val_KbarDOWN;
1598
1599 }
1600 }
1601
1602 double temp = 0.0;
1603 // G[i,i] (n_i,up + n_i,down) + 0.5 * ( J[i,i] (n_i,up + n_i,down) + K_bar_up[i,i] * n_i,up + K_bar_down[i,i] * n_i,down )
1604 for ( unsigned int i = 0; i < L; i++ ){
1605 const int num_i = bits_up[i] + bits_down[i];
1606 temp += getGmat(i, i) * num_i + 0.5 * ( Jmat[ i + L * i ] * num_i + K_bar_up[ i + L * i ] * bits_up[i]
1607 + K_bar_down[ i + L * i ] * bits_down[i] );
1608 }
1609 double myResult = temp*temp;
1610
1611 for ( unsigned int p = 0; p < L; p++ ){
1612 for ( unsigned int q = 0; q < L; q++ ){
1613 if ( getOrb2Irrep(p) == getOrb2Irrep(q) ){
1614
1615 const int special_pq = bits_up[p] * ( 1 - bits_up[q] ) + bits_down[p] * ( 1 - bits_down[q] );
1616 const double GplusJ_pq = getGmat(p, q) + Jmat[ p + L * q ];
1617 const double K_cross_pq_up = ( K_bar_up[ p + L * q ] - K_reg_up[ p + L * q ] ) * bits_up[p] * ( 1 - bits_up[q] );
1618 const double K_cross_pq_down = ( K_bar_down[ p + L * q ] - K_reg_down[ p + L * q ] ) * bits_down[p] * ( 1 - bits_down[q] );
1619
1620 myResult += ( GplusJ_pq * ( special_pq * GplusJ_pq + K_cross_pq_up + K_cross_pq_down )
1621 + 0.25 * ( K_cross_pq_up * K_cross_pq_up + K_cross_pq_down * K_cross_pq_down ) );
1622
1623 }
1624 }
1625 }
1626
1627 /*
1628
1629 Part which can be optimized most still --> For H2O/6-31G takes 82.8 % of time with Intel compiler
1630 Optimization can in principle be done with lookup tables + dgemm_ as matvec product --> bit of work :-(
1631
1632 For future reference: the quantity which is computed here:
1633
1634 0.5 * (ak|ci) * (ak|ci) * [ n_a,up * (1-n_k,up) + n_a,down * (1-n_k,down) ] * [ n_c,up * (1-n_i,up) + n_c,down * (1-n_i,down) ]
1635 - 0.5 * (ak|ci) * (ai|ck) * [ n_a,up * n_c,up * (1-n_i,up) * (1-n_k,up) + n_a,down * n_c,down * (1-n_i,down) * (1-n_k,down) ]
1636
1637 */
1638 for ( unsigned int k = 0; k < L; k++ ){
1639 if ( bits_up[k] + bits_down[k] < 2 ){
1640 for ( unsigned int a = 0; a < L; a++ ){
1641
1642 const int special_ak = ( bits_up[a] * ( 1 - bits_up[k] )
1643 + bits_down[a] * ( 1 - bits_down[k] ) );
1644 const int local_ak_up = bits_up[a] * ( 1 - bits_up[k] );
1645 const int local_ak_down = bits_down[a] * ( 1 - bits_down[k] );
1646
1647 if ( ( special_ak > 0 ) || ( local_ak_up > 0 ) || ( local_ak_down > 0 ) ){
1648
1649 const int irrep_ak = Irreps::directProd( getOrb2Irrep(a), getOrb2Irrep(k) );
1650
1651 for ( unsigned int i = 0; i < L; i++ ){
1652 if ( bits_up[i] + bits_down[i] < 2 ){
1653
1654 const int offset = Irreps::directProd( irrep_ak, getOrb2Irrep(i) ) * ( L + 1 );
1655 const int bar_i_up = 1 - bits_up[i];
1656 const int bar_i_down = 1 - bits_down[i];
1657 const int max_c_cnt = specific_orbs_irrep[ L + offset ];
1658
1659 for ( int c_cnt = 0; c_cnt < max_c_cnt; c_cnt++ ){
1660 const int c = specific_orbs_irrep[ c_cnt + offset ];
1661 const int fact_ic_up = bits_up[c] * bar_i_up;
1662 const int fact_ic_down = bits_down[c] * bar_i_down;
1663 const int prefactor1 = ( fact_ic_up + fact_ic_down ) * special_ak;
1664 const int prefactor2 = local_ak_up * fact_ic_up + local_ak_down * fact_ic_down;
1665 const double eri_akci = getERI(a, k, c, i);
1666 const double eri_aick = getERI(a, i, c, k);
1667 myResult += 0.5 * eri_akci * ( prefactor1 * eri_akci - prefactor2 * eri_aick );
1668 }
1669 }
1670 }
1671 }
1672 }
1673 }
1674 }
1675
1676 output[ counter ] = myResult;
1677
1678 }
1679
1680 delete [] bits_up;
1681 delete [] bits_down;
1682
1683 delete [] Jmat;
1684 delete [] K_reg_up;
1685 delete [] K_reg_down;
1686 delete [] K_bar_up;
1687 delete [] K_bar_down;
1688
1689 delete [] specific_orbs_irrep;
1690
1691 }
1692
1693 gettimeofday(&end, NULL);
1694 const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
1695 if ( FCIverbose > 0 ){ cout << "FCI::DiagHamSquared : Wall time = " << elapsed << " seconds" << endl; }
1696
1697 }
1698
1699
LowestEnergyDeterminant() const1700 unsigned int CheMPS2::FCI::LowestEnergyDeterminant() const{
1701
1702 const unsigned int vecLength = getVecLength( 0 );
1703 double * energies = new double[ vecLength ];
1704
1705 // Fetch the Slater determinant energies
1706 DiagHam( energies );
1707
1708 // Find the determinant with minimum energy
1709 unsigned int minEindex = 0;
1710 for ( unsigned int count = 1; count < vecLength; count++ ){
1711 if ( energies[ count ] < energies[ minEindex ] ){
1712 minEindex = count;
1713 }
1714 }
1715
1716 delete [] energies;
1717
1718 return minEindex;
1719
1720 }
1721
GetMatrixElement(int * bits_bra_up,int * bits_bra_down,int * bits_ket_up,int * bits_ket_down,int * work) const1722 double CheMPS2::FCI::GetMatrixElement(int * bits_bra_up, int * bits_bra_down, int * bits_ket_up, int * bits_ket_down, int * work) const{
1723
1724 int count_annih_up = 0;
1725 int count_creat_up = 0;
1726 int count_annih_down = 0;
1727 int count_creat_down = 0;
1728
1729 int * annih_up = work;
1730 int * creat_up = work+2;
1731 int * annih_down = work+4;
1732 int * creat_down = work+6;
1733
1734 // Find the differences between bra and ket, and store them in the arrays annih/creat
1735 for (unsigned int orb = 0; orb < L; orb++){
1736 if ( bits_bra_up[ orb ] != bits_ket_up[ orb ] ){
1737 if ( bits_ket_up[ orb ] ){ // Electron is annihilated in the ket
1738 if ( count_annih_up == 2 ){ return 0.0; }
1739 annih_up[ count_annih_up ] = orb;
1740 count_annih_up++;
1741 } else { // Electron is created in the ket
1742 if ( count_creat_up == 2 ){ return 0.0; }
1743 creat_up[ count_creat_up ] = orb;
1744 count_creat_up++;
1745 }
1746 }
1747 if ( bits_bra_down[ orb ] != bits_ket_down[ orb ] ){
1748 if ( bits_ket_down[ orb ] ){ // Electron is annihilated in the ket
1749 if ( count_annih_down == 2 ){ return 0.0; }
1750 annih_down[ count_annih_down ] = orb;
1751 count_annih_down++;
1752 } else { // Electron is created in the ket
1753 if ( count_creat_down == 2 ){ return 0.0; }
1754 creat_down[ count_creat_down ] = orb;
1755 count_creat_down++;
1756 }
1757 }
1758 }
1759
1760 // Sanity check: spin symmetry
1761 if ( count_annih_up != count_creat_up ){ return 0.0; }
1762 if ( count_annih_down != count_creat_down ){ return 0.0; }
1763
1764 // Sanity check: At most 2 annihilators and 2 creators can connect the ket and bra
1765 if ( count_annih_up + count_annih_down > 2 ){ return 0.0; }
1766 if ( count_creat_up + count_creat_down > 2 ){ return 0.0; }
1767
1768 if (( count_annih_up == 0 ) && ( count_annih_down == 0 )){ // |bra> == |ket> --> copy from DiagHam( )
1769
1770 double result = 0.0;
1771 for ( unsigned int orb1 = 0; orb1 < L; orb1++ ){
1772 const int n_tot_orb1 = bits_ket_up[ orb1 ] + bits_ket_down[ orb1 ];
1773 result += n_tot_orb1 * getGmat( orb1 , orb1 );
1774 for ( unsigned int orb2 = 0; orb2 < L; orb2++ ){
1775 result += 0.5 * n_tot_orb1 * ( bits_ket_up[ orb2 ] + bits_ket_down[ orb2 ] ) * getERI( orb1 , orb1 , orb2 , orb2 )
1776 + 0.5 * ( n_tot_orb1 - bits_ket_up[ orb1 ] * bits_ket_up[ orb2 ] - bits_ket_down[ orb1 ] * bits_ket_down[ orb2 ] ) * getERI( orb1 , orb2 , orb2 , orb1 );
1777 }
1778 }
1779 return result;
1780
1781 }
1782
1783 if (( count_annih_up == 1 ) && ( count_annih_down == 0 )){ // |bra> = a^+_j,up a_l,up |ket>
1784
1785 const int orbj = creat_up[ 0 ];
1786 const int orbl = annih_up[ 0 ];
1787
1788 double result = getGmat( orbj , orbl );
1789 for (unsigned int orb1 = 0; orb1 < L; orb1++){
1790 result += getERI(orbj, orb1, orb1, orbl) * ( 0.5 - bits_ket_up[ orb1 ] ) + getERI(orb1, orb1, orbj, orbl) * ( bits_ket_up[ orb1 ] + bits_ket_down[ orb1 ] );
1791 }
1792 int phase = 1;
1793 if ( orbj < orbl ){ for (int orbital = orbj+1; orbital < orbl; orbital++){ if ( bits_ket_up[ orbital ] ){ phase *= -1; } } }
1794 if ( orbl < orbj ){ for (int orbital = orbl+1; orbital < orbj; orbital++){ if ( bits_ket_up[ orbital ] ){ phase *= -1; } } }
1795 return ( result * phase );
1796
1797 }
1798
1799 if (( count_annih_up == 0 ) && ( count_annih_down == 1 )){ // |bra> = a^+_j,down a_l,down |ket>
1800
1801 const int orbj = creat_down[ 0 ];
1802 const int orbl = annih_down[ 0 ];
1803
1804 double result = getGmat( orbj , orbl );
1805 for (unsigned int orb1 = 0; orb1 < L; orb1++){
1806 result += getERI(orbj, orb1, orb1, orbl) * ( 0.5 - bits_ket_down[ orb1 ] ) + getERI(orb1, orb1, orbj, orbl) * ( bits_ket_up[ orb1 ] + bits_ket_down[ orb1 ] );
1807 }
1808 int phase = 1;
1809 if ( orbj < orbl ){ for (int orbital = orbj+1; orbital < orbl; orbital++){ if ( bits_ket_down[ orbital ] ){ phase *= -1; } } }
1810 if ( orbl < orbj ){ for (int orbital = orbl+1; orbital < orbj; orbital++){ if ( bits_ket_down[ orbital ] ){ phase *= -1; } } }
1811 return ( result * phase );
1812
1813 }
1814
1815 if (( count_annih_up == 2 ) && ( count_annih_down == 0 )){
1816
1817 // creat and annih are filled in increasing orbital index
1818 const int orbi = creat_up[ 0 ];
1819 const int orbj = creat_up[ 1 ];
1820 const int orbk = annih_up[ 0 ];
1821 const int orbl = annih_up[ 1 ];
1822
1823 double result = getERI(orbi, orbk, orbj, orbl) - getERI(orbi, orbl, orbj, orbk);
1824 int phase = 1;
1825 for (int orbital = orbk+1; orbital < orbl; orbital++){ if ( bits_ket_up[ orbital ] ){ phase *= -1; } } // Fermion phases orbk and orbl measured in the ket
1826 for (int orbital = orbi+1; orbital < orbj; orbital++){ if ( bits_bra_up[ orbital ] ){ phase *= -1; } } // Fermion phases orbi and orbj measured in the bra
1827 return ( result * phase );
1828
1829 }
1830
1831 if (( count_annih_up == 0 ) && ( count_annih_down == 2 )){
1832
1833 // creat and annih are filled in increasing orbital index
1834 const int orbi = creat_down[ 0 ];
1835 const int orbj = creat_down[ 1 ];
1836 const int orbk = annih_down[ 0 ];
1837 const int orbl = annih_down[ 1 ];
1838
1839 double result = getERI(orbi, orbk, orbj, orbl) - getERI(orbi, orbl, orbj, orbk);
1840 int phase = 1;
1841 for (int orbital = orbk+1; orbital < orbl; orbital++){ if ( bits_ket_down[ orbital ] ){ phase *= -1; } } // Fermion phases orbk and orbl measured in the ket
1842 for (int orbital = orbi+1; orbital < orbj; orbital++){ if ( bits_bra_down[ orbital ] ){ phase *= -1; } } // Fermion phases orbi and orbj measured in the bra
1843 return ( result * phase );
1844
1845 }
1846
1847 if (( count_annih_up == 1 ) && ( count_annih_down == 1 )){
1848
1849 const int orbi = creat_up [ 0 ];
1850 const int orbj = creat_down[ 0 ];
1851 const int orbk = annih_up [ 0 ];
1852 const int orbl = annih_down[ 0 ];
1853
1854 double result = getERI(orbi, orbk, orbj, orbl);
1855 int phase = 1;
1856 if ( orbi < orbk ){ for (int orbital = orbi+1; orbital < orbk; orbital++){ if ( bits_ket_up [ orbital ] ){ phase *= -1; } } }
1857 if ( orbk < orbi ){ for (int orbital = orbk+1; orbital < orbi; orbital++){ if ( bits_ket_up [ orbital ] ){ phase *= -1; } } }
1858 if ( orbj < orbl ){ for (int orbital = orbj+1; orbital < orbl; orbital++){ if ( bits_ket_down[ orbital ] ){ phase *= -1; } } }
1859 if ( orbl < orbj ){ for (int orbital = orbl+1; orbital < orbj; orbital++){ if ( bits_ket_down[ orbital ] ){ phase *= -1; } } }
1860 return ( result * phase );
1861
1862 }
1863
1864 return 0.0;
1865
1866 }
1867
FCIdcopy(const unsigned int vecLength,double * origin,double * target)1868 void CheMPS2::FCI::FCIdcopy(const unsigned int vecLength, double * origin, double * target){
1869
1870 int length = vecLength; // Checked "assert( max_integer >= maxVecLength );" at FCI::StartupIrrepCenter()
1871 int inc = 1;
1872 dcopy_( &length , origin , &inc , target , &inc );
1873
1874 }
1875
FCIddot(const unsigned int vecLength,double * vec1,double * vec2)1876 double CheMPS2::FCI::FCIddot(const unsigned int vecLength, double * vec1, double * vec2){
1877
1878 int length = vecLength; // Checked "assert( max_integer >= maxVecLength );" at FCI::StartupIrrepCenter()
1879 int inc = 1;
1880 return ddot_( &length , vec1 , &inc , vec2 , &inc );
1881
1882 }
1883
FCIfrobeniusnorm(const unsigned int vecLength,double * vec)1884 double CheMPS2::FCI::FCIfrobeniusnorm(const unsigned int vecLength, double * vec){
1885
1886 return sqrt( FCIddot( vecLength , vec , vec ) );
1887
1888 }
1889
FCIdaxpy(const unsigned int vecLength,const double alpha,double * vec_x,double * vec_y)1890 void CheMPS2::FCI::FCIdaxpy(const unsigned int vecLength, const double alpha, double * vec_x, double * vec_y){
1891
1892 double factor = alpha;
1893 int length = vecLength; // Checked "assert( max_integer >= maxVecLength );" at FCI::StartupIrrepCenter()
1894 int inc = 1;
1895 daxpy_( &length , &factor , vec_x , &inc , vec_y , &inc );
1896
1897 }
1898
FCIdscal(const unsigned int vecLength,const double alpha,double * vec)1899 void CheMPS2::FCI::FCIdscal(const unsigned int vecLength, const double alpha, double * vec){
1900
1901 double factor = alpha;
1902 int length = vecLength; // Checked "assert( max_integer >= maxVecLength );" at FCI::StartupIrrepCenter()
1903 int inc = 1;
1904 dscal_( &length , &factor , vec , &inc );
1905
1906 }
1907
ClearVector(const unsigned int vecLength,double * vec)1908 void CheMPS2::FCI::ClearVector(const unsigned int vecLength, double * vec){
1909
1910 for ( unsigned int cnt = 0; cnt < vecLength; cnt++ ){ vec[cnt] = 0.0; }
1911
1912 }
1913
FillRandom(const unsigned int vecLength,double * vec)1914 void CheMPS2::FCI::FillRandom(const unsigned int vecLength, double * vec){
1915
1916 for ( unsigned int cnt = 0; cnt < vecLength; cnt++ ){ vec[cnt] = ( ( 2.0 * rand() ) / RAND_MAX ) - 1.0; }
1917
1918 }
1919
GSDavidson(double * inoutput,const int DVDSN_NUM_VEC) const1920 double CheMPS2::FCI::GSDavidson(double * inoutput, const int DVDSN_NUM_VEC) const{
1921
1922 const int veclength = getVecLength( 0 ); // Checked "assert( max_integer >= maxVecLength );" at FCI::StartupIrrepCenter()
1923 Davidson deBoskabouter( veclength, DVDSN_NUM_VEC,
1924 CheMPS2::DAVIDSON_NUM_VEC_KEEP,
1925 CheMPS2::DAVIDSON_FCI_RTOL,
1926 CheMPS2::DAVIDSON_PRECOND_CUTOFF, false ); // No debug printing for FCI
1927 double ** whichpointers = new double*[2];
1928
1929 char instruction = deBoskabouter.FetchInstruction( whichpointers );
1930 assert( instruction == 'A' );
1931 if ( inoutput != NULL ){ FCIdcopy( veclength, inoutput, whichpointers[0] ); }
1932 else { FillRandom( veclength, whichpointers[0] ); }
1933 DiagHam( whichpointers[1] );
1934
1935 instruction = deBoskabouter.FetchInstruction( whichpointers );
1936 while ( instruction == 'B' ){
1937 matvec( whichpointers[0], whichpointers[1] );
1938 instruction = deBoskabouter.FetchInstruction( whichpointers );
1939 }
1940
1941 assert( instruction == 'C' );
1942 if ( inoutput != NULL ){ FCIdcopy( veclength, whichpointers[0], inoutput ); }
1943 const double FCIenergy = whichpointers[1][0] + getEconst();
1944 if ( FCIverbose > 1 ){ cout << "FCI::GSDavidson : Required number of matrix-vector multiplications = " << deBoskabouter.GetNumMultiplications() << endl; }
1945 if ( FCIverbose > 0 ){ cout << "FCI::GSDavidson : Converged ground state energy = " << FCIenergy << endl; }
1946 delete [] whichpointers;
1947 return FCIenergy;
1948
1949 }
1950
1951 /*********************************************************************************
1952 * *
1953 * Below this block all functions are for the Green's function calculations. *
1954 * *
1955 *********************************************************************************/
1956
ActWithNumberOperator(const unsigned int orbIndex,double * resultVector,double * sourceVector) const1957 void CheMPS2::FCI::ActWithNumberOperator(const unsigned int orbIndex, double * resultVector, double * sourceVector) const{
1958
1959 assert( orbIndex<L );
1960
1961 int * bits_up = new int[ L ];
1962 int * bits_down = new int[ L ];
1963
1964 const unsigned int vecLength = getVecLength( 0 );
1965 for (unsigned int counter = 0; counter < vecLength; counter++){
1966 getBitsOfCounter( 0 , counter , bits_up , bits_down );
1967 resultVector[ counter ] = ( bits_up[ orbIndex ] + bits_down[ orbIndex ] ) * sourceVector[ counter ];
1968 }
1969
1970 delete [] bits_up;
1971 delete [] bits_down;
1972
1973 }
1974
ActWithSecondQuantizedOperator(const char whichOperator,const bool isUp,const unsigned int orbIndex,double * thisVector,const FCI * otherFCI,double * otherVector) const1975 void CheMPS2::FCI::ActWithSecondQuantizedOperator(const char whichOperator, const bool isUp, const unsigned int orbIndex, double * thisVector, const FCI * otherFCI, double * otherVector) const{
1976
1977 assert( ( whichOperator=='C' ) || ( whichOperator=='A' ) ); //Operator should be a (C) Creator, or (A) Annihilator
1978 assert( orbIndex<L );
1979 assert( L==otherFCI->getL() );
1980
1981 const unsigned int vecLength = getVecLength( 0 );
1982
1983 if ( getTargetIrrep() != Irreps::directProd( otherFCI->getTargetIrrep() , getOrb2Irrep( orbIndex ) )){
1984 ClearVector( vecLength , thisVector );
1985 return;
1986 }
1987
1988 int * bits_up = new int[ L ];
1989 int * bits_down = new int[ L ];
1990
1991 if (( whichOperator=='C') && ( isUp )){
1992 for (unsigned int counter = 0; counter < vecLength; counter++){
1993
1994 getBitsOfCounter( 0 , counter , bits_up , bits_down );
1995
1996 if ( bits_up[ orbIndex ] == 1 ){ // Operator = creator_up
1997 bits_up[ orbIndex ] = 0;
1998 int phase = 1;
1999 for (unsigned int cnt = 0; cnt < orbIndex; cnt++){ if ( bits_up[ cnt ] ){ phase *= -1; } }
2000 thisVector[ counter ] = phase * otherFCI->getFCIcoeff( bits_up , bits_down , otherVector );
2001 } else {
2002 thisVector[ counter ] = 0.0;
2003 }
2004
2005 }
2006 }
2007
2008 if (( whichOperator=='C') && ( !(isUp) )){
2009 const int startphase = (( Nel_up % 2 ) == 0) ? 1 : -1;
2010 for (unsigned int counter = 0; counter < vecLength; counter++){
2011
2012 getBitsOfCounter( 0 , counter , bits_up , bits_down );
2013
2014 if ( bits_down[ orbIndex ] == 1 ){ // Operator = creator_down
2015 bits_down[ orbIndex ] = 0;
2016 int phase = startphase;
2017 for (unsigned int cnt = 0; cnt < orbIndex; cnt++){ if ( bits_down[ cnt ] ){ phase *= -1; } }
2018 thisVector[ counter ] = phase * otherFCI->getFCIcoeff( bits_up , bits_down , otherVector );
2019 } else {
2020 thisVector[ counter ] = 0.0;
2021 }
2022
2023 }
2024 }
2025
2026 if (( whichOperator=='A') && ( isUp )){
2027 for (unsigned int counter = 0; counter < vecLength; counter++){
2028
2029 getBitsOfCounter( 0 , counter , bits_up , bits_down );
2030
2031 if ( bits_up[ orbIndex ] == 0 ){ // Operator = annihilator_up
2032 bits_up[ orbIndex ] = 1;
2033 int phase = 1;
2034 for (unsigned int cnt = 0; cnt < orbIndex; cnt++){ if ( bits_up[ cnt ] ){ phase *= -1; } }
2035 thisVector[ counter ] = phase * otherFCI->getFCIcoeff( bits_up , bits_down , otherVector );
2036 } else {
2037 thisVector[ counter ] = 0.0;
2038 }
2039
2040 }
2041 }
2042
2043 if (( whichOperator=='A') && ( !(isUp) )){
2044 const int startphase = (( Nel_up % 2 ) == 0) ? 1 : -1;
2045 for (unsigned int counter = 0; counter < vecLength; counter++){
2046
2047 getBitsOfCounter( 0 , counter , bits_up , bits_down );
2048
2049 if ( bits_down[ orbIndex ] == 0 ){ // Operator = annihilator_down
2050 bits_down[ orbIndex ] = 1;
2051 int phase = startphase;
2052 for (unsigned int cnt = 0; cnt < orbIndex; cnt++){ if ( bits_down[ cnt ] ){ phase *= -1; } }
2053 thisVector[ counter ] = phase * otherFCI->getFCIcoeff( bits_up , bits_down , otherVector );
2054 } else {
2055 thisVector[ counter ] = 0.0;
2056 }
2057
2058 }
2059 }
2060
2061 delete [] bits_up;
2062 delete [] bits_down;
2063
2064 }
2065
CGSolveSystem(const double alpha,const double beta,const double eta,double * RHS,double * RealSol,double * ImagSol,const bool checkError) const2066 void CheMPS2::FCI::CGSolveSystem(const double alpha, const double beta, const double eta, double * RHS, double * RealSol, double * ImagSol, const bool checkError) const{
2067
2068 const unsigned int vecLength = getVecLength( 0 );
2069
2070 // Calculate the diagonal of the CG operator
2071 double * temp = new double[ vecLength ];
2072 double * diag = new double[ vecLength ];
2073 CGdiagonal( alpha, beta, eta, diag, temp );
2074
2075 assert( RealSol != NULL );
2076 assert( ImagSol != NULL );
2077 assert( fabs( eta ) > 0.0 );
2078
2079 /*
2080 ( alpha + beta H + I eta ) Solution = RHS
2081
2082 is solved with the conjugate gradient (CG) method. Solution = RealSol + I * ImagSol.
2083 CG requires a symmetric positive definite operator. Therefore:
2084
2085 [ ( alpha + beta H )^2 + eta^2 ] * Solution = ( alpha + beta H - I eta ) * RHS
2086
2087 Clue: Solve for ImagSol first. RealSol is then simply
2088
2089 RealSol = - ( alpha + beta H ) / eta * ImagSol
2090 */
2091
2092 /**** Solve for ImagSol ****/
2093 double ** pointers = new double*[ 3 ];
2094 double RMSerror = 0.0;
2095 {
2096 ConjugateGradient CG( vecLength, CheMPS2::CONJ_GRADIENT_RTOL, CheMPS2::CONJ_GRADIENT_PRECOND_CUTOFF, false );
2097 char instruction = CG.step( pointers );
2098 assert( instruction == 'A' );
2099 for (unsigned int cnt = 0; cnt < vecLength; cnt++){ pointers[ 0 ][ cnt ] = - eta * RHS[ cnt ] / diag[ cnt ]; } // Initial guess
2100 for (unsigned int cnt = 0; cnt < vecLength; cnt++){ pointers[ 1 ][ cnt ] = diag[ cnt ]; } // Diagonal of the operator
2101 for (unsigned int cnt = 0; cnt < vecLength; cnt++){ pointers[ 2 ][ cnt ] = - eta * RHS[ cnt ]; } // RHS of the problem
2102 instruction = CG.step( pointers );
2103 assert( instruction == 'B' );
2104 while ( instruction == 'B' ){
2105 CGoperator( alpha, beta, eta, pointers[ 0 ], temp, pointers[ 1 ] );
2106 instruction = CG.step( pointers );
2107 }
2108 assert( instruction == 'C' );
2109 FCIdcopy( vecLength, pointers[ 0 ], ImagSol );
2110 RMSerror += pointers[ 1 ][ 0 ] * pointers[ 1 ][ 0 ];
2111 }
2112
2113 /**** Solve for RealSol ****/
2114 {
2115 ConjugateGradient CG( vecLength, CheMPS2::CONJ_GRADIENT_RTOL, CheMPS2::CONJ_GRADIENT_PRECOND_CUTOFF, false );
2116 char instruction = CG.step( pointers );
2117 assert( instruction == 'A' );
2118 CGAlphaPlusBetaHAM( - alpha / eta, - beta / eta, ImagSol, pointers[ 0 ] ); // Initial guess real part can be obtained from the imaginary part
2119 for (unsigned int cnt = 0; cnt < vecLength; cnt++){ pointers[ 1 ][ cnt ] = diag[ cnt ]; } // Diagonal of the operator
2120 CGAlphaPlusBetaHAM( alpha, beta, RHS, pointers[ 2 ] ); // RHS of the problem
2121 instruction = CG.step( pointers );
2122 assert( instruction == 'B' );
2123 while ( instruction == 'B' ){
2124 CGoperator( alpha, beta, eta, pointers[ 0 ], temp, pointers[ 1 ] );
2125 instruction = CG.step( pointers );
2126 }
2127 assert( instruction == 'C' );
2128 FCIdcopy( vecLength, pointers[ 0 ], RealSol );
2129 RMSerror += pointers[ 1 ][ 0 ] * pointers[ 1 ][ 0 ];
2130 }
2131 RMSerror = sqrt( RMSerror );
2132 delete [] pointers;
2133 delete [] temp;
2134 delete [] diag;
2135
2136 if (( checkError ) && ( FCIverbose > 0 )){
2137 cout << "FCI::CGSolveSystem : RMS error when checking the solution = " << RMSerror << endl;
2138 }
2139
2140 }
2141
CGAlphaPlusBetaHAM(const double alpha,const double beta,double * in,double * out) const2142 void CheMPS2::FCI::CGAlphaPlusBetaHAM(const double alpha, const double beta, double * in, double * out) const{
2143
2144 matvec( in , out );
2145 const unsigned int vecLength = getVecLength( 0 );
2146 const double prefactor = alpha + beta * getEconst(); // matvec does only the parts with second quantized operators
2147 for (unsigned int cnt = 0; cnt < vecLength; cnt++){
2148 out[ cnt ] = prefactor * in[ cnt ] + beta * out[ cnt ]; // out = ( alpha + beta * H ) * in
2149 }
2150
2151 }
2152
CGoperator(const double alpha,const double beta,const double eta,double * in,double * temp,double * out) const2153 void CheMPS2::FCI::CGoperator(const double alpha, const double beta, const double eta, double * in, double * temp, double * out) const{
2154
2155 const unsigned int vecLength = getVecLength( 0 );
2156 CGAlphaPlusBetaHAM( alpha, beta, in, temp ); // temp = ( alpha + beta * H ) * in
2157 CGAlphaPlusBetaHAM( alpha, beta, temp, out ); // out = ( alpha + beta * H )^2 * in
2158 FCIdaxpy( vecLength, eta*eta, in, out ); // out = [ ( alpha + beta * H )^2 + eta*eta ] * in
2159
2160 }
2161
CGdiagonal(const double alpha,const double beta,const double eta,double * diagonal,double * workspace) const2162 void CheMPS2::FCI::CGdiagonal(const double alpha, const double beta, const double eta, double * diagonal, double * workspace) const{
2163
2164 // diagonal becomes diag[ ( alpha + beta * H )^2 + eta*eta ]
2165
2166 DiagHam( diagonal );
2167 DiagHamSquared( workspace );
2168
2169 const unsigned int vecLength = getVecLength( 0 );
2170 const double alpha_bis = alpha + beta * getEconst();
2171 const double factor1 = alpha_bis * alpha_bis + eta * eta;
2172 const double factor2 = 2 * alpha_bis * beta;
2173 const double factor3 = beta * beta;
2174 for (unsigned int row = 0; row < vecLength; row++){
2175 diagonal[ row ] = factor1 + factor2 * diagonal[ row ] + factor3 * workspace[ row ];
2176 }
2177
2178 if ( FCIverbose>1 ){
2179 double minval = diagonal[0];
2180 double maxval = diagonal[0];
2181 for (unsigned int cnt = 1; cnt < vecLength; cnt++){
2182 if ( diagonal[ cnt ] > maxval ){ maxval = diagonal[ cnt ]; }
2183 if ( diagonal[ cnt ] < minval ){ minval = diagonal[ cnt ]; }
2184 }
2185 cout << "FCI::CGdiagonal : Minimum value of diag[ ( alpha + beta * Ham )^2 + eta^2 ] = " << minval << endl;
2186 cout << "FCI::CGdiagonal : Maximum value of diag[ ( alpha + beta * Ham )^2 + eta^2 ] = " << maxval << endl;
2187 }
2188
2189 }
2190
RetardedGF(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const bool isUp,const double GSenergy,double * GSvector,CheMPS2::Hamiltonian * Ham,double * RePartGF,double * ImPartGF) const2191 void CheMPS2::FCI::RetardedGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF) const{
2192
2193 assert( RePartGF != NULL );
2194 assert( ImPartGF != NULL );
2195
2196 // G( omega, alpha, beta, eta ) = < 0 | a_{alpha,spin} [ omega - Ham + E_0 + I*eta ]^{-1} a^+_{beta,spin} | 0 > (addition amplitude)
2197 // + < 0 | a^+_{beta,spin} [ omega + Ham - E_0 + I*eta ]^{-1} a_{alpha,spin} | 0 > (removal amplitude)
2198
2199 double Realpart, Imagpart;
2200 RetardedGF_addition(omega, eta, orb_alpha, orb_beta, isUp, GSenergy, GSvector, Ham, &Realpart, &Imagpart);
2201 RePartGF[0] = Realpart; // Set
2202 ImPartGF[0] = Imagpart; // Set
2203
2204 RetardedGF_removal( omega, eta, orb_alpha, orb_beta, isUp, GSenergy, GSvector, Ham, &Realpart, &Imagpart);
2205 RePartGF[0] += Realpart; // Add
2206 ImPartGF[0] += Imagpart;
2207
2208 if ( FCIverbose>0 ){
2209 cout << "FCI::RetardedGF : G( omega = " << omega << " ; eta = " << eta << " ; i = " << orb_alpha << " ; j = " << orb_beta << " ) = " << RePartGF[0] << " + I * " << ImPartGF[0] << endl;
2210 cout << " Local density of states (LDOS) = " << - ImPartGF[0] / M_PI << endl;
2211 }
2212
2213 }
2214
GFmatrix_addition(const double alpha,const double beta,const double eta,int * orbsLeft,const unsigned int numLeft,int * orbsRight,const unsigned int numRight,const bool isUp,double * GSvector,CheMPS2::Hamiltonian * Ham,double * RePartsGF,double * ImPartsGF,double ** TwoRDMreal,double ** TwoRDMimag,double ** TwoRDMadd) const2215 void CheMPS2::FCI::GFmatrix_addition(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal, double ** TwoRDMimag, double ** TwoRDMadd) const{
2216
2217 /*
2218 1
2219 GF[i + numLeft * j] = < 0 | a_{orbsLeft[i], spin} -------------------------------- a^+_{orbsRight[j], spin} | 0 >
2220 [ alpha + beta * Ham + I*eta ]
2221 */
2222
2223 // Check whether some stuff is OK
2224 assert( numLeft > 0 );
2225 assert( numRight > 0 );
2226 for (unsigned int cnt = 0; cnt < numLeft; cnt++){ int orbl = orbsLeft[ cnt ]; assert((orbl < L) && (orbl >= 0)); }
2227 for (unsigned int cnt = 0; cnt < numRight; cnt++){ int orbr = orbsRight[ cnt ]; assert((orbr < L) && (orbr >= 0)); }
2228 assert( RePartsGF != NULL );
2229 assert( ImPartsGF != NULL );
2230 for ( unsigned int counter = 0; counter < numLeft * numRight; counter++ ){
2231 RePartsGF[ counter ] = 0.0;
2232 ImPartsGF[ counter ] = 0.0;
2233 }
2234 const unsigned int Lpow4 = L*L*L*L;
2235 for ( unsigned int cnt = 0; cnt < numRight; cnt++ ){
2236 if ( TwoRDMreal != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMreal[ cnt ][ elem ] = 0.0; } }
2237 if ( TwoRDMimag != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMimag[ cnt ][ elem ] = 0.0; } }
2238 if ( TwoRDMadd != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMadd[ cnt ][ elem ] = 0.0; } }
2239 }
2240
2241 const bool isOK = ( isUp ) ? ( getNel_up() < L ) : ( getNel_down() < L ); // The electron can be added
2242 for ( unsigned int cnt_right = 0; cnt_right < numRight; cnt_right++ ){
2243
2244 const int orbitalRight = orbsRight[ cnt_right ];
2245 bool matchingIrrep = false;
2246 for ( unsigned int cnt_left = 0; cnt_left < numLeft; cnt_left++ ){
2247 if ( getOrb2Irrep( orbsLeft[ cnt_left] ) == getOrb2Irrep( orbitalRight ) ){ matchingIrrep = true; }
2248 }
2249
2250 if ( isOK && matchingIrrep ){
2251
2252 const unsigned int addNelUP = getNel_up() + ((isUp) ? 1 : 0);
2253 const unsigned int addNelDOWN = getNel_down() + ((isUp) ? 0 : 1);
2254 const int addIrrep = Irreps::directProd( getTargetIrrep(), getOrb2Irrep( orbitalRight ) );
2255
2256 CheMPS2::FCI additionFCI( Ham, addNelUP, addNelDOWN, addIrrep, maxMemWorkMB, FCIverbose );
2257 const unsigned int addVecLength = additionFCI.getVecLength( 0 );
2258 double * addVector = new double[ addVecLength ];
2259 additionFCI.ActWithSecondQuantizedOperator( 'C', isUp, orbitalRight, addVector, this, GSvector ); // | addVector > = a^+_right,spin | GSvector >
2260
2261 double * RealPartSolution = new double[ addVecLength ];
2262 double * ImagPartSolution = new double[ addVecLength ];
2263 additionFCI.CGSolveSystem( alpha, beta, eta, addVector, RealPartSolution, ImagPartSolution );
2264
2265 if ( TwoRDMreal != NULL ){ additionFCI.Fill2RDM( RealPartSolution, TwoRDMreal[ cnt_right ] ); }
2266 if ( TwoRDMimag != NULL ){ additionFCI.Fill2RDM( ImagPartSolution, TwoRDMimag[ cnt_right ] ); }
2267 if ( TwoRDMadd != NULL ){ additionFCI.Fill2RDM( addVector, TwoRDMadd[ cnt_right ] ); }
2268
2269 for ( unsigned int cnt_left = 0; cnt_left < numLeft; cnt_left++ ){
2270 const int orbitalLeft = orbsLeft[ cnt_left ];
2271 if ( getOrb2Irrep( orbitalLeft ) == getOrb2Irrep( orbitalRight ) ){
2272 additionFCI.ActWithSecondQuantizedOperator( 'C', isUp, orbitalLeft, addVector, this, GSvector ); // | addVector > = a^+_left,spin | GSvector >
2273 RePartsGF[ cnt_left + numLeft * cnt_right ] = FCIddot( addVecLength, addVector, RealPartSolution );
2274 ImPartsGF[ cnt_left + numLeft * cnt_right ] = FCIddot( addVecLength, addVector, ImagPartSolution );
2275 }
2276 }
2277
2278 delete [] RealPartSolution;
2279 delete [] ImagPartSolution;
2280 delete [] addVector;
2281
2282 }
2283 }
2284
2285 }
2286
RetardedGF_addition(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const bool isUp,const double GSenergy,double * GSvector,CheMPS2::Hamiltonian * Ham,double * RePartGF,double * ImPartGF,double * TwoRDMreal,double * TwoRDMimag,double * TwoRDMadd) const2287 void CheMPS2::FCI::RetardedGF_addition(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal, double * TwoRDMimag, double * TwoRDMadd) const{
2288
2289 // Addition amplitude < 0 | a_{alpha, spin} [ omega - Ham + E_0 + I*eta ]^{-1} a^+_{beta, spin} | 0 >
2290
2291 double ** TwoRDMreal_wrap = NULL; if ( TwoRDMreal != NULL ){ TwoRDMreal_wrap = new double*[1]; TwoRDMreal_wrap[0] = TwoRDMreal; }
2292 double ** TwoRDMimag_wrap = NULL; if ( TwoRDMimag != NULL ){ TwoRDMimag_wrap = new double*[1]; TwoRDMimag_wrap[0] = TwoRDMimag; }
2293 double ** TwoRDMadd_wrap = NULL; if ( TwoRDMadd != NULL ){ TwoRDMadd_wrap = new double*[1]; TwoRDMadd_wrap[0] = TwoRDMadd; }
2294
2295 int orb_left = orb_alpha;
2296 int orb_right = orb_beta;
2297
2298 GFmatrix_addition( omega + GSenergy, -1.0, eta, &orb_left, 1, &orb_right, 1, isUp, GSvector, Ham, RePartGF, ImPartGF, TwoRDMreal_wrap, TwoRDMimag_wrap, TwoRDMadd_wrap );
2299
2300 if ( TwoRDMreal != NULL ){ delete [] TwoRDMreal_wrap; }
2301 if ( TwoRDMimag != NULL ){ delete [] TwoRDMimag_wrap; }
2302 if ( TwoRDMadd != NULL ){ delete [] TwoRDMadd_wrap; }
2303
2304 }
2305
GFmatrix_removal(const double alpha,const double beta,const double eta,int * orbsLeft,const unsigned int numLeft,int * orbsRight,const unsigned int numRight,const bool isUp,double * GSvector,CheMPS2::Hamiltonian * Ham,double * RePartsGF,double * ImPartsGF,double ** TwoRDMreal,double ** TwoRDMimag,double ** TwoRDMrem) const2306 void CheMPS2::FCI::GFmatrix_removal(const double alpha, const double beta, const double eta, int * orbsLeft, const unsigned int numLeft, int * orbsRight, const unsigned int numRight, const bool isUp, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartsGF, double * ImPartsGF, double ** TwoRDMreal, double ** TwoRDMimag, double ** TwoRDMrem) const{
2307
2308 /*
2309 1
2310 GF[i + numLeft * j] = < 0 | a^+_{orbsLeft[i], spin} -------------------------------- a_{orbsRight[j], spin} | 0 >
2311 [ alpha + beta * Ham + I*eta ]
2312 */
2313
2314 // Check whether some stuff is OK
2315 assert( numLeft > 0 );
2316 assert( numRight > 0 );
2317 for (unsigned int cnt = 0; cnt < numLeft; cnt++){ int orbl = orbsLeft [ cnt ]; assert((orbl < L) && (orbl >= 0)); }
2318 for (unsigned int cnt = 0; cnt < numRight; cnt++){ int orbr = orbsRight[ cnt ]; assert((orbr < L) && (orbr >= 0)); }
2319 assert( RePartsGF != NULL );
2320 assert( ImPartsGF != NULL );
2321 for ( unsigned int counter = 0; counter < numLeft * numRight; counter++ ){
2322 RePartsGF[ counter ] = 0.0;
2323 ImPartsGF[ counter ] = 0.0;
2324 }
2325 const unsigned int Lpow4 = L*L*L*L;
2326 for ( unsigned int cnt = 0; cnt < numRight; cnt++ ){
2327 if ( TwoRDMreal != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMreal[ cnt ][ elem ] = 0.0; } }
2328 if ( TwoRDMimag != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMimag[ cnt ][ elem ] = 0.0; } }
2329 if ( TwoRDMrem != NULL ){ for ( unsigned int elem = 0; elem < Lpow4; elem++ ){ TwoRDMrem[ cnt ][ elem ] = 0.0; } }
2330 }
2331
2332 const bool isOK = ( isUp ) ? ( getNel_up() > 0 ) : ( getNel_down() > 0 ); // The electron can be removed
2333 for ( unsigned int cnt_right = 0; cnt_right < numRight; cnt_right++ ){
2334
2335 const int orbitalRight = orbsRight[ cnt_right ];
2336 bool matchingIrrep = false;
2337 for ( unsigned int cnt_left = 0; cnt_left < numLeft; cnt_left++ ){
2338 if ( getOrb2Irrep( orbsLeft[ cnt_left] ) == getOrb2Irrep( orbitalRight ) ){ matchingIrrep = true; }
2339 }
2340
2341 if ( isOK && matchingIrrep ){
2342
2343 const unsigned int removeNelUP = getNel_up() - ((isUp) ? 1 : 0);
2344 const unsigned int removeNelDOWN = getNel_down() - ((isUp) ? 0 : 1);
2345 const int removeIrrep = Irreps::directProd( getTargetIrrep(), getOrb2Irrep( orbitalRight ) );
2346
2347 CheMPS2::FCI removalFCI( Ham, removeNelUP, removeNelDOWN, removeIrrep, maxMemWorkMB, FCIverbose );
2348 const unsigned int removeVecLength = removalFCI.getVecLength( 0 );
2349 double * removeVector = new double[ removeVecLength ];
2350 removalFCI.ActWithSecondQuantizedOperator( 'A', isUp, orbitalRight, removeVector, this, GSvector ); // | removeVector > = a_right,spin | GSvector >
2351
2352 double * RealPartSolution = new double[ removeVecLength ];
2353 double * ImagPartSolution = new double[ removeVecLength ];
2354 removalFCI.CGSolveSystem( alpha, beta, eta, removeVector, RealPartSolution, ImagPartSolution );
2355
2356 if ( TwoRDMreal != NULL ){ removalFCI.Fill2RDM( RealPartSolution, TwoRDMreal[ cnt_right ] ); }
2357 if ( TwoRDMimag != NULL ){ removalFCI.Fill2RDM( ImagPartSolution, TwoRDMimag[ cnt_right ] ); }
2358 if ( TwoRDMrem != NULL ){ removalFCI.Fill2RDM( removeVector, TwoRDMrem[ cnt_right ] ); }
2359
2360 for ( unsigned int cnt_left = 0; cnt_left < numLeft; cnt_left++ ){
2361 const int orbitalLeft = orbsLeft[ cnt_left ];
2362 if ( getOrb2Irrep( orbitalLeft ) == getOrb2Irrep( orbitalRight ) ){
2363 removalFCI.ActWithSecondQuantizedOperator( 'A', isUp, orbitalLeft, removeVector, this, GSvector ); // | removeVector > = a_left,spin | GSvector >
2364 RePartsGF[ cnt_left + numLeft * cnt_right ] = FCIddot( removeVecLength, removeVector, RealPartSolution );
2365 ImPartsGF[ cnt_left + numLeft * cnt_right ] = FCIddot( removeVecLength, removeVector, ImagPartSolution );
2366 }
2367 }
2368
2369 delete [] RealPartSolution;
2370 delete [] ImagPartSolution;
2371 delete [] removeVector;
2372
2373 }
2374 }
2375
2376 }
2377
RetardedGF_removal(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const bool isUp,const double GSenergy,double * GSvector,CheMPS2::Hamiltonian * Ham,double * RePartGF,double * ImPartGF,double * TwoRDMreal,double * TwoRDMimag,double * TwoRDMrem) const2378 void CheMPS2::FCI::RetardedGF_removal(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const bool isUp, const double GSenergy, double * GSvector, CheMPS2::Hamiltonian * Ham, double * RePartGF, double * ImPartGF, double * TwoRDMreal, double * TwoRDMimag, double * TwoRDMrem) const{
2379
2380 // Removal amplitude < 0 | a^+_{beta, spin} [ omega + Ham - E_0 + I*eta ]^{-1} a_{alpha, spin} | 0 >
2381
2382 double ** TwoRDMreal_wrap = NULL; if ( TwoRDMreal != NULL ){ TwoRDMreal_wrap = new double*[1]; TwoRDMreal_wrap[0] = TwoRDMreal; }
2383 double ** TwoRDMimag_wrap = NULL; if ( TwoRDMimag != NULL ){ TwoRDMimag_wrap = new double*[1]; TwoRDMimag_wrap[0] = TwoRDMimag; }
2384 double ** TwoRDMrem_wrap = NULL; if ( TwoRDMrem != NULL ){ TwoRDMrem_wrap = new double*[1]; TwoRDMrem_wrap[0] = TwoRDMrem; }
2385
2386 int orb_left = orb_beta;
2387 int orb_right = orb_alpha;
2388
2389 // orb_alpha = orb_right in this case !!
2390 GFmatrix_removal( omega - GSenergy, 1.0, eta, &orb_left, 1, &orb_right, 1, isUp, GSvector, Ham, RePartGF, ImPartGF, TwoRDMreal_wrap, TwoRDMimag_wrap, TwoRDMrem_wrap );
2391
2392 if ( TwoRDMreal != NULL ){ delete [] TwoRDMreal_wrap; }
2393 if ( TwoRDMimag != NULL ){ delete [] TwoRDMimag_wrap; }
2394 if ( TwoRDMrem != NULL ){ delete [] TwoRDMrem_wrap; }
2395
2396 }
2397
DensityResponseGF(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const double GSenergy,double * GSvector,double * RePartGF,double * ImPartGF) const2398 void CheMPS2::FCI::DensityResponseGF(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF) const{
2399
2400 assert( RePartGF != NULL );
2401 assert( ImPartGF != NULL );
2402
2403 // X( omega, alpha, beta, eta ) = < 0 | ( n_alpha - <0| n_alpha |0> ) [ omega - Ham + E_0 + I*eta ]^{-1} ( n_beta - <0| n_beta |0> ) | 0 > (forward amplitude)
2404 // - < 0 | ( n_beta - <0| n_beta |0> ) [ omega + Ham - E_0 + I*eta ]^{-1} ( n_alpha - <0| n_alpha |0> ) | 0 > (backward amplitude)
2405
2406 double Realpart, Imagpart;
2407 DensityResponseGF_forward( omega, eta, orb_alpha, orb_beta, GSenergy, GSvector, &Realpart, &Imagpart);
2408 RePartGF[0] = Realpart; // Set
2409 ImPartGF[0] = Imagpart; // Set
2410
2411 DensityResponseGF_backward(omega, eta, orb_alpha, orb_beta, GSenergy, GSvector, &Realpart, &Imagpart);
2412 RePartGF[0] -= Realpart; // Subtract !!!
2413 ImPartGF[0] -= Imagpart; // Subtract !!!
2414
2415 if ( FCIverbose>0 ){
2416 cout << "FCI::DensityResponseGF : X( omega = " << omega << " ; eta = " << eta << " ; i = " << orb_alpha << " ; j = " << orb_beta << " ) = " << RePartGF[0] << " + I * " << ImPartGF[0] << endl;
2417 cout << " Local density-density response (LDDR) = " << - ImPartGF[0] / M_PI << endl;
2418 }
2419
2420 }
2421
DensityResponseGF_forward(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const double GSenergy,double * GSvector,double * RePartGF,double * ImPartGF,double * TwoRDMreal,double * TwoRDMimag,double * TwoRDMdens) const2422 void CheMPS2::FCI::DensityResponseGF_forward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal, double * TwoRDMimag, double * TwoRDMdens) const{
2423
2424 // Forward amplitude: < 0 | ( n_alpha - <0| n_alpha |0> ) [ omega - Ham + E_0 + I*eta ]^{-1} ( n_beta - <0| n_beta |0> ) | 0 >
2425
2426 assert( ( orb_alpha<L ) && ( orb_beta<L ) ); // Orbital indices within bound
2427 assert( RePartGF != NULL );
2428 assert( ImPartGF != NULL );
2429
2430 const unsigned int vecLength = getVecLength( 0 );
2431 double * densityAlphaVector = new double[ vecLength ];
2432 double * densityBetaVector = ( orb_alpha == orb_beta ) ? densityAlphaVector : new double[ vecLength ];
2433 ActWithNumberOperator( orb_alpha , densityAlphaVector , GSvector ); // densityAlphaVector = n_alpha |0>
2434 const double n_alpha_0 = FCIddot( vecLength , densityAlphaVector , GSvector ); // <0| n_alpha |0>
2435 FCIdaxpy( vecLength , -n_alpha_0 , GSvector , densityAlphaVector ); // densityAlphaVector = ( n_alpha - <0| n_alpha |0> ) |0>
2436 if ( orb_alpha != orb_beta ){
2437 ActWithNumberOperator( orb_beta , densityBetaVector , GSvector ); // densityBetaVector = n_beta |0>
2438 const double n_beta_0 = FCIddot( vecLength , densityBetaVector , GSvector ); // <0| n_beta |0>
2439 FCIdaxpy( vecLength , -n_beta_0 , GSvector , densityBetaVector ); // densityBetaVector = ( n_beta - <0| n_beta |0> ) |0>
2440 }
2441
2442 double * RealPartSolution = new double[ vecLength ];
2443 double * ImagPartSolution = new double[ vecLength ];
2444 CGSolveSystem( omega + GSenergy , -1.0 , eta , densityBetaVector , RealPartSolution , ImagPartSolution );
2445 if ( TwoRDMreal != NULL ){ Fill2RDM( RealPartSolution , TwoRDMreal ); } // Sets the TwoRDMreal
2446 RePartGF[0] = FCIddot( vecLength , densityAlphaVector , RealPartSolution );
2447 delete [] RealPartSolution;
2448 if ( TwoRDMimag != NULL ){ Fill2RDM( ImagPartSolution , TwoRDMimag ); } // Sets the TwoRDMimag
2449 ImPartGF[0] = FCIddot( vecLength , densityAlphaVector , ImagPartSolution );
2450 delete [] ImagPartSolution;
2451
2452 if ( TwoRDMdens != NULL ){ Fill2RDM( densityBetaVector , TwoRDMdens ); } // Sets the TwoRDMdens
2453 if ( orb_alpha != orb_beta ){ delete [] densityBetaVector; }
2454 delete [] densityAlphaVector;
2455
2456 }
2457
DensityResponseGF_backward(const double omega,const double eta,const unsigned int orb_alpha,const unsigned int orb_beta,const double GSenergy,double * GSvector,double * RePartGF,double * ImPartGF,double * TwoRDMreal,double * TwoRDMimag,double * TwoRDMdens) const2458 void CheMPS2::FCI::DensityResponseGF_backward(const double omega, const double eta, const unsigned int orb_alpha, const unsigned int orb_beta, const double GSenergy, double * GSvector, double * RePartGF, double * ImPartGF, double * TwoRDMreal, double * TwoRDMimag, double * TwoRDMdens) const{
2459
2460 // Backward amplitude: < 0 | ( n_beta - <0| n_beta |0> ) [ omega + Ham - E_0 + I*eta ]^{-1} ( n_alpha - <0| n_alpha |0> ) | 0 >
2461
2462 assert( ( orb_alpha<L ) && ( orb_beta<L ) ); // Orbital indices within bound
2463 assert( RePartGF != NULL );
2464 assert( ImPartGF != NULL );
2465
2466 const unsigned int vecLength = getVecLength( 0 );
2467 double * densityAlphaVector = new double[ vecLength ];
2468 double * densityBetaVector = ( orb_alpha == orb_beta ) ? densityAlphaVector : new double[ vecLength ];
2469 ActWithNumberOperator( orb_alpha , densityAlphaVector , GSvector ); // densityAlphaVector = n_alpha |0>
2470 const double n_alpha_0 = FCIddot( vecLength , densityAlphaVector , GSvector ); // <0| n_alpha |0>
2471 FCIdaxpy( vecLength , -n_alpha_0 , GSvector , densityAlphaVector ); // densityAlphaVector = ( n_alpha - <0| n_alpha |0> ) |0>
2472 if ( orb_alpha != orb_beta ){
2473 ActWithNumberOperator( orb_beta , densityBetaVector , GSvector ); // densityBetaVector = n_beta |0>
2474 const double n_beta_0 = FCIddot( vecLength , densityBetaVector , GSvector ); // <0| n_beta |0>
2475 FCIdaxpy( vecLength , -n_beta_0 , GSvector , densityBetaVector ); // densityBetaVector = ( n_beta - <0| n_beta |0> ) |0>
2476 }
2477
2478 double * RealPartSolution = new double[ vecLength ];
2479 double * ImagPartSolution = new double[ vecLength ];
2480 CGSolveSystem( omega - GSenergy , 1.0 , eta , densityAlphaVector , RealPartSolution , ImagPartSolution );
2481 if ( TwoRDMreal != NULL ){ Fill2RDM( RealPartSolution , TwoRDMreal ); } // Sets the TwoRDMreal
2482 RePartGF[0] = FCIddot( vecLength , densityBetaVector , RealPartSolution );
2483 delete [] RealPartSolution;
2484 if ( TwoRDMimag != NULL ){ Fill2RDM( ImagPartSolution , TwoRDMimag ); } // Sets the TwoRDMimag
2485 ImPartGF[0] = FCIddot( vecLength , densityBetaVector , ImagPartSolution );
2486 delete [] ImagPartSolution;
2487
2488 if ( TwoRDMdens != NULL ){ Fill2RDM( densityAlphaVector , TwoRDMdens ); } // Sets the TwoRDMdens
2489 if ( orb_alpha != orb_beta ){ delete [] densityBetaVector; }
2490 delete [] densityAlphaVector;
2491
2492 }
2493
2494
2495