← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/Gciab.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/ccdensity/Gciab.cc
 
1
/*! \file
 
2
    \ingroup CCDENSITY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <strings.h>
 
7
#include <string.h>
 
8
#include <libdpd/dpd.h>
 
9
#include "MOInfo.h"
 
10
#include "Params.h"
 
11
#include "Frozen.h"
 
12
#define EXTERN
 
13
#include "globals.h"
 
14
 
 
15
namespace psi { namespace ccdensity {
 
16
 
 
17
void Gciab(void)
 
18
{
 
19
  int h, nirreps, a, b, c, i, A, B, C, I, Asym, Bsym, Csym, Isym, row, col;
 
20
  double value;
 
21
  dpdfile2 L1, T1, g;
 
22
  dpdbuf4 G, L, T, Z, Z1, Z2, V;
 
23
  double factor=0.0;
 
24
 
 
25
  nirreps = moinfo.nirreps;
 
26
 
 
27
  if(params.ref == 0) { /** RHF **/
 
28
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
29
    /* t(M,C) L(Mi,Ab) */
 
30
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LIjAb");
 
31
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
32
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, 0.0);
 
33
    dpd_file2_close(&T1);
 
34
    dpd_buf4_close(&L);
 
35
    /* l(M,C) Tau(Mi,Ab) */
 
36
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
37
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "LIA");
 
38
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
39
    dpd_file2_close(&L1);
 
40
    dpd_buf4_close(&T);
 
41
    dpd_buf4_close(&G);
 
42
    /* t(i,e) L(Mn,Ce) --> Z(Mn,Ci) */
 
43
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 11, 0, 11, 0, "Z(Mn,Ci)");
 
44
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LIjAb");
 
45
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
46
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
47
    dpd_file2_close(&T1);
 
48
    dpd_buf4_close(&L);
 
49
    /* -Z(Mn,Ci) Tau(Mn,Ab) --> G(Ci,Ab) */
 
50
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
51
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
52
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
53
    dpd_buf4_close(&T);
 
54
    dpd_buf4_close(&Z);
 
55
    dpd_buf4_close(&G);
 
56
    /* - V(iA,mC) T(m,b) --> Z(iA,bC) */
 
57
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 5, 10, 5, 0, "Z(iA,bC)");
 
58
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViAjB");
 
59
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
60
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
61
    dpd_file2_close(&T1);
 
62
    dpd_buf4_close(&V);
 
63
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 10, 5, "Z(iC,bA)");
 
64
    dpd_buf4_close(&Z);
 
65
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(iC,bA)");
 
66
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 11, 5, "Z(Ci,bA)");
 
67
    dpd_buf4_close(&Z);
 
68
    dpd_buf4_init(&Z, CC_TMP2, 0, 11, 5, 11, 5, 0, "Z(Ci,bA)");
 
69
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 11, 5, "Z(Ci,Ab)");
 
70
    dpd_buf4_close(&Z);
 
71
    /* V(ib,MC) T(M,A) --> Z(ib,AC) */
 
72
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(ib,AC)");
 
73
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViaJB");
 
74
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
75
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, 1.0, 0.0);
 
76
    dpd_file2_close(&T1);
 
77
    dpd_buf4_close(&V);
 
78
    dpd_buf4_sort(&Z, CC_TMP2, psrq, 10, 5, "Z(iC,Ab)");
 
79
    dpd_buf4_close(&Z);
 
80
    dpd_buf4_init(&Z, CC_TMP2, 0, 10, 5, 10, 5, 0, "Z(iC,Ab)");
 
81
    dpd_buf4_sort(&Z, CC_TMP1, qprs, 11, 5, "Z(Ci,Ab)");
 
82
    dpd_buf4_close(&Z);
 
83
    /* Z1(Ci,AB) + Z1(Ci,AB) --> G(Ci,AB) */
 
84
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 5, 0, "Z(Ci,Ab)");
 
85
    dpd_buf4_init(&Z2, CC_TMP1, 0, 11, 5, 11, 5, 0, "Z(Ci,Ab)");
 
86
    dpd_buf4_axpy(&Z1, &Z2, 1.0);
 
87
    dpd_buf4_close(&Z1);
 
88
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
89
    dpd_buf4_axpy(&Z2, &G, 1.0);
 
90
    dpd_buf4_close(&Z2);
 
91
    dpd_buf4_close(&G);
 
92
 
 
93
    /* g(C,A) T(i,b) --> G(Ci,Ab) */
 
94
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "GAE");
 
95
    dpd_file2_mat_init(&g);
 
96
    dpd_file2_mat_rd(&g);
 
97
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
98
    dpd_file2_mat_init(&T1);
 
99
    dpd_file2_mat_rd(&T1);
 
100
 
 
101
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
102
  
 
103
    for(h=0; h < nirreps; h++) {
 
104
      dpd_buf4_mat_irrep_init(&G, h); 0,
 
105
                                        dpd_buf4_mat_irrep_rd(&G, h);
 
106
 
 
107
      for(row=0; row < G.params->rowtot[h]; row++) {
 
108
        c = G.params->roworb[h][row][0];
 
109
        i = G.params->roworb[h][row][1];
 
110
        for(col=0; col < G.params->coltot[h]; col++) {
 
111
          a = G.params->colorb[h][col][0];
 
112
          b = G.params->colorb[h][col][1];
 
113
 
 
114
          value = 0.0;
 
115
 
 
116
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
117
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
118
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
119
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
120
              
 
121
          if((Csym==Asym) && (Isym==Bsym))
 
122
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
123
 
 
124
          G.matrix[h][row][col] -= value;
 
125
        }
 
126
      }
 
127
 
 
128
      dpd_buf4_mat_irrep_wrt(&G, h);
 
129
      dpd_buf4_mat_irrep_close(&G, h);
 
130
    }
 
131
    dpd_buf4_scm(&G, 0.5);
 
132
    dpd_buf4_close(&G);
 
133
  
 
134
    dpd_file2_mat_close(&g);
 
135
    dpd_file2_close(&g);
 
136
    dpd_file2_mat_close(&T1);
 
137
    dpd_file2_close(&T1);
 
138
  }
 
139
  else if(params.ref == 1) { /** ROHF **/
 
140
 
 
141
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "GCIAB");
 
142
    /* t(M,C) L(MI,AB) */
 
143
    dpd_buf4_init(&L, CC_GLG, 0, 0, 7, 2, 7, 0, "LIJAB");
 
144
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
145
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, 0.0);
 
146
    dpd_file2_close(&T1);
 
147
    dpd_buf4_close(&L);
 
148
    /* l(M,C) Tau(MI,AB) */
 
149
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauIJAB");
 
150
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "LIA");
 
