26 #include "Excitation.h" 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 ){
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 ),
42 assert( book_up->
gIrrep( orb1 ) == book_up->
gIrrep( orb2 ) );
46 double inproduct = 0.0;
48 #pragma omp parallel reduction(+:inproduct) 50 if ( orb1 + 1 == orb2 ){
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 );
58 double * workmem1 =
new double[ DIM * DIM ];
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 );
69 if ( orb2 == indx + 1 ){
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 );
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 );
102 void CheMPS2::Excitation::clear(
const int ikappa,
Sobject * S_up ){
106 double * storage = S_up->
gStorage();
107 for (
int cnt = start; cnt < stop; cnt++ ){ storage[ cnt ] = 0.0; }
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 ){
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 );
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 ) );
133 int size = dimLup * dimRup;
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
199 const double inproduct = ddot_( &size, block_down, &inc1, block_up, &inc1 );
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 );
218 const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
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 );
226 if (( N1 == 1 ) && ( fabs( alpha ) > 0.0 )){
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
234 * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoSR + 1 ) )
238 double * block_right = Rtrans->
gStorage( NR - 1, TwoSRdown, IRdown, NR, TwoSR, IR );
241 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRdown, &add, block_up, &dimLup );
247 if (( N1 == 2 ) && ( fabs( alpha ) > 0.0 )){
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
258 * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSR + 1 ) )
262 double * block_right = Rtrans->
gStorage( NR - 1, TwoSRdown, IRdown, NR, TwoSR, IR );
265 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRdown, &add, block_up, &dimLup );
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 );
289 const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
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 );
297 if (( N1 == 0 ) && ( fabs( beta ) > 0.0 )){
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 ){
308 * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSRdown + 1 ) )
313 double * block_right = Rregular->
gStorage( NR, TwoSR, IR, NR + 1, TwoSRdown, IRdown );
316 dgemm_( ¬rans, &trans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRup, &add, block_up, &dimLup );
324 if (( N1 == 1 ) && ( fabs( beta ) > 0.0 )){
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 );
332 * sqrt( 1.0 * ( TwoJ + 1 ) * ( TwoSRdown + 1 ) )
337 double * block_right = Rregular->
gStorage( NR, TwoSR, IR, NR + 1, TwoSRdown, IRdown );
340 dgemm_( ¬rans, &trans, &dimLup, &dimRup, &dimRdown, &factor, block_down, &dimLup, block_right, &dimRup, &add, block_up, &dimLup );
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 );
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 );
367 double inproduct = 0.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 );
375 dgemm_( ¬rans, &trans, &dimLup, &dimRup, &dimRdown, &one, block_down, &dimLup, block_right, &dimRup, &
set, workmem, &dimLup );
378 int size = dimLup * dimRup;
380 if ( fabs( gamma ) > 0.0 ){
381 double factor = gamma;
382 daxpy_( &size, &factor, workmem, &inc1, block_up, &inc1 );
384 inproduct = ddot_( &size, workmem, &inc1, block_up, &inc1 );
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 );
404 const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
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 );
412 if (( N2 == 0 ) && ( fabs( alpha ) > 0.0 )){
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
423 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoJdown + 1 ) )
428 double * block_left = Ltrans->
gStorage( NL - 1, TwoSLdown, ILdown, NL, TwoSL, IL );
431 dgemm_( &trans, ¬rans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &add, block_up, &dimLup );
439 if (( N2 == 1 ) && ( fabs( alpha ) > 0.0 )){
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
447 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoJ + 1 ) )
452 double * block_left = Ltrans->
gStorage( NL - 1, TwoSLdown, ILdown, NL, TwoSL, IL );
455 dgemm_( &trans, ¬rans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &add, block_up, &dimLup );
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 );
477 const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
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 );
485 if (( N2 == 2 ) && ( fabs( beta ) > 0.0 )){
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 ){
496 * sqrt( 1.0 * ( TwoJdown + 1 ) * ( TwoSLdown + 1 ) )
500 double * block_left = Lregular->
gStorage( NL, TwoSL, IL, NL + 1, TwoSLdown, ILdown );
503 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &add, block_up, &dimLup );
511 if (( N2 == 1 ) && ( fabs( beta ) > 0.0 )){
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 );
519 * sqrt( 1.0 * ( TwoSLdown + 1 ) * ( TwoJ + 1 ) )
523 double * block_left = Lregular->
gStorage( NL, TwoSL, IL, NL + 1, TwoSLdown, ILdown );
526 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &add, block_up, &dimLup );
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 );
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 );
553 double inproduct = 0.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 );
560 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimLdown, &one, block_left, &dimLup, block_down, &dimLdown, &
set, workmem, &dimLup );
563 int size = dimLup * dimRup;
565 if ( fabs( gamma ) > 0.0 ){
566 double factor = gamma;
567 daxpy_( &size, &factor, workmem, &inc1, block_up, &inc1 );
569 inproduct = ddot_( &size, workmem, &inc1, block_up, &inc1 );
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 );
590 const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
591 const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
593 int dimLup = book_up->
gCurrentDim( index, NL, TwoSL, IL );
594 int dimRup = book_up->
gCurrentDim( index + 2, NR, TwoSR, IR );
597 if ( fabs( alpha ) > 0.0 ){
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 ) )
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 );
617 dgemm_( &trans, ¬rans, &dimLup, &dimRdown, &dimLdown, &factor, block_left, &dimLdown, block_down, &dimLdown, &
set, workmem, &dimLup );
618 dgemm_( ¬rans, ¬rans, &dimLup, &dimRup, &dimRdown, &one, workmem, &dimLup, block_right, &dimRdown, &one, block_up, &dimLup );
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 );
642 const int TwoS1 = (( N1 == 1 ) ? 1 : 0 );
643 const int TwoS2 = (( N2 == 1 ) ? 1 : 0 );
645 int dimLup = book_up->
gCurrentDim( index, NL, TwoSL, IL );
646 int dimRup = book_up->
gCurrentDim( index + 2, NR, TwoSR, IR );
649 if ( fabs( beta ) > 0.0 ){
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 );
658 *
Special::phase( TwoSLdown + TwoSR + TwoJ + 1 + 2 * TwoS1 + 2 * TwoS2 )
659 * sqrt( 1.0 * ( TwoSLdown + 1 ) * ( TwoSRdown + 1 ) )
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 );
669 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &factor, block_left, &dimLup, block_down, &dimLdown, &
set, workmem, &dimLup );
670 dgemm_( ¬rans, &trans, &dimLup, &dimRup, &dimRdown, &one, workmem, &dimLup, block_right, &dimRup, &one, block_up, &dimLup );
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 );
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 );
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 );
706 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &one, block_left, &dimLup, block_down, &dimLdown, &
set, workmem1, &dimLup );
707 dgemm_( ¬rans, &trans, &dimLup, &dimRup, &dimRdown, &one, workmem1, &dimLup, block_right, &dimRup, &
set, workmem2, &dimLup );
710 int size = dimLup * dimRup;
712 if ( fabs( gamma ) > 0.0 ){
713 double factor = gamma;
714 daxpy_( &size, &factor, workmem2, &inc1, block_up, &inc1 );
716 inproduct = ddot_( &size, workmem2, &inc1, block_up, &inc1 );
int gTwoSL(const int ikappa) const
Get the left spin symmetry of block ikappa.
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)
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...
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.
void symm2prog()
Convert the storage from symmetric Hamiltonian convention to diagram convention.
static int phase(const int power_times_two)
Phase function.
int gTwoJ(const int ikappa) const
Get the central spin symmetry of block ikappa.
int gNKappa() const
Get the number of symmetry blocks.
double * gStorage()
Get the pointer to the storage.
int gIL(const int ikappa) const
Get the left irrep symmetry of block ikappa.
double * gStorage()
Get the pointer to the storage.
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.
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's convent...
void prog2symm()
Convert the storage from diagram convention to symmetric Hamiltonian convention.
int gNL(const int ikappa) const
Get the left particle number symmetry of block ikappa.
int gIR(const int ikappa) const
Get the right irrep symmetry of block ikappa.
int gTwoSR(const int ikappa) const
Get the right spin symmetry of block ikappa.
int gIndex() const
Get the location index.
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]) ...
int gNR(const int ikappa) const
Get the right particle number symmetry of block ikappa.
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.
int gN2(const int ikappa) const
Get the right local particle number symmetry of block ikappa.