← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/psimrcc/mrcc_t2_amps.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/psimrcc/mrcc_t2_amps.cc
 
1
/***************************************************************************
 
2
 *  PSIMRCC : Copyright (C) 2007 by Francesco Evangelista and Andrew Simmonett
 
3
 *  frank@ccc.uga.edu   andysim@ccc.uga.edu
 
4
 *  A multireference coupled cluster code
 
5
 ***************************************************************************/
 
6
#include <libmoinfo/libmoinfo.h>
 
7
#include "mrcc.h"
 
8
#include "matrix.h"
 
9
#include "blas.h"
 
10
#include "debugging.h"
 
11
#include <libutil/libutil.h>
 
12
#include "algebra_interface.h"
 
13
#include "memory_manager.h"
 
14
 
 
15
extern FILE* outfile;
 
16
 
 
17
namespace psi{ namespace psimrcc{
 
18
 
 
19
using namespace std;
 
20
 
 
21
void CCMRCC::build_t2_amplitudes()
 
22
{
 
23
  build_t2_iJaB_amplitudes();
 
24
  build_t2_ijab_amplitudes();
 
25
  build_t2_IJAB_amplitudes();
 
26
 
 
27
}
 
28
 
 
29
void CCMRCC::build_t2_ijab_amplitudes()
 
30
{
 
31
  Timer timer;
 
32
  DEBUGGING(1,
 
33
    fprintf(outfile,"\n\tBuilding the t2_ijab Amplitudes   ...");
 
34
    fflush(outfile);
 
35
  )
 
36
  if(moinfo->get_ref_size("o")==0){
 
37
    blas->append("t2_eqns[oo][vv]{c}  = t2_eqns[oO][vV]{c}");
 
38
    blas->append("t2_eqns[oo][vv]{c} += #2134# - t2_eqns[oO][vV]{c}");
 
39
  }else{
 
40
    // Closed-shell
 
41
    blas->append("t2_eqns[oo][vv]{c}  = <[oo]:[vv]>");
 
42
 
 
43
    blas->append("t2_eqns[oo][vv]{c} += #3124# - t2[v][voo]{c} 1@2 F'_ae[v][v]{c}");
 
44
    blas->append("t2_eqns[oo][vv]{c} += #4123#   t2[v][voo]{c} 1@2 F'_ae[v][v]{c}");
 
45
 
 
46
    blas->append("t2_eqns[oo][vv]{c} += #1342#   t2[o][ovv]{c} 1@1 F'_mi[o][o]{c}");
 
47
    blas->append("t2_eqns[oo][vv]{c} += #2341# - t2[o][ovv]{c} 1@1 F'_mi[o][o]{c}");
 
48
 
 
49
    blas->append("t2_eqns[oo][vv]{c} += 1/2  W_mnij[oo][oo]{c} 1@1 tau[oo][vv]{c}");
 
50
 
 
51
    blas->append("t2_eqns[oo][v>v]{c} = tau[oo][v>v]{c} 2@2 <[v>v]:[v>v]>");
 
52
 
 
53
    blas->append("t2_eqns[oo][vv]{c} +>= #1234# t2_eqns[oo][v>v]{c}");
 
54
 
 
55
    blas->append("t2_eqns[oo][vv]{c} += #1234# - Z_ijam[oov][o]{c} 2@1 t1[o][v]{c}");
 
56
    blas->append("t2_eqns[oo][vv]{c} += #1243#   Z_ijam[oov][o]{c} 2@1 t1[o][v]{c}");
 
57
 
 
58
    blas->append("t2_eqns[oo][vv]{c} += #2413#   W_jbme[ov][ov]{c} 2@2 t2[ov][ov]{c}");
 
59
    blas->append("t2_eqns[oo][vv]{c} += #2314# - W_jbme[ov][ov]{c} 2@2 t2[ov][ov]{c}");
 
60
 
 
61
    blas->append("t2_eqns[oo][vv]{c} += #1423# - W_jbme[ov][ov]{c} 2@2 t2[ov][ov]{c}");
 
62
    blas->append("t2_eqns[oo][vv]{c} += #1324#   W_jbme[ov][ov]{c} 2@2 t2[ov][ov]{c}");
 
63
 
 
64
//  blas->append("t2_eqns[oo][vv]{c} += #P-(34)P-(12)4213# - ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
65
 
 
66
    blas->append("t2_eqns[oo][vv]{c} += #4213# - ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
67
    blas->append("t2_eqns[oo][vv]{c} += #3214# + ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
68
    blas->append("t2_eqns[oo][vv]{c} += #4123# + ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
69
    blas->append("t2_eqns[oo][vv]{c} += #3124# - ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
70
 
 
71
    blas->append("t2_eqns[oo][vv]{c} += #2413#   W_jbME[ov][OV]{c} 2@2 t2[ov][OV]{c}");
 
72
    blas->append("t2_eqns[oo][vv]{c} += #2314# - W_jbME[ov][OV]{c} 2@2 t2[ov][OV]{c}");
 
73
 
 
74
    blas->append("t2_eqns[oo][vv]{c} += #1423# - W_jbME[ov][OV]{c} 2@2 t2[ov][OV]{c}");
 
75
    blas->append("t2_eqns[oo][vv]{c} += #1324#   W_jbME[ov][OV]{c} 2@2 t2[ov][OV]{c}");
 
76
 
 
77
    blas->append("t2_eqns[oo][vv]{c} += #1234#   t1[o][v]{c} 2@1 <[v]:[ovv]>");
 
78
    blas->append("t2_eqns[oo][vv]{c} += #2134# - t1[o][v]{c} 2@1 <[v]:[ovv]>");
 
79
 
 
80
    blas->append("t2_eqns[oo][vv]{c} += #3412# - t1[o][v]{c} 1@1 <[o]:[voo]>");
 
81
    blas->append("t2_eqns[oo][vv]{c} += #4312#   t1[o][v]{c} 1@1 <[o]:[voo]>");
 
82
  }
 
83
 
 
84
  // Open-shell
 
85
  blas->append("t2_eqns[oo][vv]{o}  = <[oo]:[vv]>");
 
86
 
 
87
  blas->append("t2_eqns[oo][vv]{o} += #3124# - t2[v][voo]{o} 1@2 F'_ae[v][v]{o}");
 
88
  blas->append("t2_eqns[oo][vv]{o} += #4123#   t2[v][voo]{o} 1@2 F'_ae[v][v]{o}");
 
89
 
 
90
  blas->append("t2_eqns[oo][vv]{o} += #1342#   t2[o][ovv]{o} 1@1 F'_mi[o][o]{o}");
 
91
  blas->append("t2_eqns[oo][vv]{o} += #2341# - t2[o][ovv]{o} 1@1 F'_mi[o][o]{o}");
 
92
 
 
93
  blas->append("t2_eqns[oo][vv]{o} += 1/2  W_mnij[oo][oo]{o} 1@1 tau[oo][vv]{o}");
 
94
 
 
95
  blas->append("t2_eqns[oo][v>v]{o} = tau[oo][v>v]{o} 2@2 <[v>v]:[v>v]>");
 
96
 
 
97
  blas->append("t2_eqns[oo][vv]{o} +>= #1234# t2_eqns[oo][v>v]{o}");
 
98
 
 
99
  blas->append("t2_eqns[oo][vv]{o} += #1234# - Z_ijam[oov][o]{o} 2@1 t1[o][v]{o}");
 
100
  blas->append("t2_eqns[oo][vv]{o} += #1243#   Z_ijam[oov][o]{o} 2@1 t1[o][v]{o}");
 
101
 
 
102
  blas->append("t2_eqns[oo][vv]{o} += #2413#   W_jbme[ov][ov]{o} 2@2 t2[ov][ov]{o}");
 
103
  blas->append("t2_eqns[oo][vv]{o} += #2314# - W_jbme[ov][ov]{o} 2@2 t2[ov][ov]{o}");
 
104
 
 
105
  blas->append("t2_eqns[oo][vv]{o} += #1423# - W_jbme[ov][ov]{o} 2@2 t2[ov][ov]{o}");
 
106
  blas->append("t2_eqns[oo][vv]{o} += #1324#   W_jbme[ov][ov]{o} 2@2 t2[ov][ov]{o}");
 
107
 
 
108
  blas->append("t2_eqns[oo][vv]{o} += #4213# - ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{o}");
 
109
  blas->append("t2_eqns[oo][vv]{o} += #3214# + ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{o}");
 
110
  blas->append("t2_eqns[oo][vv]{o} += #4123# + ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{o}");
 
111
  blas->append("t2_eqns[oo][vv]{o} += #3124# - ([ov]:[vo]) 1@2 t1t1_iame[ov][ov]{o}");
 
112
 
 
113
 
 
114
  blas->append("t2_eqns[oo][vv]{o} += #2413#   W_jbME[ov][OV]{o} 2@2 t2[ov][OV]{o}");
 
115
  blas->append("t2_eqns[oo][vv]{o} += #2314# - W_jbME[ov][OV]{o} 2@2 t2[ov][OV]{o}");
 
116
 
 
117
  blas->append("t2_eqns[oo][vv]{o} += #1423# - W_jbME[ov][OV]{o} 2@2 t2[ov][OV]{o}");
 
118
  blas->append("t2_eqns[oo][vv]{o} += #1324#   W_jbME[ov][OV]{o} 2@2 t2[ov][OV]{o}");
 
119
 
 
120
  blas->append("t2_eqns[oo][vv]{o} += #1234#   t1[o][v]{o} 2@1 <[v]:[ovv]>");
 
121
  blas->append("t2_eqns[oo][vv]{o} += #2134# - t1[o][v]{o} 2@1 <[v]:[ovv]>");
 
122
 
 
123
  blas->append("t2_eqns[oo][vv]{o} += #3412# - t1[o][v]{o} 1@1 <[o]:[voo]>");
 
124
  blas->append("t2_eqns[oo][vv]{o} += #4312#   t1[o][v]{o} 1@1 <[o]:[voo]>");
 
125
 
 
126
  DEBUGGING(3,
 
127
    blas->print("t2_eqns[oo][vv]{c}");
 
128
    blas->print("t2_eqns[oo][vv]{o}");
 
129
  );
 
130
 
 
131
  DEBUGGING(1,
 
132
    fprintf(outfile," done. Timing %20.6f s",timer.get());
 
133
    fflush(outfile);
 
134
  );
 
135
}
 
136
 
 
137
void CCMRCC::build_t2_iJaB_amplitudes()
 
138
{
 
139
  Timer timer;
 
140
  DEBUGGING(1,
 
141
    fprintf(outfile,"\n\tBuilding the t2_iJaB Amplitudes   ...");
 
142
    fflush(outfile);
 
143
  );
 
144
  // Closed-shell
 
145
  blas->append("t2_eqns[oO][vV]{c}  = <[oo]|[vv]>");
 
146
 
 
147
  blas->append("t2_eqns[oO][vV]{c} += #3214# t2[V][vOo]{c} 1@2 F'_ae[v][v]{c}");
 
148
  blas->append("t2_eqns[oO][vV]{c} += #4123# t2[v][VoO]{c} 1@2 F'_ae[v][v]{c}");
 
149
 
 
150
  blas->append("t2_eqns[oO][vV]{c} += #1432# - t2[O][oVv]{c} 1@1 F'_mi[o][o]{c}");
 
151
  blas->append("t2_eqns[oO][vV]{c} += #2341# - t2[o][OvV]{c} 1@1 F'_mi[o][o]{c}");
 
152
 
 
153
  blas->append("t2_eqns[oO][vV]{c} += W_mNiJ[oO][oO]{c} 1@1 tau[oO][vV]{c}");
 
154
 
 
155
  //blas->append("t2_eqns[oO][vV]{c} += tau[oO][vV]{c} 2@2 <[vv]|[vv]>");
 
156
 
 
157
  blas->append("t2_eqns[oO][vV]{c} += tau[oO][v>=V]{c} 2@2 <[vv]|[v>=v]>");
 
158
  blas->append("t2_eqns[oO][vV]{c} += #1243# tau[oO][V>=v]{c} 2@2 <[vv]|[v>=v]>");
 
159
 
 
160
 
 
161
  blas->append("t2_eqns[oO][vV]{c} += #1234#  - Z_iJaM[oOv][O]{c} 2@1 t1[O][V]{c}");
 
162
  blas->append("t2_eqns[oO][vV]{c} += #1243#    Z_iJAm[oOV][o]{c} 2@1 t1[o][v]{c}");
 
163
 
 
164
  blas->append("t2_eqns[oO][vV]{c} += #2413#   W_jbME[ov][OV]{c} 2@2 t2[ov][ov]{c}");
 
165
  blas->append("t2_eqns[oO][vV]{c} += #2413#   W_jbme[ov][ov]{c} 2@2 t2[ov][OV]{c}");
 
166
 
 
167
  blas->append("t2_eqns[oO][vV]{c} += #2314#   W_jBmE[oV][oV]{c} 2@2 t2[oV][Ov]{c}");
 
168
  blas->append("t2_eqns[oO][vV]{c} += #1423#   W_jBmE[oV][oV]{c} 2@1 t2[oV][Ov]{c}");
 
169
 
 
170
  blas->append("t2_eqns[oO][vV]{c} += #1324#   W_jbME[ov][OV]{c} 2@2 t2[OV][OV]{c}");
 
171
 
 
172
  blas->append("t2_eqns[oO][vV]{c} += #1324#   W_jbme[ov][ov]{c} 2@1 t2[ov][OV]{c}");
 
173
 
 
174
 
 
175
  blas->append("t2_eqns[oO][vV]{c} += #4213# - ([ov]|[vo]) 1@2 t1t1_iame[ov][ov]{c}");
 
176
 
 
177
  blas->append("t2_eqns[oO][vV]{c} += #2314# - <[ov]|[ov]> 1@2 t1t1_iAMe[oV][Ov]{c}");
 
178
 
 
179
  blas->append("t2_eqns[oO][vV]{c} += #1423# - <[ov]|[ov]> 1@1 t1t1_iAMe[oV][Ov]{c}");
 
180
 
 
181
  blas->append("t2_eqns[oO][vV]{c} += #3124# - ([ov]|[vo]) 1@2 t1t1_IAME[OV][OV]{c}");
 
182
 
 
183
  blas->append("t2_eqns[oO][vV]{c} += #1234#   t1[o][v]{c} 2@1 <[v]|[ovv]>");
 
184
  blas->append("t2_eqns[oO][vV]{c} += #2143#   t1[O][V]{c} 2@1 <[v]|[ovv]>");
 
185
 
 
186
  blas->append("t2_eqns[oO][vV]{c} += #3412# - t1[o][v]{c} 1@1 <[o]|[voo]>");
 
187
  blas->append("t2_eqns[oO][vV]{c} += #4321# - t1[O][V]{c} 1@1 <[o]|[voo]>");
 
188
 
 
189
 
 
190
  // Open-shell
 
191
  blas->append("t2_eqns[oO][vV]{o}  = <[oo]|[vv]>");
 
192
 
 
193
  blas->append("t2_eqns[oO][vV]{o} += #3214# t2[V][vOo]{o} 1@2 F'_AE[V][V]{o}");
 
194
  blas->append("t2_eqns[oO][vV]{o} += #4123# t2[v][VoO]{o} 1@2 F'_ae[v][v]{o}");
 
195
 
 
196
  blas->append("t2_eqns[oO][vV]{o} += #1432# - t2[O][oVv]{o} 1@1 F'_MI[O][O]{o}");
 
197
  blas->append("t2_eqns[oO][vV]{o} += #2341# - t2[o][OvV]{o} 1@1 F'_mi[o][o]{o}");
 
198
 
 
199
  blas->append("t2_eqns[oO][vV]{o} += W_mNiJ[oO][oO]{o} 1@1 tau[oO][vV]{o}");
 
200
 
 
201
  //blas->append("t2_eqns[oO][vV]{o} += tau[oO][vV]{o} 2@2 <[vv]|[vv]>");
 
202
 
 
203
 
 
204
 
 
205
  blas->append("t2_eqns[oO][vV]{o} += tau[oO][v>=V]{o} 2@2 <[vv]|[v>=v]>");
 
206
  blas->append("t2_eqns[oO][vV]{o} += #1243# tau[oO][V>=v]{o} 2@2 <[vv]|[v>=v]>");
 
207
 
 
208
 
 
209
  blas->append("t2_eqns[oO][vV]{o} += #1234#  - Z_iJaM[oOv][O]{o} 2@1 t1[O][V]{o}");
 
210
  blas->append("t2_eqns[oO][vV]{o} += #1243#    Z_iJAm[oOV][o]{o} 2@1 t1[o][v]{o}");
 
211
 
 
212
  blas->append("t2_eqns[oO][vV]{o} += #2413#   W_JBme[OV][ov]{o} 2@2 t2[ov][ov]{o}");
 
213
  blas->append("t2_eqns[oO][vV]{o} += #2413#   W_JBME[OV][OV]{o} 2@2 t2[ov][OV]{o}");
 
214
 
 
215
  blas->append("t2_eqns[oO][vV]{o} += #2314#   W_JbMe[Ov][Ov]{o} 2@2 t2[oV][Ov]{o}");
 
216
  blas->append("t2_eqns[oO][vV]{o} += #1423#   W_jBmE[oV][oV]{o} 2@1 t2[oV][Ov]{o}");
 
217
 
 
218
  blas->append("t2_eqns[oO][vV]{o} += #1324#   W_jbME[ov][OV]{o} 2@2 t2[OV][OV]{o}");
 
219
 
 
220
  blas->append("t2_eqns[oO][vV]{o} += #1324#   W_jbme[ov][ov]{o} 2@1 t2[ov][OV]{o}");
 
221
 
 
222
 
 
223
  blas->append("t2_eqns[oO][vV]{o} += #4213# - ([ov]|[vo]) 1@2 t1t1_iame[ov][ov]{o}");
 
224
 
 
225
  blas->append("t2_eqns[oO][vV]{o} += #2314# - <[ov]|[ov]> 1@2 t1t1_iAMe[oV][Ov]{o}");
 
226
 
 
227
  blas->append("t2_eqns[oO][vV]{o} += #1423# - <[ov]|[ov]> 1@1 t1t1_iAMe[oV][Ov]{o}");
 
228
 
 
229
  blas->append("t2_eqns[oO][vV]{o} += #3124# - ([ov]|[vo]) 1@2 t1t1_IAME[OV][OV]{o}");
 
230
 
 
231
  blas->append("t2_eqns[oO][vV]{o} += #1234#   t1[o][v]{o} 2@1 <[v]|[ovv]>");
 
232
  blas->append("t2_eqns[oO][vV]{o} += #2143#   t1[O][V]{o} 2@1 <[v]|[ovv]>");
 
233
 
 
234
  blas->append("t2_eqns[oO][vV]{o} += #3412# - t1[o][v]{o} 1@1 <[o]|[voo]>");
 
235
  blas->append("t2_eqns[oO][vV]{o} += #4321# - t1[O][V]{o} 1@1 <[o]|[voo]>");
 
236
 
 
237
  DEBUGGING(3,
 
238
    blas->print("t2_eqns[oO][vV]{c}");
 
239
    blas->print("t2_eqns[oO][vV]{o}");
 
240
  )
 
241
 
 
242
  DEBUGGING(1,
 
243
    fprintf(outfile," done. Timing %20.6f s",timer.get());
 
244
    fflush(outfile);
 
245
  )
 
246
}
 
247
 
 
248
void CCMRCC::build_t2_IJAB_amplitudes()
 
249
{
 
250
  Timer timer;
 
251
  DEBUGGING(1,
 
252
    fprintf(outfile,"\n\tBuilding the t2_IJAB Amplitudes   ...");
 
253
    fflush(outfile);
 
254
  )
 
255
  // Closed-shell
 
256
  blas->append("t2_eqns[OO][VV]{c}  = t2_eqns[oo][vv]{c}");
 
257
 
 
258
  // Open-shell
 
259
  blas->append("t2_eqns[OO][VV]{o}  = <[oo]:[vv]>");
 
260
 
 
261
  blas->append("t2_eqns[OO][VV]{o} += #3124# - t2[V][VOO]{o} 1@2 F'_AE[V][V]{o}");
 
262
  blas->append("t2_eqns[OO][VV]{o} += #4123#   t2[V][VOO]{o} 1@2 F'_AE[V][V]{o}");
 
263
 
 
264
  blas->append("t2_eqns[OO][VV]{o} += #1342#   t2[O][OVV]{o} 1@1 F'_MI[O][O]{o}");
 
265
  blas->append("t2_eqns[OO][VV]{o} += #2341# - t2[O][OVV]{o} 1@1 F'_MI[O][O]{o}");
 
266
 
 
267
  blas->append("t2_eqns[OO][VV]{o} += 1/2  W_MNIJ[OO][OO]{o} 1@1 tau[OO][VV]{o}");
 
268
 
 
269
  blas->append("t2_eqns[OO][V>V]{o} = tau[OO][V>V]{o} 2@2 <[v>v]:[v>v]>");
 
270
 
 
271
  blas->append("t2_eqns[OO][VV]{o} +>= #1234#  t2_eqns[OO][V>V]{o}");
 
272
 
 
273
  blas->append("t2_eqns[OO][VV]{o} += #1234# - Z_IJAM[OOV][O]{o} 2@1 t1[O][V]{o}");
 
274
  blas->append("t2_eqns[OO][VV]{o} += #1243#   Z_IJAM[OOV][O]{o} 2@1 t1[O][V]{o}");
 
275
 
 
276
  blas->append("t2_eqns[OO][VV]{o} += #2413#   W_JBME[OV][OV]{o} 2@2 t2[OV][OV]{o}");
 
277
  blas->append("t2_eqns[OO][VV]{o} += #2314# - W_JBME[OV][OV]{o} 2@2 t2[OV][OV]{o}");
 
278
 
 
279
  blas->append("t2_eqns[OO][VV]{o} += #1423# - W_JBME[OV][OV]{o} 2@2 t2[OV][OV]{o}");
 
280
  blas->append("t2_eqns[OO][VV]{o} += #1324#   W_JBME[OV][OV]{o} 2@2 t2[OV][OV]{o}");
 
281
 
 
282
  blas->append("t2_eqns[OO][VV]{o} += #4213# - ([ov]:[vo]) 1@2 t1t1_IAME[OV][OV]{o}");
 
283
  blas->append("t2_eqns[OO][VV]{o} += #3214# + ([ov]:[vo]) 1@2 t1t1_IAME[OV][OV]{o}");
 
284
  blas->append("t2_eqns[OO][VV]{o} += #4123# + ([ov]:[vo]) 1@2 t1t1_IAME[OV][OV]{o}");
 
285
  blas->append("t2_eqns[OO][VV]{o} += #3124# - ([ov]:[vo]) 1@2 t1t1_IAME[OV][OV]{o}");
 
286
 
 
287
  blas->append("t2_eqns[OO][VV]{o} += #2413#   W_JBme[OV][ov]{o} 2@1 t2[ov][OV]{o}");
 
288
  blas->append("t2_eqns[OO][VV]{o} += #2314# - W_JBme[OV][ov]{o} 2@1 t2[ov][OV]{o}");
 
289
 
 
290
  blas->append("t2_eqns[OO][VV]{o} += #1423# - W_JBme[OV][ov]{o} 2@1 t2[ov][OV]{o}");
 
291
  blas->append("t2_eqns[OO][VV]{o} += #1324#   W_JBme[OV][ov]{o} 2@1 t2[ov][OV]{o}");
 
292
 
 
293
  blas->append("t2_eqns[OO][VV]{o} += #1234#   t1[O][V]{o} 2@1 <[v]:[ovv]>");
 
294
  blas->append("t2_eqns[OO][VV]{o} += #2134# - t1[O][V]{o} 2@1 <[v]:[ovv]>");
 
295
 
 
296
  blas->append("t2_eqns[OO][VV]{o} += #3412# - t1[O][V]{o} 1@1 <[o]:[voo]>");
 
297
  blas->append("t2_eqns[OO][VV]{o} += #4312#   t1[O][V]{o} 1@1 <[o]:[voo]>");
 
298
 
 
299
  DEBUGGING(3,blas->print("t2_eqns[OO][VV]{o}"););
 
300
 
 
301
  DEBUGGING(1,
 
302
    fprintf(outfile," done. Timing %20.6f s",timer.get());
 
303
    fflush(outfile);
 
304
  )
 
305
}
 
306
 
 
307
void CCMRCC::build_t2_amplitudes_triples()
 
308
{
 
309
  Timer timer;
 
310
  DEBUGGING(1,
 
311
    fprintf(outfile,"\n\tBuilding the T3->T2 Amplitudes   ...");
 
312
    fflush(outfile);
 
313
  )
 
314
  build_t2_ijab_amplitudes_triples_diagram1();
 
315
  build_t2_iJaB_amplitudes_triples_diagram1();
 
316
  build_t2_IJAB_amplitudes_triples_diagram1();
 
317
 
 
318
  build_t2_ijab_amplitudes_triples_diagram2();
 
319
  build_t2_iJaB_amplitudes_triples_diagram2();
 
320
  build_t2_IJAB_amplitudes_triples_diagram2();
 
321
 
 
322
  build_t2_ijab_amplitudes_triples_diagram3();
 
323
  build_t2_iJaB_amplitudes_triples_diagram3();
 
324
  build_t2_IJAB_amplitudes_triples_diagram3();
 
325
  DEBUGGING(1,
 
326
    fprintf(outfile," done. Timing %20.6f s",timer.get());
 
327
    fflush(outfile);
 
328
  )
 
329
}
 
330
 
 
331
/**
 
332
 * @brief Computes the contraction
 
333
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{abe} W_{mnje} \f]
 
334
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
335
 * This algorithm hence performs
 
336
 * \f[ \frac{1}{2}\sum_{mne} t_{imn}^{abe} W_{mnje} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
337
 * \f[ \sum_{mNE} t_{imN}^{abE} W_{mNjE} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
338
 */
 
339
void CCMRCC::build_t2_ijab_amplitudes_triples_diagram1()
 
340
{
 
341
  // Loop over references
 
342
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
343
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
344
 
 
345
    // Grab the temporary matrices
 
346
    CCMatTmp  TijkabcMatTmp = blas->get_MatTmp("t3[ooo][vvv]",unique_ref,none);
 
347
    double*** Tijkabc_matrix = TijkabcMatTmp->get_matrix();
 
348
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
349
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
350
 
 
351
    CCMatTmp  WkijaMatTmp = blas->get_MatTmp("W_kija[o][oov]",unique_ref,none);
 
352
    CCMatTmp  WkiJAMatTmp = blas->get_MatTmp("W_kiJA[o][oOV]",unique_ref,none);
 
353
    double*** Wkija_matrix = WkijaMatTmp->get_matrix();
 
354
    double*** WkiJA_matrix = WkiJAMatTmp->get_matrix();
 
355
 
 
356
    CCMatTmp  HijabMatTmp = blas->get_MatTmp("t2_eqns[oo][vv]",unique_ref,none);
 
357
 
 
358
    // Grab the indexing for t3[iab][jkc]
 
359
    CCIndex* iab_indexing = blas->get_index("[ovv]");
 
360
    CCIndex* jkc_indexing = blas->get_index("[oov]");
 
361
    CCIndex*   j_indexing = blas->get_index("[o]");
 
362
 
 
363
 
 
364
    short** iab_tuples = iab_indexing->get_tuples();
 
365
    short** jkc_tuples = jkc_indexing->get_tuples();
 
366
 
 
367
    // PART A: Sort T[ijk][abc]->T[iab][jkc]
 
368
    double ***T_iabjkc;
 
369
    double ***H_iabj;
 
370
 
 
371
    allocate1(double**,T_iabjkc,moinfo->get_nirreps());
 
372
    allocate1(double**,H_iabj,moinfo->get_nirreps());
 
373
 
 
374
    for(int h =0; h < moinfo->get_nirreps();h++){
 
375
      // Allocate a block of T_iabjkc   
 
376
      allocate2(double,T_iabjkc[h],iab_indexing->get_pairpi(h),jkc_indexing->get_pairpi(h));
 
377
      allocate2(double,H_iabj[h],iab_indexing->get_pairpi(h),j_indexing->get_pairpi(h));
 
378
 
 
379
      size_t iab_offset = iab_indexing->get_first(h);
 
380
      size_t jkc_offset = jkc_indexing->get_first(h);
 
381
 
 
382
      // AAA Contribution
 
383
 
 
384
      // Sort this block
 
385
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
386
        int i = iab_tuples[iab_offset + iab][0];
 
387
        int a = iab_tuples[iab_offset + iab][1];
 
388
        int b = iab_tuples[iab_offset + iab][2];        
 
389
        for(int jkc = 0;jkc<jkc_indexing->get_pairpi(h);jkc++){
 
390
          int j = jkc_tuples[jkc_offset + jkc][0];
 
391
          int k = jkc_tuples[jkc_offset + jkc][1];
 
392
          int c = jkc_tuples[jkc_offset + jkc][2];
 
393
          T_iabjkc[h][iab][jkc] = TijkabcMatTmp->get_six_address_element(i,j,k,a,b,c);
 
394
        }
 
395
      }
 
396
      int m = iab_indexing->get_pairpi(h);
 
397
      int n = j_indexing->get_pairpi(h);
 
398
      int k = jkc_indexing->get_pairpi(h);
 
399
      if(m*n*k){
 
400
        double alpha = 1.0;
 
401
        double beta  = 0.0;
 
402
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
403
                   &(Wkija_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
404
      }
 
405
 
 
406
      // Sort the result
 
407
      size_t j_offset = j_indexing->get_first(h);
 
408
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
409
        int i = iab_tuples[iab_offset + iab][0];
 
410
        int a = iab_tuples[iab_offset + iab][1];
 
411
        int b = iab_tuples[iab_offset + iab][2];
 
412
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
413
          int abs_j = j + j_offset;
 
414
          HijabMatTmp->add_four_address_element(i,abs_j,a,b,-0.5*H_iabj[h][iab][j]);
 
415
          HijabMatTmp->add_four_address_element(abs_j,i,a,b,0.5*H_iabj[h][iab][j]);
 
416
        }
 
417
      }
 
418
 
 
419
      // AAB Contribution
 
420
 
 
421
      // Sort this block
 
422
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
423
        int i = iab_tuples[iab_offset + iab][0];
 
424
        int a = iab_tuples[iab_offset + iab][1];
 
425
        int b = iab_tuples[iab_offset + iab][2];        
 
426
        for(int jkc = 0;jkc<jkc_indexing->get_pairpi(h);jkc++){
 
427
          int j = jkc_tuples[jkc_offset + jkc][0];
 
428
          int k = jkc_tuples[jkc_offset + jkc][1];
 
429
          int c = jkc_tuples[jkc_offset + jkc][2];
 
430
          T_iabjkc[h][iab][jkc] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
431
        }
 
432
      }
 
433
      m = iab_indexing->get_pairpi(h);
 
434
      n = j_indexing->get_pairpi(h);
 
435
      k = jkc_indexing->get_pairpi(h);
 
436
      if(m*n*k){
 
437
        double alpha = 1.0;
 
438
        double beta  = 0.0;
 
439
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
440
                   &(WkiJA_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
441
      }
 
442
 
 
443
      // Sort the result
 
444
      j_offset = j_indexing->get_first(h);
 
445
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
446
        int i = iab_tuples[iab_offset + iab][0];
 
447
        int a = iab_tuples[iab_offset + iab][1];
 
448
        int b = iab_tuples[iab_offset + iab][2];
 
449
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
450
          int abs_j = j + j_offset;
 
451
          HijabMatTmp->add_four_address_element(i,abs_j,a,b,-H_iabj[h][iab][j]);
 
452
          HijabMatTmp->add_four_address_element(abs_j,i,a,b,H_iabj[h][iab][j]);
 
453
        }
 
454
      }
 
455
 
 
456
      // Deallocate the memory for the block
 
457
      release2(H_iabj[h]);
 
458
      release2(T_iabjkc[h]);
 
459
    }
 
460
    release1(H_iabj);
 
461
    release1(T_iabjkc);
 
462
  }
 
463
//   blas->print("t2_test[oo][vv]{u}");
 
464
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oo][vv]{u} . t2_test[oo][vv]{u}");
 
465
//   blas->print("ERROR{u}");
 
466
}
 
