1 /*
2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3 Copyright (C) 2013-2021 Sebastian Wouters
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <math.h>
23 #include <algorithm>
24 #include <assert.h>
25 #include <iostream>
26 #include <sstream>
27
28 #include <hdf5.h>
29
30 #include "ThreeDM.h"
31 #include "Lapack.h"
32 #include "Options.h"
33 #include "MPIchemps2.h"
34 #include "Wigner.h"
35 #include "Special.h"
36
37 using std::max;
38
ThreeDM(const SyBookkeeper * book_in,const Problem * prob_in,const bool disk_in)39 CheMPS2::ThreeDM::ThreeDM( const SyBookkeeper * book_in, const Problem * prob_in, const bool disk_in ){
40
41 book = book_in;
42 prob = prob_in;
43 disk = disk_in;
44
45 L = book->gL();
46 {
47 const long long linsize = ( long long ) L;
48 const long long size = (( disk ) ? linsize * linsize * linsize * linsize * linsize
49 : linsize * linsize * linsize * linsize * linsize * linsize );
50 assert( INT_MAX >= size );
51 array_size = size;
52 }
53
54 elements = new double[ array_size ];
55 #pragma omp simd
56 for ( int cnt = 0; cnt < array_size; cnt++ ){ elements[ cnt ] = 0.0; }
57
58 if ( disk ){
59 temp_disk_orbs = new int[ 6 * array_size ];
60 temp_disk_vals = new double [ array_size ];
61 create_file();
62 } else {
63 temp_disk_orbs = NULL;
64 temp_disk_vals = NULL;
65 }
66
67 }
68
~ThreeDM()69 CheMPS2::ThreeDM::~ThreeDM(){
70
71 delete [] elements;
72 if ( disk ){ delete [] temp_disk_orbs;
73 delete [] temp_disk_vals; }
74
75 }
76
77 #ifdef CHEMPS2_MPI_COMPILATION
mpi_allreduce()78 void CheMPS2::ThreeDM::mpi_allreduce(){
79
80 if ( disk ){
81
82 for ( int orb = 0; orb < L; orb++ ){
83 read_file( orb );
84 MPIchemps2::allreduce_array_double( elements, temp_disk_vals, array_size );
85 #pragma omp simd
86 for ( int cnt = 0; cnt < array_size; cnt++ ){ elements[ cnt ] = temp_disk_vals[ cnt ]; }
87 write_file( orb );
88 }
89
90 } else {
91
92 double * temp = new double[ array_size ];
93 MPIchemps2::allreduce_array_double( elements, temp, array_size );
94 #pragma omp simd
95 for ( int cnt = 0; cnt < array_size; cnt++ ){ elements[ cnt ] = temp[ cnt ]; }
96 delete [] temp;
97
98 }
99
100 }
101 #endif
102
set_dmrg_index(const int cnt1,const int cnt2,const int cnt3,const int cnt4,const int cnt5,const int cnt6,const double value)103 void CheMPS2::ThreeDM::set_dmrg_index( const int cnt1, const int cnt2, const int cnt3, const int cnt4, const int cnt5, const int cnt6, const double value ){
104
105 // prob assumes you use DMRG orbs, while elements is stored in hamiltonian orbitals
106 // irrep sanity checks are performed in ThreeDM::fill_site
107
108 const int orb1 = (( prob->gReorder() ) ? prob->gf2( cnt1 ) : cnt1 );
109 const int orb2 = (( prob->gReorder() ) ? prob->gf2( cnt2 ) : cnt2 );
110 const int orb3 = (( prob->gReorder() ) ? prob->gf2( cnt3 ) : cnt3 );
111 const int orb4 = (( prob->gReorder() ) ? prob->gf2( cnt4 ) : cnt4 );
112 const int orb5 = (( prob->gReorder() ) ? prob->gf2( cnt5 ) : cnt5 );
113 const int orb6 = (( prob->gReorder() ) ? prob->gf2( cnt6 ) : cnt6 );
114
115 if ( disk ){
116 int private_counter = -1;
117 #pragma omp critical
118 {
119 private_counter = temp_disk_counter;
120 temp_disk_counter++;
121 }
122 assert( private_counter < array_size );
123 temp_disk_orbs[ 6 * private_counter + 0 ] = orb1;
124 temp_disk_orbs[ 6 * private_counter + 1 ] = orb2;
125 temp_disk_orbs[ 6 * private_counter + 2 ] = orb3;
126 temp_disk_orbs[ 6 * private_counter + 3 ] = orb4;
127 temp_disk_orbs[ 6 * private_counter + 4 ] = orb5;
128 temp_disk_orbs[ 6 * private_counter + 5 ] = orb6;
129 temp_disk_vals[ private_counter ] = value;
130 return;
131 }
132
133 elements[ orb1 + L * ( orb2 + L * ( orb3 + L * ( orb4 + L * ( orb5 + L * orb6 )))) ] = value;
134 elements[ orb2 + L * ( orb3 + L * ( orb1 + L * ( orb5 + L * ( orb6 + L * orb4 )))) ] = value;
135 elements[ orb3 + L * ( orb1 + L * ( orb2 + L * ( orb6 + L * ( orb4 + L * orb5 )))) ] = value;
136 elements[ orb2 + L * ( orb1 + L * ( orb3 + L * ( orb5 + L * ( orb4 + L * orb6 )))) ] = value;
137 elements[ orb3 + L * ( orb2 + L * ( orb1 + L * ( orb6 + L * ( orb5 + L * orb4 )))) ] = value;
138 elements[ orb1 + L * ( orb3 + L * ( orb2 + L * ( orb4 + L * ( orb6 + L * orb5 )))) ] = value;
139
140 elements[ orb4 + L * ( orb5 + L * ( orb6 + L * ( orb1 + L * ( orb2 + L * orb3 )))) ] = value;
141 elements[ orb5 + L * ( orb6 + L * ( orb4 + L * ( orb2 + L * ( orb3 + L * orb1 )))) ] = value;
142 elements[ orb6 + L * ( orb4 + L * ( orb5 + L * ( orb3 + L * ( orb1 + L * orb2 )))) ] = value;
143 elements[ orb5 + L * ( orb4 + L * ( orb6 + L * ( orb2 + L * ( orb1 + L * orb3 )))) ] = value;
144 elements[ orb6 + L * ( orb5 + L * ( orb4 + L * ( orb3 + L * ( orb2 + L * orb1 )))) ] = value;
145 elements[ orb4 + L * ( orb6 + L * ( orb5 + L * ( orb1 + L * ( orb3 + L * orb2 )))) ] = value;
146
147 }
148
get_ham_index(const int cnt1,const int cnt2,const int cnt3,const int cnt4,const int cnt5,const int cnt6) const149 double CheMPS2::ThreeDM::get_ham_index( const int cnt1, const int cnt2, const int cnt3, const int cnt4, const int cnt5, const int cnt6 ) const{
150
151 assert( disk == false );
152 return elements[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * ( cnt4 + L * ( cnt5 + L * cnt6 )))) ];
153
154 }
155
fill_ham_index(const double alpha,const bool add,double * storage,const int last_orb_start,const int last_orb_num)156 void CheMPS2::ThreeDM::fill_ham_index( const double alpha, const bool add, double * storage, const int last_orb_start, const int last_orb_num ){
157
158 assert( last_orb_start >= 0 );
159 assert( last_orb_num >= 1 );
160 assert( last_orb_start + last_orb_num <= L );
161
162 if ( disk ){
163
164 for ( int ham_orb = last_orb_start; ham_orb < ( last_orb_start + last_orb_num ); ham_orb++ ){
165 read_file( ham_orb );
166 const int shift = ( ham_orb - last_orb_start ) * array_size;
167 if ( add == false ){
168 #pragma omp simd
169 for ( int cnt = 0; cnt < array_size; cnt++ ){ storage[ shift + cnt ] = alpha * elements[ cnt ]; }
170 } else {
171 #pragma omp simd
172 for ( int cnt = 0; cnt < array_size; cnt++ ){ storage[ shift + cnt ] += alpha * elements[ cnt ]; }
173 }
174 }
175
176 } else {
177
178 const int shift = last_orb_start * L * L * L * L * L;
179 const int size = last_orb_num * L * L * L * L * L;
180 if ( add == false ){
181 #pragma omp simd
182 for ( int cnt = 0; cnt < size; cnt++ ){ storage[ cnt ] = alpha * elements[ shift + cnt ]; }
183 } else {
184 #pragma omp simd
185 for ( int cnt = 0; cnt < size; cnt++ ){ storage[ cnt ] += alpha * elements[ shift + cnt ]; }
186 }
187
188 }
189
190 }
191
trace()192 double CheMPS2::ThreeDM::trace(){
193
194 double value = 0.0;
195
196 if ( disk ){
197
198 for ( int cnt3 = 0; cnt3 < L; cnt3++ ){
199 read_file( cnt3 );
200 for ( int cnt2 = 0; cnt2 < L; cnt2++ ){
201 for ( int cnt1 = 0; cnt1 < L; cnt1++ ){
202 value += elements[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * ( cnt1 + L * cnt2 ))) ];
203 }
204 }
205 }
206
207 } else {
208
209 for ( int cnt3 = 0; cnt3 < L; cnt3++ ){
210 for ( int cnt2 = 0; cnt2 < L; cnt2++ ){
211 for ( int cnt1 = 0; cnt1 < L; cnt1++ ){
212 value += get_ham_index( cnt1, cnt2, cnt3, cnt1, cnt2, cnt3 );
213 }
214 }
215 }
216
217 }
218
219 return value;
220
221 }
222
create_file() const223 void CheMPS2::ThreeDM::create_file() const{
224
225 #ifdef CHEMPS2_MPI_COMPILATION
226 const int mpi_rank = MPIchemps2::mpi_rank();
227 #else
228 const int mpi_rank = 0;
229 #endif
230
231 assert( disk == true );
232
233 std::stringstream filename;
234 filename << CheMPS2::THREE_RDM_storage_prefix << mpi_rank << ".h5";
235
236 hid_t file_id = H5Fcreate( filename.str().c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT );
237 hid_t group_id = H5Gcreate( file_id, "three_rdm", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
238
239 for ( int orb = 0; orb < L; orb++ ){
240
241 std::stringstream storagename;
242 storagename << "elements_" << orb;
243
244 hsize_t dimarray = array_size;
245 hid_t dataspace_id = H5Screate_simple( 1, &dimarray, NULL );
246 hid_t dataset_id = H5Dcreate( group_id, storagename.str().c_str(), H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
247 H5Dwrite( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, elements );
248
249 H5Dclose( dataset_id );
250 H5Sclose( dataspace_id );
251
252 }
253
254 H5Gclose( group_id );
255 H5Fclose( file_id );
256
257 }
258
write_file(const int last_ham_orb) const259 void CheMPS2::ThreeDM::write_file( const int last_ham_orb ) const{
260
261 #ifdef CHEMPS2_MPI_COMPILATION
262 const int mpi_rank = MPIchemps2::mpi_rank();
263 #else
264 const int mpi_rank = 0;
265 #endif
266
267 assert( disk == true );
268
269 std::stringstream filename;
270 filename << CheMPS2::THREE_RDM_storage_prefix << mpi_rank << ".h5";
271
272 hid_t file_id = H5Fopen( filename.str().c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
273 hid_t group_id = H5Gopen( file_id, "three_rdm", H5P_DEFAULT );
274
275 std::stringstream storagename;
276 storagename << "elements_" << last_ham_orb;
277
278 hid_t dataset_id = H5Dopen( group_id, storagename.str().c_str(), H5P_DEFAULT );
279 H5Dwrite( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, elements );
280
281 H5Dclose( dataset_id );
282
283 H5Gclose( group_id );
284 H5Fclose( file_id );
285
286 }
287
read_file(const int last_ham_orb)288 void CheMPS2::ThreeDM::read_file( const int last_ham_orb ){
289
290 #ifdef CHEMPS2_MPI_COMPILATION
291 const int mpi_rank = MPIchemps2::mpi_rank();
292 #else
293 const int mpi_rank = 0;
294 #endif
295
296 assert( disk == true );
297
298 std::stringstream filename;
299 filename << CheMPS2::THREE_RDM_storage_prefix << mpi_rank << ".h5";
300
301 hid_t file_id = H5Fopen( filename.str().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
302 hid_t group_id = H5Gopen( file_id, "three_rdm", H5P_DEFAULT );
303
304 std::stringstream storagename;
305 storagename << "elements_" << last_ham_orb;
306
307 hid_t dataset_id = H5Dopen( group_id, storagename.str().c_str(), H5P_DEFAULT );
308 H5Dread( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, elements );
309
310 H5Dclose( dataset_id );
311
312 H5Gclose( group_id );
313 H5Fclose( file_id );
314
315 }
316
correct_higher_multiplicities()317 void CheMPS2::ThreeDM::correct_higher_multiplicities(){
318
319 if ( prob->gTwoS() != 0 ){
320 double alpha = 1.0 / ( prob->gTwoS() + 1.0 );
321 int inc1 = 1;
322 if ( disk ){
323 for ( int ham_orb = 0; ham_orb < L; ham_orb++ ){
324 read_file( ham_orb );
325 dscal_( &array_size, &alpha, elements, &inc1 );
326 write_file( ham_orb );
327 }
328 } else {
329 dscal_( &array_size, &alpha, elements, &inc1 );
330 }
331 }
332
333 }
334
flush_disk()335 void CheMPS2::ThreeDM::flush_disk(){
336
337 assert( disk == true );
338 for ( int ham_orb = 0; ham_orb < L; ham_orb++ ){
339
340 read_file( ham_orb );
341
342 for ( int counter = 0; counter < temp_disk_counter; counter++ ){
343
344 const int orb1 = temp_disk_orbs[ 6 * counter + 0 ];
345 const int orb2 = temp_disk_orbs[ 6 * counter + 1 ];
346 const int orb3 = temp_disk_orbs[ 6 * counter + 2 ];
347 const int orb4 = temp_disk_orbs[ 6 * counter + 3 ];
348 const int orb5 = temp_disk_orbs[ 6 * counter + 4 ];
349 const int orb6 = temp_disk_orbs[ 6 * counter + 5 ];
350 const double value = temp_disk_vals[ counter ];
351
352 if ( orb1 == ham_orb ){ elements[ orb5 + L * ( orb6 + L * ( orb4 + L * ( orb2 + L * orb3 ))) ] = value;
353 elements[ orb6 + L * ( orb5 + L * ( orb4 + L * ( orb3 + L * orb2 ))) ] = value; }
354 if ( orb2 == ham_orb ){ elements[ orb6 + L * ( orb4 + L * ( orb5 + L * ( orb3 + L * orb1 ))) ] = value;
355 elements[ orb4 + L * ( orb6 + L * ( orb5 + L * ( orb1 + L * orb3 ))) ] = value; }
356 if ( orb3 == ham_orb ){ elements[ orb4 + L * ( orb5 + L * ( orb6 + L * ( orb1 + L * orb2 ))) ] = value;
357 elements[ orb5 + L * ( orb4 + L * ( orb6 + L * ( orb2 + L * orb1 ))) ] = value; }
358 if ( orb4 == ham_orb ){ elements[ orb2 + L * ( orb3 + L * ( orb1 + L * ( orb5 + L * orb6 ))) ] = value;
359 elements[ orb3 + L * ( orb2 + L * ( orb1 + L * ( orb6 + L * orb5 ))) ] = value; }
360 if ( orb5 == ham_orb ){ elements[ orb3 + L * ( orb1 + L * ( orb2 + L * ( orb6 + L * orb4 ))) ] = value;
361 elements[ orb1 + L * ( orb3 + L * ( orb2 + L * ( orb4 + L * orb6 ))) ] = value; }
362 if ( orb6 == ham_orb ){ elements[ orb1 + L * ( orb2 + L * ( orb3 + L * ( orb4 + L * orb5 ))) ] = value;
363 elements[ orb2 + L * ( orb1 + L * ( orb3 + L * ( orb5 + L * orb4 ))) ] = value; }
364 }
365
366 write_file( ham_orb );
367
368 }
369
370 }
371
save_HAM(const string filename) const372 void CheMPS2::ThreeDM::save_HAM( const string filename ) const{
373
374 assert( disk == false );
375 save_HAM_generic( filename, L, "3-RDM", elements );
376
377 }
378
save_HAM_generic(const string filename,const int LAS,const string tag,double * array)379 void CheMPS2::ThreeDM::save_HAM_generic( const string filename, const int LAS, const string tag, double * array ){
380
381 hid_t file_id = H5Fcreate( filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT );
382 long long linsize = ( long long ) LAS;
383 hsize_t dimarray = ( linsize * linsize * linsize * linsize * linsize * linsize );
384 hid_t group_id = H5Gcreate( file_id, tag.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
385 hid_t dataspace_id = H5Screate_simple( 1, &dimarray, NULL );
386 hid_t dataset_id = H5Dcreate( group_id, "elements", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
387
388 H5Dwrite( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, array );
389
390 H5Dclose( dataset_id );
391 H5Sclose( dataspace_id );
392 H5Gclose( group_id );
393 H5Fclose( file_id );
394
395 std::cout << "Saved the " << tag << " to the file " << filename << std::endl;
396
397 }
398
fill_site(TensorT * denT,TensorL *** Ltensors,TensorF0 **** F0tensors,TensorF1 **** F1tensors,TensorS0 **** S0tensors,TensorS1 **** S1tensors,Tensor3RDM **** dm3_a_J0_doublet,Tensor3RDM **** dm3_a_J1_doublet,Tensor3RDM **** dm3_a_J1_quartet,Tensor3RDM **** dm3_b_J0_doublet,Tensor3RDM **** dm3_b_J1_doublet,Tensor3RDM **** dm3_b_J1_quartet,Tensor3RDM **** dm3_c_J0_doublet,Tensor3RDM **** dm3_c_J1_doublet,Tensor3RDM **** dm3_c_J1_quartet,Tensor3RDM **** dm3_d_J0_doublet,Tensor3RDM **** dm3_d_J1_doublet,Tensor3RDM **** dm3_d_J1_quartet)399 void CheMPS2::ThreeDM::fill_site( TensorT * denT, TensorL *** Ltensors, TensorF0 **** F0tensors, TensorF1 **** F1tensors, TensorS0 **** S0tensors, TensorS1 **** S1tensors,
400 Tensor3RDM **** dm3_a_J0_doublet, Tensor3RDM **** dm3_a_J1_doublet, Tensor3RDM **** dm3_a_J1_quartet,
401 Tensor3RDM **** dm3_b_J0_doublet, Tensor3RDM **** dm3_b_J1_doublet, Tensor3RDM **** dm3_b_J1_quartet,
402 Tensor3RDM **** dm3_c_J0_doublet, Tensor3RDM **** dm3_c_J1_doublet, Tensor3RDM **** dm3_c_J1_quartet,
403 Tensor3RDM **** dm3_d_J0_doublet, Tensor3RDM **** dm3_d_J1_doublet, Tensor3RDM **** dm3_d_J1_quartet ){
404
405 #ifdef CHEMPS2_MPI_COMPILATION
406 const int MPIRANK = MPIchemps2::mpi_rank();
407 #endif
408
409 temp_disk_counter = 0;
410 const int orb_i = denT->gIndex();
411 const int DIM = max(book->gMaxDimAtBound( orb_i ), book->gMaxDimAtBound( orb_i+1 ));
412 const double sq3 = sqrt( 3.0 );
413
414 #pragma omp parallel
415 {
416
417 double * workmem = new double[ DIM * DIM ];
418 double * workmem2 = new double[ DIM * DIM ];
419
420 const int upperbound1 = ( orb_i * ( orb_i + 1 )) / 2;
421 int jkl[] = { 0, 0, 0 };
422 #ifdef CHEMPS2_MPI_COMPILATION
423 #pragma omp for schedule(dynamic) nowait
424 #else
425 #pragma omp for schedule(static) nowait
426 #endif
427 for ( int global = 0; global < upperbound1; global++ ){
428 Special::invert_triangle_two( global, jkl );
429 const int orb_j = jkl[ 0 ];
430 const int orb_k = jkl[ 1 ];
431 if ( book->gIrrep( orb_j ) == book->gIrrep( orb_k )){
432 #ifdef CHEMPS2_MPI_COMPILATION
433 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) )
434 #endif
435 {
436 const double d1 = diagram1( denT, F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], workmem );
437 set_dmrg_index( orb_j, orb_i, orb_i, orb_k, orb_i, orb_i, 4 * d1 );
438 set_dmrg_index( orb_j, orb_i, orb_i, orb_i, orb_k, orb_i, -2 * d1 );
439 }
440 }
441 }
442
443 const int triangle3 = L - orb_i - 1;
444 const int upperbound3 = ( triangle3 * ( triangle3 + 1 ) ) / 2;
445 #ifdef CHEMPS2_MPI_COMPILATION
446 #pragma omp for schedule(dynamic) nowait
447 #else
448 #pragma omp for schedule(static) nowait
449 #endif
450 for ( int global = 0; global < upperbound3; global++ ){
451 Special::invert_triangle_two( global, jkl );
452 const int orb_j = L - 1 - jkl[ 1 ];
453 const int orb_k = orb_j + jkl[ 0 ];
454 if ( book->gIrrep( orb_j ) == book->gIrrep( orb_k )){
455 #ifdef CHEMPS2_MPI_COMPILATION
456 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) )
457 #endif
458 {
459 const double d3 = diagram3( denT, F0tensors[orb_i][orb_k-orb_j][orb_j-1-orb_i], workmem );
460 set_dmrg_index( orb_j, orb_i, orb_i, orb_k, orb_i, orb_i, 4 * d3 );
461 set_dmrg_index( orb_j, orb_i, orb_i, orb_i, orb_k, orb_i, -2 * d3 );
462 }
463 }
464 }
465
466 const int upperbound4 = ( orb_i * ( orb_i + 1 ) * ( orb_i + 2 ) ) / 6;
467 #ifdef CHEMPS2_MPI_COMPILATION
468 #pragma omp for schedule(dynamic) nowait
469 #else
470 #pragma omp for schedule(static) nowait
471 #endif
472 for ( int global = 0; global < upperbound4; global++ ){
473 Special::invert_triangle_three( global, jkl );
474 const int orb_j = jkl[ 0 ];
475 const int orb_k = jkl[ 1 ];
476 const int orb_l = jkl[ 2 ];
477 const int recalculate_global = orb_j + ( orb_k * ( orb_k + 1 ) ) / 2 + ( orb_l * ( orb_l + 1 ) * ( orb_l + 2 ) ) / 6;
478 assert( global == recalculate_global );
479 if ( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ) == Irreps::directProd( book->gIrrep( orb_l ), book->gIrrep( orb_i ) ) ){
480 #ifdef CHEMPS2_MPI_COMPILATION
481 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
482 #endif
483 {
484 const int cnt1 = orb_k - orb_j;
485 const int cnt2 = orb_l - orb_k;
486 const int cnt3 = orb_i - 1 - orb_l;
487 if ( cnt2 > 0 ){
488 const double d4 = diagram4_5_6_7_8_9( denT, dm3_b_J0_doublet[cnt1][cnt2][cnt3], workmem, 'B' );
489 const double d5 = (( cnt1 == 0 ) ? 0.0 : diagram4_5_6_7_8_9( denT, dm3_b_J1_doublet[cnt1][cnt2][cnt3], workmem, 'B' ));
490 set_dmrg_index( orb_j, orb_k, orb_i, orb_l, orb_i, orb_i, d4 + d5 );
491 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_l, orb_i, d4 - d5 );
492 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_i, orb_l, -2 * d4 );
493 }
494 if ( cnt1 + cnt2 > 0 ){
495 const double d6 = diagram4_5_6_7_8_9( denT, dm3_c_J0_doublet[cnt1][cnt2][cnt3], workmem, 'C' );
496 const double d7 = diagram4_5_6_7_8_9( denT, dm3_c_J1_doublet[cnt1][cnt2][cnt3], workmem, 'C' );
497 set_dmrg_index( orb_j, orb_l, orb_i, orb_k, orb_i, orb_i, -2 * d6 );
498 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_k, orb_i, d6 - d7 );
499 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_i, orb_k, d6 + d7 );
500 }
501 {
502 const double d8 = diagram4_5_6_7_8_9( denT, dm3_d_J0_doublet[cnt1][cnt2][cnt3], workmem, 'D' );
503 const double d9 = diagram4_5_6_7_8_9( denT, dm3_d_J1_doublet[cnt1][cnt2][cnt3], workmem, 'D' );
504 set_dmrg_index( orb_k, orb_l, orb_i, orb_j, orb_i, orb_i, -2 * d8 );
505 set_dmrg_index( orb_k, orb_l, orb_i, orb_i, orb_j, orb_i, d8 + d9 );
506 set_dmrg_index( orb_k, orb_l, orb_i, orb_i, orb_i, orb_j, d8 - d9 );
507 }
508 }
509 }
510 }
511
512 const int upperbound10 = upperbound1 * ( L - orb_i - 1 );
513 #ifdef CHEMPS2_MPI_COMPILATION
514 #pragma omp for schedule(dynamic) nowait
515 #else
516 #pragma omp for schedule(static) nowait
517 #endif
518 for (int combined = 0; combined < upperbound10; combined++){
519 const int global = combined % upperbound1;
520 Special::invert_triangle_two( global, jkl );
521 const int orb_l = orb_i + 1 + ( combined / upperbound1 );
522 const int orb_j = jkl[ 0 ];
523 const int orb_k = jkl[ 1 ];
524 if ( Irreps::directProd(book->gIrrep( orb_j ), book->gIrrep( orb_k )) == Irreps::directProd(book->gIrrep( orb_l ), book->gIrrep( orb_i )) ){
525 #ifdef CHEMPS2_MPI_COMPILATION
526 if ( MPIRANK == MPIchemps2::owner_absigma( orb_j, orb_k ) )
527 #endif
528 {
529 const double d10 = diagram10( denT, S0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 );
530 const double d11 = (( orb_j == orb_k ) ? 0.0 :
531 diagram11( denT, S1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 ));
532 set_dmrg_index( orb_j, orb_k, orb_i, orb_l, orb_i, orb_i, d10 + d11 );
533 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_l, orb_i, d10 - d11 );
534 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_i, orb_l, - 2 * d10 );
535 }
536
537 #ifdef CHEMPS2_MPI_COMPILATION
538 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) )
539 #endif
540 {
541 {
542 const double d12 = diagram12( denT, F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 );
543 const double d13 = diagram13( denT, F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 );
544 set_dmrg_index( orb_j, orb_l, orb_i, orb_k, orb_i, orb_i, - 2 * d12 );
545 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_k, orb_i, d12 - d13 );
546 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_i, orb_k, d12 + d13 );
547 }
548 if ( orb_j < orb_k ){
549 const double d14 = diagram14( denT, F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 );
550 const double d15 = diagram15( denT, F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], Ltensors[orb_i][orb_l-1-orb_i], workmem, workmem2 );
551 set_dmrg_index( orb_k, orb_l, orb_i, orb_j, orb_i, orb_i, - 2 * d14 );
552 set_dmrg_index( orb_k, orb_l, orb_i, orb_i, orb_j, orb_i, d14 - d15 );
553 set_dmrg_index( orb_k, orb_l, orb_i, orb_i, orb_i, orb_j, d14 + d15 );
554 }
555 }
556 }
557 }
558
559 const int upperbound16 = upperbound3 * orb_i;
560 #ifdef CHEMPS2_MPI_COMPILATION
561 #pragma omp for schedule(dynamic) nowait
562 #else
563 #pragma omp for schedule(static) nowait
564 #endif
565 for ( int combined = 0; combined < upperbound16; combined++ ){
566 const int orb_j = combined % orb_i;
567 const int global = combined / orb_i;
568 Special::invert_triangle_two( global, jkl );
569 const int orb_m = L - 1 - jkl[ 1 ];
570 const int orb_n = orb_m + jkl[ 0 ];
571 if ( Irreps::directProd(book->gIrrep( orb_j ), book->gIrrep( orb_i )) == Irreps::directProd(book->gIrrep( orb_m ), book->gIrrep( orb_n )) ){
572 #ifdef CHEMPS2_MPI_COMPILATION
573 if ( MPIRANK == MPIchemps2::owner_absigma( orb_m, orb_n ) )
574 #endif
575 {
576 const double d16 = diagram16( denT, Ltensors[orb_i-1][orb_i-1-orb_j], S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
577 const double d17 = (( orb_m == orb_n ) ? 0.0 :
578 diagram17( denT, Ltensors[orb_i-1][orb_i-1-orb_j], S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 ));
579 set_dmrg_index( orb_m, orb_n, orb_i, orb_j, orb_i, orb_i, d16 - d17 );
580 set_dmrg_index( orb_m, orb_n, orb_i, orb_i, orb_j, orb_i, d16 + d17 );
581 set_dmrg_index( orb_m, orb_n, orb_i, orb_i, orb_i, orb_j, - 2 * d16 );
582 }
583
584 #ifdef CHEMPS2_MPI_COMPILATION
585 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_m, orb_n ) )
586 #endif
587 {
588 {
589 const double d18 = diagram18( denT, Ltensors[orb_i-1][orb_i-1-orb_j], F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
590 const double d19 = diagram19( denT, Ltensors[orb_i-1][orb_i-1-orb_j], F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
591 set_dmrg_index( orb_j, orb_m, orb_i, orb_n, orb_i, orb_i, d18 + d19 );
592 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_n, orb_i, - 2 * d18 );
593 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_i, orb_n, d18 - d19 );
594 }
595 if ( orb_m < orb_n ){
596 const double d20 = diagram20( denT, Ltensors[orb_i-1][orb_i-1-orb_j], F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
597 const double d21 = diagram21( denT, Ltensors[orb_i-1][orb_i-1-orb_j], F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
598 set_dmrg_index( orb_j, orb_n, orb_i, orb_m, orb_i, orb_i, d20 + d21 );
599 set_dmrg_index( orb_j, orb_n, orb_i, orb_i, orb_m, orb_i, - 2 * d20 );
600 set_dmrg_index( orb_j, orb_n, orb_i, orb_i, orb_i, orb_m, d20 - d21 );
601 }
602 }
603 }
604 }
605
606 for ( int orb_m = orb_i+1; orb_m < L; orb_m++ ){
607 for ( int orb_n = orb_m; orb_n < L; orb_n++ ){
608
609 /*
610 * Strategy for MPI of partitioning 2-2-2 and 3-1-2
611 * - owner of left renormalized operator is going to calculate
612 * - outer loop is (m,n)
613 * - in (m,n) loop make a temporary duplicate of S_mn & F_mn
614 */
615
616 const int irrep_mn = Irreps::directProd( book->gIrrep( orb_m ), book->gIrrep( orb_n ) );
617 const int irrep_imn = Irreps::directProd( book->gIrrep( orb_i ), irrep_mn );
618
619 /************************************************************
620 * Make sure every process has the relevant S_mn and F_mn *
621 ************************************************************/
622 #ifdef CHEMPS2_MPI_COMPILATION
623 #pragma omp barrier // Everyone needs to be done before tensors are created and communicated
624 #pragma omp single
625 if ( orb_m > orb_i + 1 ){ // All processes own Fx/Sx[ index ][ n - m ][ m - i - 1 == 0 ]
626 const int own_S_mn = MPIchemps2::owner_absigma( orb_m, orb_n );
627 const int own_F_mn = MPIchemps2::owner_cdf( L, orb_m, orb_n );
628 if ( MPIRANK != own_F_mn ){ F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = new TensorF0( orb_i+1, irrep_mn, false, book );
629 F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = new TensorF1( orb_i+1, irrep_mn, false, book ); }
630 if ( MPIRANK != own_S_mn ){ S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = new TensorS0( orb_i+1, irrep_mn, false, book );
631 if ( orb_m != orb_n ){ S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = new TensorS1( orb_i+1, irrep_mn, false, book ); }}
632 MPIchemps2::broadcast_tensor( F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], own_F_mn );
633 MPIchemps2::broadcast_tensor( F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], own_F_mn );
634 MPIchemps2::broadcast_tensor( S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], own_S_mn );
635 if ( orb_m != orb_n ){ MPIchemps2::broadcast_tensor( S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], own_S_mn ); }
636 }
637 #endif
638
639 /**************************************************************
640 * 2-2-2 : Find out how many contributions each process has *
641 **************************************************************/
642 int counter_Fjk = 0;
643 int counter_Sjk = 0;
644 for (int global = 0; global < upperbound1; global++){
645 Special::invert_triangle_two( global, jkl );
646 const int orb_j = jkl[ 0 ];
647 const int orb_k = jkl[ 1 ];
648 const int irrep_jk = Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) );
649 if ( irrep_jk == irrep_mn ){
650 #ifdef CHEMPS2_MPI_COMPILATION
651 if ( MPIRANK == MPIchemps2::owner_absigma( orb_j, orb_k ) ){ counter_Sjk++; }
652 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) ){ counter_Fjk++; }
653 #else
654 counter_Fjk++;
655 counter_Sjk++;
656 #endif
657 }
658 }
659
660 /*******************************
661 * 2-2-2 : contributions F-F *
662 *******************************/
663 if ( counter_Fjk > 0 ){
664
665 TensorF0 * tens_29_33 = new TensorF0( orb_i, irrep_mn, true, book );
666 TensorF1 * tens_30_32 = new TensorF1( orb_i, irrep_mn, true, book );
667 TensorF1 * tens_34_40 = new TensorF1( orb_i, irrep_mn, true, book );
668 TensorF0 * tens_35_41 = new TensorF0( orb_i, irrep_mn, true, book );
669 TensorF1 * tens_36_42 = new TensorF1( orb_i, irrep_mn, true, book );
670 TensorF1 * tens_37_44 = new TensorF1( orb_i, irrep_mn, true, book );
671 TensorF1 * tens_38_43 = new TensorF1( orb_i, irrep_mn, true, book );
672 fill_tens_29_33( denT, tens_29_33, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
673 fill_tens_30_32( denT, tens_30_32, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
674 fill_tens_36_42( denT, tens_36_42, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
675 fill_tens_34_35_37_38( denT, tens_34_40, tens_35_41, tens_37_44, tens_38_43, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
676
677 #ifdef CHEMPS2_MPI_COMPILATION
678 #pragma omp for schedule(dynamic) nowait
679 #else
680 #pragma omp for schedule(static) nowait
681 #endif
682 for ( int global = 0; global < upperbound1; global++ ){
683 Special::invert_triangle_two( global, jkl );
684 const int orb_j = jkl[ 0 ];
685 const int orb_k = jkl[ 1 ];
686 const int irrep_jk = Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) );
687 if ( irrep_jk == irrep_mn ){
688 #ifdef CHEMPS2_MPI_COMPILATION
689 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) )
690 #endif
691 {
692 if ( orb_m < orb_n ){
693 const double d31 = tens_29_33->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
694 const double d32 = tens_30_32->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
695 const double d42 = tens_36_42->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
696 const double d40 = tens_34_40->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
697 const double d41 = tens_35_41->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
698 const double d44 = tens_37_44->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
699 const double d43 = tens_38_43->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
700 set_dmrg_index( orb_j, orb_n, orb_i, orb_k, orb_m, orb_i, 4*d31 );
701 set_dmrg_index( orb_j, orb_n, orb_i, orb_m, orb_k, orb_i, -2*(d31 + d32 + d40) );
702 set_dmrg_index( orb_j, orb_n, orb_i, orb_k, orb_i, orb_m, -2*(d31 + d41) );
703 set_dmrg_index( orb_j, orb_n, orb_i, orb_m, orb_i, orb_k, d31 + d32 + d41 + d42 + d43 );
704 set_dmrg_index( orb_j, orb_n, orb_i, orb_i, orb_k, orb_m, d31 + d32 + d41 + d42 + d44 );
705 set_dmrg_index( orb_j, orb_n, orb_i, orb_i, orb_m, orb_k, -2*(d31 + d42) );
706 }
707 const double d29 = tens_29_33->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
708 const double d30 = tens_30_32->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
709 const double d36 = tens_36_42->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
710 const double d34 = tens_34_40->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
711 const double d35 = tens_35_41->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
712 const double d37 = tens_37_44->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
713 const double d38 = tens_38_43->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
714 set_dmrg_index( orb_j, orb_m, orb_i, orb_k, orb_n, orb_i, 4*d29 );
715 set_dmrg_index( orb_j, orb_m, orb_i, orb_n, orb_k, orb_i, -2*(d29 + d30 + d34) );
716 set_dmrg_index( orb_j, orb_m, orb_i, orb_k, orb_i, orb_n, -2*(d29 + d35) );
717 set_dmrg_index( orb_j, orb_m, orb_i, orb_n, orb_i, orb_k, d29 + d30 + d35 + d36 + d37 );
718 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_k, orb_n, d29 + d30 + d35 + d36 + d38 );
719 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_n, orb_k, -2*(d29 + d36) );
720 }
721 }
722 }
723
724 delete tens_29_33;
725 delete tens_30_32;
726 delete tens_34_40;
727 delete tens_35_41;
728 delete tens_36_42;
729 delete tens_37_44;
730 delete tens_38_43;
731
732 }
733
734 /***********************************
735 * 2-2-2 : contributions F-Sigma *
736 ***********************************/
737 if ( counter_Fjk > 0 ){
738
739 TensorF0 * tens_49_51 = new TensorF0( orb_i, irrep_mn, true, book );
740 TensorF1 * tens_50_52 = (( orb_m == orb_n ) ? NULL : new TensorF1( orb_i, irrep_mn, true, book ));
741 fill_tens_49_51( denT, tens_49_51, S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
742 if ( orb_m != orb_n ){ fill_tens_50_52( denT, tens_50_52, S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem ); }
743
744 #ifdef CHEMPS2_MPI_COMPILATION
745 #pragma omp for schedule(dynamic) nowait
746 #else
747 #pragma omp for schedule(static) nowait
748 #endif
749 for ( int global = 0; global < upperbound1; global++ ){
750 Special::invert_triangle_two( global, jkl );
751 const int orb_j = jkl[ 0 ];
752 const int orb_k = jkl[ 1 ];
753 const int irrep_jk = Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) );
754 if ( irrep_jk == irrep_mn ){
755 #ifdef CHEMPS2_MPI_COMPILATION
756 if ( MPIRANK == MPIchemps2::owner_cdf( L, orb_j, orb_k ) )
757 #endif
758 {
759 if ( orb_j < orb_k ){
760 const double d51 = tens_49_51->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' );
761 const double d52 = ((orb_m==orb_n)?0.0 : tens_50_52->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'N' ));
762 set_dmrg_index( orb_k, orb_m, orb_n, orb_j, orb_i, orb_i, - 2 * d51 );
763 set_dmrg_index( orb_k, orb_m, orb_n, orb_i, orb_j, orb_i, d51 + d52 );
764 set_dmrg_index( orb_k, orb_m, orb_n, orb_i, orb_i, orb_j, d51 - d52 );
765 }
766 const double d49 = tens_49_51->inproduct( F0tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' );
767 const double d50 = ((orb_m==orb_n)?0.0 : tens_50_52->inproduct( F1tensors[orb_i-1][orb_k-orb_j][orb_i-1-orb_k], 'T' ));
768 set_dmrg_index( orb_j, orb_m, orb_n, orb_k, orb_i, orb_i, - 2 * d49 );
769 set_dmrg_index( orb_j, orb_m, orb_n, orb_i, orb_k, orb_i, d49 + d50 );
770 set_dmrg_index( orb_j, orb_m, orb_n, orb_i, orb_i, orb_k, d49 - d50 );
771 }
772 }
773 }
774
775 delete tens_49_51;
776 if ( orb_m != orb_n ){ delete tens_50_52; }
777
778 }
779
780 /***************************************
781 * 2-2-2 : contributions Sigma-Sigma *
782 ***************************************/
783 if ( counter_Sjk > 0 ){
784
785 TensorS0 * tens_22 = new TensorS0( orb_i, irrep_mn, true, book );
786 TensorS1 * tens_23 = (( orb_m == orb_n ) ? NULL : new TensorS1( orb_i, irrep_mn, true, book ));
787 TensorS1 * tens_25 = (( orb_m == orb_n ) ? NULL : new TensorS1( orb_i, irrep_mn, true, book ));
788 TensorS1 * tens_26 = (( orb_m == orb_n ) ? NULL : new TensorS1( orb_i, irrep_mn, true, book ));
789 TensorS0 * tens_27 = (( orb_m == orb_n ) ? NULL : new TensorS0( orb_i, irrep_mn, true, book ));
790 TensorS1 * tens_28 = new TensorS1( orb_i, irrep_mn, true, book );
791 fill_tens_22_24( denT, tens_22, S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
792 fill_tens_28( denT, tens_28, S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
793 if ( orb_m != orb_n ){ fill_tens_23( denT, tens_23, S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
794 fill_tens_25_26_27( denT, tens_25, tens_26, tens_27, S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 ); }
795
796 #ifdef CHEMPS2_MPI_COMPILATION
797 #pragma omp for schedule(dynamic) nowait
798 #else
799 #pragma omp for schedule(static) nowait
800 #endif
801 for (int global = 0; global < upperbound1; global++){
802 Special::invert_triangle_two( global, jkl );
803 const int orb_j = jkl[ 0 ];
804 const int orb_k = jkl[ 1 ];
805 const int irrep_jk = Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) );
806 if ( irrep_jk == irrep_mn ){
807 #ifdef CHEMPS2_MPI_COMPILATION
808 if ( MPIRANK == MPIchemps2::owner_absigma( orb_j, orb_k ) )
809 #endif
810 {
811 const int cnt1 = orb_k - orb_j;
812 const double d22 = tens_22->inproduct(S0tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N');
813 const double d23 = (((cnt1==0)||(orb_m==orb_n))?0.0: tens_23->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
814 const double d25 = (((cnt1==0)||(orb_m==orb_n))?0.0: tens_25->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
815 const double d26 = (((cnt1==0)||(orb_m==orb_n))?0.0: tens_26->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
816 const double d27 = ( (orb_m==orb_n) ?0.0: tens_27->inproduct(S0tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
817 const double d28 = ( (cnt1==0) ?0.0: tens_28->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
818 set_dmrg_index( orb_j, orb_k, orb_i, orb_m, orb_n, orb_i, 2*d22 + 2*d23 + d25 );
819 set_dmrg_index( orb_j, orb_k, orb_i, orb_n, orb_m, orb_i, 2*d22 - 2*d23 - d25 );
820 set_dmrg_index( orb_j, orb_k, orb_i, orb_m, orb_i, orb_n, -d22 - d23 + d26 + d27 + d28 );
821 set_dmrg_index( orb_j, orb_k, orb_i, orb_n, orb_i, orb_m, -d22 + d23 - d26 - d27 + d28 );
822 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_m, orb_n, -d22 + d23 - d26 + d27 - d28 );
823 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_n, orb_m, -d22 - d23 + d26 - d27 - d28 );
824 }
825 }
826 }
827
828 delete tens_22;
829 if ( orb_m != orb_n ){ delete tens_23;
830 delete tens_25;
831 delete tens_26;
832 delete tens_27; }
833 delete tens_28;
834
835 }
836
837 /***********************************
838 * 2-2-2 : contributions Sigma-F *
839 ***********************************/
840 if ( counter_Sjk > 0 ){
841
842 TensorS0 * tens_45 = new TensorS0( orb_i, irrep_mn, true, book );
843 TensorS1 * tens_46 = new TensorS1( orb_i, irrep_mn, true, book );
844 TensorS0 * tens_47 = (( orb_m == orb_n ) ? NULL : new TensorS0( orb_i, irrep_mn, true, book ));
845 TensorS1 * tens_48 = (( orb_m == orb_n ) ? NULL : new TensorS1( orb_i, irrep_mn, true, book ));
846 fill_tens_45_47( denT, tens_45, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, true );
847 fill_tens_46_48( denT, tens_46, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, true );
848 if ( orb_m != orb_n ){ fill_tens_45_47( denT, tens_47, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, false );
849 fill_tens_46_48( denT, tens_48, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, false ); }
850
851
852 #ifdef CHEMPS2_MPI_COMPILATION
853 #pragma omp for schedule(dynamic) nowait
854 #else
855 #pragma omp for schedule(static) nowait
856 #endif
857 for (int global = 0; global < upperbound1; global++){
858 Special::invert_triangle_two( global, jkl );
859 const int orb_j = jkl[ 0 ];
860 const int orb_k = jkl[ 1 ];
861 const int irrep_jk = Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) );
862 if ( irrep_jk == irrep_mn ){
863 #ifdef CHEMPS2_MPI_COMPILATION
864 if ( MPIRANK == MPIchemps2::owner_absigma( orb_j, orb_k ) )
865 #endif
866 {
867 const int cnt1 = orb_k - orb_j;
868 if ( orb_m < orb_n ){
869 const double d47 = tens_47->inproduct(S0tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N');
870 const double d48 = (( cnt1 == 0 ) ? 0.0 : tens_48->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
871 set_dmrg_index( orb_j, orb_k, orb_n, orb_m, orb_i, orb_i, d47 + d48 );
872 set_dmrg_index( orb_j, orb_k, orb_n, orb_i, orb_m, orb_i, d47 - d48 );
873 set_dmrg_index( orb_j, orb_k, orb_n, orb_i, orb_i, orb_m, -2 * d47 );
874 }
875 const double d45 = tens_45->inproduct(S0tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N');
876 const double d46 = (( cnt1 == 0 ) ? 0.0: tens_46->inproduct(S1tensors[orb_i-1][cnt1][orb_i-1-orb_k],'N'));
877 set_dmrg_index( orb_j, orb_k, orb_m, orb_n, orb_i, orb_i, d45 + d46 );
878 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_n, orb_i, d45 - d46 );
879 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_i, orb_n, -2 * d45 );
880 }
881 }
882 }
883
884 delete tens_45;
885 delete tens_46;
886 if ( orb_m != orb_n ){ delete tens_47;
887 delete tens_48; }
888
889 }
890
891 /**************************************************************
892 * 3-1-2 : Find out how many contributions each process has *
893 **************************************************************/
894 int counter_jkl = 0;
895 for (int global = 0; global < upperbound4; global++){
896 Special::invert_triangle_three( global, jkl );
897 const int orb_j = jkl[ 0 ];
898 const int orb_k = jkl[ 1 ];
899 const int orb_l = jkl[ 2 ];
900 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
901 if ( irrep_jkl == irrep_imn ){
902 #ifdef CHEMPS2_MPI_COMPILATION
903 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK ){ counter_jkl++; }
904 #else
905 counter_jkl++;
906 #endif
907 }
908 }
909
910 /***********************************
911 * 3-1-2 : contributions A-Sigma *
912 ***********************************/
913 if ( counter_jkl > 0 ){
914
915 Tensor3RDM * a_S0_doublet = new Tensor3RDM( orb_i, -2, 1, 3, irrep_imn, true, book );
916 Tensor3RDM * a_S1_doublet = (( orb_m == orb_n ) ? NULL : new Tensor3RDM( orb_i, -2, 1, 3, irrep_imn, true, book ));
917 Tensor3RDM * a_S1_quartet = (( orb_m == orb_n ) ? NULL : new Tensor3RDM( orb_i, -2, 3, 3, irrep_imn, true, book ));
918 fill_a_S0( denT, a_S0_doublet, S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
919 if ( orb_m != orb_n ){ fill_a_S1( denT, a_S1_doublet, a_S1_quartet, S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 ); }
920
921 #ifdef CHEMPS2_MPI_COMPILATION
922 #pragma omp for schedule(dynamic) nowait
923 #else
924 #pragma omp for schedule(static) nowait
925 #endif
926 for (int global = 0; global < upperbound4; global++){
927 Special::invert_triangle_three( global, jkl );
928 const int orb_j = jkl[ 0 ];
929 const int orb_k = jkl[ 1 ];
930 const int orb_l = jkl[ 2 ];
931 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
932 if ( irrep_jkl == irrep_imn ){
933 #ifdef CHEMPS2_MPI_COMPILATION
934 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
935 #endif
936 {
937 const int cnt1 = orb_k - orb_j;
938 const int cnt2 = orb_l - orb_k;
939 const int cnt3 = orb_i - 1 - orb_l;
940 if ( cnt1 + cnt2 > 0 ){
941 const double d90_95 = a_S0_doublet->contract(dm3_a_J0_doublet[cnt1][cnt2][cnt3]);
942 const double d93_98 = ( (cnt1==0) ?0.0: sq3*a_S0_doublet->contract(dm3_a_J1_doublet[cnt1][cnt2][cnt3]));
943 const double d91_96 = (((orb_n==orb_m)||(cnt1==0))?0.0: a_S1_doublet->contract(dm3_a_J1_doublet[cnt1][cnt2][cnt3]));
944 const double d94_99 = ( (orb_n==orb_m) ?0.0: -sq3*a_S1_doublet->contract(dm3_a_J0_doublet[cnt1][cnt2][cnt3]));
945 const double d92_97 = (((orb_n==orb_m)||(cnt1==0))?0.0: a_S1_quartet->contract(dm3_a_J1_quartet[cnt1][cnt2][cnt3]));
946 set_dmrg_index( orb_j, orb_k, orb_l, orb_m, orb_n, orb_i, -2*d90_95 + 2*d91_96 + d92_97 );
947 set_dmrg_index( orb_j, orb_k, orb_l, orb_n, orb_m, orb_i, -2*d90_95 - 2*d91_96 - d92_97 );
948 set_dmrg_index( orb_j, orb_k, orb_l, orb_m, orb_i, orb_n, d90_95 + d91_96 - d92_97 + d93_98 + d94_99 );
949 set_dmrg_index( orb_j, orb_k, orb_l, orb_n, orb_i, orb_m, d90_95 - d91_96 + d92_97 + d93_98 - d94_99 );
950 set_dmrg_index( orb_j, orb_k, orb_l, orb_i, orb_m, orb_n, d90_95 - d91_96 + d92_97 - d93_98 + d94_99 );
951 set_dmrg_index( orb_j, orb_k, orb_l, orb_i, orb_n, orb_m, d90_95 + d91_96 - d92_97 - d93_98 - d94_99 );
952 }
953 }
954 }
955 }
956
957 delete a_S0_doublet;
958 if ( orb_m != orb_n ){ delete a_S1_doublet;
959 delete a_S1_quartet; }
960
961 }
962
963 /*************************************
964 * 3-1-2 : contributions BCD-Sigma *
965 *************************************/
966 if ( counter_jkl > 0 ){
967
968 Tensor3RDM * bcd_S0_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book );
969 Tensor3RDM * bcd_S1_doublet = (( orb_m == orb_n ) ? NULL : new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book ));
970 Tensor3RDM * bcd_S1_quartet = (( orb_m == orb_n ) ? NULL : new Tensor3RDM( orb_i, -2, 3, 1, irrep_imn, true, book ));
971 fill_bcd_S0( denT, bcd_S0_doublet, S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
972 if ( orb_m != orb_n ){ fill_bcd_S1( denT, bcd_S1_doublet, bcd_S1_quartet, S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 ); }
973
974 #ifdef CHEMPS2_MPI_COMPILATION
975 #pragma omp for schedule(dynamic) nowait
976 #else
977 #pragma omp for schedule(static) nowait
978 #endif
979 for (int global = 0; global < upperbound4; global++){
980 Special::invert_triangle_three( global, jkl );
981 const int orb_j = jkl[ 0 ];
982 const int orb_k = jkl[ 1 ];
983 const int orb_l = jkl[ 2 ];
984 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
985 if ( irrep_jkl == irrep_imn ){
986 #ifdef CHEMPS2_MPI_COMPILATION
987 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
988 #endif
989 {
990 const int cnt1 = orb_k - orb_j;
991 const int cnt2 = orb_l - orb_k;
992 const int cnt3 = orb_i - 1 - orb_l;
993 if ( cnt2 > 0 ){ // dm3_b if NOT k==l
994 const double d120_125 = bcd_S0_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]);
995 const double d124_129 = ( (cnt1==0) ?0.0: sq3*bcd_S0_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
996 const double d123_128 = ( (orb_n==orb_m) ?0.0: sq3*bcd_S1_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]));
997 const double d121_126 = (((orb_n==orb_m)||(cnt1==0))?0.0: bcd_S1_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
998 const double d122_127 = (((orb_n==orb_m)||(cnt1==0))?0.0: bcd_S1_quartet->contract(dm3_b_J1_quartet[cnt1][cnt2][cnt3]));
999 set_dmrg_index( orb_l, orb_m, orb_n, orb_j, orb_k, orb_i, -d120_125 + 3*d121_126 + d123_128 - d124_129 );
1000 set_dmrg_index( orb_l, orb_m, orb_n, orb_k, orb_j, orb_i, -d120_125 - 3*d121_126 + d123_128 + d124_129 );
1001 set_dmrg_index( orb_l, orb_m, orb_n, orb_j, orb_i, orb_k, -d120_125 - 3*d121_126 - d123_128 - d124_129 );
1002 set_dmrg_index( orb_l, orb_m, orb_n, orb_k, orb_i, orb_j, -d120_125 + 3*d121_126 - d123_128 + d124_129 );
1003 set_dmrg_index( orb_l, orb_m, orb_n, orb_i, orb_j, orb_k, 2*d120_125 + 2*d121_126 + 2*d122_127 );
1004 set_dmrg_index( orb_l, orb_m, orb_n, orb_i, orb_k, orb_j, 2*d120_125 - 2*d121_126 - 2*d122_127 );
1005 }
1006 if ( cnt1 + cnt2 > 0 ){ // dm3_c if NOT j==k==l
1007 const double d110_115 = bcd_S0_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]);
1008 const double d112_117 = -sq3*bcd_S0_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]);
1009 const double d111_116 = ((orb_n==orb_m)?0.0: -sq3*bcd_S1_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]));
1010 const double d113_118 = ((orb_n==orb_m)?0.0: bcd_S1_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]));
1011 const double d114_119 = ((orb_n==orb_m)?0.0: -2*bcd_S1_quartet->contract(dm3_c_J1_quartet[cnt1][cnt2][cnt3]));
1012 set_dmrg_index( orb_k, orb_m, orb_n, orb_j, orb_l, orb_i, 2*d110_115 + 2*d111_116 );
1013 set_dmrg_index( orb_k, orb_m, orb_n, orb_l, orb_j, orb_i, -d110_115 - d111_116 - d112_117 - 3*d113_118 );
1014 set_dmrg_index( orb_k, orb_m, orb_n, orb_j, orb_i, orb_l, 2*d110_115 - 2*d111_116 );
1015 set_dmrg_index( orb_k, orb_m, orb_n, orb_l, orb_i, orb_j, -d110_115 + d111_116 - d112_117 + 3*d113_118 );
1016 set_dmrg_index( orb_k, orb_m, orb_n, orb_i, orb_j, orb_l, -d110_115 + d111_116 + d112_117 + d113_118 + d114_119 );
1017 set_dmrg_index( orb_k, orb_m, orb_n, orb_i, orb_l, orb_j, -d110_115 - d111_116 + d112_117 - d113_118 - d114_119 );
1018 }
1019 // dm3_d for all j <= k <= l
1020 const double d100_105 = -bcd_S0_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]);
1021 const double d102_107 = -sq3*bcd_S0_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]);
1022 const double d101_106 = ( (orb_n==orb_m) ?0.0: sq3*bcd_S1_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]));
1023 const double d103_108 = ( (orb_n==orb_m) ?0.0: bcd_S1_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]));
1024 const double d104_109 = (((orb_n==orb_m)||(cnt2==0))?0.0: 2*bcd_S1_quartet->contract(dm3_d_J1_quartet[cnt1][cnt2][cnt3]));
1025 set_dmrg_index( orb_j, orb_m, orb_n, orb_k, orb_l, orb_i, 2*d100_105 + 2*d101_106 );
1026 set_dmrg_index( orb_j, orb_m, orb_n, orb_l, orb_k, orb_i, -d100_105 - d101_106 - d102_107 - 3*d103_108 );
1027 set_dmrg_index( orb_j, orb_m, orb_n, orb_k, orb_i, orb_l, 2*d100_105 - 2*d101_106 );
1028 set_dmrg_index( orb_j, orb_m, orb_n, orb_l, orb_i, orb_k, -d100_105 + d101_106 - d102_107 + 3*d103_108 );
1029 set_dmrg_index( orb_j, orb_m, orb_n, orb_i, orb_k, orb_l, -d100_105 + d101_106 + d102_107 + d103_108 + d104_109 );
1030 set_dmrg_index( orb_j, orb_m, orb_n, orb_i, orb_l, orb_k, -d100_105 - d101_106 + d102_107 - d103_108 - d104_109 );
1031 }
1032 }
1033 }
1034
1035 delete bcd_S0_doublet;
1036 if ( orb_m != orb_n ){ delete bcd_S1_doublet;
1037 delete bcd_S1_quartet; }
1038
1039 }
1040
1041 /***************************************
1042 * 3-1-2 : contributions F transpose *
1043 ***************************************/
1044 if ( counter_jkl > 0 ){
1045
1046 Tensor3RDM * F0_T_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book );
1047 Tensor3RDM * F1_T_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book );
1048 Tensor3RDM * F1_T_quartet = new Tensor3RDM( orb_i, -2, 3, 1, irrep_imn, true, book );
1049 fill_F0_T( denT, F0_T_doublet, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
1050 fill_F1_T( denT, F1_T_doublet, F1_T_quartet, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
1051
1052 #ifdef CHEMPS2_MPI_COMPILATION
1053 #pragma omp for schedule(dynamic) nowait
1054 #else
1055 #pragma omp for schedule(static) nowait
1056 #endif
1057 for (int global = 0; global < upperbound4; global++){
1058 Special::invert_triangle_three( global, jkl );
1059 const int orb_j = jkl[ 0 ];
1060 const int orb_k = jkl[ 1 ];
1061 const int orb_l = jkl[ 2 ];
1062 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
1063 if ( irrep_jkl == irrep_imn ){
1064 #ifdef CHEMPS2_MPI_COMPILATION
1065 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
1066 #endif
1067 {
1068 const int cnt1 = orb_k - orb_j;
1069 const int cnt2 = orb_l - orb_k;
1070 const int cnt3 = orb_i - 1 - orb_l;
1071 if ( cnt2 > 0 ){ // dm3_b if NOT k==l
1072 const double d130_135 = F0_T_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]);
1073 const double d131_136 = ((cnt1==0)?0.0: sq3*F0_T_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
1074 const double d132_137 = sq3*F1_T_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]);
1075 const double d133_138 = ((cnt1==0)?0.0: F1_T_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
1076 const double d134_139 = ((cnt1==0)?0.0: F1_T_quartet->contract(dm3_b_J1_quartet[cnt1][cnt2][cnt3]));
1077 set_dmrg_index( orb_j, orb_k, orb_m, orb_l, orb_n, orb_i, d130_135 + d131_136 + d132_137 + 3*d133_138 );
1078 set_dmrg_index( orb_j, orb_k, orb_m, orb_n, orb_l, orb_i, d130_135 - d131_136 + d132_137 - 3*d133_138 );
1079 set_dmrg_index( orb_j, orb_k, orb_m, orb_l, orb_i, orb_n, -2*d130_135 - 2*d131_136 );
1080 set_dmrg_index( orb_j, orb_k, orb_m, orb_n, orb_i, orb_l, d130_135 + d131_136 - d132_137 + d133_138 + 2*d134_139 );
1081 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_l, orb_n, -2*d130_135 + 2*d131_136 );
1082 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_n, orb_l, d130_135 - d131_136 - d132_137 - d133_138 - 2*d134_139 );
1083 }
1084 if ( cnt1 + cnt2 > 0 ){ // dm3_c if NOT j==k==l
1085 const double d150_155 = F0_T_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]);
1086 const double d151_156 = -sq3*F0_T_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]);
1087 const double d152_157 = sq3*F1_T_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]);
1088 const double d153_158 = F1_T_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]);
1089 const double d154_159 = -F1_T_quartet->contract(dm3_c_J1_quartet[cnt1][cnt2][cnt3]);
1090 set_dmrg_index( orb_j, orb_l, orb_m, orb_k, orb_n, orb_i, -2*d150_155 - 2*d152_157 );
1091 set_dmrg_index( orb_j, orb_l, orb_m, orb_n, orb_k, orb_i, d150_155 + d151_156 + d152_157 - 3*d153_158 );
1092 set_dmrg_index( orb_j, orb_l, orb_m, orb_k, orb_i, orb_n, 4*d150_155 );
1093 set_dmrg_index( orb_j, orb_l, orb_m, orb_n, orb_i, orb_k, -2*d150_155 + 2*d153_158 + 2*d154_159 );
1094 set_dmrg_index( orb_j, orb_l, orb_m, orb_i, orb_k, orb_n, -2*d150_155 - 2*d151_156 );
1095 set_dmrg_index( orb_j, orb_l, orb_m, orb_i, orb_n, orb_k, d150_155 + d151_156 + d152_157 + d153_158 - 2*d154_159 );
1096 }
1097 // dm3_d for all j <= k <= l
1098 const double d170_175 = -F0_T_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]);
1099 const double d171_176 = -sq3*F0_T_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]);
1100 const double d172_177 = -sq3*F1_T_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]);
1101 const double d173_178 = F1_T_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]);
1102 const double d174_179 = ((cnt2==0)?0.0: F1_T_quartet->contract(dm3_d_J1_quartet[cnt1][cnt2][cnt3]));
1103 set_dmrg_index( orb_k, orb_l, orb_m, orb_j, orb_n, orb_i, -2*d170_175 - 2*d172_177 );
1104 set_dmrg_index( orb_k, orb_l, orb_m, orb_n, orb_j, orb_i, d170_175 + d171_176 + d172_177 - 3*d173_178 );
1105 set_dmrg_index( orb_k, orb_l, orb_m, orb_j, orb_i, orb_n, 4*d170_175 );
1106 set_dmrg_index( orb_k, orb_l, orb_m, orb_n, orb_i, orb_j, -2*d170_175 + 2*d173_178 + 2*d174_179 );
1107 set_dmrg_index( orb_k, orb_l, orb_m, orb_i, orb_j, orb_n, -2*d170_175 - 2*d171_176 );
1108 set_dmrg_index( orb_k, orb_l, orb_m, orb_i, orb_n, orb_j, d170_175 + d171_176 + d172_177 + d173_178 - 2*d174_179 );
1109 }
1110 }
1111 }
1112
1113 delete F0_T_doublet;
1114 delete F1_T_doublet;
1115 delete F1_T_quartet;
1116
1117 }
1118
1119 /*************************************
1120 * 3-1-2 : contributions F regular *
1121 *************************************/
1122 if (( orb_m != orb_n ) && ( counter_jkl > 0 )){
1123
1124 Tensor3RDM * F0_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book );
1125 Tensor3RDM * F1_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_imn, true, book );
1126 Tensor3RDM * F1_quartet = new Tensor3RDM( orb_i, -2, 3, 1, irrep_imn, true, book );
1127 fill_F0( denT, F0_doublet, F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem );
1128 fill_F1( denT, F1_doublet, F1_quartet, F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i], workmem, workmem2 );
1129
1130 #ifdef CHEMPS2_MPI_COMPILATION
1131 #pragma omp for schedule(dynamic) nowait
1132 #else
1133 #pragma omp for schedule(static) nowait
1134 #endif
1135 for (int global = 0; global < upperbound4; global++){
1136 Special::invert_triangle_three( global, jkl );
1137 const int orb_j = jkl[ 0 ];
1138 const int orb_k = jkl[ 1 ];
1139 const int orb_l = jkl[ 2 ];
1140 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
1141 if ( irrep_jkl == irrep_imn ){
1142 #ifdef CHEMPS2_MPI_COMPILATION
1143 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
1144 #endif
1145 {
1146 const int cnt1 = orb_k - orb_j;
1147 const int cnt2 = orb_l - orb_k;
1148 const int cnt3 = orb_i - 1 - orb_l;
1149 if ( cnt2 > 0 ){ // dm3_b if NOT k==l
1150 const double d140_145 = F0_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]);
1151 const double d141_146 = ((cnt1==0)?0.0: sq3*F0_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
1152 const double d142_147 = sq3*F1_doublet->contract(dm3_b_J0_doublet[cnt1][cnt2][cnt3]);
1153 const double d143_148 = ((cnt1==0)?0.0: F1_doublet->contract(dm3_b_J1_doublet[cnt1][cnt2][cnt3]));
1154 const double d144_149 = ((cnt1==0)?0.0: F1_quartet->contract(dm3_b_J1_quartet[cnt1][cnt2][cnt3]));
1155 set_dmrg_index( orb_j, orb_k, orb_n, orb_l, orb_m, orb_i, d140_145 + d141_146 + d142_147 + 3*d143_148 );
1156 set_dmrg_index( orb_j, orb_k, orb_n, orb_m, orb_l, orb_i, d140_145 - d141_146 + d142_147 - 3*d143_148 );
1157 set_dmrg_index( orb_j, orb_k, orb_n, orb_l, orb_i, orb_m, -2*d140_145 - 2*d141_146 );
1158 set_dmrg_index( orb_j, orb_k, orb_n, orb_m, orb_i, orb_l, d140_145 + d141_146 - d142_147 + d143_148 + 2*d144_149 );
1159 set_dmrg_index( orb_j, orb_k, orb_n, orb_i, orb_l, orb_m, -2*d140_145 + 2*d141_146 );
1160 set_dmrg_index( orb_j, orb_k, orb_n, orb_i, orb_m, orb_l, d140_145 - d141_146 - d142_147 - d143_148 - 2*d144_149 );
1161 }
1162 if ( cnt1 + cnt2 > 0 ){ // dm3_c if NOT j==k==l
1163 const double d160_165 = F0_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]);
1164 const double d161_166 = -sq3*F0_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]);
1165 const double d162_167 = sq3*F1_doublet->contract(dm3_c_J0_doublet[cnt1][cnt2][cnt3]);
1166 const double d163_168 = F1_doublet->contract(dm3_c_J1_doublet[cnt1][cnt2][cnt3]);
1167 const double d164_169 = -F1_quartet->contract(dm3_c_J1_quartet[cnt1][cnt2][cnt3]);
1168 set_dmrg_index( orb_j, orb_l, orb_n, orb_k, orb_m, orb_i, -2*d160_165 - 2*d162_167 );
1169 set_dmrg_index( orb_j, orb_l, orb_n, orb_m, orb_k, orb_i, d160_165 + d161_166 + d162_167 - 3*d163_168 );
1170 set_dmrg_index( orb_j, orb_l, orb_n, orb_k, orb_i, orb_m, 4*d160_165 );
1171 set_dmrg_index( orb_j, orb_l, orb_n, orb_m, orb_i, orb_k, -2*d160_165 + 2*d163_168 + 2*d164_169 );
1172 set_dmrg_index( orb_j, orb_l, orb_n, orb_i, orb_k, orb_m, -2*d160_165 - 2*d161_166 );
1173 set_dmrg_index( orb_j, orb_l, orb_n, orb_i, orb_m, orb_k, d160_165 + d161_166 + d162_167 + d163_168 - 2*d164_169 );
1174 }
1175 // dm3_d for all j <= k <= l
1176 const double d180_185 = -F0_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]);
1177 const double d181_186 = -sq3*F0_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]);
1178 const double d182_187 = -sq3*F1_doublet->contract(dm3_d_J0_doublet[cnt1][cnt2][cnt3]);
1179 const double d183_188 = F1_doublet->contract(dm3_d_J1_doublet[cnt1][cnt2][cnt3]);
1180 const double d184_189 = ((cnt2==0)?0.0: F1_quartet->contract(dm3_d_J1_quartet[cnt1][cnt2][cnt3]));
1181 set_dmrg_index( orb_k, orb_l, orb_n, orb_j, orb_m, orb_i, -2*d180_185 - 2*d182_187 );
1182 set_dmrg_index( orb_k, orb_l, orb_n, orb_m, orb_j, orb_i, d180_185 + d181_186 + d182_187 - 3*d183_188 );
1183 set_dmrg_index( orb_k, orb_l, orb_n, orb_j, orb_i, orb_m, 4*d180_185 );
1184 set_dmrg_index( orb_k, orb_l, orb_n, orb_m, orb_i, orb_j, -2*d180_185 + 2*d183_188 + 2*d184_189 );
1185 set_dmrg_index( orb_k, orb_l, orb_n, orb_i, orb_j, orb_m, -2*d180_185 - 2*d181_186 );
1186 set_dmrg_index( orb_k, orb_l, orb_n, orb_i, orb_m, orb_j, d180_185 + d181_186 + d182_187 + d183_188 - 2*d184_189 );
1187 }
1188 }
1189 }
1190
1191 delete F0_doublet;
1192 delete F1_doublet;
1193 delete F1_quartet;
1194
1195 }
1196
1197 /***************************************
1198 * Clean up the double S_mn and F_mn *
1199 ***************************************/
1200 #ifdef CHEMPS2_MPI_COMPILATION
1201 #pragma omp barrier // Everyone needs to be done before tensors are deleted
1202 #pragma omp single
1203 if ( orb_m > orb_i + 1 ){ // All processes own Fx/Sx[ index ][ n - m ][ m - i - 1 == 0 ]
1204 const int own_S_mn = MPIchemps2::owner_absigma( orb_m, orb_n );
1205 const int own_F_mn = MPIchemps2::owner_cdf( L, orb_m, orb_n );
1206 if ( MPIRANK != own_F_mn ){ delete F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i]; F0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = NULL;
1207 delete F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i]; F1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = NULL; }
1208 if ( MPIRANK != own_S_mn ){ delete S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i]; S0tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = NULL;
1209 if ( orb_m != orb_n ){ delete S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i]; S1tensors[orb_i][orb_n-orb_m][orb_m-1-orb_i] = NULL; } }
1210 }
1211 #endif
1212
1213 }
1214 }
1215
1216 for ( int orb_m = orb_i+1; orb_m < L; orb_m++ ){
1217
1218 const int irrep_m = book->gIrrep( orb_m );
1219
1220 /**************************************************************
1221 * 3-2-1 : Find out how many contributions each process has *
1222 **************************************************************/
1223 int counter_jkl = 0;
1224 for (int global = 0; global < upperbound4; global++){
1225 Special::invert_triangle_three( global, jkl );
1226 const int orb_j = jkl[ 0 ];
1227 const int orb_k = jkl[ 1 ];
1228 const int orb_l = jkl[ 2 ];
1229 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
1230 if ( irrep_jkl == irrep_m ){
1231 #ifdef CHEMPS2_MPI_COMPILATION
1232 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK ){ counter_jkl++; }
1233 #else
1234 counter_jkl++;
1235 #endif
1236 }
1237 }
1238
1239 /********************
1240 * 3-2-1 : part 1 *
1241 ********************/
1242 if ( counter_jkl > 0 ){
1243
1244 Tensor3RDM * A_T_doublet = new Tensor3RDM( orb_i, -2, 1, 3, irrep_m, true, book );
1245 Tensor3RDM * BCD_T_doublet = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1246 fill_53_54( denT, A_T_doublet, Ltensors[orb_i][orb_m-1-orb_i], workmem );
1247 fill_55_to_60( denT, BCD_T_doublet, Ltensors[orb_i][orb_m-1-orb_i], workmem );
1248
1249 #ifdef CHEMPS2_MPI_COMPILATION
1250 #pragma omp for schedule(dynamic) nowait
1251 #else
1252 #pragma omp for schedule(static) nowait
1253 #endif
1254 for ( int global = 0; global < upperbound4; global++ ){
1255 Special::invert_triangle_three( global, jkl );
1256 const int orb_j = jkl[ 0 ];
1257 const int orb_k = jkl[ 1 ];
1258 const int orb_l = jkl[ 2 ];
1259 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
1260 if ( irrep_jkl == irrep_m ){
1261 #ifdef CHEMPS2_MPI_COMPILATION
1262 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
1263 #endif
1264 {
1265 const int cnt1 = orb_k - orb_j;
1266 const int cnt2 = orb_l - orb_k;
1267 const int cnt3 = orb_i - 1 - orb_l;
1268 if ( cnt1 + cnt2 > 0 ){
1269 const double d53 = A_T_doublet->contract( dm3_a_J0_doublet[cnt1][cnt2][cnt3] );
1270 const double d54 = sq3*A_T_doublet->contract( dm3_a_J1_doublet[cnt1][cnt2][cnt3] );
1271 set_dmrg_index( orb_j, orb_k, orb_l, orb_m, orb_i, orb_i, d53 - d54 );
1272 set_dmrg_index( orb_j, orb_k, orb_l, orb_i, orb_m, orb_i, d53 + d54 );
1273 set_dmrg_index( orb_j, orb_k, orb_l, orb_i, orb_i, orb_m, - 2 * d53 );
1274 }
1275 if ( cnt2 > 0 ){
1276 const double d55 = BCD_T_doublet->contract( dm3_b_J0_doublet[cnt1][cnt2][cnt3] );
1277 const double d56 = sq3*BCD_T_doublet->contract( dm3_b_J1_doublet[cnt1][cnt2][cnt3] );
1278 set_dmrg_index( orb_j, orb_k, orb_m, orb_l, orb_i, orb_i, d55 + d56 );
1279 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_l, orb_i, d55 - d56 );
1280 set_dmrg_index( orb_j, orb_k, orb_m, orb_i, orb_i, orb_l, - 2 * d55 );
1281 }
1282 if ( cnt1 + cnt2 > 0 ){
1283 const double d57 = BCD_T_doublet->contract( dm3_c_J0_doublet[cnt1][cnt2][cnt3] );
1284 const double d58 = sq3*BCD_T_doublet->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1285 set_dmrg_index( orb_j, orb_l, orb_m, orb_k, orb_i, orb_i, - 2 * d57 );
1286 set_dmrg_index( orb_j, orb_l, orb_m, orb_i, orb_k, orb_i, d57 - d58 );
1287 set_dmrg_index( orb_j, orb_l, orb_m, orb_i, orb_i, orb_k, d57 + d58 );
1288 }
1289 const double d59 = -BCD_T_doublet->contract( dm3_d_J0_doublet[cnt1][cnt2][cnt3] );
1290 const double d60 = -sq3*BCD_T_doublet->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1291 set_dmrg_index( orb_k, orb_l, orb_m, orb_j, orb_i, orb_i, - 2 * d59 );
1292 set_dmrg_index( orb_k, orb_l, orb_m, orb_i, orb_j, orb_i, d59 + d60 );
1293 set_dmrg_index( orb_k, orb_l, orb_m, orb_i, orb_i, orb_j, d59 - d60 );
1294 }
1295 }
1296 }
1297
1298 delete A_T_doublet;
1299 delete BCD_T_doublet;
1300
1301 }
1302
1303 /**************************************************************
1304 * 1-4-1 : Find out how many contributions each process has *
1305 **************************************************************/
1306 int counter_j = 0;
1307 for ( int orb_j = 0; orb_j < orb_i; orb_j++ ){
1308 const int irrep_j = book->gIrrep( orb_j );
1309 if ( irrep_j == irrep_m ){
1310 #ifdef CHEMPS2_MPI_COMPILATION
1311 if ( MPIRANK == MPIchemps2::owner_q( L, orb_j ) ){ counter_j++; }
1312 #else
1313 counter_j++;
1314 #endif
1315 }
1316 }
1317
1318 /***********************
1319 * 3-2-1 : part 2 *
1320 * 1-4-1 : L^T - L^T *
1321 ***********************/
1322 if (( counter_jkl > 0 ) || ( counter_j > 0 )){
1323
1324 Tensor3RDM * tens_61 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1325 Tensor3RDM * tens_63 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1326 Tensor3RDM * tens_65 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1327 Tensor3RDM * tens_67 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1328 Tensor3RDM * tens_68 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1329 Tensor3RDM * tens_76 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1330 Tensor3RDM * tens_77 = new Tensor3RDM( orb_i, -2, 1, 1, irrep_m, true, book );
1331 Tensor3RDM * tens_69 = new Tensor3RDM( orb_i, -2, 3, 1, irrep_m, true, book );
1332 Tensor3RDM * tens_78 = new Tensor3RDM( orb_i, -2, 3, 1, irrep_m, true, book );
1333 Tensor3RDM * tens_79 = new Tensor3RDM( orb_i, -2, 3, 1, irrep_m, true, book );
1334 fill_61( denT, tens_61, Ltensors[orb_i][orb_m-1-orb_i], workmem );
1335 fill_63_65( denT, tens_63, tens_65, tens_67, tens_68, tens_76, tens_77, Ltensors[orb_i][orb_m-1-orb_i], workmem, workmem2 );
1336 fill_69_78_79( denT, tens_69, tens_78, tens_79, Ltensors[orb_i][orb_m-1-orb_i], workmem, workmem2 );
1337
1338 #ifdef CHEMPS2_MPI_COMPILATION
1339 #pragma omp for schedule(dynamic) nowait
1340 #else
1341 #pragma omp for schedule(static) nowait
1342 #endif
1343 for ( int orb_j = 0; orb_j < orb_i; orb_j++ ){
1344 const int irrep_j = book->gIrrep( orb_j );
1345 if ( irrep_j == irrep_m ){
1346 #ifdef CHEMPS2_MPI_COMPILATION
1347 if ( MPIRANK == MPIchemps2::owner_q( L, orb_j ) ) //Everyone owns the L-tensors --> task division based on Q-tensor ownership
1348 #endif
1349 {
1350 const double d2 = sqrt( 2.0 ) * tens_61->inproduct( Ltensors[orb_i-1][orb_i-1-orb_j], 'N' );
1351 set_dmrg_index( orb_j, orb_i, orb_i, orb_m, orb_i, orb_i, 2 * d2 );
1352 set_dmrg_index( orb_j, orb_i, orb_i, orb_i, orb_m, orb_i, - d2 );
1353 }
1354 }
1355 }
1356
1357 #ifdef CHEMPS2_MPI_COMPILATION
1358 #pragma omp for schedule(dynamic) nowait
1359 #else
1360 #pragma omp for schedule(static) nowait
1361 #endif
1362 for ( int global = 0; global < upperbound4; global++ ){
1363 Special::invert_triangle_three( global, jkl );
1364 const int orb_j = jkl[ 0 ];
1365 const int orb_k = jkl[ 1 ];
1366 const int orb_l = jkl[ 2 ];
1367 const int irrep_jkl = Irreps::directProd( Irreps::directProd( book->gIrrep( orb_j ), book->gIrrep( orb_k ) ), book->gIrrep( orb_l ) );
1368 if ( irrep_jkl == irrep_m ){
1369 #ifdef CHEMPS2_MPI_COMPILATION
1370 if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
1371 #endif
1372 {
1373 const int cnt1 = orb_k - orb_j;
1374 const int cnt2 = orb_l - orb_k;
1375 const int cnt3 = orb_i - 1 - orb_l;
1376 if ( cnt2 > 0 ){
1377 const double d61 = tens_61->contract( dm3_b_J0_doublet[cnt1][cnt2][cnt3] );
1378 const double d62 = sq3*tens_61->contract( dm3_b_J1_doublet[cnt1][cnt2][cnt3] );
1379 const double d63 = tens_63->contract( dm3_b_J0_doublet[cnt1][cnt2][cnt3] );
1380 const double d64 = sq3*tens_63->contract( dm3_b_J1_doublet[cnt1][cnt2][cnt3] );
1381 const double d65 = tens_65->contract( dm3_b_J0_doublet[cnt1][cnt2][cnt3] );
1382 const double d66 = sq3*tens_65->contract( dm3_b_J1_doublet[cnt1][cnt2][cnt3] );
1383 const double d67 = tens_67->contract( dm3_b_J0_doublet[cnt1][cnt2][cnt3] );
1384 const double d68 = tens_68->contract( dm3_b_J1_doublet[cnt1][cnt2][cnt3] );
1385 const double d69 = tens_69->contract( dm3_b_J1_quartet[cnt1][cnt2][cnt3] );
1386 set_dmrg_index( orb_j, orb_k, orb_i, orb_l, orb_m, orb_i, -2*d61 - 2*d62 + d63 + d64 );
1387 set_dmrg_index( orb_j, orb_k, orb_i, orb_m, orb_l, orb_i, -2*d61 + 2*d62 + d63 - d64 );
1388 set_dmrg_index( orb_j, orb_k, orb_i, orb_l, orb_i, orb_m, d61 + d62 + d65 + d66 );
1389 set_dmrg_index( orb_j, orb_k, orb_i, orb_m, orb_i, orb_l, d61 - d62 + d67 + d68 + d69 );
1390 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_m, orb_l, d61 + d62 + d67 - d68 - d69 );
1391 set_dmrg_index( orb_j, orb_k, orb_i, orb_i, orb_l, orb_m, d61 - d62 + d65 - d66 );
1392 }
1393 if ( cnt1 + cnt2 > 0 ){
1394 const double d70 = tens_61->contract( dm3_c_J0_doublet[cnt1][cnt2][cnt3] );
1395 const double d71 = sq3*tens_61->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1396 const double d72 = -tens_63->contract( dm3_c_J0_doublet[cnt1][cnt2][cnt3] );
1397 const double d73 = -sq3*tens_63->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1398 const double d74 = -tens_65->contract( dm3_c_J0_doublet[cnt1][cnt2][cnt3] );
1399 const double d75 = -sq3*tens_65->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1400 const double d76 = tens_76->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1401 const double d77 = tens_77->contract( dm3_c_J1_doublet[cnt1][cnt2][cnt3] );
1402 const double d78 = tens_78->contract( dm3_c_J1_quartet[cnt1][cnt2][cnt3] );
1403 const double d79 = tens_79->contract( dm3_c_J1_quartet[cnt1][cnt2][cnt3] );
1404 set_dmrg_index( orb_j, orb_l, orb_i, orb_k, orb_m, orb_i, 4*d70 + 2*d72 );
1405 set_dmrg_index( orb_j, orb_l, orb_i, orb_m, orb_k, orb_i, -2*d70 + 2*d71 - d72 + d73 );
1406 set_dmrg_index( orb_j, orb_l, orb_i, orb_k, orb_i, orb_m, -2*d70 + 2*d74 );
1407 set_dmrg_index( orb_j, orb_l, orb_i, orb_m, orb_i, orb_k, d70 - d71 - d74 + d76 + d78 );
1408 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_k, orb_m, d70 - d71 - d74 + d75 );
1409 set_dmrg_index( orb_j, orb_l, orb_i, orb_i, orb_m, orb_k, -2*d70 - d72 + d77 + d79 );
1410 }
1411 const double d80 = -tens_61->contract( dm3_d_J0_doublet[cnt1][cnt2][cnt3] );
1412 const double d81 = -sq3*tens_61->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1413 const double d82 = tens_63->contract( dm3_d_J0_doublet[cnt1][cnt2][cnt3] );
1414 const double d83 = sq3*tens_63->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1415 const double d84 = tens_65->contract( dm3_d_J0_doublet[cnt1][cnt2][cnt3] );
1416 const double d85 = sq3*tens_65->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1417 const double d86 = tens_76->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1418 const double d87 = tens_77->contract( dm3_d_J1_doublet[cnt1][cnt2][cnt3] );
1419 const double d88 = -tens_78->contract( dm3_d_J1_quartet[cnt1][cnt2][cnt3] );
1420 const double d89 = -tens_79->contract( dm3_d_J1_quartet[cnt1][cnt2][cnt3] );
1421 set_dmrg_index( orb_j, orb_m, orb_i, orb_k, orb_l, orb_i, 4*d80 + 2*d82 );
1422 set_dmrg_index( orb_j, orb_m, orb_i, orb_l, orb_k, orb_i, -2*d80 - 2*d81 - d82 - d83 );
1423 set_dmrg_index( orb_j, orb_m, orb_i, orb_k, orb_i, orb_l, -2*d80 + 2*d84 );
1424 set_dmrg_index( orb_j, orb_m, orb_i, orb_l, orb_i, orb_k, d80 + d81 - d84 - d85 );
1425 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_k, orb_l, d80 + d81 - d84 + d86 + d88 );
1426 set_dmrg_index( orb_j, orb_m, orb_i, orb_i, orb_l, orb_k, -2*d80 - d82 + d87 + d89 );
1427 }
1428 }
1429 }
1430
1431 delete tens_61;
1432 delete tens_63;
1433 delete tens_65;
1434 delete tens_67;
1435 delete tens_68;
1436 delete tens_69;
1437 delete tens_76;
1438 delete tens_77;
1439 delete tens_78;
1440 delete tens_79;
1441 }
1442
1443 }
1444
1445 delete [] workmem;
1446 delete [] workmem2;
1447
1448 }
1449
1450 if (( disk ) && ( temp_disk_counter > 0 )){ flush_disk(); }
1451
1452 }
1453
diagram1(TensorT * denT,TensorF0 * denF0,double * workmem) const1454 double CheMPS2::ThreeDM::diagram1(TensorT * denT, TensorF0 * denF0, double * workmem) const{
1455
1456 assert( denF0->get_irrep() == 0 );
1457 const int orb_i = denT->gIndex();
1458
1459 double total = 0.0;
1460
1461 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1462 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1463 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1464
1465 int dimL = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1466 int dimR = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
1467
1468 if (( dimL > 0 ) && ( dimR > 0 )){
1469
1470 double * Tblock = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
1471 double * Fblock = denF0->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
1472
1473 char notrans = 'N';
1474 double alpha = 1.0;
1475 double beta = 0.0; //set
1476 dgemm_( ¬rans, ¬rans, &dimL, &dimR, &dimL, &alpha, Fblock, &dimL, Tblock, &dimL, &beta, workmem, &dimL );
1477
1478 int length = dimL * dimR;
1479 int inc = 1;
1480 total += ( TwoSL + 1 ) * ddot_( &length, workmem, &inc, Tblock, &inc );
1481
1482 }
1483 }
1484 }
1485 }
1486
1487 total *= sqrt( 0.5 );
1488 return total;
1489
1490 }
1491
diagram3(TensorT * denT,TensorF0 * denF0,double * workmem) const1492 double CheMPS2::ThreeDM::diagram3(TensorT * denT, TensorF0 * denF0, double * workmem) const{
1493
1494 assert( denF0->get_irrep() == 0 );
1495 const int orb_i = denT->gIndex();
1496
1497 double total = 0.0;
1498
1499 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1500 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1501 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1502
1503 int dimL = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1504 int dimR = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
1505
1506 if (( dimL > 0 ) && ( dimR > 0 )){
1507
1508 double * Tblock = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
1509 double * Fblock = denF0->gStorage( NL+2, TwoSL, IL, NL+2, TwoSL, IL );
1510
1511 char notrans = 'N';
1512 double alpha = 1.0;
1513 double beta = 0.0; //set
1514 dgemm_( ¬rans, ¬rans, &dimL, &dimR, &dimR, &alpha, Tblock, &dimL, Fblock, &dimR, &beta, workmem, &dimL );
1515
1516 int length = dimL * dimR;
1517 int inc = 1;
1518 total += ( TwoSL + 1 ) * ddot_( &length, workmem, &inc, Tblock, &inc );
1519
1520 }
1521 }
1522 }
1523 }
1524
1525 total *= sqrt( 0.5 );
1526 return total;
1527
1528 }
1529
diagram4_5_6_7_8_9(TensorT * denT,Tensor3RDM * d3tens,double * workmem,const char type) const1530 double CheMPS2::ThreeDM::diagram4_5_6_7_8_9(TensorT * denT, Tensor3RDM * d3tens, double * workmem, const char type) const{
1531
1532 const int orb_i = denT->gIndex();
1533 assert( d3tens->get_irrep() == book->gIrrep( orb_i ) );
1534 assert( d3tens->get_two_j2() == 1 );
1535
1536 double total = 0.0;
1537
1538 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1539 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1540 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1541
1542 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1543
1544 if ( dimLup > 0 ){
1545
1546 const int ILdown = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1547
1548 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
1549
1550 int dimR = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILdown );
1551 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILdown );
1552
1553 if (( dimLdown > 0 ) && ( dimR > 0 )){
1554
1555 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILdown );
1556 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILdown, NL+1, TwoSLprime, ILdown );
1557 double * block = d3tens->gStorage( NL-1, TwoSLprime, ILdown, NL, TwoSL, IL );
1558
1559 char notrans = 'N';
1560 double alpha = 1.0;
1561 double beta = 0.0; //set
1562 dgemm_( ¬rans, ¬rans, &dimLdown, &dimR, &dimLup, &alpha, block, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
1563
1564 int length = dimLdown * dimR;
1565 int inc = 1;
1566 const double factor = ((type =='D') ? ( sqrt( 0.5 * ( d3tens->get_two_j1() + 1 ) ) * ( TwoSLprime + 1 ) )
1567 : ( Special::phase( TwoSL + 1 - TwoSLprime )
1568 * sqrt( 0.5 * ( d3tens->get_two_j1() + 1 ) * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) ));
1569 total += factor * ddot_( &length, workmem, &inc, Tdown, &inc );
1570
1571 }
1572 }
1573 }
1574 }
1575 }
1576 }
1577
1578 return total;
1579
1580 }
1581
diagram10(TensorT * denT,TensorS0 * denS0,TensorL * denL,double * workmem,double * workmem2) const1582 double CheMPS2::ThreeDM::diagram10(TensorT * denT, TensorS0 * denS0, TensorL * denL, double * workmem, double * workmem2) const{
1583
1584 const int orb_i = denT->gIndex();
1585 assert( denS0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1586
1587 double total = 0.0;
1588
1589 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1590 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1591 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1592
1593 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1594 const int ILxIixIl = Irreps::directProd( IL, denS0->get_irrep() );
1595
1596 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1597 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSL, ILxIixIl );
1598 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxIixIl );
1599
1600 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
1601
1602 double * Tdown = denT->gStorage( NL-2, TwoSL, ILxIixIl, NL, TwoSL, ILxIixIl );
1603 double * Sblock = denS0->gStorage( NL-2, TwoSL, ILxIixIl, NL, TwoSL, IL );
1604
1605 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1606
1607 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1608
1609 if ( dimRup > 0 ){
1610
1611 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1612 double * Lblock = denL->gStorage( NL, TwoSL, ILxIixIl, NL+1, TwoSR, ILxIi );
1613
1614 char trans = 'T';
1615 char notrans = 'N';
1616 double alpha = 1.0;
1617 double beta = 0.0; //set
1618 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Sblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
1619 dgemm_( ¬rans, &trans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRdown, &beta, workmem2, &dimLdown );
1620
1621 int length = dimLdown * dimRdown;
1622 int inc = 1;
1623 total -= ( TwoSR + 1 ) * ddot_( &length, workmem2, &inc, Tdown, &inc );
1624
1625 }
1626 }
1627 }
1628 }
1629 }
1630 }
1631
1632 total *= sqrt( 0.5 );
1633 return total;
1634
1635 }
1636
diagram11(TensorT * denT,TensorS1 * denS1,TensorL * denL,double * workmem,double * workmem2) const1637 double CheMPS2::ThreeDM::diagram11(TensorT * denT, TensorS1 * denS1, TensorL * denL, double * workmem, double * workmem2) const{
1638
1639 const int orb_i = denT->gIndex();
1640 assert( denS1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1641
1642 double total = 0.0;
1643
1644 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1645 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1646 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1647
1648 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1649 const int ILxIixIl = Irreps::directProd( IL, denS1->get_irrep() );
1650
1651 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1652
1653 if ( dimLup > 0 ){
1654
1655 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
1656
1657 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSLprime, ILxIixIl );
1658 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSLprime, ILxIixIl );
1659
1660 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
1661
1662 double * Tdown = denT->gStorage( NL-2, TwoSLprime, ILxIixIl, NL, TwoSLprime, ILxIixIl );
1663 double * Sblock = denS1->gStorage( NL-2, TwoSLprime, ILxIixIl, NL, TwoSL, IL );
1664
1665 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1666
1667 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1668
1669 if (( dimRup > 0 ) && ( abs( TwoSLprime - TwoSR ) == 1 )){
1670
1671 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1672 double * Lblock = denL->gStorage( NL, TwoSLprime, ILxIixIl, NL+1, TwoSR, ILxIi );
1673
1674 char trans = 'T';
1675 char notrans = 'N';
1676 double alpha = 1.0;
1677 double beta = 0.0; //set
1678 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Sblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
1679 dgemm_( ¬rans, &trans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRdown, &beta, workmem2, &dimLdown );
1680
1681 int length = dimLdown * dimRdown;
1682 int inc = 1;
1683 total += sqrt( 3.0 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
1684 * Special::phase( TwoSR + 3 + TwoSLprime )
1685 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR )
1686 * ddot_( &length, workmem2, &inc, Tdown, &inc );
1687
1688 }
1689 }
1690 }
1691 }
1692 }
1693 }
1694 }
1695 }
1696
1697 return total;
1698
1699 }
1700
diagram12(TensorT * denT,TensorF0 * denF0,TensorL * denL,double * workmem,double * workmem2) const1701 double CheMPS2::ThreeDM::diagram12(TensorT * denT, TensorF0 * denF0, TensorL * denL, double * workmem, double * workmem2) const{
1702
1703 const int orb_i = denT->gIndex();
1704 assert( denF0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1705
1706 double total = 0.0;
1707
1708 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1709 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1710 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1711
1712 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1713 const int ILxIixIl = Irreps::directProd( IL, denF0->get_irrep() );
1714
1715 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1716 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSL, ILxIixIl );
1717 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSL, ILxIixIl );
1718
1719 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
1720
1721 double * Tdown = denT->gStorage( NL, TwoSL, ILxIixIl, NL+2, TwoSL, ILxIixIl );
1722 double * Fblock = denF0->gStorage( NL, TwoSL, ILxIixIl, NL, TwoSL, IL );
1723
1724 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1725
1726 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1727
1728 if ( dimRup > 0 ){
1729
1730 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1731 double * Lblock = denL->gStorage( NL+1, TwoSR, ILxIi, NL+2, TwoSL, ILxIixIl );
1732
1733 char notrans = 'N';
1734 double alpha = 1.0;
1735 double beta = 0.0; //set
1736 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Fblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
1737 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRup, &beta, workmem2, &dimLdown );
1738
1739 int length = dimLdown * dimRdown;
1740 int inc = 1;
1741 total += sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSR + 1 ) )
1742 * Special::phase( TwoSL + 3 - TwoSR )
1743 * ddot_( &length, workmem2, &inc, Tdown, &inc );
1744
1745 }
1746 }
1747 }
1748 }
1749 }
1750 }
1751
1752 return total;
1753
1754 }
1755
diagram13(TensorT * denT,TensorF1 * denF1,TensorL * denL,double * workmem,double * workmem2) const1756 double CheMPS2::ThreeDM::diagram13(TensorT * denT, TensorF1 * denF1, TensorL * denL, double * workmem, double * workmem2) const{
1757
1758 const int orb_i = denT->gIndex();
1759 assert( denF1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1760
1761 double total = 0.0;
1762
1763 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1764 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1765 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1766
1767 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1768 const int ILxIixIl = Irreps::directProd( IL, denF1->get_irrep() );
1769
1770 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1771
1772 if ( dimLup > 0 ){
1773
1774 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
1775
1776 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxIixIl );
1777 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSLprime, ILxIixIl );
1778
1779 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
1780
1781 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxIixIl, NL+2, TwoSLprime, ILxIixIl );
1782 double * Fblock = denF1->gStorage( NL, TwoSLprime, ILxIixIl, NL, TwoSL, IL );
1783
1784 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1785
1786 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1787
1788 if (( dimRup > 0 ) && ( abs( TwoSLprime - TwoSR ) == 1 )){
1789
1790 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1791 double * Lblock = denL->gStorage( NL+1, TwoSR, ILxIi, NL+2, TwoSLprime, ILxIixIl );
1792
1793 char notrans = 'N';
1794 double alpha = 1.0;
1795 double beta = 0.0; //set
1796 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Fblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
1797 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRup, &beta, workmem2, &dimLdown );
1798
1799 int length = dimLdown * dimRdown;
1800 int inc = 1;
1801 total += sqrt( 3.0 * ( TwoSL + 1 ) * ( TwoSR + 1 ) * ( TwoSLprime + 1 ) )
1802 * Special::phase( 2 * TwoSR + 2 )
1803 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR )
1804 * ddot_( &length, workmem2, &inc, Tdown, &inc );
1805
1806 }
1807 }
1808 }
1809 }
1810 }
1811 }
1812 }
1813 }
1814
1815 return total;
1816
1817 }
1818
diagram14(TensorT * denT,TensorF0 * denF0,TensorL * denL,double * workmem,double * workmem2) const1819 double CheMPS2::ThreeDM::diagram14(TensorT * denT, TensorF0 * denF0, TensorL * denL, double * workmem, double * workmem2) const{
1820
1821 const int orb_i = denT->gIndex();
1822 assert( denF0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1823
1824 double total = 0.0;
1825
1826 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1827 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1828 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1829
1830 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1831 const int ILxIixIl = Irreps::directProd( IL, denF0->get_irrep() );
1832
1833 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1834 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSL, ILxIixIl );
1835 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSL, ILxIixIl );
1836
1837 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
1838
1839 double * Tdown = denT->gStorage( NL, TwoSL, ILxIixIl, NL+2, TwoSL, ILxIixIl );
1840 double * Fblock = denF0->gStorage( NL, TwoSL, IL, NL, TwoSL, ILxIixIl );
1841
1842 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1843
1844 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1845
1846 if ( dimRup > 0 ){
1847
1848 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1849 double * Lblock = denL->gStorage( NL+1, TwoSR, ILxIi, NL+2, TwoSL, ILxIixIl );
1850
1851 char trans = 'T';
1852 char notrans = 'N';
1853 double alpha = 1.0;
1854 double beta = 0.0; //set
1855 dgemm_( &trans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Fblock, &dimLup, Tup, &dimLup, &beta, workmem, &dimLdown );
1856 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRup, &beta, workmem2, &dimLdown );
1857
1858 int length = dimLdown * dimRdown;
1859 int inc = 1;
1860 total += sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSR + 1 ) )
1861 * Special::phase( TwoSL + 3 - TwoSR )
1862 * ddot_( &length, workmem2, &inc, Tdown, &inc );
1863
1864 }
1865 }
1866 }
1867 }
1868 }
1869 }
1870
1871 return total;
1872
1873 }
1874
diagram15(TensorT * denT,TensorF1 * denF1,TensorL * denL,double * workmem,double * workmem2) const1875 double CheMPS2::ThreeDM::diagram15(TensorT * denT, TensorF1 * denF1, TensorL * denL, double * workmem, double * workmem2) const{
1876
1877 const int orb_i = denT->gIndex();
1878 assert( denF1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1879
1880 double total = 0.0;
1881
1882 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1883 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1884 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1885
1886 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1887 const int ILxIixIl = Irreps::directProd( IL, denF1->get_irrep() );
1888
1889 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1890
1891 if ( dimLup > 0 ){
1892
1893 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
1894
1895 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxIixIl );
1896 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSLprime, ILxIixIl );
1897
1898 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
1899
1900 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxIixIl, NL+2, TwoSLprime, ILxIixIl );
1901 double * Fblock = denF1->gStorage( NL, TwoSL, IL, NL, TwoSLprime, ILxIixIl );
1902
1903 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
1904
1905 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
1906
1907 if (( dimRup > 0 ) && ( abs( TwoSLprime - TwoSR ) == 1 )){
1908
1909 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
1910 double * Lblock = denL->gStorage( NL+1, TwoSR, ILxIi, NL+2, TwoSLprime, ILxIixIl );
1911
1912 char trans = 'T';
1913 char notrans = 'N';
1914 double alpha = 1.0;
1915 double beta = 0.0; //set
1916 dgemm_( &trans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Fblock, &dimLup, Tup, &dimLup, &beta, workmem, &dimLdown );
1917 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Lblock, &dimRup, &beta, workmem2, &dimLdown );
1918
1919 int length = dimLdown * dimRdown;
1920 int inc = 1;
1921 total += sqrt( 3.0 * ( TwoSR + 1 ) ) * ( TwoSLprime + 1 )
1922 * Special::phase( TwoSL + TwoSLprime )
1923 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR )
1924 * ddot_( &length, workmem2, &inc, Tdown, &inc );
1925
1926 }
1927 }
1928 }
1929 }
1930 }
1931 }
1932 }
1933 }
1934
1935 return total;
1936
1937 }
1938
diagram16(TensorT * denT,TensorL * denL,TensorS0 * denS0,double * workmem,double * workmem2) const1939 double CheMPS2::ThreeDM::diagram16(TensorT * denT, TensorL * denL, TensorS0 * denS0, double * workmem, double * workmem2) const{
1940
1941 const int orb_i = denT->gIndex();
1942 assert( denS0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1943
1944 double total = 0.0;
1945
1946 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
1947 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
1948 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
1949
1950 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
1951 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
1952
1953 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
1954
1955 if ( dimLup > 0 ){
1956
1957 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
1958
1959 int dimLdown = book->gCurrentDim( orb_i, NL+1, TwoSLprime, ILxIj );
1960 int dimRdown = book->gCurrentDim( orb_i+1, NL+3, TwoSLprime, ILxIj );
1961 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
1962
1963 if (( dimRup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
1964
1965 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
1966 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxIj, NL+3, TwoSLprime, ILxIj );
1967 double * Sblock = denS0->gStorage( NL+1, TwoSLprime, ILxIi, NL+3, TwoSLprime, ILxIj );
1968 double * Lblock = denL->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIj );
1969
1970 char trans = 'T';
1971 char notrans = 'N';
1972 double alpha = 1.0;
1973 double beta = 0.0; //set
1974 dgemm_( &trans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLup, Tup, &dimLup, &beta, workmem, &dimLdown );
1975 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Sblock, &dimRup, &beta, workmem2, &dimLdown );
1976
1977 int length = dimLdown * dimRdown;
1978 int inc = 1;
1979 total -= ( TwoSLprime + 1 ) * ddot_( &length, workmem2, &inc, Tdown, &inc );
1980
1981 }
1982 }
1983 }
1984 }
1985 }
1986 }
1987
1988 total *= sqrt( 0.5 );
1989 return total;
1990
1991 }
1992
diagram17(TensorT * denT,TensorL * denL,TensorS1 * denS1,double * workmem,double * workmem2) const1993 double CheMPS2::ThreeDM::diagram17(TensorT * denT, TensorL * denL, TensorS1 * denS1, double * workmem, double * workmem2) const{
1994
1995 const int orb_i = denT->gIndex();
1996 assert( denS1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
1997
1998 double total = 0.0;
1999
2000 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2001 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2002 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2003
2004 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2005 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
2006
2007 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2008
2009 if ( dimLup > 0 ){
2010
2011 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2012
2013 int dimLdown = book->gCurrentDim( orb_i, NL+1, TwoSLprime, ILxIj );
2014 int dimRdown = book->gCurrentDim( orb_i+1, NL+3, TwoSLprime, ILxIj );
2015
2016 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
2017
2018 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxIj, NL+3, TwoSLprime, ILxIj );
2019 double * Lblock = denL->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIj );
2020
2021 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2022
2023 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2024
2025 if ( dimRup > 0 ){
2026
2027 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2028 double * Sblock = denS1->gStorage( NL+1, TwoSR, ILxIi, NL+3, TwoSLprime, ILxIj );
2029
2030 char trans = 'T';
2031 char notrans = 'N';
2032 double alpha = 1.0;
2033 double beta = 0.0; //set
2034 dgemm_( &trans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLup, Tup, &dimLup, &beta, workmem, &dimLdown );
2035 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Sblock, &dimRup, &beta, workmem2, &dimLdown );
2036
2037 int length = dimLdown * dimRdown;
2038 int inc = 1;
2039 total += sqrt( 3.0 * ( TwoSR + 1 ) ) * ( TwoSLprime + 1 )
2040 * Special::phase( TwoSL + TwoSLprime + 3 )
2041 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL )
2042 * ddot_( &length, workmem2, &inc, Tdown, &inc );
2043
2044 }
2045 }
2046 }
2047 }
2048 }
2049 }
2050 }
2051 }
2052
2053 return total;
2054
2055 }
2056
diagram18(TensorT * denT,TensorL * denL,TensorF0 * denF0,double * workmem,double * workmem2) const2057 double CheMPS2::ThreeDM::diagram18(TensorT * denT, TensorL * denL, TensorF0 * denF0, double * workmem, double * workmem2) const{
2058
2059 const int orb_i = denT->gIndex();
2060 assert( denF0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
2061
2062 double total = 0.0;
2063
2064 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2065 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2066 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2067
2068 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2069 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
2070
2071 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2072
2073 if ( dimLup > 0 ){
2074
2075 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2076
2077 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIj );
2078 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIj );
2079 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2080
2081 if (( dimRup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
2082
2083 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2084 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIj, NL+1, TwoSLprime, ILxIj );
2085 double * Fblock = denF0->gStorage( NL+1, TwoSLprime, ILxIj, NL+1, TwoSLprime, ILxIi );
2086 double * Lblock = denL->gStorage( NL-1, TwoSLprime, ILxIj, NL, TwoSL, IL );
2087
2088 char trans = 'T';
2089 char notrans = 'N';
2090 double alpha = 1.0;
2091 double beta = 0.0; //set
2092 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
2093 dgemm_( ¬rans, &trans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Fblock, &dimRdown, &beta, workmem2, &dimLdown );
2094
2095 int length = dimLdown * dimRdown;
2096 int inc = 1;
2097 total += sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) )
2098 * Special::phase( TwoSL + 1 - TwoSLprime )
2099 * ddot_( &length, workmem2, &inc, Tdown, &inc );
2100
2101 }
2102 }
2103 }
2104 }
2105 }
2106 }
2107
2108 return total;
2109
2110 }
2111
diagram19(TensorT * denT,TensorL * denL,TensorF1 * denF1,double * workmem,double * workmem2) const2112 double CheMPS2::ThreeDM::diagram19(TensorT * denT, TensorL * denL, TensorF1 * denF1, double * workmem, double * workmem2) const{
2113
2114 const int orb_i = denT->gIndex();
2115 assert( denF1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
2116
2117 double total = 0.0;
2118
2119 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2120 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2121 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2122
2123 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2124 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
2125
2126 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2127
2128 if ( dimLup > 0 ){
2129
2130 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2131
2132 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIj );
2133 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIj );
2134
2135 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
2136
2137 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIj, NL+1, TwoSLprime, ILxIj );
2138 double * Lblock = denL->gStorage( NL-1, TwoSLprime, ILxIj, NL, TwoSL, IL );
2139
2140 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2141
2142 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2143
2144 if ( dimRup > 0 ){
2145
2146 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2147 double * Fblock = denF1->gStorage( NL+1, TwoSLprime, ILxIj, NL+1, TwoSR, ILxIi );
2148
2149 char trans = 'T';
2150 char notrans = 'N';
2151 double alpha = 1.0;
2152 double beta = 0.0; //set
2153 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
2154 dgemm_( ¬rans, &trans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Fblock, &dimRdown, &beta, workmem2, &dimLdown );
2155
2156 int length = dimLdown * dimRdown;
2157 int inc = 1;
2158 total += sqrt( 3.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) * ( TwoSR + 1 ) )
2159 * Special::phase( 2 * TwoSR )
2160 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL )
2161 * ddot_( &length, workmem2, &inc, Tdown, &inc );
2162
2163 }
2164 }
2165 }
2166 }
2167 }
2168 }
2169 }
2170 }
2171
2172 return total;
2173
2174 }
2175
diagram20(TensorT * denT,TensorL * denL,TensorF0 * denF0,double * workmem,double * workmem2) const2176 double CheMPS2::ThreeDM::diagram20(TensorT * denT, TensorL * denL, TensorF0 * denF0, double * workmem, double * workmem2) const{
2177
2178 const int orb_i = denT->gIndex();
2179 assert( denF0->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
2180
2181 double total = 0.0;
2182
2183 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2184 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2185 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2186
2187 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2188 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
2189
2190 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2191
2192 if ( dimLup > 0 ){
2193
2194 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2195
2196 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIj );
2197 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIj );
2198 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2199
2200 if (( dimRup > 0 ) && ( dimLdown > 0 ) && ( dimRdown > 0 )){
2201
2202 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2203 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIj, NL+1, TwoSLprime, ILxIj );
2204 double * Fblock = denF0->gStorage( NL+1, TwoSLprime, ILxIi, NL+1, TwoSLprime, ILxIj );
2205 double * Lblock = denL->gStorage( NL-1, TwoSLprime, ILxIj, NL, TwoSL, IL );
2206
2207 char notrans = 'N';
2208 double alpha = 1.0;
2209 double beta = 0.0; //set
2210 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
2211 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Fblock, &dimRup, &beta, workmem2, &dimLdown );
2212
2213 int length = dimLdown * dimRdown;
2214 int inc = 1;
2215 total += sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) )
2216 * Special::phase( TwoSL + 1 - TwoSLprime )
2217 * ddot_( &length, workmem2, &inc, Tdown, &inc );
2218
2219 }
2220 }
2221 }
2222 }
2223 }
2224 }
2225
2226 return total;
2227
2228 }
2229
diagram21(TensorT * denT,TensorL * denL,TensorF1 * denF1,double * workmem,double * workmem2) const2230 double CheMPS2::ThreeDM::diagram21(TensorT * denT, TensorL * denL, TensorF1 * denF1, double * workmem, double * workmem2) const{
2231
2232 const int orb_i = denT->gIndex();
2233 assert( denF1->get_irrep() == Irreps::directProd( book->gIrrep( orb_i ), denL->get_irrep() ) );
2234
2235 double total = 0.0;
2236
2237 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2238 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2239 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2240
2241 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2242 const int ILxIj = Irreps::directProd( IL, denL->get_irrep() );
2243
2244 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2245
2246 if ( dimLup > 0 ){
2247
2248 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2249
2250 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIj );
2251 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIj );
2252
2253 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
2254
2255 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIj, NL+1, TwoSLprime, ILxIj );
2256 double * Lblock = denL->gStorage( NL-1, TwoSLprime, ILxIj, NL, TwoSL, IL );
2257
2258 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2259
2260 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2261
2262 if ( dimRup > 0 ){
2263
2264 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2265 double * Fblock = denF1->gStorage( NL+1, TwoSR, ILxIi, NL+1, TwoSLprime, ILxIj );
2266
2267 char notrans = 'N';
2268 double alpha = 1.0;
2269 double beta = 0.0; //set
2270 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimLup, &alpha, Lblock, &dimLdown, Tup, &dimLup, &beta, workmem, &dimLdown );
2271 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRdown, &dimRup, &alpha, workmem, &dimLdown, Fblock, &dimRup, &beta, workmem2, &dimLdown );
2272
2273 int length = dimLdown * dimRdown;
2274 int inc = 1;
2275 total += sqrt( 3.0 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
2276 * Special::phase( TwoSR + TwoSLprime )
2277 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL )
2278 * ddot_( &length, workmem2, &inc, Tdown, &inc );
2279
2280 }
2281 }
2282 }
2283 }
2284 }
2285 }
2286 }
2287 }
2288
2289 return total;
2290
2291 }
2292
fill_a_S0(TensorT * denT,Tensor3RDM * tofill,TensorS0 * denS0,double * workmem) const2293 void CheMPS2::ThreeDM::fill_a_S0( TensorT * denT, Tensor3RDM * tofill, TensorS0 * denS0, double * workmem ) const{
2294
2295 const int orb_i = denT->gIndex();
2296 const int ImxInxIi = Irreps::directProd( denS0->get_irrep(), book->gIrrep( orb_i ) );
2297 assert( tofill->get_irrep() == ImxInxIi );
2298 assert( tofill->get_nelec() == 3 );
2299 assert( tofill->get_two_j2() == 1 );
2300
2301 tofill->clear();
2302
2303 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2304 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2305 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2306
2307 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2308 const int ILxImxIn = Irreps::directProd( IL, denS0->get_irrep() );
2309 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2310
2311 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2312
2313 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2314 int dimLdown = book->gCurrentDim( orb_i, NL-3, TwoSLprime, ILxImxInxIi );
2315
2316 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2317
2318 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2319 int dimRdown = book->gCurrentDim( orb_i+1, NL-2, TwoSL, ILxImxIn );
2320
2321 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2322
2323 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2324 double * Tdown = denT->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL-2, TwoSL, ILxImxIn );
2325 double * Sblock = denS0->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, IL );
2326 double * Wblock = tofill->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2327
2328 char notrans = 'N';
2329 char trans = 'T';
2330 double alpha = -0.5 * ( TwoSL + 1 );
2331 double beta = 0.0; //SET
2332 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Sblock, &dimRdown, &beta, workmem, &dimLdown );
2333 alpha = 1.0;
2334 beta = 1.0; //ADD
2335 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2336
2337 }
2338
2339 dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2340 dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSLprime, ILxImxInxIi );
2341
2342 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2343
2344 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2345 double * Tdown = denT->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL-1, TwoSLprime, ILxImxInxIi );
2346 double * Sblock = denS0->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxIi );
2347 double * Wblock = tofill->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2348
2349 char notrans = 'N';
2350 char trans = 'T';
2351 double alpha = 0.5 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSL + 1 - TwoSLprime );
2352 double beta = 0.0; //SET
2353 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Sblock, &dimRdown, &beta, workmem, &dimLdown );
2354 alpha = 1.0;
2355 beta = 1.0; //ADD
2356 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2357
2358 }
2359 }
2360 }
2361 }
2362 }
2363 }
2364
2365 }
2366
fill_bcd_S0(TensorT * denT,Tensor3RDM * tofill,TensorS0 * denS0,double * workmem) const2367 void CheMPS2::ThreeDM::fill_bcd_S0( TensorT * denT, Tensor3RDM * tofill, TensorS0 * denS0, double * workmem ) const{
2368
2369 const int orb_i = denT->gIndex();
2370 const int ImxInxIi = Irreps::directProd( denS0->get_irrep(), book->gIrrep( orb_i ) );
2371 assert( tofill->get_irrep() == ImxInxIi );
2372 assert( tofill->get_nelec() == 1 );
2373 assert( tofill->get_two_j2() == 1 );
2374
2375 tofill->clear();
2376
2377 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2378 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2379 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2380
2381 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2382 const int ILxImxIn = Irreps::directProd( IL, denS0->get_irrep() );
2383 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2384
2385 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2386
2387 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2388 int dimLdown = book->gCurrentDim( orb_i, NL+1, TwoSLprime, ILxImxInxIi );
2389
2390 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2391
2392 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2393 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSL, ILxImxIn );
2394
2395 if (( dimRup > 0 ) && ( dimRdown > 0 )){ // Type 100, 102, 110, 112, 120, 124
2396
2397 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2398 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+2, TwoSL, ILxImxIn );
2399 double * Sblock = denS0->gStorage( NL, TwoSL, IL, NL+2, TwoSL, ILxImxIn );
2400 double * Wblock = tofill->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2401
2402 char notrans = 'N';
2403 char trans = 'T';
2404 double alpha = 0.5 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSL + 1 - TwoSLprime );
2405 double beta = 0.0; //SET
2406 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimRup, &alpha, Tup, &dimLup, Sblock, &dimRup, &beta, workmem, &dimLup );
2407 alpha = 1.0;
2408 beta = 1.0; //ADD
2409 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRdown, &alpha, workmem, &dimLup, Tdown, &dimLdown, &beta, Wblock, &dimLup );
2410
2411 }
2412
2413 dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2414 dimRdown = book->gCurrentDim( orb_i+1, NL+3, TwoSLprime, ILxImxInxIi );
2415
2416 if (( dimRup > 0 ) && ( dimRdown > 0 )){ // Type 105, 107, 115, 117, 125, 129
2417
2418 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2419 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+3, TwoSLprime, ILxImxInxIi );
2420 double * Sblock = denS0->gStorage( NL+1, TwoSLprime, ILxIi , NL+3, TwoSLprime, ILxImxInxIi );
2421 double * Wblock = tofill->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2422
2423 char notrans = 'N';
2424 char trans = 'T';
2425 double alpha = - 0.5 * ( TwoSLprime + 1 );
2426 double beta = 0.0; //SET
2427 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimRup, &alpha, Tup, &dimLup, Sblock, &dimRup, &beta, workmem, &dimLup );
2428 alpha = 1.0;
2429 beta = 1.0; //ADD
2430 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRdown, &alpha, workmem, &dimLup, Tdown, &dimLdown, &beta, Wblock, &dimLup );
2431
2432 }
2433 }
2434 }
2435 }
2436 }
2437 }
2438
2439 }
2440
fill_a_S1(TensorT * denT,Tensor3RDM * doublet,Tensor3RDM * quartet,TensorS1 * denS1,double * workmem,double * workmem2) const2441 void CheMPS2::ThreeDM::fill_a_S1( TensorT * denT, Tensor3RDM * doublet, Tensor3RDM * quartet, TensorS1 * denS1, double * workmem, double * workmem2 ) const{
2442
2443 const int orb_i = denT->gIndex();
2444 const int ImxInxIi = Irreps::directProd( denS1->get_irrep(), book->gIrrep( orb_i ) );
2445 assert( doublet->get_irrep() == ImxInxIi ); assert( doublet->get_nelec() == 3 ); assert( doublet->get_two_j2() == 1 );
2446 assert( quartet->get_irrep() == ImxInxIi ); assert( quartet->get_nelec() == 3 ); assert( quartet->get_two_j2() == 3 );
2447
2448 doublet->clear();
2449 quartet->clear();
2450
2451 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2452 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2453 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2454
2455 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2456 const int ILxImxIn = Irreps::directProd( IL, denS1->get_irrep() );
2457 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2458
2459 for ( int TwoSLprime = TwoSL-3; TwoSLprime <= TwoSL+3; TwoSLprime+=2 ){
2460
2461 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2462 int dimLdown = book->gCurrentDim( orb_i, NL-3, TwoSLprime, ILxImxInxIi );
2463
2464 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2465
2466 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
2467
2468 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2469 int dimRdown = book->gCurrentDim( orb_i+1, NL-2, TwoSRprime, ILxImxIn );
2470
2471 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSL - TwoSRprime ) <= 2 )){
2472
2473 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2474 double * Tdown = denT->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL-2, TwoSRprime, ILxImxIn );
2475 double * Sblock = denS1->gStorage( NL-2, TwoSRprime, ILxImxIn, NL, TwoSL, IL );
2476
2477 char notrans = 'N';
2478 char trans = 'T';
2479 double alpha = 1.0;
2480 double beta = 0.0; //SET
2481 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Sblock, &dimRdown, &beta, workmem, &dimLdown );
2482 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
2483
2484 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2485
2486 double * Wblock = doublet->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2487 double prefactor = sqrt( 0.5 * ( TwoSRprime + 1 ) ) * ( TwoSL + 1 )
2488 * Special::phase( TwoSL + TwoSLprime + 1 )
2489 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSLprime );
2490 int length = dimLup * dimLdown;
2491 int inc = 1;
2492 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2493
2494 }
2495 {
2496
2497 double * Wblock = quartet->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2498 double prefactor = 2 * sqrt( TwoSRprime + 1.0 ) * ( TwoSL + 1 )
2499 * Special::phase( TwoSL + TwoSLprime + 3 )
2500 * Wigner::wigner6j( 1, 3, 2, TwoSL, TwoSRprime, TwoSLprime );
2501 int length = dimLup * dimLdown;
2502 int inc = 1;
2503 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2504
2505 }
2506
2507 }
2508 }
2509 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2510
2511 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2512 int dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSLprime, ILxImxInxIi );
2513
2514 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) <= 2 )){
2515
2516 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2517 double * Tdown = denT->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL-1, TwoSLprime, ILxImxInxIi );
2518 double * Sblock = denS1->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSR, ILxIi );
2519
2520 char notrans = 'N';
2521 char trans = 'T';
2522 double alpha = 1.0;
2523 double beta = 0.0; //SET
2524 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Sblock, &dimRdown, &beta, workmem, &dimLdown );
2525 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
2526
2527 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2528
2529 double * Wblock = doublet->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2530 double prefactor = sqrt( 0.5 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
2531 * Special::phase( TwoSR + TwoSLprime )
2532 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL );
2533 int length = dimLup * dimLdown;
2534 int inc = 1;
2535 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2536
2537 }
2538 {
2539
2540 double * Wblock = quartet->gStorage( NL-3, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2541 double prefactor = 2 * sqrt( TwoSL + 1.0 ) * ( TwoSR + 1 )
2542 * Special::phase( TwoSR + TwoSLprime )
2543 * Wigner::wigner6j( 1, 3, 2, TwoSLprime, TwoSR, TwoSL );
2544 int length = dimLup * dimLdown;
2545 int inc = 1;
2546 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2547
2548 }
2549
2550 }
2551 }
2552 }
2553 }
2554 }
2555 }
2556 }
2557
2558 }
2559
fill_bcd_S1(TensorT * denT,Tensor3RDM * doublet,Tensor3RDM * quartet,TensorS1 * denS1,double * workmem,double * workmem2) const2560 void CheMPS2::ThreeDM::fill_bcd_S1( TensorT * denT, Tensor3RDM * doublet, Tensor3RDM * quartet, TensorS1 * denS1, double * workmem, double * workmem2 ) const{
2561
2562 const int orb_i = denT->gIndex();
2563 const int ImxInxIi = Irreps::directProd( denS1->get_irrep(), book->gIrrep( orb_i ) );
2564 assert( doublet->get_irrep() == ImxInxIi ); assert( doublet->get_nelec() == 1 ); assert( doublet->get_two_j2() == 1 );
2565 assert( quartet->get_irrep() == ImxInxIi ); assert( quartet->get_nelec() == 1 ); assert( quartet->get_two_j2() == 3 );
2566
2567 doublet->clear();
2568 quartet->clear();
2569
2570 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2571 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2572 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2573
2574 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2575 const int ILxImxIn = Irreps::directProd( IL, denS1->get_irrep() );
2576 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2577
2578 for ( int TwoSLprime = TwoSL-3; TwoSLprime <= TwoSL+3; TwoSLprime+=2 ){
2579
2580 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2581 int dimLdown = book->gCurrentDim( orb_i, NL+1, TwoSLprime, ILxImxInxIi );
2582
2583 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2584
2585 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
2586
2587 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2588 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSRprime, ILxImxIn );
2589
2590 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSL - TwoSRprime ) <= 2 )){
2591
2592 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2593 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+2, TwoSRprime, ILxImxIn );
2594 double * Sblock = denS1->gStorage( NL, TwoSL, IL, NL+2, TwoSRprime, ILxImxIn );
2595
2596 char notrans = 'N';
2597 char trans = 'T';
2598 double alpha = 1.0;
2599 double beta = 0.0; //SET
2600 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimRup, &alpha, Tup, &dimLup, Sblock, &dimRup, &beta, workmem, &dimLup );
2601 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRdown, &alpha, workmem, &dimLup, Tdown, &dimLdown, &beta, workmem2, &dimLup );
2602
2603 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2604
2605 double * Wblock = doublet->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2606 double prefactor = sqrt( 0.5 * ( TwoSLprime + 1 ) ) * ( TwoSRprime + 1 )
2607 * Special::phase( TwoSL + TwoSRprime )
2608 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSLprime );
2609 int length = dimLup * dimLdown;
2610 int inc = 1;
2611 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2612
2613 }
2614 {
2615
2616 double * Wblock = quartet->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2617 double prefactor = sqrt( TwoSLprime + 1.0 ) * ( TwoSRprime + 1 )
2618 * Special::phase( TwoSL + TwoSRprime )
2619 * Wigner::wigner6j( 1, 2, 3, TwoSL, TwoSLprime, TwoSRprime );
2620 int length = dimLup * dimLdown;
2621 int inc = 1;
2622 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2623
2624 }
2625
2626 }
2627 }
2628 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2629
2630 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2631 int dimRdown = book->gCurrentDim( orb_i+1, NL+3, TwoSLprime, ILxImxInxIi );
2632
2633 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) <= 2 )){
2634
2635 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2636 double * Tdown = denT->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+3, TwoSLprime, ILxImxInxIi );
2637 double * Sblock = denS1->gStorage( NL+1, TwoSR, ILxIi , NL+3, TwoSLprime, ILxImxInxIi );
2638
2639 char notrans = 'N';
2640 char trans = 'T';
2641 double alpha = 1.0;
2642 double beta = 0.0; //SET
2643 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimRup, &alpha, Tup, &dimLup, Sblock, &dimRup, &beta, workmem, &dimLup );
2644 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRdown, &alpha, workmem, &dimLup, Tdown, &dimLdown, &beta, workmem2, &dimLup );
2645
2646 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2647
2648 double * Wblock = doublet->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2649 double prefactor = sqrt( 0.5 * ( TwoSR + 1 ) ) * ( TwoSLprime + 1 )
2650 * Special::phase( TwoSL + TwoSLprime + 3 )
2651 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL );
2652 int length = dimLup * dimLdown;
2653 int inc = 1;
2654 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2655
2656 }
2657 {
2658
2659 double * Wblock = quartet->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxImxInxIi );
2660 double prefactor = sqrt( TwoSR + 1.0 ) * ( TwoSLprime + 1 )
2661 * Special::phase( TwoSL + TwoSLprime + 1 )
2662 * Wigner::wigner6j( 1, 3, 2, TwoSLprime, TwoSR, TwoSL );
2663 int length = dimLup * dimLdown;
2664 int inc = 1;
2665 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2666
2667 }
2668
2669 }
2670 }
2671 }
2672 }
2673 }
2674 }
2675 }
2676
2677 }
2678
fill_F0_T(TensorT * denT,Tensor3RDM * tofill,TensorF0 * denF0,double * workmem) const2679 void CheMPS2::ThreeDM::fill_F0_T( TensorT * denT, Tensor3RDM * tofill, TensorF0 * denF0, double * workmem ) const{
2680
2681 const int orb_i = denT->gIndex();
2682 const int ImxInxIi = Irreps::directProd( denF0->get_irrep(), book->gIrrep( orb_i ) );
2683 assert( tofill->get_irrep() == ImxInxIi );
2684 assert( tofill->get_nelec() == 1 );
2685 assert( tofill->get_two_j2() == 1 );
2686
2687 tofill->clear();
2688
2689 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2690 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2691 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2692
2693 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2694 const int ILxImxIn = Irreps::directProd( IL, denF0->get_irrep() );
2695 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2696
2697 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2698
2699 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2700 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxImxInxIi );
2701
2702 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2703
2704 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2705 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxImxIn );
2706
2707 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2708
2709 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2710 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, ILxImxIn );
2711 double * Fblock = denF0->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, IL );
2712 double * Wblock = tofill->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2713
2714 char notrans = 'N';
2715 char trans = 'T';
2716 double alpha = 0.5 * ( TwoSL + 1 );
2717 double beta = 0.0; //SET
2718 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRdown, &beta, workmem, &dimLdown );
2719 alpha = 1.0;
2720 beta = 1.0; //ADD
2721 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2722
2723 }
2724
2725 dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2726 dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxImxInxIi );
2727
2728 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2729
2730 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2731 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxImxInxIi );
2732 double * Fblock = denF0->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxIi );
2733 double * Wblock = tofill->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2734
2735 char notrans = 'N';
2736 char trans = 'T';
2737 double alpha = 0.5 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSLprime + 1 - TwoSL );
2738 double beta = 0.0; //SET
2739 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRdown, &beta, workmem, &dimLdown );
2740 alpha = 1.0;
2741 beta = 1.0; //ADD
2742 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2743
2744 }
2745 }
2746 }
2747 }
2748 }
2749 }
2750
2751 }
2752
fill_F1_T(TensorT * denT,Tensor3RDM * doublet,Tensor3RDM * quartet,TensorF1 * denF1,double * workmem,double * workmem2) const2753 void CheMPS2::ThreeDM::fill_F1_T( TensorT * denT, Tensor3RDM * doublet, Tensor3RDM * quartet, TensorF1 * denF1, double * workmem, double * workmem2 ) const{
2754
2755 const int orb_i = denT->gIndex();
2756 const int ImxInxIi = Irreps::directProd( denF1->get_irrep(), book->gIrrep( orb_i ) );
2757 assert( doublet->get_irrep() == ImxInxIi ); assert( doublet->get_nelec() == 1 ); assert( doublet->get_two_j2() == 1 );
2758 assert( quartet->get_irrep() == ImxInxIi ); assert( quartet->get_nelec() == 1 ); assert( quartet->get_two_j2() == 3 );
2759
2760 doublet->clear();
2761 quartet->clear();
2762
2763 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2764 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2765 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2766
2767 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2768 const int ILxImxIn = Irreps::directProd( IL, denF1->get_irrep() );
2769 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2770
2771 for ( int TwoSLprime = TwoSL-3; TwoSLprime <= TwoSL+3; TwoSLprime+=2 ){
2772
2773 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2774 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxImxInxIi );
2775
2776 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2777
2778 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
2779
2780 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2781 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSRprime, ILxImxIn );
2782
2783 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSL - TwoSRprime ) <= 2 )){
2784
2785 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2786 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSRprime, ILxImxIn );
2787 double * Fblock = denF1->gStorage( NL, TwoSRprime, ILxImxIn, NL, TwoSL, IL );
2788
2789 char notrans = 'N';
2790 char trans = 'T';
2791 double alpha = 1.0;
2792 double beta = 0.0; //SET
2793 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRdown, &beta, workmem, &dimLdown );
2794 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
2795
2796 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2797
2798 double * Wblock = doublet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2799 double prefactor = sqrt( 0.5 * ( TwoSL + 1 ) ) * ( TwoSRprime + 1 )
2800 * Special::phase( TwoSLprime + TwoSRprime + 3 )
2801 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSLprime );
2802 int length = dimLup * dimLdown;
2803 int inc = 1;
2804 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2805
2806 }
2807 {
2808
2809 double * Wblock = quartet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2810 double prefactor = sqrt( TwoSL + 1.0 ) * ( TwoSRprime + 1 )
2811 * Special::phase( TwoSLprime + TwoSRprime + 3 )
2812 * Wigner::wigner6j( 1, 3, 2, TwoSL, TwoSRprime, TwoSLprime );
2813 int length = dimLup * dimLdown;
2814 int inc = 1;
2815 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2816
2817 }
2818
2819 }
2820 }
2821 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
2822
2823 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
2824 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxImxInxIi );
2825
2826 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) <= 2 )){
2827
2828 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
2829 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxImxInxIi );
2830 double * Fblock = denF1->gStorage( NL+1, TwoSLprime, ILxImxInxIi, NL+1, TwoSR, ILxIi );
2831
2832 char notrans = 'N';
2833 char trans = 'T';
2834 double alpha = 1.0;
2835 double beta = 0.0; //SET
2836 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRdown, &beta, workmem, &dimLdown );
2837 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
2838
2839 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2840
2841 double * Wblock = doublet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2842 double prefactor = sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) * ( TwoSR + 1 ) )
2843 * Special::phase( 2*TwoSLprime + 2 )
2844 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL );
2845 int length = dimLup * dimLdown;
2846 int inc = 1;
2847 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2848
2849 }
2850 {
2851
2852 double * Wblock = quartet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2853 double prefactor = sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) * ( TwoSR + 1 ) )
2854 * Special::phase( 2*TwoSLprime )
2855 * Wigner::wigner6j( 1, 3, 2, TwoSLprime, TwoSR, TwoSL );
2856 int length = dimLup * dimLdown;
2857 int inc = 1;
2858 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2859
2860 }
2861
2862 }
2863 }
2864 }
2865 }
2866 }
2867 }
2868 }
2869
2870 }
2871
fill_F0(TensorT * denT,Tensor3RDM * tofill,TensorF0 * denF0,double * workmem) const2872 void CheMPS2::ThreeDM::fill_F0( TensorT * denT, Tensor3RDM * tofill, TensorF0 * denF0, double * workmem ) const{
2873
2874 const int orb_i = denT->gIndex();
2875 const int ImxInxIi = Irreps::directProd( denF0->get_irrep(), book->gIrrep( orb_i ) );
2876 assert( tofill->get_irrep() == ImxInxIi );
2877 assert( tofill->get_nelec() == 1 );
2878 assert( tofill->get_two_j2() == 1 );
2879
2880 tofill->clear();
2881
2882 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2883 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2884 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2885
2886 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2887 const int ILxImxIn = Irreps::directProd( IL, denF0->get_irrep() );
2888 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2889
2890 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
2891
2892 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2893 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxImxInxIi );
2894
2895 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2896
2897 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2898 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxImxIn );
2899
2900 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2901
2902 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2903 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, ILxImxIn );
2904 double * Fblock = denF0->gStorage( NL, TwoSL, IL, NL, TwoSL, ILxImxIn );
2905 double * Wblock = tofill->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2906
2907 char notrans = 'N';
2908 char trans = 'T';
2909 double alpha = 0.5 * ( TwoSL + 1 );
2910 double beta = 0.0; //SET
2911 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRup, &beta, workmem, &dimLdown );
2912 alpha = 1.0;
2913 beta = 1.0; //ADD
2914 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2915
2916 }
2917
2918 dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIi );
2919 dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxImxInxIi );
2920
2921 if (( dimRup > 0 ) && ( dimRdown > 0 )){
2922
2923 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIi );
2924 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxImxInxIi );
2925 double * Fblock = denF0->gStorage( NL+1, TwoSLprime, ILxIi, NL+1, TwoSLprime, ILxImxInxIi );
2926 double * Wblock = tofill->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2927
2928 char notrans = 'N';
2929 char trans = 'T';
2930 double alpha = 0.5 * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSLprime + 1 - TwoSL );
2931 double beta = 0.0; //SET
2932 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRup, &beta, workmem, &dimLdown );
2933 alpha = 1.0;
2934 beta = 1.0; //ADD
2935 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, Wblock, &dimLdown );
2936
2937 }
2938 }
2939 }
2940 }
2941 }
2942 }
2943
2944 }
2945
fill_F1(TensorT * denT,Tensor3RDM * doublet,Tensor3RDM * quartet,TensorF1 * denF1,double * workmem,double * workmem2) const2946 void CheMPS2::ThreeDM::fill_F1( TensorT * denT, Tensor3RDM * doublet, Tensor3RDM * quartet, TensorF1 * denF1, double * workmem, double * workmem2 ) const{
2947
2948 const int orb_i = denT->gIndex();
2949 const int ImxInxIi = Irreps::directProd( denF1->get_irrep(), book->gIrrep( orb_i ) );
2950 assert( doublet->get_irrep() == ImxInxIi ); assert( doublet->get_nelec() == 1 ); assert( doublet->get_two_j2() == 1 );
2951 assert( quartet->get_irrep() == ImxInxIi ); assert( quartet->get_nelec() == 1 ); assert( quartet->get_two_j2() == 3 );
2952
2953 doublet->clear();
2954 quartet->clear();
2955
2956 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
2957 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
2958 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
2959
2960 const int ILxImxInxIi = Irreps::directProd( IL, ImxInxIi );
2961 const int ILxImxIn = Irreps::directProd( IL, denF1->get_irrep() );
2962 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
2963
2964 for ( int TwoSLprime = TwoSL-3; TwoSLprime <= TwoSL+3; TwoSLprime+=2 ){
2965
2966 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
2967 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxImxInxIi );
2968
2969 if (( dimLup > 0 ) && ( dimLdown > 0 )){
2970
2971 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
2972
2973 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
2974 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSRprime, ILxImxIn );
2975
2976 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSL - TwoSRprime ) <= 2 )){
2977
2978 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
2979 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSRprime, ILxImxIn );
2980 double * Fblock = denF1->gStorage( NL, TwoSL, IL, NL, TwoSRprime, ILxImxIn );
2981
2982 char notrans = 'N';
2983 char trans = 'T';
2984 double alpha = 1.0;
2985 double beta = 0.0; //SET
2986 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRup, &beta, workmem, &dimLdown );
2987 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
2988
2989 if ( abs( TwoSL - TwoSLprime ) == 1 ){
2990
2991 double * Wblock = doublet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
2992 double prefactor = sqrt( 0.5 * ( TwoSRprime + 1 ) ) * ( TwoSL + 1 )
2993 * Special::phase( TwoSLprime + TwoSL + 3 )
2994 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSLprime );
2995 int length = dimLup * dimLdown;
2996 int inc = 1;
2997 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
2998
2999 }
3000 {
3001
3002 double * Wblock = quartet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
3003 double prefactor = sqrt( TwoSRprime + 1.0 ) * ( TwoSL + 1 )
3004 * Special::phase( TwoSLprime + TwoSL + 3 )
3005 * Wigner::wigner6j( 1, 3, 2, TwoSL, TwoSRprime, TwoSLprime );
3006 int length = dimLup * dimLdown;
3007 int inc = 1;
3008 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
3009
3010 }
3011
3012 }
3013 }
3014 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3015
3016 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3017 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxImxInxIi );
3018
3019 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) <= 2 )){
3020
3021 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3022 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL+1, TwoSLprime, ILxImxInxIi );
3023 double * Fblock = denF1->gStorage( NL+1, TwoSR, ILxIi, NL+1, TwoSLprime, ILxImxInxIi );
3024
3025 char notrans = 'N';
3026 char trans = 'T';
3027 double alpha = 1.0;
3028 double beta = 0.0; //SET
3029 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &alpha, Tdown, &dimLdown, Fblock, &dimRup, &beta, workmem, &dimLdown );
3030 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &alpha, workmem, &dimLdown, Tup, &dimLup, &beta, workmem2, &dimLdown );
3031
3032 if ( abs( TwoSL - TwoSLprime ) == 1 ){
3033
3034 double * Wblock = doublet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
3035 double prefactor = sqrt( 0.5 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
3036 * Special::phase( TwoSLprime + TwoSR + 2 )
3037 * Wigner::wigner6j( 1, 1, 2, TwoSLprime, TwoSR, TwoSL );
3038 int length = dimLup * dimLdown;
3039 int inc = 1;
3040 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
3041
3042 }
3043 {
3044
3045 double * Wblock = quartet->gStorage( NL-1, TwoSLprime, ILxImxInxIi, NL, TwoSL, IL );
3046 double prefactor = sqrt( TwoSL + 1.0 ) * ( TwoSR + 1 )
3047 * Special::phase( TwoSLprime + TwoSR )
3048 * Wigner::wigner6j( 1, 3, 2, TwoSLprime, TwoSR, TwoSL );
3049 int length = dimLup * dimLdown;
3050 int inc = 1;
3051 daxpy_( &length, &prefactor, workmem2, &inc, Wblock, &inc );
3052
3053 }
3054
3055 }
3056 }
3057 }
3058 }
3059 }
3060 }
3061 }
3062
3063 }
3064
fill_tens_29_33(TensorT * denT,TensorF0 * tofill,TensorF0 * denF0,double * workmem) const3065 void CheMPS2::ThreeDM::fill_tens_29_33(TensorT * denT, TensorF0 * tofill, TensorF0 * denF0, double * workmem) const{
3066
3067 const int orb_i = denT->gIndex();
3068 assert( tofill->get_irrep() == denF0->get_irrep() );
3069 tofill->clear();
3070
3071 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3072 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3073 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3074
3075 const int ILxImxIn = Irreps::directProd( IL, denF0->get_irrep() );
3076 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3077 const int ILxImxInxIi = Irreps::directProd( ILxIi, denF0->get_irrep() );
3078
3079 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3080 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSL, ILxImxIn );
3081
3082 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3083
3084 {
3085 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3086 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSL, ILxImxIn );
3087
3088 if (( dimRup > 0 ) && ( dimRdown > 0 )){
3089
3090 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3091 double * Tdown = denT->gStorage( NL, TwoSL, ILxImxIn, NL+2, TwoSL, ILxImxIn );
3092 double * right = denF0->gStorage( NL+2, TwoSL, ILxImxIn, NL+2, TwoSL, IL );
3093 double * left = tofill->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, IL );
3094
3095 double factor = TwoSL + 1.0;
3096 char notrans = 'N';
3097 char trans = 'T';
3098 double zero = 0.0;
3099 double one = 1.0;
3100 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3101 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3102
3103 }
3104 }
3105
3106 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3107
3108 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3109 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxImxInxIi );
3110
3111 if (( dimRup > 0 ) && ( dimRdown > 0 )){
3112
3113 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3114 double * Tdown = denT->gStorage( NL, TwoSL, ILxImxIn, NL+1, TwoSR, ILxImxInxIi );
3115 double * right = denF0->gStorage( NL+1, TwoSR, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3116 double * left = tofill->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, IL );
3117
3118 double factor = 0.5 * ( TwoSR + 1 );
3119 char notrans = 'N';
3120 char trans = 'T';
3121 double zero = 0.0;
3122 double one = 1.0;
3123 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3124 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3125
3126 }
3127 }
3128 }
3129 }
3130 }
3131 }
3132
3133 }
3134
fill_tens_30_32(TensorT * denT,TensorF1 * tofill,TensorF1 * denF1,double * workmem) const3135 void CheMPS2::ThreeDM::fill_tens_30_32(TensorT * denT, TensorF1 * tofill, TensorF1 * denF1, double * workmem) const{
3136
3137 const int orb_i = denT->gIndex();
3138 assert( tofill->get_irrep() == denF1->get_irrep() );
3139 tofill->clear();
3140
3141 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3142 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3143 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3144
3145 const int ILxImxIn = Irreps::directProd( IL, denF1->get_irrep() );
3146
3147 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3148
3149 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3150 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxImxIn );
3151 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3152 int dimRdown = book->gCurrentDim( orb_i+1, NL+2, TwoSLprime, ILxImxIn );
3153
3154 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRup > 0 ) && ( dimRdown > 0 )){
3155
3156 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3157 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxImxIn, NL+2, TwoSLprime, ILxImxIn );
3158 double * right = denF1->gStorage( NL+2, TwoSLprime, ILxImxIn, NL+2, TwoSL, IL );
3159 double * left = tofill->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3160
3161 double factor = sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSL - TwoSLprime );
3162 char notrans = 'N';
3163 char trans = 'T';
3164 double zero = 0.0;
3165 double one = 1.0;
3166 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3167 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3168
3169 }
3170 }
3171 }
3172 }
3173 }
3174
3175 }
3176
fill_tens_36_42(TensorT * denT,TensorF1 * tofill,TensorF0 * denF0,double * workmem) const3177 void CheMPS2::ThreeDM::fill_tens_36_42(TensorT * denT, TensorF1 * tofill, TensorF0 * denF0, double * workmem) const{
3178
3179 const int orb_i = denT->gIndex();
3180 assert( tofill->get_irrep() == denF0->get_irrep() );
3181 tofill->clear();
3182
3183 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3184 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3185 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3186
3187 const int ILxImxIn = Irreps::directProd( IL, denF0->get_irrep() );
3188 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3189 const int ILxImxInxIi = Irreps::directProd( ILxIi, denF0->get_irrep() );
3190
3191 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3192
3193 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3194 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxImxIn );
3195
3196 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3197
3198 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3199
3200 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3201 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxImxInxIi );
3202
3203 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) == 1 )){
3204
3205 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3206 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxImxIn, NL+1, TwoSR, ILxImxInxIi );
3207 double * right = denF0->gStorage( NL+1, TwoSR, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3208 double * left = tofill->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3209
3210 double factor = 0.5 * sqrt( 6.0 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
3211 * Special::phase( TwoSR + TwoSLprime + 1 )
3212 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR );
3213 char notrans = 'N';
3214 char trans = 'T';
3215 double zero = 0.0;
3216 double one = 1.0;
3217 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3218 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3219
3220 }
3221 }
3222 }
3223 }
3224 }
3225 }
3226 }
3227
3228 }
3229
fill_tens_34_35_37_38(TensorT * denT,TensorF1 * fill34,TensorF0 * fill35,TensorF1 * fill37,TensorF1 * fill38,TensorF1 * denF1,double * workmem,double * workmem2) const3230 void CheMPS2::ThreeDM::fill_tens_34_35_37_38(TensorT * denT, TensorF1 * fill34, TensorF0 * fill35, TensorF1 * fill37, TensorF1 * fill38, TensorF1 * denF1, double * workmem, double * workmem2) const{
3231
3232 const int orb_i = denT->gIndex();
3233 fill34->clear();
3234 fill35->clear();
3235 fill37->clear();
3236 fill38->clear();
3237
3238 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3239 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3240 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3241
3242 const int ILxImxIn = Irreps::directProd( IL, denF1->get_irrep() );
3243 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3244 const int ILxImxInxIi = Irreps::directProd( ILxIi, denF1->get_irrep() );
3245
3246 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3247
3248 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3249 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxImxIn );
3250
3251 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3252
3253 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3254 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
3255
3256 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3257 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSRprime, ILxImxInxIi );
3258
3259 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSR - TwoSRprime ) <= 2 )){
3260
3261 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3262 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxImxIn, NL+1, TwoSRprime, ILxImxInxIi );
3263 double * right = denF1->gStorage( NL+1, TwoSRprime, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3264
3265 char notrans = 'N';
3266 char trans = 'T';
3267 double zero = 0.0;
3268 double one = 1.0;
3269 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &one, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3270 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &zero, workmem2, &dimLdown );
3271
3272 {
3273 double * left = fill34->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3274 double factor = 0.5 * ( TwoSRprime + 1 ) * sqrt( 1.0 * ( TwoSR + 1 ) * ( TwoSL + 1 ) )
3275 * Special::phase( TwoSL + TwoSR + 3 )
3276 * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRprime, TwoSLprime, 2 );
3277 int length = dimLup * dimLdown;
3278 int inc = 1;
3279 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3280 }
3281 if ( TwoSL == TwoSLprime ){
3282 double * left = fill35->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3283 double factor = 0.5 * ( TwoSRprime + 1 ) * sqrt( 6.0 * ( TwoSR + 1 ) )
3284 * Special::phase( TwoSL + TwoSRprime + 3 )
3285 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSRprime, TwoSL );
3286 int length = dimLup * dimLdown;
3287 int inc = 1;
3288 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3289 }
3290 {
3291 double * left = fill37->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3292 double factor = 3 * ( TwoSRprime + 1 ) * sqrt( 1.0 * ( TwoSR + 1 ) * ( TwoSL + 1 ) )
3293 * Special::phase( 2*TwoSL + TwoSR + TwoSRprime )
3294 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR )
3295 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSRprime, TwoSLprime );
3296 int length = dimLup * dimLdown;
3297 int inc = 1;
3298 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3299 }
3300 {
3301 double * left = fill38->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3302 double factor = 3 * ( TwoSRprime + 1 ) * sqrt( 1.0 * ( TwoSR + 1 ) * ( TwoSL + 1 ) )
3303 * Special::phase( 2*TwoSR + TwoSL + TwoSLprime )
3304 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSRprime )
3305 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSRprime, TwoSL );
3306 int length = dimLup * dimLdown;
3307 int inc = 1;
3308 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3309 }
3310 }
3311 }
3312 }
3313 }
3314 }
3315 }
3316 }
3317 }
3318
3319 }
3320
fill_tens_49_51(TensorT * denT,TensorF0 * tofill,TensorS0 * denS0,double * workmem) const3321 void CheMPS2::ThreeDM::fill_tens_49_51(TensorT * denT, TensorF0 * tofill, TensorS0 * denS0, double * workmem) const{
3322
3323 const int orb_i = denT->gIndex();
3324 assert( tofill->get_irrep() == denS0->get_irrep() );
3325 tofill->clear();
3326
3327 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3328 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3329 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3330
3331 const int ILxImxIn = Irreps::directProd( IL, denS0->get_irrep() );
3332
3333 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3334 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSL, ILxImxIn );
3335 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3336 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxImxIn );
3337
3338 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRup > 0 ) && ( dimRdown > 0 )){
3339
3340 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3341 double * Tdown = denT->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, ILxImxIn );
3342 double * right = denS0->gStorage( NL, TwoSL, ILxImxIn, NL+2, TwoSL, IL );
3343 double * left = tofill->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, IL );
3344
3345 double factor = - ( TwoSL + 1.0 );
3346 char notrans = 'N';
3347 char trans = 'T';
3348 double zero = 0.0;
3349 double one = 1.0;
3350 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3351 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3352
3353 }
3354 }
3355 }
3356 }
3357
3358 }
3359
fill_tens_50_52(TensorT * denT,TensorF1 * tofill,TensorS1 * denS1,double * workmem) const3360 void CheMPS2::ThreeDM::fill_tens_50_52(TensorT * denT, TensorF1 * tofill, TensorS1 * denS1, double * workmem) const{
3361
3362 const int orb_i = denT->gIndex();
3363 assert( tofill->get_irrep() == denS1->get_irrep() );
3364 tofill->clear();
3365
3366 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3367 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3368 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3369
3370 const int ILxImxIn = Irreps::directProd( IL, denS1->get_irrep() );
3371
3372 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3373
3374 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3375 int dimLdown = book->gCurrentDim( orb_i, NL, TwoSLprime, ILxImxIn );
3376 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3377 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSLprime, ILxImxIn );
3378
3379 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRup > 0 ) && ( dimRdown > 0 )){
3380
3381 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3382 double * Tdown = denT->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSLprime, ILxImxIn );
3383 double * right = denS1->gStorage( NL, TwoSLprime, ILxImxIn, NL+2, TwoSL, IL );
3384 double * left = tofill->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3385
3386 double factor = - ( TwoSL + 1.0 );
3387 char notrans = 'N';
3388 char trans = 'T';
3389 double zero = 0.0;
3390 double one = 1.0;
3391 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3392 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3393
3394 }
3395 }
3396 }
3397 }
3398 }
3399
3400 }
3401
fill_tens_22_24(TensorT * denT,TensorS0 * tofill,TensorS0 * denS0,double * workmem) const3402 void CheMPS2::ThreeDM::fill_tens_22_24(TensorT * denT, TensorS0 * tofill, TensorS0 * denS0, double * workmem) const{
3403
3404 const int orb_i = denT->gIndex();
3405 assert( tofill->get_irrep() == denS0->get_irrep() );
3406 tofill->clear();
3407
3408 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3409 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3410 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3411
3412 const int ILxImxIn = Irreps::directProd( IL, denS0->get_irrep() );
3413 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3414 const int ILxImxInxIi = Irreps::directProd( ILxIi, denS0->get_irrep() );
3415
3416 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3417 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSL, ILxImxIn );
3418
3419 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3420
3421 {
3422 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3423 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxImxIn );
3424
3425 if (( dimRup > 0 ) && ( dimRdown > 0 )){
3426
3427 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3428 double * Tdown = denT->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, ILxImxIn );
3429 double * right = denS0->gStorage( NL, TwoSL, ILxImxIn, NL+2, TwoSL, IL );
3430 double * left = tofill->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, IL );
3431
3432 double factor = TwoSL + 1.0;
3433 char notrans = 'N';
3434 char trans = 'T';
3435 double zero = 0.0;
3436 double one = 1.0;
3437 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3438 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3439
3440 }
3441 }
3442
3443 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3444
3445 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3446 int dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSR, ILxImxInxIi );
3447
3448 if (( dimRup > 0 ) && ( dimRdown > 0 )){
3449
3450 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3451 double * Tdown = denT->gStorage( NL-2, TwoSL, ILxImxIn, NL-1, TwoSR, ILxImxInxIi );
3452 double * right = denS0->gStorage( NL-1, TwoSR, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3453 double * left = tofill->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, IL );
3454
3455 double factor = 0.5 * ( TwoSR + 1 );
3456 char notrans = 'N';
3457 char trans = 'T';
3458 double zero = 0.0;
3459 double one = 1.0;
3460 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3461 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3462
3463 }
3464 }
3465 }
3466 }
3467 }
3468 }
3469
3470 }
3471
fill_tens_28(TensorT * denT,TensorS1 * tofill,TensorS0 * denS0,double * workmem) const3472 void CheMPS2::ThreeDM::fill_tens_28(TensorT * denT, TensorS1 * tofill, TensorS0 * denS0, double * workmem) const{
3473
3474 const int orb_i = denT->gIndex();
3475 assert( tofill->get_irrep() == denS0->get_irrep() );
3476 tofill->clear();
3477
3478 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3479 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3480 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3481
3482 const int ILxImxIn = Irreps::directProd( IL, denS0->get_irrep() );
3483 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3484 const int ILxImxInxIi = Irreps::directProd( ILxIi, denS0->get_irrep() );
3485
3486 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3487
3488 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3489 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSLprime, ILxImxIn );
3490
3491 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3492
3493 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3494
3495 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3496 int dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSR, ILxImxInxIi );
3497
3498 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSLprime - TwoSR ) == 1 )){
3499
3500 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3501 double * Tdown = denT->gStorage( NL-2, TwoSLprime, ILxImxIn, NL-1, TwoSR, ILxImxInxIi );
3502 double * right = denS0->gStorage( NL-1, TwoSR, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3503 double * left = tofill->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3504
3505 double factor = sqrt( 1.5 * ( TwoSL + 1 ) ) * ( TwoSR + 1 )
3506 * Special::phase( TwoSLprime + TwoSR + 1 )
3507 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSR );
3508 char notrans = 'N';
3509 char trans = 'T';
3510 double zero = 0.0;
3511 double one = 1.0;
3512 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3513 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3514
3515 }
3516 }
3517 }
3518 }
3519 }
3520 }
3521 }
3522
3523 }
3524
fill_tens_23(TensorT * denT,TensorS1 * tofill,TensorS1 * denS1,double * workmem) const3525 void CheMPS2::ThreeDM::fill_tens_23(TensorT * denT, TensorS1 * tofill, TensorS1 * denS1, double * workmem) const{
3526
3527 const int orb_i = denT->gIndex();
3528 assert( tofill->get_irrep() == denS1->get_irrep() );
3529 tofill->clear();
3530
3531 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3532 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3533 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3534
3535 const int ILxImxIn = Irreps::directProd( IL, denS1->get_irrep() );
3536
3537 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3538
3539 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3540 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSLprime, ILxImxIn );
3541 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3542 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSLprime, ILxImxIn );
3543
3544 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRup > 0 ) && ( dimRdown > 0 )){
3545
3546 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3547 double * Tdown = denT->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSLprime, ILxImxIn );
3548 double * right = denS1->gStorage( NL, TwoSLprime, ILxImxIn, NL+2, TwoSL, IL );
3549 double * left = tofill->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3550
3551 double factor = TwoSL + 1.0;
3552 char notrans = 'N';
3553 char trans = 'T';
3554 double zero = 0.0;
3555 double one = 1.0;
3556 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3557 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3558
3559 }
3560 }
3561 }
3562 }
3563 }
3564
3565 }
3566
fill_tens_25_26_27(TensorT * denT,TensorS1 * fill25,TensorS1 * fill26,TensorS0 * fill27,TensorS1 * denS1,double * workmem,double * workmem2) const3567 void CheMPS2::ThreeDM::fill_tens_25_26_27(TensorT * denT, TensorS1 * fill25, TensorS1 * fill26, TensorS0 * fill27, TensorS1 * denS1, double * workmem, double * workmem2) const{
3568
3569 const int orb_i = denT->gIndex();
3570 fill25->clear();
3571 fill26->clear();
3572 fill27->clear();
3573
3574 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3575 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3576 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3577
3578 const int ILxImxIn = Irreps::directProd( IL, denS1->get_irrep() );
3579 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3580 const int ILxImxInxIi = Irreps::directProd( ILxIi, denS1->get_irrep() );
3581
3582 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3583
3584 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3585 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSLprime, ILxImxIn );
3586
3587 if (( dimLup > 0 ) && ( dimLdown > 0 )){
3588
3589 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3590 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
3591
3592 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3593 int dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSRprime, ILxImxInxIi );
3594
3595 if (( dimRup > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSRprime - TwoSR ) <= 2 )){
3596
3597 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3598 double * Tdown = denT->gStorage( NL-2, TwoSLprime, ILxImxIn, NL-1, TwoSRprime, ILxImxInxIi );
3599 double * right = denS1->gStorage( NL-1, TwoSRprime, ILxImxInxIi, NL+1, TwoSR, ILxIi );
3600
3601 char notrans = 'N';
3602 char trans = 'T';
3603 double zero = 0.0;
3604 double one = 1.0;
3605 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &one, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3606 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &zero, workmem2, &dimLdown );
3607
3608 {
3609 double * left = fill25->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3610 double factor = ( TwoSR + 1 ) * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) )
3611 * Special::phase( TwoSL + TwoSRprime + 3 )
3612 * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRprime, TwoSLprime, 2 );
3613 int length = dimLup * dimLdown;
3614 int inc = 1;
3615 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3616 }
3617 {
3618 double * left = fill26->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3619 double factor = 3 * ( TwoSR + 1 ) * sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) )
3620 * Special::phase( TwoSR + TwoSRprime + TwoSL + TwoSLprime + 2 )
3621 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSRprime, TwoSL )
3622 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSLprime, TwoSRprime );
3623 int length = dimLup * dimLdown;
3624 int inc = 1;
3625 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3626 }
3627 if ( TwoSLprime == TwoSL ){
3628 double * left = fill27->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3629 double factor = ( TwoSR + 1 ) * sqrt( 1.5 * ( TwoSRprime + 1 ) )
3630 * Special::phase( TwoSL + TwoSR + 3 )
3631 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSRprime, TwoSL );
3632 int length = dimLup * dimLdown;
3633 int inc = 1;
3634 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3635 }
3636 }
3637 }
3638 }
3639 }
3640 }
3641 }
3642 }
3643 }
3644
3645 }
3646
fill_tens_45_47(TensorT * denT,TensorS0 * tofill,TensorF0 * denF0,double * workmem,const bool first) const3647 void CheMPS2::ThreeDM::fill_tens_45_47(TensorT * denT, TensorS0 * tofill, TensorF0 * denF0, double * workmem, const bool first) const{
3648
3649 const int orb_i = denT->gIndex();
3650 assert( tofill->get_irrep() == denF0->get_irrep() );
3651 tofill->clear();
3652
3653 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3654 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3655 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3656
3657 const int ILxImxIn = Irreps::directProd( IL, denF0->get_irrep() );
3658
3659 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3660 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSL, ILxImxIn );
3661 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
3662 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSL, ILxImxIn );
3663
3664 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( dimRup > 0 ) && ( dimRdown > 0 )){
3665
3666 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
3667 double * Tdown = denT->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, ILxImxIn );
3668 double * left = tofill->gStorage( NL-2, TwoSL, ILxImxIn, NL, TwoSL, IL );
3669
3670 double factor = -( TwoSL + 1.0 );
3671 char notrans = 'N';
3672 char trans = 'T';
3673 double zero = 0.0;
3674 double one = 1.0;
3675 if ( first ){
3676 double * right = denF0->gStorage( NL, TwoSL, ILxImxIn, NL, TwoSL, IL );
3677 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3678 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3679 } else {
3680 double * right = denF0->gStorage( NL, TwoSL, IL, NL, TwoSL, ILxImxIn );
3681 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRup, &zero, workmem, &dimLdown );
3682 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3683 }
3684 }
3685 }
3686 }
3687 }
3688
3689 }
3690
fill_tens_46_48(TensorT * denT,TensorS1 * tofill,TensorF1 * denF1,double * workmem,const bool first) const3691 void CheMPS2::ThreeDM::fill_tens_46_48(TensorT * denT, TensorS1 * tofill, TensorF1 * denF1, double * workmem, const bool first) const{
3692
3693 const int orb_i = denT->gIndex();
3694 assert( tofill->get_irrep() == denF1->get_irrep() );
3695 tofill->clear();
3696
3697 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3698 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3699 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3700
3701 const int ILxImxIn = Irreps::directProd( IL, denF1->get_irrep() );
3702
3703 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3704 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
3705
3706 if (( dimLup > 0 ) && ( dimRup > 0 )){
3707
3708 for ( int TwoSLprime = TwoSL-2; TwoSLprime <= TwoSL+2; TwoSLprime+=2 ){
3709
3710 int dimLdown = book->gCurrentDim( orb_i, NL-2, TwoSLprime, ILxImxIn );
3711 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSLprime, ILxImxIn );
3712
3713 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
3714
3715 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
3716 double * Tdown = denT->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSLprime, ILxImxIn );
3717 double * left = tofill->gStorage( NL-2, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3718
3719 char notrans = 'N';
3720 char trans = 'T';
3721 double zero = 0.0;
3722 double one = 1.0;
3723 if ( first ){
3724 double factor = sqrt( 1.0 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSL - TwoSLprime + 2 );
3725 double * right = denF1->gStorage( NL, TwoSLprime, ILxImxIn, NL, TwoSL, IL );
3726 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3727 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3728 } else {
3729 double factor = - ( TwoSL + 1.0 );
3730 double * right = denF1->gStorage( NL, TwoSL, IL, NL, TwoSLprime, ILxImxIn );
3731 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRup, &zero, workmem, &dimLdown );
3732 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3733 }
3734 }
3735 }
3736 }
3737 }
3738 }
3739 }
3740
3741 }
3742
fill_53_54(TensorT * denT,Tensor3RDM * tofill,TensorL * denL,double * workmem) const3743 void CheMPS2::ThreeDM::fill_53_54(TensorT * denT, Tensor3RDM * tofill, TensorL * denL, double * workmem) const{
3744
3745 const int orb_i = denT->gIndex();
3746 tofill->clear();
3747
3748 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3749 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3750 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3751
3752 const int ILxIm = Irreps::directProd( IL, denL->get_irrep() );
3753
3754 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3755 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
3756
3757 if (( dimLup > 0 ) && ( dimRup > 0 )){
3758
3759 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
3760
3761 int dimLdown = book->gCurrentDim( orb_i, NL-3, TwoSLprime, ILxIm );
3762 int dimRdown = book->gCurrentDim( orb_i+1, NL-1, TwoSLprime, ILxIm );
3763
3764 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
3765
3766 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
3767 double * Tdown = denT->gStorage( NL-3, TwoSLprime, ILxIm, NL-1, TwoSLprime, ILxIm );
3768 double * left = tofill->gStorage( NL-3, TwoSLprime, ILxIm, NL, TwoSL, IL );
3769 double * right = denL->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3770
3771 char notrans = 'N';
3772 char trans = 'T';
3773 double zero = 0.0;
3774 double one = 1.0;
3775 double factor = - sqrt( 0.5 ) * ( TwoSL + 1 );
3776
3777 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3778 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3779
3780 }
3781 }
3782 }
3783 }
3784 }
3785 }
3786
3787 }
3788
fill_55_to_60(TensorT * denT,Tensor3RDM * tofill,TensorL * denL,double * workmem) const3789 void CheMPS2::ThreeDM::fill_55_to_60(TensorT * denT, Tensor3RDM * tofill, TensorL * denL, double * workmem) const{
3790
3791 const int orb_i = denT->gIndex();
3792 tofill->clear();
3793
3794 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3795 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3796 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3797
3798 const int ILxIm = Irreps::directProd( IL, denL->get_irrep() );
3799
3800 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3801 int dimRup = book->gCurrentDim( orb_i+1, NL, TwoSL, IL );
3802
3803 if (( dimLup > 0 ) && ( dimRup > 0 )){
3804
3805 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
3806
3807 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIm );
3808 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIm );
3809
3810 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
3811
3812 double * Tup = denT->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
3813 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIm, NL+1, TwoSLprime, ILxIm );
3814 double * left = tofill->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3815 double * right = denL->gStorage( NL, TwoSL, IL, NL+1, TwoSLprime, ILxIm );
3816
3817 char notrans = 'N';
3818 char trans = 'T';
3819 double zero = 0.0;
3820 double one = 1.0;
3821 double factor = sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSLprime + 1 ) ) * Special::phase( TwoSL + 1 - TwoSLprime );
3822
3823 dgemm_( ¬rans, &trans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRup, &zero, workmem, &dimLdown );
3824 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3825
3826 }
3827 }
3828 }
3829 }
3830 }
3831 }
3832
3833 }
3834
fill_61(TensorT * denT,Tensor3RDM * tofill,TensorL * denL,double * workmem) const3835 void CheMPS2::ThreeDM::fill_61(TensorT * denT, Tensor3RDM * tofill, TensorL * denL, double * workmem) const{
3836
3837 const int orb_i = denT->gIndex();
3838 tofill->clear();
3839
3840 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3841 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3842 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3843
3844 const int ILxIm = Irreps::directProd( IL, denL->get_irrep() );
3845
3846 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3847 int dimRup = book->gCurrentDim( orb_i+1, NL+2, TwoSL, IL );
3848
3849 if (( dimLup > 0 ) && ( dimRup > 0 )){
3850
3851 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
3852
3853 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIm );
3854 int dimRdown = book->gCurrentDim( orb_i+1, NL+1, TwoSLprime, ILxIm );
3855
3856 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
3857
3858 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+2, TwoSL, IL );
3859 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIm, NL+1, TwoSLprime, ILxIm );
3860 double * left = tofill->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3861 double * right = denL->gStorage( NL+1, TwoSLprime, ILxIm, NL+2, TwoSL, IL );
3862
3863 char notrans = 'N';
3864 char trans = 'T';
3865 double zero = 0.0;
3866 double one = 1.0;
3867 double factor = sqrt( 0.5 ) * ( TwoSL + 1 );
3868
3869 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &factor, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3870 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &one, left, &dimLdown );
3871
3872 }
3873 }
3874 }
3875 }
3876 }
3877 }
3878
3879 }
3880
fill_63_65(TensorT * denT,Tensor3RDM * fill63,Tensor3RDM * fill65,Tensor3RDM * fill67,Tensor3RDM * fill68,Tensor3RDM * fill76,Tensor3RDM * fill77,TensorL * denL,double * workmem,double * workmem2) const3881 void CheMPS2::ThreeDM::fill_63_65(TensorT * denT, Tensor3RDM * fill63, Tensor3RDM * fill65, Tensor3RDM * fill67, Tensor3RDM * fill68, Tensor3RDM * fill76, Tensor3RDM * fill77, TensorL * denL, double * workmem, double * workmem2) const{
3882
3883 const int orb_i = denT->gIndex();
3884 fill63->clear();
3885 fill65->clear();
3886 fill67->clear();
3887 fill68->clear();
3888 fill76->clear();
3889 fill77->clear();
3890
3891 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3892 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3893 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3894
3895 const int ILxIm = Irreps::directProd( IL, denL->get_irrep() );
3896 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
3897 const int ILxImxIi = Irreps::directProd( ILxIi, denL->get_irrep() );
3898
3899 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
3900
3901 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
3902 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
3903
3904 if (( dimLup > 0 ) && ( dimRup > 0 )){
3905
3906 for ( int TwoSLprime = TwoSL-1; TwoSLprime <= TwoSL+1; TwoSLprime+=2 ){
3907 for ( int TwoSRprime = TwoSLprime-1; TwoSRprime <= TwoSLprime+1; TwoSRprime+=2 ){
3908
3909 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIm );
3910 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSRprime, ILxImxIi );
3911
3912 if (( dimLdown > 0 ) && ( dimRdown > 0 ) && ( abs( TwoSR - TwoSRprime ) == 1 )){
3913
3914 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
3915 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSRprime, ILxImxIi );
3916 double * right = denL->gStorage( NL, TwoSRprime, ILxImxIi, NL+1, TwoSR, ILxIi );
3917
3918 char notrans = 'N';
3919 char trans = 'T';
3920 double zero = 0.0;
3921 double one = 1.0;
3922 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &one, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
3923 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &zero, workmem2, &dimLdown );
3924
3925 {
3926 double factor = sqrt( 0.5 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) ) * ( TwoSR + 1 )
3927 * Special::phase( TwoSL + TwoSRprime + 2 )
3928 * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRprime, TwoSLprime, 1 );
3929 double * left = fill63->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3930 int length = dimLup * dimLdown;
3931 int inc = 1;
3932 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3933 }
3934 if ( TwoSRprime == TwoSL ){
3935 double factor = - sqrt( 0.5 ) * ( TwoSR + 1 );
3936 double * left = fill65->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3937 int length = dimLup * dimLdown;
3938 int inc = 1;
3939 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3940 }
3941 if ( TwoSR == TwoSLprime ){
3942 double factor = sqrt( 0.5 * ( TwoSRprime + 1 ) * ( TwoSL + 1 ) ) * Special::phase( TwoSL - TwoSRprime );
3943 double * left = fill67->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3944 int length = dimLup * dimLdown;
3945 int inc = 1;
3946 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3947 }
3948 {
3949 double factor = sqrt( 6.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) ) * ( TwoSR + 1 )
3950 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSLprime, TwoSRprime )
3951 * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSLprime, TwoSL );
3952 double * left = fill68->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3953 int length = dimLup * dimLdown;
3954 int inc = 1;
3955 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3956 }
3957 {
3958 double factor = sqrt( 6.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) ) * ( TwoSR + 1 )
3959 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSLprime )
3960 * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSR )
3961 * Special::phase( TwoSLprime + TwoSRprime + 2 - TwoSL - TwoSR );
3962 double * left = fill76->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3963 int length = dimLup * dimLdown;
3964 int inc = 1;
3965 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3966 }
3967 {
3968 double factor = - sqrt( 6.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) ) * ( TwoSR + 1 )
3969 * Wigner::wigner9j( 1, 1, 2, TwoSLprime, TwoSL, 1, TwoSRprime, TwoSR, 1 );
3970 double * left = fill77->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
3971 int length = dimLup * dimLdown;
3972 int inc = 1;
3973 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
3974 }
3975
3976 }
3977 }
3978 }
3979 }
3980 }
3981 }
3982 }
3983 }
3984
3985 }
3986
fill_69_78_79(TensorT * denT,Tensor3RDM * fill69,Tensor3RDM * fill78,Tensor3RDM * fill79,TensorL * denL,double * workmem,double * workmem2) const3987 void CheMPS2::ThreeDM::fill_69_78_79(TensorT * denT, Tensor3RDM * fill69, Tensor3RDM * fill78, Tensor3RDM * fill79, TensorL * denL, double * workmem, double * workmem2) const{
3988
3989 const int orb_i = denT->gIndex();
3990 fill69->clear();
3991 fill78->clear();
3992 fill79->clear();
3993
3994 for ( int NL = book->gNmin( orb_i ); NL <= book->gNmax( orb_i ); NL++ ){
3995 for ( int TwoSL = book->gTwoSmin( orb_i, NL ); TwoSL <= book->gTwoSmax( orb_i, NL ); TwoSL+=2 ){
3996 for ( int IL = 0; IL < book->getNumberOfIrreps(); IL++ ){
3997
3998 const int ILxIm = Irreps::directProd( IL, denL->get_irrep() );
3999 const int ILxIi = Irreps::directProd( IL, book->gIrrep( orb_i ) );
4000 const int ILxImxIi = Irreps::directProd( ILxIi, denL->get_irrep() );
4001
4002 for ( int TwoSR = TwoSL-1; TwoSR <= TwoSL+1; TwoSR+=2 ){
4003
4004 int dimLup = book->gCurrentDim( orb_i, NL, TwoSL, IL );
4005 int dimRup = book->gCurrentDim( orb_i+1, NL+1, TwoSR, ILxIi );
4006
4007 if (( dimLup > 0 ) && ( dimRup > 0 )){
4008
4009 for ( int TwoSRprime = TwoSR-1; TwoSRprime <= TwoSR+1; TwoSRprime+=2 ){
4010 for ( int TwoSLprime = TwoSRprime-1; TwoSLprime <= TwoSRprime+1; TwoSLprime+=2 ){
4011
4012 int dimLdown = book->gCurrentDim( orb_i, NL-1, TwoSLprime, ILxIm );
4013 int dimRdown = book->gCurrentDim( orb_i+1, NL, TwoSRprime, ILxImxIi );
4014
4015 if (( dimLdown > 0 ) && ( dimRdown > 0 )){
4016
4017 double * Tup = denT->gStorage( NL, TwoSL, IL, NL+1, TwoSR, ILxIi );
4018 double * Tdown = denT->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSRprime, ILxImxIi );
4019 double * right = denL->gStorage( NL, TwoSRprime, ILxImxIi, NL+1, TwoSR, ILxIi );
4020
4021 char notrans = 'N';
4022 char trans = 'T';
4023 double zero = 0.0;
4024 double one = 1.0;
4025 dgemm_( ¬rans, ¬rans, &dimLdown, &dimRup, &dimRdown, &one, Tdown, &dimLdown, right, &dimRdown, &zero, workmem, &dimLdown );
4026 dgemm_( ¬rans, &trans, &dimLdown, &dimLup, &dimRup, &one, workmem, &dimLdown, Tup, &dimLup, &zero, workmem2, &dimLdown );
4027
4028 const double prefactor = 2 * sqrt( 3.0 * ( TwoSL + 1 ) * ( TwoSRprime + 1 ) ) * ( TwoSR + 1 );
4029 {
4030 double factor = prefactor * Wigner::wigner6j( 1, 1, 2, TwoSR, TwoSLprime, TwoSRprime )
4031 * Wigner::wigner6j( 1, 2, 3, TwoSLprime, TwoSL, TwoSR );
4032 double * left = fill69->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
4033 int length = dimLup * dimLdown;
4034 int inc = 1;
4035 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
4036 }
4037 {
4038 double factor = prefactor * Wigner::wigner6j( 1, 1, 2, TwoSL, TwoSRprime, TwoSR )
4039 * Wigner::wigner6j( 1, 3, 2, TwoSL, TwoSRprime, TwoSLprime )
4040 * Special::phase( TwoSR + TwoSRprime + TwoSL + TwoSLprime + 2 );
4041 double * left = fill78->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
4042 int length = dimLup * dimLdown;
4043 int inc = 1;
4044 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
4045 }
4046 {
4047 double factor = - prefactor * Wigner::wigner9j( 1, 1, 2, TwoSLprime, TwoSL, 3, TwoSRprime, TwoSR, 1 );
4048 double * left = fill79->gStorage( NL-1, TwoSLprime, ILxIm, NL, TwoSL, IL );
4049 int length = dimLup * dimLdown;
4050 int inc = 1;
4051 daxpy_( &length, &factor, workmem2, &inc, left, &inc );
4052 }
4053
4054 }
4055 }
4056 }
4057 }
4058 }
4059 }
4060 }
4061 }
4062
4063 }
4064
4065
4066
4067