151
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
152
    dpd_file2_close(&L1);
 
153
    dpd_buf4_close(&T);
 
154
    dpd_buf4_close(&G);
 
155
    /* t(I,E) L(MN,CE) --> Z(MN,CI) */
 
156
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 11, 2, 11, 0, "Z(MN,CI)");
 
157
    dpd_buf4_init(&L, CC_GLG, 0, 2, 5, 2, 7, 0, "LIJAB");
 
158
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
159
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
160
    dpd_file2_close(&T1);
 
161
    dpd_buf4_close(&L);
 
162
    /* -Z(MN,CI) Tau(MN,AB) --> G(CI,AB) */
 
163
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "GCIAB");
 
164
    dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
 
165
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
166
    dpd_buf4_close(&T);
 
167
    dpd_buf4_close(&Z);
 
168
    dpd_buf4_close(&G);
 
169
    /* - V(IA,MC) T(M,B) --> Z(IA,BC) */
 
170
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 5, 10, 5, 0, "Z(IA,BC)");
 
171
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAJB");
 
172
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
173
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
174
    dpd_file2_close(&T1);
 
175
    dpd_buf4_close(&V);
 
176
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 10, 5, "Z(IC,BA)");
 
177
    dpd_buf4_close(&Z);
 
178
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(IC,BA)");
 
179
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 11, 5, "Z(CI,BA)");
 
180
    dpd_buf4_close(&Z);
 
181
    dpd_buf4_init(&Z, CC_TMP2, 0, 11, 5, 11, 5, 0, "Z(CI,BA)");
 
182
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 11, 5, "Z(CI,AB)");
 
183
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 5, 0, "Z(CI,AB)");
 
184
    dpd_buf4_axpy(&Z, &Z1, -1.0);
 
185
    dpd_buf4_close(&Z);
 
186
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 7, 0, "GCIAB");
 
187
    dpd_buf4_axpy(&Z1, &G, 1.0);
 
188
    dpd_buf4_close(&Z1);
 
189
    dpd_buf4_close(&G);
 
190
 
 
191
    /* - ( g(C,A) T(I,B) - g(C,B) T(I,A) ) --> G(CI,AB) */
 
192
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "GAE");
 
193
    dpd_file2_mat_init(&g);
 
194
    dpd_file2_mat_rd(&g);
 
195
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
196
    dpd_file2_mat_init(&T1);
 
197
    dpd_file2_mat_rd(&T1);
 
198
 
 
199
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "GCIAB");
 
200
  
 
201
    for(h=0; h < nirreps; h++) {
 
202
      dpd_buf4_mat_irrep_init(&G, h); 0,
 
203
                                        dpd_buf4_mat_irrep_rd(&G, h);
 
204
 
 
205
      for(row=0; row < G.params->rowtot[h]; row++) {
 
206
        c = G.params->roworb[h][row][0];
 
207
        i = G.params->roworb[h][row][1];
 
208
        for(col=0; col < G.params->coltot[h]; col++) {
 
209
          a = G.params->colorb[h][col][0];
 
210
          b = G.params->colorb[h][col][1];
 
211
 
 
212
          value = 0.0;
 
213
 
 
214
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
215
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
216
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
217
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
218
              
 
219
          if((Csym==Asym) && (Isym==Bsym))
 
220
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
221
 
 
222
          B = g.params->colidx[b];  A = T1.params->colidx[a];
 
223
          Bsym = g.params->qsym[b];  Asym = T1.params->qsym[a];
 
224
              
 
225
          if((Csym==Bsym) && (Isym==Asym))
 
226
            value -= g.matrix[Csym][C][B] * T1.matrix[Isym][I][A];
 
227
 
 
228
          G.matrix[h][row][col] -= value;
 
229
        }
 
230
      }
 
231
 
 
232
      dpd_buf4_mat_irrep_wrt(&G, h);
 
233
      dpd_buf4_mat_irrep_close(&G, h);
 
234
    }
 
235
    dpd_buf4_scm(&G, 0.5);
 
236
    dpd_buf4_close(&G);
 
237
  
 
238
    dpd_file2_mat_close(&g);
 
239
    dpd_file2_close(&g);
 
240
    dpd_file2_mat_close(&T1);
 
241
    dpd_file2_close(&T1);
 
242
 
 
243
 
 
244
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "Gciab");
 
245
    /* t(m,c) L(mi,ab) */
 
246
    dpd_buf4_init(&L, CC_GLG, 0, 0, 7, 2, 7, 0, "Lijab");
 
247
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
248
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, 0.0);
 
249
    dpd_file2_close(&T1);
 
250
    dpd_buf4_close(&L);
 
251
    /* l(m,c) Tau(mi,ab) */
 
252
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauijab");
 
253
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "Lia");
 
254
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
255
    dpd_file2_close(&L1);
 
256
    dpd_buf4_close(&T);
 
257
    dpd_buf4_close(&G);
 
258
    /* t(i,e) L(mn,ce) --> Z(mn,ci) */
 
259
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 11, 2, 11, 0, "Z(mn,ci)");
 
260
    dpd_buf4_init(&L, CC_GLG, 0, 2, 5, 2, 7, 0, "Lijab");
 
261
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
262
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
263
    dpd_file2_close(&T1);
 
264
    dpd_buf4_close(&L);
 
265
    /* -Z(mn,ci) Tau(mn,ab) --> G(ci,ab) */
 
266
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "Gciab");
 
267
    dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauijab");
 
268
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
269
    dpd_buf4_close(&T);
 
270
    dpd_buf4_close(&Z);
 
271
    dpd_buf4_close(&G);
 
272
    /* - V(ia,mc) T(m,b) --> Z(ia,bc) */
 
273
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 5, 10, 5, 0, "Z(ia,bc)");
 
274
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "Viajb");
 
275
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
276
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
277
    dpd_file2_close(&T1);
 
278
    dpd_buf4_close(&V);
 
279
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 10, 5, "Z(ic,ba)");
 
280
    dpd_buf4_close(&Z);
 
281
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(ic,ba)");
 
282
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 11, 5, "Z(ci,ba)");
 
283
    dpd_buf4_close(&Z);
 
284
    dpd_buf4_init(&Z, CC_TMP2, 0, 11, 5, 11, 5, 0, "Z(ci,ba)");
 
285
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 11, 5, "Z(ci,ab)");
 
286
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 5, 0, "Z(ci,ab)");
 
287
    dpd_buf4_axpy(&Z, &Z1, -1.0);
 
288
    dpd_buf4_close(&Z);
 
289
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 7, 0, "Gciab");
 
290
    dpd_buf4_axpy(&Z1, &G, 1.0);
 
291
    dpd_buf4_close(&Z1);
 