467
 
 
468
/**
 
469
 * @brief Computes the contraction
 
470
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{aBe} W_{mnJe} \f]
 
471
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
472
 * This algorithm hence performs
 
473
 * \f[ \sum_{Mne} t_{inM}^{aeB} W_{MnJe} \rightarrow  \{ -\bar{H}_{iJ}^{aB} \} \f]
 
474
 * \f[ \sum_{mNe} t_{mNJ}^{aBE} W_{mNiE} \rightarrow  \{ \bar{H}_{iJ}^{aB} \} \f]
 
475
 * \f[ \frac{1}{2} \sum_{mne} t_{mnJ}^{aeB} W_{mnie} \rightarrow  \{ -\bar{H}_{iJ}^{aB} \} \f]
 
476
 * \f[ \frac{1}{2} \sum_{NME} t_{iNM}^{aBE} W_{MNJE} \rightarrow  \{ \bar{H}_{iJ}^{aB} \} \f]
 
477
 */
 
478
void CCMRCC::build_t2_iJaB_amplitudes_triples_diagram1()
 
479
{
 
480
  // Loop over references
 
481
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
482
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
483
 
 
484
    // Grab the temporary matrices
 
485
    CCMatTmp  TijKabCMatTmp  = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
486
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
487
    CCMatTmp  TiJKaBCMatTmp  = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
488
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
489
 
 
490
    CCMatTmp  WkijaMatTmp = blas->get_MatTmp("W_kija[o][oov]",unique_ref,none);
 
491
    double*** Wkija_matrix = WkijaMatTmp->get_matrix();
 
492
    CCMatTmp  WKIjaMatTmp = blas->get_MatTmp("W_KIja[O][Oov]",unique_ref,none);
 
493
    double*** WKIja_matrix = WKIjaMatTmp->get_matrix();
 
494
    CCMatTmp  WkiJAMatTmp = blas->get_MatTmp("W_kiJA[o][oOV]",unique_ref,none);
 
495
    double*** WkiJA_matrix = WkiJAMatTmp->get_matrix();
 
496
    CCMatTmp  WKIJAMatTmp = blas->get_MatTmp("W_KIJA[O][OOV]",unique_ref,none);
 
497
    double*** WKIJA_matrix = WKIJAMatTmp->get_matrix();
 
498
 
 
499
 
 
500
    CCMatTmp  HiJaBMatTmp = blas->get_MatTmp("t2_eqns[oO][vV]",unique_ref,none);
 
501
 
 
502
    // Grab the indexing for t3[iab][jkc]
 
503
    CCIndex* jab_indexing = blas->get_index("[ovv]");
 
504
    CCIndex* kac_indexing = blas->get_index("[ovv]");
 
505
    CCIndex* iab_indexing = blas->get_index("[ovv]");
 
506
    CCIndex* iac_indexing = blas->get_index("[ovv]");
 
507
    CCIndex* kab_indexing = blas->get_index("[ovv]");
 
508
    CCIndex* kjb_indexing = blas->get_index("[oov]");
 
509
    CCIndex* kjc_indexing = blas->get_index("[oov]");
 
510
    CCIndex* ijc_indexing = blas->get_index("[oov]");
 
511
    CCIndex* ijb_indexing = blas->get_index("[oov]");
 
512
    CCIndex*   j_indexing = blas->get_index("[o]");
 
513
    CCIndex*   i_indexing = blas->get_index("[o]");
 
514
 
 
515
 
 
516
    short** iab_tuples = iab_indexing->get_tuples();
 
517
    short** jab_tuples = jab_indexing->get_tuples();
 
518
    short** iac_tuples = iac_indexing->get_tuples();
 
519
    short** kab_tuples = kab_indexing->get_tuples();
 
520
    short** kac_tuples = kac_indexing->get_tuples();
 
521
 
 
522
    short** kjb_tuples = kjb_indexing->get_tuples();
 
523
    short** kjc_tuples = kjc_indexing->get_tuples();
 
524
    short** ijc_tuples = ijc_indexing->get_tuples();
 
525
    short** ijb_tuples = ijb_indexing->get_tuples();
 
526
 
 
527
 
 
528
    double ***T_iackjb;
 
529
    double ***H_iabj;
 
530
 
 
531
    allocate1(double**,T_iackjb,moinfo->get_nirreps());
 
532
    allocate1(double**,H_iabj,moinfo->get_nirreps());
 
533
 
 
534
    for(int h =0; h < moinfo->get_nirreps();h++){
 
535
      // Allocate a block of T_iabjkc   
 
536
      allocate2(double,T_iackjb[h],iac_indexing->get_pairpi(h),kjb_indexing->get_pairpi(h));
 
537
      allocate2(double,H_iabj[h],iab_indexing->get_pairpi(h),j_indexing->get_pairpi(h));
 
538
 
 
539
      size_t iab_offset = iab_indexing->get_first(h);
 
540
      size_t iac_offset = iac_indexing->get_first(h);
 
541
      size_t kjb_offset = kjb_indexing->get_first(h);
 
542
 
 
543
      // AAA Contribution
 
544
 
 
545
      // PART A: Sort T[ijk][abc]->T[iac][kjb]
 
546
      for(int iac = 0;iac<iac_indexing->get_pairpi(h);iac++){
 
547
        int i = iac_tuples[iac_offset + iac][0];
 
548
        int a = iac_tuples[iac_offset + iac][1];
 
549
        int c = iac_tuples[iac_offset + iac][2];        
 
550
        for(int kjb = 0;kjb<kjb_indexing->get_pairpi(h);kjb++){
 
551
          int k = kjb_tuples[kjb_offset + kjb][0];
 
552
          int j = kjb_tuples[kjb_offset + kjb][1];
 
553
          int b = kjb_tuples[kjb_offset + kjb][2];
 
554
          T_iackjb[h][iac][kjb] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
555
        }
 
556
      }
 
557
      int m = iac_indexing->get_pairpi(h);
 
558
      int n = j_indexing->get_pairpi(h);
 
559
      int k = kjb_indexing->get_pairpi(h);
 
560
      if(m*n*k){
 
561
        double alpha = 1.0;
 
562
        double beta  = 0.0;
 
563
        C_DGEMM_22(m, n, k, alpha,&(T_iackjb[h][0][0]), k, 
 
564
                   &(WKIja_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
565
      }
 
566
 
 
567
      // Sort the result
 
568
      size_t j_offset = j_indexing->get_first(h);
 
569
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
570
        int i = iab_tuples[iab_offset + iab][0];
 
571
        int a = iab_tuples[iab_offset + iab][1];
 
572
        int b = iab_tuples[iab_offset + iab][2];
 
573
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
574
          int abs_j = j + j_offset;
 
575
          HiJaBMatTmp->add_four_address_element(i,abs_j,a,b,-H_iabj[h][iab][j]);
 
576
        }
 
577
      }
 
578
 
 
579
 
 
580
      size_t kab_offset = kab_indexing->get_first(h);
 
581
      size_t ijc_offset = ijc_indexing->get_first(h);
 
582
 
 
583
      // ABB Contribution
 
584
 
 
585
      // Sort this block
 
586
      for(int kab = 0;kab<kab_indexing->get_pairpi(h);kab++){
 
587
        int k = kab_tuples[kab_offset + kab][0];
 
588
        int a = kab_tuples[kab_offset + kab][1];
 
589
        int b = kab_tuples[kab_offset + kab][2];        
 
590
        for(int ijc = 0;ijc<ijc_indexing->get_pairpi(h);ijc++){
 
591
          int i = ijc_tuples[ijc_offset + ijc][0];
 
592
          int j = ijc_tuples[ijc_offset + ijc][1];
 
593
          int c = ijc_tuples[ijc_offset + ijc][2];
 
594
          T_iackjb[h][kab][ijc] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
595
        }
 
596
      }
 
597
 
 
598
      m = kab_indexing->get_pairpi(h);
 
599
      n = i_indexing->get_pairpi(h);
 
600
      k = ijc_indexing->get_pairpi(h);
 
601
      if(m*n*k){
 
602
        double alpha = 1.0;
 
603
        double beta  = 0.0;
 
604
        C_DGEMM_22(m, n, k, alpha,&(T_iackjb[h][0][0]), k, 
 
605
                   &(WkiJA_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
606
      }
 
607
 
 
608
      // Sort the result
 
609
      size_t jab_offset = jab_indexing->get_first(h);
 
610
      size_t i_offset = i_indexing->get_first(h);
 
611
      for(int jab = 0;jab<jab_indexing->get_pairpi(h);jab++){
 
612
        int j = jab_tuples[jab_offset + jab][0];
 
613
        int a = jab_tuples[jab_offset + jab][1];
 
614
        int b = jab_tuples[jab_offset + jab][2];
 
615
        for(int i = 0;i<i_indexing->get_pairpi(h);i++){
 
616
          int abs_i = i + i_offset;
 
617
          HiJaBMatTmp->add_four_address_element(abs_i,j,a,b,H_iabj[h][jab][i]);
 
618
        }
 
619
      }
 
620
 
 
621
      // t[mnJ][aeB] W[i][mne] Contribution
 
622
 
 
623
      size_t kac_offset = kac_indexing->get_first(h);
 
624
      size_t ijb_offset = ijb_indexing->get_first(h);
 
625
 
 
626
      // Sort this block
 
627
      for(int kac = 0;kac<kac_indexing->get_pairpi(h);kac++){
 
628
        int k = kac_tuples[kac_offset + kac][0];
 
629
        int a = kac_tuples[kac_offset + kac][1];
 
630
        int c = kac_tuples[kac_offset + kac][2];        
 
631
        for(int ijb = 0;ijb<ijb_indexing->get_pairpi(h);ijb++){
 
632
          int i = ijb_tuples[ijb_offset + ijb][0];
 
633
          int j = ijb_tuples[ijb_offset + ijb][1];
 
634
          int b = ijb_tuples[ijb_offset + ijb][2];
 
635
          T_iackjb[h][kac][ijb] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
636
        }
 
637
      }
 
638
 
 
639
      m = kac_indexing->get_pairpi(h);
 
640
      n = i_indexing->get_pairpi(h);
 
641
      k = ijb_indexing->get_pairpi(h);
 
642
      if(m*n*k){
 
643
        double alpha = 1.0;
 
644
        double beta  = 0.0;
 
645
        C_DGEMM_22(m, n, k, alpha,&(T_iackjb[h][0][0]), k, 
 
646
                   &(Wkija_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
647
      }
 
648
 
 
649
      // Sort the result
 
650
      for(int jab = 0;jab<jab_indexing->get_pairpi(h);jab++){
 
651
        int j = jab_tuples[jab_offset + jab][0];
 
652
        int a = jab_tuples[jab_offset + jab][1];
 
653
        int b = jab_tuples[jab_offset + jab][2];
 
654
        for(int i = 0;i<i_indexing->get_pairpi(h);i++){
 
655
          int abs_i = i + i_offset;
 
656
          HiJaBMatTmp->add_four_address_element(abs_i,j,a,b,-0.5*H_iabj[h][jab][i]);
 
657
        }
 
658
      }
 
659
 
 
660
      // t[iaB][MNE] W[J][MNE] Contribution
 
661
      size_t kjc_offset = kjc_indexing->get_first(h);
 
662
 
 
663
      // Sort this block
 
664
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
665
        int i = iab_tuples[iab_offset + iab][0];
 
666
        int a = iab_tuples[iab_offset + iab][1];
 
667
        int b = iab_tuples[iab_offset + iab][2];        
 
668
        for(int kjc = 0;kjc<kjc_indexing->get_pairpi(h);kjc++){
 
669
          int k = kjc_tuples[kjc_offset + kjc][0];
 
670
          int j = kjc_tuples[kjc_offset + kjc][1];
 
671
          int c = kjc_tuples[kjc_offset + kjc][2];
 
672
          T_iackjb[h][iab][kjc] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
673
        }
 
674
      }
 
675
 
 
676
      m = iab_indexing->get_pairpi(h);
 
677
      n = j_indexing->get_pairpi(h);
 
678
      k = kjc_indexing->get_pairpi(h);
 
679
      if(m*n*k){
 
680
        double alpha = 1.0;
 
681
        double beta  = 0.0;
 
682
        C_DGEMM_22(m, n, k, alpha,&(T_iackjb[h][0][0]), k, 
 
683
                   &(WKIJA_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
684
      }
 
685
 
 
686
      // Sort the result
 
687
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
688
        int i = iab_tuples[iab_offset + iab][0];
 
689
        int a = iab_tuples[iab_offset + iab][1];
 
690
        int b = iab_tuples[iab_offset + iab][2];
 
691
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
692
          int abs_j = j + j_offset;
 
693
          HiJaBMatTmp->add_four_address_element(i,abs_j,a,b,0.5*H_iabj[h][iab][j]);
 
694
        }
 
695
      }
 
696
      // Deallocate the memory for the block
 
697
      release2(T_iackjb[h]);
 
698
      release2(H_iabj[h]);
 
699
    }
 
700
    release1(H_iabj);
 
701
    release1(T_iackjb);
 
702
  }
 
703
//   blas->print("t2_test[oO][vV]{u}");
 
704
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oO][vV]{u} . t2_test[oO][vV]{u}");
 
705
//   blas->print("ERROR{u}");
 
706
}
 
707
 
 
708
 
 
709
/**
 
710
 * @brief Computes the contraction
 
711
 * \f[ -\frac{1}{2} P(IJ)\sum_{mne} t_{Imn}^{ABe} W_{mnJe} \f]
 
712
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
713
 * This algorithm hence performs
 
714
 * \f[ \sum_{Mne} t_{nMI}^{eAB} W_{MnJe} \rightarrow  \{\bar{H}_{IJ}^{AB},-\bar{H}_{JI}^{AB} \} \f]
 
715
 * \f[ \frac{1}{2}\sum_{MNE} t_{IMN}^{ABE} W_{MNJE} \rightarrow  \{ -\bar{H}_{IJ}^{AB},\bar{H}_{JI}^{AB} \} \f]
 
716
 */
 
717
void CCMRCC::build_t2_IJAB_amplitudes_triples_diagram1()
 
718
{
 
719
  // Loop over references
 
720
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
721
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
722
 
 
723
    // Grab the temporary matrices
 
724
    CCMatTmp  TiJKaBCMatTmp = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
725
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
726
    CCMatTmp  TIJKABCMatTmp = blas->get_MatTmp("t3[OOO][VVV]",unique_ref,none);
 
727
    double*** TIJKABC_matrix = TIJKABCMatTmp->get_matrix();
 
728
 
 
729
    CCMatTmp  WKIjaMatTmp = blas->get_MatTmp("W_KIja[O][Oov]",unique_ref,none);
 
730
    CCMatTmp  WKIJAMatTmp = blas->get_MatTmp("W_KIJA[O][OOV]",unique_ref,none);
 
731
    double*** WKIja_matrix = WKIjaMatTmp->get_matrix();
 
732
    double*** WKIJA_matrix = WKIJAMatTmp->get_matrix();
 
733
 
 
734
    CCMatTmp  HIJABMatTmp = blas->get_MatTmp("t2_eqns[OO][VV]",unique_ref,none);
 
735
 
 
736
    // Grab the indexing for t3[iab][jkc]
 
737
    CCIndex* iab_indexing = blas->get_index("[ovv]");
 
738
    CCIndex* kbc_indexing = blas->get_index("[ovv]");
 
739
    CCIndex* jia_indexing = blas->get_index("[oov]");
 
740
    CCIndex* jkc_indexing = blas->get_index("[oov]");
 
741
    CCIndex*   j_indexing = blas->get_index("[o]");
 
742
 
 
743
    short** iab_tuples = iab_indexing->get_tuples();
 
744
    short** kbc_tuples = kbc_indexing->get_tuples();
 
745
    short** jia_tuples = jia_indexing->get_tuples();
 
746
    short** jkc_tuples = jkc_indexing->get_tuples();
 
747
 
 
748
    double ***T_iabjkc;
 
749
    double ***H_iabj;
 
750
 
 
751
    allocate1(double**,T_iabjkc,moinfo->get_nirreps());
 
752
    allocate1(double**,H_iabj,moinfo->get_nirreps());
 
753
 
 
754
    for(int h =0; h < moinfo->get_nirreps();h++){
 
755
      // Allocate a block of T_iabjkc   
 
756
      allocate2(double,T_iabjkc[h],kbc_indexing->get_pairpi(h),jia_indexing->get_pairpi(h));
 
757
      allocate2(double,H_iabj[h],kbc_indexing->get_pairpi(h),j_indexing->get_pairpi(h));
 
758
 
 
759
      size_t iab_offset = iab_indexing->get_first(h);
 
760
      size_t kbc_offset = kbc_indexing->get_first(h);
 
761
      size_t jia_offset = jia_indexing->get_first(h);
 
762
      size_t jkc_offset = jkc_indexing->get_first(h);
 
763
 
 
764
      // AAA Contribution
 
765
 
 
766
      // PART A: Sort T[ijk][abc]->T[kbc][jia]
 
767
      // Sort this block
 
768
      for(int kbc = 0;kbc<kbc_indexing->get_pairpi(h);kbc++){
 
769
        int k = kbc_tuples[kbc_offset + kbc][0];
 
770
        int b = kbc_tuples[kbc_offset + kbc][1];
 
771
        int c = kbc_tuples[kbc_offset + kbc][2];        
 
772
        for(int jia = 0;jia<jia_indexing->get_pairpi(h);jia++){
 
773
          int j = jia_tuples[jia_offset + jia][0];
 
774
          int i = jia_tuples[jia_offset + jia][1];
 
775
          int a = jia_tuples[jia_offset + jia][2];
 
776
          T_iabjkc[h][kbc][jia] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
777
        }
 
778
      }
 
779
      int m = kbc_indexing->get_pairpi(h);
 
780
      int n = j_indexing->get_pairpi(h);
 
781
      int k = jia_indexing->get_pairpi(h);
 
782
      if(m*n*k){
 
783
        double alpha = 1.0;
 
784
        double beta  = 0.0;
 
785
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
786
                   &(WKIja_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
787
      }
 
788
 
 
789
      // Sort the result
 
790
      size_t j_offset = j_indexing->get_first(h);
 
791
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
792
        int i = iab_tuples[iab_offset + iab][0];
 
793
        int a = iab_tuples[iab_offset + iab][1];
 
794
        int b = iab_tuples[iab_offset + iab][2];
 
795
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
796
          int abs_j = j + j_offset;
 
797
          HIJABMatTmp->add_four_address_element(i,abs_j,a,b,H_iabj[h][iab][j]);
 
798
          HIJABMatTmp->add_four_address_element(abs_j,i,a,b,-H_iabj[h][iab][j]);
 
799
        }
 
800
      }
 
801
 
 
802
      // BBB Contribution
 
803
 
 
804
      // Sort this block
 
805
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
806
        int i = iab_tuples[iab_offset + iab][0];
 
807
        int a = iab_tuples[iab_offset + iab][1];
 
808
        int b = iab_tuples[iab_offset + iab][2];        
 
809
        for(int jkc = 0;jkc<jkc_indexing->get_pairpi(h);jkc++){
 
810
          int j = jkc_tuples[jkc_offset + jkc][0];
 
811
          int k = jkc_tuples[jkc_offset + jkc][1];
 
812
          int c = jkc_tuples[jkc_offset + jkc][2];
 
813
          T_iabjkc[h][iab][jkc] = TIJKABCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
814
        }
 
815
      }
 
816
      m = iab_indexing->get_pairpi(h);
 
817
      n = j_indexing->get_pairpi(h);
 
818
      k = jkc_indexing->get_pairpi(h);
 
819
      if(m*n*k){
 
820
        double alpha = 1.0;
 
821
        double beta  = 0.0;
 
822
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
823
                   &(WKIJA_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
824
      }
 
825
 
 
826
      // Sort the result
 
827
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
828
        int i = iab_tuples[iab_offset + iab][0];
 
829
        int a = iab_tuples[iab_offset + iab][1];
 
830
        int b = iab_tuples[iab_offset + iab][2];
 
831
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
832
          int abs_j = j + j_offset;
 
833
          HIJABMatTmp->add_four_address_element(i,abs_j,a,b,-0.5*H_iabj[h][iab][j]);
 
834
          HIJABMatTmp->add_four_address_element(abs_j,i,a,b,0.5*H_iabj[h][iab][j]);
 
835
        }
 
836
      }
 
