CheMPS2
Cumulant.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 <unistd.h>
23 #include <iostream>
24 
25 #include "Cumulant.h"
26 
27 /*void CheMPS2::Cumulant::gamma4_fock_contract_ham_slow(const Problem * prob, const ThreeDM * the3DM, const TwoDM * the2DM, double * fock, double * result){
28 
29  struct timeval start, end;
30  gettimeofday(&start, NULL);
31  const int L = prob->gL();
32 
33  for ( int cnt = 0; cnt < L*L*L*L*L*L; cnt++ ){ result[ cnt ] = 0.0; }
34 
35  int * irreps = new int[ L ];
36  for ( int orb = 0; orb < L; orb++ ){ irreps[ orb ] = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( orb ) : orb ); }
37 
38  #pragma omp parallel for schedule(dynamic)
39  for ( int i = 0; i < L; i++ ){
40  for ( int j = i; j < L; j++ ){
41  for ( int k = i; k < L; k++ ){
42  const int irrep_ijk = Irreps::directProd( Irreps::directProd( irreps[ i ], irreps[ j ] ), irreps[ k ] );
43  for ( int p = i; p < L; p++ ){
44  for ( int q = i; q < L; q++ ){
45  for ( int r = q; r < L; r++ ){
46  const int irrep_pqr = Irreps::directProd( Irreps::directProd( irreps[ p ], irreps[ q ] ), irreps[ r ] );
47  if ( irrep_ijk == irrep_pqr ){
48  double value = 0.0;
49  for ( int l = 0; l < L; l++ ){
50  for ( int s = 0; s < L; s++ ){
51  if ( irreps[ l ] == irreps[ s ] ){
52  value += fock[ l + L * s ] * gamma4_ham(prob, the3DM, the2DM, i, j, k, l, p, q, r, s);
53  }
54  }
55  }
56  result[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] = value;
57  result[ i + L * ( k + L * ( j + L * ( p + L * ( r + L * q )))) ] = value;
58  result[ j + L * ( i + L * ( k + L * ( q + L * ( p + L * r )))) ] = value;
59  result[ k + L * ( i + L * ( j + L * ( r + L * ( p + L * q )))) ] = value;
60  result[ j + L * ( k + L * ( i + L * ( q + L * ( r + L * p )))) ] = value;
61  result[ k + L * ( j + L * ( i + L * ( r + L * ( q + L * p )))) ] = value;
62  result[ p + L * ( q + L * ( r + L * ( i + L * ( j + L * k )))) ] = value;
63  result[ p + L * ( r + L * ( q + L * ( i + L * ( k + L * j )))) ] = value;
64  result[ q + L * ( p + L * ( r + L * ( j + L * ( i + L * k )))) ] = value;
65  result[ r + L * ( p + L * ( q + L * ( k + L * ( i + L * j )))) ] = value;
66  result[ q + L * ( r + L * ( p + L * ( j + L * ( k + L * i )))) ] = value;
67  result[ r + L * ( q + L * ( p + L * ( k + L * ( j + L * i )))) ] = value;
68  }
69  }
70  }
71  }
72  }
73  }
74  }
75 
76  delete [] irreps;
77 
78  gettimeofday(&end, NULL);
79  const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
80  std::cout << "Cumulant :: Contraction of cu(4)-4RDM with CASPT2 Fock operator took " << elapsed << " seconds." << std::endl;
81 
82 }*/
83 
84 void CheMPS2::Cumulant::gamma4_fock_contract_ham(const Problem * prob, const ThreeDM * the3DM, const TwoDM * the2DM, double * fock, double * result){
85 
86  struct timeval start, end;
87  gettimeofday(&start, NULL);
88  const int L = prob->gL();
89 
90  /* Clear result */
91  for ( int cnt = 0; cnt < L*L*L*L*L*L; cnt++ ){ result[ cnt ] = 0.0; }
92 
93  /* Construct an array with the orbital irreps in Hamiltonian indices */
94  int * irreps = new int[ L ];
95  for ( int orb = 0; orb < L; orb++ ){ irreps[ orb ] = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( orb ) : orb ); }
96 
97  /* Helper arrays with partial (mult) or full (dot) contractions of objects with the CASPT2 Fock operator */
98  double * G3dotF = new double[ L*L*L*L ];
99  double * lambda2 = new double[ L*L*L*L ];
100  double * G2multF = new double[ L*L*L*L ];
101  double * L2multF = new double[ L*L*L*L ];
102  double * gamma1 = new double[ L*L ];
103  double * G1multF = new double[ L*L ];
104  double * G2dotF = new double[ L*L ];
105  double * L2dotF = new double[ L*L ];
106  for ( int cnt = 0; cnt < L*L*L*L; cnt++ ){ G3dotF[ cnt ] = 0.0; }
107  for ( int cnt = 0; cnt < L*L*L*L; cnt++ ){ lambda2[ cnt ] = 0.0; }
108  for ( int cnt = 0; cnt < L*L*L*L; cnt++ ){ G2multF[ cnt ] = 0.0; }
109  for ( int cnt = 0; cnt < L*L*L*L; cnt++ ){ L2multF[ cnt ] = 0.0; }
110  for ( int cnt = 0; cnt < L*L; cnt++ ){ gamma1[ cnt ] = 0.0; }
111  for ( int cnt = 0; cnt < L*L; cnt++ ){ G1multF[ cnt ] = 0.0; }
112  for ( int cnt = 0; cnt < L*L; cnt++ ){ G2dotF[ cnt ] = 0.0; }
113  for ( int cnt = 0; cnt < L*L; cnt++ ){ L2dotF[ cnt ] = 0.0; }
114  double G1dotF = 0.0;
115 
116  /* Fill gamma1 */
117  for ( int i = 0; i < L; i++ ){
118  for ( int j = i; j < L; j++ ){
119  const double value = the2DM->get1RDM_HAM( i, j );
120  gamma1[ i + L * j ] = value;
121  gamma1[ j + L * i ] = value;
122  }
123  }
124 
125  /* Build G3dotF[i,j,p,q] = sum_[l,s] Gamma3[i,j,l,p,q,s] F[l,s]
126  Build lambda2[i,j,p,q] = Lambda2[i,j,p,q] */
127  for ( int i = 0; i < L; i++ ){
128  for ( int j = i; j < L; j++ ){
129  const int irrep_ij = Irreps::directProd( irreps[ i ], irreps[ j ] );
130  for ( int p = i; p < L; p++ ){
131  for ( int q = i; q < L; q++ ){
132  const int irrep_pq = Irreps::directProd( irreps[ p ], irreps[ q ] );
133  if ( irrep_ij == irrep_pq ){
134  { // sum_[l,s] Gamma3[i,j,l,p,q,s] F[l,s]
135  double value = 0.0;
136  for ( int l = 0; l < L; l++ ){
137  for ( int s = 0; s < L; s++ ){
138  if ( irreps[ l ] == irreps[ s ] ){
139  value += the3DM->get_ham_index(i,j,l,p,q,s) * fock[ l + L * s ];
140  }
141  }
142  }
143  G3dotF[ i + L * ( j + L * ( p + L * q )) ] = value;
144  G3dotF[ j + L * ( i + L * ( q + L * p )) ] = value;
145  G3dotF[ p + L * ( q + L * ( i + L * j )) ] = value;
146  G3dotF[ q + L * ( p + L * ( j + L * i )) ] = value;
147  }
148  { // Lambda2[i,j,p,q]
149  const double val_lambda = the2DM->getTwoDMA_HAM( i, j, p, q )
150  - gamma1[ i + L * p ] * gamma1[ j + L * q ]
151  + gamma1[ i + L * q ] * gamma1[ j + L * p ] * 0.5;
152  lambda2[ i + L * ( j + L * ( p + L * q )) ] = val_lambda;
153  lambda2[ j + L * ( i + L * ( q + L * p )) ] = val_lambda;
154  lambda2[ p + L * ( q + L * ( i + L * j )) ] = val_lambda;
155  lambda2[ q + L * ( p + L * ( j + L * i )) ] = val_lambda;
156  }
157  }
158  }
159  }
160  }
161  }
162 
163  /* Build G1multF[i,j] = sum_[p] Gamma1[i,p] F[p,j]
164  Build G1dotF = sum_[i,j] Gamma1[i,j] F[j,i] */
165  for ( int i = 0; i < L; i++ ){
166  for ( int j = 0; j < L; j++ ){
167  if ( irreps[ i ] == irreps[ j ] ){
168  double value = 0.0;
169  for ( int p = 0; p < L; p++ ){
170  if ( irreps[ j ] == irreps[ p ] ){
171  value += gamma1[ i + L * p ] * fock[ p + L * j ];
172  }
173  }
174  G1multF[ i + L * j ] = value;
175  }
176  }
177  G1dotF += G1multF[ i * ( L + 1 ) ];
178  }
179 
180  /* Build G2dotF[i,j] = sum_[p,q] Gamma2[i,p,j,q] F[p,q]
181  Build L2dotF[i,j] = sum_[p,q] Lambda2[i,p,j,q] F[p,q] */
182  for ( int i = 0; i < L; i++ ){
183  for ( int j = i; j < L; j++ ){
184  if ( irreps[ i ] == irreps[ j ] ){
185  double val_gamma = 0.0;
186  double val_lambda = 0.0;
187  for ( int p = 0; p < L; p++ ){
188  for ( int q = 0; q < L; q++ ){
189  if ( irreps[ p ] == irreps[ q ] ){
190  val_gamma += fock[ p + L * q ] * the2DM->getTwoDMA_HAM( i, p, j, q );
191  val_lambda += fock[ p + L * q ] * lambda2[ i + L * ( p + L * ( j + L * q )) ];
192  }
193  }
194  }
195  G2dotF[ i + L * j ] = val_gamma;
196  G2dotF[ j + L * i ] = val_gamma;
197  L2dotF[ i + L * j ] = val_lambda;
198  L2dotF[ j + L * i ] = val_lambda;
199  }
200  }
201  }
202 
203  /* Build G2multF[i,s,j,k] = sum_[l] Gamma2[i,l,j,k] F[l,s]
204  Build L2multF[i,s,j,k] = sum_[l] Lambda2[i,l,j,k] F[l,s] */
205  for ( int i = 0; i < L; i++ ){
206  for ( int j = 0; j < L; j++ ){
207  const int irrep_ij = Irreps::directProd( irreps[ i ], irreps[ j ] );
208  for ( int k = 0; k < L; k++ ){
209  for ( int s = 0; s < L; s++ ){
210  const int irrep_ks = Irreps::directProd( irreps[ k ], irreps[ s ] );
211  if ( irrep_ij == irrep_ks ){
212  double val_gamma = 0.0;
213  double val_lambda = 0.0;
214  for ( int l = 0; l < L; l++ ){
215  if ( irreps[ l ] == irreps[ s ] ){
216  val_gamma += fock[ l + L * s ] * the2DM->getTwoDMA_HAM( i, l, j, k );
217  val_lambda += fock[ l + L * s ] * lambda2[ i + L * ( l + L * ( j + L * k ) ) ];
218  }
219  }
220  G2multF[ i + L * ( s + L * ( j + L * k ) ) ] = val_gamma;
221  L2multF[ i + L * ( s + L * ( j + L * k ) ) ] = val_lambda;
222  }
223  }
224  }
225  }
226  }
227 
228  /* Fill result */
229  #pragma omp parallel for schedule(dynamic)
230  for ( int i = 0; i < L; i++ ){
231  for ( int j = i; j < L; j++ ){
232  for ( int k = i; k < L; k++ ){
233  const int irrep_ijk = Irreps::directProd( Irreps::directProd( irreps[ i ], irreps[ j ] ), irreps[ k ] );
234  for ( int p = i; p < L; p++ ){
235  for ( int q = i; q < L; q++ ){
236  for ( int r = q; r < L; r++ ){
237  const int irrep_pqr = Irreps::directProd( Irreps::directProd( irreps[ p ], irreps[ q ] ), irreps[ r ] );
238  if ( irrep_ijk == irrep_pqr ){
239 
240  double dm3_contribution = 0.0;
241  double gamma2_part1 = 0.0;
242  double gamma2_part2 = 0.0;
243  double lambda2_part1 = 0.0;
244  double lambda2_part2 = 0.0;
245 
246  for ( int ls = 0; ls < L; ls++ ){
247 
248  dm3_contribution += ( the3DM->get_ham_index( ls, j, k, p, q, r ) * G1multF[ i + L * ls ]
249  + the3DM->get_ham_index( i, ls, k, p, q, r ) * G1multF[ j + L * ls ]
250  + the3DM->get_ham_index( i, j, ls, p, q, r ) * G1multF[ k + L * ls ]
251  + the3DM->get_ham_index( i, j, k, ls, q, r ) * G1multF[ p + L * ls ]
252  + the3DM->get_ham_index( i, j, k, p, ls, r ) * G1multF[ q + L * ls ]
253  + the3DM->get_ham_index( i, j, k, p, q, ls ) * G1multF[ r + L * ls ] );
254 
255  gamma2_part1 += ( the2DM->getTwoDMA_HAM( i, j, p, ls ) * G2multF[ k + L * ( ls + L * ( r + L * q )) ]
256  + the2DM->getTwoDMA_HAM( i, j, ls, q ) * G2multF[ k + L * ( ls + L * ( r + L * p )) ]
257  + the2DM->getTwoDMA_HAM( i, k, p, ls ) * G2multF[ j + L * ( ls + L * ( q + L * r )) ]
258  + the2DM->getTwoDMA_HAM( i, k, ls, r ) * G2multF[ j + L * ( ls + L * ( q + L * p )) ]
259  + the2DM->getTwoDMA_HAM( k, j, ls, q ) * G2multF[ i + L * ( ls + L * ( p + L * r )) ]
260  + the2DM->getTwoDMA_HAM( k, j, r, ls ) * G2multF[ i + L * ( ls + L * ( p + L * q )) ] );
261 
262  gamma2_part2 += ( the2DM->getTwoDMA_HAM( i, j, r, ls ) * ( G2multF[ k + L * ( ls + L * ( p + L * q )) ]
263  + 0.5 * G2multF[ k + L * ( ls + L * ( q + L * p )) ] )
264  + the2DM->getTwoDMA_HAM( i, j, ls, r ) * ( G2multF[ k + L * ( ls + L * ( q + L * p )) ]
265  + 0.5 * G2multF[ k + L * ( ls + L * ( p + L * q )) ] )
266  + the2DM->getTwoDMA_HAM( i, k, q, ls ) * ( G2multF[ j + L * ( ls + L * ( p + L * r )) ]
267  + 0.5 * G2multF[ j + L * ( ls + L * ( r + L * p )) ] )
268  + the2DM->getTwoDMA_HAM( i, k, ls, q ) * ( G2multF[ j + L * ( ls + L * ( r + L * p )) ]
269  + 0.5 * G2multF[ j + L * ( ls + L * ( p + L * r )) ] )
270  + the2DM->getTwoDMA_HAM( k, j, p, ls ) * ( G2multF[ i + L * ( ls + L * ( r + L * q )) ]
271  + 0.5 * G2multF[ i + L * ( ls + L * ( q + L * r )) ] )
272  + the2DM->getTwoDMA_HAM( k, j, ls, p ) * ( G2multF[ i + L * ( ls + L * ( q + L * r )) ]
273  + 0.5 * G2multF[ i + L * ( ls + L * ( r + L * q )) ] ) );
274 
275  lambda2_part1 += ( lambda2[ i + L * ( j + L * ( p + L * ls )) ] * L2multF[ k + L * ( ls + L * ( r + L * q )) ]
276  + lambda2[ i + L * ( j + L * ( ls + L * q )) ] * L2multF[ k + L * ( ls + L * ( r + L * p )) ]
277  + lambda2[ i + L * ( k + L * ( p + L * ls )) ] * L2multF[ j + L * ( ls + L * ( q + L * r )) ]
278  + lambda2[ i + L * ( k + L * ( ls + L * r )) ] * L2multF[ j + L * ( ls + L * ( q + L * p )) ]
279  + lambda2[ k + L * ( j + L * ( ls + L * q )) ] * L2multF[ i + L * ( ls + L * ( p + L * r )) ]
280  + lambda2[ k + L * ( j + L * ( r + L * ls )) ] * L2multF[ i + L * ( ls + L * ( p + L * q )) ] );
281 
282  lambda2_part2 += ( lambda2[ i + L * ( j + L * ( r + L * ls )) ] * ( L2multF[ k + L * ( ls + L * ( p + L * q )) ]
283  + 0.5 * L2multF[ k + L * ( ls + L * ( q + L * p )) ] )
284  + lambda2[ i + L * ( j + L * ( ls + L * r )) ] * ( L2multF[ k + L * ( ls + L * ( q + L * p )) ]
285  + 0.5 * L2multF[ k + L * ( ls + L * ( p + L * q )) ] )
286  + lambda2[ i + L * ( k + L * ( q + L * ls )) ] * ( L2multF[ j + L * ( ls + L * ( p + L * r )) ]
287  + 0.5 * L2multF[ j + L * ( ls + L * ( r + L * p )) ] )
288  + lambda2[ i + L * ( k + L * ( ls + L * q )) ] * ( L2multF[ j + L * ( ls + L * ( r + L * p )) ]
289  + 0.5 * L2multF[ j + L * ( ls + L * ( p + L * r )) ] )
290  + lambda2[ k + L * ( j + L * ( p + L * ls )) ] * ( L2multF[ i + L * ( ls + L * ( r + L * q )) ]
291  + 0.5 * L2multF[ i + L * ( ls + L * ( q + L * r )) ] )
292  + lambda2[ k + L * ( j + L * ( ls + L * p )) ] * ( L2multF[ i + L * ( ls + L * ( q + L * r )) ]
293  + 0.5 * L2multF[ i + L * ( ls + L * ( r + L * q )) ] ) );
294  }
295 
296  const double contracted_value = ( the3DM->get_ham_index( i, j, k, p, q, r ) * G1dotF
297 
298  + G3dotF[ i + L * ( j + L * ( p + L * q )) ] * gamma1[ k + L * r ]
299  - 0.5 * G3dotF[ i + L * ( j + L * ( r + L * q )) ] * gamma1[ k + L * p ]
300  - 0.5 * G3dotF[ i + L * ( j + L * ( p + L * r )) ] * gamma1[ k + L * q ]
301  + G3dotF[ i + L * ( k + L * ( p + L * r )) ] * gamma1[ j + L * q ]
302  - 0.5 * G3dotF[ i + L * ( k + L * ( q + L * r )) ] * gamma1[ j + L * p ]
303  - 0.5 * G3dotF[ i + L * ( k + L * ( p + L * q )) ] * gamma1[ j + L * r ]
304  + G3dotF[ j + L * ( k + L * ( q + L * r )) ] * gamma1[ i + L * p ]
305  - 0.5 * G3dotF[ j + L * ( k + L * ( p + L * r )) ] * gamma1[ i + L * q ]
306  - 0.5 * G3dotF[ j + L * ( k + L * ( q + L * p )) ] * gamma1[ i + L * r ]
307 
308  - 0.5 * dm3_contribution
309 
310  - the2DM->getTwoDMA_HAM( i, j, p, q ) * G2dotF[ k + L * r ]
311  + 0.5 * the2DM->getTwoDMA_HAM( i, j, p, r ) * G2dotF[ k + L * q ]
312  + 0.5 * the2DM->getTwoDMA_HAM( i, j, r, q ) * G2dotF[ k + L * p ]
313  - the2DM->getTwoDMA_HAM( i, k, p, r ) * G2dotF[ j + L * q ]
314  + 0.5 * the2DM->getTwoDMA_HAM( i, k, p, q ) * G2dotF[ j + L * r ]
315  + 0.5 * the2DM->getTwoDMA_HAM( i, k, q, r ) * G2dotF[ j + L * p ]
316  - the2DM->getTwoDMA_HAM( k, j, r, q ) * G2dotF[ i + L * p ]
317  + 0.5 * the2DM->getTwoDMA_HAM( k, j, p, q ) * G2dotF[ i + L * r ]
318  + 0.5 * the2DM->getTwoDMA_HAM( k, j, r, p ) * G2dotF[ i + L * q ]
319 
320  + 0.5 * gamma2_part1
321  - gamma2_part2 / 3.0
322 
323  + 2 * lambda2[ i + L * ( j + L * ( p + L * q )) ] * L2dotF[ k + L * r ]
324  - lambda2[ i + L * ( j + L * ( p + L * r )) ] * L2dotF[ k + L * q ]
325  - lambda2[ i + L * ( j + L * ( r + L * q )) ] * L2dotF[ k + L * p ]
326  + 2 * lambda2[ i + L * ( k + L * ( p + L * r )) ] * L2dotF[ j + L * q ]
327  - lambda2[ i + L * ( k + L * ( p + L * q )) ] * L2dotF[ j + L * r ]
328  - lambda2[ i + L * ( k + L * ( q + L * r )) ] * L2dotF[ j + L * p ]
329  + 2 * lambda2[ k + L * ( j + L * ( r + L * q )) ] * L2dotF[ i + L * p ]
330  - lambda2[ k + L * ( j + L * ( p + L * q )) ] * L2dotF[ i + L * r ]
331  - lambda2[ k + L * ( j + L * ( r + L * p )) ] * L2dotF[ i + L * q ]
332 
333  - lambda2_part1
334  + lambda2_part2 / 1.5 );
335 
336  result[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] = contracted_value;
337  result[ i + L * ( k + L * ( j + L * ( p + L * ( r + L * q )))) ] = contracted_value;
338  result[ j + L * ( i + L * ( k + L * ( q + L * ( p + L * r )))) ] = contracted_value;
339  result[ k + L * ( i + L * ( j + L * ( r + L * ( p + L * q )))) ] = contracted_value;
340  result[ j + L * ( k + L * ( i + L * ( q + L * ( r + L * p )))) ] = contracted_value;
341  result[ k + L * ( j + L * ( i + L * ( r + L * ( q + L * p )))) ] = contracted_value;
342  result[ p + L * ( q + L * ( r + L * ( i + L * ( j + L * k )))) ] = contracted_value;
343  result[ p + L * ( r + L * ( q + L * ( i + L * ( k + L * j )))) ] = contracted_value;
344  result[ q + L * ( p + L * ( r + L * ( j + L * ( i + L * k )))) ] = contracted_value;
345  result[ r + L * ( p + L * ( q + L * ( k + L * ( i + L * j )))) ] = contracted_value;
346  result[ q + L * ( r + L * ( p + L * ( j + L * ( k + L * i )))) ] = contracted_value;
347  result[ r + L * ( q + L * ( p + L * ( k + L * ( j + L * i )))) ] = contracted_value;
348  }
349  }
350  }
351  }
352  }
353  }
354  }
355 
356  delete [] G3dotF;
357  delete [] lambda2;
358  delete [] G2multF;
359  delete [] L2multF;
360  delete [] gamma1;
361  delete [] G1multF;
362  delete [] G2dotF;
363  delete [] L2dotF;
364  delete [] irreps;
365 
366  gettimeofday(&end, NULL);
367  const double elapsed = (end.tv_sec - start.tv_sec) + 1e-6 * (end.tv_usec - start.tv_usec);
368  std::cout << "Cumulant :: Contraction of cu(4)-4RDM with CASPT2 Fock operator took " << elapsed << " seconds." << std::endl;
369 
370 }
371 
372 double CheMPS2::Cumulant::lambda2_ham(const TwoDM * the2DM, const int i, const int j, const int p, const int q){
373 
374  const double value = the2DM->getTwoDMA_HAM( i, j, p, q )
375  - the2DM->get1RDM_HAM( i, p ) * the2DM->get1RDM_HAM( j, q )
376  + the2DM->get1RDM_HAM( i, q ) * the2DM->get1RDM_HAM( j, p ) * 0.5;
377  return value;
378 
379 }
380 
381 double CheMPS2::Cumulant::gamma4_ham(const Problem * prob, const ThreeDM * the3DM, const TwoDM * the2DM, const int i, const int j, const int k, const int l,
382  const int p, const int q, const int r, const int s){
383 
384  //Prob assumes you use DMRG orbs... f1 converts HAM orbs to DMRG orbs
385  const int irrep_i = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( i ) : i );
386  const int irrep_j = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( j ) : j );
387  const int irrep_k = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( k ) : k );
388  const int irrep_l = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( l ) : l );
389  const int irrep_p = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( p ) : p );
390  const int irrep_q = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( q ) : q );
391  const int irrep_r = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( r ) : r );
392  const int irrep_s = prob->gIrrep(( prob->gReorder() ) ? prob->gf1( s ) : s );
393 
394  const int irrep_ij = Irreps::directProd( irrep_i, irrep_j );
395  const int irrep_kl = Irreps::directProd( irrep_k, irrep_l );
396  const int irrep_pq = Irreps::directProd( irrep_p, irrep_q );
397  const int irrep_rs = Irreps::directProd( irrep_r, irrep_s );
398 
399  if ( Irreps::directProd( irrep_ij, irrep_kl ) != Irreps::directProd( irrep_pq, irrep_rs ) ){ return 0.0; }
400 
401  const double part1 = ( the3DM->get_ham_index(i,j,k,p,q,r) * the2DM->get1RDM_HAM(l,s)
402  - 0.5 * the3DM->get_ham_index(i,j,k,s,q,r) * the2DM->get1RDM_HAM(l,p)
403  - 0.5 * the3DM->get_ham_index(i,j,k,p,s,r) * the2DM->get1RDM_HAM(l,q)
404  - 0.5 * the3DM->get_ham_index(i,j,k,p,q,s) * the2DM->get1RDM_HAM(l,r)
405  + the3DM->get_ham_index(i,j,l,p,q,s) * the2DM->get1RDM_HAM(k,r)
406  - 0.5 * the3DM->get_ham_index(i,j,l,r,q,s) * the2DM->get1RDM_HAM(k,p)
407  - 0.5 * the3DM->get_ham_index(i,j,l,p,r,s) * the2DM->get1RDM_HAM(k,q)
408  - 0.5 * the3DM->get_ham_index(i,j,l,p,q,r) * the2DM->get1RDM_HAM(k,s)
409  + the3DM->get_ham_index(i,k,l,p,r,s) * the2DM->get1RDM_HAM(j,q)
410  - 0.5 * the3DM->get_ham_index(i,k,l,q,r,s) * the2DM->get1RDM_HAM(j,p)
411  - 0.5 * the3DM->get_ham_index(i,k,l,p,q,s) * the2DM->get1RDM_HAM(j,r)
412  - 0.5 * the3DM->get_ham_index(i,k,l,p,r,q) * the2DM->get1RDM_HAM(j,s)
413  + the3DM->get_ham_index(j,k,l,q,r,s) * the2DM->get1RDM_HAM(i,p)
414  - 0.5 * the3DM->get_ham_index(j,k,l,p,r,s) * the2DM->get1RDM_HAM(i,q)
415  - 0.5 * the3DM->get_ham_index(j,k,l,q,p,s) * the2DM->get1RDM_HAM(i,r)
416  - 0.5 * the3DM->get_ham_index(j,k,l,q,r,p) * the2DM->get1RDM_HAM(i,s) );
417 
418  const double part2 = ( the2DM->getTwoDMA_HAM(i,j,p,q) * the2DM->getTwoDMA_HAM(k,l,r,s)
419  - the2DM->getTwoDMA_HAM(i,j,p,r) * the2DM->getTwoDMA_HAM(k,l,q,s) * 0.5
420  - the2DM->getTwoDMA_HAM(i,j,p,s) * the2DM->getTwoDMA_HAM(k,l,r,q) * 0.5
421  - the2DM->getTwoDMA_HAM(i,j,r,q) * the2DM->getTwoDMA_HAM(k,l,p,s) * 0.5
422  - the2DM->getTwoDMA_HAM(i,j,s,q) * the2DM->getTwoDMA_HAM(k,l,r,p) * 0.5
423  + the2DM->getTwoDMA_HAM(i,j,r,s) * the2DM->getTwoDMA_HAM(k,l,p,q) / 3.0
424  + the2DM->getTwoDMA_HAM(i,j,r,s) * the2DM->getTwoDMA_HAM(k,l,q,p) / 6.0
425  + the2DM->getTwoDMA_HAM(i,j,s,r) * the2DM->getTwoDMA_HAM(k,l,p,q) / 6.0
426  + the2DM->getTwoDMA_HAM(i,j,s,r) * the2DM->getTwoDMA_HAM(k,l,q,p) / 3.0
427  + the2DM->getTwoDMA_HAM(i,k,p,r) * the2DM->getTwoDMA_HAM(j,l,q,s)
428  - the2DM->getTwoDMA_HAM(i,k,p,q) * the2DM->getTwoDMA_HAM(j,l,r,s) * 0.5
429  - the2DM->getTwoDMA_HAM(i,k,p,s) * the2DM->getTwoDMA_HAM(j,l,q,r) * 0.5
430  - the2DM->getTwoDMA_HAM(i,k,q,r) * the2DM->getTwoDMA_HAM(j,l,p,s) * 0.5
431  - the2DM->getTwoDMA_HAM(i,k,s,r) * the2DM->getTwoDMA_HAM(j,l,q,p) * 0.5
432  + the2DM->getTwoDMA_HAM(i,k,q,s) * the2DM->getTwoDMA_HAM(j,l,p,r) / 3.0
433  + the2DM->getTwoDMA_HAM(i,k,s,q) * the2DM->getTwoDMA_HAM(j,l,p,r) / 6.0
434  + the2DM->getTwoDMA_HAM(i,k,q,s) * the2DM->getTwoDMA_HAM(j,l,r,p) / 6.0
435  + the2DM->getTwoDMA_HAM(i,k,s,q) * the2DM->getTwoDMA_HAM(j,l,r,p) / 3.0
436  + the2DM->getTwoDMA_HAM(i,l,p,s) * the2DM->getTwoDMA_HAM(k,j,r,q)
437  - the2DM->getTwoDMA_HAM(i,l,p,r) * the2DM->getTwoDMA_HAM(k,j,s,q) * 0.5
438  - the2DM->getTwoDMA_HAM(i,l,p,q) * the2DM->getTwoDMA_HAM(k,j,r,s) * 0.5
439  - the2DM->getTwoDMA_HAM(i,l,r,s) * the2DM->getTwoDMA_HAM(k,j,p,q) * 0.5
440  - the2DM->getTwoDMA_HAM(i,l,q,s) * the2DM->getTwoDMA_HAM(k,j,r,p) * 0.5
441  + the2DM->getTwoDMA_HAM(i,l,r,q) * the2DM->getTwoDMA_HAM(k,j,p,s) / 3.0
442  + the2DM->getTwoDMA_HAM(i,l,q,r) * the2DM->getTwoDMA_HAM(k,j,p,s) / 6.0
443  + the2DM->getTwoDMA_HAM(i,l,r,q) * the2DM->getTwoDMA_HAM(k,j,s,p) / 6.0
444  + the2DM->getTwoDMA_HAM(i,l,q,r) * the2DM->getTwoDMA_HAM(k,j,s,p) / 3.0 );
445 
446  const double part3 = ( lambda2_ham(the2DM,i,j,p,q) * lambda2_ham(the2DM,k,l,r,s)
447  - lambda2_ham(the2DM,i,j,p,r) * lambda2_ham(the2DM,k,l,q,s) * 0.5
448  - lambda2_ham(the2DM,i,j,p,s) * lambda2_ham(the2DM,k,l,r,q) * 0.5
449  - lambda2_ham(the2DM,i,j,r,q) * lambda2_ham(the2DM,k,l,p,s) * 0.5
450  - lambda2_ham(the2DM,i,j,s,q) * lambda2_ham(the2DM,k,l,r,p) * 0.5
451  + lambda2_ham(the2DM,i,j,r,s) * lambda2_ham(the2DM,k,l,p,q) / 3.0
452  + lambda2_ham(the2DM,i,j,r,s) * lambda2_ham(the2DM,k,l,q,p) / 6.0
453  + lambda2_ham(the2DM,i,j,s,r) * lambda2_ham(the2DM,k,l,p,q) / 6.0
454  + lambda2_ham(the2DM,i,j,s,r) * lambda2_ham(the2DM,k,l,q,p) / 3.0
455  + lambda2_ham(the2DM,i,k,p,r) * lambda2_ham(the2DM,j,l,q,s)
456  - lambda2_ham(the2DM,i,k,p,q) * lambda2_ham(the2DM,j,l,r,s) * 0.5
457  - lambda2_ham(the2DM,i,k,p,s) * lambda2_ham(the2DM,j,l,q,r) * 0.5
458  - lambda2_ham(the2DM,i,k,q,r) * lambda2_ham(the2DM,j,l,p,s) * 0.5
459  - lambda2_ham(the2DM,i,k,s,r) * lambda2_ham(the2DM,j,l,q,p) * 0.5
460  + lambda2_ham(the2DM,i,k,q,s) * lambda2_ham(the2DM,j,l,p,r) / 3.0
461  + lambda2_ham(the2DM,i,k,s,q) * lambda2_ham(the2DM,j,l,p,r) / 6.0
462  + lambda2_ham(the2DM,i,k,q,s) * lambda2_ham(the2DM,j,l,r,p) / 6.0
463  + lambda2_ham(the2DM,i,k,s,q) * lambda2_ham(the2DM,j,l,r,p) / 3.0
464  + lambda2_ham(the2DM,i,l,p,s) * lambda2_ham(the2DM,k,j,r,q)
465  - lambda2_ham(the2DM,i,l,p,r) * lambda2_ham(the2DM,k,j,s,q) * 0.5
466  - lambda2_ham(the2DM,i,l,p,q) * lambda2_ham(the2DM,k,j,r,s) * 0.5
467  - lambda2_ham(the2DM,i,l,r,s) * lambda2_ham(the2DM,k,j,p,q) * 0.5
468  - lambda2_ham(the2DM,i,l,q,s) * lambda2_ham(the2DM,k,j,r,p) * 0.5
469  + lambda2_ham(the2DM,i,l,r,q) * lambda2_ham(the2DM,k,j,p,s) / 3.0
470  + lambda2_ham(the2DM,i,l,q,r) * lambda2_ham(the2DM,k,j,p,s) / 6.0
471  + lambda2_ham(the2DM,i,l,r,q) * lambda2_ham(the2DM,k,j,s,p) / 6.0
472  + lambda2_ham(the2DM,i,l,q,r) * lambda2_ham(the2DM,k,j,s,p) / 3.0 );
473 
474  return ( part1 - part2 + 2 * part3 );
475 
476 }
477 
478 
479 
480 
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::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::gIrrep
int gIrrep(const int nOrb) const
Get an orbital irrep.
Definition: Problem.cpp:336
CheMPS2::Problem::gf1
int gf1(const int HamOrb) const
Get the DMRG index corresponding to a Ham index.
Definition: Problem.cpp:347
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::Cumulant::gamma4_fock_contract_ham
static void gamma4_fock_contract_ham(const Problem *prob, const ThreeDM *the3DM, const TwoDM *the2DM, double *fock, double *result)
Contract the CASPT2 Fock operator with the cumulant approximation of in time, using HAM indices...
Definition: Cumulant.cpp:84
CheMPS2::ThreeDM::get_ham_index
double get_ham_index(const int cnt1, const int cnt2, const int cnt3, const int cnt4, const int cnt5, const int cnt6) const
Get a 3-RDM, using the HAM indices.
Definition: ThreeDM.cpp:148
CheMPS2::Cumulant::gamma4_ham
static double gamma4_ham(const Problem *prob, const ThreeDM *the3DM, const TwoDM *the2DM, const int i, const int j, const int k, const int l, const int p, const int q, const int r, const int s)
Get the cumulant approximation of , using HAM indices.
Definition: Cumulant.cpp:381