292
    dpd_buf4_close(&G);
 
293
 
 
294
    /* - ( g(c,a) T(i,b) - g(c,b) T(i,a) ) --> G(ci,ab) */
 
295
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "Gae");
 
296
    dpd_file2_mat_init(&g);
 
297
    dpd_file2_mat_rd(&g);
 
298
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
299
    dpd_file2_mat_init(&T1);
 
300
    dpd_file2_mat_rd(&T1);
 
301
 
 
302
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 7, 11, 7, 0, "Gciab");
 
303
  
 
304
    for(h=0; h < nirreps; h++) {
 
305
      dpd_buf4_mat_irrep_init(&G, h); 0,
 
306
                                        dpd_buf4_mat_irrep_rd(&G, h);
 
307
 
 
308
      for(row=0; row < G.params->rowtot[h]; row++) {
 
309
        c = G.params->roworb[h][row][0];
 
310
        i = G.params->roworb[h][row][1];
 
311
        for(col=0; col < G.params->coltot[h]; col++) {
 
312
          a = G.params->colorb[h][col][0];
 
313
          b = G.params->colorb[h][col][1];
 
314
 
 
315
          value = 0.0;
 
316
 
 
317
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
318
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
319
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
320
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
321
              
 
322
          if((Csym==Asym) && (Isym==Bsym))
 
323
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
324
 
 
325
          B = g.params->colidx[b];  A = T1.params->colidx[a];
 
326
          Bsym = g.params->qsym[b];  Asym = T1.params->qsym[a];
 
327
              
 
328
          if((Csym==Bsym) && (Isym==Asym))
 
329
            value -= g.matrix[Csym][C][B] * T1.matrix[Isym][I][A];
 
330
 
 
331
          G.matrix[h][row][col] -= value;
 
332
        }
 
333
      }
 
334
 
 
335
      dpd_buf4_mat_irrep_wrt(&G, h);
 
336
      dpd_buf4_mat_irrep_close(&G, h);
 
337
    }
 
338
    dpd_buf4_scm(&G, 0.5);
 
339
    dpd_buf4_close(&G);
 
340
  
 
341
    dpd_file2_mat_close(&g);
 
342
    dpd_file2_close(&g);
 
343
    dpd_file2_mat_close(&T1);
 
344
    dpd_file2_close(&T1);
 
345
 
 
346
 
 
347
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
348
    /* t(M,C) L(Mi,Ab) */
 
349
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LIjAb");
 
350
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
351
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, 0.0);
 
352
    dpd_file2_close(&T1);
 
353
    dpd_buf4_close(&L);
 
354
    /* l(M,C) Tau(Mi,Ab) */
 
355
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
356
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "LIA");
 
357
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
358
    dpd_file2_close(&L1);
 
359
    dpd_buf4_close(&T);
 
360
    dpd_buf4_close(&G);
 
361
    /* t(i,e) L(Mn,Ce) --> Z(Mn,Ci) */
 
362
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 11, 0, 11, 0, "Z(Mn,Ci)");
 
363
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LIjAb");
 
364
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
365
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
366
    dpd_file2_close(&T1);
 
367
    dpd_buf4_close(&L);
 
368
    /* -Z(Mn,Ci) Tau(Mn,Ab) --> G(Ci,Ab) */
 
369
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
370
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
371
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
372
    dpd_buf4_close(&T);
 
373
    dpd_buf4_close(&Z);
 
374
    dpd_buf4_close(&G);
 
375
    /* - V(iA,mC) T(m,b) --> Z(iA,bC) */
 
376
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 5, 10, 5, 0, "Z(iA,bC)");
 
377
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViAjB");
 
378
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
379
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
380
    dpd_file2_close(&T1);
 
381
    dpd_buf4_close(&V);
 
382
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 10, 5, "Z(iC,bA)");
 
383
    dpd_buf4_close(&Z);
 
384
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(iC,bA)");
 
385
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 11, 5, "Z(Ci,bA)");
 
386
    dpd_buf4_close(&Z);
 
387
    dpd_buf4_init(&Z, CC_TMP2, 0, 11, 5, 11, 5, 0, "Z(Ci,bA)");
 
388
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 11, 5, "Z(Ci,Ab)");
 
389
    dpd_buf4_close(&Z);
 
390
    /* V(ib,MC) T(M,A) --> Z(ib,AC) */
 
391
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(ib,AC)");
 
392
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViaJB");
 
393
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
394
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, 1.0, 0.0);
 
395
    dpd_file2_close(&T1);
 
396
    dpd_buf4_close(&V);
 
397
    dpd_buf4_sort(&Z, CC_TMP2, psrq, 10, 5, "Z(iC,Ab)");
 
398
    dpd_buf4_close(&Z);
 
399
    dpd_buf4_init(&Z, CC_TMP2, 0, 10, 5, 10, 5, 0, "Z(iC,Ab)");
 
400
    dpd_buf4_sort(&Z, CC_TMP1, qprs, 11, 5, "Z(Ci,Ab)");
 
401
    dpd_buf4_close(&Z);
 
402
    /* Z1(Ci,AB) + Z1(Ci,AB) --> G(Ci,AB) */
 
403
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 5, 0, "Z(Ci,Ab)");
 
404
    dpd_buf4_init(&Z2, CC_TMP1, 0, 11, 5, 11, 5, 0, "Z(Ci,Ab)");
 
405
    dpd_buf4_axpy(&Z1, &Z2, 1.0);
 
406
    dpd_buf4_close(&Z1);
 
407
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
408
    dpd_buf4_axpy(&Z2, &G, 1.0);
 
409
    dpd_buf4_close(&Z2);
 
410
    dpd_buf4_close(&G);
 
411
 
 
412
    /* g(C,A) T(i,b) --> G(Ci,Ab) */
 
413
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "GAE");
 
414
    dpd_file2_mat_init(&g);
 
415
    dpd_file2_mat_rd(&g);
 
416
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
417
    dpd_file2_mat_init(&T1);
 
418
    dpd_file2_mat_rd(&T1);
 
419
 
 
420
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
 
421
  
 
422
    for(h=0; h < nirreps; h++) {
 
423
      dpd_buf4_mat_irrep_init(&G, h); 0,
 
424
                                        dpd_buf4_mat_irrep_rd(&G, h);
 
425
 
 
426
      for(row=0; row < G.params->rowtot[h]; row++) {
 
427
        c = G.params->roworb[h][row][0];
 
428
        i = G.params->roworb[h][row][1];
 
429
        for(col=0; col < G.params->coltot[h]; col++) {
 
430
          a = G.params->colorb[h][col][0];
 
431
          b = G.params->colorb[h][col][1];
 
432
 
 
433
          value = 0.0;
 
434
 
 
435
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
436
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
437
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
438
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
439
              
 
440
          if((Csym==Asym) && (Isym==Bsym))
 
441
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
442
 
 
443
          G.matrix[h][row][col] -= value;
 
444
        }
 
445
      }
 
446
 
 
447
      dpd_buf4_mat_irrep_wrt(&G, h);
 
448
      dpd_buf4_mat_irrep_close(&G, h);
 
449
    }
 