837
 
 
838
      // Deallocate the memory for the block
 
839
      release2(H_iabj[h]);
 
840
      release2(T_iabjkc[h]);
 
841
    }
 
842
    release1(H_iabj);
 
843
    release1(T_iabjkc);
 
844
  }
 
845
//   blas->print("t2_test[OO][VV]{u}");
 
846
//   blas->solve("ERROR{u} = 1000000.0 t2_test[OO][VV]{u} . t2_test[OO][VV]{u}");
 
847
//   blas->print("ERROR{u}");
 
848
}
 
849
 
 
850
 
 
851
 
 
852
/**
 
853
 * @brief Computes the contraction
 
854
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{abe} W_{mnje} \f]
 
855
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
856
 * This algorithm hence performs
 
857
 * \f[ \frac{1}{2}\sum_{mne} t_{imn}^{abe} W_{mnje} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
858
 * \f[ \sum_{mNE} t_{imN}^{abE} W_{mNjE} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
859
 */
 
860
void CCMRCC::build_t2_ijab_amplitudes_triples_diagram2()
 
861
{
 
862
  // Loop over references
 
863
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
864
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
865
 
 
866
    // Grab the temporary matrices
 
867
    CCMatTmp  TijkabcMatTmp = blas->get_MatTmp("t3[ooo][vvv]",unique_ref,none);
 
868
    double*** Tijkabc_matrix = TijkabcMatTmp->get_matrix();
 
869
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
870
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
871
 
 
872
    CCMatTmp  WaibcMatTmp  = blas->get_MatTmp("W_aibc[v][ovv]",unique_ref,none);
 
873
    double*** Waibc_matrix = WaibcMatTmp->get_matrix();
 
874
    CCMatTmp  WaIbCMatTmp  = blas->get_MatTmp("W_aIbC[v][OvV]",unique_ref,none);
 
875
    double*** WaIbC_matrix = WaIbCMatTmp->get_matrix();
 
876
 
 
877
    CCMatTmp  HijabMatTmp = blas->get_MatTmp("t2_eqns[oo][vv]",unique_ref,none);
 
878
 
 
879
    // Grab the indexing for t3[iab][jkc]
 
880
    CCIndex* ovv_indexing = blas->get_index("[ovv]");
 
881
    CCIndex* oov_indexing = blas->get_index("[oov]");
 
882
    CCIndex*   v_indexing = blas->get_index("[v]");
 
883
 
 
884
    short** ovv_tuples = ovv_indexing->get_tuples();
 
885
    short** oov_tuples = oov_indexing->get_tuples();
 
886
 
 
887
    double ***T_oovovv;
 
888
    double ***H_ijab;
 
889
 
 
890
    allocate1(double**,T_oovovv,moinfo->get_nirreps());
 
891
    allocate1(double**,H_ijab,moinfo->get_nirreps());
 
892
 
 
893
    for(int h =0; h < moinfo->get_nirreps();h++){
 
894
      // Allocate a block of T_iabjkc   
 
895
      allocate2(double,T_oovovv[h],oov_indexing->get_pairpi(h),ovv_indexing->get_pairpi(h));
 
896
      allocate2(double,H_ijab[h],oov_indexing->get_pairpi(h),v_indexing->get_pairpi(h));
 
897
 
 
898
      size_t ovv_offset = ovv_indexing->get_first(h);
 
899
      size_t oov_offset = oov_indexing->get_first(h);
 
900
 
 
901
      // AAA Contribution
 
902
 
 
903
      // Sort this block
 
904
      for(int kbc = 0;kbc<ovv_indexing->get_pairpi(h);kbc++){
 
905
        int k = ovv_tuples[ovv_offset + kbc][0];
 
906
        int b = ovv_tuples[ovv_offset + kbc][1];
 
907
        int c = ovv_tuples[ovv_offset + kbc][2];        
 
908
        for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
909
          int i = oov_tuples[oov_offset + ija][0];
 
910
          int j = oov_tuples[oov_offset + ija][1];
 
911
          int a = oov_tuples[oov_offset + ija][2];        
 
912
          T_oovovv[h][ija][kbc] = TijkabcMatTmp->get_six_address_element(i,j,k,a,b,c);
 
913
        }
 
914
      }
 
915
      int m = oov_indexing->get_pairpi(h);
 
916
      int n =   v_indexing->get_pairpi(h);
 
917
      int k = ovv_indexing->get_pairpi(h);
 
918
      if(m*n*k){
 
919
        double alpha = 1.0;
 
920
        double beta  = 0.0;
 
921
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
922
                   &(Waibc_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
923
      }
 
924
 
 
925
      // Sort the result
 
926
      size_t v_offset = v_indexing->get_first(h);
 
927
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
928
        int i = oov_tuples[oov_offset + ija][0];
 
929
        int j = oov_tuples[oov_offset + ija][1];
 
930
        int a = oov_tuples[oov_offset + ija][2];
 
931
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
932
          int abs_b = b + v_offset;
 
933
          HijabMatTmp->add_four_address_element(i,j,a,abs_b,0.5*H_ijab[h][ija][b]);
 
934
          HijabMatTmp->add_four_address_element(i,j,abs_b,a,-0.5*H_ijab[h][ija][b]);
 
935
        }
 
936
      }
 
