CheMPS2
TwoDM.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 <iostream>
23 #include <math.h>
24 #include <algorithm>
25 #include <assert.h>
26 
27 #include "TwoDM.h"
28 #include "Lapack.h"
29 #include "MyHDF5.h"
30 #include "Options.h"
31 #include "MPIchemps2.h"
32 #include "Wigner.h"
33 #include "Special.h"
34 
35 using std::max;
36 using std::cout;
37 using std::endl;
38 
39 CheMPS2::TwoDM::TwoDM(const SyBookkeeper * denBKIn, const Problem * ProbIn){
40 
41  denBK = denBKIn;
42  Prob = ProbIn;
43  L = denBK->gL();
44 
45  const long long size = ((long long) L ) * ((long long) L ) * ((long long) L ) * ((long long) L );
46  assert( INT_MAX >= size );
47  two_rdm_A = new double[ size ];
48  two_rdm_B = new double[ size ];
49 
50  //Clear the storage so that an allreduce can be performed in the end
51  for (int cnt = 0; cnt < size; cnt++){ two_rdm_A[ cnt ] = 0.0; }
52  for (int cnt = 0; cnt < size; cnt++){ two_rdm_B[ cnt ] = 0.0; }
53 
54 }
55 
56 CheMPS2::TwoDM::~TwoDM(){
57 
58  delete [] two_rdm_A;
59  delete [] two_rdm_B;
60 
61 }
62 
63 #ifdef CHEMPS2_MPI_COMPILATION
64 void CheMPS2::TwoDM::mpi_allreduce(){
65 
66  const int size = L*L*L*L; // Tested OK in creator TwoDM
67  double * temp = new double[ size ];
68  MPIchemps2::allreduce_array_double( two_rdm_A, temp, size ); for (int cnt = 0; cnt < size; cnt++){ two_rdm_A[ cnt ] = temp[ cnt ]; }
69  MPIchemps2::allreduce_array_double( two_rdm_B, temp, size ); for (int cnt = 0; cnt < size; cnt++){ two_rdm_B[ cnt ] = temp[ cnt ]; }
70  delete [] temp;
71 
72 }
73 #endif
74 
75 void CheMPS2::TwoDM::set_2rdm_A_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value){
76 
77  //Prob assumes you use DMRG orbs...
78  //Irrep sanity checks are performed in TwoDM::FillSite
79  two_rdm_A[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * cnt4 ) ) ] = value;
80  two_rdm_A[ cnt2 + L * ( cnt1 + L * ( cnt4 + L * cnt3 ) ) ] = value;
81  two_rdm_A[ cnt3 + L * ( cnt4 + L * ( cnt1 + L * cnt2 ) ) ] = value;
82  two_rdm_A[ cnt4 + L * ( cnt3 + L * ( cnt2 + L * cnt1 ) ) ] = value;
83 
84 }
85 
86 void CheMPS2::TwoDM::set_2rdm_B_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const double value){
87 
88  //Prob assumes you use DMRG orbs...
89  //Irrep sanity checks are performed in TwoDM::FillSite
90  two_rdm_B[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * cnt4 ) ) ] = value;
91  two_rdm_B[ cnt2 + L * ( cnt1 + L * ( cnt4 + L * cnt3 ) ) ] = value;
92  two_rdm_B[ cnt3 + L * ( cnt4 + L * ( cnt1 + L * cnt2 ) ) ] = value;
93  two_rdm_B[ cnt4 + L * ( cnt3 + L * ( cnt2 + L * cnt1 ) ) ] = value;
94 
95 }
96 
97 double CheMPS2::TwoDM::getTwoDMA_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const{
98 
99  //Prob assumes you use DMRG orbs...
100  const int irrep1 = Prob->gIrrep(cnt1);
101  const int irrep2 = Prob->gIrrep(cnt2);
102  const int irrep3 = Prob->gIrrep(cnt3);
103  const int irrep4 = Prob->gIrrep(cnt4);
104  if ( Irreps::directProd(irrep1, irrep2) == Irreps::directProd(irrep3, irrep4) ){
105  return two_rdm_A[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * cnt4 ) ) ];
106  }
107 
108  return 0.0;
109 
110 }
111 
112 double CheMPS2::TwoDM::getTwoDMB_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const{
113 
114  //Prob assumes you use DMRG orbs...
115  const int irrep1 = Prob->gIrrep(cnt1);
116  const int irrep2 = Prob->gIrrep(cnt2);
117  const int irrep3 = Prob->gIrrep(cnt3);
118  const int irrep4 = Prob->gIrrep(cnt4);
119  if ( Irreps::directProd(irrep1, irrep2) == Irreps::directProd(irrep3, irrep4) ){
120  return two_rdm_B[ cnt1 + L * ( cnt2 + L * ( cnt3 + L * cnt4 ) ) ];
121  }
122 
123  return 0.0;
124 
125 }
126 
127 double CheMPS2::TwoDM::get1RDM_DMRG(const int cnt1, const int cnt2) const{
128 
129  //Prob assumes you use DMRG orbs...
130  const int irrep1 = Prob->gIrrep(cnt1);
131  const int irrep2 = Prob->gIrrep(cnt2);
132  if ( irrep1 == irrep2 ){
133  double value = 0.0;
134  for ( int orbsum = 0; orbsum < L; orbsum++ ){ value += getTwoDMA_DMRG( cnt1, orbsum, cnt2, orbsum ); }
135  value = value / ( Prob->gN() - 1.0 );
136  return value;
137  }
138 
139  return 0.0;
140 
141 }
142 
143 double CheMPS2::TwoDM::spin_density_dmrg( const int cnt1, const int cnt2 ) const{
144 
145  //Prob assumes you use DMRG orbs...
146  const int irrep1 = Prob->gIrrep(cnt1);
147  const int irrep2 = Prob->gIrrep(cnt2);
148  if ( irrep1 == irrep2 ){
149  const int two_s = Prob->gTwoS();
150  if ( two_s > 0 ){
151  double value = ( 2 - Prob->gN() ) * get1RDM_DMRG( cnt1, cnt2 );
152  for ( int orb = 0; orb < Prob->gL(); orb++ ){
153  value -= ( getTwoDMA_DMRG( cnt1, orb, orb, cnt2 ) + getTwoDMB_DMRG( cnt1, orb, orb, cnt2 ) );
154  }
155  value = 1.5 * value / ( 0.5 * two_s + 1 );
156  return value;
157  }
158  }
159 
160  return 0.0;
161 
162 }
163 
164 double CheMPS2::TwoDM::getTwoDMA_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const{
165 
166  //Prob assumes you use DMRG orbs... f1 converts HAM orbs to DMRG orbs
167  if ( Prob->gReorder() ){
168  return getTwoDMA_DMRG( Prob->gf1(cnt1), Prob->gf1(cnt2), Prob->gf1(cnt3), Prob->gf1(cnt4) );
169  }
170  return getTwoDMA_DMRG( cnt1, cnt2, cnt3, cnt4 );
171 
172 }
173 
174 double CheMPS2::TwoDM::getTwoDMB_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const{
175 
176  //Prob assumes you use DMRG orbs... f1 converts HAM orbs to DMRG orbs
177  if ( Prob->gReorder() ){
178  return getTwoDMB_DMRG( Prob->gf1(cnt1), Prob->gf1(cnt2), Prob->gf1(cnt3), Prob->gf1(cnt4) );
179  }
180  return getTwoDMB_DMRG( cnt1, cnt2, cnt3, cnt4 );
181 
182 }
183 
184 double CheMPS2::TwoDM::get1RDM_HAM(const int cnt1, const int cnt2) const{
185 
186  //Prob assumes you use DMRG orbs... f1 converts HAM orbs to DMRG orbs
187  if ( Prob->gReorder() ){
188  return get1RDM_DMRG( Prob->gf1(cnt1), Prob->gf1(cnt2) );
189  }
190  return get1RDM_DMRG( cnt1, cnt2 );
191 
192 }
193 
194 double CheMPS2::TwoDM::spin_density_ham( const int cnt1, const int cnt2 ) const{
195 
196  //Prob assumes you use DMRG orbs... f1 converts HAM orbs to DMRG orbs
197  if ( Prob->gReorder() ){
198  return spin_density_dmrg( Prob->gf1(cnt1), Prob->gf1(cnt2) );
199  }
200  return spin_density_dmrg( cnt1, cnt2 );
201 
202 }
203 
204 double CheMPS2::TwoDM::trace() const{
205 
206  double val = 0.0;
207  for (int cnt1=0; cnt1<L; cnt1++){
208  for (int cnt2=0; cnt2<L; cnt2++){
209  val += getTwoDMA_DMRG(cnt1,cnt2,cnt1,cnt2);
210  }
211  }
212  return val;
213 
214 }
215 
216 double CheMPS2::TwoDM::energy() const{
217 
218  double val = 0.0;
219  for (int cnt1=0; cnt1<L; cnt1++){
220  for (int cnt2=0; cnt2<L; cnt2++){
221  for (int cnt3=0; cnt3<L; cnt3++){
222  for (int cnt4=0; cnt4<L; cnt4++){
223  val += getTwoDMA_DMRG(cnt1,cnt2,cnt3,cnt4) * Prob->gMxElement(cnt1,cnt2,cnt3,cnt4);
224  }
225  }
226  }
227  }
228  val *= 0.5;
229  return val + Prob->gEconst();
230 
231 }
232 
233 void CheMPS2::TwoDM::print_noon() const{
234 
235  int lwork = 3 * L;
236  double * OneRDM = new double[ L * L ];
237  double * work = new double[ lwork ];
238  double * eigs = new double[ L ];
239 
240  Irreps my_irreps( Prob->gSy() );
241 
242  for ( int irrep = 0; irrep < denBK->getNumberOfIrreps(); irrep++ ){
243 
244  int jump1 = 0;
245  for ( int orb1 = 0; orb1 < L; orb1++ ){
246  if ( Prob->gIrrep( orb1 ) == irrep ){
247  int jump2 = jump1;
248  for ( int orb2 = orb1; orb2 < L; orb2++ ){
249  if ( Prob->gIrrep( orb2 ) == irrep ){
250  const double value = get1RDM_DMRG( orb1, orb2 );
251  OneRDM[ jump1 + L * jump2 ] = value;
252  OneRDM[ jump2 + L * jump1 ] = value;
253  jump2 += 1;
254  }
255  }
256  jump1 += 1;
257  }
258  }
259 
260  if ( jump1 > 0 ){
261  char jobz = 'N'; // Eigenvalues only
262  char uplo = 'U';
263  int lda = L;
264  int info;
265  dsyev_(&jobz, &uplo, &jump1, OneRDM, &lda, eigs, work, &lwork, &info);
266  cout << " NOON of irrep " << my_irreps.getIrrepName( irrep ) << " = [ ";
267  for ( int cnt = 0; cnt < jump1 - 1; cnt++ ){ cout << eigs[ jump1 - 1 - cnt ] << " , "; } // Print from large to small
268  cout << eigs[ 0 ] << " ]." << endl;
269  }
270 
271  }
272  delete [] OneRDM;
273  delete [] work;
274  delete [] eigs;
275 
276 }
277 
278 void CheMPS2::TwoDM::save_HAM( const string filename ) const{
279 
280  // Create an array with the 2-RDM in the ORIGINAL HAM indices
281  const int total_size = L * L * L * L;
282  double * local_array = new double[ total_size ];
283  for ( int ham4 = 0; ham4 < L; ham4++ ){
284  for ( int ham3 = 0; ham3 < L; ham3++ ){
285  for ( int ham2 = 0; ham2 < L; ham2++ ){
286  for ( int ham1 = 0; ham1 < L; ham1++ ){
287  local_array[ ham1 + L * ( ham2 + L * ( ham3 + L * ham4 ) ) ] = getTwoDMA_HAM( ham1, ham2, ham3, ham4 );
288  }
289  }
290  }
291  }
292 
293  // (Re)create the HDF5 file with the 2-RDM
294  hid_t file_id = H5Fcreate( filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT );
295  hsize_t dimarray = total_size;
296  hid_t group_id = H5Gcreate( file_id, "2-RDM", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
297  hid_t dataspace_id = H5Screate_simple( 1, &dimarray, NULL );
298  hid_t dataset_id = H5Dcreate( group_id, "elements", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
299 
300  H5Dwrite( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, local_array );
301 
302  H5Dclose( dataset_id );
303  H5Sclose( dataspace_id );
304  H5Gclose( group_id );
305  H5Fclose( file_id );
306 
307  // Deallocate the array
308  delete [] local_array;
309 
310  cout << "Saved the 2-RDM to the file " << filename << endl;
311 
312 }
313 
314 void CheMPS2::TwoDM::save() const{
315 
316  hid_t file_id = H5Fcreate(CheMPS2::TWO_RDM_storagename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
317  hsize_t dimarray = L*L*L*L;
318  {
319  hid_t group_id = H5Gcreate(file_id, "two_rdm_A", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
320 
321  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
322  hid_t dataset_id = H5Dcreate(group_id, "elements", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
323  H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, two_rdm_A);
324 
325  H5Dclose(dataset_id);
326  H5Sclose(dataspace_id);
327 
328  H5Gclose(group_id);
329  }
330  {
331  hid_t group_id = H5Gcreate(file_id, "two_rdm_B", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
332 
333  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
334  hid_t dataset_id = H5Dcreate(group_id, "elements", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
335  H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, two_rdm_B);
336  H5Dclose(dataset_id);
337  H5Sclose(dataspace_id);
338 
339  H5Gclose(group_id);
340  }
341  H5Fclose(file_id);
342 
343 }
344 
345 void CheMPS2::TwoDM::read(){
346 
347  hid_t file_id = H5Fopen(CheMPS2::TWO_RDM_storagename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
348  {
349  hid_t group_id = H5Gopen(file_id, "two_rdm_A", H5P_DEFAULT);
350 
351  hid_t dataset_id = H5Dopen(group_id, "elements", H5P_DEFAULT);
352  H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, two_rdm_A);
353  H5Dclose(dataset_id);
354 
355  H5Gclose(group_id);
356  }
357  {
358  hid_t group_id = H5Gopen(file_id, "two_rdm_B", H5P_DEFAULT);
359 
360  hid_t dataset_id = H5Dopen(group_id, "elements", H5P_DEFAULT);
361  H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, two_rdm_B);
362  H5Dclose(dataset_id);
363 
364  H5Gclose(group_id);
365  }
366  H5Fclose(file_id);
367 
368  std::cout << "TwoDM::read : Everything loaded!" << std::endl;
369 
370 }
371 
372 void CheMPS2::TwoDM::write2DMAfile(const string filename) const{
373 
374  int * psi2molpro = new int[ denBK->getNumberOfIrreps() ];
375  Irreps my_irreps( Prob->gSy() );
376  my_irreps.symm_psi2molpro( psi2molpro );
377 
378  FILE * capturing;
379  capturing = fopen( filename.c_str(), "w" ); // "w" with fopen means truncate file
380  fprintf( capturing, " &2-RDM NORB= %d,NELEC= %d,MS2= %d,\n", L, Prob->gN(), Prob->gTwoS() );
381  fprintf( capturing, " ORBSYM=" );
382  for (int HamOrb=0; HamOrb<L; HamOrb++){
383  const int DMRGOrb = ( Prob->gReorder() ) ? Prob->gf1( HamOrb ) : HamOrb;
384  fprintf( capturing, "%d,", psi2molpro[ Prob->gIrrep( DMRGOrb ) ] );
385  }
386  fprintf( capturing, "\n ISYM=%d,\n /\n", psi2molpro[ Prob->gIrrep() ] );
387  delete [] psi2molpro;
388 
389  for (int ham_p=0; ham_p<L; ham_p++){
390  const int dmrg_p = (Prob->gReorder())?Prob->gf1(ham_p):ham_p;
391  for (int ham_q=0; ham_q<=ham_p; ham_q++){ // p >= q
392  const int dmrg_q = (Prob->gReorder())?Prob->gf1(ham_q):ham_q;
393  const int irrep_pq = Irreps::directProd( Prob->gIrrep(dmrg_p), Prob->gIrrep(dmrg_q) );
394  for (int ham_r=0; ham_r<=ham_p; ham_r++){ // p >= r
395  const int dmrg_r = (Prob->gReorder())?Prob->gf1(ham_r):ham_r;
396  for (int ham_s=0; ham_s<=ham_p; ham_s++){ // p >= s
397  const int dmrg_s = (Prob->gReorder())?Prob->gf1(ham_s):ham_s;
398  const int irrep_rs = Irreps::directProd( Prob->gIrrep(dmrg_r), Prob->gIrrep(dmrg_s) );
399  if ( irrep_pq == irrep_rs ){
400  const int num_equal = (( ham_q == ham_p ) ? 1 : 0 )
401  + (( ham_r == ham_p ) ? 1 : 0 )
402  + (( ham_s == ham_p ) ? 1 : 0 );
403  /* 1. p > q,r,s
404  2. p==q > r,s
405  3. p==r > q,s
406  4. p==s > q,r
407  5. p==q==r > s
408  6. p==q==s > r
409  7. p==r==s > q
410  8. p==r==s==q
411  While 2-4 are inequivalent ( num_equal == 1 ), 5-7 are equivalent ( num_equal == 2 ). Hence: */
412  if ( ( num_equal != 2 ) || ( ham_p > ham_s ) ){
413  const double value = getTwoDMA_DMRG(dmrg_p, dmrg_r, dmrg_q, dmrg_s);
414  fprintf( capturing, " % 23.16E %3d %3d %3d %3d\n", value, ham_p+1, ham_q+1, ham_r+1, ham_s+1 );
415  }
416  }
417  }
418  }
419  }
420  }
421 
422  // 1-RDM in Hamiltonian indices with p >= q
423  const double prefactor = 1.0 / ( Prob->gN() - 1.0 );
424  for (int ham_p=0; ham_p<L; ham_p++){
425  const int dmrg_p = (Prob->gReorder())?Prob->gf1(ham_p):ham_p;
426  for (int ham_q=0; ham_q<=ham_p; ham_q++){
427  const int dmrg_q = (Prob->gReorder())?Prob->gf1(ham_q):ham_q;
428  if ( Prob->gIrrep(dmrg_p) == Prob->gIrrep(dmrg_q) ){
429  double value = 0.0;
430  for ( int orbsum = 0; orbsum < L; orbsum++ ){ value += getTwoDMA_DMRG( dmrg_p, orbsum, dmrg_q, orbsum ); }
431  value *= prefactor;
432  fprintf( capturing, " % 23.16E %3d %3d %3d %3d\n", value, ham_p+1, ham_q+1, 0, 0 );
433  }
434  }
435  }
436 
437  // 0-RDM == Norm of the wavefunction?
438  fprintf( capturing, " % 23.16E %3d %3d %3d %3d", 1.0, 0, 0, 0, 0 );
439  fclose( capturing );
440  cout << "Created the file " << filename << "." << endl;
441 
442 }
443 
444 void CheMPS2::TwoDM::FillSite(TensorT * denT, TensorL *** Ltens, TensorF0 **** F0tens, TensorF1 **** F1tens, TensorS0 **** S0tens, TensorS1 **** S1tens){
445 
446  #ifdef CHEMPS2_MPI_COMPILATION
447  const int MPIRANK = MPIchemps2::mpi_rank();
448  #endif
449 
450  const int theindex = denT->gIndex();
451  const int DIM = max(denBK->gMaxDimAtBound(theindex), denBK->gMaxDimAtBound(theindex+1));
452 
453  #ifdef CHEMPS2_MPI_COMPILATION
454  if ( MPIRANK == MPI_CHEMPS2_MASTER )
455  #endif
456  {
457  //Diagram 1
458  const double d1 = doD1(denT);
459  set_2rdm_A_DMRG(theindex,theindex,theindex,theindex, 2*d1);
460  set_2rdm_B_DMRG(theindex,theindex,theindex,theindex,-2*d1);
461  }
462 
463  #pragma omp parallel
464  {
465 
466  double * workmem = new double[DIM*DIM];
467  double * workmem2 = new double[DIM*DIM];
468 
469  #pragma omp for schedule(static) nowait
470  for (int j_index=theindex+1; j_index<L; j_index++){
471  if (denBK->gIrrep(j_index) == denBK->gIrrep(theindex)){
472  #ifdef CHEMPS2_MPI_COMPILATION
473  if ( MPIRANK == MPIchemps2::owner_q( L, j_index ) ) //Everyone owns the L-tensors --> task division based on Q-tensor ownership
474  #endif
475  {
476  //Diagram 2
477  const double d2 = doD2(denT, Ltens[theindex][j_index-theindex-1], workmem);
478  set_2rdm_A_DMRG(theindex,j_index,theindex,theindex, 2*d2);
479  set_2rdm_B_DMRG(theindex,j_index,theindex,theindex,-2*d2);
480  }
481  }
482  }
483 
484  const int dimTriangle = L - theindex - 1;
485  const int upperboundTriangle = ( dimTriangle * ( dimTriangle + 1 ) ) / 2;
486  int result[ 2 ];
487  #pragma omp for schedule(static) nowait
488  for ( int global = 0; global < upperboundTriangle; global++ ){
489  Special::invert_triangle_two( global, result );
490  const int j_index = L - 1 - result[ 1 ];
491  const int k_index = j_index + result[ 0 ];
492  if ( denBK->gIrrep( j_index ) == denBK->gIrrep( k_index )){
493 
494  #ifdef CHEMPS2_MPI_COMPILATION
495  if ( MPIRANK == MPIchemps2::owner_absigma( j_index, k_index ) )
496  #endif
497  {
498  //Diagram 3
499  const double d3 = doD3(denT, S0tens[theindex][k_index-j_index][j_index-theindex-1], workmem);
500  set_2rdm_A_DMRG(theindex,theindex,j_index,k_index, 2*d3);
501  set_2rdm_B_DMRG(theindex,theindex,j_index,k_index,-2*d3);
502  }
503 
504  #ifdef CHEMPS2_MPI_COMPILATION
505  if ( MPIRANK == MPIchemps2::owner_cdf( L, j_index, k_index ) )
506  #endif
507  {
508  //Diagrams 4,5 & 6
509  const double d4 = doD4(denT, F0tens[theindex][k_index-j_index][j_index-theindex-1], workmem);
510  const double d5 = doD5(denT, F0tens[theindex][k_index-j_index][j_index-theindex-1], workmem);
511  const double d6 = doD6(denT, F1tens[theindex][k_index-j_index][j_index-theindex-1], workmem);
512  set_2rdm_A_DMRG(theindex,j_index,k_index,theindex, -2*d4 - 2*d5 - 3*d6);
513  set_2rdm_B_DMRG(theindex,j_index,k_index,theindex, -2*d4 - 2*d5 + d6);
514  set_2rdm_A_DMRG(theindex,j_index,theindex,k_index, 4*d4 + 4*d5);
515  set_2rdm_B_DMRG(theindex,j_index,theindex,k_index, 2*d6);
516  }
517  }
518  }
519 
520  #pragma omp for schedule(static) nowait
521  for (int g_index=0; g_index<theindex; g_index++){
522  if (denBK->gIrrep(g_index) == denBK->gIrrep(theindex)){
523  #ifdef CHEMPS2_MPI_COMPILATION
524  if ( MPIRANK == MPIchemps2::owner_q( L, g_index ) ) //Everyone owns the L-tensors --> task division based on Q-tensor ownership
525  #endif
526  {
527  //Diagram 7
528  const double d7 = doD7(denT, Ltens[theindex-1][theindex-g_index-1], workmem);
529  set_2rdm_A_DMRG(g_index,theindex,theindex,theindex, 2*d7);
530  set_2rdm_B_DMRG(g_index,theindex,theindex,theindex,-2*d7);
531  }
532  }
533  }
534 
535  const int globalsize8to12 = theindex * ( L - 1 - theindex );
536  #pragma omp for schedule(static) nowait
537  for (int gj_index=0; gj_index<globalsize8to12; gj_index++){
538  const int g_index = gj_index % theindex;
539  const int j_index = ( gj_index / theindex ) + theindex + 1;
540  const int I_g = denBK->gIrrep(g_index);
541  if (denBK->gIrrep(g_index) == denBK->gIrrep(j_index)){
542  #ifdef CHEMPS2_MPI_COMPILATION
543  if ( MPIRANK == MPIchemps2::owner_absigma( g_index, j_index ) ) //Everyone owns the L-tensors --> task division based on ABSigma-tensor ownership
544  #endif
545  {
546  //Diagrams 8,9,10 & 11
547  const double d8 = doD8(denT, Ltens[theindex-1][theindex-g_index-1], Ltens[theindex][j_index-theindex-1], workmem, workmem2, I_g);
548  double d9, d10, d11;
549  doD9andD10andD11(denT, Ltens[theindex-1][theindex-g_index-1], Ltens[theindex][j_index-theindex-1], workmem, workmem2, &d9, &d10, &d11, I_g);
550  set_2rdm_A_DMRG(g_index,theindex,j_index,theindex, -4*d8-d9);
551  set_2rdm_A_DMRG(g_index,theindex,theindex,j_index, 2*d8 + d11);
552  set_2rdm_B_DMRG(g_index,theindex,j_index,theindex, d9 - 2*d10);
553  set_2rdm_B_DMRG(g_index,theindex,theindex,j_index, 2*d8 + 2*d10 - d11);
554 
555  //Diagram 12
556  const double d12 = doD12(denT, Ltens[theindex-1][theindex-g_index-1], Ltens[theindex][j_index-theindex-1], workmem, workmem2, I_g);
557  set_2rdm_A_DMRG(g_index,j_index,theindex,theindex, 2*d12);
558  set_2rdm_B_DMRG(g_index,j_index,theindex,theindex,-2*d12);
559  }
560  }
561  }
562 
563  const int globalsize = theindex * upperboundTriangle;
564  #pragma omp for schedule(static) nowait
565  for ( int gjk_index = 0; gjk_index < globalsize; gjk_index++ ){
566  const int g_index = gjk_index % theindex;
567  const int global = gjk_index / theindex;
568  Special::invert_triangle_two( global, result );
569  const int j_index = L - 1 - result[ 1 ];
570  const int k_index = j_index + result[ 0 ];
571  const int I_g = denBK->gIrrep( g_index );
572  const int cnt1 = k_index - j_index;
573  const int cnt2 = j_index - theindex - 1;
574 
575  if (Irreps::directProd(I_g, denBK->gIrrep(theindex)) == Irreps::directProd(denBK->gIrrep(j_index), denBK->gIrrep(k_index))){
576  #ifdef CHEMPS2_MPI_COMPILATION
577  if ( MPIRANK == MPIchemps2::owner_absigma( j_index, k_index ) )
578  #endif
579  {
580  //Diagrams 13,14,15 & 16
581  const double d13 = doD13(denT, Ltens[theindex-1][theindex-g_index-1], S0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g);
582  const double d14 = doD14(denT, Ltens[theindex-1][theindex-g_index-1], S0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g);
583  double d15 = 0.0;
584  double d16 = 0.0;
585  if (k_index>j_index){
586  d15 = doD15(denT, Ltens[theindex-1][theindex-g_index-1], S1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g);
587  d16 = doD16(denT, Ltens[theindex-1][theindex-g_index-1], S1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g);
588  }
589  set_2rdm_A_DMRG(g_index,theindex,j_index,k_index, 2*d13 + 2*d14 + 3*d15 + 3*d16);
590  set_2rdm_A_DMRG(g_index,theindex,k_index,j_index, 2*d13 + 2*d14 - 3*d15 - 3*d16);
591  set_2rdm_B_DMRG(g_index,theindex,j_index,k_index,-2*d13 - 2*d14 + d15 + d16);
592  set_2rdm_B_DMRG(g_index,theindex,k_index,j_index,-2*d13 - 2*d14 - d15 - d16);
593  }
594 
595  #ifdef CHEMPS2_MPI_COMPILATION
596  if ( MPIRANK == MPIchemps2::owner_cdf( L, j_index, k_index ) )
597  #endif
598  {
599  //Diagrams 17,18,19 & 20
600  const double d17 = doD17orD21(denT, Ltens[theindex-1][theindex-g_index-1], F0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, true);
601  const double d18 = doD18orD22(denT, Ltens[theindex-1][theindex-g_index-1], F0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, true);
602  const double d19 = doD19orD23(denT, Ltens[theindex-1][theindex-g_index-1], F1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, true);
603  const double d20 = doD20orD24(denT, Ltens[theindex-1][theindex-g_index-1], F1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, true);
604  set_2rdm_A_DMRG(g_index,j_index,k_index,theindex, -2*d17 - 2*d18 - 3*d19 - 3*d20);
605  set_2rdm_A_DMRG(g_index,j_index,theindex,k_index, 4*d17 + 4*d18 );
606  set_2rdm_B_DMRG(g_index,j_index,k_index,theindex, -2*d17 - 2*d18 + d19 + d20);
607  set_2rdm_B_DMRG(g_index,j_index,theindex,k_index, 2*d19 + 2*d20);
608 
609  //Diagrams 21,22,23 & 24
610  const double d21 = doD17orD21(denT, Ltens[theindex-1][theindex-g_index-1], F0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, false);
611  const double d22 = doD18orD22(denT, Ltens[theindex-1][theindex-g_index-1], F0tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, false);
612  const double d23 = doD19orD23(denT, Ltens[theindex-1][theindex-g_index-1], F1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, false);
613  const double d24 = doD20orD24(denT, Ltens[theindex-1][theindex-g_index-1], F1tens[theindex][cnt1][cnt2], workmem, workmem2, I_g, false);
614  set_2rdm_A_DMRG(g_index,k_index,j_index,theindex, -2*d21 - 2*d22 - 3*d23 - 3*d24);
615  set_2rdm_A_DMRG(g_index,k_index,theindex,j_index, 4*d21 + 4*d22 );
616  set_2rdm_B_DMRG(g_index,k_index,j_index,theindex, -2*d21 - 2*d22 + d23 + d24);
617  set_2rdm_B_DMRG(g_index,k_index,theindex,j_index, 2*d23 + 2*d24);
618  }
619  }
620  }
621 
622  delete [] workmem;
623  delete [] workmem2;
624 
625  }
626 
627 }
628 
629 void CheMPS2::TwoDM::correct_higher_multiplicities(){
630 
631  if ( Prob->gTwoS() != 0 ){
632  double alpha = 1.0 / ( Prob->gTwoS() + 1.0 );
633  int length = L*L*L*L;
634  int inc = 1;
635  dscal_( &length, &alpha, two_rdm_A, &inc );
636  dscal_( &length, &alpha, two_rdm_B, &inc );
637  }
638 
639 }
640 
641 double CheMPS2::TwoDM::doD1(TensorT * denT){
642 
643  int theindex = denT->gIndex();
644 
645  double total = 0.0;
646 
647  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
648  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
649  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
650  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
651  int dimR = denBK->gCurrentDim(theindex+1,NL+2,TwoSL,IL);
652  if ((dimL>0) && (dimR>0)){
653 
654  double * Tblock = denT->gStorage(NL,TwoSL,IL,NL+2,TwoSL,IL);
655 
656  int length = dimL*dimR;
657  int inc = 1;
658  total += (TwoSL+1) * ddot_(&length, Tblock, &inc, Tblock, &inc);
659 
660  }
661  }
662  }
663  }
664 
665  return total;
666 
667 }
668 
669 double CheMPS2::TwoDM::doD2(TensorT * denT, TensorL * Lright, double * workmem){
670 
671  int theindex = denT->gIndex();
672 
673  double total = 0.0;
674 
675  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
676  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
677  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
678  for (int TwoSR = TwoSL-1; TwoSR<=TwoSL+1; TwoSR+=2){
679 
680  int IRup = Irreps::directProd(IL, denBK->gIrrep(theindex));
681 
682  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
683  int dimRdown = denBK->gCurrentDim(theindex+1,NL+2,TwoSL,IL);
684  int dimRup = denBK->gCurrentDim(theindex+1,NL+1,TwoSR,IRup);
685 
686  if ((dimL>0) && (dimRup>0) && (dimRdown>0)){
687 
688  double * Tdown = denT->gStorage(NL,TwoSL,IL,NL+2,TwoSL,IL );
689  double * Tup = denT->gStorage(NL,TwoSL,IL,NL+1,TwoSR,IRup);
690  double * Lblock = Lright->gStorage(NL+1,TwoSR,IRup,NL+2,TwoSL,IL);
691 
692  char trans = 'T';
693  char notrans = 'N';
694  double alpha = 1.0;
695  double beta = 0.0; //set
696  dgemm_(&notrans,&trans,&dimL,&dimRup,&dimRdown,&alpha,Tdown,&dimL,Lblock,&dimRup,&beta,workmem,&dimL);
697 
698  const double factor = Special::phase( TwoSL + 1 - TwoSR ) * 0.5 * sqrt((TwoSL+1)*(TwoSR+1.0));
699 
700  int length = dimL * dimRup;
701  int inc = 1;
702  total += factor * ddot_(&length, workmem, &inc, Tup, &inc);
703 
704  }
705  }
706  }
707  }
708  }
709 
710  return total;
711 
712 }
713 
714 double CheMPS2::TwoDM::doD3(TensorT * denT, TensorS0 * S0right, double * workmem){
715 
716  int theindex = denT->gIndex();
717 
718  double total = 0.0;
719 
720  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
721  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
722  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
723 
724  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
725  int dimRdown = denBK->gCurrentDim(theindex+1,NL, TwoSL,IL);
726  int dimRup = denBK->gCurrentDim(theindex+1,NL+2,TwoSL,IL);
727 
728  if ((dimL>0) && (dimRup>0) && (dimRdown>0)){
729 
730  double * Tdown = denT->gStorage(NL,TwoSL,IL,NL, TwoSL,IL);
731  double * Tup = denT->gStorage(NL,TwoSL,IL,NL+2,TwoSL,IL);
732  double * S0block = S0right->gStorage(NL, TwoSL, IL, NL+2, TwoSL, IL);
733 
734  char notrans = 'N';
735  double alpha = 1.0;
736  double beta = 0.0; //set
737  dgemm_(&notrans,&notrans,&dimL,&dimRup,&dimRdown,&alpha,Tdown,&dimL,S0block,&dimRdown,&beta,workmem,&dimL);
738 
739  double factor = sqrt(0.5) * (TwoSL+1);
740 
741  int length = dimL * dimRup;
742  int inc = 1;
743  total += factor * ddot_(&length, workmem, &inc, Tup, &inc);
744 
745  }
746  }
747  }
748  }
749 
750  return total;
751 
752 }
753 
754 double CheMPS2::TwoDM::doD4(TensorT * denT, TensorF0 * F0right, double * workmem){
755 
756  int theindex = denT->gIndex();
757 
758  double total = 0.0;
759 
760  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
761  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
762  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
763 
764  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
765  int dimR = denBK->gCurrentDim(theindex+1,NL+2,TwoSL,IL);
766 
767  if ((dimL>0) && (dimR>0)){
768 
769  double * Tblock = denT->gStorage(NL, TwoSL, IL, NL+2, TwoSL, IL);
770  double * F0block = F0right->gStorage(NL+2, TwoSL, IL, NL+2, TwoSL, IL);
771 
772  char notrans = 'N';
773  double alpha = 1.0;
774  double beta = 0.0; //set
775  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimR,&alpha,Tblock,&dimL,F0block,&dimR,&beta,workmem,&dimL);
776 
777  double factor = sqrt(0.5) * (TwoSL+1);
778 
779  int length = dimL * dimR;
780  int inc = 1;
781  total += factor * ddot_(&length, workmem, &inc, Tblock, &inc);
782 
783  }
784  }
785  }
786  }
787 
788  return total;
789 
790 }
791 
792 double CheMPS2::TwoDM::doD5(TensorT * denT, TensorF0 * F0right, double * workmem){
793 
794  int theindex = denT->gIndex();
795 
796  double total = 0.0;
797 
798  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
799  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
800  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
801  for (int TwoSR=TwoSL-1; TwoSR<=TwoSL+1; TwoSR+=2){
802 
803  int IR = Irreps::directProd(IL,denBK->gIrrep(theindex));
804  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
805  int dimR = denBK->gCurrentDim(theindex+1,NL+1,TwoSR,IR);
806 
807  if ((dimL>0) && (dimR>0)){
808 
809  double * Tblock = denT->gStorage(NL, TwoSL, IL, NL+1, TwoSR, IR);
810  double * F0block = F0right->gStorage(NL+1, TwoSR, IR, NL+1, TwoSR, IR);
811 
812  char notrans = 'N';
813  double alpha = 1.0;
814  double beta = 0.0; //set
815  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimR,&alpha,Tblock,&dimL,F0block,&dimR,&beta,workmem,&dimL);
816 
817  double factor = 0.5 * sqrt(0.5) * (TwoSR+1);
818 
819  int length = dimL * dimR;
820  int inc = 1;
821  total += factor * ddot_(&length, workmem, &inc, Tblock, &inc);
822 
823  }
824  }
825  }
826  }
827  }
828 
829  return total;
830 
831 }
832 
833 double CheMPS2::TwoDM::doD6(TensorT * denT, TensorF1 * F1right, double * workmem){
834 
835  int theindex = denT->gIndex();
836 
837  double total = 0.0;
838 
839  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
840  for (int TwoSL = denBK->gTwoSmin(theindex,NL); TwoSL<= denBK->gTwoSmax(theindex,NL); TwoSL+=2){
841  for (int IL = 0; IL<denBK->getNumberOfIrreps(); IL++){
842  for (int TwoSRup=TwoSL-1; TwoSRup<=TwoSL+1; TwoSRup+=2){
843  for (int TwoSRdown=TwoSL-1; TwoSRdown<=TwoSL+1; TwoSRdown+=2){
844 
845  int IR = Irreps::directProd(IL,denBK->gIrrep(theindex));
846  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL, IL);
847  int dimRup = denBK->gCurrentDim(theindex+1,NL+1,TwoSRup, IR);
848  int dimRdown = denBK->gCurrentDim(theindex+1,NL+1,TwoSRdown,IR);
849 
850  if ((dimL>0) && (dimRup>0) && (dimRdown>0)){
851 
852  double * Tup = denT->gStorage(NL, TwoSL, IL, NL+1, TwoSRup, IR);
853  double * Tdown = denT->gStorage(NL, TwoSL, IL, NL+1, TwoSRdown, IR);
854  double * F1block = F1right->gStorage(NL+1, TwoSRdown, IR, NL+1, TwoSRup, IR);
855 
856  char notrans = 'N';
857  double alpha = 1.0;
858  double beta = 0.0; //set
859  dgemm_(&notrans,&notrans,&dimL,&dimRup,&dimRdown,&alpha,Tdown,&dimL,F1block,&dimRdown,&beta,workmem,&dimL);
860 
861  const double factor = sqrt((TwoSRup+1)/3.0) * ( TwoSRdown + 1 )
862  * Special::phase( TwoSL + TwoSRdown - 1 )
863  * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSRdown, TwoSL );
864 
865  int length = dimL * dimRup;
866  int inc = 1;
867  total += factor * ddot_(&length, workmem, &inc, Tup, &inc);
868 
869  }
870  }
871  }
872  }
873  }
874  }
875 
876  return total;
877 
878 }
879 
880 double CheMPS2::TwoDM::doD7(TensorT * denT, TensorL * Lleft, double * workmem){
881 
882  int theindex = denT->gIndex();
883 
884  double total = 0.0;
885 
886  for (int NLup = denBK->gNmin(theindex); NLup<=denBK->gNmax(theindex); NLup++){
887  for (int TwoSLup = denBK->gTwoSmin(theindex,NLup); TwoSLup<= denBK->gTwoSmax(theindex,NLup); TwoSLup+=2){
888  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
889 
890  int dimLup = denBK->gCurrentDim(theindex, NLup, TwoSLup, ILup);
891 
892  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
893 
894  int IR = Irreps::directProd(ILup,denBK->gIrrep(theindex));
895  int dimLdown = denBK->gCurrentDim(theindex, NLup-1, TwoSLdown, IR);
896  int dimR = denBK->gCurrentDim(theindex+1, NLup+1, TwoSLdown, IR);
897 
898  if ((dimLup>0) && (dimLdown>0) && (dimR>0)){
899 
900  double * Tup = denT->gStorage(NLup, TwoSLup, ILup, NLup+1, TwoSLdown, IR);
901  double * Tdown = denT->gStorage(NLup-1, TwoSLdown, IR, NLup+1, TwoSLdown, IR);
902  double * Lblock = Lleft->gStorage(NLup-1, TwoSLdown, IR, NLup, TwoSLup, ILup);
903 
904  char trans = 'T';
905  char notrans = 'N';
906  double alpha = 1.0;
907  double beta = 0.0; //set
908  dgemm_(&trans,&notrans,&dimLup,&dimR,&dimLdown,&alpha,Lblock,&dimLdown,Tdown,&dimLdown,&beta,workmem,&dimLup);
909 
910  const double factor = 0.5 * sqrt((TwoSLdown+1)*(TwoSLup+1.0)) * Special::phase( TwoSLup - TwoSLdown + 3 );
911 
912  int length = dimLup * dimR;
913  int inc = 1;
914  total += factor * ddot_(&length, workmem, &inc, Tup, &inc);
915 
916  }
917  }
918  }
919  }
920  }
921 
922  return total;
923 
924 }
925 
926 double CheMPS2::TwoDM::doD8(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g){
927 
928  int theindex = denT->gIndex();
929 
930  double total = 0.0;
931 
932  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
933  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
934  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
935 
936  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
937  int dimRup = denBK->gCurrentDim(theindex+1, NL+2, TwoSLup, ILup);
938 
939  if ((dimLup>0) && (dimRup>0)){
940 
941  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
942 
943  int Idown = Irreps::directProd(ILup,Irrep_g);
944 
945  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
946  int dimRdown = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, Idown);
947 
948  if ((dimLdown>0) && (dimRdown>0)){
949 
950  double * Tup = denT->gStorage(NL, TwoSLup, ILup, NL+2, TwoSLup, ILup);
951  double * Tdown = denT->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSLdown, Idown);
952  double * LleftBlk = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
953  double * LrightBlk = Lright->gStorage(NL+1, TwoSLdown, Idown, NL+2, TwoSLup, ILup);
954 
955  char trans = 'T';
956  char notrans = 'N';
957  double alpha = 1.0;
958  double beta = 0.0; //set
959  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,LleftBlk,&dimLdown,Tdown,&dimLdown,&beta,workmem,&dimLup);
960 
961  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,LrightBlk,&dimRdown,&beta,workmem2,&dimLup);
962 
963  double factor = -0.5 * (TwoSLup+1);
964 
965  int length = dimLup * dimRup;
966  int inc = 1;
967  total += factor * ddot_(&length, workmem2, &inc, Tup, &inc);
968 
969  }
970  }
971  }
972  }
973  }
974  }
975 
976  return total;
977 
978 }
979 
980 void CheMPS2::TwoDM::doD9andD10andD11(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, double * d9, double * d10, double * d11, int Irrep_g){
981 
982  d9[0] = 0.0;
983  d10[0] = 0.0;
984  d11[0] = 0.0;
985 
986  int theindex = denT->gIndex();
987 
988  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
989  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
990  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
991 
992  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
993  if (dimLup>0){
994 
995  int IRup = Irreps::directProd(ILup, denBK->gIrrep(theindex));
996  int ILdown = Irreps::directProd(ILup, Irrep_g);
997  int IRdown = Irreps::directProd(ILdown, denBK->gIrrep(theindex));
998 
999  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1000  for (int TwoSRup=TwoSLup-1; TwoSRup<=TwoSLup+1; TwoSRup+=2){
1001  for (int TwoSRdown=TwoSRup-1; TwoSRdown<=TwoSRup+1; TwoSRdown+=2){
1002  if ((TwoSLdown>=0) && (TwoSRup>=0) && (TwoSRdown>=0) && (abs(TwoSLdown - TwoSRdown)<=1)){
1003 
1004  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1005  int dimRup = denBK->gCurrentDim(theindex+1, NL+1, TwoSRup, IRup);
1006  int dimRdown = denBK->gCurrentDim(theindex+1, NL, TwoSRdown, IRdown);
1007 
1008  if ((dimLdown>0) && (dimRup>0) && (dimRdown>0)){
1009 
1010  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+1, TwoSRup, IRup);
1011  double * T_down = denT->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSRdown, IRdown);
1012  double * LleftBlk = Lleft->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, ILup);
1013  double * LrightBlk = Lright->gStorage(NL, TwoSRdown, IRdown, NL+1, TwoSRup, IRup);
1014 
1015  char trans = 'T';
1016  char notrans = 'N';
1017  double alpha = 1.0;
1018  double beta = 0.0; //SET
1019 
1020  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,LleftBlk,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1021 
1022  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,LrightBlk,&dimRdown,&beta,workmem2,&dimLup);
1023 
1024  int length = dimLup * dimRup;
1025  int inc = 1;
1026  double value = ddot_(&length, workmem2, &inc, T_up, &inc);
1027 
1028  const double fact1 = Special::phase( TwoSLup + TwoSRdown + 2 )
1029  * ( TwoSRup + 1 ) * sqrt( ( TwoSRdown + 1 ) * ( TwoSLup + 1.0 ) )
1030  * Wigner::wigner6j( TwoSRup, 1, TwoSLup, TwoSLdown, 1, TwoSRdown );
1031  const double fact2 = 2 * ( TwoSRup + 1 ) * sqrt( ( TwoSRdown + 1 ) * ( TwoSLup + 1.0 ) )
1032  * Wigner::wigner6j( TwoSRup, TwoSLdown, 2, 1, 1, TwoSLup )
1033  * Wigner::wigner6j( TwoSRup, TwoSLdown, 2, 1, 1, TwoSRdown );
1034  const int fact3 = (( TwoSRdown == TwoSLup ) ? ( TwoSRup + 1 ) : 0 );
1035 
1036  d9[0] += fact1 * value;
1037  d10[0] += fact2 * value;
1038  d11[0] += fact3 * value;
1039 
1040  }
1041  }
1042  }
1043  }
1044  }
1045  }
1046  }
1047  }
1048  }
1049 
1050 }
1051 
1052 double CheMPS2::TwoDM::doD12(TensorT * denT, TensorL * Lleft, TensorL * Lright, double * workmem, double * workmem2, int Irrep_g){
1053 
1054  double d12 = 0.0;
1055 
1056  int theindex = denT->gIndex();
1057 
1058  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1059  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1060  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1061 
1062  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1063  int dimRup = denBK->gCurrentDim(theindex+1, NL, TwoSLup, ILup);
1064  if ((dimLup>0) && (dimRup>0)){
1065 
1066  int Idown = Irreps::directProd(ILup, Irrep_g);
1067 
1068  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1069 
1070  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
1071  int dimRdown = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, Idown);
1072 
1073  if ((dimLdown>0) && (dimRdown>0)){
1074 
1075  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL, TwoSLup, ILup);
1076  double * T_down = denT->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSLdown, Idown);
1077  double * LleftBlk = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
1078  double * LrightBlk = Lright->gStorage(NL, TwoSLup, ILup, NL+1, TwoSLdown, Idown);
1079 
1080  char trans = 'T';
1081  char notrans = 'N';
1082  double alpha = 1.0;
1083  double beta = 0.0; //SET
1084 
1085  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,LleftBlk,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1086 
1087  dgemm_(&notrans,&trans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,LrightBlk,&dimRup,&beta,workmem2,&dimLup);
1088 
1089  const double factor = Special::phase( TwoSLdown + 1 - TwoSLup )
1090  * 0.5 * sqrt( ( TwoSLup + 1 ) * ( TwoSLdown + 1.0 ) );
1091 
1092  int length = dimLup * dimRup;
1093  int inc = 1;
1094  d12 += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1095 
1096  }
1097  }
1098  }
1099  }
1100  }
1101  }
1102 
1103  return d12;
1104 
1105 }
1106 
1107 double CheMPS2::TwoDM::doD13(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g){
1108 
1109  double d13 = 0.0;
1110 
1111  int theindex = denT->gIndex();
1112 
1113  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1114  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1115  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1116 
1117  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1118  int dimRup = denBK->gCurrentDim(theindex+1, NL+2, TwoSLup, ILup);
1119 
1120  if ((dimLup>0) && (dimRup>0)){
1121 
1122  int ILdown = Irreps::directProd(ILup, Irrep_g);
1123  int IRdown = Irreps::directProd(ILdown, denBK->gIrrep(theindex));
1124 
1125  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1126 
1127  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1128  int dimRdown = denBK->gCurrentDim(theindex+1, NL, TwoSLup, IRdown);
1129 
1130  if ((dimLdown>0) && (dimRdown>0)){
1131 
1132  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+2, TwoSLup, ILup);
1133  double * T_down = denT->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, IRdown);
1134  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, ILup);
1135  double * S0block = S0right->gStorage(NL, TwoSLup, IRdown, NL+2, TwoSLup, ILup);
1136 
1137  char trans = 'T';
1138  char notrans = 'N';
1139  double alpha = 1.0;
1140  double beta = 0.0; //SET
1141 
1142  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1143 
1144  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,S0block,&dimRdown,&beta,workmem2,&dimLup);
1145 
1146  double factor = -0.5 * sqrt(0.5) * (TwoSLup+1);
1147 
1148  int length = dimLup * dimRup;
1149  int inc = 1;
1150  d13 += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1151 
1152  }
1153  }
1154  }
1155  }
1156  }
1157  }
1158 
1159  return d13;
1160 
1161 }
1162 
1163 double CheMPS2::TwoDM::doD14(TensorT * denT, TensorL * Lleft, TensorS0 * S0right, double * workmem, double * workmem2, int Irrep_g){
1164 
1165  double d14 = 0.0;
1166 
1167  int theindex = denT->gIndex();
1168 
1169  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1170  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1171  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1172 
1173  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1174 
1175  if (dimLup>0){
1176 
1177  int Idown = Irreps::directProd(ILup, Irrep_g);
1178  int IRup = Irreps::directProd(ILup, denBK->gIrrep(theindex));
1179 
1180  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1181 
1182  int dimRup = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, IRup);
1183  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
1184  int dimRdown = denBK->gCurrentDim(theindex+1, NL-1, TwoSLdown, Idown);
1185 
1186  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1187 
1188  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+1, TwoSLdown, IRup);
1189  double * T_down = denT->gStorage(NL-1, TwoSLdown, Idown, NL-1, TwoSLdown, Idown);
1190  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
1191  double * S0block = S0right->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSLdown, IRup);
1192 
1193  char trans = 'T';
1194  char notrans = 'N';
1195  double alpha = 1.0;
1196  double beta = 0.0; //SET
1197 
1198  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1199 
1200  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,S0block,&dimRdown,&beta,workmem2,&dimLup);
1201 
1202  const double factor = Special::phase( TwoSLdown + 1 - TwoSLup ) * 0.5 * sqrt(0.5 * (TwoSLup+1) * (TwoSLdown+1));
1203 
1204  int length = dimLup * dimRup;
1205  int inc = 1;
1206  d14 += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1207 
1208  }
1209  }
1210  }
1211  }
1212  }
1213  }
1214 
1215  return d14;
1216 
1217 }
1218 
1219 double CheMPS2::TwoDM::doD15(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g){
1220 
1221  double d15 = 0.0;
1222 
1223  int theindex = denT->gIndex();
1224 
1225  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1226  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1227  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1228 
1229  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1230  int dimRup = denBK->gCurrentDim(theindex+1, NL+2, TwoSLup, ILup);
1231 
1232  if ((dimLup>0) && (dimRup>0)){
1233 
1234  int ILdown = Irreps::directProd(ILup, Irrep_g);
1235  int IRdown = Irreps::directProd(ILdown, denBK->gIrrep(theindex));
1236 
1237  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1238  for (int TwoSRdown=TwoSLdown-1; TwoSRdown<=TwoSLdown+1; TwoSRdown+=2){
1239 
1240  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1241  int dimRdown = denBK->gCurrentDim(theindex+1, NL, TwoSRdown, IRdown);
1242 
1243  if ((dimLdown>0) && (dimRdown>0)){
1244 
1245  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+2, TwoSLup, ILup);
1246  double * T_down = denT->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSRdown, IRdown);
1247  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, ILup);
1248  double * S1block = S1right->gStorage(NL, TwoSRdown, IRdown, NL+2, TwoSLup, ILup);
1249 
1250  char trans = 'T';
1251  char notrans = 'N';
1252  double alpha = 1.0;
1253  double beta = 0.0; //SET
1254 
1255  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1256 
1257  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,S1block,&dimRdown,&beta,workmem2,&dimLup);
1258 
1259  const double factor = Special::phase( TwoSLdown + TwoSLup + 1 )
1260  * ( TwoSLup + 1 ) * sqrt((TwoSRdown+1)/3.0) * Wigner::wigner6j( 1, 1, 2, TwoSLup, TwoSRdown, TwoSLdown );
1261 
1262  int length = dimLup * dimRup;
1263  int inc = 1;
1264  d15 += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1265 
1266  }
1267  }
1268  }
1269  }
1270  }
1271  }
1272  }
1273 
1274  return d15;
1275 
1276 }
1277 
1278 double CheMPS2::TwoDM::doD16(TensorT * denT, TensorL * Lleft, TensorS1 * S1right, double * workmem, double * workmem2, int Irrep_g){
1279 
1280  double d16 = 0.0;
1281 
1282  int theindex = denT->gIndex();
1283 
1284  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1285  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1286  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1287 
1288  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1289 
1290  if (dimLup>0){
1291 
1292  int Idown = Irreps::directProd(ILup, Irrep_g);
1293  int IRup = Irreps::directProd(ILup, denBK->gIrrep(theindex));
1294 
1295  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1296  for (int TwoSRup=TwoSLup-1; TwoSRup<=TwoSLup+1; TwoSRup+=2){
1297 
1298  int dimRup = denBK->gCurrentDim(theindex+1, NL+1, TwoSRup, IRup);
1299  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
1300  int dimRdown = denBK->gCurrentDim(theindex+1, NL-1, TwoSLdown, Idown);
1301 
1302  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1303 
1304  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+1, TwoSRup, IRup);
1305  double * T_down = denT->gStorage(NL-1, TwoSLdown, Idown, NL-1, TwoSLdown, Idown);
1306  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
1307  double * S1block = S1right->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSRup, IRup);
1308 
1309  char trans = 'T';
1310  char notrans = 'N';
1311  double alpha = 1.0;
1312  double beta = 0.0; //SET
1313 
1314  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1315 
1316  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,S1block,&dimRdown,&beta,workmem2,&dimLup);
1317 
1318  const double factor = Special::phase( TwoSRup + TwoSLdown + 2 ) * (TwoSRup+1) * sqrt((TwoSLup+1)/3.0)
1319  * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSLdown, TwoSLup );
1320 
1321  int length = dimLup * dimRup;
1322  int inc = 1;
1323  d16 += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1324 
1325  }
1326  }
1327  }
1328  }
1329  }
1330  }
1331  }
1332 
1333  return d16;
1334 
1335 }
1336 
1337 double CheMPS2::TwoDM::doD17orD21(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD17){
1338 
1339  double total = 0.0;
1340 
1341  int theindex = denT->gIndex();
1342 
1343  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1344  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1345  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1346 
1347  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1348 
1349  if (dimLup>0){
1350 
1351  int ILdown = Irreps::directProd(ILup, Irrep_g);
1352  int IRdown = Irreps::directProd(ILdown, denBK->gIrrep(theindex));
1353 
1354  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1355 
1356  int dimRup = denBK->gCurrentDim(theindex+1, NL, TwoSLup, ILup);
1357  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1358  int dimRdown = denBK->gCurrentDim(theindex+1, NL, TwoSLup, IRdown);
1359 
1360  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1361 
1362  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL, TwoSLup, ILup);
1363  double * T_down = denT->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, IRdown);
1364  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, ILup);
1365  double * F0block = (shouldIdoD17) ? F0right->gStorage(NL, TwoSLup, IRdown, NL, TwoSLup, ILup)
1366  : F0right->gStorage(NL, TwoSLup, ILup, NL, TwoSLup, IRdown) ;
1367 
1368  char trans = 'T';
1369  char notrans = 'N';
1370  char var = (shouldIdoD17) ? notrans : trans;
1371  double alpha = 1.0;
1372  double beta = 0.0; //SET
1373  int dimvar = (shouldIdoD17) ? dimRdown : dimRup ;
1374 
1375  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1376 
1377  dgemm_(&notrans,&var,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,F0block,&dimvar,&beta,workmem2,&dimLup);
1378 
1379  int length = dimLup * dimRup;
1380  int inc = 1;
1381  total += sqrt(0.5) * 0.5 * (TwoSLup+1) * ddot_(&length, workmem2, &inc, T_up, &inc);
1382 
1383  }
1384  }
1385  }
1386  }
1387  }
1388  }
1389 
1390  return total;
1391 
1392 }
1393 
1394 double CheMPS2::TwoDM::doD18orD22(TensorT * denT, TensorL * Lleft, TensorF0 * F0right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD18){
1395 
1396  double total = 0.0;
1397 
1398  int theindex = denT->gIndex();
1399 
1400  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1401  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1402  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1403 
1404  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1405 
1406  if (dimLup>0){
1407 
1408  int Idown = Irreps::directProd(ILup, Irrep_g);
1409  int IRup = Irreps::directProd(ILup, denBK->gIrrep(theindex));
1410 
1411  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1412 
1413  int dimRup = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, IRup);
1414  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
1415  int dimRdown = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, Idown);
1416 
1417  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1418 
1419  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+1, TwoSLdown, IRup);
1420  double * T_down = denT->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSLdown, Idown);
1421  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
1422  double * F0block = (shouldIdoD18) ? F0right->gStorage(NL+1, TwoSLdown, Idown, NL+1, TwoSLdown, IRup)
1423  : F0right->gStorage(NL+1, TwoSLdown, IRup, NL+1, TwoSLdown, Idown) ;
1424 
1425  char trans = 'T';
1426  char notrans = 'N';
1427  char var = (shouldIdoD18) ? notrans : trans;
1428  double alpha = 1.0;
1429  double beta = 0.0; //SET
1430  int dimvar = (shouldIdoD18) ? dimRdown : dimRup ;
1431 
1432  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1433 
1434  dgemm_(&notrans,&var,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,F0block,&dimvar,&beta,workmem2,&dimLup);
1435 
1436  const double factor = Special::phase( TwoSLdown + 1 - TwoSLup ) * 0.5 * sqrt(0.5*(TwoSLup+1)*(TwoSLdown+1));
1437 
1438  int length = dimLup * dimRup;
1439  int inc = 1;
1440  total += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1441 
1442  }
1443  }
1444  }
1445  }
1446  }
1447  }
1448 
1449  return total;
1450 
1451 }
1452 
1453 double CheMPS2::TwoDM::doD19orD23(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD19){
1454 
1455  double total = 0.0;
1456 
1457  int theindex = denT->gIndex();
1458 
1459  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1460  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1461  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1462 
1463  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1464 
1465  if (dimLup>0){
1466 
1467  int ILdown = Irreps::directProd(ILup, Irrep_g);
1468  int IRdown = Irreps::directProd(ILdown, denBK->gIrrep(theindex));
1469 
1470  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1471  for (int TwoSRdown=TwoSLdown-1; TwoSRdown<=TwoSLdown+1; TwoSRdown+=2){
1472 
1473  int dimRup = denBK->gCurrentDim(theindex+1, NL, TwoSLup, ILup);
1474  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1475  int dimRdown = denBK->gCurrentDim(theindex+1, NL, TwoSRdown, IRdown);
1476 
1477  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1478 
1479  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL, TwoSLup, ILup);
1480  double * T_down = denT->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSRdown, IRdown);
1481  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, ILdown, NL, TwoSLup, ILup);
1482  double * F1block = (shouldIdoD19) ? F1right->gStorage(NL, TwoSRdown, IRdown, NL, TwoSLup, ILup)
1483  : F1right->gStorage(NL, TwoSLup, ILup, NL, TwoSRdown, IRdown) ;
1484 
1485  char trans = 'T';
1486  char notrans = 'N';
1487  char var = (shouldIdoD19) ? notrans : trans;
1488  double alpha = 1.0;
1489  double beta = 0.0; //SET
1490  int dimvar = (shouldIdoD19) ? dimRdown : dimRup ;
1491 
1492  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1493 
1494  dgemm_(&notrans,&var,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,F1block,&dimvar,&beta,workmem2,&dimLup);
1495 
1496  double factor = 0.0;
1497  if (shouldIdoD19){
1498  factor = Special::phase( TwoSLdown + TwoSRdown - 1 )
1499  * ( TwoSRdown + 1 ) * sqrt((TwoSLup+1)/3.0) * Wigner::wigner6j( 1, 1, 2, TwoSLup, TwoSRdown, TwoSLdown );
1500  } else {
1501  factor = Special::phase( TwoSLdown + TwoSLup - 1 )
1502  * ( TwoSLup + 1 ) * sqrt((TwoSRdown+1)/3.0) * Wigner::wigner6j( 1, 1, 2, TwoSLup, TwoSRdown, TwoSLdown );
1503  }
1504 
1505  int length = dimLup * dimRup;
1506  int inc = 1;
1507  total += factor * ddot_(&length,workmem2,&inc,T_up,&inc);
1508 
1509  }
1510  }
1511  }
1512  }
1513  }
1514  }
1515  }
1516 
1517  return total;
1518 
1519 }
1520 
1521 double CheMPS2::TwoDM::doD20orD24(TensorT * denT, TensorL * Lleft, TensorF1 * F1right, double * workmem, double * workmem2, int Irrep_g, bool shouldIdoD20){
1522 
1523  double total = 0.0;
1524 
1525  int theindex = denT->gIndex();
1526 
1527  for (int NL = denBK->gNmin(theindex); NL<=denBK->gNmax(theindex); NL++){
1528  for (int TwoSLup = denBK->gTwoSmin(theindex,NL); TwoSLup<= denBK->gTwoSmax(theindex,NL); TwoSLup+=2){
1529  for (int ILup = 0; ILup<denBK->getNumberOfIrreps(); ILup++){
1530 
1531  int dimLup = denBK->gCurrentDim(theindex, NL, TwoSLup, ILup);
1532 
1533  if (dimLup>0){
1534 
1535  int Idown = Irreps::directProd(ILup, Irrep_g);
1536  int IRup = Irreps::directProd(ILup, denBK->gIrrep(theindex));
1537 
1538  for (int TwoSLdown=TwoSLup-1; TwoSLdown<=TwoSLup+1; TwoSLdown+=2){
1539  for (int TwoSRup=TwoSLup-1; TwoSRup<=TwoSLup+1; TwoSRup+=2){
1540 
1541  int dimRup = denBK->gCurrentDim(theindex+1, NL+1, TwoSRup, IRup);
1542  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, Idown);
1543  int dimRdown = denBK->gCurrentDim(theindex+1, NL+1, TwoSLdown, Idown);
1544 
1545  if ((dimLdown>0) && (dimRdown>0) && (dimRup>0)){
1546 
1547  double * T_up = denT->gStorage(NL, TwoSLup, ILup, NL+1, TwoSRup, IRup);
1548  double * T_down = denT->gStorage(NL-1, TwoSLdown, Idown, NL+1, TwoSLdown, Idown);
1549  double * Lblock = Lleft->gStorage(NL-1, TwoSLdown, Idown, NL, TwoSLup, ILup);
1550  double * F1block = (shouldIdoD20) ? F1right->gStorage(NL+1, TwoSLdown, Idown, NL+1, TwoSRup, IRup)
1551  : F1right->gStorage(NL+1, TwoSRup, IRup, NL+1, TwoSLdown, Idown) ;
1552 
1553  char trans = 'T';
1554  char notrans = 'N';
1555  char var = (shouldIdoD20) ? notrans : trans;
1556  double alpha = 1.0;
1557  double beta = 0.0; //SET
1558  int dimvar = (shouldIdoD20) ? dimRdown : dimRup ;
1559 
1560  dgemm_(&trans,&notrans,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,T_down,&dimLdown,&beta,workmem,&dimLup);
1561 
1562  dgemm_(&notrans,&var,&dimLup,&dimRup,&dimRdown,&alpha,workmem,&dimLup,F1block,&dimvar,&beta,workmem2,&dimLup);
1563 
1564  double factor = 0.0;
1565  if (shouldIdoD20){
1566  factor = Special::phase( 2 * TwoSLup )
1567  * sqrt((TwoSLup+1)*(TwoSRup+1)*(TwoSLdown+1)/3.0) * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSLdown, TwoSLup );
1568  } else {
1569  factor = Special::phase( 2 * TwoSLup + TwoSRup - TwoSLdown )
1570  * (TwoSRup+1) * sqrt((TwoSLup+1)/3.0) * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSLdown, TwoSLup );
1571  }
1572 
1573  int length = dimLup * dimRup;
1574  int inc = 1;
1575  total += factor * ddot_(&length, workmem2, &inc, T_up, &inc);
1576 
1577  }
1578  }
1579  }
1580  }
1581  }
1582  }
1583  }
1584 
1585  return total;
1586 
1587 }
1588 
1589 
1590 
1591 
CheMPS2::TwoDM::FillSite
void FillSite(TensorT *denT, TensorL ***Ltens, TensorF0 ****F0tens, TensorF1 ****F1tens, TensorS0 ****S0tens, TensorS1 ****S1tens)
Fill the 2DM terms with as second site index denT->gIndex()
Definition: TwoDM.cpp:444
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::TwoDM::~TwoDM
virtual ~TwoDM()
Destructor.
Definition: TwoDM.cpp:56
CheMPS2::TwoDM::get1RDM_HAM
double get1RDM_HAM(const int cnt1, const int cnt2) const
Get a 1-RDM term, using the HAM indices.
Definition: TwoDM.cpp:184
CheMPS2::Problem::gReorder
bool gReorder() const
Get whether the Hamiltonian orbitals are reordered for the DMRG calculation.
Definition: Problem.cpp:346
CheMPS2::TensorT::gIndex
int gIndex() const
Get the location index.
Definition: TensorT.cpp:161
CheMPS2::TwoDM::save
void save() const
Save the TwoDMs to disk.
Definition: TwoDM.cpp:314
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::MPIchemps2::mpi_rank
static int mpi_rank()
Get the rank of this MPI process.
Definition: MPIchemps2.h:131
CheMPS2::TwoDM::getTwoDMB_DMRG
double getTwoDMB_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const
Get a 2DM_B term, using the DMRG indices.
Definition: TwoDM.cpp:112
CheMPS2::Irreps::symm_psi2molpro
void symm_psi2molpro(int *psi2molpro) const
Fill the array psi2molpro with the irrep conventions of molpro for the currently activated group...
Definition: Irreps.cpp:177
CheMPS2::Special::phase
static int phase(const int power_times_two)
Phase function.
Definition: Special.h:36
CheMPS2::TwoDM::energy
double energy() const
Calculate the energy based on the 2DM-A.
Definition: TwoDM.cpp:216
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::TwoDM::spin_density_dmrg
double spin_density_dmrg(const int cnt1, const int cnt2) const
Get a spin-density term, using the DMRG indices.
Definition: TwoDM.cpp:143
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::TwoDM::correct_higher_multiplicities
void correct_higher_multiplicities()
After the whole 2-RDM is filled, a prefactor for higher multiplicities should be applied.
Definition: TwoDM.cpp:629
CheMPS2::Problem::gN
int gN() const
Get the targeted particle number.
Definition: Problem.h:68
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::TwoDM::getTwoDMA_HAM
double getTwoDMA_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const
Get a 2DM_A term, using the HAM indices.
Definition: TwoDM.cpp:164
CheMPS2::Problem::gMxElement
double gMxElement(const int alpha, const int beta, const int gamma, const int delta) const
Get a specific interaction matrix element.
Definition: Problem.cpp:350
CheMPS2::TwoDM::get1RDM_DMRG
double get1RDM_DMRG(const int cnt1, const int cnt2) const
Get a 1-RDM term, using the DMRG indices.
Definition: TwoDM.cpp:127
CheMPS2::Problem::gIrrep
int gIrrep(const int nOrb) const
Get an orbital irrep.
Definition: Problem.cpp:336
CheMPS2::TwoDM::getTwoDMA_DMRG
double getTwoDMA_DMRG(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const
Get a 2DM_A term, using the DMRG indices.
Definition: TwoDM.cpp:97
CheMPS2::TwoDM::trace
double trace() const
Return the double trace of 2DM-A (should be N(N-1))
Definition: TwoDM.cpp:204
CheMPS2::TensorT::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorT.cpp:134
CheMPS2::Problem::gf1
int gf1(const int HamOrb) const
Get the DMRG index corresponding to a Ham index.
Definition: Problem.cpp:347
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::Problem::gL
int gL() const
Get the number of orbitals.
Definition: Problem.h:51
CheMPS2::TwoDM::spin_density_ham
double spin_density_ham(const int cnt1, const int cnt2) const
Get a spin-density term, using the HAM indices.
Definition: TwoDM.cpp:194
CheMPS2::TwoDM::mpi_allreduce
void mpi_allreduce()
Add the 2-RDM elements of all MPI processes.
CheMPS2::SyBookkeeper::getNumberOfIrreps
int getNumberOfIrreps() const
Get the total number of irreps.
Definition: SyBookkeeper.cpp:153
CheMPS2::TwoDM::getTwoDMB_HAM
double getTwoDMB_HAM(const int cnt1, const int cnt2, const int cnt3, const int cnt4) const
Get a 2DM_B term, using the HAM indices.
Definition: TwoDM.cpp:174
CheMPS2::TwoDM::print_noon
void print_noon() const
Print the natural orbital occupation numbers.
Definition: TwoDM.cpp:233
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::TwoDM::read
void read()
Load the TwoDMs from disk.
Definition: TwoDM.cpp:345
CheMPS2::TwoDM::save_HAM
void save_HAM(const string filename) const
Save the 2-RDM-A to disk in Hamiltonian indices.
Definition: TwoDM.cpp:278
CheMPS2::SyBookkeeper::gNmax
int gNmax(const int boundary) const
Get the max. possible particle number for a certain boundary.
Definition: SyBookkeeper.cpp:157
CheMPS2::TwoDM::TwoDM
TwoDM(const SyBookkeeper *denBKIn, const Problem *ProbIn)
Constructor.
Definition: TwoDM.cpp:39
CheMPS2::TwoDM::write2DMAfile
void write2DMAfile(const string filename) const
Write the 2-RDM-A to a file.
Definition: TwoDM.cpp:372
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145
CheMPS2::Problem::gEconst
double gEconst() const
Get the constant part of the Hamiltonian.
Definition: Problem.h:76
CheMPS2::Problem::gSy
int gSy() const
Get the point group symmetry.
Definition: Problem.h:55