30 #include "ConjugateGradient.h" 46 f_dot_4dm = contract_4dm;
49 struct timeval start, end;
50 gettimeofday( &start, NULL );
55 make_AA_CC(
true, 0.0 );
57 make_EE_GG(
true, 0.0 );
58 make_BB_FF_singlet(
true, 0.0 );
59 make_BB_FF_triplet(
true, 0.0 );
61 construct_rhs( oei, ints );
63 make_AA_CC(
false, IPEA );
64 make_DD(
false, IPEA );
65 make_EE_GG(
false, IPEA );
66 make_BB_FF_singlet(
false, IPEA );
67 make_BB_FF_triplet(
false, IPEA );
70 make_FAB_FCF_singlet();
71 make_FAB_FCF_triplet();
72 make_FBE_FFG_singlet();
73 make_FBE_FFG_triplet();
79 gettimeofday( &end, NULL );
80 double elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
81 cout <<
"CASPT2 : Wall time tensors = " << elapsed <<
" seconds" << endl;
82 gettimeofday( &start, NULL );
86 gettimeofday( &end, NULL );
87 elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
88 cout <<
"CASPT2 : Wall time diag(ovlp) = " << elapsed <<
" seconds" << endl;
94 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
95 delete [] FAA[ irrep ];
96 delete [] FCC[ irrep ];
97 delete [] FDD[ irrep ];
98 delete [] FEE[ irrep ];
99 delete [] FGG[ irrep ];
100 delete [] FBB_singlet[ irrep ];
101 delete [] FBB_triplet[ irrep ];
102 delete [] FFF_singlet[ irrep ];
103 delete [] FFF_triplet[ irrep ];
110 delete [] FBB_singlet;
111 delete [] FBB_triplet;
112 delete [] FFF_singlet;
113 delete [] FFF_triplet;
115 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
116 for (
int irrep_right = 0; irrep_right < num_irreps; irrep_right++ ){
118 const int num_w = indices->
getNDMRG( irrep_w );
119 for (
int w = 0; w < num_w; w++ ){
120 delete [] FAD[ irrep_left ][ irrep_right ][ w ];
121 delete [] FCD[ irrep_left ][ irrep_right ][ w ];
122 delete [] FAB_singlet[ irrep_left ][ irrep_right ][ w ];
123 delete [] FAB_triplet[ irrep_left ][ irrep_right ][ w ];
124 delete [] FCF_singlet[ irrep_left ][ irrep_right ][ w ];
125 delete [] FCF_triplet[ irrep_left ][ irrep_right ][ w ];
126 delete [] FBE_singlet[ irrep_left ][ irrep_right ][ w ];
127 delete [] FBE_triplet[ irrep_left ][ irrep_right ][ w ];
128 delete [] FFG_singlet[ irrep_left ][ irrep_right ][ w ];
129 delete [] FFG_triplet[ irrep_left ][ irrep_right ][ w ];
130 delete [] FDE_singlet[ irrep_left ][ irrep_right ][ w ];
131 delete [] FDE_triplet[ irrep_left ][ irrep_right ][ w ];
132 delete [] FDG_singlet[ irrep_left ][ irrep_right ][ w ];
133 delete [] FDG_triplet[ irrep_left ][ irrep_right ][ w ];
135 delete [] FAD[ irrep_left ][ irrep_right ];
136 delete [] FCD[ irrep_left ][ irrep_right ];
137 delete [] FAB_singlet[ irrep_left ][ irrep_right ];
138 delete [] FAB_triplet[ irrep_left ][ irrep_right ];
139 delete [] FCF_singlet[ irrep_left ][ irrep_right ];
140 delete [] FCF_triplet[ irrep_left ][ irrep_right ];
141 delete [] FBE_singlet[ irrep_left ][ irrep_right ];
142 delete [] FBE_triplet[ irrep_left ][ irrep_right ];
143 delete [] FFG_singlet[ irrep_left ][ irrep_right ];
144 delete [] FFG_triplet[ irrep_left ][ irrep_right ];
145 delete [] FDE_singlet[ irrep_left ][ irrep_right ];
146 delete [] FDE_triplet[ irrep_left ][ irrep_right ];
147 delete [] FDG_singlet[ irrep_left ][ irrep_right ];
148 delete [] FDG_triplet[ irrep_left ][ irrep_right ];
150 delete [] FAD[ irrep_left ];
151 delete [] FCD[ irrep_left ];
152 delete [] FAB_singlet[ irrep_left ];
153 delete [] FAB_triplet[ irrep_left ];
154 delete [] FCF_singlet[ irrep_left ];
155 delete [] FCF_triplet[ irrep_left ];
156 delete [] FBE_singlet[ irrep_left ];
157 delete [] FBE_triplet[ irrep_left ];
158 delete [] FFG_singlet[ irrep_left ];
159 delete [] FFG_triplet[ irrep_left ];
160 delete [] FDE_singlet[ irrep_left ];
161 delete [] FDE_triplet[ irrep_left ];
162 delete [] FDG_singlet[ irrep_left ];
163 delete [] FDG_triplet[ irrep_left ];
167 delete [] FAB_singlet;
168 delete [] FAB_triplet;
169 delete [] FCF_singlet;
170 delete [] FCF_triplet;
171 delete [] FBE_singlet;
172 delete [] FBE_triplet;
173 delete [] FFG_singlet;
174 delete [] FFG_triplet;
175 delete [] FDE_singlet;
176 delete [] FDE_triplet;
177 delete [] FDG_singlet;
178 delete [] FDG_triplet;
180 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
181 const int num_w = indices->
getNDMRG( irrep );
182 for (
int w = 0; w < num_w; w++ ){
183 delete [] FEH[ irrep ][ w ];
184 delete [] FGH[ irrep ][ w ];
186 delete [] FEH[ irrep ];
187 delete [] FGH[ irrep ];
197 delete [] size_B_singlet;
198 delete [] size_B_triplet;
199 delete [] size_F_singlet;
200 delete [] size_F_triplet;
202 delete [] vector_rhs;
208 struct timeval start, end;
209 gettimeofday( &start, NULL );
212 const int normalizations[] = { 1, 2, 2, 1, 1, 2, 6, 2, 2, 2, 6, 4, 12 };
213 const bool apply_shift = (( fabs( imag_shift ) > 0.0 ) ?
true :
false );
215 int total_size = jump[ CHEMPS2_CASPT2_NUM_CASES * num_irreps ];
216 double * diag_fock =
new double[ total_size ];
217 diagonal( diag_fock );
218 double min_eig = diag_fock[ 0 ];
219 for (
int elem = 1; elem < total_size; elem++ ){ min_eig = min( min_eig, diag_fock[ elem ] ); }
220 cout <<
"CASPT2 : Solution algorithm = " << (( CONJUGATE_GRADIENT ) ?
"Conjugate Gradient" :
"Davidson" ) << endl;
221 cout <<
"CASPT2 : Minimum(diagonal) = " << min_eig << endl;
223 ConjugateGradient * CG = (( CONJUGATE_GRADIENT ) ?
new ConjugateGradient( total_size, CheMPS2::CONJ_GRADIENT_RTOL, CheMPS2::CONJ_GRADIENT_PRECOND_CUTOFF,
false ) : NULL );
224 Davidson * DAVID = (( CONJUGATE_GRADIENT ) ? NULL :
new Davidson( total_size,
225 CheMPS2::DAVIDSON_NUM_VEC,
226 CheMPS2::DAVIDSON_NUM_VEC_KEEP,
227 CheMPS2::CONJ_GRADIENT_RTOL,
228 CheMPS2::CONJ_GRADIENT_PRECOND_CUTOFF,
231 double ** pointers =
new double*[ 3 ];
232 char instruction = (( CONJUGATE_GRADIENT ) ? CG->
step( pointers ) : DAVID->
FetchInstruction( pointers ));
233 assert( instruction ==
'A' );
234 for (
int elem = 0; elem < total_size; elem++ ){ pointers[ 0 ][ elem ] = vector_rhs[ elem ] / diag_fock[ elem ]; }
235 for (
int elem = 0; elem < total_size; elem++ ){ pointers[ 1 ][ elem ] = diag_fock[ elem ]; }
236 for (
int elem = 0; elem < total_size; elem++ ){ pointers[ 2 ][ elem ] = vector_rhs[ elem ]; }
238 const double E2_DIAGONAL = - ddot_( &total_size, pointers[ 0 ], &inc1, pointers[ 2 ], &inc1 );
239 instruction = (( CONJUGATE_GRADIENT ) ? CG->
step( pointers ) : DAVID->
FetchInstruction( pointers ));
240 assert( instruction ==
'B' );
241 while ( instruction ==
'B' ){
242 matvec( pointers[ 0 ], pointers[ 1 ], diag_fock );
243 if ( apply_shift ){ add_shift( pointers[ 0 ], pointers[ 1 ], diag_fock, imag_shift, normalizations ); }
244 instruction = (( CONJUGATE_GRADIENT ) ? CG->
step( pointers ) : DAVID->
FetchInstruction( pointers ));
246 assert( instruction ==
'C' );
247 const double E2_NONVARIATIONAL = - ddot_( &total_size, pointers[ 0 ], &inc1, vector_rhs, &inc1 );
248 const double rnorm = pointers[ 1 ][ 0 ];
250 cout <<
"CASPT2 : Residual norm = " << rnorm << endl;
251 matvec( pointers[ 0 ], pointers[ 1 ], diag_fock );
252 const double E2_VARIATIONAL = 2 * E2_NONVARIATIONAL + ddot_( &total_size, pointers[ 0 ], &inc1, pointers[ 1 ], &inc1 );
255 const double inproduct = inproduct_vectors( pointers[ 0 ], pointers[ 0 ], normalizations );
256 const double reference_weight = 1.0 / ( 1.0 + inproduct );
257 cout <<
"CASPT2 : Reference weight = " << reference_weight << endl;
260 if ( CG != NULL ){
delete CG; }
261 if ( DAVID != NULL ){
delete DAVID; }
263 gettimeofday( &end, NULL );
264 double elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
265 cout <<
"CASPT2 : Wall time solution = " << elapsed <<
" seconds" << endl;
267 cout <<
"CASPT2 : E2 [DIAGONAL] = " << E2_DIAGONAL << endl;
268 cout <<
"CASPT2 : E2 [NON-VARIATIONAL] = " << E2_NONVARIATIONAL << endl;
269 cout <<
"CASPT2 : E2 [VARIATIONAL] = " << E2_VARIATIONAL << endl;
270 return E2_VARIATIONAL;
274 void CheMPS2::CASPT2::add_shift(
double * vector,
double * result,
double * diag_fock,
const double imag_shift,
const int * normalizations )
const{
276 for (
int sector = 0; sector < CHEMPS2_CASPT2_NUM_CASES; sector++ ){
277 const int start = jump[ num_irreps * sector ];
278 const int stop = jump[ num_irreps * ( sector + 1 ) ];
279 const double factor = imag_shift * imag_shift * normalizations[ sector ] * normalizations[ sector ];
280 for (
int elem = start; elem < stop; elem ++ ){
281 result[ elem ] += factor * vector[ elem ] / diag_fock[ elem ];
287 double CheMPS2::CASPT2::inproduct_vectors(
double * first,
double * second,
const int * normalizations )
const{
291 for (
int sector = 0; sector < CHEMPS2_CASPT2_NUM_CASES; sector++ ){
292 int pointer = jump[ num_irreps * sector ];
293 int size = jump[ num_irreps * ( sector + 1 ) ] - pointer;
294 value += normalizations[ sector ] * ddot_( &size, first + pointer, &inc1, second + pointer, &inc1 );
300 void CheMPS2::CASPT2::energy_per_sector(
double * solution )
const{
303 double energies[ CHEMPS2_CASPT2_NUM_CASES ];
304 for (
int sector = 0; sector < CHEMPS2_CASPT2_NUM_CASES; sector++ ){
305 int pointer = jump[ num_irreps * sector ];
306 int size = jump[ num_irreps * ( sector + 1 ) ] - pointer;
307 energies[ sector ] = - ddot_( &size, solution + pointer, &inc1, vector_rhs + pointer, &inc1 );
309 cout <<
"************************************************" << endl;
310 cout <<
"* CASPT2 non-variational energy per sector *" << endl;
311 cout <<
"************************************************" << endl;
312 cout <<
" A or VJTU = " << energies[ CHEMPS2_CASPT2_A ] << endl;
313 cout <<
" B or VJTI = " << energies[ CHEMPS2_CASPT2_B_SINGLET ] + energies[ CHEMPS2_CASPT2_B_TRIPLET ] << endl;
314 cout <<
" C or ATVX = " << energies[ CHEMPS2_CASPT2_C ] << endl;
315 cout <<
" D or AIVX = " << energies[ CHEMPS2_CASPT2_D ] << endl;
316 cout <<
" E or VJAI = " << energies[ CHEMPS2_CASPT2_E_SINGLET ] + energies[ CHEMPS2_CASPT2_E_TRIPLET ] << endl;
317 cout <<
" F or BVAT = " << energies[ CHEMPS2_CASPT2_F_SINGLET ] + energies[ CHEMPS2_CASPT2_F_TRIPLET ] << endl;
318 cout <<
" G or BJAT = " << energies[ CHEMPS2_CASPT2_G_SINGLET ] + energies[ CHEMPS2_CASPT2_G_TRIPLET ] << endl;
319 cout <<
" H or BJAI = " << energies[ CHEMPS2_CASPT2_H_SINGLET ] + energies[ CHEMPS2_CASPT2_H_TRIPLET ] << endl;
320 cout <<
"************************************************" << endl;
324 void CheMPS2::CASPT2::create_f_dots(){
327 f_dot_3dm =
new double[ LAS * LAS * LAS * LAS ];
328 f_dot_2dm =
new double[ LAS * LAS ];
331 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
332 const int NDMRG = indices->
getNDMRG( irrep );
333 const int NOCC = indices->
getNOCC( irrep );
337 for (
int orb = 0; orb < NDMRG; orb++ ){
338 value += fock->
get( irrep, NOCC + orb, NOCC + orb ) * one_rdm[ jumpx + orb + LAS * ( jumpx + orb ) ];
343 for (
int irrep1 = 0; irrep1 < num_irreps; irrep1++ ){
345 const int upper1 = lower1 + indices->
getNDMRG( irrep1 );
347 for (
int i1 = lower1; i1 < upper1; i1++ ){
348 for (
int i2 = lower1; i2 < upper1; i2++ ){
351 for (
int irrepx = 0; irrepx < num_irreps; irrepx++ ){
352 const int NOCCx = indices->
getNOCC( irrepx );
353 const int NDMRGx = indices->
getNDMRG( irrepx );
356 for (
int orbx = 0; orbx < NDMRGx; orbx++ ){
357 value += fock->
get( irrepx, NOCCx + orbx, NOCCx + orbx ) * two_rdm[ i1 + LAS * ( jumpx + orbx + LAS * ( i2 + LAS * ( jumpx + orbx ))) ];
360 f_dot_2dm[ i1 + LAS * i2 ] = value;
365 for (
int irrep1 = 0; irrep1 < num_irreps; irrep1++ ){
367 const int upper1 = lower1 + indices->
getNDMRG( irrep1 );
368 for (
int irrep2 = 0; irrep2 < num_irreps; irrep2++ ){
370 const int upper2 = lower2 + indices->
getNDMRG( irrep2 );
372 for (
int irrep3 = 0; irrep3 < num_irreps; irrep3++ ){
376 const int upper3 = lower3 + indices->
getNDMRG( irrep3 );
377 const int upper4 = lower4 + indices->
getNDMRG( irrep4 );
379 for (
int i1 = lower1; i1 < upper1; i1++ ){
380 for (
int i2 = lower2; i2 < upper2; i2++ ){
381 for (
int i3 = lower3; i3 < upper3; i3++ ){
382 for (
int i4 = lower4; i4 < upper4; i4++ ){
385 for (
int irrepx = 0; irrepx < num_irreps; irrepx++ ){
387 const int NOCCx = indices->
getNOCC( irrepx );
388 const int NDMRGx = indices->
getNDMRG( irrepx );
390 for (
int orbx = 0; orbx < NDMRGx; orbx++ ){
391 value += ( fock->
get( irrepx, NOCCx + orbx, NOCCx + orbx )
392 * three_rdm[ i1 + LAS*( i2 + LAS*( jumpx + orbx + LAS*( i3 + LAS*( i4 + LAS*( jumpx + orbx ))))) ] );
395 f_dot_3dm[ i1 + LAS * ( i2 + LAS * ( i3 + LAS * i4 )) ] = value;
419 int CheMPS2::CASPT2::vector_helper(){
421 int * helper =
new int[ CHEMPS2_CASPT2_NUM_CASES * num_irreps ];
438 size_A =
new int[ num_irreps ];
439 size_C =
new int[ num_irreps ];
440 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
442 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
443 for (
int irrep_u = 0; irrep_u < num_irreps; irrep_u++ ){
448 size_A[ irrep ] = linsize_AC;
449 size_C[ irrep ] = linsize_AC;
450 helper[ irrep + num_irreps * CHEMPS2_CASPT2_A ] = size_A[ irrep ] * indices->
getNOCC( irrep );
451 helper[ irrep + num_irreps * CHEMPS2_CASPT2_C ] = size_C[ irrep ] * indices->
getNVIRT( irrep );
474 size_D =
new int[ num_irreps ];
475 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
477 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
479 const int nact_t = indices->
getNDMRG( irrep_t );
480 const int nact_u = indices->
getNDMRG( irrep_u );
481 jump_tu += nact_t * nact_u;
483 size_D[ irrep ] = 2 * jump_tu;
485 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
487 const int nocc_i = indices->
getNOCC( irrep_i );
488 const int nvirt_a = indices->
getNVIRT( irrep_a );
489 jump_ai += nocc_i * nvirt_a;
491 helper[ irrep + num_irreps * CHEMPS2_CASPT2_D ] = jump_ai * size_D[ irrep ];
546 size_B_singlet =
new int[ num_irreps ];
547 size_B_triplet =
new int[ num_irreps ];
548 size_F_singlet =
new int[ num_irreps ];
549 size_F_triplet =
new int[ num_irreps ];
550 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
551 int jump_tu_singlet = 0;
552 int jump_tu_triplet = 0;
554 for (
int irrep_tu = 0; irrep_tu < num_irreps; irrep_tu++ ){
555 const int nact_tu = indices->
getNDMRG( irrep_tu );
556 jump_tu_singlet += ( nact_tu * ( nact_tu + 1 )) / 2;
557 jump_tu_triplet += ( nact_tu * ( nact_tu - 1 )) / 2;
560 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
562 if ( irrep_t < irrep_u ){
563 const int nact_t = indices->
getNDMRG( irrep_t );
564 const int nact_u = indices->
getNDMRG( irrep_u );
565 jump_tu_singlet += nact_t * nact_u;
566 jump_tu_triplet += nact_t * nact_u;
570 size_B_singlet[ irrep ] = jump_tu_singlet;
571 size_B_triplet[ irrep ] = jump_tu_triplet;
572 size_F_singlet[ irrep ] = jump_tu_singlet;
573 size_F_triplet[ irrep ] = jump_tu_triplet;
575 int linsize_B_singlet = 0;
576 int linsize_B_triplet = 0;
577 int linsize_F_singlet = 0;
578 int linsize_F_triplet = 0;
580 for (
int irrep_ijab = 0; irrep_ijab < num_irreps; irrep_ijab++ ){
581 const int nocc_ij = indices->
getNOCC( irrep_ijab );
582 linsize_B_singlet += ( nocc_ij * ( nocc_ij + 1 ))/2;
583 linsize_B_triplet += ( nocc_ij * ( nocc_ij - 1 ))/2;
584 const int nvirt_ab = indices->
getNVIRT( irrep_ijab );
585 linsize_F_singlet += ( nvirt_ab * ( nvirt_ab + 1 ))/2;
586 linsize_F_triplet += ( nvirt_ab * ( nvirt_ab - 1 ))/2;
589 for (
int irrep_ai = 0; irrep_ai < num_irreps; irrep_ai++ ){
591 if ( irrep_ai < irrep_bj ){
592 const int nocc_i = indices->
getNOCC( irrep_ai );
593 const int nocc_j = indices->
getNOCC( irrep_bj );
594 linsize_B_singlet += nocc_i * nocc_j;
595 linsize_B_triplet += nocc_i * nocc_j;
596 const int nvirt_a = indices->
getNVIRT( irrep_ai );
597 const int nvirt_b = indices->
getNVIRT( irrep_bj );
598 linsize_F_singlet += nvirt_a * nvirt_b;
599 linsize_F_triplet += nvirt_a * nvirt_b;
603 helper[ irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] = linsize_B_singlet * size_B_singlet[ irrep ];
604 helper[ irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] = linsize_B_triplet * size_B_triplet[ irrep ];
605 helper[ irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] = linsize_F_singlet * size_F_singlet[ irrep ];
606 helper[ irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] = linsize_F_triplet * size_F_triplet[ irrep ];
632 size_E =
new int[ num_irreps ];
633 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
634 size_E[ irrep ] = indices->
getNDMRG( irrep );
635 int linsize_E_singlet = 0;
636 int linsize_E_triplet = 0;
637 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
638 const int nvirt_a = indices->
getNVIRT( irrep_a );
640 if ( irrep_occ == 0 ){
641 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
642 const int nocc_ij = indices->
getNOCC( irrep_ij );
643 linsize_E_singlet += ( nvirt_a * nocc_ij * ( nocc_ij + 1 )) / 2;
644 linsize_E_triplet += ( nvirt_a * nocc_ij * ( nocc_ij - 1 )) / 2;
647 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
649 if ( irrep_i < irrep_j ){
650 const int nocc_i = indices->
getNOCC( irrep_i );
651 const int nocc_j = indices->
getNOCC( irrep_j );
652 linsize_E_singlet += nvirt_a * nocc_i * nocc_j;
653 linsize_E_triplet += nvirt_a * nocc_i * nocc_j;
658 helper[ irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] = linsize_E_singlet * size_E[ irrep ];
659 helper[ irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] = linsize_E_triplet * size_E[ irrep ];
685 size_G =
new int[ num_irreps ];
686 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
687 size_G[ irrep ] = indices->
getNDMRG( irrep );
688 int linsize_G_singlet = 0;
689 int linsize_G_triplet = 0;
690 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
691 const int nocc_i = indices->
getNOCC( irrep_i );
693 if ( irrep_virt == 0 ){
694 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
695 const int nvirt_ab = indices->
getNVIRT( irrep_ab );
696 linsize_G_singlet += ( nocc_i * nvirt_ab * ( nvirt_ab + 1 )) / 2;
697 linsize_G_triplet += ( nocc_i * nvirt_ab * ( nvirt_ab - 1 )) / 2;
700 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
702 if ( irrep_a < irrep_b ){
703 const int nvirt_a = indices->
getNVIRT( irrep_a );
704 const int nvirt_b = indices->
getNVIRT( irrep_b );
705 linsize_G_singlet += nocc_i * nvirt_a * nvirt_b;
706 linsize_G_triplet += nocc_i * nvirt_a * nvirt_b;
711 helper[ irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] = linsize_G_singlet * size_G[ irrep ];
712 helper[ irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] = linsize_G_triplet * size_G[ irrep ];
744 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
745 int linsize_H_singlet = 0;
746 int linsize_H_triplet = 0;
748 int linsize_ij_singlet = 0;
749 int linsize_ij_triplet = 0;
750 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
751 const int nocc_ij = indices->
getNOCC( irrep_ij );
752 linsize_ij_singlet += ( nocc_ij * ( nocc_ij + 1 )) / 2;
753 linsize_ij_triplet += ( nocc_ij * ( nocc_ij - 1 )) / 2;
755 int linsize_ab_singlet = 0;
756 int linsize_ab_triplet = 0;
757 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
758 const int nvirt_ab = indices->
getNVIRT( irrep_ab );
759 linsize_ab_singlet += ( nvirt_ab * ( nvirt_ab + 1 )) / 2;
760 linsize_ab_triplet += ( nvirt_ab * ( nvirt_ab - 1 )) / 2;
762 linsize_H_singlet = linsize_ij_singlet * linsize_ab_singlet;
763 linsize_H_triplet = linsize_ij_triplet * linsize_ab_triplet;
766 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
768 if ( irrep_i < irrep_j ){ linsize_ij += indices->
getNOCC( irrep_i ) * indices->
getNOCC( irrep_j ); }
771 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
773 if ( irrep_a < irrep_b ){ linsize_ab += indices->
getNVIRT( irrep_a ) * indices->
getNVIRT( irrep_b ); }
775 linsize_H_singlet = linsize_ij * linsize_ab;
776 linsize_H_triplet = linsize_ij * linsize_ab;
778 helper[ irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] = linsize_H_singlet;
779 helper[ irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] = linsize_H_triplet;
782 jump =
new int[ CHEMPS2_CASPT2_NUM_CASES * num_irreps + 1 ];
784 for (
int cnt = 0; cnt < CHEMPS2_CASPT2_NUM_CASES * num_irreps; cnt++ ){ jump[ cnt+1 ] = jump[ cnt ] + helper[ cnt ]; }
786 const int total_size = jump[ CHEMPS2_CASPT2_NUM_CASES * num_irreps ];
788 cout <<
"CASPT2 : Old size V_SD space = " << total_size << endl;
796 long long length = 0;
797 for (
int i1 = 0; i1 < NUM_IRREPS; i1++ ){
798 const long long nocc1 = idx->
getNOCC( i1 );
799 const long long nact1 = idx->
getNDMRG( i1 );
800 const long long nvir1 = idx->
getNVIRT( i1 );
801 for (
int i2 = 0; i2 < NUM_IRREPS; i2++ ){
802 const long long nocc2 = idx->
getNOCC( i2 );
803 const long long nact2 = idx->
getNDMRG( i2 );
804 for (
int i3 = 0; i3 < NUM_IRREPS; i3++ ){
806 const long long nact3 = idx->
getNDMRG( i3 );
807 const long long nvir3 = idx->
getNVIRT( i3 );
808 const long long nocc4 = idx->
getNOCC( i4 );
809 const long long nact4 = idx->
getNDMRG( i4 );
810 const long long nvir4 = idx->
getNVIRT( i4 );
812 length += nocc1 * nact2 * nact3 * nact4;
813 length += nact1 * nact2 * nact3 * nvir4;
814 length += 2 * nocc1 * nact2 * nact3 * nvir4;
815 length += nocc1 * nocc2 * nact3 * nvir4;
816 length += nocc1 * nact2 * nvir3 * nvir4;
818 length += nact1 * nact3 * nocc2 * nocc4;
819 length += nvir1 * nvir3 * nact2 * nact4;
820 length += nvir1 * nvir3 * nocc2 * nocc4;
824 length += ( nact1 * nact3 * nocc2 * ( nocc2 - 1 ) ) / 2;
825 length += ( nvir1 * nvir3 * nact2 * ( nact2 - 1 ) ) / 2;
826 length += ( nvir1 * nvir3 * nocc2 * ( nocc2 - 1 ) ) / 2;
828 length += ( nact1 * ( nact3 + 1 ) * nocc2 ) / 2;
829 length += ( nvir1 * ( nvir3 + 1 ) * nact2 ) / 2;
830 length += ( nvir1 * ( nvir3 + 1 ) * nocc2 ) / 2;
840 int CheMPS2::CASPT2::recreatehelper1(
double * FOCK,
double * OVLP,
int SIZE,
double * work,
double * eigs,
int lwork ){
842 if ( SIZE == 0 ){
return SIZE; }
848 dsyev_( &jobz, &uplo, &SIZE, OVLP, &SIZE, eigs, work, &lwork, &info );
853 while (( skip < SIZE ) && ( ctu )){
854 if ( eigs[ skip ] < CheMPS2::CASPT2_OVLP_CUTOFF ){
860 int NEWSIZE = SIZE - skip;
862 if ( NEWSIZE == 0 ){
return NEWSIZE; }
865 for (
int col = skip; col < SIZE; col++ ){
866 const double prefactor = 1.0 / sqrt( eigs[ col ] );
867 for (
int row = 0; row < SIZE; row++ ){
868 OVLP[ row + SIZE * col ] *= prefactor;
877 dgemm_( ¬rans, ¬rans, &SIZE, &NEWSIZE, &SIZE, &one, FOCK, &SIZE, OVLP + skip * SIZE, &SIZE, &
set, work, &SIZE );
878 dgemm_( &trans, ¬rans, &NEWSIZE, &NEWSIZE, &SIZE, &one, OVLP + skip * SIZE, &SIZE, work, &SIZE, &
set, FOCK, &NEWSIZE );
881 dsyev_( &jobz, &uplo, &NEWSIZE, FOCK, &NEWSIZE, eigs, work, &lwork, &info );
884 dgemm_( ¬rans, ¬rans, &SIZE, &NEWSIZE, &NEWSIZE, &one, OVLP + skip * SIZE, &SIZE, FOCK, &NEWSIZE, &
set, work, &SIZE );
885 int size_copy = SIZE * NEWSIZE;
887 dcopy_( &size_copy, work, &inc1, OVLP, &inc1 );
890 dcopy_( &NEWSIZE, eigs, &inc1, FOCK, &inc1 );
896 void CheMPS2::CASPT2::recreatehelper2(
double * LEFT,
double * RIGHT,
double ** matrix,
double * work,
int OLD_LEFT,
int NEW_LEFT,
int OLD_RIGHT,
int NEW_RIGHT,
const int number ){
902 for (
int count = 0; count < number; count++ ){
904 dgemm_( &trans, ¬rans, &NEW_LEFT, &OLD_RIGHT, &OLD_LEFT, &one, LEFT, &OLD_LEFT, matrix[ count ], &OLD_LEFT, &
set, work, &NEW_LEFT );
906 dgemm_( ¬rans, ¬rans, &NEW_LEFT, &NEW_RIGHT, &OLD_RIGHT, &one, work, &NEW_LEFT, RIGHT, &OLD_RIGHT, &
set, matrix[ count ], &NEW_LEFT );
911 void CheMPS2::CASPT2::recreatehelper3(
double * OVLP,
int OLDSIZE,
int NEWSIZE,
double * rhs_old,
double * rhs_new,
const int num_rhs ){
918 for (
int sector = 0; sector < num_rhs; sector++ ){
919 dgemv_( &trans, &OLDSIZE, &NEWSIZE, &one, OVLP, &OLDSIZE, rhs_old + OLDSIZE * sector, &inc1, &
set, rhs_new + NEWSIZE * sector, &inc1 );
924 void CheMPS2::CASPT2::recreate(){
926 const int maxsize = get_maxsize();
927 const int lwork = maxsize * maxsize;
928 double * work =
new double[ lwork ];
929 double * eigs =
new double[ maxsize ];
931 int * newsize_A =
new int[ num_irreps ];
932 int * newsize_C =
new int[ num_irreps ];
933 int * newsize_D =
new int[ num_irreps ];
934 int * newsize_E =
new int[ num_irreps ];
935 int * newsize_G =
new int[ num_irreps ];
936 int * newsize_B_singlet =
new int[ num_irreps ];
937 int * newsize_B_triplet =
new int[ num_irreps ];
938 int * newsize_F_singlet =
new int[ num_irreps ];
939 int * newsize_F_triplet =
new int[ num_irreps ];
941 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
942 newsize_A[ irrep ] = recreatehelper1( FAA[ irrep ], SAA[ irrep ], size_A[ irrep ], work, eigs, lwork );
943 newsize_C[ irrep ] = recreatehelper1( FCC[ irrep ], SCC[ irrep ], size_C[ irrep ], work, eigs, lwork );
944 newsize_D[ irrep ] = recreatehelper1( FDD[ irrep ], SDD[ irrep ], size_D[ irrep ], work, eigs, lwork );
945 newsize_E[ irrep ] = recreatehelper1( FEE[ irrep ], SEE[ irrep ], size_E[ irrep ], work, eigs, lwork );
946 newsize_G[ irrep ] = recreatehelper1( FGG[ irrep ], SGG[ irrep ], size_G[ irrep ], work, eigs, lwork );
947 newsize_B_singlet[ irrep ] = recreatehelper1( FBB_singlet[ irrep ], SBB_singlet[ irrep ], size_B_singlet[ irrep ], work, eigs, lwork );
948 newsize_B_triplet[ irrep ] = recreatehelper1( FBB_triplet[ irrep ], SBB_triplet[ irrep ], size_B_triplet[ irrep ], work, eigs, lwork );
949 newsize_F_singlet[ irrep ] = recreatehelper1( FFF_singlet[ irrep ], SFF_singlet[ irrep ], size_F_singlet[ irrep ], work, eigs, lwork );
950 newsize_F_triplet[ irrep ] = recreatehelper1( FFF_triplet[ irrep ], SFF_triplet[ irrep ], size_F_triplet[ irrep ], work, eigs, lwork );
953 for (
int IL = 0; IL < num_irreps; IL++ ){
954 for (
int IR = 0; IR < num_irreps; IR++ ){
957 const int num_w = indices->
getNDMRG( irrep_w );
959 if ( newsize_A[ IL ] * newsize_D[ IR ] * num_w > 0 ){
960 recreatehelper2( SAA[ IL ], SDD[ IR ], FAD[ IL ][ IR ], work, size_A[ IL ], newsize_A[ IL ], size_D[ IR ], newsize_D[ IR ], num_w );
963 if ( newsize_C[ IL ] * newsize_D[ IR ] * num_w > 0 ){
964 recreatehelper2( SCC[ IL ], SDD[ IR ], FCD[ IL ][ IR ], work, size_C[ IL ], newsize_C[ IL ], size_D[ IR ], newsize_D[ IR ], num_w );
967 if ( newsize_A[ IL ] * newsize_B_singlet[ IR ] * num_w > 0 ){
968 recreatehelper2( SAA[ IL ], SBB_singlet[ IR ], FAB_singlet[ IL ][ IR ], work, size_A[ IL ], newsize_A[ IL ], size_B_singlet[ IR ], newsize_B_singlet[ IR ], num_w );
971 if ( newsize_A[ IL ] * newsize_B_triplet[ IR ] * num_w > 0 ){
972 recreatehelper2( SAA[ IL ], SBB_triplet[ IR ], FAB_triplet[ IL ][ IR ], work, size_A[ IL ], newsize_A[ IL ], size_B_triplet[ IR ], newsize_B_triplet[ IR ], num_w );
975 if ( newsize_C[ IL ] * newsize_F_singlet[ IR ] * num_w > 0 ){
976 recreatehelper2( SCC[ IL ], SFF_singlet[ IR ], FCF_singlet[ IL ][ IR ], work, size_C[ IL ], newsize_C[ IL ], size_F_singlet[ IR ], newsize_F_singlet[ IR ], num_w );
979 if ( newsize_C[ IL ] * newsize_F_triplet[ IR ] * num_w > 0 ){
980 recreatehelper2( SCC[ IL ], SFF_triplet[ IR ], FCF_triplet[ IL ][ IR ], work, size_C[ IL ], newsize_C[ IL ], size_F_triplet[ IR ], newsize_F_triplet[ IR ], num_w );
983 if ( newsize_B_singlet[ IL ] * newsize_E[ IR ] * num_w > 0 ){
984 recreatehelper2( SBB_singlet[ IL ], SEE[ IR ], FBE_singlet[ IL ][ IR ], work, size_B_singlet[ IL ], newsize_B_singlet[ IL ], size_E[ IR ], newsize_E[ IR ], num_w );
987 if ( newsize_B_triplet[ IL ] * newsize_E[ IR ] * num_w > 0 ){
988 recreatehelper2( SBB_triplet[ IL ], SEE[ IR ], FBE_triplet[ IL ][ IR ], work, size_B_triplet[ IL ], newsize_B_triplet[ IL ], size_E[ IR ], newsize_E[ IR ], num_w );
991 if ( newsize_F_singlet[ IL ] * newsize_G[ IR ] * num_w > 0 ){
992 recreatehelper2( SFF_singlet[ IL ], SGG[ IR ], FFG_singlet[ IL ][ IR ], work, size_F_singlet[ IL ], newsize_F_singlet[ IL ], size_G[ IR ], newsize_G[ IR ], num_w );
995 if ( newsize_F_triplet[ IL ] * newsize_G[ IR ] * num_w > 0 ){
996 recreatehelper2( SFF_triplet[ IL ], SGG[ IR ], FFG_triplet[ IL ][ IR ], work, size_F_triplet[ IL ], newsize_F_triplet[ IL ], size_G[ IR ], newsize_G[ IR ], num_w );
999 if ( newsize_D[ IL ] * newsize_E[ IR ] * num_w > 0 ){
1000 recreatehelper2( SDD[ IL ], SEE[ IR ], FDE_singlet[ IL ][ IR ], work, size_D[ IL ], newsize_D[ IL ], size_E[ IR ], newsize_E[ IR ], num_w );
1001 recreatehelper2( SDD[ IL ], SEE[ IR ], FDE_triplet[ IL ][ IR ], work, size_D[ IL ], newsize_D[ IL ], size_E[ IR ], newsize_E[ IR ], num_w );
1004 if ( newsize_D[ IL ] * newsize_G[ IR ] * num_w > 0 ){
1005 recreatehelper2( SDD[ IL ], SGG[ IR ], FDG_singlet[ IL ][ IR ], work, size_D[ IL ], newsize_D[ IL ], size_G[ IR ], newsize_G[ IR ], num_w );
1006 recreatehelper2( SDD[ IL ], SGG[ IR ], FDG_triplet[ IL ][ IR ], work, size_D[ IL ], newsize_D[ IL ], size_G[ IR ], newsize_G[ IR ], num_w );
1010 if ( newsize_E[ IL ] * num_w > 0 ){
1012 recreatehelper2( SEE[ IL ], &one, FEH[ IL ], work, size_E[ IL ], newsize_E[ IL ], 1, 1, num_w );
1015 if ( newsize_G[ IL ] * num_w > 0 ){
1017 recreatehelper2( SGG[ IL ], &one, FGH[ IL ], work, size_G[ IL ], newsize_G[ IL ], 1, 1, num_w );
1026 double * tempvector_rhs =
new double[ jump[ num_irreps * CHEMPS2_CASPT2_NUM_CASES ] ];
1027 int * helper =
new int[ num_irreps * CHEMPS2_CASPT2_NUM_CASES ];
1028 for (
int ptr = 0; ptr < num_irreps * CHEMPS2_CASPT2_NUM_CASES; ptr++ ){ helper[ ptr ] = 0; }
1030 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
1032 if ( newsize_A[ irrep ] > 0 ){
1033 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_A;
1034 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_A[ irrep ];
1035 assert( num_rhs * size_A[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1036 helper[ ptr ] = num_rhs * newsize_A[ irrep ];
1037 recreatehelper3( SAA[ irrep ], size_A[ irrep ], newsize_A[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1040 if ( newsize_B_singlet[ irrep ] > 0 ){
1041 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET;
1042 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_B_singlet[ irrep ];
1043 assert( num_rhs * size_B_singlet[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1044 helper[ ptr ] = num_rhs * newsize_B_singlet[ irrep ];
1045 recreatehelper3( SBB_singlet[ irrep ], size_B_singlet[ irrep ], newsize_B_singlet[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1048 if ( newsize_B_triplet[ irrep ] > 0 ){
1049 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET;
1050 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_B_triplet[ irrep ];
1051 assert( num_rhs * size_B_triplet[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1052 helper[ ptr ] = num_rhs * newsize_B_triplet[ irrep ];
1053 recreatehelper3( SBB_triplet[ irrep ], size_B_triplet[ irrep ], newsize_B_triplet[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1056 if ( newsize_C[ irrep ] > 0 ){
1057 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_C;
1058 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_C[ irrep ];
1059 assert( num_rhs * size_C[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1060 helper[ ptr ] = num_rhs * newsize_C[ irrep ];
1061 recreatehelper3( SCC[ irrep ], size_C[ irrep ], newsize_C[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1064 if ( newsize_D[ irrep ] > 0 ){
1065 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_D;
1066 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_D[ irrep ];
1067 assert( num_rhs * size_D[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1068 helper[ ptr ] = num_rhs * newsize_D[ irrep ];
1069 recreatehelper3( SDD[ irrep ], size_D[ irrep ], newsize_D[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1072 if ( newsize_E[ irrep ] > 0 ){
1073 const int ptr1 = irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET;
1074 const int num_rhs1 = ( jump[ ptr1 + 1 ] - jump[ ptr1 ] ) / size_E[ irrep ];
1075 assert( num_rhs1 * size_E[ irrep ] == jump[ ptr1 + 1 ] - jump[ ptr1 ] );
1076 helper[ ptr1 ] = num_rhs1 * newsize_E[ irrep ];
1077 recreatehelper3( SEE[ irrep ], size_E[ irrep ], newsize_E[ irrep ], vector_rhs + jump[ ptr1 ], tempvector_rhs + jump[ ptr1 ], num_rhs1 );
1078 const int ptr2 = irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET;
1079 const int num_rhs2 = ( jump[ ptr2 + 1 ] - jump[ ptr2 ] ) / size_E[ irrep ];
1080 assert( num_rhs2 * size_E[ irrep ] == jump[ ptr2 + 1 ] - jump[ ptr2 ] );
1081 helper[ ptr2 ] = num_rhs2 * newsize_E[ irrep ];
1082 recreatehelper3( SEE[ irrep ], size_E[ irrep ], newsize_E[ irrep ], vector_rhs + jump[ ptr2 ], tempvector_rhs + jump[ ptr2 ], num_rhs2 );
1085 if ( newsize_F_singlet[ irrep ] > 0 ){
1086 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET;
1087 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_F_singlet[ irrep ];
1088 assert( num_rhs * size_F_singlet[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1089 helper[ ptr ] = num_rhs * newsize_F_singlet[ irrep ];
1090 recreatehelper3( SFF_singlet[ irrep ], size_F_singlet[ irrep ], newsize_F_singlet[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1093 if ( newsize_F_triplet[ irrep ] > 0 ){
1094 const int ptr = irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET;
1095 const int num_rhs = ( jump[ ptr + 1 ] - jump[ ptr ] ) / size_F_triplet[ irrep ];
1096 assert( num_rhs * size_F_triplet[ irrep ] == jump[ ptr + 1 ] - jump[ ptr ] );
1097 helper[ ptr ] = num_rhs * newsize_F_triplet[ irrep ];
1098 recreatehelper3( SFF_triplet[ irrep ], size_F_triplet[ irrep ], newsize_F_triplet[ irrep ], vector_rhs + jump[ ptr ], tempvector_rhs + jump[ ptr ], num_rhs );
1101 if ( newsize_G[ irrep ] > 0 ){
1102 const int ptr1 = irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET;
1103 const int num_rhs1 = ( jump[ ptr1 + 1 ] - jump[ ptr1 ] ) / size_G[ irrep ];
1104 assert( num_rhs1 * size_G[ irrep ] == jump[ ptr1 + 1 ] - jump[ ptr1 ] );
1105 helper[ ptr1 ] = num_rhs1 * newsize_G[ irrep ];
1106 recreatehelper3( SGG[ irrep ], size_G[ irrep ], newsize_G[ irrep ], vector_rhs + jump[ ptr1 ], tempvector_rhs + jump[ ptr1 ], num_rhs1 );
1107 const int ptr2 = irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET;
1108 const int num_rhs2 = ( jump[ ptr2 + 1 ] - jump[ ptr2 ] ) / size_G[ irrep ];
1109 assert( num_rhs2 * size_G[ irrep ] == jump[ ptr2 + 1 ] - jump[ ptr2 ] );
1110 helper[ ptr2 ] = num_rhs2 * newsize_G[ irrep ];
1111 recreatehelper3( SGG[ irrep ], size_G[ irrep ], newsize_G[ irrep ], vector_rhs + jump[ ptr2 ], tempvector_rhs + jump[ ptr2 ], num_rhs2 );
1114 const int ptr1 = irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET;
1115 helper[ ptr1 ] = jump[ ptr1 + 1 ] - jump[ ptr1 ];
1117 dcopy_( helper + ptr1, vector_rhs + jump[ ptr1 ], &inc1, tempvector_rhs + jump[ ptr1 ], &inc1 );
1118 const int ptr2 = irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET;
1119 helper[ ptr2 ] = jump[ ptr2 + 1 ] - jump[ ptr2 ];
1120 dcopy_( helper + ptr2, vector_rhs + jump[ ptr2 ], &inc1, tempvector_rhs + jump[ ptr2 ], &inc1 );
1124 delete [] vector_rhs;
1125 delete [] size_A; size_A = newsize_A;
1126 delete [] size_C; size_C = newsize_C;
1127 delete [] size_D; size_D = newsize_D;
1128 delete [] size_E; size_E = newsize_E;
1129 delete [] size_G; size_G = newsize_G;
1130 delete [] size_B_singlet; size_B_singlet = newsize_B_singlet;
1131 delete [] size_B_triplet; size_B_triplet = newsize_B_triplet;
1132 delete [] size_F_singlet; size_F_singlet = newsize_F_singlet;
1133 delete [] size_F_triplet; size_F_triplet = newsize_F_triplet;
1135 int * newjump =
new int[ num_irreps * CHEMPS2_CASPT2_NUM_CASES + 1 ];
1137 for (
int cnt = 0; cnt < num_irreps * CHEMPS2_CASPT2_NUM_CASES; cnt++ ){
1138 newjump[ cnt + 1 ] = newjump[ cnt ] + helper[ cnt ];
1140 vector_rhs =
new double[ newjump[ num_irreps * CHEMPS2_CASPT2_NUM_CASES ] ];
1142 for (
int cnt = 0; cnt < num_irreps * CHEMPS2_CASPT2_NUM_CASES; cnt++ ){
1143 dcopy_( helper + cnt, tempvector_rhs + jump[ cnt ], &inc1, vector_rhs + newjump[ cnt ], &inc1 );
1145 delete [] tempvector_rhs;
1150 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
1151 delete [] SAA[ irrep ];
1152 delete [] SCC[ irrep ];
1153 delete [] SDD[ irrep ];
1154 delete [] SEE[ irrep ];
1155 delete [] SGG[ irrep ];
1156 delete [] SBB_singlet[ irrep ];
1157 delete [] SBB_triplet[ irrep ];
1158 delete [] SFF_singlet[ irrep ];
1159 delete [] SFF_triplet[ irrep ];
1166 delete [] SBB_singlet;
1167 delete [] SBB_triplet;
1168 delete [] SFF_singlet;
1169 delete [] SFF_triplet;
1171 double * temp = NULL;
1172 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
1173 temp =
new double[ size_A[ irrep ] ]; dcopy_( size_A + irrep, FAA[ irrep ], &inc1, temp, &inc1 );
delete [] FAA[ irrep ]; FAA[ irrep ] = temp;
1174 temp =
new double[ size_C[ irrep ] ]; dcopy_( size_C + irrep, FCC[ irrep ], &inc1, temp, &inc1 );
delete [] FCC[ irrep ]; FCC[ irrep ] = temp;
1175 temp =
new double[ size_D[ irrep ] ]; dcopy_( size_D + irrep, FDD[ irrep ], &inc1, temp, &inc1 );
delete [] FDD[ irrep ]; FDD[ irrep ] = temp;
1176 temp =
new double[ size_E[ irrep ] ]; dcopy_( size_E + irrep, FEE[ irrep ], &inc1, temp, &inc1 );
delete [] FEE[ irrep ]; FEE[ irrep ] = temp;
1177 temp =
new double[ size_G[ irrep ] ]; dcopy_( size_G + irrep, FGG[ irrep ], &inc1, temp, &inc1 );
delete [] FGG[ irrep ]; FGG[ irrep ] = temp;
1178 temp =
new double[ size_B_singlet[ irrep ] ]; dcopy_( size_B_singlet + irrep, FBB_singlet[ irrep ], &inc1, temp, &inc1 );
delete [] FBB_singlet[ irrep ]; FBB_singlet[ irrep ] = temp;
1179 temp =
new double[ size_B_triplet[ irrep ] ]; dcopy_( size_B_triplet + irrep, FBB_triplet[ irrep ], &inc1, temp, &inc1 );
delete [] FBB_triplet[ irrep ]; FBB_triplet[ irrep ] = temp;
1180 temp =
new double[ size_F_singlet[ irrep ] ]; dcopy_( size_F_singlet + irrep, FFF_singlet[ irrep ], &inc1, temp, &inc1 );
delete [] FFF_singlet[ irrep ]; FFF_singlet[ irrep ] = temp;
1181 temp =
new double[ size_F_triplet[ irrep ] ]; dcopy_( size_F_triplet + irrep, FFF_triplet[ irrep ], &inc1, temp, &inc1 );
delete [] FFF_triplet[ irrep ]; FFF_triplet[ irrep ] = temp;
1184 const int total_size = jump[ num_irreps * CHEMPS2_CASPT2_NUM_CASES ];
1185 cout <<
"CASPT2 : New size V_SD space = " << total_size << endl;
1189 int CheMPS2::CASPT2::get_maxsize()
const{
1192 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
1193 maxsize = max( max( max( max( max( max( max( max( max(
1199 size_B_singlet[irrep] ),
1200 size_B_triplet[irrep] ),
1201 size_F_singlet[irrep] ),
1202 size_F_triplet[irrep] ),
1205 if ( maxsize <= 2 ){ maxsize = 3; }
1210 void CheMPS2::CASPT2::matmat(
char totrans,
int rowdim,
int coldim,
int sumdim,
double alpha,
double * matrix,
int ldaM,
double * origin,
int ldaO,
double * target,
int ldaT ){
1214 dgemm_( &totrans, ¬rans, &rowdim, &coldim, &sumdim, &alpha, matrix, &ldaM, origin, &ldaO, &add, target, &ldaT );
1218 void CheMPS2::CASPT2::matvec(
double * vector,
double * result,
double * diag_fock )
const{
1240 const int vectorlength = jump[ CHEMPS2_CASPT2_NUM_CASES * num_irreps ];
1242 for (
int elem = 0; elem < vectorlength; elem++ ){ result[ elem ] = diag_fock[ elem ] * vector[ elem ]; }
1243 const int maxlinsize = get_maxsize();
1244 double * workspace =
new double[ maxlinsize * maxlinsize ];
1245 const double SQRT2 = sqrt( 2.0 );
1248 for (
int IL = 0; IL < num_irreps; IL++ ){
1249 const int SIZE_L = size_A[ IL ];
1250 const int nocc_ij = indices->
getNOCC( IL );
1251 for (
int IR = 0; IR < num_irreps; IR++ ){
1252 const int SIZE_R = size_D[ IR ];
1254 const int shift = shift_D_nonactive( indices, IL, Iw );
1255 const int nocc_w = indices->
getNOCC( Iw );
1256 const int nact_w = indices->
getNDMRG( Iw );
1257 const int n_oa_w = nocc_w + nact_w;
1258 const int nvir_w = indices->
getNVIRT( Iw );
1259 int total_size = SIZE_L * SIZE_R;
1260 if ( total_size * nocc_ij * nact_w * nvir_w > 0 ){
1261 for (
int ac = 0; ac < nvir_w; ac++ ){
1262 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1263 for (
int w = 0; w < nact_w; w++ ){
1264 double f_wc = fock->
get( Iw, nocc_w + w, n_oa_w + ac );
1266 daxpy_( &total_size, &f_wc, FAD[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1268 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ];
1269 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_R * ( shift + nocc_ij * ac );
1270 matmat(
'N', SIZE_L, nocc_ij, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1271 matmat(
'T', SIZE_R, nocc_ij, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1278 for (
int IL = 0; IL < num_irreps; IL++ ){
1279 const int SIZE_L = size_C[ IL ];
1280 const int nvir_ab = indices->
getNVIRT( IL );
1281 for (
int IR = 0; IR < num_irreps; IR++ ){
1282 const int SIZE_R = size_D[ IR ];
1284 const int shift = shift_D_nonactive( indices, Iw, IL );
1285 const int nocc_w = indices->
getNOCC( Iw );
1286 const int nact_w = indices->
getNDMRG( Iw );
1287 int total_size = SIZE_L * SIZE_R;
1288 if ( total_size * nvir_ab * nact_w * nocc_w > 0 ){
1289 for (
int ik = 0; ik < nocc_w; ik++ ){
1290 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1291 for (
int w = 0; w < nact_w; w++ ){
1292 double f_kw = fock->
get( Iw, ik, nocc_w + w );
1294 daxpy_( &total_size, &f_kw, FCD[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1296 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ];
1297 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_R * ( shift + ik );
1298 const int LDA_R = SIZE_R * nocc_w;
1299 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1300 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1307 for (
int IL = 0; IL < num_irreps; IL++ ){
1308 const int SIZE_L = size_A[ IL ];
1309 const int nocc_l = indices->
getNOCC( IL );
1310 for (
int IR = 0; IR < num_irreps; IR++ ){
1311 const int SIZE_R = size_B_singlet[ IR ];
1313 const int shift = (( Iw < IL ) ? shift_B_nonactive( indices, Iw, IL, +1 ) : shift_B_nonactive( indices, IL, Iw, +1 ));
1314 const int nocc_w = indices->
getNOCC( Iw );
1315 const int nact_w = indices->
getNDMRG( Iw );
1316 int total_size = SIZE_L * SIZE_R;
1317 if ( total_size * nocc_l * nact_w * nocc_w > 0 ){
1318 for (
int k = 0; k < nocc_w; k++ ){
1319 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1320 for (
int w = 0; w < nact_w; w++ ){
1321 double f_kw = fock->
get( Iw, k, nocc_w + w );
1323 daxpy_( &total_size, &f_kw, FAB_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1327 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ];
1328 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + SIZE_R * ( shift + ( k * ( k + 1 ) ) / 2 );
1329 matmat(
'N', SIZE_L, k, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1330 matmat(
'T', SIZE_R, k, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1332 #pragma omp parallel for schedule(static) 1333 for (
int l = k; l < nocc_l; l++ ){
1334 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ] + SIZE_L * l;
1335 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + SIZE_R * ( shift + k + ( l * ( l + 1 ) ) / 2 );
1336 const double factor = (( k == l ) ? SQRT2 : 1.0 );
1337 matmat(
'N', SIZE_L, 1, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1338 matmat(
'T', SIZE_R, 1, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1341 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ];
1342 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + SIZE_R * ( shift + (( Iw < IL ) ? k : nocc_l * k ));
1343 const int LDA_R = (( Iw < IL ) ? SIZE_R * nocc_w : SIZE_R );
1344 matmat(
'N', SIZE_L, nocc_l, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1345 matmat(
'T', SIZE_R, nocc_l, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1353 for (
int IL = 0; IL < num_irreps; IL++ ){
1354 const int SIZE_L = size_A[ IL ];
1355 const int nocc_l = indices->
getNOCC( IL );
1356 for (
int IR = 0; IR < num_irreps; IR++ ){
1357 const int SIZE_R = size_B_triplet[ IR ];
1359 const int shift = (( Iw < IL ) ? shift_B_nonactive( indices, Iw, IL, -1 ) : shift_B_nonactive( indices, IL, Iw, -1 ));
1360 const int nocc_w = indices->
getNOCC( Iw );
1361 const int nact_w = indices->
getNDMRG( Iw );
1362 int total_size = SIZE_L * SIZE_R;
1363 if ( total_size * nocc_l * nact_w * nocc_w > 0 ){
1364 for (
int k = 0; k < nocc_w; k++ ){
1365 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1366 for (
int w = 0; w < nact_w; w++ ){
1367 double f_kw = fock->
get( Iw, k, nocc_w + w );
1369 daxpy_( &total_size, &f_kw, FAB_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1373 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ];
1374 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + SIZE_R * ( shift + ( k * ( k - 1 ) ) / 2 );
1375 matmat(
'N', SIZE_L, k, SIZE_R, -1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1376 matmat(
'T', SIZE_R, k, SIZE_L, -1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1378 #pragma omp parallel for schedule(static) 1379 for (
int l = k+1; l < nocc_l; l++ ){
1380 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ] + SIZE_L * l;
1381 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + SIZE_R * ( shift + k + ( l * ( l - 1 ) ) / 2 );
1382 matmat(
'N', SIZE_L, 1, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1383 matmat(
'T', SIZE_R, 1, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1386 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_A ];
1387 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + SIZE_R * ( shift + (( Iw < IL ) ? k : nocc_l * k ));
1388 const int LDA_R = (( Iw < IL ) ? SIZE_R * nocc_w : SIZE_R );
1389 const double factor = (( Iw < IL ) ? 1.0 : -1.0 );
1390 matmat(
'N', SIZE_L, nocc_l, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1391 matmat(
'T', SIZE_R, nocc_l, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1399 for (
int IL = 0; IL < num_irreps; IL++ ){
1400 const int SIZE_L = size_C[ IL ];
1401 const int nvir_d = indices->
getNVIRT( IL );
1402 for (
int IR = 0; IR < num_irreps; IR++ ){
1403 const int SIZE_R = size_F_singlet[ IR ];
1405 const int shift = (( Iw < IL ) ? shift_F_nonactive( indices, Iw, IL, +1 ) : shift_F_nonactive( indices, IL, Iw, +1 ));
1406 const int nocc_w = indices->
getNOCC( Iw );
1407 const int nact_w = indices->
getNDMRG( Iw );
1408 const int n_oa_w = nocc_w + nact_w;
1409 const int nvir_w = indices->
getNVIRT( Iw );
1410 int total_size = SIZE_L * SIZE_R;
1411 if ( total_size * nvir_d * nact_w * nvir_w > 0 ){
1412 for (
int c = 0; c < nvir_w; c++ ){
1413 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1414 for (
int w = 0; w < nact_w; w++ ){
1415 double f_wc = fock->
get( Iw, nocc_w + w, n_oa_w + c );
1417 daxpy_( &total_size, &f_wc, FCF_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1421 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ];
1422 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + SIZE_R * ( shift + ( c * ( c + 1 ) ) / 2 );
1423 matmat(
'N', SIZE_L, c, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1424 matmat(
'T', SIZE_R, c, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1426 #pragma omp parallel for schedule(static) 1427 for (
int d = c; d < nvir_d; d++ ){
1428 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ] + SIZE_L * d;
1429 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + SIZE_R * ( shift + c + ( d * ( d + 1 ) ) / 2 );
1430 const double factor = (( c == d ) ? SQRT2 : 1.0 );
1431 matmat(
'N', SIZE_L, 1, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1432 matmat(
'T', SIZE_R, 1, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1435 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ];
1436 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + SIZE_R * ( shift + (( Iw < IL ) ? c : nvir_d * c ));
1437 const int LDA_R = (( Iw < IL ) ? SIZE_R * nvir_w : SIZE_R );
1438 matmat(
'N', SIZE_L, nvir_d, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1439 matmat(
'T', SIZE_R, nvir_d, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1447 for (
int IL = 0; IL < num_irreps; IL++ ){
1448 const int SIZE_L = size_C[ IL ];
1449 const int nvir_d = indices->
getNVIRT( IL );
1450 for (
int IR = 0; IR < num_irreps; IR++ ){
1451 const int SIZE_R = size_F_triplet[ IR ];
1453 const int shift = (( Iw < IL ) ? shift_F_nonactive( indices, Iw, IL, -1 ) : shift_F_nonactive( indices, IL, Iw, -1 ));
1454 const int nocc_w = indices->
getNOCC( Iw );
1455 const int nact_w = indices->
getNDMRG( Iw );
1456 const int n_oa_w = nocc_w + nact_w;
1457 const int nvir_w = indices->
getNVIRT( Iw );
1458 int total_size = SIZE_L * SIZE_R;
1459 if ( total_size * nvir_d * nact_w * nvir_w > 0 ){
1460 for (
int c = 0; c < nvir_w; c++ ){
1461 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1462 for (
int w = 0; w < nact_w; w++ ){
1463 double f_wc = fock->
get( Iw, nocc_w + w, n_oa_w + c );
1465 daxpy_( &total_size, &f_wc, FCF_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1469 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ];
1470 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + SIZE_R * ( shift + ( c * ( c - 1 ) ) / 2 );
1471 matmat(
'N', SIZE_L, c, SIZE_R, -1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1472 matmat(
'T', SIZE_R, c, SIZE_L, -1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1474 #pragma omp parallel for schedule(static) 1475 for (
int d = c+1; d < nvir_d; d++ ){
1476 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ] + SIZE_L * d;
1477 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + SIZE_R * ( shift + c + ( d * ( d - 1 ) ) / 2 );
1478 matmat(
'N', SIZE_L, 1, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
1479 matmat(
'T', SIZE_R, 1, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
1482 const double factor = (( Iw < IL ) ? 1.0 : -1.0 );
1483 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_C ];
1484 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + SIZE_R * ( shift + (( Iw < IL ) ? c : nvir_d * c ));
1485 const int LDA_R = (( Iw < IL ) ? SIZE_R * nvir_w : SIZE_R );
1486 matmat(
'N', SIZE_L, nvir_d, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1487 matmat(
'T', SIZE_R, nvir_d, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1495 for (
int IL = 0; IL < num_irreps; IL++ ){
1496 const int SIZE_L = size_B_singlet[ IL ];
1497 for (
int IR = 0; IR < num_irreps; IR++ ){
1498 const int SIZE_R = size_E[ IR ];
1500 const int nocc_w = indices->
getNOCC( Iw );
1501 const int nact_w = indices->
getNDMRG( Iw );
1502 const int n_oa_w = nocc_w + nact_w;
1503 const int nvir_w = indices->
getNVIRT( Iw );
1505 for (
int Iik = 0; Iik < num_irreps; Iik++ ){
1508 const int nocc_ik = indices->
getNOCC( Iik );
1509 linsize += (( IL == 0 ) ? ( nocc_ik * ( nocc_ik + 1 ) ) / 2 : nocc_ik * indices->
getNOCC( Ijl ) );
1512 const int size_ij = linsize;
1513 int total_size = SIZE_L * SIZE_R;
1514 if ( total_size * nact_w * nvir_w * size_ij > 0 ){
1515 const int shift_E = shift_E_nonactive( indices, Iw, 0, IL, +1 );
1516 for (
int ac = 0; ac < nvir_w; ac++ ){
1517 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1518 for (
int w = 0; w < nact_w; w++ ){
1519 double f_wc = fock->
get( Iw, nocc_w + w, n_oa_w + ac );
1521 daxpy_( &total_size, &f_wc, FBE_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1523 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_B_SINGLET ];
1524 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE_R * ( shift_E + ac );
1525 const int LDA_R = SIZE_R * nvir_w;
1526 matmat(
'N', SIZE_L, size_ij, SIZE_R, 2.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1527 matmat(
'T', SIZE_R, size_ij, SIZE_L, 2.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1534 for (
int IL = 0; IL < num_irreps; IL++ ){
1535 const int SIZE_L = size_B_triplet[ IL ];
1536 for (
int IR = 0; IR < num_irreps; IR++ ){
1537 const int SIZE_R = size_E[ IR ];
1539 const int nocc_w = indices->
getNOCC( Iw );
1540 const int nact_w = indices->
getNDMRG( Iw );
1541 const int n_oa_w = nocc_w + nact_w;
1542 const int nvir_w = indices->
getNVIRT( Iw );
1544 for (
int Iik = 0; Iik < num_irreps; Iik++ ){
1547 const int nocc_ik = indices->
getNOCC( Iik );
1548 linsize += (( IL == 0 ) ? ( nocc_ik * ( nocc_ik - 1 ) ) / 2 : nocc_ik * indices->
getNOCC( Ijl ) );
1551 const int size_ij = linsize;
1552 int total_size = SIZE_L * SIZE_R;
1553 if ( total_size * nact_w * nvir_w * size_ij > 0 ){
1554 const int shift_E = shift_E_nonactive( indices, Iw, 0, IL, -1 );
1555 for (
int ac = 0; ac < nvir_w; ac++ ){
1556 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1557 for (
int w = 0; w < nact_w; w++ ){
1558 double f_wc = fock->
get( Iw, nocc_w + w, n_oa_w + ac );
1560 daxpy_( &total_size, &f_wc, FBE_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1562 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ];
1563 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE_R * ( shift_E + ac );
1564 const int LDA_R = SIZE_R * nvir_w;
1565 matmat(
'N', SIZE_L, size_ij, SIZE_R, 2.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1566 matmat(
'T', SIZE_R, size_ij, SIZE_L, 2.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1573 for (
int IL = 0; IL < num_irreps; IL++ ){
1574 const int SIZE_L = size_F_singlet[ IL ];
1575 for (
int IR = 0; IR < num_irreps; IR++ ){
1576 const int SIZE_R = size_G[ IR ];
1578 const int nocc_w = indices->
getNOCC( Iw );
1579 const int nact_w = indices->
getNDMRG( Iw );
1581 for (
int Iac = 0; Iac < num_irreps; Iac++ ){
1584 const int nvir_ac = indices->
getNVIRT( Iac );
1585 linsize += (( IL == 0 ) ? ( nvir_ac * ( nvir_ac + 1 ) ) / 2 : nvir_ac * indices->
getNVIRT( Ibd ) );
1588 const int size_ab = linsize;
1589 int total_size = SIZE_L * SIZE_R;
1590 if ( total_size * nact_w * nocc_w * size_ab > 0 ){
1591 const int shift_G = shift_G_nonactive( indices, Iw, 0, IL, +1 );
1592 for (
int ik = 0; ik < nocc_w; ik++ ){
1593 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1594 for (
int w = 0; w < nact_w; w++ ){
1595 double f_kw = fock->
get( Iw, ik, nocc_w + w );
1597 daxpy_( &total_size, &f_kw, FFG_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1599 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_F_SINGLET ];
1600 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE_R * ( shift_G + ik );
1601 const int LDA_R = SIZE_R * nocc_w;
1602 matmat(
'N', SIZE_L, size_ab, SIZE_R, 2.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1603 matmat(
'T', SIZE_R, size_ab, SIZE_L, 2.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1610 for (
int IL = 0; IL < num_irreps; IL++ ){
1611 const int SIZE_L = size_F_triplet[ IL ];
1612 for (
int IR = 0; IR < num_irreps; IR++ ){
1613 const int SIZE_R = size_G[ IR ];
1615 const int nocc_w = indices->
getNOCC( Iw );
1616 const int nact_w = indices->
getNDMRG( Iw );
1618 for (
int Iac = 0; Iac < num_irreps; Iac++ ){
1621 const int nvir_ac = indices->
getNVIRT( Iac );
1622 linsize += (( IL == 0 ) ? ( nvir_ac * ( nvir_ac - 1 ) ) / 2 : nvir_ac * indices->
getNVIRT( Ibd ) );
1625 const int size_ab = linsize;
1626 int total_size = SIZE_L * SIZE_R;
1627 if ( total_size * nact_w * nocc_w * size_ab > 0 ){
1628 const int shift_G = shift_G_nonactive( indices, Iw, 0, IL, -1 );
1629 for (
int ik = 0; ik < nocc_w; ik++ ){
1630 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1631 for (
int w = 0; w < nact_w; w++ ){
1632 double f_kw = fock->
get( Iw, ik, nocc_w + w );
1634 daxpy_( &total_size, &f_kw, FFG_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
1636 const int ptr_L = jump[ IL + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ];
1637 const int ptr_R = jump[ IR + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE_R * ( shift_G + ik );
1638 const int LDA_R = SIZE_R * nocc_w;
1639 matmat(
'N', SIZE_L, size_ab, SIZE_R, 2.0, workspace, SIZE_L, vector + ptr_R, LDA_R, result + ptr_L, SIZE_L );
1640 matmat(
'T', SIZE_R, size_ab, SIZE_L, 2.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, LDA_R );
1647 for (
int IL = 0; IL < num_irreps; IL++ ){
1648 int SIZE = size_E[ IL ];
1649 const int nocc_w = indices->
getNOCC( IL );
1650 const int nact_w = indices->
getNDMRG( IL );
1651 const int n_oa_w = nocc_w + nact_w;
1652 const int nvir_w = indices->
getNVIRT( IL );
1653 if ( SIZE * nact_w * nvir_w > 0 ){
1654 for (
int c = 0; c < nvir_w; c++ ){
1657 for (
int cnt = 0; cnt < SIZE; cnt++ ){ workspace[ cnt ] = 0.0; }
1658 for (
int w = 0; w < nact_w; w++ ){
1659 double f_wc = fock->
get( IL, nocc_w + w, n_oa_w + c );
1661 daxpy_( &SIZE, &f_wc, FEH[ IL ][ w ], &inc1, workspace, &inc1 );
1664 for (
int Id = 0; Id < num_irreps; Id++ ){
1667 const int nvir_d = indices->
getNVIRT( Id );
1668 const int LDA_E = SIZE * nvir_d;
1671 for (
int Ii = 0; Ii < num_irreps; Ii++ ){
1674 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE * shift_E_nonactive( indices, Id, Ii, Ij, +1 );
1675 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, Ii, Ij, IL, Id, +1 );
1676 const int size_ij = indices->
getNOCC( Ii ) * indices->
getNOCC( Ij );
1677 #pragma omp parallel for schedule(static) 1678 for (
int d = 0; d < nvir_d; d++ ){
1679 const int ptr_H = jump_H + size_ij * ( c + nvir_w * d );
1680 const int ptr_E = jump_E + SIZE * d;
1681 matmat(
'N', SIZE, size_ij, 1, 2.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1682 matmat(
'T', 1, size_ij, SIZE, 2.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1689 for (
int Iij = 0; Iij < num_irreps; Iij++ ){
1690 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE * shift_E_nonactive( indices, Id, Iij, Iij, +1 );
1691 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, Iij, Iij, IL, Id, +1 );
1692 const int size_ij = ( indices->
getNOCC( Iij ) * ( indices->
getNOCC( Iij ) + 1 ) ) / 2;
1693 #pragma omp parallel for schedule(static) 1694 for (
int d = 0; d < c; d++ ){
1695 const int ptr_H = jump_H + size_ij * ( d + ( c * ( c + 1 ) ) / 2 );
1696 const int ptr_E = jump_E + SIZE * d;
1697 matmat(
'N', SIZE, size_ij, 1, 2.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1698 matmat(
'T', 1, size_ij, SIZE, 2.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1700 #pragma omp parallel for schedule(static) 1701 for (
int d = c; d < nvir_d; d++ ){
1702 const int ptr_H = jump_H + size_ij * ( c + ( d * ( d + 1 ) ) / 2 );
1703 const int ptr_E = jump_E + SIZE * d;
1704 const double factor = 2 * (( c == d ) ? SQRT2 : 1.0 );
1705 matmat(
'N', SIZE, size_ij, 1, factor, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1706 matmat(
'T', 1, size_ij, SIZE, factor, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1712 for (
int Ii = 0; Ii < num_irreps; Ii++ ){
1715 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE * shift_E_nonactive( indices, Id, Ii, Ij, +1 );
1716 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, Ii, Ij, Id, IL, +1 );
1717 const int size_ij = indices->
getNOCC( Ii ) * indices->
getNOCC( Ij );
1718 #pragma omp parallel for schedule(static) 1719 for (
int d = 0; d < nvir_d; d++ ){
1720 const int ptr_H = jump_H + size_ij * ( d + nvir_d * c );
1721 const int ptr_E = jump_E + SIZE * d;
1722 matmat(
'N', SIZE, size_ij, 1, 2.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1723 matmat(
'T', 1, size_ij, SIZE, 2.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1734 for (
int IL = 0; IL < num_irreps; IL++ ){
1735 int SIZE = size_E[ IL ];
1736 const int nocc_w = indices->
getNOCC( IL );
1737 const int nact_w = indices->
getNDMRG( IL );
1738 const int n_oa_w = nocc_w + nact_w;
1739 const int nvir_w = indices->
getNVIRT( IL );
1740 if ( SIZE * nact_w * nvir_w > 0 ){
1741 for (
int c = 0; c < nvir_w; c++ ){
1744 for (
int cnt = 0; cnt < SIZE; cnt++ ){ workspace[ cnt ] = 0.0; }
1745 for (
int w = 0; w < nact_w; w++ ){
1746 double f_wc = fock->
get( IL, nocc_w + w, n_oa_w + c );
1748 daxpy_( &SIZE, &f_wc, FEH[ IL ][ w ], &inc1, workspace, &inc1 );
1751 for (
int Id = 0; Id < num_irreps; Id++ ){
1754 const int nvir_d = indices->
getNVIRT( Id );
1755 const int LDA_E = SIZE * nvir_d;
1758 for (
int Ii = 0; Ii < num_irreps; Ii++ ){
1761 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE * shift_E_nonactive( indices, Id, Ii, Ij, -1 );
1762 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, Ii, Ij, IL, Id, -1 );
1763 const int size_ij = indices->
getNOCC( Ii ) * indices->
getNOCC( Ij );
1764 #pragma omp parallel for schedule(static) 1765 for (
int d = 0; d < nvir_d; d++ ){
1766 const int ptr_H = jump_H + size_ij * ( c + nvir_w * d );
1767 const int ptr_E = jump_E + SIZE * d;
1768 matmat(
'N', SIZE, size_ij, 1, 6.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1769 matmat(
'T', 1, size_ij, SIZE, 6.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1776 for (
int Iij = 0; Iij < num_irreps; Iij++ ){
1777 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE * shift_E_nonactive( indices, Id, Iij, Iij, -1 );
1778 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, Iij, Iij, IL, Id, -1 );
1779 const int size_ij = ( indices->
getNOCC( Iij ) * ( indices->
getNOCC( Iij ) - 1 ) ) / 2;
1780 #pragma omp parallel for schedule(static) 1781 for (
int d = 0; d < c; d++ ){
1782 const int ptr_H = jump_H + size_ij * ( d + ( c * ( c - 1 ) ) / 2 );
1783 const int ptr_E = jump_E + SIZE * d;
1784 matmat(
'N', SIZE, size_ij, 1, -6.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1785 matmat(
'T', 1, size_ij, SIZE, -6.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1787 #pragma omp parallel for schedule(static) 1788 for (
int d = c+1; d < nvir_d; d++ ){
1789 const int ptr_H = jump_H + size_ij * ( c + ( d * ( d - 1 ) ) / 2 );
1790 const int ptr_E = jump_E + SIZE * d;
1791 matmat(
'N', SIZE, size_ij, 1, 6.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1792 matmat(
'T', 1, size_ij, SIZE, 6.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1798 for (
int Ii = 0; Ii < num_irreps; Ii++ ){
1801 const int jump_E = jump[ IL + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE * shift_E_nonactive( indices, Id, Ii, Ij, -1 );
1802 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, Ii, Ij, Id, IL, -1 );
1803 const int size_ij = indices->
getNOCC( Ii ) * indices->
getNOCC( Ij );
1804 #pragma omp parallel for schedule(static) 1805 for (
int d = 0; d < nvir_d; d++ ){
1806 const int ptr_H = jump_H + size_ij * ( d + nvir_d * c );
1807 const int ptr_E = jump_E + SIZE * d;
1808 matmat(
'N', SIZE, size_ij, 1, -6.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_E, LDA_E );
1809 matmat(
'T', 1, size_ij, SIZE, -6.0, workspace, SIZE, vector + ptr_E, LDA_E, result + ptr_H, 1 );
1820 for (
int IL = 0; IL < num_irreps; IL++ ){
1821 int SIZE = size_G[ IL ];
1822 const int nocc_w = indices->
getNOCC( IL );
1823 const int nact_w = indices->
getNDMRG( IL );
1824 if ( SIZE * nocc_w * nact_w > 0 ){
1825 for (
int k = 0; k < nocc_w; k++ ){
1828 for (
int cnt = 0; cnt < SIZE; cnt++ ){ workspace[ cnt ] = 0.0; }
1829 for (
int w = 0; w < nact_w; w++ ){
1830 double f_kw = fock->
get( IL, k, nocc_w + w );
1832 daxpy_( &SIZE, &f_kw, FGH[ IL ][ w ], &inc1, workspace, &inc1 );
1835 for (
int Il = 0; Il < num_irreps; Il++ ){
1838 const int nocc_l = indices->
getNOCC( Il );
1841 for (
int Ia = 0; Ia < num_irreps; Ia++ ){
1844 const int colsize = nocc_l * indices->
getNVIRT( Ia ) * indices->
getNVIRT( Ib );
1846 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE * shift_G_nonactive( indices, Il, Ia, Ib, +1 );
1847 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, IL, Il, Ia, Ib, +1 ) + k;
1848 matmat(
'N', SIZE, colsize, 1, 2.0, workspace, SIZE, vector + jump_H, nocc_w, result + jump_G, SIZE );
1849 matmat(
'T', 1, colsize, SIZE, 2.0, workspace, SIZE, vector + jump_G, SIZE, result + jump_H, nocc_w );
1856 for (
int Iab = 0; Iab < num_irreps; Iab++ ){
1857 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE * shift_G_nonactive( indices, Il, Iab, Iab, +1 );
1858 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, IL, Il, Iab, Iab, +1 );
1859 const int size_ij = ( nocc_w * ( nocc_w + 1 ) ) / 2;
1860 const int size_ab = ( indices->
getNVIRT( Iab ) * ( indices->
getNVIRT( Iab ) + 1 ) ) / 2;
1861 const int LDA_G = SIZE * nocc_l;
1862 #pragma omp parallel for schedule(static) 1863 for (
int l = 0; l < k; l++ ){
1864 const int ptr_H = jump_H + ( l + ( k * ( k + 1 ) ) / 2 );
1865 const int ptr_G = jump_G + SIZE * l;
1866 matmat(
'N', SIZE, size_ab, 1, 2.0, workspace, SIZE, vector + ptr_H, size_ij, result + ptr_G, LDA_G );
1867 matmat(
'T', 1, size_ab, SIZE, 2.0, workspace, SIZE, vector + ptr_G, LDA_G, result + ptr_H, size_ij );
1869 #pragma omp parallel for schedule(static) 1870 for (
int l = k; l < nocc_l; l++ ){
1871 const int ptr_H = jump_H + ( k + ( l * ( l + 1 ) ) / 2 );
1872 const int ptr_G = jump_G + SIZE * l;
1873 const double factor = 2 * (( k == l ) ? SQRT2 : 1.0 );
1874 matmat(
'N', SIZE, size_ab, 1, factor, workspace, SIZE, vector + ptr_H, size_ij, result + ptr_G, LDA_G );
1875 matmat(
'T', 1, size_ab, SIZE, factor, workspace, SIZE, vector + ptr_G, LDA_G, result + ptr_H, size_ij );
1881 for (
int Ia = 0; Ia < num_irreps; Ia++ ){
1884 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE * shift_G_nonactive( indices, Il, Ia, Ib, +1 );
1885 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift_H_nonactive( indices, Il, IL, Ia, Ib, +1 );
1887 #pragma omp parallel for schedule(static) 1888 for (
int ab = 0; ab < size_ab; ab++ ){
1889 const int ptr_H = jump_H + nocc_l * ( k + nocc_w * ab );
1890 const int ptr_G = jump_G + SIZE * nocc_l * ab;
1891 matmat(
'N', SIZE, nocc_l, 1, 2.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_G, SIZE );
1892 matmat(
'T', 1, nocc_l, SIZE, 2.0, workspace, SIZE, vector + ptr_G, SIZE, result + ptr_H, 1 );
1903 for (
int IL = 0; IL < num_irreps; IL++ ){
1904 int SIZE = size_G[ IL ];
1905 const int nocc_w = indices->
getNOCC( IL );
1906 const int nact_w = indices->
getNDMRG( IL );
1907 if ( SIZE * nocc_w * nact_w > 0 ){
1908 for (
int k = 0; k < nocc_w; k++ ){
1911 for (
int cnt = 0; cnt < SIZE; cnt++ ){ workspace[ cnt ] = 0.0; }
1912 for (
int w = 0; w < nact_w; w++ ){
1913 double f_kw = fock->
get( IL, k, nocc_w + w );
1915 daxpy_( &SIZE, &f_kw, FGH[ IL ][ w ], &inc1, workspace, &inc1 );
1918 for (
int Il = 0; Il < num_irreps; Il++ ){
1921 const int nocc_l = indices->
getNOCC( Il );
1924 for (
int Ia = 0; Ia < num_irreps; Ia++ ){
1927 const int colsize = nocc_l * indices->
getNVIRT( Ia ) * indices->
getNVIRT( Ib );
1929 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE * shift_G_nonactive( indices, Il, Ia, Ib, -1 );
1930 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, IL, Il, Ia, Ib, -1 ) + k;
1931 matmat(
'N', SIZE, colsize, 1, -6.0, workspace, SIZE, vector + jump_H, nocc_w, result + jump_G, SIZE );
1932 matmat(
'T', 1, colsize, SIZE, -6.0, workspace, SIZE, vector + jump_G, SIZE, result + jump_H, nocc_w );
1939 for (
int Iab = 0; Iab < num_irreps; Iab++ ){
1940 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE * shift_G_nonactive( indices, Il, Iab, Iab, -1 );
1941 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, IL, Il, Iab, Iab, -1 );
1942 const int size_ij = ( nocc_w * ( nocc_w - 1 ) ) / 2;
1943 const int size_ab = ( indices->
getNVIRT( Iab ) * ( indices->
getNVIRT( Iab ) - 1 ) ) / 2;
1944 const int LDA_G = SIZE * nocc_l;
1945 #pragma omp parallel for schedule(static) 1946 for (
int l = 0; l < k; l++ ){
1947 const int ptr_H = jump_H + ( l + ( k * ( k - 1 ) ) / 2 );
1948 const int ptr_G = jump_G + SIZE * l;
1949 matmat(
'N', SIZE, size_ab, 1, 6.0, workspace, SIZE, vector + ptr_H, size_ij, result + ptr_G, LDA_G );
1950 matmat(
'T', 1, size_ab, SIZE, 6.0, workspace, SIZE, vector + ptr_G, LDA_G, result + ptr_H, size_ij );
1952 #pragma omp parallel for schedule(static) 1953 for (
int l = k+1; l < nocc_l; l++ ){
1954 const int ptr_H = jump_H + ( k + ( l * ( l - 1 ) ) / 2 );
1955 const int ptr_G = jump_G + SIZE * l;
1956 matmat(
'N', SIZE, size_ab, 1, -6.0, workspace, SIZE, vector + ptr_H, size_ij, result + ptr_G, LDA_G );
1957 matmat(
'T', 1, size_ab, SIZE, -6.0, workspace, SIZE, vector + ptr_G, LDA_G, result + ptr_H, size_ij );
1963 for (
int Ia = 0; Ia < num_irreps; Ia++ ){
1966 const int jump_G = jump[ IL + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE * shift_G_nonactive( indices, Il, Ia, Ib, -1 );
1967 const int jump_H = jump[ Icenter + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift_H_nonactive( indices, Il, IL, Ia, Ib, -1 );
1969 #pragma omp parallel for schedule(static) 1970 for (
int ab = 0; ab < size_ab; ab++ ){
1971 const int ptr_H = jump_H + nocc_l * ( k + nocc_w * ab );
1972 const int ptr_G = jump_G + SIZE * nocc_l * ab;
1973 matmat(
'N', SIZE, nocc_l, 1, 6.0, workspace, SIZE, vector + ptr_H, 1, result + ptr_G, SIZE );
1974 matmat(
'T', 1, nocc_l, SIZE, 6.0, workspace, SIZE, vector + ptr_G, SIZE, result + ptr_H, 1 );
1985 for (
int IL = 0; IL < num_irreps; IL++ ){
1986 const int SIZE_L = size_D[ IL ];
1987 for (
int IR = 0; IR < num_irreps; IR++ ){
1988 const int SIZE_R = size_E[ IR ];
1990 const int nocc_kw = indices->
getNOCC( Ikw );
1991 const int nact_kw = indices->
getNDMRG( Ikw );
1992 int total_size = SIZE_L * SIZE_R;
1993 if ( total_size * nocc_kw * nact_kw > 0 ){
1994 for (
int k = 0; k < nocc_kw; k++ ){
1995 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
1996 for (
int w = 0; w < nact_kw; w++ ){
1997 double f_kw = fock->
get( Ikw, k, nocc_kw + w );
1999 daxpy_( &total_size, &f_kw, FDE_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
2001 for (
int Iab = 0; Iab < num_irreps; Iab++ ){
2002 const int nvir_ab = indices->
getNVIRT( Iab );
2004 const int nocc_l = indices->
getNOCC( Il );
2005 if ( nvir_ab * nocc_l > 0 ){
2006 const int jump_D = jump[ IL + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_L * shift_D_nonactive( indices, Il, Iab );
2007 const int jump_E = jump[ IR + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE_R * (( Ikw <= Il ) ? shift_E_nonactive( indices, Iab, Ikw, Il, +1 )
2008 : shift_E_nonactive( indices, Iab, Il, Ikw, +1 ));
2010 #pragma omp parallel for schedule(static) 2011 for (
int l = 0; l < nocc_l; l++ ){
2012 const double factor = (( k == l ) ? SQRT2 : 1.0 );
2013 const int ptr_L = jump_D + SIZE_L * l;
2014 const int ptr_R = jump_E + SIZE_R * nvir_ab * (( k < l ) ? ( k + ( l * ( l + 1 ) ) / 2 ) : ( l + ( k * ( k + 1 ) ) / 2 ));
2015 const int LDA_L = SIZE_L * nocc_l;
2016 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, LDA_L );
2017 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, LDA_L, result + ptr_R, SIZE_R );
2020 #pragma omp parallel for schedule(static) 2021 for (
int l = 0; l < nocc_l; l++ ){
2022 const int ptr_L = jump_D + SIZE_L * l;
2023 const int ptr_R = jump_E + SIZE_R * nvir_ab * (( Ikw < Il ) ? ( k + nocc_kw * l ) : ( l + nocc_l * k ));
2024 const int LDA_L = SIZE_L * nocc_l;
2025 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, LDA_L );
2026 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, LDA_L, result + ptr_R, SIZE_R );
2037 for (
int IL = 0; IL < num_irreps; IL++ ){
2038 const int SIZE_L = size_D[ IL ];
2039 for (
int IR = 0; IR < num_irreps; IR++ ){
2040 const int SIZE_R = size_E[ IR ];
2042 const int nocc_kw = indices->
getNOCC( Ikw );
2043 const int nact_kw = indices->
getNDMRG( Ikw );
2044 int total_size = SIZE_L * SIZE_R;
2045 if ( total_size * nocc_kw * nact_kw > 0 ){
2046 for (
int k = 0; k < nocc_kw; k++ ){
2047 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
2048 for (
int w = 0; w < nact_kw; w++ ){
2049 double f_kw = fock->
get( Ikw, k, nocc_kw + w );
2051 daxpy_( &total_size, &f_kw, FDE_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
2053 for (
int Iab = 0; Iab < num_irreps; Iab++ ){
2054 const int nvir_ab = indices->
getNVIRT( Iab );
2056 const int nocc_l = indices->
getNOCC( Il );
2057 if ( nvir_ab * nocc_l > 0 ){
2058 const int jump_D = jump[ IL + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_L * shift_D_nonactive( indices, Il, Iab );
2059 const int jump_E = jump[ IR + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE_R * (( Ikw <= Il ) ? shift_E_nonactive( indices, Iab, Ikw, Il, -1 )
2060 : shift_E_nonactive( indices, Iab, Il, Ikw, -1 ));
2062 #pragma omp parallel for schedule(static) 2063 for (
int l = 0; l < k; l++ ){
2064 const int ptr_L = jump_D + SIZE_L * l;
2065 const int ptr_R = jump_E + SIZE_R * nvir_ab * ( l + ( k * ( k - 1 ) ) / 2 );
2066 const int LDA_L = SIZE_L * nocc_l;
2067 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, -3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, LDA_L );
2068 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, -3.0, workspace, SIZE_L, vector + ptr_L, LDA_L, result + ptr_R, SIZE_R );
2070 #pragma omp parallel for schedule(static) 2071 for (
int l = k+1; l < nocc_l; l++ ){
2072 const int ptr_L = jump_D + SIZE_L * l;
2073 const int ptr_R = jump_E + SIZE_R * nvir_ab * ( k + ( l * ( l - 1 ) ) / 2 );
2074 const int LDA_L = SIZE_L * nocc_l;
2075 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, 3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, LDA_L );
2076 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, 3.0, workspace, SIZE_L, vector + ptr_L, LDA_L, result + ptr_R, SIZE_R );
2079 #pragma omp parallel for schedule(static) 2080 for (
int l = 0; l < nocc_l; l++ ){
2081 const double factor = (( Ikw < Il ) ? 3.0 : -3.0 );
2082 const int ptr_L = jump_D + SIZE_L * l;
2083 const int ptr_R = jump_E + SIZE_R * nvir_ab * (( Ikw < Il ) ? ( k + nocc_kw * l ) : ( l + nocc_l * k ));
2084 const int LDA_L = SIZE_L * nocc_l;
2085 matmat(
'N', SIZE_L, nvir_ab, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, LDA_L );
2086 matmat(
'T', SIZE_R, nvir_ab, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, LDA_L, result + ptr_R, SIZE_R );
2097 for (
int IL = 0; IL < num_irreps; IL++ ){
2098 const int SIZE_L = size_D[ IL ];
2099 for (
int IR = 0; IR < num_irreps; IR++ ){
2100 const int SIZE_R = size_E[ IR ];
2102 const int nocc_wc = indices->
getNOCC( Iwc );
2103 const int nact_wc = indices->
getNDMRG( Iwc );
2104 const int n_oa_wc = nocc_wc + nact_wc;
2105 const int nvir_wc = indices->
getNVIRT( Iwc );
2106 int total_size = SIZE_L * SIZE_R;
2107 if ( total_size * nvir_wc * nact_wc > 0 ){
2108 for (
int c = 0; c < nvir_wc; c++ ){
2109 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
2110 for (
int w = 0; w < nact_wc; w++ ){
2111 double f_wc = fock->
get( Iwc, nocc_wc + w, n_oa_wc + c );
2113 daxpy_( &total_size, &f_wc, FDG_singlet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
2115 for (
int Iij = 0; Iij < num_irreps; Iij++ ){
2116 const int nocc_ij = indices->
getNOCC( Iij );
2118 const int nvir_d = indices->
getNVIRT( Id );
2119 if ( nvir_d * nocc_ij > 0 ){
2120 const int jump_D = jump[ IL + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_L * shift_D_nonactive( indices, Iij, Id );
2121 const int jump_G = jump[ IR + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE_R * (( Iwc <= Id ) ? shift_G_nonactive( indices, Iij, Iwc, Id, +1 )
2122 : shift_G_nonactive( indices, Iij, Id, Iwc, +1 ));
2125 const int ptr_L = jump_D;
2126 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c * ( c + 1 ) ) / 2;
2127 matmat(
'N', SIZE_L, c * nocc_ij, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2128 matmat(
'T', SIZE_R, c * nocc_ij, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2130 #pragma omp parallel for schedule(static) 2131 for (
int d = c; d < nvir_d; d++ ){
2132 const double factor = (( c == d ) ? SQRT2 : 1.0 );
2133 const int ptr_L = jump_D + SIZE_L * nocc_ij * d;
2134 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c + ( d * ( d + 1 ) ) / 2 );
2135 matmat(
'N', SIZE_L, nocc_ij, SIZE_R, factor, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2136 matmat(
'T', SIZE_R, nocc_ij, SIZE_L, factor, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2140 #pragma omp parallel for schedule(static) 2141 for (
int d = 0; d < nvir_d; d++ ){
2142 const int ptr_L = jump_D + SIZE_L * nocc_ij * d;
2143 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c + nvir_wc * d );
2144 matmat(
'N', SIZE_L, nocc_ij, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2145 matmat(
'T', SIZE_R, nocc_ij, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2148 const int ptr_L = jump_D;
2149 const int ptr_R = jump_G + SIZE_R * nocc_ij * nvir_d * c;
2150 matmat(
'N', SIZE_L, nvir_d * nocc_ij, SIZE_R, 1.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2151 matmat(
'T', SIZE_R, nvir_d * nocc_ij, SIZE_L, 1.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2162 for (
int IL = 0; IL < num_irreps; IL++ ){
2163 const int SIZE_L = size_D[ IL ];
2164 for (
int IR = 0; IR < num_irreps; IR++ ){
2165 const int SIZE_R = size_E[ IR ];
2167 const int nocc_wc = indices->
getNOCC( Iwc );
2168 const int nact_wc = indices->
getNDMRG( Iwc );
2169 const int n_oa_wc = nocc_wc + nact_wc;
2170 const int nvir_wc = indices->
getNVIRT( Iwc );
2171 int total_size = SIZE_L * SIZE_R;
2172 if ( total_size * nvir_wc * nact_wc > 0 ){
2173 for (
int c = 0; c < nvir_wc; c++ ){
2174 for (
int cnt = 0; cnt < total_size; cnt++ ){ workspace[ cnt ] = 0.0; }
2175 for (
int w = 0; w < nact_wc; w++ ){
2176 double f_wc = fock->
get( Iwc, nocc_wc + w, n_oa_wc + c );
2178 daxpy_( &total_size, &f_wc, FDG_triplet[ IL ][ IR ][ w ], &inc1, workspace, &inc1 );
2180 for (
int Iij = 0; Iij < num_irreps; Iij++ ){
2181 const int nocc_ij = indices->
getNOCC( Iij );
2183 const int nvir_d = indices->
getNVIRT( Id );
2184 if ( nvir_d * nocc_ij > 0 ){
2185 const int jump_D = jump[ IL + num_irreps * CHEMPS2_CASPT2_D ] + SIZE_L * shift_D_nonactive( indices, Iij, Id );
2186 const int jump_G = jump[ IR + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE_R * (( Iwc <= Id ) ? shift_G_nonactive( indices, Iij, Iwc, Id, -1 )
2187 : shift_G_nonactive( indices, Iij, Id, Iwc, -1 ));
2190 const int ptr_L = jump_D;
2191 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c * ( c - 1 ) ) / 2;
2192 matmat(
'N', SIZE_L, c * nocc_ij, SIZE_R, -3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2193 matmat(
'T', SIZE_R, c * nocc_ij, SIZE_L, -3.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2195 #pragma omp parallel for schedule(static) 2196 for (
int d = c+1; d < nvir_d; d++ ){
2197 const int ptr_L = jump_D + SIZE_L * nocc_ij * d;
2198 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c + ( d * ( d - 1 ) ) / 2 );
2199 matmat(
'N', SIZE_L, nocc_ij, SIZE_R, 3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2200 matmat(
'T', SIZE_R, nocc_ij, SIZE_L, 3.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2204 #pragma omp parallel for schedule(static) 2205 for (
int d = 0; d < nvir_d; d++ ){
2206 const int ptr_L = jump_D + SIZE_L * nocc_ij * d;
2207 const int ptr_R = jump_G + SIZE_R * nocc_ij * ( c + nvir_wc * d );
2208 matmat(
'N', SIZE_L, nocc_ij, SIZE_R, 3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2209 matmat(
'T', SIZE_R, nocc_ij, SIZE_L, 3.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2212 const int ptr_L = jump_D;
2213 const int ptr_R = jump_G + SIZE_R * nocc_ij * nvir_d * c;
2214 matmat(
'N', SIZE_L, nvir_d * nocc_ij, SIZE_R, -3.0, workspace, SIZE_L, vector + ptr_R, SIZE_R, result + ptr_L, SIZE_L );
2215 matmat(
'T', SIZE_R, nvir_d * nocc_ij, SIZE_L, -3.0, workspace, SIZE_L, vector + ptr_L, SIZE_L, result + ptr_R, SIZE_R );
2225 delete [] workspace;
2229 int CheMPS2::CASPT2::shift_H_nonactive(
const DMRGSCFindices * idx,
const int irrep_i,
const int irrep_j,
const int irrep_a,
const int irrep_b,
const int ST ){
2231 assert( irrep_i <= irrep_j );
2232 assert( irrep_a <= irrep_b );
2235 assert( irrep_prod == irrep_virt );
2239 if ( irrep_prod == 0 ){
2240 for (
int Iij = 0; Iij < n_irreps; Iij++ ){
2241 for (
int Iab = 0; Iab < n_irreps; Iab++ ){
2242 if (( irrep_i == Iij ) && ( irrep_j == Iij ) && ( irrep_a == Iab ) && ( irrep_b == Iab )){
2251 for (
int Ii = 0; Ii < n_irreps; Ii++ ){
2254 for (
int Ia = 0; Ia < n_irreps; Ia++ ){
2257 if (( irrep_i == Ii ) && ( irrep_j == Ij ) && ( irrep_a == Ia ) && ( irrep_b == Ib )){
2272 int CheMPS2::CASPT2::shift_G_nonactive(
const DMRGSCFindices * idx,
const int irrep_i,
const int irrep_a,
const int irrep_b,
const int ST ){
2274 assert( irrep_a <= irrep_b );
2280 for (
int Ii = 0; Ii < n_irreps; Ii++ ){
2283 for (
int Iab = 0; Iab < n_irreps; Iab++ ){
2284 if (( irrep_i == Ii ) && ( irrep_a == Iab ) && ( irrep_b == Iab )){
2292 for (
int Ia = 0; Ia < n_irreps; Ia++ ){
2295 if (( irrep_i == Ii ) && ( irrep_a == Ia ) && ( irrep_b == Ib )){
2309 int CheMPS2::CASPT2::shift_E_nonactive(
const DMRGSCFindices * idx,
const int irrep_a,
const int irrep_i,
const int irrep_j,
const int ST ){
2311 assert( irrep_i <= irrep_j );
2317 for (
int Ia = 0; Ia < n_irreps; Ia++ ){
2320 for (
int Iij = 0; Iij < n_irreps; Iij++ ){
2321 if (( irrep_a == Ia ) && ( irrep_i == Iij ) && ( irrep_j == Iij )){
2329 for (
int Ii = 0; Ii < n_irreps; Ii++ ){
2332 if (( irrep_a == Ia ) && ( irrep_i == Ii ) && ( irrep_j == Ij )){
2346 int CheMPS2::CASPT2::shift_F_nonactive(
const DMRGSCFindices * idx,
const int irrep_a,
const int irrep_b,
const int ST ){
2348 assert( irrep_a <= irrep_b );
2353 if ( irr_prod == 0 ){
2354 for (
int Iab = 0; Iab < n_irreps; Iab++ ){
2355 if (( irrep_a == Iab ) && ( irrep_b == Iab )){
2362 for (
int Ia = 0; Ia < n_irreps; Ia++ ){
2365 if (( irrep_a == Ia ) && ( irrep_b == Ib )){
2377 int CheMPS2::CASPT2::shift_B_nonactive(
const DMRGSCFindices * idx,
const int irrep_i,
const int irrep_j,
const int ST ){
2379 assert( irrep_i <= irrep_j );
2384 if ( irr_prod == 0 ){
2385 for (
int Iij = 0; Iij < n_irreps; Iij++ ){
2386 if (( Iij == irrep_i ) && ( Iij == irrep_j )){
2389 shift += ( idx->
getNOCC( Iij ) * ( idx->
getNOCC( Iij ) + ST ) ) / 2;
2393 for (
int Ii = 0; Ii < n_irreps; Ii++ ){
2396 if (( Ii == irrep_i ) && ( Ij == irrep_j )){
2408 int CheMPS2::CASPT2::shift_D_nonactive(
const DMRGSCFindices * idx,
const int irrep_i,
const int irrep_a ){
2414 for (
int Ii = 0; Ii < n_irreps; Ii++ ){
2416 if (( Ii == irrep_i ) && ( Ia == irrep_a )){
2426 void CheMPS2::CASPT2::diagonal(
double * result )
const{
2428 #pragma omp parallel 2432 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2433 const int SIZE = size_A[ irrep ];
2435 const int NOCC = indices->
getNOCC( irrep );
2436 #pragma omp for schedule(static) 2437 for (
int count = 0; count < NOCC; count++ ){
2438 const double beta = - f_dot_1dm - fock->
get( irrep, count, count );
2439 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_A ] + SIZE * count;
2441 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = FAA[ irrep ][ elem ] + beta; }
2447 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2448 const int SIZE = size_C[ irrep ];
2450 const int NVIR = indices->
getNVIRT( irrep );
2451 const int N_OA = indices->
getNOCC( irrep ) + indices->
getNDMRG( irrep );
2452 #pragma omp for schedule(static) 2453 for (
int count = 0; count < NVIR; count++ ){
2454 const double beta = - f_dot_1dm + fock->
get( irrep, N_OA + count, N_OA + count );
2455 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_C ] + SIZE * count;
2457 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = FCC[ irrep ][ elem ] + beta; }
2463 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2464 const int SIZE = size_D[ irrep ];
2467 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2469 const int NOCC_i = indices->
getNOCC( irrep_i );
2470 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2471 const int SIZE_ia = NOCC_i * indices->
getNVIRT( irrep_a );
2472 #pragma omp for schedule(static) 2473 for (
int combined = 0; combined < SIZE_ia; combined++ ){
2474 const int i = combined % NOCC_i;
2475 const int a = combined / NOCC_i;
2476 const double f_ii = fock->
get( irrep_i, i, i );
2477 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2478 const double beta = - f_dot_1dm + f_aa - f_ii;
2479 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_D ] + SIZE * ( shift + combined );
2481 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = FDD[ irrep ][ elem ] + beta; }
2488 int triangle_idx[] = { 0, 0 };
2492 const int SIZE = size_B_singlet[ 0 ];
2495 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
2496 const int nocc_ij = indices->
getNOCC( irrep_ij );
2497 const int size_ij = ( nocc_ij * ( nocc_ij + 1 ) ) / 2;
2498 #pragma omp for schedule(static) 2499 for (
int combined = 0; combined < size_ij; combined++ ){
2501 const int i = triangle_idx[ 0 ];
2502 const int j = triangle_idx[ 1 ];
2503 const double f_ii = fock->
get( irrep_ij, i, i );
2504 const double f_jj = fock->
get( irrep_ij, j, j );
2505 const double beta = - f_dot_1dm - f_ii - f_jj;
2506 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + SIZE * ( shift + combined );
2508 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FBB_singlet[ 0 ][ elem ] + beta ); }
2514 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
2515 const int SIZE = size_B_singlet[ irrep ];
2518 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2520 if ( irrep_i < irrep_j ){
2521 const int nocc_i = indices->
getNOCC( irrep_i );
2522 const int size_ij = nocc_i * indices->
getNOCC( irrep_j );
2523 #pragma omp for schedule(static) 2524 for (
int combined = 0; combined < size_ij; combined++ ){
2525 const int i = combined % nocc_i;
2526 const int j = combined / nocc_i;
2527 const double f_ii = fock->
get( irrep_i, i, i );
2528 const double f_jj = fock->
get( irrep_j, j, j );
2529 const double beta = - f_dot_1dm - f_ii - f_jj;
2530 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + SIZE * ( shift + combined );
2532 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FBB_singlet[ irrep ][ elem ] + beta ); }
2542 const int SIZE = size_B_triplet[ 0 ];
2545 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
2546 const int nocc_ij = indices->
getNOCC( irrep_ij );
2547 const int size_ij = ( nocc_ij * ( nocc_ij - 1 ) ) / 2;
2548 #pragma omp for schedule(static) 2549 for (
int combined = 0; combined < size_ij; combined++ ){
2551 const int i = triangle_idx[ 0 ];
2552 const int j = triangle_idx[ 1 ];
2553 const double f_ii = fock->
get( irrep_ij, i, i );
2554 const double f_jj = fock->
get( irrep_ij, j, j );
2555 const double beta = - f_dot_1dm - f_ii - f_jj;
2556 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + SIZE * ( shift + combined );
2558 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FBB_triplet[ 0 ][ elem ] + beta ); }
2564 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
2565 const int SIZE = size_B_triplet[ irrep ];
2568 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2570 if ( irrep_i < irrep_j ){
2571 const int nocc_i = indices->
getNOCC( irrep_i );
2572 const int size_ij = nocc_i * indices->
getNOCC( irrep_j );
2573 #pragma omp for schedule(static) 2574 for (
int combined = 0; combined < size_ij; combined++ ){
2575 const int i = combined % nocc_i;
2576 const int j = combined / nocc_i;
2577 const double f_ii = fock->
get( irrep_i, i, i );
2578 const double f_jj = fock->
get( irrep_j, j, j );
2579 const double beta = - f_dot_1dm - f_ii - f_jj;
2580 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + SIZE * ( shift + combined );
2582 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FBB_triplet[ irrep ][ elem ] + beta ); }
2592 const int SIZE = size_F_singlet[ 0 ];
2595 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
2596 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
2597 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
2598 const int SIZE_ab = ( NVIR_ab * ( NVIR_ab + 1 ) ) / 2;
2599 #pragma omp for schedule(static) 2600 for (
int combined = 0; combined < SIZE_ab; combined++ ){
2602 const int a = triangle_idx[ 0 ];
2603 const int b = triangle_idx[ 1 ];
2604 const double f_aa = fock->
get( irrep_ab, N_OA_ab + a, N_OA_ab + a );
2605 const double f_bb = fock->
get( irrep_ab, N_OA_ab + b, N_OA_ab + b );
2606 const double beta = - f_dot_1dm + f_aa + f_bb;
2607 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + SIZE * ( shift + combined );
2609 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FFF_singlet[ 0 ][ elem ] + beta ); }
2615 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
2616 const int SIZE = size_F_singlet[ irrep ];
2619 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2621 if ( irrep_a < irrep_b ){
2622 const int NVIR_a = indices->
getNVIRT( irrep_a );
2623 const int SIZE_ab = NVIR_a * indices->
getNVIRT( irrep_b );
2624 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2625 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
2626 #pragma omp for schedule(static) 2627 for (
int combined = 0; combined < SIZE_ab; combined++ ){
2628 const int a = combined % NVIR_a;
2629 const int b = combined / NVIR_a;
2630 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2631 const double f_bb = fock->
get( irrep_b, N_OA_b + b, N_OA_b + b );
2632 const double beta = - f_dot_1dm + f_aa + f_bb;
2633 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + SIZE * ( shift + combined );
2635 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FFF_singlet[ irrep ][ elem ] + beta ); }
2645 const int SIZE = size_F_triplet[ 0 ];
2648 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
2649 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
2650 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
2651 const int SIZE_ab = ( NVIR_ab * ( NVIR_ab - 1 ) ) / 2;
2652 #pragma omp for schedule(static) 2653 for (
int combined = 0; combined < SIZE_ab; combined++ ){
2655 const int a = triangle_idx[ 0 ];
2656 const int b = triangle_idx[ 1 ];
2657 const double f_aa = fock->
get( irrep_ab, N_OA_ab + a, N_OA_ab + a );
2658 const double f_bb = fock->
get( irrep_ab, N_OA_ab + b, N_OA_ab + b );
2659 const double beta = - f_dot_1dm + f_aa + f_bb;
2660 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + SIZE * ( shift + combined );
2662 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FFF_triplet[ 0 ][ elem ] + beta ); }
2668 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
2669 const int SIZE = size_F_triplet[ irrep ];
2672 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2674 if ( irrep_a < irrep_b ){
2675 const int NVIR_a = indices->
getNVIRT( irrep_a );
2676 const int SIZE_ab = NVIR_a * indices->
getNVIRT( irrep_b );
2677 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2678 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
2679 #pragma omp for schedule(static) 2680 for (
int combined = 0; combined < SIZE_ab; combined++ ){
2681 const int a = combined % NVIR_a;
2682 const int b = combined / NVIR_a;
2683 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2684 const double f_bb = fock->
get( irrep_b, N_OA_b + b, N_OA_b + b );
2685 const double beta = - f_dot_1dm + f_aa + f_bb;
2686 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + SIZE * ( shift + combined );
2688 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FFF_triplet[ irrep ][ elem ] + beta ); }
2697 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2698 const int SIZE = size_E[ irrep ];
2701 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2702 const int NVIR_a = indices->
getNVIRT( irrep_a );
2703 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2705 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2706 const int NOCC_i = indices->
getNOCC( irrep_i );
2708 if ( irrep_i == irrep_j ){
2709 const int SIZE_aij = ( NVIR_a * NOCC_i * ( NOCC_i + 1 ) ) / 2;
2710 #pragma omp for schedule(static) 2711 for (
int combined = 0; combined < SIZE_aij; combined++ ){
2712 const int a = combined % NVIR_a;
2714 const int i = triangle_idx[ 0 ];
2715 const int j = triangle_idx[ 1 ];
2716 const double f_ii = fock->
get( irrep_i, i, i );
2717 const double f_jj = fock->
get( irrep_j, j, j );
2718 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2719 const double beta = - f_dot_1dm + f_aa - f_ii - f_jj;
2720 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE * ( shift + combined );
2722 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FEE[ irrep ][ elem ] + beta ); }
2726 if ( irrep_i < irrep_j ){
2727 const int SIZE_aij = NVIR_a * NOCC_i * indices->
getNOCC( irrep_j );
2728 #pragma omp for schedule(static) 2729 for (
int combined = 0; combined < SIZE_aij; combined++ ){
2730 const int a = combined % NVIR_a;
2731 const int temp = combined / NVIR_a;
2732 const int i = temp % NOCC_i;
2733 const int j = temp / NOCC_i;
2734 const double f_ii = fock->
get( irrep_i, i, i );
2735 const double f_jj = fock->
get( irrep_j, j, j );
2736 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2737 const double beta = - f_dot_1dm + f_aa - f_ii - f_jj;
2738 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] + SIZE * ( shift + combined );
2740 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FEE[ irrep ][ elem ] + beta ); }
2750 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2751 const int SIZE = size_E[ irrep ];
2754 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2755 const int NVIR_a = indices->
getNVIRT( irrep_a );
2756 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2758 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2759 const int NOCC_i = indices->
getNOCC( irrep_i );
2761 if ( irrep_i == irrep_j ){
2762 const int SIZE_aij = ( NVIR_a * NOCC_i * ( NOCC_i - 1 ) ) / 2;
2763 #pragma omp for schedule(static) 2764 for (
int combined = 0; combined < SIZE_aij; combined++ ){
2765 const int a = combined % NVIR_a;
2767 const int i = triangle_idx[ 0 ];
2768 const int j = triangle_idx[ 1 ];
2769 const double f_ii = fock->
get( irrep_i, i, i );
2770 const double f_jj = fock->
get( irrep_j, j, j );
2771 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2772 const double beta = - f_dot_1dm + f_aa - f_ii - f_jj;
2773 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE * ( shift + combined );
2775 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 6 * ( FEE[ irrep ][ elem ] + beta ); }
2779 if ( irrep_i < irrep_j ){
2780 const int SIZE_aij = NVIR_a * NOCC_i * indices->
getNOCC( irrep_j );
2781 #pragma omp for schedule(static) 2782 for (
int combined = 0; combined < SIZE_aij; combined++ ){
2783 const int a = combined % NVIR_a;
2784 const int temp = combined / NVIR_a;
2785 const int i = temp % NOCC_i;
2786 const int j = temp / NOCC_i;
2787 const double f_ii = fock->
get( irrep_i, i, i );
2788 const double f_jj = fock->
get( irrep_j, j, j );
2789 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2790 const double beta = - f_dot_1dm + f_aa - f_ii - f_jj;
2791 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] + SIZE * ( shift + combined );
2793 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 6 * ( FEE[ irrep ][ elem ] + beta ); }
2803 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2804 const int SIZE = size_G[ irrep ];
2807 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2808 const int NOCC_i = indices->
getNOCC( irrep_i );
2810 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2811 const int NVIR_a = indices->
getNVIRT( irrep_a );
2812 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2814 if ( irrep_a == irrep_b ){
2815 const int SIZE_iab = ( NOCC_i * NVIR_a * ( NVIR_a + 1 ) ) / 2;
2816 #pragma omp for schedule(static) 2817 for (
int combined = 0; combined < SIZE_iab; combined++ ){
2818 const int i = combined % NOCC_i;
2820 const int a = triangle_idx[ 0 ];
2821 const int b = triangle_idx[ 1 ];
2822 const double f_ii = fock->
get( irrep_i, i, i );
2823 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2824 const double f_bb = fock->
get( irrep_b, N_OA_a + b, N_OA_a + b );
2825 const double beta = - f_dot_1dm + f_aa + f_bb - f_ii;
2826 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE * ( shift + combined );
2828 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FGG[ irrep ][ elem ] + beta ); }
2832 if ( irrep_a < irrep_b ){
2833 const int SIZE_iab = NOCC_i * NVIR_a * indices->
getNVIRT( irrep_b );
2834 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
2835 #pragma omp for schedule(static) 2836 for (
int combined = 0; combined < SIZE_iab; combined++ ){
2837 const int i = combined % NOCC_i;
2838 const int temp = combined / NOCC_i;
2839 const int a = temp % NVIR_a;
2840 const int b = temp / NVIR_a;
2841 const double f_ii = fock->
get( irrep_i, i, i );
2842 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2843 const double f_bb = fock->
get( irrep_b, N_OA_b + b, N_OA_b + b );
2844 const double beta = - f_dot_1dm + f_aa + f_bb - f_ii;
2845 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] + SIZE * ( shift + combined );
2847 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 2 * ( FGG[ irrep ][ elem ] + beta ); }
2857 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
2858 const int SIZE = size_G[ irrep ];
2861 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2862 const int NOCC_i = indices->
getNOCC( irrep_i );
2864 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2865 const int NVIR_a = indices->
getNVIRT( irrep_a );
2866 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2868 if ( irrep_a == irrep_b ){
2869 const int SIZE_iab = ( NOCC_i * NVIR_a * ( NVIR_a - 1 ) ) / 2;
2870 #pragma omp for schedule(static) 2871 for (
int combined = 0; combined < SIZE_iab; combined++ ){
2872 const int i = combined % NOCC_i;
2874 const int a = triangle_idx[ 0 ];
2875 const int b = triangle_idx[ 1 ];
2876 const double f_ii = fock->
get( irrep_i, i, i );
2877 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2878 const double f_bb = fock->
get( irrep_b, N_OA_a + b, N_OA_a + b );
2879 const double beta = - f_dot_1dm + f_aa + f_bb - f_ii;
2880 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE * ( shift + combined );
2882 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 6 * ( FGG[ irrep ][ elem ] + beta ); }
2886 if ( irrep_a < irrep_b ){
2887 const int SIZE_iab = NOCC_i * NVIR_a * indices->
getNVIRT( irrep_b );
2888 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
2889 #pragma omp for schedule(static) 2890 for (
int combined = 0; combined < SIZE_iab; combined++ ){
2891 const int i = combined % NOCC_i;
2892 const int temp = combined / NOCC_i;
2893 const int a = temp % NVIR_a;
2894 const int b = temp / NVIR_a;
2895 const double f_ii = fock->
get( irrep_i, i, i );
2896 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2897 const double f_bb = fock->
get( irrep_b, N_OA_b + b, N_OA_b + b );
2898 const double beta = - f_dot_1dm + f_aa + f_bb - f_ii;
2899 double * target = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] + SIZE * ( shift + combined );
2901 for (
int elem = 0; elem < SIZE; elem++ ){ target[ elem ] = 6 * ( FGG[ irrep ][ elem ] + beta ); }
2913 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
2914 const int NOCC_ij = indices->
getNOCC( irrep_ij );
2915 const int size_ij = ( NOCC_ij * ( NOCC_ij + 1 ) ) / 2;
2916 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
2917 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
2918 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
2919 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift;
2920 const int size_ijab = ( size_ij * NVIR_ab * ( NVIR_ab + 1 ) ) / 2;
2921 #pragma omp for schedule(static) 2922 for (
int combined = 0; combined < size_ijab; combined++ ){
2924 const int i = triangle_idx[ 0 ];
2925 const int j = triangle_idx[ 1 ];
2927 const int a = triangle_idx[ 0 ];
2928 const int b = triangle_idx[ 1 ];
2929 const double f_ii = fock->
get( irrep_ij, i, i );
2930 const double f_jj = fock->
get( irrep_ij, j, j );
2931 const double f_aa = fock->
get( irrep_ab, N_OA_ab + a, N_OA_ab + a );
2932 const double f_bb = fock->
get( irrep_ab, N_OA_ab + b, N_OA_ab + b );
2933 const double term = f_aa + f_bb - f_ii - f_jj;
2934 target[ combined ] = 4 * term;
2942 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
2943 const int NOCC_ij = indices->
getNOCC( irrep_ij );
2944 const int size_ij = ( NOCC_ij * ( NOCC_ij - 1 ) ) / 2;
2945 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
2946 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
2947 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
2948 double * target = result + jump[ num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift;
2949 const int size_ijab = ( size_ij * NVIR_ab * ( NVIR_ab - 1 ) ) / 2;
2950 #pragma omp for schedule(static) 2951 for (
int combined = 0; combined < size_ijab; combined++ ){
2953 const int i = triangle_idx[ 0 ];
2954 const int j = triangle_idx[ 1 ];
2956 const int a = triangle_idx[ 0 ];
2957 const int b = triangle_idx[ 1 ];
2958 const double f_ii = fock->
get( irrep_ij, i, i );
2959 const double f_jj = fock->
get( irrep_ij, j, j );
2960 const double f_aa = fock->
get( irrep_ab, N_OA_ab + a, N_OA_ab + a );
2961 const double f_bb = fock->
get( irrep_ab, N_OA_ab + b, N_OA_ab + b );
2962 const double term = f_aa + f_bb - f_ii - f_jj;
2963 target[ combined ] = 12 * term;
2969 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
2971 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
2973 if ( irrep_i < irrep_j ){
2974 const int NOCC_i = indices->
getNOCC( irrep_i );
2975 const int NOCC_j = indices->
getNOCC( irrep_j );
2976 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
2978 if ( irrep_a < irrep_b ){
2979 const int NVIR_a = indices->
getNVIRT( irrep_a );
2980 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
2981 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
2982 double * target_singlet = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] + shift;
2983 double * target_triplet = result + jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] + shift;
2984 const int size_ijab = NOCC_i * NOCC_j * NVIR_a * indices->
getNVIRT( irrep_b );
2985 #pragma omp for schedule(static) 2986 for (
int combined = 0; combined < size_ijab; combined++ ){
2987 const int i = combined % NOCC_i;
2988 const int temp1 = combined / NOCC_i;
2989 const int j = temp1 % NOCC_j;
2990 const int temp2 = temp1 / NOCC_j;
2991 const int a = temp2 % NVIR_a;
2992 const int b = temp2 / NVIR_a;
2993 const double f_ii = fock->
get( irrep_i, i, i );
2994 const double f_jj = fock->
get( irrep_j, j, j );
2995 const double f_aa = fock->
get( irrep_a, N_OA_a + a, N_OA_a + a );
2996 const double f_bb = fock->
get( irrep_b, N_OA_b + b, N_OA_b + b );
2997 const double term = f_aa + f_bb - f_ii - f_jj;
2998 target_singlet[ combined ] = 4 * term;
2999 target_triplet[ combined ] = 12 * term;
3066 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3067 const int NORB = indices->
getNORB( irrep );
3068 for (
int row = 0; row < NORB; row++ ){
3069 for (
int col = row; col < NORB; col++ ){
3070 double value = oei->
get( irrep, row, col );
3071 for (
int irrep_occ = 0; irrep_occ < num_irreps; irrep_occ++ ){
3072 const int NOCC = indices->
getNOCC( irrep_occ );
3073 for (
int occ = 0; occ < NOCC; occ++ ){
3074 value += ( 2 * integrals->
FourIndexAPI( irrep, irrep_occ, irrep, irrep_occ, row, occ, col, occ )
3075 - integrals->
FourIndexAPI( irrep, irrep_occ, irrep_occ, irrep, row, occ, occ, col ) );
3078 MAT->
set( irrep, row, col, value );
3079 MAT->
set( irrep, col, row, value );
3086 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3087 const int NORB = indices->
getNORB( irrep );
3088 for (
int row = 0; row < NORB; row++ ){
3089 for (
int col = row; col < NORB; col++ ){
3091 for (
int irrep_act = 0; irrep_act < num_irreps; irrep_act++ ){
3092 const int NOCC = indices->
getNOCC( irrep_act );
3093 const int NACT = indices->
getNDMRG( irrep_act );
3094 for (
int act = 0; act < NACT; act++ ){
3095 value += integrals->
FourIndexAPI( irrep, irrep_act, irrep_act, irrep, row, NOCC + act, NOCC + act, col );
3098 MAT2->
set( irrep, row, col, value );
3099 MAT2->
set( irrep, col, row, value );
3104 vector_rhs =
new double[ jump[ num_irreps * CHEMPS2_CASPT2_NUM_CASES ] ];
3106 #pragma omp parallel 3108 const int max_size = get_maxsize();
3109 double * workspace =
new double[ max_size ];
3112 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3113 if ( size_A[ irrep ] > 0 ){
3114 const int NOCC = indices->
getNOCC( irrep );
3115 const int NACT = indices->
getNDMRG( irrep );
3117 #pragma omp for schedule(static) 3118 for (
int count_i = 0; count_i < NOCC; count_i++ ){
3120 double * target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_A ] + size_A[ irrep ] * count_i;
3126 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3127 const int occ_x = indices->
getNOCC( irrep_x );
3128 const int num_x = indices->
getNDMRG( irrep_x );
3130 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
3132 const int occ_y = indices->
getNOCC( irrep_y );
3133 const int occ_z = indices->
getNOCC( irrep_z );
3134 const int num_y = indices->
getNDMRG( irrep_y );
3135 const int num_z = indices->
getNDMRG( irrep_z );
3140 for (
int z = 0; z < num_z; z++ ){
3141 for (
int y = 0; y < num_y; y++ ){
3142 for (
int x = 0; x < num_x; x++ ){
3143 const double ix_zy = integrals->
get_coulomb( irrep, irrep_x, irrep_z, irrep_y, count_i, occ_x + x, occ_z + z, occ_y + y );
3144 workspace[ jump_xyz + x + num_x * ( y + num_y * z ) ] = ix_zy;
3150 for (
int z = 0; z < num_z; z++ ){
3151 for (
int y = 0; y < num_y; y++ ){
3152 for (
int x = 0; x < num_x; x++ ){
3154 for (
int w = 0; w < NACT; w++ ){
3155 value += MAT->
get( irrep, count_i, NOCC + w ) * two_rdm[ d_x + x + LAS * ( d_y + y + LAS * ( d_w + w + LAS * ( d_z + z ) ) ) ];
3157 target[ jump_xyz + x + num_x * ( y + num_y * z ) ] = - value;
3163 if ( irrep_x == irrep ){
3164 for (
int z = 0; z < num_z; z++ ){
3165 for (
int y = 0; y < num_y; y++ ){
3166 for (
int x = 0; x < num_x; x++ ){
3167 target[ jump_xyz + x + num_x * ( y + num_y * z ) ] += 2 * MAT->
get( irrep, count_i, occ_x + x ) * one_rdm[ d_y + y + LAS * ( d_z + z ) ];
3174 if ( irrep_y == irrep ){
3175 for (
int z = 0; z < num_z; z++ ){
3176 for (
int y = 0; y < num_y; y++ ){
3177 for (
int x = 0; x < num_x; x++ ){
3178 target[ jump_xyz + x + num_x * ( y + num_y * z ) ] -= MAT->
get( irrep, count_i, occ_y + y ) * one_rdm[ d_x + x + LAS * ( d_z + z ) ];
3183 jump_xyz += num_x * num_y * num_z;
3186 assert( jump_xyz == size_A[ irrep ] );
3192 dgemv_( ¬rans, &jump_xyz, &jump_xyz, &one, SAA[ irrep ], &jump_xyz, workspace, &inc1, &one, target, &inc1 );
3194 assert( NOCC * size_A[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_A ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_A ] );
3199 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3200 if ( size_C[ irrep ] > 0 ){
3201 const int NOCC = indices->
getNOCC( irrep );
3202 const int NVIR = indices->
getNVIRT( irrep );
3203 const int NACT = indices->
getNDMRG( irrep );
3204 const int N_OA = NOCC + NACT;
3206 #pragma omp for schedule(static) 3207 for (
int count_a = 0; count_a < NVIR; count_a++ ){
3209 double * target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_C ] + size_C[ irrep ] * count_a;
3214 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3215 const int occ_x = indices->
getNOCC( irrep_x );
3216 const int num_x = indices->
getNDMRG( irrep_x );
3218 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
3220 const int occ_y = indices->
getNOCC( irrep_y );
3221 const int occ_z = indices->
getNOCC( irrep_z );
3222 const int num_y = indices->
getNDMRG( irrep_y );
3223 const int num_z = indices->
getNDMRG( irrep_z );
3228 for (
int z = 0; z < num_z; z++ ){
3229 for (
int y = 0; y < num_y; y++ ){
3230 for (
int x = 0; x < num_x; x++ ){
3231 const double zy_xa = integrals->
get_coulomb( irrep_z, irrep_y, irrep_x, irrep, occ_z + z, occ_y + y, occ_x + x, N_OA + count_a );
3232 workspace[ jump_xyz + x + num_x * ( y + num_y * z ) ] = zy_xa;
3238 for (
int z = 0; z < num_z; z++ ){
3239 for (
int y = 0; y < num_y; y++ ){
3240 for (
int x = 0; x < num_x; x++ ){
3242 for (
int w = 0; w < NACT; w++ ){
3243 value += ( ( MAT->
get( irrep, NOCC + w, N_OA + count_a ) - MAT2->
get( irrep, NOCC + w, N_OA + count_a ) )
3244 * two_rdm[ d_w + w + LAS * ( d_y + y + LAS * ( d_x + x + LAS * ( d_z + z ) ) ) ] );
3246 target[ jump_xyz + x + num_x * ( y + num_y * z ) ] = value;
3252 if (( irrep_z == irrep ) && ( irrep_x == irrep_y )){
3253 for (
int z = 0; z < num_z; z++ ){
3255 for (
int w = 0; w < NACT; w++ ){
3256 value += ( ( MAT->
get( irrep, NOCC + w, N_OA + count_a ) - MAT2->
get( irrep, NOCC + w, N_OA + count_a ) )
3257 * one_rdm[ d_w + w + LAS * ( d_z + z ) ] );
3259 for (
int xy = 0; xy < num_x; xy++ ){
3260 target[ jump_xyz + xy + num_x * ( xy + num_x * z ) ] += value;
3264 jump_xyz += num_x * num_y * num_z;
3267 assert( jump_xyz == size_C[ irrep ] );
3273 dgemv_( ¬rans, &jump_xyz, &jump_xyz, &one, SCC[ irrep ], &jump_xyz, workspace, &inc1, &one, target, &inc1 );
3275 assert( NVIR * size_C[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_C ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_C ] );
3284 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3285 if ( size_D[ irrep ] > 0 ){
3287 const int D2JUMP = size_D[ irrep ] / 2;
3288 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
3290 assert( shift == shift_D_nonactive( indices, irrep_i, irrep_a ) );
3291 const int NOCC_i = indices->
getNOCC( irrep_i );
3292 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
3293 const int NVIR_a = indices->
getNVIRT( irrep_a );
3294 const int loopsize = NOCC_i * NVIR_a;
3295 #pragma omp for schedule(static) 3296 for (
int combined = 0; combined < loopsize; combined++ ){
3298 const int count_i = combined % NOCC_i;
3299 const int count_a = combined / NOCC_i;
3300 double * target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_D ] + size_D[ irrep ] * ( shift + combined );
3301 const double MAT_ia = ( ( irrep_i == irrep_a ) ? MAT->
get( irrep_i, count_i, N_OA_a + count_a ) : 0.0 );
3312 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3314 const int occ_x = indices->
getNOCC( irrep_x );
3315 const int occ_y = indices->
getNOCC( irrep_y );
3316 const int num_x = indices->
getNDMRG( irrep_x );
3317 const int num_y = indices->
getNDMRG( irrep_y );
3320 for (
int y = 0; y < num_y; y++ ){
3321 for (
int x = 0; x < num_x; x++ ){
3322 const double ia_yx = integrals->
get_coulomb( irrep_y, irrep_x, irrep_i, irrep_a, occ_y + y, occ_x + x, count_i, N_OA_a + count_a );
3323 workspace[ jump_xy + x + num_x * y ] = ia_yx;
3328 for (
int y = 0; y < num_y; y++ ){
3329 for (
int x = 0; x < num_x; x++ ){
3330 const double ix_ya = integrals->
get_coulomb( irrep_i, irrep_x, irrep_y, irrep_a, count_i, occ_x + x, occ_y + y, N_OA_a + count_a );
3331 workspace[ D2JUMP + jump_xy + x + num_x * y ] = ix_ya;
3339 for (
int y = 0; y < num_y; y++ ){
3340 for (
int x = 0; x < num_x; x++ ){
3341 const double value = MAT_ia * one_rdm[ d_xy + x + LAS * ( d_xy + y ) ];
3342 target[ jump_xy + x + num_x * y ] = 2 * value;
3343 target[ D2JUMP + jump_xy + x + num_x * y ] = - value;
3347 for (
int y = 0; y < num_y; y++ ){
3348 for (
int x = 0; x < num_x; x++ ){
3349 target[ jump_xy + x + num_x * y ] = 0.0;
3350 target[ D2JUMP + jump_xy + x + num_x * y ] = 0.0;
3354 jump_xy += num_x * num_y;
3356 jump_xy = 2 * jump_xy;
3357 assert( jump_xy == size_D[ irrep ] );
3363 dgemv_( ¬rans, &jump_xy, &jump_xy, &one, SDD[ irrep ], &jump_xy, workspace, &inc1, &one, target, &inc1 );
3367 assert( shift * size_D[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_D ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_D ] );
3374 const double SQRT_0p5 = sqrt( 0.5 );
3375 int triangle_idx[] = { 0, 0 };
3378 if ( size_B_singlet[ 0 ] > 0 ){
3380 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
3381 assert( shift == shift_B_nonactive( indices, irrep_ij, irrep_ij, +1 ) );
3382 const int NOCC_ij = indices->
getNOCC( irrep_ij );
3383 const int loopsize = ( NOCC_ij * ( NOCC_ij + 1 ) ) / 2;
3384 #pragma omp for schedule(static) 3385 for (
int combined = 0; combined < loopsize; combined++ ){
3388 const int i = triangle_idx[ 0 ];
3389 const int j = triangle_idx[ 1 ];
3390 assert( combined == i + ( j * ( j + 1 ) ) / 2 );
3394 for (
int irrep_xy = 0; irrep_xy < num_irreps; irrep_xy++ ){
3395 const int occ_xy = indices->
getNOCC( irrep_xy );
3396 const int num_xy = indices->
getNDMRG( irrep_xy );
3398 for (
int x = 0; x < num_xy; x++ ){
3399 for (
int y = x; y < num_xy; y++ ){
3400 const double ix_jy = integrals->
get_coulomb( irrep_ij, irrep_xy, irrep_ij, irrep_xy, i, occ_xy + x, j, occ_xy + y );
3401 const double iy_jx = integrals->
get_coulomb( irrep_ij, irrep_xy, irrep_ij, irrep_xy, i, occ_xy + y, j, occ_xy + x );
3402 workspace[ jump_xy + x + (y*(y+1))/2 ] = ( ix_jy + iy_jx ) * (( x==y ) ? 0.5 : 1.0 );
3405 jump_xy += ( num_xy * ( num_xy + 1 ) ) / 2;
3407 assert( jump_xy == size_B_singlet[ 0 ] );
3412 double alpha = (( i == j ) ? SQRT_0p5 : 1.0 );
3414 double * target = vector_rhs + jump[ num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + size_B_singlet[ 0 ] * ( shift + combined );
3415 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SBB_singlet[ 0 ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3419 assert( shift * size_B_singlet[ 0 ] == jump[ 1 + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] - jump[ num_irreps * CHEMPS2_CASPT2_B_SINGLET ] );
3421 if ( size_B_triplet[ 0 ] > 0 ){
3423 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
3424 assert( shift == shift_B_nonactive( indices, irrep_ij, irrep_ij, -1 ) );
3425 const int NOCC_ij = indices->
getNOCC( irrep_ij );
3426 const int loopsize = ( NOCC_ij * ( NOCC_ij - 1 ) ) / 2;
3427 #pragma omp for schedule(static) 3428 for (
int combined = 0; combined < loopsize; combined++ ){
3431 const int i = triangle_idx[ 0 ];
3432 const int j = triangle_idx[ 1 ];
3433 assert( combined == i + ( j * ( j - 1 ) ) / 2 );
3437 for (
int irrep_xy = 0; irrep_xy < num_irreps; irrep_xy++ ){
3438 const int occ_xy = indices->
getNOCC( irrep_xy );
3439 const int num_xy = indices->
getNDMRG( irrep_xy );
3441 for (
int x = 0; x < num_xy; x++ ){
3442 for (
int y = x+1; y < num_xy; y++ ){
3443 const double ix_jy = integrals->
get_coulomb( irrep_ij, irrep_xy, irrep_ij, irrep_xy, i, occ_xy + x, j, occ_xy + y );
3444 const double iy_jx = integrals->
get_coulomb( irrep_ij, irrep_xy, irrep_ij, irrep_xy, i, occ_xy + y, j, occ_xy + x );
3445 workspace[ jump_xy + x + (y*(y-1))/2 ] = ix_jy - iy_jx;
3448 jump_xy += ( num_xy * ( num_xy - 1 ) ) / 2;
3450 assert( jump_xy == size_B_triplet[ 0 ] );
3457 double * target = vector_rhs + jump[ num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + size_B_triplet[ 0 ] * ( shift + combined );
3458 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SBB_triplet[ 0 ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3462 assert( shift * size_B_triplet[ 0 ] == jump[ 1 + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] - jump[ num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] );
3464 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
3465 assert( size_B_singlet[ irrep ] == size_B_triplet[ irrep ] );
3466 if ( size_B_singlet[ irrep ] > 0 ){
3468 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
3470 if ( irrep_i < irrep_j ){
3471 assert( shift == shift_B_nonactive( indices, irrep_i, irrep_j, 0 ) );
3472 const int NOCC_i = indices->
getNOCC( irrep_i );
3473 const int NOCC_j = indices->
getNOCC( irrep_j );
3474 const int loopsize = NOCC_i * NOCC_j;
3475 #pragma omp for schedule(static) 3476 for (
int combined = 0; combined < loopsize; combined++ ){
3478 const int i = combined % NOCC_i;
3479 const int j = combined / NOCC_i;
3483 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3485 if ( irrep_x < irrep_y ){
3486 const int occ_x = indices->
getNOCC( irrep_x );
3487 const int occ_y = indices->
getNOCC( irrep_y );
3488 const int num_x = indices->
getNDMRG( irrep_x );
3489 const int num_y = indices->
getNDMRG( irrep_y );
3491 for (
int y = 0; y < num_y; y++ ){
3492 for (
int x = 0; x < num_x; x++ ){
3493 const double ix_jy = integrals->
get_coulomb( irrep_i, irrep_x, irrep_j, irrep_y, i, occ_x + x, j, occ_y + y );
3494 const double iy_jx = integrals->
get_coulomb( irrep_i, irrep_y, irrep_j, irrep_x, i, occ_y + y, j, occ_x + x );
3495 workspace[ jump_xy + x + num_x * y ] = ix_jy + iy_jx;
3498 jump_xy += num_x * num_y;
3501 assert( jump_xy == size_B_singlet[ irrep ] );
3508 double * target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] + size_B_singlet[ irrep ] * ( shift + combined );
3509 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SBB_singlet[ irrep ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3513 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3515 if ( irrep_x < irrep_y ){
3516 const int occ_x = indices->
getNOCC( irrep_x );
3517 const int occ_y = indices->
getNOCC( irrep_y );
3518 const int num_x = indices->
getNDMRG( irrep_x );
3519 const int num_y = indices->
getNDMRG( irrep_y );
3521 for (
int y = 0; y < num_y; y++ ){
3522 for (
int x = 0; x < num_x; x++ ){
3523 const double ix_jy = integrals->
get_coulomb( irrep_i, irrep_x, irrep_j, irrep_y, i, occ_x + x, j, occ_y + y );
3524 const double iy_jx = integrals->
get_coulomb( irrep_i, irrep_y, irrep_j, irrep_x, i, occ_y + y, j, occ_x + x );
3525 workspace[ jump_xy + x + num_x * y ] = ix_jy - iy_jx;
3528 jump_xy += num_x * num_y;
3531 assert( jump_xy == size_B_triplet[ irrep ] );
3534 target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] + size_B_triplet[ irrep ] * ( shift + combined );
3535 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SBB_triplet[ irrep ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3540 assert( shift * size_B_singlet[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_SINGLET ] );
3541 assert( shift * size_B_triplet[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_B_TRIPLET ] );
3546 if ( size_F_singlet[ 0 ] > 0 ){
3548 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
3549 assert( shift == shift_F_nonactive( indices, irrep_ab, irrep_ab, +1 ) );
3550 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
3551 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
3552 const int loopsize = ( NVIR_ab * ( NVIR_ab + 1 ) ) / 2;
3553 #pragma omp for schedule(static) 3554 for (
int combined = 0; combined < loopsize; combined++ ){
3557 const int a = triangle_idx[ 0 ];
3558 const int b = triangle_idx[ 1 ];
3559 assert( combined == a + ( b * ( b + 1 ) ) / 2 );
3563 for (
int irrep_xy = 0; irrep_xy < num_irreps; irrep_xy++ ){
3564 const int occ_xy = indices->
getNOCC( irrep_xy );
3565 const int num_xy = indices->
getNDMRG( irrep_xy );
3567 for (
int x = 0; x < num_xy; x++ ){
3568 for (
int y = x; y < num_xy; y++ ){
3569 const double ax_by = integrals->
get_exchange( irrep_xy, irrep_xy, irrep_ab, irrep_ab, occ_xy + x, occ_xy + y, N_OA_ab + a, N_OA_ab + b );
3570 const double ay_bx = integrals->
get_exchange( irrep_xy, irrep_xy, irrep_ab, irrep_ab, occ_xy + y, occ_xy + x, N_OA_ab + a, N_OA_ab + b );
3571 workspace[ jump_xy + x + (y*(y+1))/2 ] = ( ax_by + ay_bx ) * (( x==y ) ? 0.5 : 1.0 );
3574 jump_xy += ( num_xy * ( num_xy + 1 ) ) / 2;
3576 assert( jump_xy == size_F_singlet[ 0 ] );
3581 double alpha = (( a == b ) ? SQRT_0p5 : 1.0 );
3583 double * target = vector_rhs + jump[ 0 + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + size_F_singlet[ 0 ] * ( shift + combined );
3584 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SFF_singlet[ 0 ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3588 assert( shift * size_F_singlet[ 0 ] == jump[ 1 + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] - jump[ num_irreps * CHEMPS2_CASPT2_F_SINGLET ] );
3590 if ( size_F_triplet[ 0 ] > 0 ){
3592 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
3593 assert( shift == shift_F_nonactive( indices, irrep_ab, irrep_ab, -1 ) );
3594 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
3595 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
3596 const int loopsize = ( NVIR_ab * ( NVIR_ab - 1 ) ) / 2;
3597 #pragma omp for schedule(static) 3598 for (
int combined = 0; combined < loopsize; combined++ ){
3601 const int a = triangle_idx[ 0 ];
3602 const int b = triangle_idx[ 1 ];
3603 assert( combined == a + ( b * ( b - 1 ) ) / 2 );
3607 for (
int irrep_xy = 0; irrep_xy < num_irreps; irrep_xy++ ){
3608 const int occ_xy = indices->
getNOCC( irrep_xy );
3609 const int num_xy = indices->
getNDMRG( irrep_xy );
3611 for (
int x = 0; x < num_xy; x++ ){
3612 for (
int y = x+1; y < num_xy; y++ ){
3613 const double ax_by = integrals->
get_exchange( irrep_xy, irrep_xy, irrep_ab, irrep_ab, occ_xy + x, occ_xy + y, N_OA_ab + a, N_OA_ab + b );
3614 const double ay_bx = integrals->
get_exchange( irrep_xy, irrep_xy, irrep_ab, irrep_ab, occ_xy + y, occ_xy + x, N_OA_ab + a, N_OA_ab + b );
3615 workspace[ jump_xy + x + (y*(y-1))/2 ] = ax_by - ay_bx;
3618 jump_xy += ( num_xy * ( num_xy - 1 ) ) / 2;
3620 assert( jump_xy == size_F_triplet[ 0 ] );
3627 double * target = vector_rhs + jump[ 0 + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + size_F_triplet[ 0 ] * ( shift + combined );
3628 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SFF_triplet[ 0 ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3632 assert( shift * size_F_triplet[ 0 ] == jump[ 1 + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] - jump[ num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] );
3634 for (
int irrep = 1; irrep < num_irreps; irrep++ ){
3635 assert( size_F_singlet[ irrep ] == size_F_triplet[ irrep ] );
3636 if ( size_F_singlet[ irrep ] > 0 ){
3638 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
3640 if ( irrep_a < irrep_b ){
3641 assert( shift == shift_F_nonactive( indices, irrep_a, irrep_b, 0 ) );
3642 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
3643 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
3644 const int NVIR_a = indices->
getNVIRT( irrep_a );
3645 const int NVIR_b = indices->
getNVIRT( irrep_b );
3646 const int loopsize = NVIR_a * NVIR_b;
3647 #pragma omp for schedule(static) 3648 for (
int combined = 0; combined < loopsize; combined++ ){
3650 const int a = combined % NVIR_a;
3651 const int b = combined / NVIR_a;
3655 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3657 if ( irrep_x < irrep_y ){
3658 const int occ_x = indices->
getNOCC( irrep_x );
3659 const int occ_y = indices->
getNOCC( irrep_y );
3660 const int num_x = indices->
getNDMRG( irrep_x );
3661 const int num_y = indices->
getNDMRG( irrep_y );
3663 for (
int y = 0; y < num_y; y++ ){
3664 for (
int x = 0; x < num_x; x++ ){
3665 const double ax_by = integrals->
get_exchange( irrep_x, irrep_y, irrep_a, irrep_b, occ_x + x, occ_y + y, N_OA_a + a, N_OA_b + b );
3666 const double ay_bx = integrals->
get_exchange( irrep_y, irrep_x, irrep_a, irrep_b, occ_y + y, occ_x + x, N_OA_a + a, N_OA_b + b );
3667 workspace[ jump_xy + x + num_x * y ] = ax_by + ay_bx;
3670 jump_xy += num_x * num_y;
3673 assert( jump_xy == size_F_singlet[ irrep ] );
3680 double * target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] + size_F_singlet[ irrep ] * ( shift + combined );
3681 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SFF_singlet[ irrep ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3685 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
3687 if ( irrep_x < irrep_y ){
3688 const int occ_x = indices->
getNOCC( irrep_x );
3689 const int occ_y = indices->
getNOCC( irrep_y );
3690 const int num_x = indices->
getNDMRG( irrep_x );
3691 const int num_y = indices->
getNDMRG( irrep_y );
3693 for (
int y = 0; y < num_y; y++ ){
3694 for (
int x = 0; x < num_x; x++ ){
3695 const double ax_by = integrals->
get_exchange( irrep_x, irrep_y, irrep_a, irrep_b, occ_x + x, occ_y + y, N_OA_a + a, N_OA_b + b );
3696 const double ay_bx = integrals->
get_exchange( irrep_y, irrep_x, irrep_a, irrep_b, occ_y + y, occ_x + x, N_OA_a + a, N_OA_b + b );
3697 workspace[ jump_xy + x + num_x * y ] = ax_by - ay_bx;
3700 jump_xy += num_x * num_y;
3703 assert( jump_xy == size_F_triplet[ irrep ] );
3706 target = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] + size_F_triplet[ irrep ] * ( shift + combined );
3707 dgemv_( ¬rans, &jump_xy, &jump_xy, &alpha, SFF_triplet[ irrep ], &jump_xy, workspace, &inc1, &
set, target, &inc1 );
3712 assert( shift * size_F_singlet[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_SINGLET ] );
3713 assert( shift * size_F_triplet[ irrep ] == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_F_TRIPLET ] );
3716 delete [] workspace;
3719 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3720 const int occ_t = indices->
getNOCC( irrep );
3721 const int num_t = indices->
getNDMRG( irrep );
3723 double * target_singlet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ];
3724 double * target_triplet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ];
3725 int shift_singlet = 0;
3726 int shift_triplet = 0;
3727 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
3728 const int NVIR_a = indices->
getNVIRT( irrep_a );
3729 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
3731 if ( irrep_occ == 0 ){
3732 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
3733 assert( shift_singlet == shift_E_nonactive( indices, irrep_a, irrep_ij, irrep_ij, +1 ) );
3734 assert( shift_triplet == shift_E_nonactive( indices, irrep_a, irrep_ij, irrep_ij, -1 ) );
3735 const int NOCC_ij = indices->
getNOCC( irrep_ij );
3736 const int loopsize_singlet = ( NVIR_a * NOCC_ij * ( NOCC_ij + 1 ) ) / 2;
3737 const int loopsize_triplet = ( NVIR_a * NOCC_ij * ( NOCC_ij - 1 ) ) / 2;
3738 #pragma omp for schedule(static) 3739 for (
int combined_singlet = 0; combined_singlet < loopsize_singlet; combined_singlet++ ){
3741 const int a = combined_singlet % NVIR_a;
3742 const int remainder = combined_singlet / NVIR_a;
3744 const int i = triangle_idx[ 0 ];
3745 const int j = triangle_idx[ 1 ];
3746 const int combined_triplet = a + NVIR_a * ( i + ( j * ( j - 1 ) ) / 2 );
3747 const double ij_factor = (( i == j ) ? SQRT_0p5 : 1.0 );
3749 const int count_aij_singlet = shift_singlet + combined_singlet;
3750 const int count_aij_triplet = shift_triplet + combined_triplet;
3751 for (
int t = 0; t < num_t; t++ ){
3752 double value_singlet = 0.0;
3753 double value_triplet = 0.0;
3754 for (
int w = 0; w < num_t; w++ ){
3755 const double SEE_wt = SEE[ irrep ][ w + num_t * t ];
3756 const double aj_wi = integrals->
get_coulomb( irrep_ij, irrep, irrep_ij, irrep_a, i, occ_t + w, j, N_OA_a + a );
3757 const double ai_wj = integrals->
get_coulomb( irrep_ij, irrep, irrep_ij, irrep_a, j, occ_t + w, i, N_OA_a + a );
3758 value_singlet += SEE_wt * ( aj_wi + ai_wj );
3759 value_triplet += 3 * SEE_wt * ( aj_wi - ai_wj );
3761 target_singlet[ t + num_t * count_aij_singlet ] = value_singlet * ij_factor;
3762 if ( j > i ){ target_triplet[ t + num_t * count_aij_triplet ] = value_triplet; }
3765 shift_singlet += loopsize_singlet;
3766 shift_triplet += loopsize_triplet;
3769 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
3771 if ( irrep_i < irrep_j ){
3772 assert( shift_singlet == shift_E_nonactive( indices, irrep_a, irrep_i, irrep_j, +1 ) );
3773 assert( shift_triplet == shift_E_nonactive( indices, irrep_a, irrep_i, irrep_j, -1 ) );
3774 const int NOCC_i = indices->
getNOCC( irrep_i );
3775 const int NOCC_j = indices->
getNOCC( irrep_j );
3776 const int loopsize = NOCC_i * NOCC_j * NVIR_a;
3777 #pragma omp for schedule(static) 3778 for (
int combined = 0; combined < loopsize; combined++ ){
3780 const int a = combined % NVIR_a;
3781 const int remainder = combined / NVIR_a;
3782 const int i = remainder % NOCC_i;
3783 const int j = remainder / NOCC_i;
3785 const int count_aij_singlet = shift_singlet + combined;
3786 const int count_aij_triplet = shift_triplet + combined;
3787 for (
int t = 0; t < num_t; t++ ){
3788 double value_singlet = 0.0;
3789 double value_triplet = 0.0;
3790 for (
int w = 0; w < num_t; w++ ){
3791 const double SEE_wt = SEE[ irrep ][ w + num_t * t ];
3792 const double aj_wi = integrals->
get_coulomb( irrep_i, irrep, irrep_j, irrep_a, i, occ_t + w, j, N_OA_a + a );
3793 const double ai_wj = integrals->
get_coulomb( irrep_j, irrep, irrep_i, irrep_a, j, occ_t + w, i, N_OA_a + a );
3794 value_singlet += SEE_wt * ( aj_wi + ai_wj );
3795 value_triplet += 3 * SEE_wt * ( aj_wi - ai_wj );
3797 target_singlet[ t + num_t * count_aij_singlet ] = value_singlet;
3798 target_triplet[ t + num_t * count_aij_triplet ] = value_triplet;
3801 shift_singlet += loopsize;
3802 shift_triplet += loopsize;
3807 assert( shift_singlet * num_t == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_SINGLET ] );
3808 assert( shift_triplet * num_t == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_E_TRIPLET ] );
3813 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3814 const int occ_t = indices->
getNOCC( irrep );
3815 const int num_t = indices->
getNDMRG( irrep );
3817 double * target_singlet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ];
3818 double * target_triplet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ];
3819 int shift_singlet = 0;
3820 int shift_triplet = 0;
3821 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
3822 const int NOCC_i = indices->
getNOCC( irrep_i );
3824 if ( irrep_virt == 0 ){
3825 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
3826 assert( shift_singlet == shift_G_nonactive( indices, irrep_i, irrep_ab, irrep_ab, +1 ) );
3827 assert( shift_triplet == shift_G_nonactive( indices, irrep_i, irrep_ab, irrep_ab, -1 ) );
3828 const int N_OA_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
3829 const int NVIR_ab = indices->
getNVIRT( irrep_ab );
3830 const int loopsize_singlet = ( NOCC_i * NVIR_ab * ( NVIR_ab + 1 ) ) / 2;
3831 const int loopsize_triplet = ( NOCC_i * NVIR_ab * ( NVIR_ab - 1 ) ) / 2;
3832 #pragma omp for schedule(static) 3833 for (
int combined_singlet = 0; combined_singlet < loopsize_singlet; combined_singlet++ ){
3835 const int i = combined_singlet % NOCC_i;
3836 const int remainder = combined_singlet / NOCC_i;
3838 const int a = triangle_idx[ 0 ];
3839 const int b = triangle_idx[ 1 ];
3840 const int combined_triplet = i + NOCC_i * ( a + ( b * ( b - 1 ) ) / 2 );
3841 const double ab_factor = (( a == b ) ? SQRT_0p5 : 1.0 );
3843 const int count_abi_singlet = shift_singlet + combined_singlet;
3844 const int count_abi_triplet = shift_triplet + combined_triplet;
3845 for (
int t = 0; t < num_t; t++ ){
3846 double value_singlet = 0.0;
3847 double value_triplet = 0.0;
3848 for (
int u = 0; u < num_t; u++ ){
3849 const double SGG_ut = SGG[ irrep ][ u + num_t * t ];
3850 const double ai_bu = integrals->
get_exchange( irrep_i, irrep, irrep_ab, irrep_ab, i, occ_t + u, N_OA_ab + a, N_OA_ab + b );
3851 const double bi_au = integrals->
get_exchange( irrep_i, irrep, irrep_ab, irrep_ab, i, occ_t + u, N_OA_ab + b, N_OA_ab + a );
3852 value_singlet += SGG_ut * ( ai_bu + bi_au );
3853 value_triplet += 3 * SGG_ut * ( ai_bu - bi_au );
3855 target_singlet[ t + num_t * count_abi_singlet ] = value_singlet * ab_factor;
3856 if ( b > a ){ target_triplet[ t + num_t * count_abi_triplet ] = value_triplet; }
3859 shift_singlet += loopsize_singlet;
3860 shift_triplet += loopsize_triplet;
3863 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
3865 if ( irrep_a < irrep_b ){
3866 assert( shift_singlet == shift_G_nonactive( indices, irrep_i, irrep_a, irrep_b, +1 ) );
3867 assert( shift_triplet == shift_G_nonactive( indices, irrep_i, irrep_a, irrep_b, -1 ) );
3868 const int N_OA_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
3869 const int N_OA_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
3870 const int NVIR_a = indices->
getNVIRT( irrep_a );
3871 const int NVIR_b = indices->
getNVIRT( irrep_b );
3872 const int loopsize = NOCC_i * NVIR_a * NVIR_b;
3873 #pragma omp for schedule(static) 3874 for (
int combined = 0; combined < loopsize; combined++ ){
3876 const int i = combined % NOCC_i;
3877 const int remainder = combined / NOCC_i;
3878 const int a = remainder % NVIR_a;
3879 const int b = remainder / NVIR_a;
3881 const int count_abi_singlet = shift_singlet + combined;
3882 const int count_abi_triplet = shift_triplet + combined;
3883 for (
int t = 0; t < num_t; t++ ){
3884 double value_singlet = 0.0;
3885 double value_triplet = 0.0;
3886 for (
int u = 0; u < num_t; u++ ){
3887 const double SGG_ut = SGG[ irrep ][ u + num_t * t ];
3888 const double ai_bu = integrals->
get_exchange( irrep_i, irrep, irrep_a, irrep_b, i, occ_t + u, N_OA_a + a, N_OA_b + b );
3889 const double bi_au = integrals->
get_exchange( irrep_i, irrep, irrep_b, irrep_a, i, occ_t + u, N_OA_b + b, N_OA_a + a );
3890 value_singlet += SGG_ut * ( ai_bu + bi_au );
3891 value_triplet += 3 * SGG_ut * ( ai_bu - bi_au );
3893 target_singlet[ t + num_t * count_abi_singlet ] = value_singlet;
3894 target_triplet[ t + num_t * count_abi_triplet ] = value_triplet;
3897 shift_singlet += loopsize;
3898 shift_triplet += loopsize;
3903 assert( shift_singlet * num_t == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_SINGLET ] );
3904 assert( shift_triplet * num_t == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_G_TRIPLET ] );
3909 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
3910 int shift_singlet = 0;
3911 int shift_triplet = 0;
3912 double * target_singlet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET ];
3913 double * target_triplet = vector_rhs + jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ];
3915 for (
int irrep_ij = 0; irrep_ij < num_irreps; irrep_ij++ ){
3916 const int nocc_ij = indices->
getNOCC( irrep_ij );
3917 const int linsize_ij_singlet = ( nocc_ij * ( nocc_ij + 1 ) ) / 2;
3918 const int linsize_ij_triplet = ( nocc_ij * ( nocc_ij - 1 ) ) / 2;
3919 for (
int irrep_ab = 0; irrep_ab < num_irreps; irrep_ab++ ){
3920 assert( shift_singlet == shift_H_nonactive( indices, irrep_ij, irrep_ij, irrep_ab, irrep_ab, +1 ) );
3921 assert( shift_triplet == shift_H_nonactive( indices, irrep_ij, irrep_ij, irrep_ab, irrep_ab, -1 ) );
3922 const int nvirt_ab = indices->
getNVIRT( irrep_ab );
3923 const int noa_ab = indices->
getNOCC( irrep_ab ) + indices->
getNDMRG( irrep_ab );
3924 const int linsize_ab_singlet = ( nvirt_ab * ( nvirt_ab + 1 ) ) / 2;
3925 const int linsize_ab_triplet = ( nvirt_ab * ( nvirt_ab - 1 ) ) / 2;
3926 #pragma omp for schedule(static) 3927 for (
int combined_ab_singlet = 0; combined_ab_singlet < linsize_ab_singlet; combined_ab_singlet++ ){
3930 const int a = triangle_idx[ 0 ];
3931 const int b = triangle_idx[ 1 ];
3932 const double ab_factor = (( a == b ) ? SQRT_0p5 : 1.0 );
3933 const int combined_ab_triplet = a + ( b * ( b - 1 ) ) / 2;
3935 for (
int i = 0; i < nocc_ij; i++ ){
3936 for (
int j = i; j < nocc_ij; j++ ){
3937 const double ij_factor = (( i == j ) ? SQRT_0p5 : 1.0 );
3938 const int counter_singlet = shift_singlet + i + ( j * ( j + 1 ) ) / 2 + linsize_ij_singlet * combined_ab_singlet;
3939 const int counter_triplet = shift_triplet + i + ( j * ( j - 1 ) ) / 2 + linsize_ij_triplet * combined_ab_triplet;
3941 const double ai_bj = integrals->
get_exchange( irrep_ij, irrep_ij, irrep_ab, irrep_ab, i, j, noa_ab + a, noa_ab + b );
3942 const double aj_bi = integrals->
get_exchange( irrep_ij, irrep_ij, irrep_ab, irrep_ab, j, i, noa_ab + a, noa_ab + b );
3943 target_singlet[ counter_singlet ] = 2 * ( ai_bj + aj_bi ) * ij_factor * ab_factor;
3944 if ((a<b) && (i<j)){ target_triplet[ counter_triplet ] = 6 * ( ai_bj - aj_bi ); }
3948 shift_singlet += linsize_ij_singlet * linsize_ab_singlet;
3949 shift_triplet += linsize_ij_triplet * linsize_ab_triplet;
3953 for (
int irrep_i = 0; irrep_i < num_irreps; irrep_i++ ){
3955 if ( irrep_i < irrep_j ){
3956 const int nocc_i = indices->
getNOCC( irrep_i );
3957 const int nocc_j = indices->
getNOCC( irrep_j );
3958 const int linsize_ij = nocc_i * nocc_j;
3959 for (
int irrep_a = 0; irrep_a < num_irreps; irrep_a++ ){
3961 if ( irrep_a < irrep_b ){
3962 assert( shift_singlet == shift_H_nonactive( indices, irrep_i, irrep_j, irrep_a, irrep_b, +1 ) );
3963 assert( shift_triplet == shift_H_nonactive( indices, irrep_i, irrep_j, irrep_a, irrep_b, -1 ) );
3964 const int nvir_a = indices->
getNVIRT( irrep_a );
3965 const int nvir_b = indices->
getNVIRT( irrep_b );
3966 const int noa_a = indices->
getNOCC( irrep_a ) + indices->
getNDMRG( irrep_a );
3967 const int noa_b = indices->
getNOCC( irrep_b ) + indices->
getNDMRG( irrep_b );
3968 const int linsize_ab = nvir_a * nvir_b;
3969 #pragma omp for schedule(static) 3970 for (
int combined_ab = 0; combined_ab < linsize_ab; combined_ab++ ){
3972 const int a = combined_ab % nvir_a;
3973 const int b = combined_ab / nvir_a;
3975 for (
int j = 0; j < nocc_j; j++ ){
3976 for (
int i = 0; i < nocc_i; i++ ){
3977 const int count_singlet = shift_singlet + i + nocc_i * ( j + nocc_j * combined_ab );
3978 const int count_triplet = shift_triplet + i + nocc_i * ( j + nocc_j * combined_ab );
3979 const double ai_bj = integrals->
get_exchange( irrep_i, irrep_j, irrep_a, irrep_b, i, j, noa_a + a, noa_b + b );
3980 const double aj_bi = integrals->
get_exchange( irrep_j, irrep_i, irrep_a, irrep_b, j, i, noa_a + a, noa_b + b );
3981 target_singlet[ count_singlet ] = 2 * ( ai_bj + aj_bi );
3982 target_triplet[ count_triplet ] = 6 * ( ai_bj - aj_bi );
3986 shift_singlet += linsize_ij * linsize_ab;
3987 shift_triplet += linsize_ij * linsize_ab;
3993 assert( shift_singlet == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_SINGLET ] );
3994 assert( shift_triplet == jump[ 1 + irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] - jump[ irrep + num_irreps * CHEMPS2_CASPT2_H_TRIPLET ] );
4000 int CheMPS2::CASPT2::jump_AC_active(
const DMRGSCFindices * idx,
const int irrep_t,
const int irrep_u,
const int irrep_v ){
4006 for (
int It = 0; It < n_irreps; It++ ){
4007 for (
int Iu = 0; Iu < n_irreps; Iu++ ){
4009 if (( It == irrep_t ) && ( Iu == irrep_u ) && ( Iv == irrep_v )){
4021 int CheMPS2::CASPT2::jump_BF_active(
const DMRGSCFindices * idx,
const int irrep_t,
const int irrep_u,
const int ST ){
4023 assert( irrep_t <= irrep_u );
4028 if ( irrep_tu == 0 ){
4029 for (
int Itu = 0; Itu < n_irreps; Itu++ ){
4030 if (( Itu == irrep_t ) && ( Itu == irrep_u )){
4037 for (
int It = 0; It < n_irreps; It++ ){
4040 if (( It == irrep_t ) && ( Iu == irrep_u )){
4052 void CheMPS2::CASPT2::make_FDE_FDG(){
4106 FDE_singlet =
new double***[ num_irreps ];
4107 FDE_triplet =
new double***[ num_irreps ];
4108 FDG_singlet =
new double***[ num_irreps ];
4109 FDG_triplet =
new double***[ num_irreps ];
4113 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
4115 const int SIZE_left = size_D[ irrep_left ];
4116 const int D2JUMP = SIZE_left / 2;
4117 FDE_singlet[ irrep_left ] =
new double**[ num_irreps ];
4118 FDE_triplet[ irrep_left ] =
new double**[ num_irreps ];
4119 FDG_singlet[ irrep_left ] =
new double**[ num_irreps ];
4120 FDG_triplet[ irrep_left ] =
new double**[ num_irreps ];
4122 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
4124 const int SIZE_right = indices->
getNDMRG( irrep_t );
4126 const int num_t = indices->
getNDMRG( irrep_t );
4129 const int num_w = indices->
getNDMRG( irrep_w );
4130 FDE_singlet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4131 FDE_triplet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4132 FDG_singlet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4133 FDG_triplet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4135 for (
int w = 0; w < num_w; w++ ){
4137 FDE_singlet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4138 FDE_triplet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4139 FDG_singlet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4140 FDG_triplet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4142 double * FDE_sing = FDE_singlet[ irrep_left ][ irrep_t ][ w ];
4143 double * FDE_trip = FDE_triplet[ irrep_left ][ irrep_t ][ w ];
4144 double * FDG_sing = FDG_singlet[ irrep_left ][ irrep_t ][ w ];
4145 double * FDG_trip = FDG_triplet[ irrep_left ][ irrep_t ][ w ];
4148 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
4150 const int num_x = indices->
getNDMRG( irrep_x );
4153 const int num_y = indices->
getNDMRG( irrep_y );
4155 for (
int t = 0; t < num_t; t++ ){
4156 for (
int y = 0; y < num_y; y++ ){
4157 for (
int x = 0; x < num_x; x++ ){
4158 const double gamma_ytxw = two_rdm[ d_y + y + LAS * ( d_t + t + LAS * ( d_x + x + LAS * ( d_w + w ))) ];
4159 const double gamma_ytwx = two_rdm[ d_y + y + LAS * ( d_t + t + LAS * ( d_w + w + LAS * ( d_x + x ))) ];
4160 FDE_sing[ jump_row + x + num_x * y + SIZE_left * t ] = - gamma_ytxw;
4161 FDE_sing[ D2JUMP + jump_row + x + num_x * y + SIZE_left * t ] = + gamma_ytxw + gamma_ytwx;
4162 FDE_trip[ jump_row + x + num_x * y + SIZE_left * t ] = - gamma_ytxw;
4163 FDE_trip[ D2JUMP + jump_row + x + num_x * y + SIZE_left * t ] = ( gamma_ytxw - gamma_ytwx ) / 3.0;
4164 const double gamma_ywtx = two_rdm[ d_y + y + LAS * ( d_w + w + LAS * ( d_t + t + LAS * ( d_x + x ))) ];
4165 const double gamma_ywxt = two_rdm[ d_y + y + LAS * ( d_w + w + LAS * ( d_x + x + LAS * ( d_t + t ))) ];
4166 FDG_sing[ jump_row + x + num_x * y + SIZE_left * t ] = + gamma_ywxt;
4167 FDG_sing[ D2JUMP + jump_row + x + num_x * y + SIZE_left * t ] = - gamma_ywxt - gamma_ywtx;
4168 FDG_trip[ jump_row + x + num_x * y + SIZE_left * t ] = - gamma_ywxt;
4169 FDG_trip[ D2JUMP + jump_row + x + num_x * y + SIZE_left * t ] = ( gamma_ywxt - gamma_ywtx ) / 3.0;
4174 if (( irrep_t == irrep_w ) && ( irrep_x == irrep_y )){
4175 for (
int x = 0; x < num_x; x++ ){
4176 for (
int y = 0; y < num_x; y++ ){
4177 const double gamma_yx = one_rdm[ d_y + y + LAS * ( d_x + x ) ];
4178 FDE_sing[ jump_row + x + num_x * y + SIZE_left * w ] += 2 * gamma_yx;
4179 FDE_sing[ D2JUMP + jump_row + x + num_x * y + SIZE_left * w ] -= gamma_yx;
4180 FDE_trip[ jump_row + x + num_x * y + SIZE_left * w ] += 2 * gamma_yx;
4181 FDE_trip[ D2JUMP + jump_row + x + num_x * y + SIZE_left * w ] -= gamma_yx;
4186 if (( irrep_t == irrep_x ) && ( irrep_y == irrep_w )){
4187 for (
int y = 0; y < num_y; y++ ){
4188 const double gamma_yw = one_rdm[ d_y + y + LAS * ( d_w + w ) ];
4189 for (
int tx = 0; tx < num_x; tx++ ){
4190 FDE_sing[ jump_row + tx + num_x * y + SIZE_left * tx ] -= gamma_yw;
4191 FDE_sing[ D2JUMP + jump_row + tx + num_x * y + SIZE_left * tx ] -= gamma_yw;
4192 FDE_trip[ jump_row + tx + num_x * y + SIZE_left * tx ] -= gamma_yw;
4193 FDE_trip[ D2JUMP + jump_row + tx + num_x * y + SIZE_left * tx ] += gamma_yw;
4198 if (( irrep_t == irrep_y ) && ( irrep_w == irrep_x )){
4199 for (
int t = 0; t < num_t; t++ ){
4200 for (
int y = 0; y < num_y; y++ ){
4201 const double gamma_yt = one_rdm[ d_y + y + LAS * ( d_t + t ) ];
4202 FDG_sing[ jump_row + w + num_x * y + SIZE_left * t ] += gamma_yt;
4203 FDG_sing[ D2JUMP + jump_row + w + num_x * y + SIZE_left * t ] += gamma_yt;
4204 FDG_trip[ jump_row + w + num_x * y + SIZE_left * t ] -= gamma_yt;
4205 FDG_trip[ D2JUMP + jump_row + w + num_x * y + SIZE_left * t ] += gamma_yt;
4209 jump_row += num_x * num_y;
4217 void CheMPS2::CASPT2::make_FEH_FGH(){
4229 FEH =
new double**[ num_irreps ];
4230 FGH =
new double**[ num_irreps ];
4232 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
4234 const int NACT = indices->
getNDMRG( irrep_x );
4235 FEH[ irrep_x ] =
new double*[ NACT ];
4236 FGH[ irrep_x ] =
new double*[ NACT ];
4238 for (
int w = 0; w < NACT; w++ ){
4239 FEH[ irrep_x ][ w ] =
new double[ NACT ];
4240 FGH[ irrep_x ][ w ] =
new double[ NACT ];
4242 for (
int x = 0; x < NACT; x++ ){
4243 FEH[ irrep_x ][ w ][ x ] = + SEE[ irrep_x ][ x + NACT * w ];
4244 FGH[ irrep_x ][ w ][ x ] = - SGG[ irrep_x ][ x + NACT * w ];
4251 void CheMPS2::CASPT2::make_FBE_FFG_singlet(){
4261 FBE_singlet =
new double***[ num_irreps ];
4262 FFG_singlet =
new double***[ num_irreps ];
4264 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
4266 assert( size_B_singlet[ irrep_left ] == size_F_singlet[ irrep_left ] );
4267 const int SIZE_left = size_B_singlet[ irrep_left ];
4268 FBE_singlet[ irrep_left ] =
new double**[ num_irreps ];
4269 FFG_singlet[ irrep_left ] =
new double**[ num_irreps ];
4271 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
4273 const int SIZE_right = indices->
getNDMRG( irrep_t );
4274 const int num_t = indices->
getNDMRG( irrep_t );
4276 const int num_w = indices->
getNDMRG( irrep_w );
4277 FBE_singlet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4278 FFG_singlet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4280 for (
int w = 0; w < num_w; w++ ){
4282 FBE_singlet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4283 FFG_singlet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4285 double * BEptr = FBE_singlet[ irrep_left ][ irrep_t ][ w ];
4286 double * FGptr = FFG_singlet[ irrep_left ][ irrep_t ][ w ];
4288 if ( irrep_left == 0 ){
4289 const int jump_singlet = jump_BF_active( indices, irrep_t, irrep_w, +1 );
4290 for (
int t = 0; t < w; t++ ){
4291 for (
int xy = 0; xy < SIZE_left; xy++ ){
4292 BEptr[ xy + SIZE_left * t ] = + SBB_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + t + ( w * ( w + 1 ) ) / 2 ) ];
4293 FGptr[ xy + SIZE_left * t ] = - SFF_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + t + ( w * ( w + 1 ) ) / 2 ) ];
4296 for (
int t = w; t < num_t; t++ ){
4297 for (
int xy = 0; xy < SIZE_left; xy++ ){
4298 BEptr[ xy + SIZE_left * t ] = + SBB_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + w + ( t * ( t + 1 ) ) / 2 ) ];
4299 FGptr[ xy + SIZE_left * t ] = - SFF_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + w + ( t * ( t + 1 ) ) / 2 ) ];
4303 if ( irrep_t < irrep_w ){
4304 const int jump_singlet = jump_BF_active( indices, irrep_t, irrep_w, +1 );
4305 for (
int t = 0; t < num_t; t++ ){
4306 for (
int xy = 0; xy < SIZE_left; xy++ ){
4307 BEptr[ xy + SIZE_left * t ] = + SBB_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + t + num_t * w ) ];
4308 FGptr[ xy + SIZE_left * t ] = - SFF_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + t + num_t * w ) ];
4312 const int jump_singlet = jump_BF_active( indices, irrep_w, irrep_t, +1 );
4313 for (
int t = 0; t < num_t; t++ ){
4314 for (
int xy = 0; xy < SIZE_left; xy++ ){
4315 BEptr[ xy + SIZE_left * t ] = + SBB_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + w + num_w * t ) ];
4316 FGptr[ xy + SIZE_left * t ] = - SFF_singlet[ irrep_left ][ xy + SIZE_left * ( jump_singlet + w + num_w * t ) ];
4327 void CheMPS2::CASPT2::make_FBE_FFG_triplet(){
4337 FBE_triplet =
new double***[ num_irreps ];
4338 FFG_triplet =
new double***[ num_irreps ];
4340 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
4342 assert( size_B_triplet[ irrep_left ] == size_F_triplet[ irrep_left ] );
4343 const int SIZE_left = size_B_triplet[ irrep_left ];
4344 FBE_triplet[ irrep_left ] =
new double**[ num_irreps ];
4345 FFG_triplet[ irrep_left ] =
new double**[ num_irreps ];
4347 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
4349 const int SIZE_right = indices->
getNDMRG( irrep_t );
4350 const int num_t = indices->
getNDMRG( irrep_t );
4352 const int num_w = indices->
getNDMRG( irrep_w );
4353 FBE_triplet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4354 FFG_triplet[ irrep_left ][ irrep_t ] =
new double*[ num_w ];
4356 for (
int w = 0; w < num_w; w++ ){
4358 FBE_triplet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4359 FFG_triplet[ irrep_left ][ irrep_t ][ w ] =
new double[ SIZE_left * SIZE_right ];
4361 double * BEptr = FBE_triplet[ irrep_left ][ irrep_t ][ w ];
4362 double * FGptr = FFG_triplet[ irrep_left ][ irrep_t ][ w ];
4364 if ( irrep_left == 0 ){
4365 const int jump_triplet = jump_BF_active( indices, irrep_t, irrep_w, -1 );
4366 for (
int t = 0; t < w; t++ ){
4367 for (
int xy = 0; xy < SIZE_left; xy++ ){
4368 BEptr[ xy + SIZE_left * t ] = + SBB_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + t + ( w * ( w - 1 ) ) / 2 ) ];
4369 FGptr[ xy + SIZE_left * t ] = + SFF_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + t + ( w * ( w - 1 ) ) / 2 ) ];
4372 for (
int xy = 0; xy < SIZE_left; xy++ ){
4373 BEptr[ xy + SIZE_left * w ] = 0.0;
4374 FGptr[ xy + SIZE_left * w ] = 0.0;
4376 for (
int t = w+1; t < num_t; t++ ){
4377 for (
int xy = 0; xy < SIZE_left; xy++ ){
4378 BEptr[ xy + SIZE_left * t ] = - SBB_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + w + ( t * ( t - 1 ) ) / 2 ) ];
4379 FGptr[ xy + SIZE_left * t ] = - SFF_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + w + ( t * ( t - 1 ) ) / 2 ) ];
4383 if ( irrep_t < irrep_w ){
4384 const int jump_triplet = jump_BF_active( indices, irrep_t, irrep_w, -1 );
4385 for (
int t = 0; t < num_t; t++ ){
4386 for (
int xy = 0; xy < SIZE_left; xy++ ){
4387 BEptr[ xy + SIZE_left * t ] = + SBB_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + t + num_t * w ) ];
4388 FGptr[ xy + SIZE_left * t ] = + SFF_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + t + num_t * w ) ];
4392 const int jump_triplet = jump_BF_active( indices, irrep_w, irrep_t, -1 );
4393 for (
int t = 0; t < num_t; t++ ){
4394 for (
int xy = 0; xy < SIZE_left; xy++ ){
4395 BEptr[ xy + SIZE_left * t ] = - SBB_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + w + num_w * t ) ];
4396 FGptr[ xy + SIZE_left * t ] = - SFF_triplet[ irrep_left ][ xy + SIZE_left * ( jump_triplet + w + num_w * t ) ];
4407 void CheMPS2::CASPT2::make_FAB_FCF_singlet(){
4433 FAB_singlet =
new double***[ num_irreps ];
4434 FCF_singlet =
new double***[ num_irreps ];
4438 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
4440 assert( size_A[ irrep_left ] == size_C[ irrep_left ] );
4441 const int SIZE_left = size_A[ irrep_left ];
4442 FAB_singlet[ irrep_left ] =
new double**[ num_irreps ];
4443 FCF_singlet[ irrep_left ] =
new double**[ num_irreps ];
4447 assert( size_B_singlet[ 0 ] == size_F_singlet[ 0 ] );
4448 const int SIZE_right = size_B_singlet[ 0 ];
4449 const int num_w = indices->
getNDMRG( irrep_left );
4450 FAB_singlet[ irrep_left ][ 0 ] =
new double*[ num_w ];
4451 FCF_singlet[ irrep_left ][ 0 ] =
new double*[ num_w ];
4453 for (
int w = 0; w < num_w; w++ ){
4455 FAB_singlet[ irrep_left ][ 0 ][ w ] =
new double[ SIZE_left * SIZE_right ];
4456 FCF_singlet[ irrep_left ][ 0 ][ w ] =
new double[ SIZE_left * SIZE_right ];
4458 double * ABptr = FAB_singlet[ irrep_left ][ 0 ][ w ];
4459 double * CFptr = FCF_singlet[ irrep_left ][ 0 ][ w ];
4462 for (
int irrep_ut = 0; irrep_ut < num_irreps; irrep_ut++ ){
4464 const int num_ut = indices->
getNDMRG( irrep_ut );
4465 assert( jump_col == jump_BF_active( indices, irrep_ut, irrep_ut, +1 ) );
4467 const int jump_AB = jump_AC_active( indices, irrep_ut, irrep_ut, irrep_left );
4468 const int jump_CF = jump_AC_active( indices, irrep_ut, irrep_left, irrep_ut );
4470 for (
int t = 0; t < num_ut; t++ ){
4471 for (
int u = t; u < num_ut; u++ ){
4472 for (
int xyz = 0; xyz < SIZE_left; xyz++ ){
4473 ABptr[ xyz + SIZE_left * ( jump_col + t + ( u * ( u + 1 ) ) / 2 ) ] = ( - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB + u + num_ut * ( t + num_ut * w )) ]
4474 - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB + t + num_ut * ( u + num_ut * w )) ] );
4475 CFptr[ xyz + SIZE_left * ( jump_col + t + ( u * ( u + 1 ) ) / 2 ) ] = ( + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF + u + num_ut * ( w + num_w * t )) ]
4476 + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF + t + num_ut * ( w + num_w * u )) ] );
4482 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
4484 const int num_x = indices->
getNDMRG( irrep_x );
4485 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
4488 const int num_y = indices->
getNDMRG( irrep_y );
4490 const int num_z = indices->
getNDMRG( irrep_z );
4491 assert( jump_row == jump_AC_active( indices, irrep_x, irrep_y, irrep_z ) );
4493 if ( irrep_ut == irrep_left ){
4496 for (
int u = w; u < num_ut; u++ ){
4497 for (
int x = 0; x < num_x; x++ ){
4498 for (
int y = 0; y < num_y; y++ ){
4499 for (
int z = 0; z < num_z; z++ ){
4500 const double gamma_zuyx = two_rdm[ d_z + z + LAS * ( d_ut + u + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4501 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( u * ( u + 1 ) ) / 2 ) ] -= gamma_zuyx;
4508 if ( irrep_ut == irrep_y ){
4509 for (
int x = 0; x < num_x; x++ ){
4510 for (
int z = 0; z < num_z; z++ ){
4511 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4512 for (
int uy = w; uy < num_y; uy++ ){
4513 ABptr[ jump_row + x + num_x * ( uy + num_y * z ) + SIZE_left * ( jump_col + w + ( uy * ( uy + 1 ) ) / 2 ) ] -= gamma_zx;
4520 if ( irrep_ut == irrep_x ){
4521 for (
int y = 0; y < num_y; y++ ){
4522 for (
int z = 0; z < num_z; z++ ){
4523 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4524 for (
int ux = w; ux < num_x; ux++ ){
4525 ABptr[ jump_row + ux + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( ux * ( ux + 1 ) ) / 2 ) ] += 2 * gamma_zy;
4532 for (
int u = w; u < num_ut; u++ ){
4533 for (
int x = 0; x < num_x; x++ ){
4534 for (
int y = 0; y < num_y; y++ ){
4535 for (
int z = 0; z < num_z; z++ ){
4536 const double gamma_zxyu = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_ut + u ))) ];
4537 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( u * ( u + 1 ) ) / 2 ) ] -= gamma_zxyu;
4544 if ( irrep_x == irrep_y ){
4545 for (
int u = w; u < num_ut; u++ ){
4546 for (
int z = 0; z < num_z; z++ ){
4547 const double gamma_zu = one_rdm[ d_z + z + LAS * ( d_ut + u ) ];
4548 for (
int xy = 0; xy < num_x; xy++ ){
4549 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + w + ( u * ( u + 1 ) ) / 2 ) ] -= gamma_zu;
4556 for (
int t = 0; t <= w; t++ ){
4557 for (
int x = 0; x < num_x; x++ ){
4558 for (
int y = 0; y < num_y; y++ ){
4559 for (
int z = 0; z < num_z; z++ ){
4560 const double gamma_ztyx = two_rdm[ d_z + z + LAS * ( d_ut + t + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4561 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + ( w * ( w + 1 ) ) / 2 ) ] -= gamma_ztyx;
4568 if ( irrep_ut == irrep_y ){
4569 for (
int x = 0; x < num_x; x++ ){
4570 for (
int z = 0; z < num_z; z++ ){
4571 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4572 for (
int ty = 0; ty <= w; ty++ ){
4573 ABptr[ jump_row + x + num_x * ( ty + num_y * z ) + SIZE_left * ( jump_col + ty + ( w * ( w + 1 ) ) / 2 ) ] -= gamma_zx;
4580 if ( irrep_ut == irrep_x ){
4581 for (
int y = 0; y < num_y; y++ ){
4582 for (
int z = 0; z < num_z; z++ ){
4583 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4584 for (
int tx = 0; tx <= w; tx++ ){
4585 ABptr[ jump_row + tx + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + tx + ( w * ( w + 1 ) ) / 2 ) ] += 2 * gamma_zy;
4592 for (
int t = 0; t <= w; t++ ){
4593 for (
int x = 0; x < num_x; x++ ){
4594 for (
int y = 0; y < num_y; y++ ){
4595 for (
int z = 0; z < num_z; z++ ){
4596 const double gamma_zxyt = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_ut + t ))) ];
4597 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + ( w * ( w + 1 ) ) / 2 ) ] -= gamma_zxyt;
4604 if ( irrep_x == irrep_y ){
4605 for (
int t = 0; t <= w; t++ ){
4606 for (
int z = 0; z < num_z; z++ ){
4607 const double gamma_zt = one_rdm[ d_z + z + LAS * ( d_ut + t ) ];
4608 for (
int xy = 0; xy < num_x; xy++ ){
4609 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + t + ( w * ( w + 1 ) ) / 2 ) ] -= gamma_zt;
4615 jump_row += num_x * num_y * num_z;
4618 jump_col += ( num_ut * ( num_ut + 1 ) ) / 2;
4623 for (
int irrep_right = 1; irrep_right < num_irreps; irrep_right++ ){
4625 assert( size_B_singlet[ irrep_right ] == size_F_singlet[ irrep_right ] );
4626 const int SIZE_right = size_B_singlet[ irrep_right ];
4628 const int num_w = indices->
getNDMRG( irrep_w );
4629 FAB_singlet[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
4630 FCF_singlet[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
4632 for (
int w = 0; w < num_w; w++ ){
4634 FAB_singlet[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
4635 FCF_singlet[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
4637 double * ABptr = FAB_singlet[ irrep_left ][ irrep_right ][ w ];
4638 double * CFptr = FCF_singlet[ irrep_left ][ irrep_right ][ w ];
4641 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
4643 if ( irrep_t < irrep_u ){
4645 const int num_t = indices->
getNDMRG( irrep_t );
4647 const int num_u = indices->
getNDMRG( irrep_u );
4648 assert( jump_col == jump_BF_active( indices, irrep_t, irrep_u, +1 ) );
4650 const int jump_AB1 = jump_AC_active( indices, irrep_u, irrep_t, irrep_w );
4651 const int jump_AB2 = jump_AC_active( indices, irrep_t, irrep_u, irrep_w );
4652 const int jump_CF1 = jump_AC_active( indices, irrep_u, irrep_w, irrep_t );
4653 const int jump_CF2 = jump_AC_active( indices, irrep_t, irrep_w, irrep_u );
4655 for (
int t = 0; t < num_t; t++ ){
4656 for (
int u = 0; u < num_u; u++ ){
4657 for (
int xyz = 0; xyz < SIZE_left; xyz++ ){
4658 ABptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = ( - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB1 + u + num_u * ( t + num_t * w )) ]
4659 - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB2 + t + num_t * ( u + num_u * w )) ] );
4660 CFptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = ( + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF1 + u + num_u * ( w + num_w * t )) ]
4661 + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF2 + t + num_t * ( w + num_w * u )) ] );
4667 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
4669 const int num_x = indices->
getNDMRG( irrep_x );
4670 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
4673 const int num_y = indices->
getNDMRG( irrep_y );
4675 const int num_z = indices->
getNDMRG( irrep_z );
4676 assert( jump_row == jump_AC_active( indices, irrep_x, irrep_y, irrep_z ) );
4678 if ( irrep_t == irrep_w ){
4681 for (
int u = 0; u < num_u; u++ ){
4682 for (
int x = 0; x < num_x; x++ ){
4683 for (
int y = 0; y < num_y; y++ ){
4684 for (
int z = 0; z < num_z; z++ ){
4685 const double gamma_zuyx = two_rdm[ d_z + z + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4686 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] -= gamma_zuyx;
4693 if ( irrep_u == irrep_y ){
4694 for (
int x = 0; x < num_x; x++ ){
4695 for (
int z = 0; z < num_z; z++ ){
4696 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4697 for (
int uy = 0; uy < num_y; uy++ ){
4698 ABptr[ jump_row + x + num_x * ( uy + num_y * z ) + SIZE_left * ( jump_col + w + num_w * uy ) ] -= gamma_zx;
4705 if ( irrep_u == irrep_x ){
4706 for (
int y = 0; y < num_y; y++ ){
4707 for (
int z = 0; z < num_z; z++ ){
4708 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4709 for (
int ux = 0; ux < num_x; ux++ ){
4710 ABptr[ jump_row + ux + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * ux ) ] += 2 * gamma_zy;
4717 for (
int u = 0; u < num_u; u++ ){
4718 for (
int x = 0; x < num_x; x++ ){
4719 for (
int y = 0; y < num_y; y++ ){
4720 for (
int z = 0; z < num_z; z++ ){
4721 const double gamma_zxyu = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_u + u ))) ];
4722 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] -= gamma_zxyu;
4729 if ( irrep_x == irrep_y ){
4730 for (
int u = 0; u < num_u; u++ ){
4731 for (
int z = 0; z < num_z; z++ ){
4732 const double gamma_zu = one_rdm[ d_z + z + LAS * ( d_u + u ) ];
4733 for (
int xy = 0; xy < num_x; xy++ ){
4734 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] -= gamma_zu;
4741 if ( irrep_u == irrep_w ){
4744 for (
int t = 0; t < num_t; t++ ){
4745 for (
int x = 0; x < num_x; x++ ){
4746 for (
int y = 0; y < num_y; y++ ){
4747 for (
int z = 0; z < num_z; z++ ){
4748 const double gamma_ztyx = two_rdm[ d_z + z + LAS * ( d_t + t + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4749 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] -= gamma_ztyx;
4756 if ( irrep_t == irrep_y ){
4757 for (
int x = 0; x < num_x; x++ ){
4758 for (
int z = 0; z < num_z; z++ ){
4759 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4760 for (
int ty = 0; ty < num_y; ty++ ){
4761 ABptr[ jump_row + x + num_x * ( ty + num_y * z ) + SIZE_left * ( jump_col + ty + num_y * w ) ] -= gamma_zx;
4768 if ( irrep_t == irrep_x ){
4769 for (
int y = 0; y < num_y; y++ ){
4770 for (
int z = 0; z < num_z; z++ ){
4771 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4772 for (
int tx = 0; tx < num_x; tx++ ){
4773 ABptr[ jump_row + tx + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + tx + num_x * w ) ] += 2 * gamma_zy;
4780 for (
int t = 0; t < num_t; t++ ){
4781 for (
int x = 0; x < num_x; x++ ){
4782 for (
int y = 0; y < num_y; y++ ){
4783 for (
int z = 0; z < num_z; z++ ){
4784 const double gamma_zxyt = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_t + t ))) ];
4785 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] -= gamma_zxyt;
4792 if ( irrep_x == irrep_y ){
4793 for (
int t = 0; t < num_t; t++ ){
4794 for (
int z = 0; z < num_z; z++ ){
4795 const double gamma_zt = one_rdm[ d_z + z + LAS * ( d_t + t ) ];
4796 for (
int xy = 0; xy < num_x; xy++ ){
4797 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] -= gamma_zt;
4803 jump_row += num_x * num_y * num_z;
4806 jump_col += num_t * num_u;
4815 void CheMPS2::CASPT2::make_FAB_FCF_triplet(){
4841 FAB_triplet =
new double***[ num_irreps ];
4842 FCF_triplet =
new double***[ num_irreps ];
4846 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
4848 assert( size_A[ irrep_left ] == size_C[ irrep_left ] );
4849 const int SIZE_left = size_A[ irrep_left ];
4850 FAB_triplet[ irrep_left ] =
new double**[ num_irreps ];
4851 FCF_triplet[ irrep_left ] =
new double**[ num_irreps ];
4855 assert( size_B_triplet[ 0 ] == size_F_triplet[ 0 ] );
4856 const int SIZE_right = size_B_triplet[ 0 ];
4857 const int num_w = indices->
getNDMRG( irrep_left );
4858 FAB_triplet[ irrep_left ][ 0 ] =
new double*[ num_w ];
4859 FCF_triplet[ irrep_left ][ 0 ] =
new double*[ num_w ];
4861 for (
int w = 0; w < num_w; w++ ){
4863 FAB_triplet[ irrep_left ][ 0 ][ w ] =
new double[ SIZE_left * SIZE_right ];
4864 FCF_triplet[ irrep_left ][ 0 ][ w ] =
new double[ SIZE_left * SIZE_right ];
4866 double * ABptr = FAB_triplet[ irrep_left ][ 0 ][ w ];
4867 double * CFptr = FCF_triplet[ irrep_left ][ 0 ][ w ];
4870 for (
int irrep_ut = 0; irrep_ut < num_irreps; irrep_ut++ ){
4872 const int num_ut = indices->
getNDMRG( irrep_ut );
4873 assert( jump_col == jump_BF_active( indices, irrep_ut, irrep_ut, -1 ) );
4875 const int jump_AB = jump_AC_active( indices, irrep_ut, irrep_ut, irrep_left );
4876 const int jump_CF = jump_AC_active( indices, irrep_ut, irrep_left, irrep_ut );
4878 for (
int t = 0; t < num_ut; t++ ){
4879 for (
int u = t+1; u < num_ut; u++ ){
4880 for (
int xyz = 0; xyz < SIZE_left; xyz++ ){
4881 ABptr[ xyz + SIZE_left * ( jump_col + t + ( u * ( u - 1 ) ) / 2 ) ] = ( - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB + u + num_ut * ( t + num_ut * w )) ]
4882 + SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB + t + num_ut * ( u + num_ut * w )) ] );
4883 CFptr[ xyz + SIZE_left * ( jump_col + t + ( u * ( u - 1 ) ) / 2 ) ] = ( + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF + u + num_ut * ( w + num_w * t )) ]
4884 - SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF + t + num_ut * ( w + num_w * u )) ] );
4890 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
4892 const int num_x = indices->
getNDMRG( irrep_x );
4893 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
4896 const int num_y = indices->
getNDMRG( irrep_y );
4898 const int num_z = indices->
getNDMRG( irrep_z );
4899 assert( jump_row == jump_AC_active( indices, irrep_x, irrep_y, irrep_z ) );
4901 if ( irrep_ut == irrep_left ){
4904 for (
int u = w+1; u < num_ut; u++ ){
4905 for (
int x = 0; x < num_x; x++ ){
4906 for (
int y = 0; y < num_y; y++ ){
4907 for (
int z = 0; z < num_z; z++ ){
4908 const double gamma_zuyx = two_rdm[ d_z + z + LAS * ( d_ut + u + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4909 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( u * ( u - 1 ) ) / 2 ) ] -= 3 * gamma_zuyx;
4916 if ( irrep_ut == irrep_y ){
4917 for (
int x = 0; x < num_x; x++ ){
4918 for (
int z = 0; z < num_z; z++ ){
4919 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4920 for (
int uy = w+1; uy < num_y; uy++ ){
4921 ABptr[ jump_row + x + num_x * ( uy + num_y * z ) + SIZE_left * ( jump_col + w + ( uy * ( uy - 1 ) ) / 2 ) ] -= 3 * gamma_zx;
4928 if ( irrep_ut == irrep_x ){
4929 for (
int y = 0; y < num_y; y++ ){
4930 for (
int z = 0; z < num_z; z++ ){
4931 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4932 for (
int ux = w+1; ux < num_x; ux++ ){
4933 ABptr[ jump_row + ux + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( ux * ( ux - 1 ) ) / 2 ) ] += 6 * gamma_zy;
4940 for (
int u = w+1; u < num_ut; u++ ){
4941 for (
int x = 0; x < num_x; x++ ){
4942 for (
int y = 0; y < num_y; y++ ){
4943 for (
int z = 0; z < num_z; z++ ){
4944 const double gamma_zxyu = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_ut + u ))) ];
4945 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + ( u * ( u - 1 ) ) / 2 ) ] += gamma_zxyu;
4952 if ( irrep_x == irrep_y ){
4953 for (
int u = w+1; u < num_ut; u++ ){
4954 for (
int z = 0; z < num_z; z++ ){
4955 const double gamma_zu = one_rdm[ d_z + z + LAS * ( d_ut + u ) ];
4956 for (
int xy = 0; xy < num_x; xy++ ){
4957 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + w + ( u * ( u - 1 ) ) / 2 ) ] += gamma_zu;
4964 for (
int t = 0; t < w; t++ ){
4965 for (
int x = 0; x < num_x; x++ ){
4966 for (
int y = 0; y < num_y; y++ ){
4967 for (
int z = 0; z < num_z; z++ ){
4968 const double gamma_ztyx = two_rdm[ d_z + z + LAS * ( d_ut + t + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
4969 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + ( w * ( w - 1 ) ) / 2 ) ] += 3 * gamma_ztyx;
4976 if ( irrep_ut == irrep_y ){
4977 for (
int x = 0; x < num_x; x++ ){
4978 for (
int z = 0; z < num_z; z++ ){
4979 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
4980 for (
int ty = 0; ty < w; ty++ ){
4981 ABptr[ jump_row + x + num_x * ( ty + num_y * z ) + SIZE_left * ( jump_col + ty + ( w * ( w - 1 ) ) / 2 ) ] += 3 * gamma_zx;
4988 if ( irrep_ut == irrep_x ){
4989 for (
int y = 0; y < num_y; y++ ){
4990 for (
int z = 0; z < num_z; z++ ){
4991 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
4992 for (
int tx = 0; tx < w; tx++ ){
4993 ABptr[ jump_row + tx + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + tx + ( w * ( w - 1 ) ) / 2 ) ] -= 6 * gamma_zy;
5000 for (
int t = 0; t < w; t++ ){
5001 for (
int x = 0; x < num_x; x++ ){
5002 for (
int y = 0; y < num_y; y++ ){
5003 for (
int z = 0; z < num_z; z++ ){
5004 const double gamma_zxyt = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_ut + t ))) ];
5005 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + ( w * ( w - 1 ) ) / 2 ) ] -= gamma_zxyt;
5012 if ( irrep_x == irrep_y ){
5013 for (
int t = 0; t < w; t++ ){
5014 for (
int z = 0; z < num_z; z++ ){
5015 const double gamma_zt = one_rdm[ d_z + z + LAS * ( d_ut + t ) ];
5016 for (
int xy = 0; xy < num_x; xy++ ){
5017 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + t + ( w * ( w - 1 ) ) / 2 ) ] -= gamma_zt;
5023 jump_row += num_x * num_y * num_z;
5026 jump_col += ( num_ut * ( num_ut - 1 ) ) / 2;
5031 for (
int irrep_right = 1; irrep_right < num_irreps; irrep_right++ ){
5033 assert( size_B_triplet[ irrep_right ] == size_F_triplet[ irrep_right ] );
5034 const int SIZE_right = size_B_triplet[ irrep_right ];
5036 const int num_w = indices->
getNDMRG( irrep_w );
5037 FAB_triplet[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
5038 FCF_triplet[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
5040 for (
int w = 0; w < num_w; w++ ){
5042 FAB_triplet[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
5043 FCF_triplet[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
5045 double * ABptr = FAB_triplet[ irrep_left ][ irrep_right ][ w ];
5046 double * CFptr = FCF_triplet[ irrep_left ][ irrep_right ][ w ];
5049 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
5051 if ( irrep_t < irrep_u ){
5053 const int num_t = indices->
getNDMRG( irrep_t );
5055 const int num_u = indices->
getNDMRG( irrep_u );
5056 assert( jump_col == jump_BF_active( indices, irrep_t, irrep_u, -1 ) );
5058 const int jump_AB1 = jump_AC_active( indices, irrep_u, irrep_t, irrep_w );
5059 const int jump_AB2 = jump_AC_active( indices, irrep_t, irrep_u, irrep_w );
5060 const int jump_CF1 = jump_AC_active( indices, irrep_u, irrep_w, irrep_t );
5061 const int jump_CF2 = jump_AC_active( indices, irrep_t, irrep_w, irrep_u );
5063 for (
int t = 0; t < num_t; t++ ){
5064 for (
int u = 0; u < num_u; u++ ){
5065 for (
int xyz = 0; xyz < SIZE_left; xyz++ ){
5066 ABptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = ( - SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB1 + u + num_u * ( t + num_t * w )) ]
5067 + SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AB2 + t + num_t * ( u + num_u * w )) ] );
5068 CFptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = ( + SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF1 + u + num_u * ( w + num_w * t )) ]
5069 - SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CF2 + t + num_t * ( w + num_w * u )) ] );
5075 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
5077 const int num_x = indices->
getNDMRG( irrep_x );
5078 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
5081 const int num_y = indices->
getNDMRG( irrep_y );
5083 const int num_z = indices->
getNDMRG( irrep_z );
5084 assert( jump_row == jump_AC_active( indices, irrep_x, irrep_y, irrep_z ) );
5086 if ( irrep_t == irrep_w ){
5089 for (
int u = 0; u < num_u; u++ ){
5090 for (
int x = 0; x < num_x; x++ ){
5091 for (
int y = 0; y < num_y; y++ ){
5092 for (
int z = 0; z < num_z; z++ ){
5093 const double gamma_zuyx = two_rdm[ d_z + z + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
5094 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] -= 3 * gamma_zuyx;
5101 if ( irrep_u == irrep_y ){
5102 for (
int x = 0; x < num_x; x++ ){
5103 for (
int z = 0; z < num_z; z++ ){
5104 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
5105 for (
int uy = 0; uy < num_y; uy++ ){
5106 ABptr[ jump_row + x + num_x * ( uy + num_y * z ) + SIZE_left * ( jump_col + w + num_w * uy ) ] -= 3 * gamma_zx;
5113 if ( irrep_u == irrep_x ){
5114 for (
int y = 0; y < num_y; y++ ){
5115 for (
int z = 0; z < num_z; z++ ){
5116 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
5117 for (
int ux = 0; ux < num_x; ux++ ){
5118 ABptr[ jump_row + ux + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * ux ) ] += 6 * gamma_zy;
5125 for (
int u = 0; u < num_u; u++ ){
5126 for (
int x = 0; x < num_x; x++ ){
5127 for (
int y = 0; y < num_y; y++ ){
5128 for (
int z = 0; z < num_z; z++ ){
5129 const double gamma_zxyu = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_u + u ))) ];
5130 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] += gamma_zxyu;
5137 if ( irrep_x == irrep_y ){
5138 for (
int u = 0; u < num_u; u++ ){
5139 for (
int z = 0; z < num_z; z++ ){
5140 const double gamma_zu = one_rdm[ d_z + z + LAS * ( d_u + u ) ];
5141 for (
int xy = 0; xy < num_x; xy++ ){
5142 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + w + num_w * u ) ] += gamma_zu;
5149 if ( irrep_u == irrep_w ){
5152 for (
int t = 0; t < num_t; t++ ){
5153 for (
int x = 0; x < num_x; x++ ){
5154 for (
int y = 0; y < num_y; y++ ){
5155 for (
int z = 0; z < num_z; z++ ){
5156 const double gamma_ztyx = two_rdm[ d_z + z + LAS * ( d_t + t + LAS * ( d_y + y + LAS * ( d_x + x ))) ];
5157 ABptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] += 3 * gamma_ztyx;
5164 if ( irrep_t == irrep_y ){
5165 for (
int x = 0; x < num_x; x++ ){
5166 for (
int z = 0; z < num_z; z++ ){
5167 const double gamma_zx = one_rdm[ d_z + z + LAS * ( d_x + x ) ];
5168 for (
int ty = 0; ty < num_y; ty++ ){
5169 ABptr[ jump_row + x + num_x * ( ty + num_y * z ) + SIZE_left * ( jump_col + ty + num_y * w ) ] += 3 * gamma_zx;
5176 if ( irrep_t == irrep_x ){
5177 for (
int y = 0; y < num_y; y++ ){
5178 for (
int z = 0; z < num_z; z++ ){
5179 const double gamma_zy = one_rdm[ d_z + z + LAS * ( d_y + y ) ];
5180 for (
int tx = 0; tx < num_x; tx++ ){
5181 ABptr[ jump_row + tx + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + tx + num_x * w ) ] -= 6 * gamma_zy;
5188 for (
int t = 0; t < num_t; t++ ){
5189 for (
int x = 0; x < num_x; x++ ){
5190 for (
int y = 0; y < num_y; y++ ){
5191 for (
int z = 0; z < num_z; z++ ){
5192 const double gamma_zxyt = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_t + t ))) ];
5193 CFptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] -= gamma_zxyt;
5200 if ( irrep_x == irrep_y ){
5201 for (
int t = 0; t < num_t; t++ ){
5202 for (
int z = 0; z < num_z; z++ ){
5203 const double gamma_zt = one_rdm[ d_z + z + LAS * ( d_t + t ) ];
5204 for (
int xy = 0; xy < num_x; xy++ ){
5205 CFptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + t + num_t * w ) ] -= gamma_zt;
5211 jump_row += num_x * num_y * num_z;
5214 jump_col += num_t * num_u;
5223 void CheMPS2::CASPT2::make_FAD_FCD(){
5245 FAD =
new double***[ num_irreps ];
5246 FCD =
new double***[ num_irreps ];
5250 for (
int irrep_left = 0; irrep_left < num_irreps; irrep_left++ ){
5252 assert( size_A[ irrep_left ] == size_C[ irrep_left ] );
5253 const int SIZE_left = size_A[ irrep_left ];
5254 FAD[ irrep_left ] =
new double**[ num_irreps ];
5255 FCD[ irrep_left ] =
new double**[ num_irreps ];
5257 for (
int irrep_right = 0; irrep_right < num_irreps; irrep_right++ ){
5259 const int SIZE_right = size_D[ irrep_right ];
5260 const int D2JUMP = SIZE_right / 2;
5263 const int num_w = indices->
getNDMRG( irrep_w );
5264 FAD[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
5265 FCD[ irrep_left ][ irrep_right ] =
new double*[ num_w ];
5267 for (
int w = 0; w < num_w; w++ ){
5269 FAD[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
5270 FCD[ irrep_left ][ irrep_right ][ w ] =
new double[ SIZE_left * SIZE_right ];
5272 double * AD1ptr = FAD[ irrep_left ][ irrep_right ][ w ];
5273 double * AD2ptr = FAD[ irrep_left ][ irrep_right ][ w ] + SIZE_left * D2JUMP;
5274 double * CD1ptr = FCD[ irrep_left ][ irrep_right ][ w ];
5275 double * CD2ptr = FCD[ irrep_left ][ irrep_right ][ w ] + SIZE_left * D2JUMP;
5278 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
5280 const int num_t = indices->
getNDMRG( irrep_t );
5282 const int num_u = indices->
getNDMRG( irrep_u );
5284 const int jump_AD1 = jump_AC_active( indices, irrep_w, irrep_t, irrep_u );
5285 const int jump_AD2 = jump_AC_active( indices, irrep_t, irrep_w, irrep_u );
5286 const int jump_CD2 = jump_AC_active( indices, irrep_u, irrep_t, irrep_w );
5288 for (
int t = 0; t < num_t; t++ ){
5289 for (
int u = 0; u < num_u; u++ ){
5290 for (
int xyz = 0; xyz < SIZE_left; xyz++ ){
5291 AD1ptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = + SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AD1 + w + num_w * ( t + num_t * u )) ];
5292 AD2ptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = + SAA[ irrep_left ][ xyz + SIZE_left * ( jump_AD2 + t + num_t * ( w + num_w * u )) ];
5293 CD1ptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = - SCC[ irrep_left ][ xyz + SIZE_left * ( jump_AD1 + w + num_w * ( t + num_t * u )) ];
5294 CD2ptr[ xyz + SIZE_left * ( jump_col + t + num_t * u ) ] = - SCC[ irrep_left ][ xyz + SIZE_left * ( jump_CD2 + u + num_u * ( t + num_t * w )) ];
5300 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
5302 const int num_x = indices->
getNDMRG( irrep_x );
5303 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
5306 const int num_y = indices->
getNDMRG( irrep_y );
5308 const int num_z = indices->
getNDMRG( irrep_z );
5309 assert( jump_row == jump_AC_active( indices, irrep_x, irrep_y, irrep_z ) );
5311 if ( irrep_t == irrep_w ){
5313 for (
int u = 0; u < num_u; u++ ){
5314 for (
int x = 0; x < num_x; x++ ){
5315 for (
int y = 0; y < num_y; y++ ){
5316 for (
int z = 0; z < num_z; z++ ){
5317 const double gamma_zxyu = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_u + u ))) ];
5318 CD2ptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + w + num_t * u ) ] += 2 * gamma_zxyu;
5325 if ( irrep_x == irrep_y ){
5326 for (
int u = 0; u < num_u; u++ ){
5327 for (
int z = 0; z < num_z; z++ ){
5328 const double gamma_zu = one_rdm[ d_z + z + LAS * ( d_u + u ) ];
5329 for (
int xy = 0; xy < num_x; xy++ ){
5330 CD2ptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + w + num_t * u ) ] += 2 * gamma_zu;
5337 if ( irrep_t == irrep_u ){
5339 for (
int tu = 0; tu < num_u; tu++ ){
5340 for (
int x = 0; x < num_x; x++ ){
5341 for (
int y = 0; y < num_y; y++ ){
5342 for (
int z = 0; z < num_z; z++ ){
5343 const double gamma_zxyw = two_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_y + y + LAS * ( d_w + w ))) ];
5344 CD2ptr[ jump_row + x + num_x * ( y + num_y * z ) + SIZE_left * ( jump_col + tu + num_u * tu ) ] += gamma_zxyw;
5351 if ( irrep_x == irrep_y ){
5352 for (
int z = 0; z < num_z; z++ ){
5353 const double gamma_zw = one_rdm[ d_z + z + LAS * ( d_w + w ) ];
5354 for (
int tu = 0; tu < num_u; tu++ ){
5355 for (
int xy = 0; xy < num_x; xy++ ){
5356 CD2ptr[ jump_row + xy + num_x * ( xy + num_x * z ) + SIZE_left * ( jump_col + tu + num_u * tu ) ] += gamma_zw;
5362 jump_row += num_x * num_y * num_z;
5365 jump_col += num_t * num_u;
5373 void CheMPS2::CASPT2::make_AA_CC(
const bool OVLP,
const double IPEA ){
5413 if ( OVLP ){ SAA =
new double*[ num_irreps ];
5414 SCC =
new double*[ num_irreps ]; }
5415 else { FAA =
new double*[ num_irreps ];
5416 FCC =
new double*[ num_irreps ]; }
5420 for (
int irrep = 0; irrep < num_irreps; irrep++ ){
5422 assert( size_A[ irrep ] == size_C[ irrep ] );
5423 const int SIZE = size_A[ irrep ];
5424 if ( OVLP ){ SAA[ irrep ] =
new double[ SIZE * SIZE ];
5425 SCC[ irrep ] =
new double[ SIZE * SIZE ]; }
5426 else { FAA[ irrep ] =
new double[ SIZE * SIZE ];
5427 FCC[ irrep ] =
new double[ SIZE * SIZE ]; }
5430 for (
int irrep_t = 0; irrep_t < num_irreps; irrep_t++ ){
5432 const int num_t = indices->
getNDMRG( irrep_t );
5433 const int nocc_t = indices->
getNOCC( irrep_t );
5434 for (
int irrep_u = 0; irrep_u < num_irreps; irrep_u++ ){
5436 const int num_u = indices->
getNDMRG( irrep_u );
5437 const int nocc_u = indices->
getNOCC( irrep_u );
5440 const int num_v = indices->
getNDMRG( irrep_v );
5442 for (
int irrep_x = 0; irrep_x < num_irreps; irrep_x++ ){
5444 const int num_x = indices->
getNDMRG( irrep_x );
5445 const int nocc_x = indices->
getNOCC( irrep_x );
5446 for (
int irrep_y = 0; irrep_y < num_irreps; irrep_y++ ){
5448 const int num_y = indices->
getNDMRG( irrep_y );
5449 const int nocc_y = indices->
getNOCC( irrep_y );
5452 const int num_z = indices->
getNDMRG( irrep_z );
5457 #pragma omp parallel for schedule(static) 5458 for (
int v = 0; v < num_v; v++ ){
5459 for (
int u = 0; u < num_u; u++ ){
5460 for (
int t = 0; t < num_t; t++ ){
5461 for (
int z = 0; z < num_z; z++ ){
5462 for (
int y = 0; y < num_y; y++ ){
5463 for (
int x = 0; x < num_x; x++ ){
5464 const int ptr = jump_row + x + num_x * ( y + num_y * z ) + SIZE * ( jump_col + t + num_t * ( u + num_u * v ) );
5465 SAA[ irrep ][ ptr ] = - three_rdm[ d_z + z + LAS * ( d_t + t + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_x + x + LAS * ( d_v + v ))))) ];
5473 #pragma omp parallel for schedule(static) 5474 for (
int v = 0; v < num_v; v++ ){
5475 for (
int u = 0; u < num_u; u++ ){
5476 const double f_uu = fock->
get( irrep_u, nocc_u + u, nocc_u + u );
5477 for (
int t = 0; t < num_t; t++ ){
5478 const double f_tt = fock->
get( irrep_t, nocc_t + t, nocc_t + t );
5479 for (
int z = 0; z < num_z; z++ ){
5480 for (
int y = 0; y < num_y; y++ ){
5481 const double f_yy = fock->
get( irrep_y, nocc_y + y, nocc_y + y );
5482 for (
int x = 0; x < num_x; x++ ){
5483 const double f_xx = fock->
get( irrep_x, nocc_x + x, nocc_x + x );
5484 const int ptr = jump_row + x + num_x * ( y + num_y * z ) + SIZE * ( jump_col + t + num_t * ( u + num_u * v ) );
5485 FAA[ irrep ][ ptr ] = - f_dot_4dm[ d_z + z + LAS * ( d_t + t + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_x + x + LAS * ( d_v + v ))))) ]
5486 + ( f_tt + f_uu + f_xx + f_yy ) * SAA[ irrep ][ ptr ];
5498 #pragma omp parallel for schedule(static) 5499 for (
int v = 0; v < num_v; v++ ){
5500 for (
int u = 0; u < num_u; u++ ){
5501 for (
int t = 0; t < num_t; t++ ){
5502 for (
int z = 0; z < num_z; z++ ){
5503 for (
int y = 0; y < num_y; y++ ){
5504 for (
int x = 0; x < num_x; x++ ){
5505 const int ptr = jump_row + x + num_x * ( y + num_y * z ) + SIZE * ( jump_col + t + num_t * ( u + num_u * v ) );
5506 SCC[ irrep ][ ptr ] = three_rdm[ d_z + z + LAS * ( d_x + x + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_t + t + LAS * ( d_v + v ))))) ];
5514 #pragma omp parallel for schedule(static) 5515 for (
int v = 0; v < num_v; v++ ){
5516 for (
int u = 0; u < num_u; u++ ){
5517 const double f_uu = fock->
get( irrep_u, nocc_u + u, nocc_u + u );
5518 for (
int t = 0; t < num_t; t++ ){
5519 for (
int z = 0; z < num_z; z++ ){
5520 for (
int y = 0; y < num_y; y++ ){
5521 const double f_yy = fock->
get( irrep_y, nocc_y + y, nocc_y + y );
5522 for (
int x = 0; x < num_x; x++ ){
5523 const int ptr = jump_row + x + num_x * ( y + num_y * z ) + SIZE * ( jump_col + t + num_t * ( u + num_u * v ) );
5524 FCC[ irrep ][ ptr ] = f_dot_4dm[ d_z + z + LAS * ( d_x + x + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_t + t + LAS * ( d_v + v ))))) ]
5525 + ( f_yy + f_uu ) * SCC[ irrep ][ ptr ];
5536 if ( irrep_t == irrep_x ){
5538 #pragma omp parallel for schedule(static) 5539 for (
int v = 0; v < num_v; v++ ){
5540 for (
int u = 0; u < num_u; u++ ){
5541 for (
int xt = 0; xt < num_t; xt++ ){
5542 for (
int z = 0; z < num_z; z++ ){
5543 for (
int y = 0; y < num_y; y++ ){
5544 SAA[ irrep ][ jump_row + xt + num_x * ( y + num_y * z ) + SIZE * ( jump_col + xt + num_t * ( u + num_u * v ) ) ]
5545 += 2 * two_rdm[ d_z + z + LAS * ( d_u + u + LAS * ( d_y + y + LAS * ( d_v + v ))) ];
5552 #pragma omp parallel for schedule(static) 5553 for (
int v = 0; v < num_v; v++ ){