937
 
 
938
      // Sort this block
 
939
      for(int kbc = 0;kbc<ovv_indexing->get_pairpi(h);kbc++){
 
940
        int k = ovv_tuples[ovv_offset + kbc][0];
 
941
        int b = ovv_tuples[ovv_offset + kbc][1];
 
942
        int c = ovv_tuples[ovv_offset + kbc][2];        
 
943
        for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
944
          int i = oov_tuples[oov_offset + ija][0];
 
945
          int j = oov_tuples[oov_offset + ija][1];
 
946
          int a = oov_tuples[oov_offset + ija][2];        
 
947
          T_oovovv[h][ija][kbc] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
948
        }
 
949
      }
 
950
 
 
951
      m = oov_indexing->get_pairpi(h);
 
952
      n =   v_indexing->get_pairpi(h);
 
953
      k = ovv_indexing->get_pairpi(h);
 
954
      if(m*n*k){
 
955
        double alpha = 1.0;
 
956
        double beta  = 0.0;
 
957
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
958
                   &(WaIbC_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
959
      }
 
960
 
 
961
      // Sort the result
 
962
      v_offset = v_indexing->get_first(h);
 
963
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
964
        int i = oov_tuples[oov_offset + ija][0];
 
965
        int j = oov_tuples[oov_offset + ija][1];
 
966
        int a = oov_tuples[oov_offset + ija][2];
 
967
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
968
          int abs_b = b + v_offset;
 
969
          HijabMatTmp->add_four_address_element(i,j,a,abs_b,H_ijab[h][ija][b]);
 
970
          HijabMatTmp->add_four_address_element(i,j,abs_b,a,-H_ijab[h][ija][b]);
 
971
        }
 