450
    dpd_buf4_scm(&G, 0.5);
 
451
    dpd_buf4_close(&G);
 
452
  
 
453
    dpd_file2_mat_close(&g);
 
454
    dpd_file2_close(&g);
 
455
    dpd_file2_mat_close(&T1);
 
456
    dpd_file2_close(&T1);
 
457
 
 
458
 
 
459
 
 
460
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
 
461
    /* t(m,c) L(mI,aB) */
 
462
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LiJaB");
 
463
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
464
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, 0.0);
 
465
    dpd_file2_close(&T1);
 
466
    dpd_buf4_close(&L);
 
467
    /* l(m,c) Tau(mI,aB) */
 
468
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
 
469
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "Lia");
 
470
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
471
    dpd_file2_close(&L1);
 
472
    dpd_buf4_close(&T);
 
473
    dpd_buf4_close(&G);
 
474
    /* t(I,E) L(mN,cE) --> Z(mN,cI) */
 
475
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 11, 0, 11, 0, "Z(mN,cI)");
 
476
    dpd_buf4_init(&L, CC_GLG, 0, 0, 5, 0, 5, 0, "LiJaB");
 
477
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
478
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
479
    dpd_file2_close(&T1);
 
480
    dpd_buf4_close(&L);
 
481
    /* -Z(mN,cI) Tau(mN,aB) --> G(cI,aB) */
 
482
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
 
483
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
 
484
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
485
    dpd_buf4_close(&T);
 
486
    dpd_buf4_close(&Z);
 
487
    dpd_buf4_close(&G);
 
488
    /* - V(Ia,Mc) T(M,B) --> Z(Ia,Bc) */
 
489
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 5, 10, 5, 0, "Z(Ia,Bc)");
 
490
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIaJb");
 
491
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
492
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
493
    dpd_file2_close(&T1);
 
494
    dpd_buf4_close(&V);
 
495
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 10, 5, "Z(Ic,Ba)");
 
496
    dpd_buf4_close(&Z);
 
497
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(Ic,Ba)");
 
498
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 11, 5, "Z(cI,Ba)");
 
499
    dpd_buf4_close(&Z);
 
500
    dpd_buf4_init(&Z, CC_TMP2, 0, 11, 5, 11, 5, 0, "Z(cI,Ba)");
 
501
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 11, 5, "Z(cI,aB)");
 
502
    dpd_buf4_close(&Z);
 
503
    /* V(IB,mc) T(m,a) --> Z(IB,ac) */
 
504
    dpd_buf4_init(&Z, CC_TMP1, 0, 10, 5, 10, 5, 0, "Z(IB,ac)");
 
505
    dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAjb");
 
506
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
507
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, 1.0, 0.0);
 
508
    dpd_file2_close(&T1);
 
509
    dpd_buf4_close(&V);
 
510
    dpd_buf4_sort(&Z, CC_TMP2, psrq, 10, 5, "Z(Ic,aB)");
 
511
    dpd_buf4_close(&Z);
 
512
    dpd_buf4_init(&Z, CC_TMP2, 0, 10, 5, 10, 5, 0, "Z(Ic,aB)");
 
513
    dpd_buf4_sort(&Z, CC_TMP1, qprs, 11, 5, "Z(cI,aB)");
 
514
    dpd_buf4_close(&Z);
 
515
    /* Z1(cI,aB) + Z2(cI,aB) --> G(cI,aB) */
 
516
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 5, 0, "Z(cI,aB)");
 
517
    dpd_buf4_init(&Z2, CC_TMP1, 0, 11, 5, 11, 5, 0, "Z(cI,aB)");
 
518
    dpd_buf4_axpy(&Z1, &Z2, 1.0);
 
519
    dpd_buf4_close(&Z1);
 
520
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
 
521
    dpd_buf4_axpy(&Z2, &G, 1.0);
 
522
    dpd_buf4_close(&Z2);
 
523
    dpd_buf4_close(&G);
 
524
 
 
525
    /* g(c,a) T(I,B) --> G(cI,aB) */
 
526
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "Gae");
 
527
    dpd_file2_mat_init(&g);
 
528
    dpd_file2_mat_rd(&g);
 
529
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
530
    dpd_file2_mat_init(&T1);
 
531
    dpd_file2_mat_rd(&T1);
 
532
 
 
533
    dpd_buf4_init(&G, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
 
534
  
 
535
    for(h=0; h < nirreps; h++) {
 
536
      dpd_buf4_mat_irrep_init(&G, h); 0,
 
537
                                        dpd_buf4_mat_irrep_rd(&G, h);
 
538
 
 
539
      for(row=0; row < G.params->rowtot[h]; row++) {
 
540
        c = G.params->roworb[h][row][0];
 
541
        i = G.params->roworb[h][row][1];
 
542
        for(col=0; col < G.params->coltot[h]; col++) {
 
543
          a = G.params->colorb[h][col][0];
 
544
          b = G.params->colorb[h][col][1];
 
545
 
 
546
          value = 0.0;
 
547
 
 
548
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
549
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
550
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
551
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
552
              
 
553
          if((Csym==Asym) && (Isym==Bsym))
 
554
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
555
 
 
556
          G.matrix[h][row][col] -= value;
 
557
        }
 
558
      }
 
559
 
 
560
      dpd_buf4_mat_irrep_wrt(&G, h);
 
561
      dpd_buf4_mat_irrep_close(&G, h);
 
562
    }
 
563
    dpd_buf4_scm(&G, 0.5);
 
564
    dpd_buf4_close(&G);
 
565
  
 
566
    dpd_file2_mat_close(&g);
 
567
    dpd_file2_close(&g);
 
568
    dpd_file2_mat_close(&T1);
 
569
    dpd_file2_close(&T1);
 
570
  }
 
