CheMPS2
Excitation.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 <math.h>
21 #include <stdlib.h>
22 #include <iostream>
23 #include <algorithm>
24 #include <assert.h>
25 
26 #include "Excitation.h"
27 #include "Lapack.h"
28 #include "Special.h"
29 #include "Wigner.h"
30 
31 double CheMPS2::Excitation::matvec( const SyBookkeeper * book_up, const SyBookkeeper * book_down, const int orb1, const int orb2, const double alpha, const double beta, const double gamma, Sobject * S_up, Sobject * S_down, TensorO ** overlaps, TensorL ** regular, TensorL ** trans ){
32 
33  const int indx = S_up->gIndex();
34  assert( orb1 < orb2 );
35  assert( indx == S_down->gIndex() );
36  assert( indx >= orb1 );
37  assert( indx < orb2 );
38  const int DIM = std::max( std::max( book_up->gMaxDimAtBound( indx ),
39  book_up->gMaxDimAtBound( indx + 2 ) ),
40  std::max( book_down->gMaxDimAtBound( indx ),
41  book_down->gMaxDimAtBound( indx + 2 ) ) );
42  assert( book_up->gIrrep( orb1 ) == book_up->gIrrep( orb2 ) );
43 
44  S_down->prog2symm();
45 
46  double inproduct = 0.0;
47 
48  #pragma omp parallel reduction(+:inproduct)
49  {
50  if ( orb1 + 1 == orb2 ){ // The second quantized operators are neighbours
51  #pragma omp for schedule(dynamic)
52  for ( int dummy = 0; dummy < S_up->gNKappa(); dummy++ ){
53  const int ikappa = S_up->gReorder( dummy );
54  clear( ikappa, S_up );
55  inproduct += neighbours( ikappa, book_up, book_down, alpha, beta, gamma, S_up, S_down );
56  }
57  } else {
58  double * workmem1 = new double[ DIM * DIM ];
59  if ( orb1 == indx ){ // All the way at the left
60  #pragma omp for schedule(dynamic)
61  for ( int dummy = 0; dummy < S_up->gNKappa(); dummy++ ){
62  const int ikappa = S_up->gReorder( dummy );
63  clear( ikappa, S_up );
64  first_left( ikappa, book_up, book_down, alpha, S_up, S_down, trans[ indx + 1 ] );
65  second_left( ikappa, book_up, book_down, beta, S_up, S_down, regular[ indx + 1 ] );
66  inproduct += third_left( ikappa, book_up, book_down, gamma, S_up, S_down, overlaps[ indx + 1 ], workmem1 );
67  }
68  } else {
69  if ( orb2 == indx + 1 ){ // All the way at the right
70  #pragma omp for schedule(dynamic)
71  for ( int dummy = 0; dummy < S_up->gNKappa(); dummy++ ){
72  const int ikappa = S_up->gReorder( dummy );
73  clear( ikappa, S_up );
74  first_right( ikappa, book_up, book_down, alpha, S_up, S_down, trans[ indx - 1 ] );
75  second_right( ikappa, book_up, book_down, beta, S_up, S_down, regular[ indx - 1 ] );
76  inproduct += third_right( ikappa, book_up, book_down, gamma, S_up, S_down, overlaps[ indx - 1 ], workmem1 );
77  }
78  } else { // Somewhere in the middle
79  double * workmem2 = new double[ DIM * DIM ];
80  #pragma omp for schedule(dynamic)
81  for ( int dummy = 0; dummy < S_up->gNKappa(); dummy++ ){
82  const int ikappa = S_up->gReorder( dummy );
83  clear( ikappa, S_up );
84  first_middle( ikappa, book_up, book_down, alpha, S_up, S_down, trans[ indx - 1 ], trans[ indx + 1 ], workmem1 );
85  second_middle( ikappa, book_up, book_down, beta, S_up, S_down, regular[ indx - 1 ], regular[ indx + 1 ], workmem1 );
86  inproduct += third_middle( ikappa, book_up, book_down, gamma, S_up, S_down, overlaps[ indx - 1 ], overlaps[ indx + 1 ], workmem1, workmem2 );
87  }
88  delete [] workmem2;
89  }
90  }
91  delete [] workmem1;
92  }
93  }
94 
95  S_up ->symm2prog();
96  // S_down->symm2prog(); S_down is not used anymore afterwards
97 
98  return inproduct;
99 
100 }
101 
102 void CheMPS2::Excitation::clear( const int ikappa, Sobject * S_up ){
103 
104  const int start = S_up->gKappa2index( ikappa );
105  const int stop = S_up->gKappa2index( ikappa + 1 );
106  double * storage = S_up->gStorage();
107  for ( int cnt = start; cnt < stop; cnt++ ){ storage[ cnt ] = 0.0; }
108 
109 }
110 
111 double CheMPS2::Excitation::neighbours( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double alpha, const double beta, const double gamma, Sobject * S_up, Sobject * S_down ){
112 
113  const int index = S_up->gIndex();
114  const int TwoSL = S_up->gTwoSL( ikappa );
115  const int TwoSR = S_up->gTwoSR( ikappa );
116  const int TwoJ = S_up->gTwoJ( ikappa );
117  const int NL = S_up->gNL( ikappa );
118  const int NR = S_up->gNR( ikappa );
119  const int IL = S_up->gIL( ikappa );
120  const int IR = S_up->gIR( ikappa );
121  const int N1 = S_up->gN1( ikappa );
122  const int N2 = S_up->gN2( ikappa );
123 
124  const int dimLup = book_up ->gCurrentDim( index, NL, TwoSL, IL );
125  const int dimRup = book_up ->gCurrentDim( index + 2, NR, TwoSR, IR );
126  const int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
127  const int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
128  assert( dimLup == dimLdown );
129  assert( dimRup == dimRdown );
130  assert( book_up->gIrrep( index ) == book_up->gIrrep( index + 1 ) );
131 
132  double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
133  int size = dimLup * dimRup;
134  int inc1 = 1;
135 
136  // Add a^+ a
137  if ( fabs( alpha ) > 0.0 ){
138  if (( TwoJ == 0 ) && ( N1 == 1 ) && ( N2 == 1 )){
139  double * block_down = S_down->gStorage( NL, TwoSL, IL, 0, 2, TwoJ, NR, TwoSR, IR );
140  assert( block_down != NULL );
141  double factor = sqrt( 2.0 ) * alpha;
142  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
143  }
144  if (( TwoJ == 1 ) && ( N1 == 2 ) && ( N2 == 1 )){
145  double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 2, TwoJ, NR, TwoSR, IR );
146  assert( block_down != NULL );
147  double factor = -alpha;
148  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
149  }
150  if (( TwoJ == 1 ) && ( N1 == 1 ) && ( N2 == 0 )){
151  double * block_down = S_down->gStorage( NL, TwoSL, IL, 0, 1, TwoJ, NR, TwoSR, IR );
152  assert( block_down != NULL );
153  double factor = alpha;
154  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
155  }
156  if (( TwoJ == 0 ) && ( N1 == 2 ) && ( N2 == 0 )){
157  double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 1, TwoJ, NR, TwoSR, IR );
158  assert( block_down != NULL );
159  double factor = sqrt( 2.0 ) * alpha;
160  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
161  }
162  }
163 
164  // Add a a^+
165  if ( fabs( beta ) > 0.0 ){
166  if (( TwoJ == 0 ) && ( N1 == 1 ) && ( N2 == 1 )){
167  double * block_down = S_down->gStorage( NL, TwoSL, IL, 2, 0, TwoJ, NR, TwoSR, IR );
168  assert( block_down != NULL );
169  double factor = sqrt( 2.0 ) * beta;
170  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
171  }
172  if (( TwoJ == 1 ) && ( N1 == 0 ) && ( N2 == 1 )){
173  double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 0, TwoJ, NR, TwoSR, IR );
174  assert( block_down != NULL );
175  double factor = beta;
176  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
177  }
178  if (( TwoJ == 1 ) && ( N1 == 1 ) && ( N2 == 2 )){
179  double * block_down = S_down->gStorage( NL, TwoSL, IL, 2, 1, TwoJ, NR, TwoSR, IR );
180  assert( block_down != NULL );
181  double factor = -beta;
182  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
183  }
184  if (( TwoJ == 0 ) && ( N1 == 0 ) && ( N2 == 2 )){
185  double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 1, TwoJ, NR, TwoSR, IR );
186  assert( block_down != NULL );
187  double factor = sqrt( 2.0 ) * beta;
188  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
189  }
190  }
191 
192  // Add the constant part
193  double * block_down = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
194  assert( block_down != NULL );
195  if ( fabs( gamma ) > 0.0 ){
196  double factor = gamma;
197  daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
198  }
199  const double inproduct = ddot_( &size, block_down, &inc1, block_up, &inc1 );
200  return inproduct;
201 
202 }
203 
204 void CheMPS2::Excitation::first_left( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double alpha, Sobject * S_up, Sobject * S_down, TensorL * Rtrans ){
205 
206  const int index = S_up->gIndex();
207  const int TwoSL = S_up->gTwoSL( ikappa );
208  const int TwoSR = S_up->gTwoSR( ikappa );
209  const int TwoJ = S_up->gTwoJ( ikappa );
210  const int NL = S_up->gNL( ikappa );
211  const int NR = S_up->gNR( ikappa );
212  const int IL = S_up->gIL( ikappa );
213  const int IR = S_up->gIR( ikappa );
214  const int N1 = S_up->gN1( ikappa );
215  const int N2 = S_up->gN2( ikappa );
216 
217  const int IRdown = Irreps::directProd( IR, book_up->gIrrep( index ) );
218  const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
219 
220  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
221  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
222  const int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
223  assert( dimLup == dimLdown );
224  assert( book_up->gIrrep( index ) == Rtrans->get_irrep() );
225 
226  if (( N1 == 1 ) && ( fabs( alpha ) > 0.0 )){ // Cfr. 3K1A in HeffDiagrams3.cpp
227  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
228  if (( abs( TwoSL - TwoSRdown ) <= TwoS2 ) && ( TwoSRdown >= 0 )){
229  const int memSkappa = S_down->gKappa( NL, TwoSL, IL, 0, N2, TwoS2, NR - 1, TwoSRdown, IRdown );
230  if ( memSkappa != -1 ){
231  int dimRdown = book_down->gCurrentDim( index + 2, NR - 1, TwoSRdown, IRdown );
232  double factor = alpha
233  * Special::phase( TwoSL + TwoSR + TwoJ + 2 * TwoS2 )
234  * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoSR + 1 ) )
235  * Wigner::wigner6j( TwoS2, TwoJ, 1, TwoSR, TwoSRdown, TwoSL );
236  double add = 1.0;
237  char notrans = 'N';
238  double * block_right = Rtrans->gStorage( NR - 1, TwoSRdown, IRdown, NR, TwoSR, IR );
239  double * block_down = S_down->gStorage() + S_down->gKappa2index( memSkappa );
240  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
241  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRdown, &add, block_up, &dimLup );
242  }
243  }
244  }
245  }
246 
247  if (( N1 == 2 ) && ( fabs( alpha ) > 0.0 )){ // Cfr. 3K1B in HeffDiagrams3.cpp
248  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
249  int dimRdown = book_down->gCurrentDim( index + 2, NR - 1, TwoSRdown, IRdown );
250  if (( dimRdown > 0 ) && ( TwoSRdown >= 0 )){
251  const int TwoJstart = ((( TwoSRdown != TwoSL ) || ( TwoS2 == 0 )) ? 1 + TwoS2 : 0 );
252  for ( int TwoJdown = TwoJstart; TwoJdown <= 1 + TwoS2; TwoJdown += 2 ){
253  if ( abs( TwoSL - TwoSRdown ) <= TwoJdown ){
254  const int memSkappa = S_down->gKappa( NL, TwoSL, IL, 1, N2, TwoJdown, NR - 1, TwoSRdown, IRdown );
255  if ( memSkappa != -1 ){
256  double factor = alpha
257  * Special::phase( TwoSL + TwoSR + TwoJdown + 1 + 2 * TwoS2 )
258  * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSR + 1 ) )
259  * Wigner::wigner6j( TwoJdown, TwoS2, 1, TwoSR, TwoSRdown, TwoSL );
260  double add = 1.0;
261  char notrans = 'N';
262  double * block_right = Rtrans->gStorage( NR - 1, TwoSRdown, IRdown, NR, TwoSR, IR );
263  double * block_down = S_down->gStorage() + S_down->gKappa2index( memSkappa );
264  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
265  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRdown, &add, block_up, &dimLup );
266  }
267  }
268  }
269  }
270  }
271  }
272 
273 }
274 
275 void CheMPS2::Excitation::second_left( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double beta, Sobject * S_up, Sobject * S_down, TensorL * Rregular ){
276 
277  const int index = S_up->gIndex();
278  const int TwoSL = S_up->gTwoSL( ikappa );
279  const int TwoSR = S_up->gTwoSR( ikappa );
280  const int TwoJ = S_up->gTwoJ( ikappa );
281  const int NL = S_up->gNL( ikappa );
282  const int NR = S_up->gNR( ikappa );
283  const int IL = S_up->gIL( ikappa );
284  const int IR = S_up->gIR( ikappa );
285  const int N1 = S_up->gN1( ikappa );
286  const int N2 = S_up->gN2( ikappa );
287 
288  const int IRdown = Irreps::directProd( IR, book_up->gIrrep( index ) );
289  const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
290 
291  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
292  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
293  const int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
294  assert( dimLup == dimLdown );
295  assert( book_up->gIrrep( index ) == Rregular->get_irrep() );
296 
297  if (( N1 == 0 ) && ( fabs( beta ) > 0.0 )){ // Cfr. 3K2A in HeffDiagrams3.cpp
298  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
299  int dimRdown = book_down->gCurrentDim( index + 2, NR + 1, TwoSRdown, IRdown );
300  if (( dimRdown > 0 ) && ( TwoSRdown >= 0 )){
301  const int TwoJstart = ((( TwoSRdown != TwoSL ) || ( TwoS2 == 0 )) ? 1 + TwoS2 : 0 );
302  for ( int TwoJdown = TwoJstart; TwoJdown <= 1 + TwoS2; TwoJdown += 2 ){
303  if ( abs( TwoSL - TwoSRdown ) <= TwoJdown ){
304  const int memSkappa = S_down->gKappa( NL, TwoSL, IL, 1, N2, TwoJdown, NR + 1, TwoSRdown, IRdown );
305  if ( memSkappa != -1 ){
306  double factor = beta
307  * Special::phase( TwoSL + TwoSRdown + TwoJdown + 2 * TwoS2 )
308  * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSRdown + 1 ) )
309  * Wigner::wigner6j( TwoJdown, TwoS2, 1, TwoSR, TwoSRdown, TwoSL );
310  double add = 1.0;
311  char notrans = 'N';
312  char trans = 'T';
313  double * block_right = Rregular->gStorage( NR, TwoSR, IR, NR + 1, TwoSRdown, IRdown );
314  double * block_down = S_down ->gStorage() + S_down->gKappa2index( memSkappa );
315  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
316  dgemm_( &notrans, &trans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRup, &add, block_up, &dimLup );
317  }
318  }
319  }
320  }
321  }
322  }
323 
324  if (( N1 == 1 ) && ( fabs( beta ) > 0.0 )){ // Cfr. 3K2B in HeffDiagrams3.cpp
325  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
326  if (( abs( TwoSL - TwoSRdown ) <= TwoS2 ) && ( TwoSRdown >= 0 )){
327  const int memSkappa = S_down->gKappa( NL, TwoSL, IL, 2, N2, TwoS2, NR + 1, TwoSRdown, IRdown );
328  if ( memSkappa != -1 ){
329  int dimRdown = book_down->gCurrentDim( index + 2, NR + 1, TwoSRdown, IRdown );
330  double factor = beta
331  * Special::phase( TwoSL + TwoSRdown + TwoJ + 1 + 2 * TwoS2 )
332  * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoSRdown + 1 ) )
333  * Wigner::wigner6j( TwoS2, TwoJ, 1, TwoSR, TwoSRdown, TwoSL );
334  double add = 1.0;
335  char notrans = 'N';
336  char trans = 'T';
337  double * block_right = Rregular->gStorage( NR, TwoSR, IR, NR + 1, TwoSRdown, IRdown );
338  double * block_down = S_down ->gStorage() + S_down->gKappa2index( memSkappa );
339  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
340  dgemm_( &notrans, &trans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRup, &add, block_up, &dimLup );
341  }
342  }
343  }
344  }
345 
346 }
347 
348 double CheMPS2::Excitation::third_left( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double gamma, Sobject * S_up, Sobject * S_down, TensorO * Rovlp, double * workmem ){
349 
350  const int index = S_up->gIndex();
351  const int TwoSL = S_up->gTwoSL( ikappa );
352  const int TwoSR = S_up->gTwoSR( ikappa );
353  const int TwoJ = S_up->gTwoJ( ikappa );
354  const int NL = S_up->gNL( ikappa );
355  const int NR = S_up->gNR( ikappa );
356  const int IL = S_up->gIL( ikappa );
357  const int IR = S_up->gIR( ikappa );
358  const int N1 = S_up->gN1( ikappa );
359  const int N2 = S_up->gN2( ikappa );
360 
361  int dimLup = book_up ->gCurrentDim( index, NL, TwoSL, IL );
362  int dimRup = book_up ->gCurrentDim( index + 2, NR, TwoSR, IR );
363  int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
364  int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
365  assert( dimLup == dimLdown );
366 
367  double inproduct = 0.0;
368  if ( dimRdown > 0 ){
369  double * block_down = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
370  double * block_right = Rovlp ->gStorage( NR, TwoSR, IR, NR, TwoSR, IR );
371  char trans = 'T';
372  char notrans = 'N';
373  double one = 1.0;
374  double set = 0.0;
375  dgemm_( &notrans, &trans, &dimLup, &dimRup, &dimRdown, &one, block_down, &dimLup, block_right, &dimRup, &set, workmem, &dimLup );
376 
377  double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
378  int size = dimLup * dimRup;
379  int inc1 = 1;
380  if ( fabs( gamma ) > 0.0 ){
381  double factor = gamma;
382  daxpy_( &size, &factor, workmem, &inc1, block_up, &inc1 );
383  }
384  inproduct = ddot_( &size, workmem, &inc1, block_up, &inc1 );
385  }
386  return inproduct;
387 
388 }
389 
390 void CheMPS2::Excitation::first_right( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double alpha, Sobject * S_up, Sobject * S_down, TensorL * Ltrans ){
391 
392  const int index = S_up->gIndex();
393  const int TwoSL = S_up->gTwoSL( ikappa );
394  const int TwoSR = S_up->gTwoSR( ikappa );
395  const int TwoJ = S_up->gTwoJ( ikappa );
396  const int NL = S_up->gNL( ikappa );
397  const int NR = S_up->gNR( ikappa );
398  const int IL = S_up->gIL( ikappa );
399  const int IR = S_up->gIR( ikappa );
400  const int N1 = S_up->gN1( ikappa );
401  const int N2 = S_up->gN2( ikappa );
402 
403  const int ILdown = Irreps::directProd( IL, book_up->gIrrep( index + 1 ) );
404  const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
405 
406  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
407  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
408  const int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
409  assert( dimRup == dimRdown );
410  assert( book_up->gIrrep( index + 1 ) == Ltrans->get_irrep() );
411 
412  if (( N2 == 0 ) && ( fabs( alpha ) > 0.0 )){ // Cfr. 3B2A in HeffDiagrams3.cpp
413  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
414  int dimLdown = book_down->gCurrentDim( index, NL - 1, TwoSLdown, ILdown );
415  if (( dimLdown > 0 ) && ( TwoSLdown >= 0 )){
416  const int TwoJstart = ((( TwoSR != TwoSLdown ) || ( TwoS1 == 0 )) ? 1 + TwoS1 : 0 );
417  for ( int TwoJdown = TwoJstart; TwoJdown <= 1 + TwoS1; TwoJdown += 2 ){
418  if ( abs( TwoSLdown - TwoSR ) <= TwoJdown ){
419  const int memSkappa = S_down->gKappa( NL - 1, TwoSLdown, ILdown, N1, 1, TwoJdown, NR, TwoSR, IR );
420  if ( memSkappa != -1 ){
421  double factor = alpha
422  * Special::phase( TwoSLdown + TwoSR + 2 - TwoJdown )
423  * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoJdown + 1 ) )
424  * Wigner::wigner6j( TwoJdown, TwoS1, 1, TwoSL, TwoSLdown, TwoSR );
425  double add = 1.0;
426  char notrans = 'N';
427  char trans = 'T';
428  double * block_left = Ltrans->gStorage( NL - 1, TwoSLdown, ILdown, NL, TwoSL, IL );
429  double * block_down = S_down->gStorage() + S_down->gKappa2index( memSkappa );
430  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
431  dgemm_( &trans, &notrans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &add, block_up, &dimLup );
432  }
433  }
434  }
435  }
436  }
437  }
438 
439  if (( N2 == 1 ) && ( fabs( alpha ) > 0.0 )){ // Cfr. 3B2B in HeffDiagrams3.cpp
440  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
441  if (( abs( TwoSLdown - TwoSR ) <= TwoS1 ) && ( TwoSLdown >= 0 )){
442  const int memSkappa = S_down->gKappa( NL - 1, TwoSLdown, ILdown, N1, 2, TwoS1, NR, TwoSR, IR );
443  if ( memSkappa != -1 ){
444  int dimLdown = book_down->gCurrentDim( index, NL - 1, TwoSLdown, ILdown );
445  double factor = alpha
446  * Special::phase( TwoSLdown + TwoSR + 3 - TwoJ )
447  * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoJ + 1 ) )
448  * Wigner::wigner6j( TwoS1, TwoJ, 1, TwoSL, TwoSLdown, TwoSR );
449  double add = 1.0;
450  char notrans = 'N';
451  char trans = 'T';
452  double * block_left = Ltrans->gStorage( NL - 1, TwoSLdown, ILdown, NL, TwoSL, IL );
453  double * block_down = S_down->gStorage() + S_down->gKappa2index( memSkappa );
454  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
455  dgemm_( &trans, &notrans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &add, block_up, &dimLup );
456  }
457  }
458  }
459  }
460 
461 }
462 
463 void CheMPS2::Excitation::second_right( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double beta, Sobject * S_up, Sobject * S_down, TensorL * Lregular ){
464 
465  const int index = S_up->gIndex();
466  const int TwoSL = S_up->gTwoSL( ikappa );
467  const int TwoSR = S_up->gTwoSR( ikappa );
468  const int TwoJ = S_up->gTwoJ( ikappa );
469  const int NL = S_up->gNL( ikappa );
470  const int NR = S_up->gNR( ikappa );
471  const int IL = S_up->gIL( ikappa );
472  const int IR = S_up->gIR( ikappa );
473  const int N1 = S_up->gN1( ikappa );
474  const int N2 = S_up->gN2( ikappa );
475 
476  const int ILdown = Irreps::directProd( IL, book_up->gIrrep( index + 1 ) );
477  const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
478 
479  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
480  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
481  const int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
482  assert( dimRup == dimRdown );
483  assert( book_up->gIrrep( index + 1 ) == Lregular->get_irrep() );
484 
485  if (( N2 == 2 ) && ( fabs( beta ) > 0.0 )){ // Cfr. 3B1A in HeffDiagrams3.cpp
486  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
487  int dimLdown = book_down->gCurrentDim( index, NL + 1, TwoSLdown, ILdown );
488  if (( dimLdown > 0 ) && ( TwoSLdown >= 0 )){
489  const int TwoJstart = ((( TwoSR != TwoSLdown ) || ( TwoS1 == 0 )) ? 1 + TwoS1 : 0 );
490  for ( int TwoJdown = TwoJstart; TwoJdown <= 1 + TwoS1; TwoJdown += 2 ){
491  if ( abs( TwoSLdown - TwoSR ) <= TwoJdown ){
492  const int memSkappa = S_down->gKappa( NL + 1, TwoSLdown, ILdown, N1, 1, TwoJdown, NR, TwoSR, IR );
493  if ( memSkappa != -1 ){
494  double factor = beta
495  * Special::phase( TwoSL + TwoSR + 3 - TwoJdown )
496  * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSLdown + 1 ) )
497  * Wigner::wigner6j( TwoJdown, TwoS1, 1, TwoSL, TwoSLdown, TwoSR );
498  double add = 1.0;
499  char notrans = 'N';
500  double * block_left = Lregular->gStorage( NL, TwoSL, IL, NL + 1, TwoSLdown, ILdown );
501  double * block_down = S_down ->gStorage() + S_down->gKappa2index( memSkappa );
502  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
503  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &add, block_up, &dimLup );
504  }
505  }
506  }
507  }
508  }
509  }
510 
511  if (( N2 == 1 ) && ( fabs( beta ) > 0.0 )){ // Cfr. 3B1B in HeffDiagrams3.cpp
512  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
513  if (( abs( TwoSLdown - TwoSR ) <= TwoS1 ) && ( TwoSLdown >= 0 )){
514  const int memSkappa = S_down->gKappa( NL + 1, TwoSLdown, ILdown, N1, 0, TwoS1, NR, TwoSR, IR );
515  if ( memSkappa != -1 ){
516  int dimLdown = book_down->gCurrentDim( index, NL + 1, TwoSLdown, ILdown );
517  double factor = beta
518  * Special::phase( TwoSL + TwoSR + 2 - TwoJ )
519  * sqrt( 1.0 * ( TwoSLdown + 1 ) * ( TwoJ + 1 ) )
520  * Wigner::wigner6j( TwoS1, TwoJ, 1, TwoSL, TwoSLdown, TwoSR );
521  double add = 1.0;
522  char notrans = 'N';
523  double * block_left = Lregular->gStorage( NL, TwoSL, IL, NL + 1, TwoSLdown, ILdown );
524  double * block_down = S_down ->gStorage() + S_down->gKappa2index( memSkappa );
525  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
526  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &add, block_up, &dimLup );
527  }
528  }
529  }
530  }
531 
532 }
533 
534 double CheMPS2::Excitation::third_right( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double gamma, Sobject * S_up, Sobject * S_down, TensorO * Lovlp, double * workmem ){
535 
536  const int index = S_up->gIndex();
537  const int TwoSL = S_up->gTwoSL( ikappa );
538  const int TwoSR = S_up->gTwoSR( ikappa );
539  const int TwoJ = S_up->gTwoJ( ikappa );
540  const int NL = S_up->gNL( ikappa );
541  const int NR = S_up->gNR( ikappa );
542  const int IL = S_up->gIL( ikappa );
543  const int IR = S_up->gIR( ikappa );
544  const int N1 = S_up->gN1( ikappa );
545  const int N2 = S_up->gN2( ikappa );
546 
547  int dimLup = book_up ->gCurrentDim( index, NL, TwoSL, IL );
548  int dimRup = book_up ->gCurrentDim( index + 2, NR, TwoSR, IR );
549  int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
550  int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
551  assert( dimRup == dimRdown );
552 
553  double inproduct = 0.0;
554  if ( dimLdown > 0 ){
555  double * block_down = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
556  double * block_left = Lovlp ->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
557  char notrans = 'N';
558  double one = 1.0;
559  double set = 0.0;
560  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimLdown, &one, block_left, &dimLup, block_down, &dimLdown, &set, workmem, &dimLup );
561 
562  double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
563  int size = dimLup * dimRup;
564  int inc1 = 1;
565  if ( fabs( gamma ) > 0.0 ){
566  double factor = gamma;
567  daxpy_( &size, &factor, workmem, &inc1, block_up, &inc1 );
568  }
569  inproduct = ddot_( &size, workmem, &inc1, block_up, &inc1 );
570  }
571  return inproduct;
572 
573 }
574 
575 void CheMPS2::Excitation::first_middle( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double alpha, Sobject * S_up, Sobject * S_down, TensorL * Ltrans, TensorL * Rtrans, double * workmem ){
576 
577  const int index = S_up->gIndex();
578  const int TwoSL = S_up->gTwoSL( ikappa );
579  const int TwoSR = S_up->gTwoSR( ikappa );
580  const int TwoJ = S_up->gTwoJ( ikappa );
581  const int NL = S_up->gNL( ikappa );
582  const int NR = S_up->gNR( ikappa );
583  const int IL = S_up->gIL( ikappa );
584  const int IR = S_up->gIR( ikappa );
585  const int N1 = S_up->gN1( ikappa );
586  const int N2 = S_up->gN2( ikappa );
587 
588  const int ILdown = Irreps::directProd( IL, Ltrans->get_irrep() );
589  const int IRdown = Irreps::directProd( IR, Rtrans->get_irrep() );
590  const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
591  const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
592 
593  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
594  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
595  assert( Ltrans->get_irrep() == Rtrans->get_irrep() );
596 
597  if ( fabs( alpha ) > 0.0 ){ // Cfr. 3C2 in HeffDiagrams3.cpp
598  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
599  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
600  if (( abs( TwoSLdown - TwoSRdown ) <= TwoJ ) && ( TwoSLdown >= 0 ) && ( TwoSRdown >= 0 )){
601  const int memSkappa = S_down->gKappa( NL - 1, TwoSLdown, ILdown, N1, N2, TwoJ, NR - 1, TwoSRdown, IRdown );
602  if ( memSkappa != -1 ){
603  int dimLdown = book_down->gCurrentDim( index, NL - 1, TwoSLdown, ILdown );
604  int dimRdown = book_down->gCurrentDim( index + 2, NR - 1, TwoSRdown, IRdown );
605  double factor = alpha
606  * Special::phase( TwoSL + TwoSRdown + TwoJ + 1 + 2 * TwoS1 + 2 * TwoS2 )
607  * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSR + 1 ) )
608  * Wigner::wigner6j( TwoSL, TwoSR, TwoJ, TwoSRdown, TwoSLdown, 1 );
609  char trans = 'T';
610  char notrans = 'N';
611  double set = 0.0;
612  double one = 1.0;
613  double * block_left = Ltrans->gStorage( NL - 1, TwoSLdown, ILdown, NL, TwoSL, IL );
614  double * block_right = Rtrans->gStorage( NR - 1, TwoSRdown, IRdown, NR, TwoSR, IR );
615  double * block_down = S_down->gStorage() + S_down->gKappa2index( memSkappa );
616  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
617  dgemm_( &trans, &notrans, &dimLup, &dimRdown, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &set, workmem, &dimLup );
618  dgemm_( &notrans, &notrans, &dimLup, &dimRup, &dimRdown, &one, workmem, &dimLup, block_right, &dimRdown, &one, block_up, &dimLup );
619  }
620  }
621  }
622  }
623  }
624 
625 }
626 
627 void CheMPS2::Excitation::second_middle( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double beta, Sobject * S_up, Sobject * S_down, TensorL * Lregular, TensorL * Rregular, double * workmem ){
628 
629  const int index = S_up->gIndex();
630  const int TwoSL = S_up->gTwoSL( ikappa );
631  const int TwoSR = S_up->gTwoSR( ikappa );
632  const int TwoJ = S_up->gTwoJ( ikappa );
633  const int NL = S_up->gNL( ikappa );
634  const int NR = S_up->gNR( ikappa );
635  const int IL = S_up->gIL( ikappa );
636  const int IR = S_up->gIR( ikappa );
637  const int N1 = S_up->gN1( ikappa );
638  const int N2 = S_up->gN2( ikappa );
639 
640  const int ILdown = Irreps::directProd( IL, Lregular->get_irrep() );
641  const int IRdown = Irreps::directProd( IR, Rregular->get_irrep() );
642  const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
643  const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
644 
645  int dimLup = book_up->gCurrentDim( index, NL, TwoSL, IL );
646  int dimRup = book_up->gCurrentDim( index + 2, NR, TwoSR, IR );
647  assert( Lregular->get_irrep() == Rregular->get_irrep() );
648 
649  if ( fabs( beta ) > 0.0 ){ // Cfr. 3C1 in HeffDiagrams3.cpp
650  for ( int TwoSLdown = TwoSL - 1; TwoSLdown <= TwoSL + 1; TwoSLdown += 2 ){
651  for ( int TwoSRdown = TwoSR - 1; TwoSRdown <= TwoSR + 1; TwoSRdown += 2 ){
652  if (( abs( TwoSLdown - TwoSRdown ) <= TwoJ ) && ( TwoSLdown >= 0 ) && ( TwoSRdown >= 0 )){
653  const int memSkappa = S_down->gKappa( NL + 1, TwoSLdown, ILdown, N1, N2, TwoJ, NR + 1, TwoSRdown, IRdown );
654  if ( memSkappa != -1 ){
655  int dimLdown = book_down->gCurrentDim( index, NL + 1, TwoSLdown, ILdown );
656  int dimRdown = book_down->gCurrentDim( index + 2, NR + 1, TwoSRdown, IRdown );
657  double factor = beta
658  * Special::phase( TwoSLdown + TwoSR + TwoJ + 1 + 2 * TwoS1 + 2 * TwoS2 )
659  * sqrt( 1.0 * ( TwoSLdown + 1 ) * ( TwoSRdown + 1 ) )
660  * Wigner::wigner6j( TwoSL, TwoSR, TwoJ, TwoSRdown, TwoSLdown, 1 );
661  char trans = 'T';
662  char notrans = 'N';
663  double set = 0.0;
664  double one = 1.0;
665  double * block_left = Lregular->gStorage( NL, TwoSL, IL, NL + 1, TwoSLdown, ILdown );
666  double * block_right = Rregular->gStorage( NR, TwoSR, IR, NR + 1, TwoSRdown, IRdown );
667  double * block_down = S_down ->gStorage() + S_down->gKappa2index( memSkappa );
668  double * block_up = S_up ->gStorage() + S_up ->gKappa2index( ikappa );
669  dgemm_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &set, workmem, &dimLup );
670  dgemm_( &notrans, &trans, &dimLup, &dimRup, &dimRdown, &one, workmem, &dimLup, block_right, &dimRup, &one, block_up, &dimLup );
671  }
672  }
673  }
674  }
675  }
676 
677 }
678 
679 double CheMPS2::Excitation::third_middle( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double gamma, Sobject * S_up, Sobject * S_down, TensorO * Lovlp, TensorO * Rovlp, double * workmem1, double * workmem2 ){
680 
681  const int index = S_up->gIndex();
682  const int TwoSL = S_up->gTwoSL( ikappa );
683  const int TwoSR = S_up->gTwoSR( ikappa );
684  const int TwoJ = S_up->gTwoJ( ikappa );
685  const int NL = S_up->gNL( ikappa );
686  const int NR = S_up->gNR( ikappa );
687  const int IL = S_up->gIL( ikappa );
688  const int IR = S_up->gIR( ikappa );
689  const int N1 = S_up->gN1( ikappa );
690  const int N2 = S_up->gN2( ikappa );
691 
692  int dimLup = book_up ->gCurrentDim( index, NL, TwoSL, IL );
693  int dimRup = book_up ->gCurrentDim( index + 2, NR, TwoSR, IR );
694  int dimLdown = book_down->gCurrentDim( index, NL, TwoSL, IL );
695  int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
696 
697  double inproduct = 0.0;
698  if (( dimLdown > 0 ) && ( dimRdown > 0 )){
699  double * block_down = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
700  double * block_left = Lovlp ->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
701  double * block_right = Rovlp ->gStorage( NR, TwoSR, IR, NR, TwoSR, IR );
702  char trans = 'T';
703  char notrans = 'N';
704  double one = 1.0;
705  double set = 0.0;
706  dgemm_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &one, block_left, &dimLup, block_down, &dimLdown, &set, workmem1, &dimLup );
707  dgemm_( &notrans, &trans, &dimLup, &dimRup, &dimRdown, &one, workmem1, &dimLup, block_right, &dimRup, &set, workmem2, &dimLup );
708 
709  double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
710  int size = dimLup * dimRup;
711  int inc1 = 1;
712  if ( fabs( gamma ) > 0.0 ){
713  double factor = gamma;
714  daxpy_( &size, &factor, workmem2, &inc1, block_up, &inc1 );
715  }
716  inproduct = ddot_( &size, workmem2, &inc1, block_up, &inc1 );
717  }
718  return inproduct;
719 
720 }
721 
722 
int gTwoSL(const int ikappa) const
Get the left spin symmetry of block ikappa.
Definition: Sobject.cpp:203
static double wigner6j(const int two_ja, const int two_jb, const int two_jc, const int two_jd, const int two_je, const int two_jf)
Wigner-6j symbol (gsl api)
Definition: Wigner.cpp:294
static double matvec(const SyBookkeeper *book_up, const SyBookkeeper *book_down, const int orb1, const int orb2, const double alpha, const double beta, const double gamma, Sobject *S_up, Sobject *S_down, TensorO **overlaps, TensorL **regular, TensorL **trans)
Matrix-vector multiplication S_up = ( alpha * E_{orb1,orb2} + beta * E_{orb2,orb1} + gamma ) x S_down...
Definition: Excitation.cpp:31
int gCurrentDim(const int boundary, const int N, const int TwoS, const int irrep) const
Get the current virtual dimensions.
int gKappa2index(const int kappa) const
Get the storage jump corresponding to a certain symmetry block.
Definition: Sobject.cpp:190
void symm2prog()
Convert the storage from symmetric Hamiltonian convention to diagram convention.
Definition: Sobject.cpp:636
static int phase(const int power_times_two)
Phase function.
Definition: Special.h:36
int gTwoJ(const int ikappa) const
Get the central spin symmetry of block ikappa.
Definition: Sobject.cpp:207
int gNKappa() const
Get the number of symmetry blocks.
Definition: Sobject.cpp:166
double * gStorage()
Get the pointer to the storage.
int gIL(const int ikappa) const
Get the left irrep symmetry of block ikappa.
Definition: Sobject.cpp:204
double * gStorage()
Get the pointer to the storage.
Definition: Sobject.cpp:168
int gKappa(const int NL, const int TwoSL, const int IL, const int N1, const int N2, const int TwoJ, const int NR, const int TwoSR, const int IR) const
Get the index corresponding to a certain symmetry block.
Definition: Sobject.cpp:172
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
void prog2symm()
Convert the storage from diagram convention to symmetric Hamiltonian convention.
Definition: Sobject.cpp:622
int gNL(const int ikappa) const
Get the left particle number symmetry of block ikappa.
Definition: Sobject.cpp:202
int gIR(const int ikappa) const
Get the right irrep symmetry of block ikappa.
Definition: Sobject.cpp:210
int gTwoSR(const int ikappa) const
Get the right spin symmetry of block ikappa.
Definition: Sobject.cpp:209
int gIndex() const
Get the location index.
Definition: Sobject.cpp:200
int gMaxDimAtBound(const int boundary) const
Get the maximum virtual dimension at a certain boundary.
int gReorder(const int ikappa) const
Get the blocks from large to small: blocksize(reorder[i])>=blocksize(reorder[i+1]) ...
Definition: Sobject.cpp:170
int gNR(const int ikappa) const
Get the right particle number symmetry of block ikappa.
Definition: Sobject.cpp:208
int get_irrep() const
Get the (real-valued abelian) point group irrep difference between the symmetry sectors of the lower ...
int gIrrep(const int orbital) const
Get an orbital irrep.
int gN1(const int ikappa) const
Get the left local particle number symmetry of block ikappa.
Definition: Sobject.cpp:205
int gN2(const int ikappa) const
Get the right local particle number symmetry of block ikappa.
Definition: Sobject.cpp:206