972
      }
 
973
      // Deallocate the memory for the block
 
974
      release2(H_ijab[h]);
 
975
      release2(T_oovovv[h]);
 
976
    }
 
977
    release1(H_ijab);
 
978
    release1(T_oovovv);
 
979
  }
 
980
//   blas->print("t2_test[oo][vv]{u}");
 
981
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oo][vv]{u} . t2_test[oo][vv]{u}");
 
982
//   blas->print("ERROR{u}");
 
983
}
 
984
 
 
985
 
 
986
 
 
987
/**
 
988
 * @brief Computes the contraction
 
989
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{abe} W_{mnje} \f]
 
990
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
991
 * This algorithm hence performs
 
992
 * \f[ \frac{1}{2}\sum_{mne} t_{imn}^{abe} W_{mnje} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
993
 * \f[ \sum_{mNE} t_{imN}^{abE} W_{mNjE} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
994
 */
 
995
void CCMRCC::build_t2_iJaB_amplitudes_triples_diagram2()
 
996
{
 
997
  // Loop over references
 
998
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
999
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1000
 
 
1001
    // Grab the temporary matrices
 
1002
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
1003
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
1004
    CCMatTmp  TiJKaBCMatTmp = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
1005
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
1006
 
 
1007
    CCMatTmp  WaibcMatTmp  = blas->get_MatTmp("W_aibc[v][ovv]",unique_ref,none);
 
1008
    double*** Waibc_matrix = WaibcMatTmp->get_matrix();
 
1009
    CCMatTmp  WaIbCMatTmp  = blas->get_MatTmp("W_aIbC[v][OvV]",unique_ref,none);
 
1010
    double*** WaIbC_matrix = WaIbCMatTmp->get_matrix();
 
1011
    CCMatTmp  WAiBcMatTmp  = blas->get_MatTmp("W_AiBc[V][oVv]",unique_ref,none);
 
1012
    double*** WAiBc_matrix = WAiBcMatTmp->get_matrix();
 
1013
    CCMatTmp  WAIBCMatTmp  = blas->get_MatTmp("W_AIBC[V][OVV]",unique_ref,none);
 
1014
    double*** WAIBC_matrix = WAIBCMatTmp->get_matrix();
 
1015
 
 
1016
    CCMatTmp  HiJaBMatTmp = blas->get_MatTmp("t2_eqns[oO][vV]",unique_ref,none);
 
1017
 
 
1018
    // Grab the indexing for t3[iab][jkc]
 
1019
    CCIndex* ovv_indexing = blas->get_index("[ovv]");
 
1020
    CCIndex* oov_indexing = blas->get_index("[oov]");
 
1021
    CCIndex*   v_indexing = blas->get_index("[v]");
 
1022
 
 
1023
    short** ovv_tuples = ovv_indexing->get_tuples();
 
1024
    short** oov_tuples = oov_indexing->get_tuples();
 
1025
 
 
1026
    double ***T_oovovv;
 
1027
    double ***H_ijab;
 
1028
 
 
1029
    allocate1(double**,T_oovovv,moinfo->get_nirreps());
 
1030
    allocate1(double**,H_ijab,moinfo->get_nirreps());
 
1031
 
 
1032
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1033
      // Allocate a block of T_iabjkc   
 
1034
      allocate2(double,T_oovovv[h],oov_indexing->get_pairpi(h),ovv_indexing->get_pairpi(h));
 
1035
      allocate2(double,H_ijab[h],oov_indexing->get_pairpi(h),v_indexing->get_pairpi(h));
 
1036
 
 
1037
      size_t ovv_offset = ovv_indexing->get_first(h);
 
1038
      size_t oov_offset = oov_indexing->get_first(h);
 
1039
 
 
1040
      // AAA Contribution
 
1041
 
 
1042
      // Sort this block
 
1043
      for(int jab = 0;jab<ovv_indexing->get_pairpi(h);jab++){
 
1044
        int j = ovv_tuples[ovv_offset + jab][0];
 
1045
        int a = ovv_tuples[ovv_offset + jab][1];
 
1046
        int b = ovv_tuples[ovv_offset + jab][2];        
 
1047
        for(int ikc = 0;ikc<oov_indexing->get_pairpi(h);ikc++){
 
1048
          int i = oov_tuples[oov_offset + ikc][0];
 
1049
          int k = oov_tuples[oov_offset + ikc][1];
 
1050
          int c = oov_tuples[oov_offset + ikc][2];        
 
1051
          T_oovovv[h][ikc][jab] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1052
        }
 
1053
      }
 
1054
      int m = oov_indexing->get_pairpi(h);
 
1055
      int n =   v_indexing->get_pairpi(h);
 
1056
      int k = ovv_indexing->get_pairpi(h);
 
1057
      if(m*n*k){
 
1058
        double alpha = 1.0;
 
1059
        double beta  = 0.0;
 
1060
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1061
                   &(Waibc_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1062
      }
 
1063
 
 
1064
      // Sort the result
 
1065
      size_t v_offset = v_indexing->get_first(h);
 
1066
      for(int ijb = 0;ijb<oov_indexing->get_pairpi(h);ijb++){
 
1067
        int i = oov_tuples[oov_offset + ijb][0];
 
1068
        int j = oov_tuples[oov_offset + ijb][1];
 
1069
        int b = oov_tuples[oov_offset + ijb][2];
 
1070
        for(int a = 0;a<v_indexing->get_pairpi(h);a++){
 
1071
          int abs_a = a + v_offset;
 
1072
          HiJaBMatTmp->add_four_address_element(i,j,abs_a,b,0.5*H_ijab[h][ijb][a]);
 
1073
        }
 
1074
      }
 
1075
 
 
1076
      // Sort this block
 
1077
      for(int kac = 0;kac<ovv_indexing->get_pairpi(h);kac++){
 
1078
        int k = ovv_tuples[ovv_offset + kac][0];
 
1079
        int a = ovv_tuples[ovv_offset + kac][1];
 
1080
        int c = ovv_tuples[ovv_offset + kac][2];        
 
1081
        for(int ijb = 0;ijb<oov_indexing->get_pairpi(h);ijb++){
 
1082
          int i = oov_tuples[oov_offset + ijb][0];
 
1083
          int j = oov_tuples[oov_offset + ijb][1];
 
1084
          int b = oov_tuples[oov_offset + ijb][2];        
 
1085
          T_oovovv[h][ijb][kac] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1086
        }
 
1087
      }
 
1088
      m = oov_indexing->get_pairpi(h);
 
1089
      n =   v_indexing->get_pairpi(h);
 
1090
      k = ovv_indexing->get_pairpi(h);
 
1091
      if(m*n*k){
 
1092
        double alpha = 1.0;
 
1093
        double beta  = 0.0;
 
1094
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1095
                   &(WaIbC_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1096
      }
 
1097
 
 
1098
      // Sort the result
 
1099
      v_offset = v_indexing->get_first(h);
 
1100
      for(int ijb = 0;ijb<oov_indexing->get_pairpi(h);ijb++){
 
1101
        int i = oov_tuples[oov_offset + ijb][0];
 
1102
        int j = oov_tuples[oov_offset + ijb][1];
 
1103
        int b = oov_tuples[oov_offset + ijb][2];
 
1104
        for(int a = 0;a<v_indexing->get_pairpi(h);a++){
 
1105
          int abs_a = a + v_offset;
 
1106
          HiJaBMatTmp->add_four_address_element(i,j,abs_a,b,H_ijab[h][ijb][a]);
 
1107
        }
 
1108
      }
 
1109
 
 
1110
      // Sort this block
 
1111
      for(int jcb = 0;jcb<ovv_indexing->get_pairpi(h);jcb++){
 
1112
        int j = ovv_tuples[ovv_offset + jcb][0];
 
1113
        int c = ovv_tuples[ovv_offset + jcb][1];
 
1114
        int b = ovv_tuples[ovv_offset + jcb][2];        
 
1115
        for(int ika = 0;ika<oov_indexing->get_pairpi(h);ika++){
 
1116
          int i = oov_tuples[oov_offset + ika][0];
 
1117
          int k = oov_tuples[oov_offset + ika][1];
 
1118
          int a = oov_tuples[oov_offset + ika][2];        
 
1119
          T_oovovv[h][ika][jcb] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1120
        }
 
1121
      }
 
1122
      m = oov_indexing->get_pairpi(h);
 
1123
      n =   v_indexing->get_pairpi(h);
 
1124
      k = ovv_indexing->get_pairpi(h);
 
1125
      if(m*n*k){
 
1126
        double alpha = 1.0;
 
1127
        double beta  = 0.0;
 
1128
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1129
                   &(WAiBc_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1130
      }
 
1131
 
 
1132
      // Sort the result
 
1133
      v_offset = v_indexing->get_first(h);
 
1134
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1135
        int i = oov_tuples[oov_offset + ija][0];
 
1136
        int j = oov_tuples[oov_offset + ija][1];
 
1137
        int a = oov_tuples[oov_offset + ija][2];
 
1138
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
1139
          int abs_b = b + v_offset;
 
1140
          HiJaBMatTmp->add_four_address_element(i,j,a,abs_b,H_ijab[h][ija][b]);
 
1141
        }
 
1142
      }
 
1143
 
 
1144
      // Sort this block
 
1145
      for(int kbc = 0;kbc<ovv_indexing->get_pairpi(h);kbc++){
 
1146
        int k = ovv_tuples[ovv_offset + kbc][0];
 
1147
        int b = ovv_tuples[ovv_offset + kbc][1];
 
1148
        int c = ovv_tuples[ovv_offset + kbc][2];        
 
1149
        for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1150
          int i = oov_tuples[oov_offset + ija][0];
 
1151
          int j = oov_tuples[oov_offset + ija][1];
 
1152
          int a = oov_tuples[oov_offset + ija][2];        
 
1153
          T_oovovv[h][ija][kbc] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1154
        }
 
1155
      }
 
1156
      m = oov_indexing->get_pairpi(h);
 
1157
      n =   v_indexing->get_pairpi(h);
 
1158
      k = ovv_indexing->get_pairpi(h);
 
1159
      if(m*n*k){
 
1160
        double alpha = 1.0;
 
1161
        double beta  = 0.0;
 
1162
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1163
                   &(WAIBC_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1164
      }
 
1165
 
 
1166
      // Sort the result
 
1167
      v_offset = v_indexing->get_first(h);
 
1168
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1169
        int i = oov_tuples[oov_offset + ija][0];
 
1170
        int j = oov_tuples[oov_offset + ija][1];
 
1171
        int a = oov_tuples[oov_offset + ija][2];
 
1172
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
1173
          int abs_b = b + v_offset;
 
1174
          HiJaBMatTmp->add_four_address_element(i,j,a,abs_b,0.5*H_ijab[h][ija][b]);
 
1175
        }
 
1176
      }
 