571
  else if(params.ref == 2) { /** UHF **/
 
572
 
 
573
    if(!strcmp(params.wfn,"CCSD_T") && params.dertype==1) {
 
574
      /* For CCSD(T) gradients, some density contributions are
 
575
         calculated in cctriples */
 
576
      dpd_buf4_init(&G, CC_GAMMA, 0, 20, 7, 20, 7, 0, "GIDAB");
 
577
      dpd_buf4_sort(&G, CC_GAMMA, qprs, 21, 7, "GCIAB");
 
578
      dpd_buf4_close(&G);
 
579
      dpd_buf4_init(&G, CC_GAMMA, 0, 30, 17, 30, 17, 0, "Gidab");
 
580
      dpd_buf4_sort(&G, CC_GAMMA, qprs, 31, 17, "Gciab");
 
581
      dpd_buf4_close(&G);
 
582
      dpd_buf4_init(&G, CC_GAMMA, 0, 24, 28, 24, 28, 0, "GIdAb");
 
583
      dpd_buf4_sort(&G, CC_GAMMA, qpsr, 25, 29, "GcIaB");
 
584
      dpd_buf4_close(&G);
 
585
      dpd_buf4_init(&G, CC_GAMMA, 0, 27, 29, 27, 29, 0, "GiDaB");
 
586
      dpd_buf4_sort(&G, CC_GAMMA, qpsr, 26, 28, "GCiAb");
 
587
      dpd_buf4_close(&G);
 
588
 
 
589
      factor = 1.0;
 
590
    }
 
591
 
 
592
    dpd_buf4_init(&G, CC_GAMMA, 0, 21, 7, 21, 7, 0, "GCIAB");
 
593
    /* t(M,C) L(MI,AB) */
 
594
    dpd_buf4_init(&L, CC_GLG, 0, 0, 7, 2, 7, 0, "LIJAB");
 
595
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
596
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, factor);
 
597
    dpd_file2_close(&T1);
 
598
    dpd_buf4_close(&L);
 
599
    /* l(M,C) Tau(MI,AB) */
 
600
    dpd_buf4_init(&T, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauIJAB");
 
601
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "LIA");
 
602
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
603
    dpd_file2_close(&L1);
 
604
    dpd_buf4_close(&T);
 
605
    dpd_buf4_close(&G);
 
606
    /* t(I,E) L(MN,CE) --> Z(MN,CI) */
 
607
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 21, 2, 21, 0, "Z(MN,CI)");
 
608
    dpd_buf4_init(&L, CC_GLG, 0, 2, 5, 2, 7, 0, "LIJAB");
 
609
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
610
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
611
    dpd_file2_close(&T1);
 
612
    dpd_buf4_close(&L);
 
613
    /* -Z(MN,CI) Tau(MN,AB) --> G(CI,AB) */
 
614
    dpd_buf4_init(&G, CC_GAMMA, 0, 21, 7, 21, 7, 0, "GCIAB");
 
615
    dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
 
616
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
617
    dpd_buf4_close(&T);
 
618
    dpd_buf4_close(&Z);
 
619
    dpd_buf4_close(&G);
 
620
    /* - V(IA,MC) T(M,B) --> Z(IA,BC) */
 
621
    dpd_buf4_init(&Z, CC_TMP0, 0, 20, 5, 20, 5, 0, "Z(IA,BC)");
 
622
    dpd_buf4_init(&V, CC_MISC, 0, 20, 20, 20, 20, 0, "VIAJB");
 
623
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
624
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
625
    dpd_file2_close(&T1);
 
626
    dpd_buf4_close(&V);
 
627
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 20, 5, "Z(IC,BA)");
 
628
    dpd_buf4_close(&Z);
 
629
    dpd_buf4_init(&Z, CC_TMP1, 0, 20, 5, 20, 5, 0, "Z(IC,BA)");
 
630
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 21, 5, "Z(CI,BA)");
 
631
    dpd_buf4_close(&Z);
 
632
    dpd_buf4_init(&Z, CC_TMP2, 0, 21, 5, 21, 5, 0, "Z(CI,BA)");
 
633
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 21, 5, "Z(CI,AB)");
 
634
    dpd_buf4_init(&Z1, CC_TMP0, 0, 21, 5, 21, 5, 0, "Z(CI,AB)");
 
635
    dpd_buf4_axpy(&Z, &Z1, -1.0);
 
636
    dpd_buf4_close(&Z);
 
637
    dpd_buf4_init(&G, CC_GAMMA, 0, 21, 5, 21, 7, 0, "GCIAB");
 
638
    dpd_buf4_axpy(&Z1, &G, 1.0);
 
639
    dpd_buf4_close(&Z1);
 
640
    dpd_buf4_close(&G);
 
641
 
 
642
    /* - ( g(C,A) T(I,B) - g(C,B) T(I,A) ) --> G(CI,AB) */
 
643
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "GAE");
 
644
    dpd_file2_mat_init(&g);
 
645
    dpd_file2_mat_rd(&g);
 
646
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
647
    dpd_file2_mat_init(&T1);
 
648
    dpd_file2_mat_rd(&T1);
 
649
 
 
650
    dpd_buf4_init(&G, CC_GAMMA, 0, 21, 7, 21, 7, 0, "GCIAB");
 
651
  
 
652
    for(h=0; h < nirreps; h++) {
 
653
      dpd_buf4_mat_irrep_init(&G, h);
 
654
      dpd_buf4_mat_irrep_rd(&G, h);
 
655
 
 
656
      for(row=0; row < G.params->rowtot[h]; row++) {
 
657
        c = G.params->roworb[h][row][0];
 
658
        i = G.params->roworb[h][row][1];
 
659
        for(col=0; col < G.params->coltot[h]; col++) {
 
660
          a = G.params->colorb[h][col][0];
 
661
          b = G.params->colorb[h][col][1];
 
662
 
 
663
          value = 0.0;
 
664
 
 
665
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
666
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
667
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
668
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
669
              
 
670
          if((Csym==Asym) && (Isym==Bsym))
 
671
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
672
 
 
673
          B = g.params->colidx[b];  A = T1.params->colidx[a];
 
674
          Bsym = g.params->qsym[b];  Asym = T1.params->qsym[a];
 
675
              
 
676
          if((Csym==Bsym) && (Isym==Asym))
 
677
            value -= g.matrix[Csym][C][B] * T1.matrix[Isym][I][A];
 
678
 
 
679
          G.matrix[h][row][col] -= value;
 
680
        }
 
681
      }
 
682
 
 
683
      dpd_buf4_mat_irrep_wrt(&G, h);
 
684
      dpd_buf4_mat_irrep_close(&G, h);
 
685
    }
 
686
    dpd_buf4_scm(&G, 0.5);
 
687
    dpd_buf4_close(&G);
 
688
  
 
689
    dpd_file2_mat_close(&g);
 
690
    dpd_file2_close(&g);
 
691
    dpd_file2_mat_close(&T1);
 
692
    dpd_file2_close(&T1);
 
693
 
 
694
 
 
695
    dpd_buf4_init(&G, CC_GAMMA, 0, 31, 17, 31, 17, 0, "Gciab");
 
696
    /* t(m,c) L(mi,ab) */
 
