CheMPS2
DMRGoperators3RDM.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 <stdlib.h>
21 #include <sys/time.h>
22 #include <assert.h>
23 
24 #include "DMRG.h"
25 #include "MPIchemps2.h"
26 #include "Special.h"
27 
28 void CheMPS2::DMRG::update_safe_3rdm_operators(const int boundary){
29 
30  /*
31  indices 0 <= j <= k <= l < boundary
32 
33  tensor_3rdm[ boundary - 1 == index ][ k - j ][ l - k ][ boundary - 1 - l ]
34 
35  **************************
36  * anni / anni / anni *
37  **************************
38 
39  1/ j == k == l is forbidden ( no three annihilators on the same site )
40  2/ if j == k, J1 must be zero ( Sigma_J1 does not exist then )
41  3/ if k == l, J2 must be 1/2
42 
43  tensor_3rdm_a_J0_doublet --> j <= k < l & j < k == l ( or NOT j == k == l )
44  tensor_3rdm_a_J1_doublet --> j < k <= l
45  tensor_3rdm_a_J1_quartet --> j < k < l
46 
47  **************************
48  * anni / anni / crea *
49  **************************
50 
51  1/ j <= k < l is allowed ( k == l is part of tensor_3rdm_c )
52  2/ if j == k, J1 must be zero ( Sigma_J1 does not exist then )
53 
54  tensor_3rdm_b_J0_doublet --> j <= k < l
55  tensor_3rdm_b_J1_doublet --> j < k < l
56  tensor_3rdm_b_J1_quartet --> j < k < l
57 
58  **************************
59  * anni / crea / anni *
60  **************************
61 
62  1/ j == k == l is forbidden ( j == k == l is part of tensor_3rdm_d )
63 
64  tensor_3rdm_c_J0_doublet --> j <= k < l & j < k == l ( or NOT j == k == l )
65  tensor_3rdm_c_J1_doublet --> j <= k < l & j < k == l ( or NOT j == k == l )
66  tensor_3rdm_c_J1_quartet --> j <= k < l & j < k == l ( or NOT j == k == l )
67 
68  **************************
69  * crea / anni / anni *
70  **************************
71 
72  1/ j <= k <= l is allowed
73  2/ if k == l, J2 must be 1/2
74 
75  tensor_3rdm_d_J0_doublet --> j <= k <= l
76  tensor_3rdm_d_J1_doublet --> j <= k <= l
77  tensor_3rdm_d_J1_quartet --> j <= k < l
78  */
79 
80  allocate_3rdm_operators( boundary );
81  update_3rdm_operators( boundary );
82  if ( boundary >= 2 ){ delete_3rdm_operators( boundary - 1 ); }
83 
84 }
85 
86 void CheMPS2::DMRG::update_3rdm_operators(const int boundary){
87 
88  struct timeval start, end;
89  gettimeofday(&start, NULL);
90 
91  const int index = boundary - 1;
92  const int dimL = denBK->gMaxDimAtBound(boundary-1);
93  const int dimR = denBK->gMaxDimAtBound(boundary);
94  #ifdef CHEMPS2_MPI_COMPILATION
95  const int MPIRANK = MPIchemps2::mpi_rank();
96  #endif
97 
98  #pragma omp parallel
99  {
100  double * workmem = new double[dimL*dimR];
101 
102 #ifdef CHEMPS2_MPI_COMPILATION //######( loop j<=k<=l MPI )######//
103 
104 /* Strategy for MPI:
105  - outer loop is (j,k)
106  - everyone has a temporary duplicate of S_jk and F_jk
107 */
108  for ( int orb_j = 0; orb_j < boundary; orb_j++ ){
109  for ( int orb_k = orb_j; orb_k < boundary; orb_k++ ){
110 
111  const int irrjk = Irreps::directProd( denBK->gIrrep( orb_j ), denBK->gIrrep( orb_k ) );
112  const int cnt1 = orb_k - orb_j;
113 
114  #pragma omp single
115  if ( orb_k < index-1 ){ // All processes own Fx/Sx[ index - 1 ][ k - j ][ index - 1 - k == 0 ]
116  const int own_S_jk = MPIchemps2::owner_absigma( orb_j, orb_k );
117  const int own_F_jk = MPIchemps2::owner_cdf( L, orb_j, orb_k );
118  if ( MPIRANK != own_F_jk ){ F0tensors[index-1][cnt1][index-orb_k-1] = new TensorF0( index, irrjk, true, denBK );
119  F1tensors[index-1][cnt1][index-orb_k-1] = new TensorF1( index, irrjk, true, denBK ); }
120  if ( MPIRANK != own_S_jk ){ S0tensors[index-1][cnt1][index-orb_k-1] = new TensorS0( index, irrjk, true, denBK );
121  if ( cnt1 > 0 ){ S1tensors[index-1][cnt1][index-orb_k-1] = new TensorS1( index, irrjk, true, denBK ); }}
122  MPIchemps2::broadcast_tensor( F0tensors[index-1][cnt1][index-orb_k-1], own_F_jk );
123  MPIchemps2::broadcast_tensor( F1tensors[index-1][cnt1][index-orb_k-1], own_F_jk );
124  MPIchemps2::broadcast_tensor( S0tensors[index-1][cnt1][index-orb_k-1], own_S_jk );
125  if ( cnt1 > 0 ){ MPIchemps2::broadcast_tensor( S1tensors[index-1][cnt1][index-orb_k-1], own_S_jk ); }
126  }
127 
128  #pragma omp for schedule(dynamic)
129  for ( int orb_l = orb_k; orb_l < boundary; orb_l++ ){
130  if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK ){
131  const int cnt2 = orb_l - orb_k;
132  const int cnt3 = index - orb_l;
133 
134 #else //######( loop j<=k<=l MPI )######//
135 
136  const int upperbound = (boundary*(boundary+1)*(boundary+2))/6;
137  int jkl[] = { 0, 0, 0 };
138  #pragma omp for schedule(static)
139  for ( int global = 0; global < upperbound; global++ ){
140  Special::invert_triangle_three( global, jkl );
141  const int orb_j = jkl[ 0 ];
142  const int orb_k = jkl[ 1 ];
143  const int orb_l = jkl[ 2 ];
144  const int recalculate_global = orb_j + (orb_k*(orb_k+1))/2 + (orb_l*(orb_l+1)*(orb_l+2))/6;
145  assert( global == recalculate_global );
146  const int cnt1 = orb_k - orb_j;
147  const int cnt2 = orb_l - orb_k;
148  const int cnt3 = index - orb_l;
149 
150 #endif //######( loop j<=k<=l MPI )######//
151 
152  /* PERFORM THE UPDATES */
153  if ( cnt3 == 0 ){ // Create tensors
154  if ( cnt2 > 0 ){
155  tensor_3rdm_a_J0_doublet[index][cnt1][cnt2][0]->a1(S0tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
156  if (cnt1>0){ tensor_3rdm_a_J1_doublet[index][cnt1][cnt2][0]->a1(S1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
157  tensor_3rdm_a_J1_quartet[index][cnt1][cnt2][0]->a1(S1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem); }
158  tensor_3rdm_b_J0_doublet[index][cnt1][cnt2][0]->b1(S0tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
159  if (cnt1>0){ tensor_3rdm_b_J1_doublet[index][cnt1][cnt2][0]->b1(S1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
160  tensor_3rdm_b_J1_quartet[index][cnt1][cnt2][0]->b1(S1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem); }
161  tensor_3rdm_c_J0_doublet[index][cnt1][cnt2][0]->c1(F0tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
162  tensor_3rdm_c_J1_doublet[index][cnt1][cnt2][0]->c1(F1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
163  tensor_3rdm_c_J1_quartet[index][cnt1][cnt2][0]->c1(F1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
164  tensor_3rdm_d_J0_doublet[index][cnt1][cnt2][0]->d1(F0tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
165  tensor_3rdm_d_J1_doublet[index][cnt1][cnt2][0]->d1(F1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
166  tensor_3rdm_d_J1_quartet[index][cnt1][cnt2][0]->d1(F1tensors[index-1][cnt1][cnt2-1], MPS[index], workmem);
167  } else {
168  if ( cnt1 > 0 ){
169  tensor_3rdm_a_J0_doublet[index][cnt1][0][0]->extra2(Ltensors[index-1][cnt1-1], MPS[index], workmem);
170  tensor_3rdm_a_J1_doublet[index][cnt1][0][0]->extra2(Ltensors[index-1][cnt1-1], MPS[index], workmem);
171  tensor_3rdm_c_J0_doublet[index][cnt1][0][0]->extra4(Ltensors[index-1][cnt1-1], MPS[index], workmem);
172  tensor_3rdm_c_J1_doublet[index][cnt1][0][0]->extra4(Ltensors[index-1][cnt1-1], MPS[index], workmem);
173  tensor_3rdm_c_J1_quartet[index][cnt1][0][0]->extra4(Ltensors[index-1][cnt1-1], MPS[index], workmem);
174  tensor_3rdm_d_J0_doublet[index][cnt1][0][0]->extra3(Ltensors[index-1][cnt1-1], MPS[index], workmem);
175  tensor_3rdm_d_J1_doublet[index][cnt1][0][0]->extra3(Ltensors[index-1][cnt1-1], MPS[index], workmem);
176  } else {
177  tensor_3rdm_d_J0_doublet[index][0][0][0]->extra1(MPS[index]);
178  tensor_3rdm_d_J1_doublet[index][0][0][0]->extra1(MPS[index]);
179  }
180  }
181  } else { // Update tensors
182  if (cnt1+cnt2>0){ tensor_3rdm_a_J0_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_a_J0_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
183  if (cnt1>0) { tensor_3rdm_a_J1_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_a_J1_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
184  if (cnt1*cnt2>0){ tensor_3rdm_a_J1_quartet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_a_J1_quartet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
185  if (cnt2>0) { tensor_3rdm_b_J0_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_b_J0_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
186  if (cnt1*cnt2>0){ tensor_3rdm_b_J1_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_b_J1_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
187  if (cnt1*cnt2>0){ tensor_3rdm_b_J1_quartet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_b_J1_quartet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
188  if (cnt1+cnt2>0){ tensor_3rdm_c_J0_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_c_J0_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
189  if (cnt1+cnt2>0){ tensor_3rdm_c_J1_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_c_J1_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
190  if (cnt1+cnt2>0){ tensor_3rdm_c_J1_quartet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_c_J1_quartet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
191  tensor_3rdm_d_J0_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_d_J0_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem);
192  tensor_3rdm_d_J1_doublet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_d_J1_doublet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem);
193  if (cnt2>0) { tensor_3rdm_d_J1_quartet[index][cnt1][cnt2][cnt3]->update(tensor_3rdm_d_J1_quartet[index-1][cnt1][cnt2][cnt3-1], MPS[index], MPS[index], workmem); }
194  }
195 
196 #ifdef CHEMPS2_MPI_COMPILATION //######( close loop j<=k<=l MPI )######//
197 
198  }
199  }
200 
201  #pragma omp single
202  if ( orb_k < index - 1 ){ // All processes own Fx/Sx[ index - 1 ][ k - j ][ index - 1 - k == 0 ]
203  const int own_S_jk = MPIchemps2::owner_absigma( orb_j, orb_k );
204  const int own_F_jk = MPIchemps2::owner_cdf( L, orb_j, orb_k );
205  if ( MPIRANK != own_F_jk ){ delete F0tensors[index-1][cnt1][index-orb_k-1]; F0tensors[index-1][cnt1][index-orb_k-1] = NULL;
206  delete F1tensors[index-1][cnt1][index-orb_k-1]; F1tensors[index-1][cnt1][index-orb_k-1] = NULL; }
207  if ( MPIRANK != own_S_jk ){ delete S0tensors[index-1][cnt1][index-orb_k-1]; S0tensors[index-1][cnt1][index-orb_k-1] = NULL;
208  if ( cnt1 > 0 ){ delete S1tensors[index-1][cnt1][index-orb_k-1]; S1tensors[index-1][cnt1][index-orb_k-1] = NULL; }}
209  }
210  }
211  }
212 
213 #else //######( close loop j<=k<=l MPI )######//
214 
215  }
216 
217 #endif //######( close loop j<=k<=l MPI )######//
218 
219  delete [] workmem;
220  }
221 
222  gettimeofday(&end, NULL);
223  timings[ CHEMPS2_TIME_TENS_CALC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
224 
225 }
226 
227 void CheMPS2::DMRG::allocate_3rdm_operators(const int boundary){
228 
229  struct timeval start, end;
230  gettimeofday(&start, NULL);
231 
232  #ifdef CHEMPS2_MPI_COMPILATION
233  const int MPIRANK = MPIchemps2::mpi_rank();
234  #endif
235  const int index = boundary - 1;
236 
237  tensor_3rdm_a_J0_doublet[ index ] = new Tensor3RDM***[ boundary ];
238  tensor_3rdm_a_J1_doublet[ index ] = new Tensor3RDM***[ boundary ];
239  tensor_3rdm_a_J1_quartet[ index ] = new Tensor3RDM***[ boundary ];
240  tensor_3rdm_b_J0_doublet[ index ] = new Tensor3RDM***[ boundary ];
241  tensor_3rdm_b_J1_doublet[ index ] = new Tensor3RDM***[ boundary ];
242  tensor_3rdm_b_J1_quartet[ index ] = new Tensor3RDM***[ boundary ];
243  tensor_3rdm_c_J0_doublet[ index ] = new Tensor3RDM***[ boundary ];
244  tensor_3rdm_c_J1_doublet[ index ] = new Tensor3RDM***[ boundary ];
245  tensor_3rdm_c_J1_quartet[ index ] = new Tensor3RDM***[ boundary ];
246  tensor_3rdm_d_J0_doublet[ index ] = new Tensor3RDM***[ boundary ];
247  tensor_3rdm_d_J1_doublet[ index ] = new Tensor3RDM***[ boundary ];
248  tensor_3rdm_d_J1_quartet[ index ] = new Tensor3RDM***[ boundary ];
249 
250  for ( int cnt1 = 0; cnt1 < boundary; cnt1++ ){ // cnt1 = k - j < boundary
251 
252  tensor_3rdm_a_J0_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
253  tensor_3rdm_a_J1_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
254  tensor_3rdm_a_J1_quartet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
255  tensor_3rdm_b_J0_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
256  tensor_3rdm_b_J1_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
257  tensor_3rdm_b_J1_quartet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
258  tensor_3rdm_c_J0_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
259  tensor_3rdm_c_J1_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
260  tensor_3rdm_c_J1_quartet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
261  tensor_3rdm_d_J0_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
262  tensor_3rdm_d_J1_doublet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
263  tensor_3rdm_d_J1_quartet[ index ][ cnt1 ] = new Tensor3RDM**[ boundary - cnt1 ];
264 
265  for ( int cnt2 = 0; cnt2 < boundary - cnt1; cnt2++ ){ // cnt2 = l - k < boundary - k = boundary - cnt1 - j <= boundary - cnt1
266 
267  tensor_3rdm_a_J0_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
268  tensor_3rdm_a_J1_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
269  tensor_3rdm_a_J1_quartet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
270  tensor_3rdm_b_J0_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
271  tensor_3rdm_b_J1_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
272  tensor_3rdm_b_J1_quartet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
273  tensor_3rdm_c_J0_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
274  tensor_3rdm_c_J1_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
275  tensor_3rdm_c_J1_quartet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
276  tensor_3rdm_d_J0_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
277  tensor_3rdm_d_J1_doublet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
278  tensor_3rdm_d_J1_quartet[ index ][ cnt1 ][ cnt2 ] = new Tensor3RDM*[ boundary - cnt1 - cnt2 ];
279 
280  for ( int cnt3 = 0; cnt3 < boundary - cnt1 - cnt2; cnt3++ ){ // cnt3 = boundary - 1 - l < boundary - ( l - k ) - ( k - j ) = boundary - cnt1 - cnt2
281 
282  const int orb_l = boundary - 1 - cnt3;
283  const int orb_k = orb_l - cnt2;
284  const int orb_j = orb_k - cnt1;
285  const int irr = Irreps::directProd( Irreps::directProd( denBK->gIrrep( orb_j ), denBK->gIrrep( orb_k ) ), denBK->gIrrep( orb_l ) );
286 
287  #ifdef CHEMPS2_MPI_COMPILATION
288  if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK ){
289  #endif
290  tensor_3rdm_a_J0_doublet[index][cnt1][cnt2][cnt3] = (cnt1+cnt2>0) ? new Tensor3RDM(boundary, 0, 1, 3, irr, true, denBK) : NULL; // NOT j == k == l
291  tensor_3rdm_a_J1_doublet[index][cnt1][cnt2][cnt3] = (cnt1>0) ? new Tensor3RDM(boundary, 2, 1, 3, irr, true, denBK) : NULL; // j < k <= l
292  tensor_3rdm_a_J1_quartet[index][cnt1][cnt2][cnt3] = (cnt1*cnt2>0) ? new Tensor3RDM(boundary, 2, 3, 3, irr, true, denBK) : NULL; // j < k < l
293  tensor_3rdm_b_J0_doublet[index][cnt1][cnt2][cnt3] = (cnt2>0) ? new Tensor3RDM(boundary, 0, 1, 1, irr, true, denBK) : NULL; // j <= k < l
294  tensor_3rdm_b_J1_doublet[index][cnt1][cnt2][cnt3] = (cnt1*cnt2>0) ? new Tensor3RDM(boundary, 2, 1, 1, irr, true, denBK) : NULL; // j < k < l
295  tensor_3rdm_b_J1_quartet[index][cnt1][cnt2][cnt3] = (cnt1*cnt2>0) ? new Tensor3RDM(boundary, 2, 3, 1, irr, true, denBK) : NULL; // j < k < l
296  tensor_3rdm_c_J0_doublet[index][cnt1][cnt2][cnt3] = (cnt1+cnt2>0) ? new Tensor3RDM(boundary, 0, 1, 1, irr, true, denBK) : NULL; // NOT j == k == l
297  tensor_3rdm_c_J1_doublet[index][cnt1][cnt2][cnt3] = (cnt1+cnt2>0) ? new Tensor3RDM(boundary, 2, 1, 1, irr, true, denBK) : NULL; // NOT j == k == l
298  tensor_3rdm_c_J1_quartet[index][cnt1][cnt2][cnt3] = (cnt1+cnt2>0) ? new Tensor3RDM(boundary, 2, 3, 1, irr, true, denBK) : NULL; // NOT j == k == l
299  tensor_3rdm_d_J0_doublet[index][cnt1][cnt2][cnt3] = new Tensor3RDM(boundary, 0, 1, 1, irr, false, denBK); // j <= k <= l
300  tensor_3rdm_d_J1_doublet[index][cnt1][cnt2][cnt3] = new Tensor3RDM(boundary, 2, 1, 1, irr, false, denBK); // j <= k <= l
301  tensor_3rdm_d_J1_quartet[index][cnt1][cnt2][cnt3] = (cnt2>0) ? new Tensor3RDM(boundary, 2, 3, 1, irr, false, denBK) : NULL; // j <= k < l
302  #ifdef CHEMPS2_MPI_COMPILATION
303  } else {
304  tensor_3rdm_a_J0_doublet[index][cnt1][cnt2][cnt3] = NULL;
305  tensor_3rdm_a_J1_doublet[index][cnt1][cnt2][cnt3] = NULL;
306  tensor_3rdm_a_J1_quartet[index][cnt1][cnt2][cnt3] = NULL;
307  tensor_3rdm_b_J0_doublet[index][cnt1][cnt2][cnt3] = NULL;
308  tensor_3rdm_b_J1_doublet[index][cnt1][cnt2][cnt3] = NULL;
309  tensor_3rdm_b_J1_quartet[index][cnt1][cnt2][cnt3] = NULL;
310  tensor_3rdm_c_J0_doublet[index][cnt1][cnt2][cnt3] = NULL;
311  tensor_3rdm_c_J1_doublet[index][cnt1][cnt2][cnt3] = NULL;
312  tensor_3rdm_c_J1_quartet[index][cnt1][cnt2][cnt3] = NULL;
313  tensor_3rdm_d_J0_doublet[index][cnt1][cnt2][cnt3] = NULL;
314  tensor_3rdm_d_J1_doublet[index][cnt1][cnt2][cnt3] = NULL;
315  tensor_3rdm_d_J1_quartet[index][cnt1][cnt2][cnt3] = NULL;
316  }
317  #endif
318 
319  }
320  }
321  }
322 
323  gettimeofday(&end, NULL);
324  timings[ CHEMPS2_TIME_TENS_ALLOC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
325 
326 }
327 
328 void CheMPS2::DMRG::delete_3rdm_operators(const int boundary){
329 
330  struct timeval start, end;
331  gettimeofday(&start, NULL);
332 
333  #ifdef CHEMPS2_MPI_COMPILATION
334  const int MPIRANK = MPIchemps2::mpi_rank();
335  #endif
336  const int index = boundary - 1;
337 
338  for ( int cnt1 = 0; cnt1 < boundary; cnt1++ ){ // cnt1 = k - j < boundary
339 
340  for ( int cnt2 = 0; cnt2 < boundary - cnt1; cnt2++ ){ // cnt2 = l - k < boundary - k = boundary - cnt1 - j <= boundary - cnt1
341 
342  for ( int cnt3 = 0; cnt3 < boundary - cnt1 - cnt2; cnt3++ ){ // cnt3 = boundary - 1 - l < boundary - ( l - k ) - ( k - j ) = boundary - cnt1 - cnt2
343 
344  const int orb_l = boundary - 1 - cnt3;
345  const int orb_k = orb_l - cnt2;
346  const int orb_j = orb_k - cnt1;
347 
348  #ifdef CHEMPS2_MPI_COMPILATION
349  if ( MPIchemps2::owner_3rdm_diagram( L, orb_j, orb_k, orb_l ) == MPIRANK )
350  #endif
351  {
352  if (cnt1+cnt2>0){ delete tensor_3rdm_a_J0_doublet[index][cnt1][cnt2][cnt3]; }
353  if (cnt1>0) { delete tensor_3rdm_a_J1_doublet[index][cnt1][cnt2][cnt3]; }
354  if (cnt1*cnt2>0){ delete tensor_3rdm_a_J1_quartet[index][cnt1][cnt2][cnt3]; }
355  if (cnt2>0) { delete tensor_3rdm_b_J0_doublet[index][cnt1][cnt2][cnt3]; }
356  if (cnt1*cnt2>0){ delete tensor_3rdm_b_J1_doublet[index][cnt1][cnt2][cnt3]; }
357  if (cnt1*cnt2>0){ delete tensor_3rdm_b_J1_quartet[index][cnt1][cnt2][cnt3]; }
358  if (cnt1+cnt2>0){ delete tensor_3rdm_c_J0_doublet[index][cnt1][cnt2][cnt3]; }
359  if (cnt1+cnt2>0){ delete tensor_3rdm_c_J1_doublet[index][cnt1][cnt2][cnt3]; }
360  if (cnt1+cnt2>0){ delete tensor_3rdm_c_J1_quartet[index][cnt1][cnt2][cnt3]; }
361  delete tensor_3rdm_d_J0_doublet[index][cnt1][cnt2][cnt3];
362  delete tensor_3rdm_d_J1_doublet[index][cnt1][cnt2][cnt3];
363  if (cnt2>0) { delete tensor_3rdm_d_J1_quartet[index][cnt1][cnt2][cnt3]; }
364  }
365  }
366 
367  delete [] tensor_3rdm_a_J0_doublet[ index ][ cnt1 ][ cnt2 ];
368  delete [] tensor_3rdm_a_J1_doublet[ index ][ cnt1 ][ cnt2 ];
369  delete [] tensor_3rdm_a_J1_quartet[ index ][ cnt1 ][ cnt2 ];
370  delete [] tensor_3rdm_b_J0_doublet[ index ][ cnt1 ][ cnt2 ];
371  delete [] tensor_3rdm_b_J1_doublet[ index ][ cnt1 ][ cnt2 ];
372  delete [] tensor_3rdm_b_J1_quartet[ index ][ cnt1 ][ cnt2 ];
373  delete [] tensor_3rdm_c_J0_doublet[ index ][ cnt1 ][ cnt2 ];
374  delete [] tensor_3rdm_c_J1_doublet[ index ][ cnt1 ][ cnt2 ];
375  delete [] tensor_3rdm_c_J1_quartet[ index ][ cnt1 ][ cnt2 ];
376  delete [] tensor_3rdm_d_J0_doublet[ index ][ cnt1 ][ cnt2 ];
377  delete [] tensor_3rdm_d_J1_doublet[ index ][ cnt1 ][ cnt2 ];
378  delete [] tensor_3rdm_d_J1_quartet[ index ][ cnt1 ][ cnt2 ];
379 
380  }
381 
382  delete [] tensor_3rdm_a_J0_doublet[ index ][ cnt1 ];
383  delete [] tensor_3rdm_a_J1_doublet[ index ][ cnt1 ];
384  delete [] tensor_3rdm_a_J1_quartet[ index ][ cnt1 ];
385  delete [] tensor_3rdm_b_J0_doublet[ index ][ cnt1 ];
386  delete [] tensor_3rdm_b_J1_doublet[ index ][ cnt1 ];
387  delete [] tensor_3rdm_b_J1_quartet[ index ][ cnt1 ];
388  delete [] tensor_3rdm_c_J0_doublet[ index ][ cnt1 ];
389  delete [] tensor_3rdm_c_J1_doublet[ index ][ cnt1 ];
390  delete [] tensor_3rdm_c_J1_quartet[ index ][ cnt1 ];
391  delete [] tensor_3rdm_d_J0_doublet[ index ][ cnt1 ];
392  delete [] tensor_3rdm_d_J1_doublet[ index ][ cnt1 ];
393  delete [] tensor_3rdm_d_J1_quartet[ index ][ cnt1 ];
394 
395  }
396 
397  delete [] tensor_3rdm_a_J0_doublet[ index ];
398  delete [] tensor_3rdm_a_J1_doublet[ index ];
399  delete [] tensor_3rdm_a_J1_quartet[ index ];
400  delete [] tensor_3rdm_b_J0_doublet[ index ];
401  delete [] tensor_3rdm_b_J1_doublet[ index ];
402  delete [] tensor_3rdm_b_J1_quartet[ index ];
403  delete [] tensor_3rdm_c_J0_doublet[ index ];
404  delete [] tensor_3rdm_c_J1_doublet[ index ];
405  delete [] tensor_3rdm_c_J1_quartet[ index ];
406  delete [] tensor_3rdm_d_J0_doublet[ index ];
407  delete [] tensor_3rdm_d_J1_doublet[ index ];
408  delete [] tensor_3rdm_d_J1_quartet[ index ];
409 
410  gettimeofday(&end, NULL);
411  timings[ CHEMPS2_TIME_TENS_FREE ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
412 
413 }
414 
415 void CheMPS2::DMRG::update_correlations_tensors(const int siteindex){
416 
417  struct timeval start, end;
418 
419  const int dimL = denBK->gMaxDimAtBound(siteindex-1);
420  const int dimR = denBK->gMaxDimAtBound(siteindex);
421  double * workmemLR = new double[dimL*dimR];
422 
423  for ( int previousindex = 0; previousindex < siteindex-1; previousindex++ ){
424 
425  gettimeofday(&start, NULL);
426  TensorGYZ * newG = new TensorGYZ(siteindex, 'G', denBK);
427  TensorGYZ * newY = new TensorGYZ(siteindex, 'Y', denBK);
428  TensorGYZ * newZ = new TensorGYZ(siteindex, 'Z', denBK);
429  TensorKM * newK = new TensorKM( siteindex, 'K', denBK->gIrrep(previousindex), denBK );
430  TensorKM * newM = new TensorKM( siteindex, 'M', denBK->gIrrep(previousindex), denBK );
431  gettimeofday(&end, NULL);
432  timings[ CHEMPS2_TIME_TENS_ALLOC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
433 
434  gettimeofday(&start, NULL);
435  newG->update(Gtensors[previousindex], MPS[siteindex-1], MPS[siteindex-1], workmemLR);
436  newY->update(Ytensors[previousindex], MPS[siteindex-1], MPS[siteindex-1], workmemLR);
437  newZ->update(Ztensors[previousindex], MPS[siteindex-1], MPS[siteindex-1], workmemLR);
438  newK->update(Ktensors[previousindex], MPS[siteindex-1], MPS[siteindex-1], workmemLR);
439  newM->update(Mtensors[previousindex], MPS[siteindex-1], MPS[siteindex-1], workmemLR);
440  gettimeofday(&end, NULL);
441  timings[ CHEMPS2_TIME_TENS_CALC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
442 
443  gettimeofday(&start, NULL);
444  delete Gtensors[previousindex];
445  delete Ytensors[previousindex];
446  delete Ztensors[previousindex];
447  delete Ktensors[previousindex];
448  delete Mtensors[previousindex];
449  gettimeofday(&end, NULL);
450  timings[ CHEMPS2_TIME_TENS_FREE ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
451 
452  Gtensors[previousindex] = newG;
453  Ytensors[previousindex] = newY;
454  Ztensors[previousindex] = newZ;
455  Ktensors[previousindex] = newK;
456  Mtensors[previousindex] = newM;
457 
458  }
459  delete [] workmemLR;
460 
461  gettimeofday(&start, NULL);
462  Gtensors[siteindex-1] = new TensorGYZ(siteindex, 'G', denBK);
463  Ytensors[siteindex-1] = new TensorGYZ(siteindex, 'Y', denBK);
464  Ztensors[siteindex-1] = new TensorGYZ(siteindex, 'Z', denBK);
465  Ktensors[siteindex-1] = new TensorKM( siteindex, 'K', denBK->gIrrep(siteindex-1), denBK );
466  Mtensors[siteindex-1] = new TensorKM( siteindex, 'M', denBK->gIrrep(siteindex-1), denBK );
467  gettimeofday(&end, NULL);
468  timings[ CHEMPS2_TIME_TENS_ALLOC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
469 
470  gettimeofday(&start, NULL);
471  Gtensors[siteindex-1]->construct(MPS[siteindex-1]);
472  Ytensors[siteindex-1]->construct(MPS[siteindex-1]);
473  Ztensors[siteindex-1]->construct(MPS[siteindex-1]);
474  Ktensors[siteindex-1]->construct(MPS[siteindex-1]);
475  Mtensors[siteindex-1]->construct(MPS[siteindex-1]);
476  gettimeofday(&end, NULL);
477  timings[ CHEMPS2_TIME_TENS_CALC ] += (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
478 
479 }
480 
481 
CheMPS2::TensorOperator::update
void update(TensorOperator *previous, TensorT *mps_tensor_up, TensorT *mps_tensor_down, double *workmem)
Clear and update.
Definition: TensorOperator.cpp:163
CheMPS2::Tensor3RDM::extra1
void extra1(TensorT *denT)
Make diagram extra1.
Definition: Tensor3RDM.cpp:356
CheMPS2::MPIchemps2::mpi_rank
static int mpi_rank()
Get the rank of this MPI process.
Definition: MPIchemps2.h:131
CheMPS2::Tensor3RDM::c1
void c1(TensorOperator *denF, TensorT *denT, double *workmem)
Make diagram c1.
Definition: Tensor3RDM.cpp:204
CheMPS2::Tensor3RDM::extra3
void extra3(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra3.
Definition: Tensor3RDM.cpp:430
CheMPS2::Tensor3RDM::a1
void a1(TensorOperator *Sigma, TensorT *denT, double *workmem)
Make diagram a1.
Definition: Tensor3RDM.cpp:52
CheMPS2::Tensor3RDM::extra2
void extra2(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra2.
Definition: Tensor3RDM.cpp:391
CheMPS2::Tensor3RDM::d1
void d1(TensorOperator *denF, TensorT *denT, double *workmem)
Make diagram d1.
Definition: Tensor3RDM.cpp:280
CheMPS2::TensorGYZ::construct
void construct(TensorT *denT)
Construct a new TensorGYZ.
Definition: TensorGYZ.cpp:42
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::Tensor3RDM::extra4
void extra4(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra4.
Definition: Tensor3RDM.cpp:469
CheMPS2::SyBookkeeper::gMaxDimAtBound
int gMaxDimAtBound(const int boundary) const
Get the maximum virtual dimension at a certain boundary.
Definition: SyBookkeeper.cpp:282
CheMPS2::Tensor3RDM::b1
void b1(TensorOperator *Sigma, TensorT *denT, double *workmem)
Make diagram b1.
Definition: Tensor3RDM.cpp:128
CheMPS2::TensorKM::construct
void construct(TensorT *denT)
Definition: TensorKM.cpp:42
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