1177
 
 
1178
      // Deallocate the memory for the block
 
1179
      release2(H_ijab[h]);
 
1180
      release2(T_oovovv[h]);
 
1181
    }
 
1182
    release1(H_ijab);
 
1183
    release1(T_oovovv);
 
1184
  }
 
1185
//   blas->print("t2_test[oO][vV]{u}");
 
1186
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oO][vV]{u} . t2_test[oO][vV]{u}");
 
1187
//   blas->print("ERROR{u}");
 
1188
}
 
1189
 
 
1190
/**
 
1191
 * @brief Computes the contraction
 
1192
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{abe} W_{mnje} \f]
 
1193
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
1194
 * This algorithm hence performs
 
1195
 * \f[ \frac{1}{2}\sum_{mne} t_{imn}^{abe} W_{mnje} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
1196
 * \f[ \sum_{mNE} t_{imN}^{abE} W_{mNjE} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
1197
 */
 
1198
void CCMRCC::build_t2_IJAB_amplitudes_triples_diagram2()
 
1199
{
 
1200
  // Loop over references
 
1201
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
1202
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1203
 
 
1204
    // Grab the temporary matrices
 
1205
    CCMatTmp  TiJKaBCMatTmp = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
1206
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
1207
    CCMatTmp  TIJKABCMatTmp = blas->get_MatTmp("t3[OOO][VVV]",unique_ref,none);
 
1208
    double*** TIJKABC_matrix = TIJKABCMatTmp->get_matrix();
 
1209
 
 
1210
    CCMatTmp  WAiBcMatTmp  = blas->get_MatTmp("W_AiBc[V][oVv]",unique_ref,none);
 
1211
    double*** WAiBc_matrix = WAiBcMatTmp->get_matrix();
 
1212
    CCMatTmp  WAIBCMatTmp  = blas->get_MatTmp("W_AIBC[V][OVV]",unique_ref,none);
 
1213
    double*** WAIBC_matrix = WAIBCMatTmp->get_matrix();
 
1214
 
 
1215
    CCMatTmp  HIJABMatTmp = blas->get_MatTmp("t2_eqns[OO][VV]",unique_ref,none);
 
1216
 
 
1217
    // Grab the indexing for t3[iab][jkc]
 
1218
    CCIndex* ovv_indexing = blas->get_index("[ovv]");
 
1219
    CCIndex* oov_indexing = blas->get_index("[oov]");
 
1220
    CCIndex*   v_indexing = blas->get_index("[v]");
 
1221
 
 
1222
    short** ovv_tuples = ovv_indexing->get_tuples();
 
1223
    short** oov_tuples = oov_indexing->get_tuples();
 
1224
 
 
1225
    double ***T_oovovv;
 
1226
    double ***H_ijab;
 
1227
 
 
1228
    allocate1(double**,T_oovovv,moinfo->get_nirreps());
 
1229
    allocate1(double**,H_ijab,moinfo->get_nirreps());
 
1230
 
 
1231
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1232
      // Allocate a block of T_iabjkc   
 
1233
      allocate2(double,T_oovovv[h],oov_indexing->get_pairpi(h),ovv_indexing->get_pairpi(h));
 
1234
      allocate2(double,H_ijab[h],oov_indexing->get_pairpi(h),v_indexing->get_pairpi(h));
 
1235
 
 
1236
      size_t ovv_offset = ovv_indexing->get_first(h);
 
1237
      size_t oov_offset = oov_indexing->get_first(h);
 
1238
 
 
1239
      // AAA Contribution
 
1240
 
 
1241
      // Sort this block
 
1242
      for(int kbc = 0;kbc<ovv_indexing->get_pairpi(h);kbc++){
 
1243
        int k = ovv_tuples[ovv_offset + kbc][0];
 
1244
        int b = ovv_tuples[ovv_offset + kbc][1];
 
1245
        int c = ovv_tuples[ovv_offset + kbc][2];        
 
1246
        for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1247
          int i = oov_tuples[oov_offset + ija][0];
 
1248
          int j = oov_tuples[oov_offset + ija][1];
 
1249
          int a = oov_tuples[oov_offset + ija][2];        
 
1250
          T_oovovv[h][ija][kbc] = TIJKABCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1251
        }
 
1252
      }
 
1253
      int m = oov_indexing->get_pairpi(h);
 
1254
      int n =   v_indexing->get_pairpi(h);
 
1255
      int k = ovv_indexing->get_pairpi(h);
 
1256
      if(m*n*k){
 
1257
        double alpha = 1.0;
 
1258
        double beta  = 0.0;
 
1259
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1260
                   &(WAIBC_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1261
      }
 
1262
 
 
1263
      // Sort the result
 
1264
      size_t v_offset = v_indexing->get_first(h);
 
1265
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1266
        int i = oov_tuples[oov_offset + ija][0];
 
1267
        int j = oov_tuples[oov_offset + ija][1];
 
1268
        int a = oov_tuples[oov_offset + ija][2];
 
1269
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
1270
          int abs_b = b + v_offset;
 
1271
          HIJABMatTmp->add_four_address_element(i,j,a,abs_b,0.5*H_ijab[h][ija][b]);
 
1272
          HIJABMatTmp->add_four_address_element(i,j,abs_b,a,-0.5*H_ijab[h][ija][b]);
 
1273
        }
 