697
    dpd_buf4_init(&L, CC_GLG, 0, 10, 17, 12, 17, 0, "Lijab");
 
698
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
699
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, factor);
 
700
    dpd_file2_close(&T1);
 
701
    dpd_buf4_close(&L);
 
702
    /* l(m,c) Tau(mi,ab) */
 
703
    dpd_buf4_init(&T, CC_TAMPS, 0, 10, 17, 12, 17, 0, "tauijab");
 
704
    dpd_file2_init(&L1, CC_GLG, 0, 2, 3, "Lia");
 
705
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
706
    dpd_file2_close(&L1);
 
707
    dpd_buf4_close(&T);
 
708
    dpd_buf4_close(&G);
 
709
    /* t(i,e) L(mn,ce) --> Z(mn,ci) */
 
710
    dpd_buf4_init(&Z, CC_TMP0, 0, 12, 31, 12, 31, 0, "Z(mn,ci)");
 
711
    dpd_buf4_init(&L, CC_GLG, 0, 12, 15, 12, 17, 0, "Lijab");
 
712
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
713
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
714
    dpd_file2_close(&T1);
 
715
    dpd_buf4_close(&L);
 
716
    /* -Z(mn,ci) Tau(mn,ab) --> G(ci,ab) */
 
717
    dpd_buf4_init(&G, CC_GAMMA, 0, 31, 17, 31, 17, 0, "Gciab");
 
718
    dpd_buf4_init(&T, CC_TAMPS, 0, 12, 17, 12, 17, 0, "tauijab");
 
719
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
720
    dpd_buf4_close(&T);
 
721
    dpd_buf4_close(&Z);
 
722
    dpd_buf4_close(&G);
 
723
    /* - V(ia,mc) T(m,b) --> Z(ia,bc) */
 
724
    dpd_buf4_init(&Z, CC_TMP0, 0, 30, 15, 30, 15, 0, "Z(ia,bc)");
 
725
    dpd_buf4_init(&V, CC_MISC, 0, 30, 30, 30, 30, 0, "Viajb");
 
726
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
727
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
728
    dpd_file2_close(&T1);
 
729
    dpd_buf4_close(&V);
 
730
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 30, 15, "Z(ic,ba)");
 
731
    dpd_buf4_close(&Z);
 
732
    dpd_buf4_init(&Z, CC_TMP1, 0, 30, 15, 30, 15, 0, "Z(ic,ba)");
 
733
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 31, 15, "Z(ci,ba)");
 
734
    dpd_buf4_close(&Z);
 
735
    dpd_buf4_init(&Z, CC_TMP2, 0, 31, 15, 31, 15, 0, "Z(ci,ba)");
 
736
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 31, 15, "Z(ci,ab)");
 
737
    dpd_buf4_init(&Z1, CC_TMP0, 0, 31, 15, 31, 15, 0, "Z(ci,ab)");
 
738
    dpd_buf4_axpy(&Z, &Z1, -1.0);
 
739
    dpd_buf4_close(&Z);
 
740
    dpd_buf4_init(&G, CC_GAMMA, 0, 31, 15, 31, 17, 0, "Gciab");
 
741
    dpd_buf4_axpy(&Z1, &G, 1.0);
 
742
    dpd_buf4_close(&Z1);
 
743
    dpd_buf4_close(&G);
 
744
 
 
745
    /* - ( g(c,a) T(i,b) - g(c,b) T(i,a) ) --> G(ci,ab) */
 
746
    dpd_file2_init(&g, CC_GLG, 0, 3, 3, "Gae");
 
747
    dpd_file2_mat_init(&g);
 
748
    dpd_file2_mat_rd(&g);
 
749
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
750
    dpd_file2_mat_init(&T1);
 
751
    dpd_file2_mat_rd(&T1);
 
752
 
 
753
    dpd_buf4_init(&G, CC_GAMMA, 0, 31, 17, 31, 17, 0, "Gciab");
 
754
  
 
755
    for(h=0; h < nirreps; h++) {
 
756
      dpd_buf4_mat_irrep_init(&G, h);
 
757
      dpd_buf4_mat_irrep_rd(&G, h);
 
758
 
 
759
      for(row=0; row < G.params->rowtot[h]; row++) {
 
760
        c = G.params->roworb[h][row][0];
 
761
        i = G.params->roworb[h][row][1];
 
762
        for(col=0; col < G.params->coltot[h]; col++) {
 
763
          a = G.params->colorb[h][col][0];
 
764
          b = G.params->colorb[h][col][1];
 
765
 
 
766
          value = 0.0;
 
767
 
 
768
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
769
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
770
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
771
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
772
              
 
773
          if((Csym==Asym) && (Isym==Bsym))
 
774
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
775
 
 
776
          B = g.params->colidx[b];  A = T1.params->colidx[a];
 
777
          Bsym = g.params->qsym[b];  Asym = T1.params->qsym[a];
 
778
              
 
779
          if((Csym==Bsym) && (Isym==Asym))
 
780
            value -= g.matrix[Csym][C][B] * T1.matrix[Isym][I][A];
 
781
 
 
782
          G.matrix[h][row][col] -= value;
 
783
        }
 
784
      }
 
785
 
 
786
      dpd_buf4_mat_irrep_wrt(&G, h);
 
787
      dpd_buf4_mat_irrep_close(&G, h);
 
788
    }
 
789
    dpd_buf4_scm(&G, 0.5);
 
790
    dpd_buf4_close(&G);
 
791
  
 
792
    dpd_file2_mat_close(&g);
 
793
    dpd_file2_close(&g);
 
794
    dpd_file2_mat_close(&T1);
 
795
    dpd_file2_close(&T1);
 
796
 
 
797
 
 
798
    dpd_buf4_init(&G, CC_GAMMA, 0, 26, 28, 26, 28, 0, "GCiAb");
 
799
    /* t(M,C) L(Mi,Ab) */
 
800
    dpd_buf4_init(&L, CC_GLG, 0, 22, 28, 22, 28, 0, "LIjAb");
 
801
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
802
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, factor);
 
803
    dpd_file2_close(&T1);
 
804
    dpd_buf4_close(&L);
 
805
    /* l(M,C) Tau(Mi,Ab) */
 
