CheMPS2
FCI.cpp
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013-2016 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 
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 
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 
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 
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 
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 
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 
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 
409 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 
417 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 
435 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 
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 
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 
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 
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 
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 
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 
618 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_( &notrans, &notrans, &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_( &notrans, &notrans, &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 
772 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 
805 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 
896 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 
1160 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 
1168 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 
1176 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 
1184 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 
1405 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 
1467 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 
1504 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 
1700 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 
1722 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 
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 
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 
1884 double CheMPS2::FCI::FCIfrobeniusnorm(const unsigned int vecLength, double * vec){
1885 
1886  return sqrt( FCIddot( vecLength , vec , vec ) );
1887 
1888 }
1889 
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 
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 
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 
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 
1920 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 
1957 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 
1975 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 
2066 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 
2142 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 
2153 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 
2162 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 
2191 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 
2215 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 
2287 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 
2306 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 
2378 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 
2398 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 
2422 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 
2458 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 
CheMPS2::FCI::FCIddot
static double FCIddot(const unsigned int vecLength, double *vec1, double *vec2)
Take the inproduct of two vectors.
Definition: FCI.cpp:1876
CheMPS2::FCI::Diag4RDM
void Diag4RDM(double *vector, double *three_rdm, const unsigned int orbz, double *output) const
Construct part of the 4-RDM: output(i,j,k,p,q,r) = Gamma^4(i,j,k,z,p,q,r,z)
Definition: FCI.cpp:1176
CheMPS2::FCI::apply_excitation
void apply_excitation(double *orig_vector, double *result_vector, const int crea, const int anni, const int orig_target_irrep) const
Apply E_{crea,anni} to |orig_vector> and store the result in |result_vector>
Definition: FCI.cpp:772
CheMPS2::FCI::getEconst
double getEconst() const
Function which returns the nuclear repulsion energy.
Definition: FCI.h:92
CheMPS2::FCI::DensityResponseGF
void 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
Calculate the density response Green&#39;s function (= forward - backward propagating part) ...
Definition: FCI.cpp:2398
CheMPS2::FCI::DensityResponseGF_backward
void 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=NULL, double *TwoRDMimag=NULL, double *TwoRDMdens=NULL) const
Calculate the backward propagating part of the density response Green&#39;s function: <GSvector| ( n_beta...
Definition: FCI.cpp:2458
CheMPS2::FCI::getERI
double getERI(const int orb1, const int orb2, const int orb3, const int orb4) const
Function which returns the electron repulsion integral (in chemists notation: integral dr1 dr2 orb1(r...
Definition: FCI.h:82
CheMPS2::FCI::DiagHam
void DiagHam(double *diag) const
Function which returns the diagonal elements of the FCI Hamiltonian (without Econstant!!) = Slater de...
Definition: FCI.cpp:1467
CheMPS2::FCI::FCIdcopy
static void FCIdcopy(const unsigned int vecLength, double *origin, double *target)
Copy a vector.
Definition: FCI.cpp:1868
CheMPS2::FCI::CGdiagonal
void CGdiagonal(const double alpha, const double beta, const double eta, double *diagonal, double *workspace) const
Calculate (without approximation) diagonal = diag[ (alpha + beta * Hamiltonian)^2 + eta^2 ] (Econstan...
Definition: FCI.cpp:2162
CheMPS2::Davidson::GetNumMultiplications
int GetNumMultiplications() const
Get the number of matrix vector multiplications which have been performed.
Definition: Davidson.cpp:103
CheMPS2::FCI::GSDavidson
double GSDavidson(double *inoutput=NULL, const int DVDSN_NUM_VEC=CheMPS2::DAVIDSON_NUM_VEC) const
Calculates the FCI ground state with Davidson&#39;s algorithm.
Definition: FCI.cpp:1920
CheMPS2::Irreps::getNumberOfIrreps
int getNumberOfIrreps() const
Get the number of irreps for the currently activated group.
Definition: Irreps.cpp:94
CheMPS2::ConjugateGradient::step
char step(double **pointers)
The iterator to converge the ground state vector.
Definition: ConjugateGradient.cpp:63
CheMPS2::FCI::ActWithSecondQuantizedOperator
void ActWithSecondQuantizedOperator(const char whichOperator, const bool isUp, const unsigned int orbIndex, double *thisVector, const FCI *otherFCI, double *otherVector) const
Set thisVector to a creator/annihilator acting on otherVector.
Definition: FCI.cpp:1975
CheMPS2::FCI::Fill4RDM
void Fill4RDM(double *vector, double *FourRDM) const
Construct the (spin-summed) 4-RDM of a FCI vector: Gamma^4(i,j,k,l,p,q,r,t) = sum_sigma,tau,s,z < a^+_{i,sigma} a^+_{j,tau} a^+_{k,s} a^+_{l,z} a_{t,z} a_{r,s} a_{q,tau} a_{p,sigma} > = FourRDM[ i + L * ( j + L * ( k + L * ( l + L * ( p + L * ( q + L * ( r + L * t ) ) ) ) ) ) ].
Definition: FCI.cpp:896
CheMPS2::Hamiltonian::getOrbitalIrrep
int getOrbitalIrrep(const int nOrb) const
Get an orbital irrep number.
Definition: Hamiltonian.cpp:97
CheMPS2::FCI::FCIdaxpy
static void FCIdaxpy(const unsigned int vecLength, const double alpha, double *vec_x, double *vec_y)
Do lapack&#39;s daxpy vec_y += alpha * vec_x.
Definition: FCI.cpp:1890
CheMPS2::FCI::DensityResponseGF_forward
void 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=NULL, double *TwoRDMimag=NULL, double *TwoRDMdens=NULL) const
Calculate the forward propagating part of the density response Green&#39;s function: <GSvector| ( n_alpha...
Definition: FCI.cpp:2422
CheMPS2::Hamiltonian::getVmat
double getVmat(const int index1, const int index2, const int index3, const int index4) const
Get a Vmat element.
Definition: Hamiltonian.cpp:136
CheMPS2::FCI::GFmatrix_removal
void 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=NULL, double **TwoRDMimag=NULL, double **TwoRDMrem=NULL) const
Calculate the removal Green&#39;s function: GF[i+numLeft*j] = <GSvector| a^+_{orbsLeft[i], spin} [ alpha + beta * Ham + I*eta ]^{-1} a_{orbsRight[j], spin} |GSvector>
Definition: FCI.cpp:2306
CheMPS2::FCI::DiagHamSquared
void DiagHamSquared(double *output) const
Function which returns the diagonal elements of the FCI Hamiltonian squared (without Econstant!!) ...
Definition: FCI.cpp:1504
CheMPS2::FCI::FillRandom
static void FillRandom(const unsigned int vecLength, double *vec)
Fill a vector with random numbers in the interval [-1,1[; used when output for GSDavidson is desired ...
Definition: FCI.cpp:1914
CheMPS2::FCI::RetardedGF_removal
void 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=NULL, double *TwoRDMimag=NULL, double *TwoRDMrem=NULL) const
Calculate the removal part of the retarded Green&#39;s function: <GSvector| a^+_{beta, spin(isUp)} [ omega + Ham - GSenergy + I*eta ]^{-1} a_{alpha, spin(isUp)} |GSvector>
Definition: FCI.cpp:2378
CheMPS2::FCI::FCIdscal
static void FCIdscal(const unsigned int vecLength, const double alpha, double *vec)
Do lapack&#39;s dscal vec *= alpha.
Definition: FCI.cpp:1899
CheMPS2::FCI::LowestEnergyDeterminant
unsigned int LowestEnergyDeterminant() const
Return the global counter of the Slater determinant with the lowest energy.
Definition: FCI.cpp:1700
CheMPS2::FCI::CGoperator
void CGoperator(const double alpha, const double beta, const double eta, double *in, double *temp, double *out) const
Calculate out = [(alpha + beta * Hamiltonian)^2 + eta^2] * in (Econstant is taken into account!!) ...
Definition: FCI.cpp:2153
CheMPS2::FCI::getNel_down
unsigned int getNel_down() const
Getter for the number of down or beta electrons.
Definition: FCI.h:60
CheMPS2::Hamiltonian::getTmat
double getTmat(const int index1, const int index2) const
Get a Tmat element.
Definition: Hamiltonian.cpp:110
CheMPS2::FCI::getGmat
double getGmat(const int orb1, const int orb2) const
Function which returns the AUGMENTED one-body matrix elements ( G_{ij} = T_{ij} - 0...
Definition: FCI.h:88
CheMPS2::FCI::CGSolveSystem
void CGSolveSystem(const double alpha, const double beta, const double eta, double *RHS, double *RealSol, double *ImagSol, const bool checkError=true) const
Calculate the solution of the equation ( alpha + beta * Hamiltonian + I * eta ) Solution = RHS with c...
Definition: FCI.cpp:2066
CheMPS2::FCI::RetardedGF
void 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
Calculate the retarded Green&#39;s function (= addition + removal amplitude)
Definition: FCI.cpp:2191
CheMPS2::FCI::getFCIcoeff
double getFCIcoeff(int *bits_up, int *bits_down, double *vector) const
Function which returns a FCI coefficient.
Definition: FCI.cpp:435
CheMPS2::FCI::GFmatrix_addition
void 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=NULL, double **TwoRDMimag=NULL, double **TwoRDMadd=NULL) const
Calculate the addition Green&#39;s function: GF[i+numLeft*j] = <GSvector| a_{orbsLeft[i], spin} [ alpha + beta * Ham + I*eta ]^{-1} a^+_{orbsRight[j], spin} |GSvector>
Definition: FCI.cpp:2215
CheMPS2::FCI::~FCI
virtual ~FCI()
Destructor.
Definition: FCI.cpp:82
CheMPS2::Irreps::directProd
static int directProd(const int Irrep1, const int Irrep2)
Get the direct product of the irreps with numbers Irrep1 and Irrep2: a bitwise XOR for psi4&#39;s convent...
Definition: Irreps.h:123
CheMPS2::FCI::getVecLength
unsigned int getVecLength(const int irrep_center) const
Getter for the number of variables in the vector " E_ij | FCI vector > " ; where irrep_center = I_i x...
Definition: FCI.h:65
CheMPS2::FCI::Fill3RDM
void Fill3RDM(double *vector, double *ThreeRDM) const
Construct the (spin-summed) 3-RDM of a FCI vector: Gamma^3(i,j,k,l,m,n) = sum_sigma,tau,s < a^+_{i,sigma} a^+_{j,tau} a^+_{k,s} a_{n,s} a_{m,tau} a_{l,sigma} > = ThreeRDM[ i + L * ( j + L * ( k + L * ( l + L * ( m + L * n ) ) ) ) ].
Definition: FCI.cpp:1168
CheMPS2::FCI::CGAlphaPlusBetaHAM
void CGAlphaPlusBetaHAM(const double alpha, const double beta, double *in, double *out) const
Calculate out = (alpha + beta * Hamiltonian) * in (Econstant is taken into account!!) ...
Definition: FCI.cpp:2142
CheMPS2::FCI::GetMatrixElement
double GetMatrixElement(int *bits_bra_up, int *bits_bra_down, int *bits_ket_up, int *bits_ket_down, int *work) const
Sandwich the Hamiltonian between two Slater determinants (return a specific element) (without Econsta...
Definition: FCI.cpp:1722
CheMPS2::FCI::ClearVector
static void ClearVector(const unsigned int vecLength, double *vec)
Set the entries of a vector to zero.
Definition: FCI.cpp:1908
CheMPS2::FCI::ActWithNumberOperator
void ActWithNumberOperator(const unsigned int orbIndex, double *resultVector, double *sourceVector) const
Set resultVector to the number operator of a specific site acting on sourceVector.
Definition: FCI.cpp:1957
CheMPS2::FCI::RetardedGF_addition
void 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=NULL, double *TwoRDMimag=NULL, double *TwoRDMadd=NULL) const
Calculate the addition part of the retarded Green&#39;s function: <GSvector| a_{alpha, spin(isUp)} [ omega - Ham + GSenergy + I*eta ]^{-1} a^+_{beta, spin(isUp)} |GSvector>
Definition: FCI.cpp:2287
CheMPS2::FCI::str2bits
static void str2bits(const unsigned int Lvalue, const unsigned int bitstring, int *bits)
Convertor between two representations of a same spin-projection Slater determinant.
Definition: FCI.cpp:391
CheMPS2::FCI::bits2str
static unsigned int bits2str(const unsigned int Lvalue, int *bits)
Convertor between two representations of a same spin-projection Slater determinant.
Definition: FCI.cpp:397
CheMPS2::Hamiltonian::getNGroup
int getNGroup() const
Get the group number.
Definition: Hamiltonian.cpp:95
CheMPS2::FCI::getNel_up
unsigned int getNel_up() const
Getter for the number of up or alpha electrons.
Definition: FCI.h:56
CheMPS2::Davidson::FetchInstruction
char FetchInstruction(double **pointers)
The iterator to converge the ground state vector.
Definition: Davidson.cpp:105
CheMPS2::FCI::matvec
void matvec(double *input, double *output) const
Function which performs the Hamiltonian times Vector product (without Econstant!!) making use of the ...
Definition: FCI.cpp:618
CheMPS2::Hamiltonian::getEconst
double getEconst() const
Get the constant energy.
Definition: Hamiltonian.cpp:101
CheMPS2::FCI::getBitsOfCounter
void getBitsOfCounter(const int irrep_center, const unsigned int counter, int *bits_up, int *bits_down) const
Find the bit representation of a global counter corresponding to " E_ij | FCI vector > " ; where irre...
Definition: FCI.cpp:417
CheMPS2::FCI::CalcSpinSquared
double CalcSpinSquared(double *vector) const
Measure S(S+1) (spin squared)
Definition: FCI.cpp:1405
CheMPS2::FCI::getUpIrrepOfCounter
int getUpIrrepOfCounter(const int irrep_center, const unsigned int counter) const
Find the irrep of the up Slater determinant of a global counter corresponding to " E_ij | FCI vector ...
Definition: FCI.cpp:409
CheMPS2::FCI::getL
unsigned int getL() const
Getter for the number of orbitals.
Definition: FCI.h:52
CheMPS2::FCI::getTargetIrrep
int getTargetIrrep() const
Get the target irrep.
Definition: FCI.h:69
CheMPS2::FCI::FCIfrobeniusnorm
static double FCIfrobeniusnorm(const unsigned int vecLength, double *vec)
Calculate the 2-norm of a vector.
Definition: FCI.cpp:1884
CheMPS2::FCI::Fill2RDM
double Fill2RDM(double *vector, double *TwoRDM) const
Construct the (spin-summed) 2-RDM of a FCI vector: Gamma^2(i,j,k,l) = sum_sigma,tau < a^+_i...
Definition: FCI.cpp:805
CheMPS2::FCI::Fock4RDM
void Fock4RDM(double *vector, double *ThreeRDM, double *Fock, double *output) const
Construct the (spin-summed) contraction of the 4-RDM with the Fock operator: output(i,j,k,p,q,r) = sum_{l,t} Fock(l,t) * Gamma^4(i,j,k,l,p,q,r,t)
Definition: FCI.cpp:1160
CheMPS2::FCI::getOrb2Irrep
int getOrb2Irrep(const int orb) const
Get an orbital irrep.
Definition: FCI.h:74
CheMPS2::Hamiltonian::getL
int getL() const
Get the number of orbitals.
Definition: Hamiltonian.cpp:93
CheMPS2::FCI::FCI
FCI(CheMPS2::Hamiltonian *Ham, const unsigned int Nel_up, const unsigned int Nel_down, const int TargetIrrep, const double maxMemWorkMB=100.0, const int FCIverbose=2)
Constructor.
Definition: FCI.cpp:37