1274
      }
 
1275
 
 
1276
      // Sort this block
 
1277
      for(int ica = 0;ica<ovv_indexing->get_pairpi(h);ica++){
 
1278
        int i = ovv_tuples[ovv_offset + ica][0];
 
1279
        int c = ovv_tuples[ovv_offset + ica][1];
 
1280
        int a = ovv_tuples[ovv_offset + ica][2];        
 
1281
        for(int jkb = 0;jkb<oov_indexing->get_pairpi(h);jkb++){
 
1282
          int j = oov_tuples[oov_offset + jkb][0];
 
1283
          int k = oov_tuples[oov_offset + jkb][1];
 
1284
          int b = oov_tuples[oov_offset + jkb][2];        
 
1285
          T_oovovv[h][jkb][ica] = TiJKaBCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1286
        }
 
1287
      }
 
1288
 
 
1289
      m = oov_indexing->get_pairpi(h);
 
1290
      n =   v_indexing->get_pairpi(h);
 
1291
      k = ovv_indexing->get_pairpi(h);
 
1292
      if(m*n*k){
 
1293
        double alpha = 1.0;
 
1294
        double beta  = 0.0;
 
1295
        C_DGEMM_22(m, n, k, alpha,&(T_oovovv[h][0][0]), k, 
 
1296
                   &(WAiBc_matrix[h][0][0]), k, beta, &(H_ijab[h][0][0]),n);
 
1297
      }
 
1298
 
 
1299
      // Sort the result
 
1300
      v_offset = v_indexing->get_first(h);
 
1301
      for(int ija = 0;ija<oov_indexing->get_pairpi(h);ija++){
 
1302
        int i = oov_tuples[oov_offset + ija][0];
 
1303
        int j = oov_tuples[oov_offset + ija][1];
 
1304
        int a = oov_tuples[oov_offset + ija][2];
 
1305
        for(int b = 0;b<v_indexing->get_pairpi(h);b++){
 
1306
          int abs_b = b + v_offset;
 
1307
          HIJABMatTmp->add_four_address_element(i,j,a,abs_b,H_ijab[h][ija][b]);
 
1308
          HIJABMatTmp->add_four_address_element(i,j,abs_b,a,-H_ijab[h][ija][b]);
 
1309
        }
 
1310
      }
 
1311
      // Deallocate the memory for the block
 
1312
      release2(H_ijab[h]);
 
1313
      release2(T_oovovv[h]);
 
1314
    }
 
1315
    release1(H_ijab);
 
1316
    release1(T_oovovv);
 
1317
  }
 
1318
//   blas->print("t2_test[OO][VV]{u}");
 
1319
//   blas->solve("ERROR{u} = 1000000.0 t2_test[OO][VV]{u} . t2_test[OO][VV]{u}");
 
1320
//   blas->print("ERROR{u}");
 
1321
}
 
1322
 
 
1323
/**
 
1324
 * @brief Computes the contraction
 
1325
 * \f[ \sum_{me} t_{ijm}^{abe} F''_{me} \f]
 
1326
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
1327
 * This algorithm hence performs
 
1328
 * \f[ \sum_{me} t_{ijm}^{abe} F''_{me} + \sum_{ME} t_{ijM}^{abE} F''_{ME} \rightarrow  \{ \bar{H}_{ij}^{ab} \} \f]
 
1329
 */
 
1330
void CCMRCC::build_t2_ijab_amplitudes_triples_diagram3()
 
1331
{
 
1332
  // Loop over references
 
1333
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
1334
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1335
 
 
1336
    // Grab the temporary matrices
 
1337
    CCMatTmp  HijabMatTmp   = blas->get_MatTmp("t2_eqns[oo][vv]",unique_ref,none);
 
1338
    CCMatTmp  TijkabcMatTmp = blas->get_MatTmp("t3[ooo][vvv]",unique_ref,none);
 
1339
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
1340
    CCMatTmp  FmeMatTmp     = blas->get_MatTmp("F2_me[o][v]",unique_ref,none);
 
1341
    CCMatTmp  FMEMatTmp     = blas->get_MatTmp("F2_ME[O][V]",unique_ref,none);
 
1342
 
 
1343
    // Grab the indexing for t3[ijk][abc]
 
1344
    short**   ij_tuples = HijabMatTmp->get_left()->get_tuples();
 
1345
    short**   ab_tuples = HijabMatTmp->get_right()->get_tuples();
 
1346
    short**   m_tuples  = FmeMatTmp->get_left()->get_tuples();
 
1347
    short**   e_tuples  = FmeMatTmp->get_right()->get_tuples();
 
1348
 
 
1349
    double*** Tijkabc_matrix = TijkabcMatTmp->get_matrix();
 
1350
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
1351
    double*** Hijab_matrix   = HijabMatTmp->get_matrix();
 
1352
    double*** Fme_matrix     = FmeMatTmp->get_matrix();
 
1353
    double*** FME_matrix     = FMEMatTmp->get_matrix();
 
1354
    CCIndex*  ijkIndex       = blas->get_index("[ooo]");
 
1355
    CCIndex*  abcIndex       = blas->get_index("[vvv]");
 
1356
 
 
1357
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1358
      size_t ij_offset = HijabMatTmp->get_left()->get_first(h);
 
1359
      size_t ab_offset = HijabMatTmp->get_right()->get_first(h);
 
1360
      for(int ab = 0;ab <HijabMatTmp->get_right_pairpi(h);ab++){
 
1361
        int a = ab_tuples[ab_offset + ab][0];
 
1362
        int b = ab_tuples[ab_offset + ab][1];
 
1363
        for(int ij = 0;ij<HijabMatTmp->get_left_pairpi(h);ij++){
 
1364
          int i = ij_tuples[ij_offset + ij][0];
 
1365
          int j = ij_tuples[ij_offset + ij][1];
 
1366
          for(int m_sym =0; m_sym < moinfo->get_nirreps();m_sym++){
 
1367
            size_t m_offset = FmeMatTmp->get_left()->get_first(m_sym);
 
1368
            size_t e_offset = FmeMatTmp->get_right()->get_first(m_sym);
 
1369
            for(int e = 0;e <FmeMatTmp->get_right_pairpi(m_sym);e++){
 
1370
              int e_abs = e + e_offset;
 
1371
              size_t abe  = abcIndex->get_tuple_index(a,b,e_abs);
 
1372
              int abe_sym = abcIndex->get_tuple_irrep(a,b,e_abs);
 
1373
              for(int m = 0;m <FmeMatTmp->get_left_pairpi(m_sym);m++){
 
1374
                int m_abs = m + m_offset;
 
1375
                size_t ijm  = ijkIndex->get_tuple_index(i,j,m_abs);
 
1376
                Hijab_matrix[h][ij][ab] += Tijkabc_matrix[abe_sym][ijm][abe] * Fme_matrix[m_sym][m][e];
 
1377
                Hijab_matrix[h][ij][ab] += TijKabC_matrix[abe_sym][ijm][abe] * FME_matrix[m_sym][m][e];
 
1378
              }
 
1379
            }
 
1380
          }
 
1381
        }
 
1382
      }
 
1383
    }
 
1384
  }
 
1385
//   blas->print("t2_eqns[oo][vv]{u}");
 
1386
//   blas->solve("ERROR{u} = 1000000.0 t2_eqns[oo][vv]{u} . t2_eqns[oo][vv]{u}");
 
1387
//   blas->print("ERROR{u}");
 
1388
}
 
1389
 
 
1390
/**
 
1391
 * @brief Computes the contraction
 
1392
 * \f[ \sum_{me} t_{ijm}^{abe} F''_{me} \f]
 
1393
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
1394
 * This algorithm hence performs
 
1395
 * \f[ \sum_{me} t_{imJ}^{aeB} F''_{me} + \sum_{ME} t_{iMJ}^{aEB} F''_{ME} \rightarrow  \{ \bar{H}_{iJ}^{aB} \} \f]
 
1396
 */
 
1397
void CCMRCC::build_t2_iJaB_amplitudes_triples_diagram3()
 
1398
{
 
1399
  // Loop over references
 
1400
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
1401
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1402
 
 
1403
    // Grab the temporary matrices
 
1404
    CCMatTmp  HiJaBMatTmp   = blas->get_MatTmp("t2_eqns[oO][vV]",unique_ref,none);
 
1405
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
1406
    CCMatTmp  TiJKaBCMatTmp = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
1407
    CCMatTmp  FmeMatTmp     = blas->get_MatTmp("F2_me[o][v]",unique_ref,none);
 
1408
    CCMatTmp  FMEMatTmp     = blas->get_MatTmp("F2_ME[O][V]",unique_ref,none);
 
1409
 
 
1410
    // Grab the indexing for t3[ijk][abc]
 
1411
    short**   ij_tuples = HiJaBMatTmp->get_left()->get_tuples();
 
1412
    short**   ab_tuples = HiJaBMatTmp->get_right()->get_tuples();
 
1413
    short**   m_tuples  = FmeMatTmp->get_left()->get_tuples();
 
1414
    short**   e_tuples  = FmeMatTmp->get_right()->get_tuples();
 
1415
 
 
1416
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
1417
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
1418
    double*** HiJaB_matrix   = HiJaBMatTmp->get_matrix();
 
1419
    double*** Fme_matrix     = FmeMatTmp->get_matrix();
 
1420
    double*** FME_matrix     = FMEMatTmp->get_matrix();
 
1421
    CCIndex*  ijkIndex       = blas->get_index("[ooo]");
 
1422
    CCIndex*  abcIndex       = blas->get_index("[vvv]");
 
1423
 
 
1424
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1425
      size_t ij_offset = HiJaBMatTmp->get_left()->get_first(h);
 
1426
      size_t ab_offset = HiJaBMatTmp->get_right()->get_first(h);
 
1427
      for(int ab = 0;ab <HiJaBMatTmp->get_right_pairpi(h);ab++){
 
1428
        int a = ab_tuples[ab_offset + ab][0];
 
1429
        int b = ab_tuples[ab_offset + ab][1];
 
1430
        for(int ij = 0;ij<HiJaBMatTmp->get_left_pairpi(h);ij++){
 
1431
          int i = ij_tuples[ij_offset + ij][0];
 
1432
          int j = ij_tuples[ij_offset + ij][1];
 
1433
          for(int m_sym =0; m_sym < moinfo->get_nirreps();m_sym++){
 
1434
            size_t m_offset = FmeMatTmp->get_left()->get_first(m_sym);
 
1435
            size_t e_offset = FmeMatTmp->get_right()->get_first(m_sym);
 
1436
            for(int e = 0;e <FmeMatTmp->get_right_pairpi(m_sym);e++){
 
1437
              int e_abs = e + e_offset;
 
1438
              size_t aeb  = abcIndex->get_tuple_index(a,e_abs,b);
 
1439
              int aeb_sym = abcIndex->get_tuple_irrep(a,e_abs,b);
 
1440
              for(int m = 0;m <FmeMatTmp->get_left_pairpi(m_sym);m++){
 
1441
                int m_abs = m + m_offset;
 
1442
                size_t imj  = ijkIndex->get_tuple_index(i,m_abs,j);
 
1443
                HiJaB_matrix[h][ij][ab] += TijKabC_matrix[aeb_sym][imj][aeb] * Fme_matrix[m_sym][m][e];
 
1444
                HiJaB_matrix[h][ij][ab] += TiJKaBC_matrix[aeb_sym][imj][aeb] * FME_matrix[m_sym][m][e];
 
1445
              }
 
1446
            }
 
1447
          }
 
1448
        }
 
1449
      }
 
1450
    }
 
1451
  }
 
1452
//   blas->print("t2_test[oO][vV]{u}");
 
1453
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oO][vV]{u} . t2_test[oO][vV]{u}");
 
1454
//   blas->print("ERROR{u}");
 
1455
}
 
1456
 
 
1457
/**
 
1458
 * @brief Computes the contraction
 
1459
 * \f[ \sum_{me} t_{ijm}^{abe} F''_{me} \f]
 
1460
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
1461
 * This algorithm hence performs
 
1462
 * \f[ \sum_{me} t_{mIJ}^{eAB} F''_{me} + \sum_{ME} t_{MIJ}^{EAB} F''_{ME} \rightarrow  \{ \bar{H}_{IJ}^{AB} \} \f]
 
1463
 */
 
1464
void CCMRCC::build_t2_IJAB_amplitudes_triples_diagram3()
 