806
    dpd_buf4_init(&T, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
 
807
    dpd_file2_init(&L1, CC_GLG, 0, 0, 1, "LIA");
 
808
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
809
    dpd_file2_close(&L1);
 
810
    dpd_buf4_close(&T);
 
811
    dpd_buf4_close(&G);
 
812
    /* t(i,e) L(Mn,Ce) --> Z(Mn,Ci) */
 
813
    dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z(Mn,Ci)");
 
814
    dpd_buf4_init(&L, CC_GLG, 0, 22, 28, 22, 28, 0, "LIjAb");
 
815
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
816
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
817
    dpd_file2_close(&T1);
 
818
    dpd_buf4_close(&L);
 
819
    /* -Z(Mn,Ci) Tau(Mn,Ab) --> G(Ci,Ab) */
 
820
    dpd_buf4_init(&G, CC_GAMMA, 0, 26, 28, 26, 28, 0, "GCiAb");
 
821
    dpd_buf4_init(&T, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
 
822
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
823
    dpd_buf4_close(&T);
 
824
    dpd_buf4_close(&Z);
 
825
    dpd_buf4_close(&G);
 
826
    /* - V(iA,mC) T(m,b) --> Z(iA,bC) */
 
827
    dpd_buf4_init(&Z, CC_TMP0, 0, 27, 29, 27, 29, 0, "Z(iA,bC)");
 
828
    dpd_buf4_init(&V, CC_MISC, 0, 27, 27, 27, 27, 0, "ViAjB");
 
829
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
830
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
831
    dpd_file2_close(&T1);
 
832
    dpd_buf4_close(&V);
 
833
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 27, 29, "Z(iC,bA)");
 
834
    dpd_buf4_close(&Z);
 
835
    dpd_buf4_init(&Z, CC_TMP1, 0, 27, 29, 27, 29, 0, "Z(iC,bA)");
 
836
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 26, 29, "Z(Ci,bA)");
 
837
    dpd_buf4_close(&Z);
 
838
    dpd_buf4_init(&Z, CC_TMP2, 0, 26, 29, 26, 29, 0, "Z(Ci,bA)");
 
839
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 26, 28, "Z(Ci,Ab)");
 
840
    dpd_buf4_close(&Z);
 
841
    /* V(ib,MC) T(M,A) --> Z(ib,AC) */
 
842
    dpd_buf4_init(&Z, CC_TMP1, 0, 30, 5, 30, 5, 0, "Z(ib,AC)");
 
843
    dpd_buf4_init(&V, CC_MISC, 0, 30, 20, 30, 20, 0, "ViaJB");
 
844
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
845
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, 1.0, 0.0);
 
846
    dpd_file2_close(&T1);
 
847
    dpd_buf4_close(&V);
 
848
    dpd_buf4_sort(&Z, CC_TMP2, psrq, 27, 28, "Z(iC,Ab)");
 
849
    dpd_buf4_close(&Z);
 
850
    dpd_buf4_init(&Z, CC_TMP2, 0, 27, 28, 27, 28, 0, "Z(iC,Ab)");
 
851
    dpd_buf4_sort(&Z, CC_TMP1, qprs, 26, 28, "Z(Ci,Ab)");
 
852
    dpd_buf4_close(&Z);
 
853
    /* Z1(Ci,AB) + Z1(Ci,AB) --> G(Ci,AB) */
 
854
    dpd_buf4_init(&Z1, CC_TMP0, 0, 26, 28, 26, 28, 0, "Z(Ci,Ab)");
 
855
    dpd_buf4_init(&Z2, CC_TMP1, 0, 26, 28, 26, 28, 0, "Z(Ci,Ab)");
 
856
    dpd_buf4_axpy(&Z1, &Z2, 1.0);
 
857
    dpd_buf4_close(&Z1);
 
858
    dpd_buf4_init(&G, CC_GAMMA, 0, 26, 28, 26, 28, 0, "GCiAb");
 
859
    dpd_buf4_axpy(&Z2, &G, 1.0);
 
860
    dpd_buf4_close(&Z2);
 
861
    dpd_buf4_close(&G);
 
862
 
 
863
    /* g(C,A) T(i,b) --> G(Ci,Ab) */
 
864
    dpd_file2_init(&g, CC_GLG, 0, 1, 1, "GAE");
 
865
    dpd_file2_mat_init(&g);
 
866
    dpd_file2_mat_rd(&g);
 
867
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
868
    dpd_file2_mat_init(&T1);
 
869
    dpd_file2_mat_rd(&T1);
 
870
 
 
871
    dpd_buf4_init(&G, CC_GAMMA, 0, 26, 28, 26, 28, 0, "GCiAb");
 
872
  
 
873
    for(h=0; h < nirreps; h++) {
 
874
      dpd_buf4_mat_irrep_init(&G, h);
 
875
      dpd_buf4_mat_irrep_rd(&G, h);
 
876
 
 
877
      for(row=0; row < G.params->rowtot[h]; row++) {
 
878
        c = G.params->roworb[h][row][0];
 
879
        i = G.params->roworb[h][row][1];
 
880
        for(col=0; col < G.params->coltot[h]; col++) {
 
881
          a = G.params->colorb[h][col][0];
 
882
          b = G.params->colorb[h][col][1];
 
883
 
 
884
          value = 0.0;
 
885
 
 
886
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
887
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
888
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
889
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
890
              
 
891
          if((Csym==Asym) && (Isym==Bsym))
 
892
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
893
 
 
894
          G.matrix[h][row][col] -= value;
 
895
        }
 
896
      }
 
897
 
 
898
      dpd_buf4_mat_irrep_wrt(&G, h);
 
899
      dpd_buf4_mat_irrep_close(&G, h);
 
900
    }
 
901
    dpd_buf4_scm(&G, 0.5);
 
902
    dpd_buf4_close(&G);
 
903
  
 
904
    dpd_file2_mat_close(&g);
 
905
    dpd_file2_close(&g);
 
906
    dpd_file2_mat_close(&T1);
 
907
    dpd_file2_close(&T1);
 
908
 
 
909
 
 
910
 
 
911
    dpd_buf4_init(&G, CC_GAMMA, 0, 25, 29, 25, 29, 0, "GcIaB");
 
912
    /* t(m,c) L(mI,aB) */
 
913
    dpd_buf4_init(&L, CC_GLG, 0, 23, 29, 23, 29, 0, "LiJaB");
 
914
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
915
    dpd_contract244(&T1, &L, &G, 0, 0, 0, 1.0, factor);
 
916
    dpd_file2_close(&T1);
 
917
    dpd_buf4_close(&L);
 
918
    /* l(m,c) Tau(mI,aB) */
 
919
    dpd_buf4_init(&T, CC_TAMPS, 0, 23, 29, 23, 29, 0, "tauiJaB");
 
920
    dpd_file2_init(&L1, CC_GLG, 0, 2, 3, "Lia");
 
921
    dpd_contract244(&L1, &T, &G, 0, 0, 0, 1.0, 1.0);
 
922
    dpd_file2_close(&L1);
 
923
    dpd_buf4_close(&T);
 
924
    dpd_buf4_close(&G);
 
925
    /* t(I,E) L(mN,cE) --> Z(mN,cI) */
 
