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