1465
{
 
1466
  // Loop over references
 
1467
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
1468
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1469
 
 
1470
    // Grab the temporary matrices
 
1471
    CCMatTmp  HIJABMatTmp   = blas->get_MatTmp("t2_eqns[OO][VV]",unique_ref,none);
 
1472
    CCMatTmp  TiJKaBCMatTmp = blas->get_MatTmp("t3[oOO][vVV]",unique_ref,none);
 
1473
    CCMatTmp  TIJKABCMatTmp = blas->get_MatTmp("t3[OOO][VVV]",unique_ref,none);
 
1474
    CCMatTmp  FmeMatTmp     = blas->get_MatTmp("F2_me[o][v]",unique_ref,none);
 
1475
    CCMatTmp  FMEMatTmp     = blas->get_MatTmp("F2_ME[O][V]",unique_ref,none);
 
1476
 
 
1477
    // Grab the indexing for t3[ijk][abc]
 
1478
    short**   ij_tuples = HIJABMatTmp->get_left()->get_tuples();
 
1479
    short**   ab_tuples = HIJABMatTmp->get_right()->get_tuples();
 
1480
    short**   m_tuples  = FmeMatTmp->get_left()->get_tuples();
 
1481
    short**   e_tuples  = FmeMatTmp->get_right()->get_tuples();
 
1482
 
 
1483
    double*** TiJKaBC_matrix = TiJKaBCMatTmp->get_matrix();
 
1484
    double*** TIJKABC_matrix = TIJKABCMatTmp->get_matrix();
 
1485
    double*** HIJAB_matrix   = HIJABMatTmp->get_matrix();
 
1486
    double*** Fme_matrix     = FmeMatTmp->get_matrix();
 
1487
    double*** FME_matrix     = FMEMatTmp->get_matrix();
 
1488
    CCIndex*  ijkIndex       = blas->get_index("[ooo]");
 
1489
    CCIndex*  abcIndex       = blas->get_index("[vvv]");
 
1490
 
 
1491
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1492
      size_t ij_offset = HIJABMatTmp->get_left()->get_first(h);
 
1493
      size_t ab_offset = HIJABMatTmp->get_right()->get_first(h);
 
1494
      for(int ab = 0;ab <HIJABMatTmp->get_right_pairpi(h);ab++){
 
1495
        int a = ab_tuples[ab_offset + ab][0];
 
1496
        int b = ab_tuples[ab_offset + ab][1];
 
1497
        for(int ij = 0;ij<HIJABMatTmp->get_left_pairpi(h);ij++){
 
1498
          int i = ij_tuples[ij_offset + ij][0];
 
1499
          int j = ij_tuples[ij_offset + ij][1];
 
1500
          for(int m_sym =0; m_sym < moinfo->get_nirreps();m_sym++){
 
1501
            size_t m_offset = FmeMatTmp->get_left()->get_first(m_sym);
 
1502
            size_t e_offset = FmeMatTmp->get_right()->get_first(m_sym);
 
1503
            for(int e = 0;e <FmeMatTmp->get_right_pairpi(m_sym);e++){
 
1504
              int e_abs = e + e_offset;
 
1505
              size_t eab  = abcIndex->get_tuple_index(e_abs,a,b);
 
1506
              int eab_sym = abcIndex->get_tuple_irrep(e_abs,a,b);
 
1507
              for(int m = 0;m <FmeMatTmp->get_left_pairpi(m_sym);m++){
 
1508
                int m_abs = m + m_offset;
 
1509
                size_t mij  = ijkIndex->get_tuple_index(m_abs,i,j);
 
1510
                HIJAB_matrix[h][ij][ab] += TiJKaBC_matrix[eab_sym][mij][eab] * Fme_matrix[m_sym][m][e];
 
1511
                HIJAB_matrix[h][ij][ab] += TIJKABC_matrix[eab_sym][mij][eab] * FME_matrix[m_sym][m][e];
 
1512
              }
 
1513
            }
 
1514
          }
 
1515
        }
 
1516
      }
 
1517
    }
 
1518
  }
 
1519
//   blas->print("t2_test[oO][vV]{u}");
 
1520
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oO][vV]{u} . t2_test[oO][vV]{u}");
 
1521
//   blas->print("ERROR{u}");
 
1522
}
 
1523
 
 
1524
 
 
1525
/**
 
1526
 * @brief Computes the contraction
 
1527
 * \f[ -\frac{1}{2} P(ij)\sum_{mne} t_{imn}^{abe} W_{mnje} \f]
 
1528
 * as described in Fig. 1 of Phys. Chem. Chem. Phys. vol. 2, pg. 2047 (2000).
 
1529
 * This algorithm hence performs
 
1530
 * \f[ \frac{1}{2}\sum_{mne} t_{imn}^{abe} W_{mnje} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
1531
 * \f[ \sum_{mNE} t_{imN}^{abE} W_{mNjE} \rightarrow  \{ -\bar{H}_{ij}^{ab},\bar{H}_{ji}^{ab} \} \f]
 
1532
 */
 
1533
/*
 
1534
void CCMRCC::build_t2_ijab_amplitudes_triples_diagram1()
 
1535
{
 
1536
  // Loop over references
 
1537
  for(int ref=0;ref<moinfo->get_nunique();ref++){
 
1538
    int unique_ref  = moinfo->get_ref_number("u",ref);
 
1539
 
 
1540
    // Grab the temporary matrices
 
1541
    CCMatTmp  TijkabcMatTmp = blas->get_MatTmp("t3[ooo][vvv]",unique_ref,none);
 
1542
    double*** Tijkabc_matrix = TijkabcMatTmp->get_matrix();
 
1543
    CCMatTmp  TijKabCMatTmp = blas->get_MatTmp("t3[ooO][vvV]",unique_ref,none);
 
1544
    double*** TijKabC_matrix = TijKabCMatTmp->get_matrix();
 
1545
 
 
1546
    CCMatTmp  WkijaMatTmp = blas->get_MatTmp("W_kija[o][oov]",unique_ref,none);
 
1547
    CCMatTmp  WkiJAMatTmp = blas->get_MatTmp("W_kiJA[o][oOV]",unique_ref,none);
 
1548
    double*** Wkija_matrix = WkijaMatTmp->get_matrix();
 
1549
    double*** WkiJA_matrix = WkiJAMatTmp->get_matrix();
 
1550
 
 
1551
    CCMatTmp  HijabMatTmp = blas->get_MatTmp("t2_eqns[oo][vv]",unique_ref,none);
 
1552
 
 
1553
    // Grab the indexing for t3[iab][jkc]
 
1554
    CCIndex* iab_indexing = blas->get_index("[ovv]");
 
1555
    CCIndex* jkc_indexing = blas->get_index("[oov]");
 
1556
    CCIndex*   j_indexing = blas->get_index("[o]");
 
1557
 
 
1558
 
 
1559
    short** iab_tuples = iab_indexing->get_tuples();
 
1560
    short** jkc_tuples = jkc_indexing->get_tuples();
 
1561
 
 
1562
    // PART A: Sort T[ijk][abc]->T[iab][jkc]
 
1563
    double ***T_iabjkc;
 
1564
    double ***H_iabj;
 
1565
 
 
1566
    init_matrix<double**>(T_iabjkc,moinfo->get_nirreps());
 
1567
    init_matrix<double**>(H_iabj,moinfo->get_nirreps());
 
1568
 
 
1569
    for(int h =0; h < moinfo->get_nirreps();h++){
 
1570
      // Allocate a block of T_iabjkc   
 
1571
      init_matrix<double>(T_iabjkc[h],iab_indexing->get_pairpi(h),jkc_indexing->get_pairpi(h));
 
1572
      init_matrix<double>(H_iabj[h],iab_indexing->get_pairpi(h),j_indexing->get_pairpi(h));
 
1573
 
 
1574
      size_t iab_offset = iab_indexing->get_first(h);
 
1575
      size_t jkc_offset = jkc_indexing->get_first(h);
 
1576
 
 
1577
      // AAA Contribution
 
1578
 
 
1579
      // Sort this block
 
1580
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
1581
        int i = iab_tuples[iab_offset + iab][0];
 
1582
        int a = iab_tuples[iab_offset + iab][1];
 
1583
        int b = iab_tuples[iab_offset + iab][2];        
 
1584
        for(int jkc = 0;jkc<jkc_indexing->get_pairpi(h);jkc++){
 
1585
          int j = jkc_tuples[jkc_offset + jkc][0];
 
1586
          int k = jkc_tuples[jkc_offset + jkc][1];
 
1587
          int c = jkc_tuples[jkc_offset + jkc][2];
 
1588
          T_iabjkc[h][iab][jkc] = TijkabcMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1589
        }
 
1590
      }
 
1591
      int m = iab_indexing->get_pairpi(h);
 
1592
      int n = j_indexing->get_pairpi(h);
 
1593
      int k = jkc_indexing->get_pairpi(h);
 
1594
      if(m*n*k){
 
1595
        double alpha = 1.0;
 
1596
        double beta  = 0.0;
 
1597
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
1598
                   &(Wkija_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
1599
      }
 
1600
 
 
1601
      // Sort the result
 
1602
      size_t j_offset = j_indexing->get_first(h);
 
1603
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
1604
        int i = iab_tuples[iab_offset + iab][0];
 
1605
        int a = iab_tuples[iab_offset + iab][1];
 
1606
        int b = iab_tuples[iab_offset + iab][2];
 
1607
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
1608
          int abs_j = j + j_offset;
 
1609
          HijabMatTmp->add_four_address_element(i,abs_j,a,b,-0.5*H_iabj[h][iab][j]);
 
1610
          HijabMatTmp->add_four_address_element(abs_j,i,a,b,0.5*H_iabj[h][iab][j]);
 
1611
        }
 
1612
      }
 
1613
 
 
1614
      // AAB Contribution
 
1615
 
 
1616
      // Sort this block
 
1617
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
1618
        int i = iab_tuples[iab_offset + iab][0];
 
1619
        int a = iab_tuples[iab_offset + iab][1];
 
1620
        int b = iab_tuples[iab_offset + iab][2];        
 
1621
        for(int jkc = 0;jkc<jkc_indexing->get_pairpi(h);jkc++){
 
1622
          int j = jkc_tuples[jkc_offset + jkc][0];
 
1623
          int k = jkc_tuples[jkc_offset + jkc][1];
 
1624
          int c = jkc_tuples[jkc_offset + jkc][2];
 
1625
          T_iabjkc[h][iab][jkc] = TijKabCMatTmp->get_six_address_element(i,j,k,a,b,c);
 
1626
        }
 
1627
      }
 
1628
      m = iab_indexing->get_pairpi(h);
 
1629
      n = j_indexing->get_pairpi(h);
 
1630
      k = jkc_indexing->get_pairpi(h);
 
1631
      if(m*n*k){
 
1632
        double alpha = 1.0;
 
1633
        double beta  = 0.0;
 
1634
        C_DGEMM_22(m, n, k, alpha,&(T_iabjkc[h][0][0]), k, 
 
1635
                   &(WkiJA_matrix[h][0][0]), k, beta, &(H_iabj[h][0][0]),n);
 
1636
      }
 
1637
 
 
1638
      // Sort the result
 
1639
      j_offset = j_indexing->get_first(h);
 
1640
      for(int iab = 0;iab<iab_indexing->get_pairpi(h);iab++){
 
1641
        int i = iab_tuples[iab_offset + iab][0];
 
1642
        int a = iab_tuples[iab_offset + iab][1];
 
1643
        int b = iab_tuples[iab_offset + iab][2];
 
1644
        for(int j = 0;j<j_indexing->get_pairpi(h);j++){
 
1645
          int abs_j = j + j_offset;
 
1646
          HijabMatTmp->add_four_address_element(i,abs_j,a,b,-H_iabj[h][iab][j]);
 
1647
          HijabMatTmp->add_four_address_element(abs_j,i,a,b,H_iabj[h][iab][j]);
 
1648
        }
 
1649
      }
 
1650
 
 
1651
      // Deallocate the memory for the block
 
1652
      free_matrix<double>(H_iabj[h],iab_indexing->get_pairpi(h),j_indexing->get_pairpi(h));
 
1653
      free_matrix<double>(T_iabjkc[h],iab_indexing->get_pairpi(h),jkc_indexing->get_pairpi(h));
 
1654
    }
 
1655
    free_matrix<double**>(H_iabj,moinfo->get_nirreps());      
 
1656
    free_matrix<double**>(T_iabjkc,moinfo->get_nirreps());
 
1657
  }
 
1658
//   blas->print("t2_test[oo][vv]{u}");
 
1659
//   blas->solve("ERROR{u} = 1000000.0 t2_test[oo][vv]{u} . t2_test[oo][vv]{u}");
 
1660
//   blas->print("ERROR{u}");
 
1661
}
 
1662
*/
 
1663
 
 
1664
}} /* End Namespaces */