926
    dpd_buf4_init(&Z, CC_TMP0, 0, 23, 25, 23, 25, 0, "Z(mN,cI)");
 
927
    dpd_buf4_init(&L, CC_GLG, 0, 23, 29, 23, 29, 0, "LiJaB");
 
928
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
929
    dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
 
930
    dpd_file2_close(&T1);
 
931
    dpd_buf4_close(&L);
 
932
    /* -Z(mN,cI) Tau(mN,aB) --> G(cI,aB) */
 
933
    dpd_buf4_init(&G, CC_GAMMA, 0, 25, 29, 25, 29, 0, "GcIaB");
 
934
    dpd_buf4_init(&T, CC_TAMPS, 0, 23, 29, 23, 29, 0, "tauiJaB");
 
935
    dpd_contract444(&Z, &T, &G, 1, 1, -1.0, 1.0);
 
936
    dpd_buf4_close(&T);
 
937
    dpd_buf4_close(&Z);
 
938
    dpd_buf4_close(&G);
 
939
    /* - V(Ia,Mc) T(M,B) --> Z(Ia,Bc) */
 
940
    dpd_buf4_init(&Z, CC_TMP0, 0, 24, 28, 24, 28, 0, "Z(Ia,Bc)");
 
941
    dpd_buf4_init(&V, CC_MISC, 0, 24, 24, 24, 24, 0, "VIaJb");
 
942
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
943
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, -1.0, 0.0);
 
944
    dpd_file2_close(&T1);
 
945
    dpd_buf4_close(&V);
 
946
    dpd_buf4_sort(&Z, CC_TMP1, psrq, 24, 28, "Z(Ic,Ba)");
 
947
    dpd_buf4_close(&Z);
 
948
    dpd_buf4_init(&Z, CC_TMP1, 0, 24, 28, 24, 28, 0, "Z(Ic,Ba)");
 
949
    dpd_buf4_sort(&Z, CC_TMP2, qprs, 25, 28, "Z(cI,Ba)");
 
950
    dpd_buf4_close(&Z);
 
951
    dpd_buf4_init(&Z, CC_TMP2, 0, 25, 28, 25, 28, 0, "Z(cI,Ba)");
 
952
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 25, 29, "Z(cI,aB)");
 
953
    dpd_buf4_close(&Z);
 
954
    /* V(IB,mc) T(m,a) --> Z(IB,ac) */
 
955
    dpd_buf4_init(&Z, CC_TMP1, 0, 20, 15, 20, 15, 0, "Z(IB,ac)");
 
956
    dpd_buf4_init(&V, CC_MISC, 0, 20, 30, 20, 30, 0, "VIAjb");
 
957
    dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
958
    dpd_contract244(&T1, &V, &Z, 0, 2, 1, 1.0, 0.0);
 
959
    dpd_file2_close(&T1);
 
960
    dpd_buf4_close(&V);
 
961
    dpd_buf4_sort(&Z, CC_TMP2, psrq, 24, 29, "Z(Ic,aB)");
 
962
    dpd_buf4_close(&Z);
 
963
    dpd_buf4_init(&Z, CC_TMP2, 0, 24, 29, 24, 29, 0, "Z(Ic,aB)");
 
964
    dpd_buf4_sort(&Z, CC_TMP1, qprs, 25, 29, "Z(cI,aB)");
 
965
    dpd_buf4_close(&Z);
 
966
    /* Z1(cI,aB) + Z2(cI,aB) --> G(cI,aB) */
 
967
    dpd_buf4_init(&Z1, CC_TMP0, 0, 25, 29, 25, 29, 0, "Z(cI,aB)");
 
968
    dpd_buf4_init(&Z2, CC_TMP1, 0, 25, 29, 25, 29, 0, "Z(cI,aB)");
 
969
    dpd_buf4_axpy(&Z1, &Z2, 1.0);
 
970
    dpd_buf4_close(&Z1);
 
971
    dpd_buf4_init(&G, CC_GAMMA, 0, 25, 29, 25, 29, 0, "GcIaB");
 
972
    dpd_buf4_axpy(&Z2, &G, 1.0);
 
973
    dpd_buf4_close(&Z2);
 
974
    dpd_buf4_close(&G);
 
975
 
 
976
    /* g(c,a) T(I,B) --> G(cI,aB) */
 
977
    dpd_file2_init(&g, CC_GLG, 0, 3, 3, "Gae");
 
978
    dpd_file2_mat_init(&g);
 
979
    dpd_file2_mat_rd(&g);
 
980
    dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
981
    dpd_file2_mat_init(&T1);
 
982
    dpd_file2_mat_rd(&T1);
 
983
 
 
984
    dpd_buf4_init(&G, CC_GAMMA, 0, 25, 29, 25, 29, 0, "GcIaB");
 
985
  
 
986
    for(h=0; h < nirreps; h++) {
 
987
      dpd_buf4_mat_irrep_init(&G, h);
 
988
      dpd_buf4_mat_irrep_rd(&G, h);
 
989
 
 
990
      for(row=0; row < G.params->rowtot[h]; row++) {
 
991
        c = G.params->roworb[h][row][0];
 
992
        i = G.params->roworb[h][row][1];
 
993
        for(col=0; col < G.params->coltot[h]; col++) {
 
994
          a = G.params->colorb[h][col][0];
 
995
          b = G.params->colorb[h][col][1];
 
996
 
 
997
          value = 0.0;
 
998
 
 
999
          C = g.params->rowidx[c];  I = T1.params->rowidx[i];
 
1000
          Csym = g.params->psym[c]; Isym = T1.params->psym[i];
 
1001
          A = g.params->colidx[a];  B = T1.params->colidx[b];
 
1002
          Asym = g.params->qsym[a];  Bsym = T1.params->qsym[b];
 
1003
              
 
1004
          if((Csym==Asym) && (Isym==Bsym))
 
1005
            value += g.matrix[Csym][C][A] * T1.matrix[Isym][I][B];
 
1006
 
 
1007
          G.matrix[h][row][col] -= value;
 
1008
        }
 
1009
      }
 
1010
 
 
1011
      dpd_buf4_mat_irrep_wrt(&G, h);
 
1012
      dpd_buf4_mat_irrep_close(&G, h);
 
1013
    }
 
1014
    dpd_buf4_scm(&G, 0.5);
 
1015
    dpd_buf4_close(&G);
 
1016
  
 
1017
    dpd_file2_mat_close(&g);
 
1018
    dpd_file2_close(&g);
 
1019
    dpd_file2_mat_close(&T1);
 
1020
    dpd_file2_close(&T1);
 
1021
 
 
1022
  }
 
1023
}
 
1024
 
 
1025
}} // namespace psi::ccdensity