← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/lib/libdpd/cc3_sigma_UHF.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/lib/libdpd/cc3_sigma_UHF.cc
 
1
/*! \file
 
2
    \ingroup DPD
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cstring>
 
8
#include <cmath>
 
9
#include <libciomr/libciomr.h>
 
10
#include <libqt/qt.h>
 
11
#include <libdpd/dpd.h>
 
12
#include <ccfiles.h>
 
13
 
 
14
extern "C" {
 
15
 
 
16
/*
 
17
  This function computes contributions to singles and doubles of
 
18
  matrix elements of triples:
 
19
    SIA   <-- <S|(Dints)           <T|(Wmbij,Wabei) CMNEF |0> |T> / (w-wt)
 
20
    SIjAb <-- <D|(FME,WAmEf,WMnIe) <T|(Wmbij,Wabei) CMNEF |0> |T> / (w-wt)
 
21
  These are used to make X3 quantity in T3_RHF:
 
22
    CIjAb, WAbEi, WMbIj, fIJ, fAB, omega
 
23
*/
 
24
 
 
25
void cc3_sigma_UHF_AAA(dpdbuf4 *CMNEF, dpdbuf4 *WABEI, dpdbuf4 *WMBIJ,
 
26
    int do_singles, dpdbuf4 *Dints_anti, dpdfile2 *SIA, int do_doubles, dpdfile2 *FME,
 
27
    dpdbuf4 *WMAFE, dpdbuf4 *WMNIE, dpdbuf4 *SIJAB, int *aoccpi, int *aocc_off,
 
28
    int *avirtpi, int *avir_off, double omega, FILE *outfile)
 
29
{
 
30
  int h, nirreps;
 
31
  int *occ_off, *occpi, *vir_off, *virtpi;
 
32
  int Gi, Gj, Gk, Gijk, Ga, Gb, Gc;
 
33
  int i, j, k, I, J, K;
 
34
  int a, b, c, A, B, C;
 
35
  int Gij, ij, Gab, ab, Gjk, jk;
 
36
  double ***W1, **Z;
 
37
  dpdfile2 fIJ, fAB;
 
38
  int nrows, ncols, nlinks;
 
39
  int Gd, d, cd, dc, Gid, id, DD;
 
40
  int Gm, m, Gmi, mi, im, mc, M;
 
41
  int GS, GH, GU, GC, GX3;
 
42
 
 
43
  GC = CMNEF->file.my_irrep;
 
44
 
 
45
  GU = WMBIJ->file.my_irrep;
 
46
 
 
47
  if (do_singles)
 
48
    GH = Dints_anti->file.my_irrep;
 
49
  else if (do_doubles)
 
50
    GH = WMAFE->file.my_irrep;
 
51
 
 
52
  GS = SIJAB->file.my_irrep;
 
53
 
 
54
  GX3 = GU^GC;
 
55
 
 
56
  nirreps = CMNEF->params->nirreps;
 
57
  occpi = aoccpi;
 
58
  occ_off = aocc_off;
 
59
  virtpi = avirtpi;
 
60
  vir_off = avir_off;
 
61
 
 
62
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
63
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
64
 
 
65
  if (do_singles) {
 
66
    dpd_file2_mat_init(SIA);
 
67
    dpd_file2_mat_rd(SIA);
 
68
 
 
69
    for(h=0; h < nirreps; h++) {
 
70
      dpd_buf4_mat_irrep_init(Dints_anti, h);
 
71
      dpd_buf4_mat_irrep_rd(Dints_anti, h);
 
72
    }
 
73
  }
 
74
 
 
75
  if (do_doubles) {
 
76
    dpd_file2_mat_init(FME);
 
77
    dpd_file2_mat_rd(FME);
 
78
 
 
79
    for(h=0; h < nirreps; h++) {
 
80
      dpd_buf4_mat_irrep_init(SIJAB, h);
 
81
      dpd_buf4_mat_irrep_rd(SIJAB, h);
 
82
    }
 
83
    for(h=0; h < nirreps; h++) {
 
84
      dpd_buf4_mat_irrep_init(WMNIE, h);
 
85
      dpd_buf4_mat_irrep_rd(WMNIE, h);
 
86
    }
 
87
  }
 
88
 
 
89
  /* target T3 amplitudes go in here */
 
90
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
91
 
 
92
  for(Gi=0; Gi < nirreps; Gi++) {
 
93
    for(Gj=0; Gj < nirreps; Gj++) {
 
94
      Gij = Gi ^ Gj;
 
95
      for(Gk=0; Gk < nirreps; Gk++) {
 
96
        Gijk = Gi ^ Gj ^ Gk;
 
97
        Gjk = Gj ^ Gk;
 
98
 
 
99
        for(Gab=0; Gab < nirreps; Gab++) {
 
100
          /* changed */
 
101
          Gc = Gab ^ Gijk ^ GX3;
 
102
          W1[Gab] = dpd_block_matrix(WABEI->params->coltot[Gab], virtpi[Gc]);
 
103
        }
 
104
 
 
105
        for(i=0; i < occpi[Gi]; i++) {
 
106
          I = occ_off[Gi] + i;
 
107
          for(j=0; j < occpi[Gj]; j++) {
 
108
            J = occ_off[Gj] + j;
 
109
            for(k=0; k < occpi[Gk]; k++) {
 
110
              K = occ_off[Gk] + k;
 
111
 
 
112
              T3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, CMNEF, WABEI, WMBIJ, &fIJ, &fAB, 
 
113
                     occpi, occ_off, virtpi, vir_off, omega);
 
114
 
 
115
          if (do_singles) {
 
116
              /* t_KC <-- 1/4 t_IJKABC <IJ||AB> */
 
117
 
 
118
              Gc = Gk ^ GS;  /* changed */
 
119
              Gab = Gij ^ GH ;  /* changed */
 
120
 
 
121
              ij = Dints_anti->params->rowidx[I][J];
 
122
 
 
123
              nrows = Dints_anti->params->coltot[Gij^GH]; /* changed */
 
124
              ncols = virtpi[Gc];
 
125
 
 
126
              if(nrows && ncols)
 
127
                C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, Dints_anti->matrix[Gij][ij], 1,
 
128
                        1.0, SIA->matrix[Gk][k], 1);
 
129
          }
 
130
 
 
131
          if (do_doubles) {
 
132
              /* t_IJAB <-- t_IJKABC F_KC */
 
133
              Gc = Gk ^ GH;    
 
134
              Gab = Gij ^ GS; 
 
135
 
 
136
              nrows = SIJAB->params->coltot[Gij^GS];
 
137
              ncols = virtpi[Gc];
 
138
              ij = SIJAB->params->rowidx[I][J];
 
139
 
 
140
              if(nrows && ncols)
 
141
                C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][0], ncols, FME->matrix[Gk][k], 1,
 
142
                        1.0, SIJAB->matrix[Gij][ij], 1);
 
143
 
 
144
              /* t_JKDC <-- +1/2 t_IJKABC W_IDAB */
 
145
              /* t_JKCD <-- -1/2 t_IJKABC W_IDAB */
 
146
              jk = SIJAB->params->rowidx[J][K];
 
147
              for(Gd=0; Gd < nirreps; Gd++) {
 
148
                Gid = Gi ^ Gd;  
 
149
                Gab = Gid ^ GH;
 
150
                Gc = Gab ^ Gijk ^ GX3; 
 
151
 
 
152
                id = WMAFE->row_offset[Gid][I];
 
153
 
 
154
                Z = block_matrix(virtpi[Gc],virtpi[Gd]);
 
155
                WMAFE->matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMAFE->params->coltot[Gid^GH]);
 
156
                dpd_buf4_mat_irrep_rd_block(WMAFE, Gid, id, virtpi[Gd]);
 
157
 
 
158
                nrows = virtpi[Gc];
 
159
                ncols = virtpi[Gd];
 
160
                nlinks = WMAFE->params->coltot[Gid^GH];
 
161
 
 
162
                if(nrows && ncols && nlinks)
 
163
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows, 
 
164
                          WMAFE->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
165
 
 
166
                for(c=0; c < virtpi[Gc]; c++) {
 
167
                  C = vir_off[Gc] + c;
 
168
                  for(d=0; d < virtpi[Gd]; d++) {
 
169
                    DD = vir_off[Gd] + d;
 
170
                    cd = SIJAB->params->colidx[C][DD];
 
171
                    dc = SIJAB->params->colidx[DD][C];
 
172
                    SIJAB->matrix[Gjk][jk][dc] += Z[c][d];
 
173
                    SIJAB->matrix[Gjk][jk][cd] += -Z[c][d];
 
174
                  }
 
175
                }
 
176
                dpd_free_block(WMAFE->matrix[Gid], virtpi[Gd], WMAFE->params->coltot[Gid^GH]);
 
177
                free_block(Z);
 
178
              }
 
179
 
 
180
 
 
181
              /* S_MIAB <-- +1/2 t_IJKABC W_JKMC */
 
182
              /* S_IMAB <-- -1/2 t_IJKABC W_JKMC */
 
183
              jk = WMNIE->params->rowidx[J][K];
 
184
              for(Gm=0; Gm < nirreps; Gm++) {
 
185
                Gmi = Gm ^ Gi;
 
186
                Gab = Gmi ^ GS;
 
187
                Gc = Gab ^ Gijk ^ GX3;   
 
188
 
 
189
                mc = WMNIE->col_offset[Gjk][Gm];
 
190
 
 
191
                nrows = WABEI->params->coltot[Gab];
 
192
                ncols = occpi[Gm];
 
193
                nlinks = virtpi[Gc];
 
194
 
 
195
                Z = dpd_block_matrix(nrows, ncols);
 
196
 
 
197
                if(nrows && ncols && nlinks) 
 
198
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nlinks,
 
199
                          &(WMNIE->matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
200
 
 
201
                for(m=0; m < ncols; m++) {
 
202
                  M = occ_off[Gm] + m;
 
203
                  mi = SIJAB->params->rowidx[M][I];
 
204
                  im = SIJAB->params->rowidx[I][M];
 
205
                  for(ab=0; ab < nrows; ab++) {
 
206
                    SIJAB->matrix[Gmi][mi][ab] += Z[ab][m];
 
207
                    SIJAB->matrix[Gmi][im][ab] -= Z[ab][m];
 
208
                  }
 
209
                }
 
210
 
 
211
                dpd_free_block(Z, nrows, ncols);
 
212
              }
 
213
            } /* end do_doubles */
 
214
 
 
215
            } /* k */
 
216
          } /* j */
 
217
        } /* i */
 
218
 
 
219
        for(Gab=0; Gab < nirreps; Gab++) {
 
220
          /* This will need to change for non-totally-symmetric cases */
 
221
          Gc = Gab ^ Gijk;
 
222
          dpd_free_block(W1[Gab], WABEI->params->coltot[Gab], virtpi[Gc]);
 
223
        }
 
224
 
 
225
      } /* Gk */
 
226
    } /* Gj */
 
227
  } /* Gi */
 
228
 
 
229
  free(W1);
 
230
 
 
231
  dpd_file2_close(&fIJ);
 
232
  dpd_file2_close(&fAB);
 
233
 
 
234
  if (do_singles) {
 
235
 
 
236
    dpd_file2_mat_wrt(SIA);
 
237
    dpd_file2_mat_close(SIA);
 
238
 
 
239
    for(h=0; h < nirreps; h++)
 
240
      dpd_buf4_mat_irrep_close(Dints_anti, h);
 
241
  }
 
242
 
 
243
  if (do_doubles) {
 
244
 
 
245
    dpd_file2_mat_close(FME);
 
246
 
 
247
    for(h=0; h < nirreps; h++)
 
248
      dpd_buf4_mat_irrep_close(WMNIE, h);
 
249
 
 
250
    for(h=0; h < nirreps; h++) {
 
251
      dpd_buf4_mat_irrep_wrt(SIJAB, h);
 
252
      dpd_buf4_mat_irrep_close(SIJAB, h);
 
253
    }
 
254
  }
 
255
}
 
256
 
 
257
void cc3_sigma_UHF_BBB(dpdbuf4 *Cmnef, dpdbuf4 *Wabei, dpdbuf4 *Wmbij,
 
258
     int do_singles, dpdbuf4 *Dijab_anti, dpdfile2 *Sia, int do_doubles,
 
259
     dpdfile2 *Fme, dpdbuf4 *Wmafe, dpdbuf4 *Wmnie, dpdbuf4 *Sijab,
 
260
     int *boccpi, int *bocc_off, int *bvirtpi, int *bvir_off,
 
261
     double omega, FILE *outfile)
 
262
{
 
263
  int h, nirreps;
 
264
  int *occ_off, *occpi, *vir_off, *virtpi;
 
265
  int Gi, Gj, Gk, Gijk, Ga, Gb, Gc;
 
266
  int i, j, k, I, J, K;
 
267
  int a, b, c, A, B, C;
 
268
  int Gij, ij, Gab, ab, Gjk, jk;
 
269
  double ***W1, **Z;
 
270
  dpdfile2 fIJ, fAB;
 
271
  int nrows, ncols, nlinks;
 
272
  int Gd, d, cd, dc, Gid, id, DD;
 
273
  int Gm, m, Gmi, mi, im, mc, M;
 
274
  int GS, GH, GU, GC, GX3;
 
275
  
 
276
  GC = Cmnef->file.my_irrep;
 
277
  GU = Wmbij->file.my_irrep; 
 
278
 
 
279
  if (do_singles)
 
280
    GH = Dijab_anti->file.my_irrep;
 
281
  else if (do_doubles)
 
282
    GH = Wmafe->file.my_irrep;
 
283
 
 
284
  GS = Sijab->file.my_irrep;
 
285
  GX3 = GU^GC;
 
286
 
 
287
  nirreps = Cmnef->params->nirreps;
 
288
  occpi = boccpi;
 
289
  occ_off = bocc_off;
 
290
  virtpi = bvirtpi;
 
291
  vir_off = bvir_off;
 
292
 
 
293
  dpd_file2_init(&fIJ, CC_OEI, 0, 2, 2, "fij");
 
294
  dpd_file2_init(&fAB, CC_OEI, 0, 3, 3, "fab");
 
295
 
 
296
  if (do_singles) {
 
297
 
 
298
    dpd_file2_mat_init(Sia);
 
299
    dpd_file2_mat_rd(Sia);
 
300
  
 
301
    for(h=0; h < nirreps; h++) {
 
302
      dpd_buf4_mat_irrep_init(Dijab_anti, h);
 
303
      dpd_buf4_mat_irrep_rd(Dijab_anti, h);
 
304
    }
 
305
  }
 
306
 
 
307
  if (do_doubles) {
 
308
 
 
309
    dpd_file2_mat_init(Fme);
 
310
    dpd_file2_mat_rd(Fme);
 
311
 
 
312
    for(h=0; h < nirreps; h++) {
 
313
      dpd_buf4_mat_irrep_init(Sijab, h);
 
314
      dpd_buf4_mat_irrep_rd(Sijab, h);
 
315
    }
 
316
 
 
317
    for(h=0; h < nirreps; h++) {
 
318
      dpd_buf4_mat_irrep_init(Wmnie, h);
 
319
      dpd_buf4_mat_irrep_rd(Wmnie, h);
 
320
    }
 
321
  }
 
322
 
 
323
  /* target T3 amplitudes go in here */
 
324
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
325
 
 
326
  for(Gi=0; Gi < nirreps; Gi++) {
 
327
    for(Gj=0; Gj < nirreps; Gj++) {
 
328
      Gij = Gi ^ Gj;
 
329
      for(Gk=0; Gk < nirreps; Gk++) {
 
330
        Gijk = Gi ^ Gj ^ Gk;
 
331
        Gjk = Gj ^ Gk;
 
332
 
 
333
        for(Gab=0; Gab < nirreps; Gab++) {
 
334
          Gc = Gab ^ Gijk ^ GX3;
 
335
          W1[Gab] = dpd_block_matrix(Wabei->params->coltot[Gab], virtpi[Gc]);
 
336
        }
 
337
 
 
338
        for(i=0; i < occpi[Gi]; i++) {
 
339
          I = occ_off[Gi] + i;
 
340
          for(j=0; j < occpi[Gj]; j++) {
 
341
            J = occ_off[Gj] + j;
 
342
            for(k=0; k < occpi[Gk]; k++) {
 
343
              K = occ_off[Gk] + k;
 
344
 
 
345
              T3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, Cmnef, Wabei, Wmbij, &fIJ, &fAB, 
 
346
                     occpi, occ_off, virtpi, vir_off, omega);
 
347
 
 
348
          if (do_singles) {
 
349
 
 
350
              /* S_kc <-- 1/4 t_ijkabc <ij||ab> */
 
351
              Gc = Gk ^ GS;   
 
352
              Gab = Gij ^ GH;
 
353
 
 
354
              ij = Dijab_anti->params->rowidx[I][J];
 
355
 
 
356
              nrows = Dijab_anti->params->coltot[Gij^GH];
 
357
              ncols = virtpi[Gc];
 
358
 
 
359
              if(nrows && ncols)
 
360
                C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, Dijab_anti->matrix[Gij][ij], 1,
 
361
                        1.0, Sia->matrix[Gk][k], 1);
 
362
          }
 
363
 
 
364
          if (do_doubles) {
 
365
 
 
366
              /* S_ijab <-- t_ijkabc F_kc */
 
367
              Gc = Gk ^ GH;   
 
368
              Gab = Gij ^ GS;
 
369
 
 
370
              nrows = Sijab->params->coltot[Gij^GS];
 
371
              ncols = virtpi[Gc];
 
372
              ij = Sijab->params->rowidx[I][J];
 
373
 
 
374
              if(nrows && ncols)
 
375
                C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][0], ncols, Fme->matrix[Gk][k], 1,
 
376
                        1.0, Sijab->matrix[Gij][ij], 1);
 
377
 
 
378
              /* S_jkdc <-- 1/2 t_ijkabc W_idab */
 
379
              /* S_jkcd <-- -1/2 t_ijkabc W_idab */
 
380
              jk = Sijab->params->rowidx[J][K];
 
381
              for(Gd=0; Gd < nirreps; Gd++) {
 
382
 
 
383
                Gid = Gi ^ Gd; 
 
384
                Gab = Gid ^ GH; 
 
385
                Gc = Gab ^ Gijk ^ GX3; 
 
386
 
 
387
                id = Wmafe->row_offset[Gid][I];
 
388
 
 
389
                Z = block_matrix(virtpi[Gc],virtpi[Gd]);
 
390
                Wmafe->matrix[Gid] = dpd_block_matrix(virtpi[Gd], Wmafe->params->coltot[Gid^GH]);
 
391
                dpd_buf4_mat_irrep_rd_block(Wmafe, Gid, id, virtpi[Gd]);
 
392
 
 
393
                nrows = virtpi[Gc];
 
394
                ncols = virtpi[Gd];
 
395
                nlinks = Wmafe->params->coltot[Gid^GH];
 
396
 
 
397
                if(nrows && ncols && nlinks)
 
398
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows, 
 
399
                          Wmafe->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
400
 
 
401
                for(c=0; c < virtpi[Gc]; c++) {
 
402
                  C = vir_off[Gc] + c;
 
403
                  for(d=0; d < virtpi[Gd]; d++) {
 
404
                    DD = vir_off[Gd] + d;
 
405
                    cd = Sijab->params->colidx[C][DD];
 
406
                    dc = Sijab->params->colidx[DD][C];
 
407
                    Sijab->matrix[Gjk][jk][dc] += Z[c][d];
 
408
                    Sijab->matrix[Gjk][jk][cd] += -Z[c][d];
 
409
                  }
 
410
                }
 
411
                dpd_free_block(Wmafe->matrix[Gid], virtpi[Gd], Wmafe->params->coltot[Gid^GH]);
 
412
                free_block(Z);
 
413
              }
 
414
 
 
415
              /* S_miab <-- +1/2 t_ijkabc W_jkmc */
 
416
              /* S_imab <-- -1/2 t_ijkabc W_jkmc */
 
417
              jk = Wmnie->params->rowidx[J][K];
 
418
              for(Gm=0; Gm < nirreps; Gm++) {
 
419
                Gmi = Gm ^ Gi; 
 
420
                Gab = Gmi ^ GS; 
 
421
                Gc = Gab ^ Gijk ^ GX3;    
 
422
 
 
423
                mc = Wmnie->col_offset[Gjk][Gm];
 
424
 
 
425
                nrows = Wabei->params->coltot[Gab];
 
426
                ncols = occpi[Gm];
 
427
                nlinks = virtpi[Gc];
 
428
 
 
429
                Z = dpd_block_matrix(nrows, ncols);
 
430
 
 
431
                if(nrows && ncols && nlinks) 
 
432
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nlinks,
 
433
                          &(Wmnie->matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
434
 
 
435
                for(m=0; m < ncols; m++) {
 
436
                  M = occ_off[Gm] + m;
 
437
                  mi = Sijab->params->rowidx[M][I];
 
438
                  im = Sijab->params->rowidx[I][M];
 
439
                  for(ab=0; ab < nrows; ab++) {
 
440
                    Sijab->matrix[Gmi][mi][ab] += Z[ab][m];
 
441
                    Sijab->matrix[Gmi][im][ab] -= Z[ab][m];
 
442
                  }
 
443
                }
 
444
 
 
445
                dpd_free_block(Z, nrows, ncols);
 
446
              }
 
447
          } /* end do_doubles */
 
448
 
 
449
            } /* k */
 
450
          } /* j */
 
451
        } /* i */
 
452
 
 
453
        for(Gab=0; Gab < nirreps; Gab++) {
 
454
          /* This will need to change for non-totally-symmetric cases */
 
455
          Gc = Gab ^ Gijk;
 
456
          dpd_free_block(W1[Gab], Wabei->params->coltot[Gab], virtpi[Gc]);
 
457
        }
 
458
 
 
459
      } /* Gk */
 
460
    } /* Gj */
 
461
  } /* Gi */
 
462
 
 
463
  free(W1);
 
464
 
 
465
  dpd_file2_close(&fIJ);
 
466
  dpd_file2_close(&fAB);
 
467
 
 
468
  if (do_singles) {
 
469
 
 
470
    dpd_file2_mat_wrt(Sia);
 
471
    dpd_file2_mat_close(Sia);
 
472
 
 
473
    for(h=0; h < nirreps; h++)
 
474
      dpd_buf4_mat_irrep_close(Dijab_anti, h);
 
475
  }
 
476
 
 
477
  if (do_doubles) {
 
478
 
 
479
    dpd_file2_mat_close(Fme);
 
480
 
 
481
    for(h=0; h < nirreps; h++)
 
482
      dpd_buf4_mat_irrep_close(Wmnie, h);
 
483
 
 
484
    for(h=0; h < nirreps; h++) {
 
485
      dpd_buf4_mat_irrep_wrt(Sijab, h);
 
486
      dpd_buf4_mat_irrep_close(Sijab, h);
 
487
    }
 
488
  }
 
489
}
 
490
 
 
491
void cc3_sigma_UHF_AAB(dpdbuf4 *C2AA, dpdbuf4 *C2AB, dpdbuf4 *C2BA,
 
492
    dpdbuf4 *FAA, dpdbuf4 *FAB, dpdbuf4 *FBA,
 
493
    dpdbuf4 *EAA, dpdbuf4 *EAB, dpdbuf4 *EBA, 
 
494
    int do_singles, dpdbuf4 *DAA, dpdbuf4 *DAB, dpdfile2 *SIA, dpdfile2 *Sia,
 
495
    int do_doubles, dpdfile2 *FME, dpdfile2 *Fme,
 
496
    dpdbuf4 *WMAFE, dpdbuf4 *WMaFe, dpdbuf4 *WmAfE,
 
497
    dpdbuf4 *WMNIE, dpdbuf4 *WMnIe, dpdbuf4 *WmNiE,
 
498
    dpdbuf4 *SIJAB, dpdbuf4 *SIjAb, int *aoccpi, int *aocc_off, int *boccpi,
 
499
    int *bocc_off, int *avirtpi, int *avir_off, int *bvirtpi, int *bvir_off,
 
500
    double omega, FILE *outfile)
 
501
{
 
502
  int h, nirreps;
 
503
  int Gi, Gj, Gk, Gijk;
 
504
  int Ga, Gb, Gc, Gab;
 
505
  int Gij, ij, ji, Gjk, jk, Gbc, bc, Gcb, cb;
 
506
  int Gd, Gkd, kd, d, DD, ad, da, dc, Gid, id, Gac, ac, bd;
 
507
  int i, j, k, I, J, K;
 
508
  int a, b, c, A, B, C;
 
509
  int ab;
 
510
  double ***W1, ***W2, ***W3;
 
511
  dpdfile2 fIJ, fAB, fij, fab;
 
512
  int nrows, ncols, nlinks;
 
513
  int **W_offset, offset;
 
514
  double **Z;
 
515
  int Gm, m, Gmi, mi, im, mc, M;
 
516
  int Gmk, mk, ma, kj, Gim;
 
517
  int GS, GH, GU, GC, GX3;
 
518
 
 
519
  GC = C2AA->file.my_irrep;
 
520
  GU = EAA->file.my_irrep;
 
521
 
 
522
  if (do_singles)
 
523
    GH = DAA->file.my_irrep;
 
524
  else if (do_doubles)
 
525
    GH = WMAFE->file.my_irrep;
 
526
 
 
527
  GS = SIJAB->file.my_irrep;
 
528
  GX3 = GU^GC;
 
529
 
 
530
  nirreps = C2AA->params->nirreps;
 
531
 
 
532
  W_offset = init_int_matrix(nirreps, nirreps);
 
533
  for(Gab=0; Gab < nirreps; Gab++) {
 
534
    for(Ga=0,offset=0; Ga < nirreps; Ga++) {
 
535
      Gb = Ga ^ Gab;
 
536
      W_offset[Gab][Ga] = offset;
 
537
      offset += avirtpi[Ga] * avirtpi[Gb];
 
538
    }
 
539
  }
 
540
 
 
541
  if (do_singles) {
 
542
    dpd_file2_mat_init(SIA);
 
543
    dpd_file2_mat_rd(SIA);
 
544
    dpd_file2_mat_init(Sia);
 
545
    dpd_file2_mat_rd(Sia);
 
546
 
 
547
    for(h=0; h < nirreps; h++) {
 
548
      dpd_buf4_mat_irrep_init(DAA, h);
 
549
      dpd_buf4_mat_irrep_rd(DAA, h);
 
550
      dpd_buf4_mat_irrep_init(DAB, h);
 
551
      dpd_buf4_mat_irrep_rd(DAB, h);
 
552
    }
 
553
  }
 
554
 
 
555
  if (do_doubles) {
 
556
    dpd_file2_mat_init(FME);
 
557
    dpd_file2_mat_rd(FME);
 
558
    dpd_file2_mat_init(Fme);
 
559
    dpd_file2_mat_rd(Fme);
 
560
 
 
561
    for(h=0; h < nirreps; h++) {
 
562
      dpd_buf4_mat_irrep_init(WMnIe, h);
 
563
      dpd_buf4_mat_irrep_rd(WMnIe, h);
 
564
 
 
565
      dpd_buf4_mat_irrep_init(WMNIE, h);
 
566
      dpd_buf4_mat_irrep_rd(WMNIE, h);
 
567
 
 
568
      dpd_buf4_mat_irrep_init(WmNiE, h);
 
569
      dpd_buf4_mat_irrep_rd(WmNiE, h);
 
570
    }
 
571
    for(h=0; h < nirreps; h++) {
 
572
      dpd_buf4_mat_irrep_init(SIJAB, h);
 
573
      dpd_buf4_mat_irrep_rd(SIJAB, h);
 
574
    }
 
575
    for(h=0; h < nirreps; h++) {
 
576
      dpd_buf4_mat_irrep_init(SIjAb, h);
 
577
      dpd_buf4_mat_irrep_rd(SIjAb, h);
 
578
    }
 
579
  }
 
580
 
 
581
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
582
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
583
  dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
 
584
  dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
 
585
 
 
586
  /* target T3 amplitudes go in here */
 
587
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
588
  W2 = (double ***) malloc(nirreps * sizeof(double **));
 
589
  W3 = (double ***) malloc(nirreps * sizeof(double **));
 
590
 
 
591
  for(Gi=0; Gi < nirreps; Gi++) {
 
592
    for(Gj=0; Gj < nirreps; Gj++) {
 
593
      Gij = Gi ^ Gj;
 
594
      for(Gk=0; Gk < nirreps; Gk++) {
 
595
        Gijk = Gi ^ Gj ^ Gk;
 
596
        Gjk = Gj ^ Gk;
 
597
 
 
598
        for(Gab=0; Gab < nirreps; Gab++) {
 
599
          Gc = Gab ^ Gijk ^ GX3;
 
600
          W1[Gab] = dpd_block_matrix(FAA->params->coltot[Gab], bvirtpi[Gc]);
 
601
        }
 
602
        for(Ga=0; Ga < nirreps; Ga++) {
 
603
          Gcb = Ga ^ Gijk ^ GX3;
 
604
          W2[Ga] = dpd_block_matrix(avirtpi[Ga], WmAfE->params->coltot[Gcb]);  /* alpha-beta-alpha */
 
605
        }
 
606
        for(Gb=0; Gb < nirreps; Gb++) {
 
607
          Gac = Gb ^ Gijk ^ GX3;
 
608
          W3[Gb] = dpd_block_matrix(avirtpi[Gb], WMaFe->params->coltot[Gac]);  /* alpha-alpha-beta */
 
609
        }
 
610
 
 
611
        for(i=0; i < aoccpi[Gi]; i++) {
 
612
          I = aocc_off[Gi] + i;
 
613
          for(j=0; j < aoccpi[Gj]; j++) {
 
614
            J = aocc_off[Gj] + j;
 
615
            for(k=0; k < boccpi[Gk]; k++) {
 
616
              K = bocc_off[Gk] + k;
 
617
 
 
618
              T3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, C2AA, C2AB, C2BA, 
 
619
                     FAA, FAB, FBA, EAA, EAB, EBA, &fIJ, &fij, &fAB, &fab,
 
620
                     aoccpi, aocc_off, boccpi, bocc_off, avirtpi, avir_off, bvirtpi, bvir_off, omega);
 
621
 
 
622
          if (do_singles) {
 
623
 
 
624
              /* S_kc <-- 1/4 t_IJkABc <IJ||AB> */
 
625
 
 
626
              Gc = Gk ^ GS;  
 
627
              Gab = Gij ^ GH;
 
628
 
 
629
              ij = DAA->params->rowidx[I][J];
 
630
 
 
631
              nrows = DAA->params->coltot[Gij^GH];
 
632
              ncols = bvirtpi[Gc];
 
633
 
 
634
              if(nrows && ncols)
 
635
                C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, DAA->matrix[Gij][ij], 1,
 
636
                        1.0, Sia->matrix[Gk][k], 1);
 
637
 
 
638
              /* S_IA <-- t_IJkABc <Jk|Bc> */
 
639
 
 
640
              Ga = Gi ^ GS;   
 
641
              Gbc = Gjk ^ GH;
 
642
 
 
643
              jk = DAB->params->rowidx[J][K];
 
644
 
 
645
              for(Gab=0; Gab < nirreps; Gab++) {
 
646
                Gb = Ga ^ Gab;
 
647
                Gc = Gb ^ Gbc;
 
648
 
 
649
                ab = W_offset[Gab][Ga];
 
650
                bc = DAB->col_offset[Gjk][Gb];
 
651
 
 
652
                nrows = avirtpi[Ga];
 
653
                ncols = avirtpi[Gb] * bvirtpi[Gc];
 
654
 
 
655
                if(nrows && ncols)
 
656
                  C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][ab], ncols, &(DAB->matrix[Gjk][jk][bc]), 1,
 
657
                          1.0, SIA->matrix[Gi][i], 1);
 
658
              }
 
659
          }
 
660
 
 
661
          if (do_doubles) {
 
662
              /* S_IJAB <-- t_IJkABc F_kc */
 
663
              Gc = Gk ^ GH;   
 
664
              Gab = Gij ^ GS;
 
665
 
 
666
              nrows = SIJAB->params->coltot[Gij^GS];
 
667
              ncols = bvirtpi[Gc];
 
668
              ij = SIJAB->params->rowidx[I][J];
 
669
 
 
670
              if(nrows && ncols)
 
671
                C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][0], ncols, Fme->matrix[Gk][k], 1,
 
672
                        1.0, SIJAB->matrix[Gij][ij], 1);
 
673
 
 
674
              /* S_JkBc <-- t_IJkABc F_IA */
 
675
              Ga = Gi ^ GH;   
 
676
              Gbc = Gjk ^ GS;
 
677
 
 
678
              jk = C2AB->params->rowidx[J][K];
 
679
 
 
680
              for(Gab=0; Gab < nirreps; Gab++) {
 
681
                Gb = Ga ^ Gab;
 
682
                Gc = Gb ^ Gbc;
 
683
 
 
684
                ab = W_offset[Gab][Ga];
 
685
                bc = SIjAb->col_offset[Gjk][Gb];
 
686
 
 
687
                nrows = avirtpi[Ga];
 
688
                ncols = avirtpi[Gb] * bvirtpi[Gc];
 
689
 
 
690
                if(nrows && ncols)
 
691
                  C_DGEMV('t', nrows, ncols, 1.0, W1[Gab][ab], ncols, FME->matrix[Gi][i], 1,
 
692
                          1.0, &(SIjAb->matrix[Gjk][jk][bc]), 1);
 
693
              }
 
694
 
 
695
              /* S_JIDA <-- t_IJkABc W_kDcB */
 
696
              /* sort W1(AB,c) to W2(A,cB) */
 
697
              for(Gab=0; Gab < nirreps; Gab++) {
 
698
                Gc = Gab ^ Gijk ^ GX3;
 
699
                for(ab=0; ab < FAA->params->coltot[Gab]; ab++) {
 
700
                  A = FAA->params->colorb[Gab][ab][0];
 
701
                  B = FAA->params->colorb[Gab][ab][1];
 
702
                  Ga = FAA->params->rsym[A];
 
703
                  a = A - avir_off[Ga];
 
704
                  for(c=0; c < bvirtpi[Gc]; c++) {
 
705
                    C = bvir_off[Gc] + c;
 
706
                    cb = WmAfE->params->colidx[C][B];
 
707
                    W2[Ga][a][cb] = W1[Gab][ab][c];
 
708
                  }
 
709
                }
 
710
              }
 
711
 
 
712
              ji = SIJAB->params->rowidx[J][I];
 
713
 
 
714
              for(Gd=0; Gd < nirreps; Gd++) {
 
715
                Gkd = Gk ^ Gd; 
 
716
                Gcb = Gkd ^ GH; 
 
717
                Ga = Gd ^ Gij ^ GS;      
 
718
 
 
719
                kd = WmAfE->row_offset[Gkd][K];
 
720
                WmAfE->matrix[Gkd] = dpd_block_matrix(avirtpi[Gd], WmAfE->params->coltot[Gkd^GH]);
 
721
                dpd_buf4_mat_irrep_rd_block(WmAfE, Gkd, kd, avirtpi[Gd]);
 
722
                Z = block_matrix(avirtpi[Ga], avirtpi[Gd]);
 
723
 
 
724
                nrows = avirtpi[Ga];
 
725
                ncols = avirtpi[Gd];
 
726
                nlinks = WmAfE->params->coltot[Gkd^GH];
 
727
 
 
728
                if(nrows && ncols && nlinks)
 
729
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nlinks, 
 
730
                          WmAfE->matrix[Gkd][0], nlinks, 0.0, Z[0], ncols);
 
731
 
 
732
                for(a=0; a < avirtpi[Ga]; a++) {
 
733
                  A = avir_off[Ga] + a;
 
734
                  for(d=0; d < avirtpi[Gd]; d++) {
 
735
                    DD = avir_off[Gd] + d;
 
736
                    ad = SIJAB->params->colidx[A][DD];
 
737
                    da = SIJAB->params->colidx[DD][A];
 
738
                    SIJAB->matrix[Gij][ji][ad] += -Z[a][d];
 
739
                    SIJAB->matrix[Gij][ji][da] += Z[a][d];
 
740
                  }
 
741
                }
 
742
 
 
743
                dpd_free_block(WmAfE->matrix[Gkd], avirtpi[Gd], WmAfE->params->coltot[Gkd^GH]);
 
744
                free_block(Z);
 
745
              }
 
746
 
 
747
              /* S_JkDc <-- 1/2 t_IJkABc W_IDAB */
 
748
 
 
749
              jk = SIjAb->params->rowidx[J][K];
 
750
 
 
751
              for(Gd=0; Gd < nirreps; Gd++) {
 
752
                Gid = Gi ^ Gd; 
 
753
                Gab = Gid ^ GH; 
 
754
                Gc = Gab ^ Gijk ^ GX3;    
 
755
 
 
756
                id = WMAFE->row_offset[Gid][I];
 
757
                WMAFE->matrix[Gid] = dpd_block_matrix(avirtpi[Gd], WMAFE->params->coltot[Gid^GH]);
 
758
                dpd_buf4_mat_irrep_rd_block(WMAFE, Gid, id, avirtpi[Gd]);
 
759
                Z = block_matrix(bvirtpi[Gc], avirtpi[Gd]);
 
760
 
 
761
                nrows = bvirtpi[Gc];
 
762
                ncols = avirtpi[Gd];
 
763
                nlinks = WMAFE->params->coltot[Gid^GH];
 
764
 
 
765
                if(nrows && ncols && nlinks)
 
766
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows,
 
767
                          WMAFE->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
768
 
 
769
                for(c=0; c < bvirtpi[Gc]; c++) {
 
770
                  C = bvir_off[Gc] + c;
 
771
                  for(d=0; d < avirtpi[Gd]; d++) {
 
772
                    DD = avir_off[Gd] + d;
 
773
                    dc = SIjAb->params->colidx[DD][C];
 
774
                    SIjAb->matrix[Gjk][jk][dc] += Z[c][d];
 
775
                  }
 
776
                }
 
777
 
 
778
                free_block(Z);
 
779
                dpd_free_block(WMAFE->matrix[Gid], avirtpi[Gd], WMAFE->params->coltot[Gid^GH]);
 
780
              }
 
781
 
 
782
              /* t_JkBd <-- t_IJkABc W_IdAc */
 
783
              /* sort W1(AB,c) to W3(B,Ac) */
 
784
              for(Gab=0; Gab < nirreps; Gab++) {
 
785
                Gc = Gab ^ Gijk ^ GX3;
 
786
                for(ab=0; ab < FAA->params->coltot[Gab]; ab++) {
 
787
                  A = FAA->params->colorb[Gab][ab][0];
 
788
                  B = FAA->params->colorb[Gab][ab][1];
 
789
                  Gb = FAA->params->ssym[B];
 
790
                  b = B - avir_off[Gb];
 
791
                  for(c=0; c < bvirtpi[Gc]; c++) {
 
792
                    C = bvir_off[Gc] + c;
 
793
                    ac = WMaFe->params->colidx[A][C];
 
794
                    W3[Gb][b][ac] = W1[Gab][ab][c];
 
795
                  }
 
796
                }
 
797
              }
 
798
 
 
799
              jk = SIjAb->params->rowidx[J][K];
 
800
 
 
801
              for(Gd=0; Gd < nirreps; Gd++) {
 
802
                Gid = Gi ^ Gd; 
 
803
                Gac = Gid ^ GH; 
 
804
                Gb = Gac ^ Gijk ^ GX3;    
 
805
 
 
806
                id = WMaFe->row_offset[Gid][I]; 
 
807
                WMaFe->matrix[Gid] = dpd_block_matrix(bvirtpi[Gd], WMaFe->params->coltot[Gid^GH]);
 
808
                dpd_buf4_mat_irrep_rd_block(WMaFe, Gid, id, bvirtpi[Gd]);
 
809
                Z = block_matrix(avirtpi[Gb], bvirtpi[Gd]);
 
810
 
 
811
                nrows = avirtpi[Gb];
 
812
                ncols = bvirtpi[Gd];
 
813
                nlinks = WMaFe->params->coltot[Gid^GH];
 
814
 
 
815
                if(nrows && ncols && nlinks)
 
816
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3[Gb][0], nlinks,
 
817
                          WMaFe->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
818
 
 
819
                for(b=0; b < avirtpi[Gb]; b++) {
 
820
                  B = avir_off[Gb] + b;
 
821
                  for(d=0; d < bvirtpi[Gd]; d++) {
 
822
                    DD = bvir_off[Gd] + d;
 
823
                    bd = SIjAb->params->colidx[B][DD];
 
824
                    SIjAb->matrix[Gjk][jk][bd] += Z[b][d];
 
825
                  }
 
826
                }
 
827
 
 
828
                dpd_free_block(WMaFe->matrix[Gid], bvirtpi[Gd], WMaFe->params->coltot[Gid^GH]);
 
829
                free_block(Z);
 
830
              }
 
831
 
 
832
              /* S_MIAB <--- +t_IJkABc W_JkMc */
 
833
              /* S_IMAB <--- -t_IJkABc W_JkMc */
 
834
 
 
835
              jk = WMnIe->params->rowidx[J][K];
 
836
 
 
837
              for(Gm=0; Gm < nirreps; Gm++) {
 
838
                Gmi = Gm ^ Gi;  
 
839
                Gab = Gmi ^ GS;  
 
840
                Gc = Gab ^ Gijk ^ GX3;     
 
841
 
 
842
                mc = WMnIe->col_offset[Gjk][Gm];
 
843
 
 
844
                nrows = FAA->params->coltot[Gab];
 
845
                ncols = aoccpi[Gm];
 
846
                nlinks = bvirtpi[Gc];
 
847
 
 
848
                Z = dpd_block_matrix(nrows, ncols);
 
849
 
 
850
                if(nrows && ncols && nlinks)
 
851
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W1[Gab][0], nlinks,
 
852
                          &(WMnIe->matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
853
 
 
854
                for(m=0; m < ncols; m++) {
 
855
                  M = aocc_off[Gm] + m;
 
856
                  mi = SIJAB->params->rowidx[M][I];
 
857
                  im = SIJAB->params->rowidx[I][M];
 
858
                  for(ab=0; ab < nrows; ab++) {
 
859
                    SIJAB->matrix[Gmi][mi][ab] += Z[ab][m];
 
860
                    SIJAB->matrix[Gmi][im][ab] -= Z[ab][m];
 
861
                  }
 
862
                }
 
863
 
 
864
                dpd_free_block(Z, nrows, ncols);
 
865
              }
 
866
 
 
867
              /* t_MkBc <-- 1/2 t_IJkABc W_IJMA */
 
868
              /* sort W(AB,c) to W(A,Bc) */
 
869
              for(Gab=0; Gab < nirreps; Gab++) {
 
870
                Gc = Gab ^ Gijk ^ GX3;  
 
871
                for(ab=0; ab < FAA->params->coltot[Gab]; ab++ ){
 
872
                  A = FAA->params->colorb[Gab][ab][0];
 
873
                  B = FAA->params->colorb[Gab][ab][1];
 
874
                  Ga = FAA->params->rsym[A];
 
875
                  a = A - avir_off[Ga];
 
876
                  for(c=0; c < bvirtpi[Gc]; c++) {
 
877
                    C = bvir_off[Gc] + c;
 
878
                    bc = SIjAb->params->colidx[B][C];
 
879
                    W3[Ga][a][bc] = W1[Gab][ab][c];
 
880
                  }
 
881
                }
 
882
              }
 
883
 
 
884
              ij = WMNIE->params->rowidx[I][J];
 
885
 
 
886
              for(Gm=0; Gm < nirreps; Gm++) {
 
887
                Gmk = Gm ^ Gk; 
 
888
                Gbc = Gmk ^ GS; 
 
889
                Ga = Gbc ^ Gijk ^ GX3;    
 
890
 
 
891
                ma = WMNIE->col_offset[Gij][Gm];
 
892
 
 
893
                nrows = SIjAb->params->coltot[Gmk^GS];
 
894
                ncols = aoccpi[Gm];
 
895
                nlinks = avirtpi[Ga];
 
896
 
 
897
                Z = dpd_block_matrix(nrows, ncols);
 
898
 
 
899
                if(nrows && ncols && nlinks)
 
900
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W3[Ga][0], nrows,
 
901
                          &(WMNIE->matrix[Gij][ij][ma]), nlinks, 0.0, Z[0], ncols);
 
902
 
 
903
                for(m=0; m < aoccpi[Gm]; m++) {
 
904
                  M = aocc_off[Gm] + m;
 
905
                  mk = SIjAb->params->rowidx[M][K];
 
906
                  for(bc=0; bc < nrows; bc++) {
 
907
                    SIjAb->matrix[Gmk][mk][bc] += Z[bc][m];
 
908
                  }
 
909
                }
 
910
 
 
911
                dpd_free_block(Z, nrows, ncols);
 
912
              }
 
913
 
 
914
              /* S_ImBc <-- t_IJkABc W_kJmA */
 
915
              /* sort W(AB,c) to W(A,Bc) */
 
916
              for(Gab=0; Gab < nirreps; Gab++) {
 
917
                Gc = Gab ^ Gijk ^ GX3;  /* assumes totally symmetric */
 
918
                for(ab=0; ab < FAA->params->coltot[Gab]; ab++ ){
 
919
                  A = FAA->params->colorb[Gab][ab][0];
 
920
                  B = FAA->params->colorb[Gab][ab][1];
 
921
                  Ga = FAA->params->rsym[A];
 
922
                  a = A - avir_off[Ga];
 
923
                  for(c=0; c < bvirtpi[Gc]; c++) {
 
924
                    C = bvir_off[Gc] + c;
 
925
                    bc = SIjAb->params->colidx[B][C];
 
926
                    W3[Ga][a][bc] = W1[Gab][ab][c];
 
927
                  }
 
928
                }
 
929
              }
 
930
 
 
931
              kj = WmNiE->params->rowidx[K][J];
 
932
 
 
933
              for(Gm=0; Gm < nirreps; Gm++) {
 
934
                Gim = Gi ^ Gm; 
 
935
                Gbc = Gim ^ GS; 
 
936
                Ga = Gbc ^ Gijk ^ GX3;    
 
937
 
 
938
                ma = WmNiE->col_offset[Gjk][Gm];
 
939
 
 
940
                nrows = SIjAb->params->coltot[Gim^GS];
 
941
                ncols = boccpi[Gm];
 
942
                nlinks = avirtpi[Ga];
 
943
 
 
944
                Z = dpd_block_matrix(nrows, ncols);
 
945
 
 
946
                if(nrows && ncols && nlinks)
 
947
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 1.0, W3[Ga][0], nrows,
 
948
                          &(WmNiE->matrix[Gjk][kj][ma]), nlinks, 0.0, Z[0], ncols);
 
949
 
 
950
                for(m=0; m < boccpi[Gm]; m++) {
 
951
                  M = bocc_off[Gm] + m;
 
952
                  im = SIjAb->params->rowidx[I][M];
 
953
                  for(bc=0; bc < nrows; bc++) {
 
954
                    SIjAb->matrix[Gim][im][bc] += Z[bc][m];
 
955
                  }
 
956
                }
 
957
 
 
958
                dpd_free_block(Z, nrows, ncols);
 
959
              }
 
960
          }
 
961
 
 
962
            } /* k */
 
963
          } /* j */
 
964
        } /* i */
 
965
 
 
966
        for(Gab=0; Gab < nirreps; Gab++) {
 
967
          /* This will need to change for non-totally-symmetric cases */
 
968
          Gc = Gab ^ Gijk ^ GX3;
 
969
          dpd_free_block(W1[Gab], FAA->params->coltot[Gab], bvirtpi[Gc]);
 
970
        }
 
971
        for(Ga=0; Ga < nirreps; Ga++) {
 
972
          Gcb = Ga ^ Gijk ^ GX3; /* assumes totally symmetric */
 
973
          dpd_free_block(W2[Ga], avirtpi[Ga], WmAfE->params->coltot[Gcb]);
 
974
        }
 
975
        for(Gb=0; Gb < nirreps; Gb++) {
 
976
          Gac = Gb ^ Gijk ^ GX3; /* assumes totally symmetric */
 
977
          dpd_free_block(W3[Gb], avirtpi[Gb], WMaFe->params->coltot[Gac]);
 
978
        }
 
979
 
 
980
      } /* Gk */
 
981
    } /* Gj */
 
982
  } /* Gi */
 
983
 
 
984
  free(W1);
 
985
  free(W2);
 
986
  free(W3);
 
987
  free_int_matrix(W_offset);
 
988
 
 
989
  dpd_file2_close(&fIJ);
 
990
  dpd_file2_close(&fAB);
 
991
  dpd_file2_close(&fij);
 
992
  dpd_file2_close(&fab);
 
993
 
 
994
  if (do_singles) {
 
995
 
 
996
    dpd_file2_mat_wrt(SIA);
 
997
    dpd_file2_mat_close(SIA);
 
998
    dpd_file2_mat_wrt(Sia);
 
999
    dpd_file2_mat_close(Sia);
 
1000
 
 
1001
    for(h=0; h < nirreps; h++) {
 
1002
      dpd_buf4_mat_irrep_close(DAA, h);
 
1003
      dpd_buf4_mat_irrep_close(DAB, h);
 
1004
    }
 
1005
  }
 
1006
 
 
1007
  if (do_doubles) {
 
1008
 
 
1009
    dpd_file2_mat_close(FME);
 
1010
    dpd_file2_mat_close(Fme);
 
1011
 
 
1012
    for(h=0; h < nirreps; h++) {
 
1013
      dpd_buf4_mat_irrep_close(WMNIE, h);
 
1014
      dpd_buf4_mat_irrep_close(WMnIe, h);
 
1015
      dpd_buf4_mat_irrep_close(WmNiE, h);
 
1016
    }
 
1017
    for(h=0; h < nirreps; h++) {
 
1018
      dpd_buf4_mat_irrep_wrt(SIJAB, h);
 
1019
      dpd_buf4_mat_irrep_close(SIJAB, h);
 
1020
    }
 
1021
    for(h=0; h < nirreps; h++) {
 
1022
      dpd_buf4_mat_irrep_wrt(SIjAb, h);
 
1023
      dpd_buf4_mat_irrep_close(SIjAb, h);
 
1024
    }
 
1025
  }
 
1026
}
 
1027
 
 
1028
void cc3_sigma_UHF_BBA(dpdbuf4 *C2BB, dpdbuf4 *C2AB, dpdbuf4 *C2BA,
 
1029
    dpdbuf4 *FBB, dpdbuf4 *FAB, dpdbuf4 *FBA,  
 
1030
    dpdbuf4 *EBB, dpdbuf4 *EAB, dpdbuf4 *EBA,
 
1031
    int do_singles, dpdbuf4 *DBB, dpdbuf4 *DBA, dpdfile2 *SIA, dpdfile2 *Sia,
 
1032
    int do_doubles, dpdfile2 *FME, dpdfile2 *Fme,
 
1033
    dpdbuf4 *Wmafe, dpdbuf4 *WMaFe, dpdbuf4 *WmAfE,
 
1034
    dpdbuf4 *Wmnie, dpdbuf4 *WMnIe, dpdbuf4 *WmNiE,
 
1035
    dpdbuf4 *Sijab, dpdbuf4 *SIjAb, int *aoccpi, int *aocc_off, int *boccpi,
 
1036
    int *bocc_off, int *avirtpi, int *avir_off, int *bvirtpi, int *bvir_off,
 
1037
    double omega, FILE *outfile)
 
1038
{
 
1039
  int h, nirreps, S_irr;
 
1040
  int Gi, Gj, Gk, Gijk;
 
1041
  int Ga, Gb, Gc, Gab;
 
1042
  int Gij, ij, ji, Gjk, jk, Gbc, bc, Gcb, cb;
 
1043
  int Gd, Gkd, kd, d, DD, ad, da, dc, Gid, id, Gac, ac, bd;
 
1044
  int i, j, k, I, J, K;
 
1045
  int a, b, c, A, B, C;
 
1046
  int ab;
 
1047
  double ***W1, ***W2, ***W3;
 
1048
  dpdfile2 fIJ, fAB, fij, fab;
 
1049
  int nrows, ncols, nlinks;
 
1050
  int **W_offset, offset;
 
1051
  dpdbuf4 SiJaB;
 
1052
  double **Z;
 
1053
  int Gm, m, Gmi, mi, im, mc, M;
 
1054
  int Gmk, mk, ma, kj, Gim;
 
1055
  int GS, GH, GU, GC, GX3;
 
1056
    
 
1057
  GC = C2BB->file.my_irrep;
 
1058
  GU = EBB->file.my_irrep;
 
1059
    
 
1060
  if (do_singles)
 
1061
    GH = DBB->file.my_irrep;
 
1062
  else if (do_doubles)
 
1063
    GH = Wmafe->file.my_irrep;
 
1064
  
 
1065
  GS = Sijab->file.my_irrep;
 
1066
  GX3 = GU^GC;
 
1067
 
 
1068
  nirreps = C2BB->params->nirreps;
 
1069
 
 
1070
  W_offset = init_int_matrix(nirreps, nirreps);
 
1071
  for(Gab=0; Gab < nirreps; Gab++) {
 
1072
    for(Ga=0,offset=0; Ga < nirreps; Ga++) {
 
1073
      Gb = Ga ^ Gab;
 
1074
      W_offset[Gab][Ga] = offset;
 
1075
      offset += bvirtpi[Ga] * bvirtpi[Gb];
 
1076
    }
 
1077
  }
 
1078
 
 
1079
  if (do_singles) {
 
1080
 
 
1081
    dpd_file2_mat_init(SIA);
 
1082
    dpd_file2_mat_rd(SIA);
 
1083
    dpd_file2_mat_init(Sia);
 
1084
    dpd_file2_mat_rd(Sia);
 
1085
 
 
1086
    for(h=0; h < nirreps; h++) {
 
1087
      dpd_buf4_mat_irrep_init(DBB, h);
 
1088
      dpd_buf4_mat_irrep_rd(DBB, h);
 
1089
      dpd_buf4_mat_irrep_init(DBA, h);
 
1090
      dpd_buf4_mat_irrep_rd(DBA, h);
 
1091
    }
 
1092
  }
 
1093
 
 
1094
  if (do_doubles) {
 
1095
 
 
1096
    dpd_file2_mat_init(FME);
 
1097
    dpd_file2_mat_rd(FME);
 
1098
    dpd_file2_mat_init(Fme);
 
1099
    dpd_file2_mat_rd(Fme);
 
1100
 
 
1101
    for(h=0; h < nirreps; h++) {
 
1102
      dpd_buf4_mat_irrep_init(WmNiE, h);
 
1103
      dpd_buf4_mat_irrep_rd(WmNiE, h);
 
1104
  
 
1105
      dpd_buf4_mat_irrep_init(Wmnie, h);
 
1106
      dpd_buf4_mat_irrep_rd(Wmnie, h);
 
1107
  
 
1108
      dpd_buf4_mat_irrep_init(WMnIe, h);
 
1109
      dpd_buf4_mat_irrep_rd(WMnIe, h);
 
1110
    }
 
1111
    for(h=0; h < nirreps; h++) {
 
1112
      dpd_buf4_mat_irrep_init(Sijab, h);
 
1113
      dpd_buf4_mat_irrep_rd(Sijab, h);
 
1114
    }
 
1115
 
 
1116
    /* put result here until end of function */
 
1117
    S_irr = SIjAb->file.my_irrep;
 
1118
    dpd_buf4_init(&SiJaB, CC_TMP0, S_irr, 23, 29, 23, 29, 0, "CC3 SiJaB");
 
1119
    for(h=0; h < nirreps; h++) {
 
1120
      dpd_buf4_mat_irrep_init(&SiJaB, h);
 
1121
    }
 
1122
  }
 
1123
 
 
1124
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
1125
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
1126
  dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
 
1127
  dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
 
1128
 
 
1129
  /* target T3 amplitudes go in here */
 
1130
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
1131
  W2 = (double ***) malloc(nirreps * sizeof(double **));
 
1132
  W3 = (double ***) malloc(nirreps * sizeof(double **));
 
1133
 
 
1134
  for(Gi=0; Gi < nirreps; Gi++) {
 
1135
    for(Gj=0; Gj < nirreps; Gj++) {
 
1136
      Gij = Gi ^ Gj;
 
1137
      for(Gk=0; Gk < nirreps; Gk++) {
 
1138
        Gijk = Gi ^ Gj ^ Gk;
 
1139
        Gjk = Gj ^ Gk;
 
1140
 
 
1141
        for(Gab=0; Gab < nirreps; Gab++) {
 
1142
          Gc = Gab ^ Gijk ^ GX3;
 
1143
          W1[Gab] = dpd_block_matrix(FBB->params->coltot[Gab], avirtpi[Gc]);
 
1144
        }
 
1145
        for(Ga=0; Ga < nirreps; Ga++) {
 
1146
          Gcb = Ga ^ Gijk ^ GX3;
 
1147
          W2[Ga] = dpd_block_matrix(bvirtpi[Ga], WMaFe->params->coltot[Gcb]);
 
1148
        }
 
1149
        for(Gb=0; Gb < nirreps; Gb++) {
 
1150
          Gac = Gb ^ Gijk ^ GX3; 
 
1151
          W3[Gb] = dpd_block_matrix(bvirtpi[Gb], WmAfE->params->coltot[Gac]);
 
1152
        }
 
1153
 
 
1154
        for(i=0; i < boccpi[Gi]; i++) {
 
1155
          I = bocc_off[Gi] + i;
 
1156
          for(j=0; j < boccpi[Gj]; j++) {
 
1157
            J = bocc_off[Gj] + j;
 
1158
            for(k=0; k < aoccpi[Gk]; k++) {
 
1159
              K = aocc_off[Gk] + k;
 
1160
 
 
1161
              T3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, C2BB, C2BA, C2AB, 
 
1162
                     FBB, FBA, FAB, EBB, EBA, EAB, &fij, &fIJ, &fab, &fAB,
 
1163
                     boccpi, bocc_off, aoccpi, aocc_off, bvirtpi, bvir_off, avirtpi, avir_off, omega);
 
1164
 
 
1165
          if (do_singles) {
 
1166
              /* S_KC <-- 1/4 t_ijKabC <ij||ab> */
 
1167
 
 
1168
              Gc = Gk ^ GS; 
 
1169
              Gab = Gij ^ GH;
 
1170
 
 
1171
              ij = DBB->params->rowidx[I][J];
 
1172
 
 
1173
              nrows = DBB->params->coltot[Gij^GH];
 
1174
              ncols = avirtpi[Gc];
 
1175
 
 
1176
              if(nrows && ncols)
 
1177
                C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, DBB->matrix[Gij][ij], 1,
 
1178
                        1.0, SIA->matrix[Gk][k], 1);
 
1179
 
 
1180
              /* S_ia <-- t_ijKabC <jK|bC> */
 
1181
 
 
1182
              Ga = Gi ^ GS;  
 
1183
              Gbc = Gjk ^ GH;
 
1184
 
 
1185
              jk = DBA->params->rowidx[J][K];
 
1186
 
 
1187
              for(Gab=0; Gab < nirreps; Gab++) {
 
1188
                Gb = Ga ^ Gab;
 
1189
                Gc = Gb ^ Gbc;
 
1190
 
 
1191
                ab = W_offset[Gab][Ga];
 
1192
                bc = DBA->col_offset[Gjk][Gb];
 
1193
 
 
1194
                nrows = bvirtpi[Ga];
 
1195
                ncols = bvirtpi[Gb] * avirtpi[Gc];
 
1196
 
 
1197
                if(nrows && ncols)
 
1198
                  C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][ab], ncols, &(DBA->matrix[Gjk][jk][bc]), 1,
 
1199
                          1.0, Sia->matrix[Gi][i], 1);
 
1200
              }
 
1201
          } /* end do_singles */
 
1202
 
 
1203
          if (do_doubles ) {
 
1204
              /* S_ijab <-- t_ijKabC F_KC */
 
1205
              Gc = Gk ^ GH;  
 
1206
              Gab = Gij ^ GS;
 
1207
 
 
1208
              nrows = Sijab->params->coltot[Gij^GS];
 
1209
              ncols = avirtpi[Gc];
 
1210
              ij = Sijab->params->rowidx[I][J];
 
1211
 
 
1212
              if(nrows && ncols)
 
1213
                C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][0], ncols, FME->matrix[Gk][k], 1,
 
1214
                        1.0, Sijab->matrix[Gij][ij], 1);
 
1215
 
 
1216
              /* S_jKbC <-- t_ijKabC F_ia */
 
1217
              Ga = Gi ^ GH;
 
1218
              Gbc = Gjk ^ GS;
 
1219
 
 
1220
              jk = C2BA->params->rowidx[J][K];
 
1221
 
 
1222
              for(Gab=0; Gab < nirreps; Gab++) {
 
1223
                Gb = Ga ^ Gab;
 
1224
                Gc = Gb ^ Gbc;
 
1225
 
 
1226
                ab = W_offset[Gab][Ga];
 
1227
                bc = SiJaB.col_offset[Gjk][Gb];
 
1228
 
 
1229
                nrows = bvirtpi[Ga];
 
1230
                ncols = bvirtpi[Gb] * avirtpi[Gc];
 
1231
 
 
1232
                if(nrows && ncols)
 
1233
                  C_DGEMV('t', nrows, ncols, 1.0, W1[Gab][ab], ncols, Fme->matrix[Gi][i], 1,
 
1234
                          1.0, &(SiJaB.matrix[Gjk][jk][bc]), 1);
 
1235
              }
 
1236
 
 
1237
              /* S_jida <-- t_ijKabC W_KdCb */
 
1238
              /* sort W1(ab,C) to W2(a,Cb) */
 
1239
              for(Gab=0; Gab < nirreps; Gab++) {
 
1240
                Gc = Gab ^ Gijk ^ GX3;
 
1241
                for(ab=0; ab < FBB->params->coltot[Gab]; ab++) {
 
1242
                  A = FBB->params->colorb[Gab][ab][0];
 
1243
                  B = FBB->params->colorb[Gab][ab][1];
 
1244
                  Ga = FBB->params->rsym[A];
 
1245
                  a = A - bvir_off[Ga];
 
1246
                  for(c=0; c < avirtpi[Gc]; c++) {
 
1247
                    C = avir_off[Gc] + c;
 
1248
                    cb = WMaFe->params->colidx[C][B];
 
1249
                    W2[Ga][a][cb] = W1[Gab][ab][c];
 
1250
                  }
 
1251
                }
 
1252
              }
 
1253
 
 
1254
              ji = Sijab->params->rowidx[J][I];
 
1255
 
 
1256
              for(Gd=0; Gd < nirreps; Gd++) {
 
1257
                Gkd = Gk ^ Gd; 
 
1258
                Gcb = Gkd ^ GH; 
 
1259
                Ga = Gd ^ Gij ^ GS;      
 
1260
 
 
1261
                kd = WMaFe->row_offset[Gkd][K];
 
1262
                WMaFe->matrix[Gkd] = dpd_block_matrix(bvirtpi[Gd], WMaFe->params->coltot[Gkd^GH]);
 
1263
                dpd_buf4_mat_irrep_rd_block(WMaFe, Gkd, kd, bvirtpi[Gd]);
 
1264
                Z = block_matrix(bvirtpi[Ga], bvirtpi[Gd]);
 
1265
 
 
1266
                nrows = bvirtpi[Ga];
 
1267
                ncols = bvirtpi[Gd];
 
1268
                nlinks = WMaFe->params->coltot[Gkd^GH];
 
1269
 
 
1270
                if(nrows && ncols && nlinks)
 
1271
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nlinks, 
 
1272
                          WMaFe->matrix[Gkd][0], nlinks, 0.0, Z[0], ncols);
 
1273
 
 
1274
                for(a=0; a < bvirtpi[Ga]; a++) {
 
1275
                  A = bvir_off[Ga] + a;
 
1276
                  for(d=0; d < bvirtpi[Gd]; d++) {
 
1277
                    DD = bvir_off[Gd] + d;
 
1278
                    ad = Sijab->params->colidx[A][DD];
 
1279
                    da = Sijab->params->colidx[DD][A];
 
1280
                    Sijab->matrix[Gij][ji][ad] += -Z[a][d];
 
1281
                    Sijab->matrix[Gij][ji][da] += Z[a][d];
 
1282
                  }
 
1283
                }
 
1284
 
 
1285
                dpd_free_block(WMaFe->matrix[Gkd], bvirtpi[Gd], WMaFe->params->coltot[Gkd^GH]);
 
1286
                free_block(Z);
 
1287
              }
 
1288
 
 
1289
              /* t_jKcD <-- 1/2 t_ijKabC W_idab */
 
1290
 
 
1291
              jk = SiJaB.params->rowidx[J][K];
 
1292
 
 
1293
              for(Gd=0; Gd < nirreps; Gd++) {
 
1294
                Gid = Gi ^ Gd;
 
1295
                Gab = Gid ^ GH;
 
1296
                Gc = Gab ^ Gijk ^ GX3;   
 
1297
 
 
1298
                id = Wmafe->row_offset[Gid][I];
 
1299
                Wmafe->matrix[Gid] = dpd_block_matrix(bvirtpi[Gd], Wmafe->params->coltot[Gid^GH]);
 
1300
                dpd_buf4_mat_irrep_rd_block(Wmafe, Gid, id, bvirtpi[Gd]);
 
1301
                Z = block_matrix(avirtpi[Gc], bvirtpi[Gd]);
 
1302
 
 
1303
                nrows = avirtpi[Gc];
 
1304
                ncols = bvirtpi[Gd];
 
1305
                nlinks = Wmafe->params->coltot[Gid^GH];
 
1306
 
 
1307
                if(nrows && ncols && nlinks)
 
1308
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows,
 
1309
                          Wmafe->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
1310
 
 
1311
                for(c=0; c < avirtpi[Gc]; c++) {
 
1312
                  C = avir_off[Gc] + c;
 
1313
                  for(d=0; d < bvirtpi[Gd]; d++) {
 
1314
                    DD = bvir_off[Gd] + d;
 
1315
                    dc = SiJaB.params->colidx[DD][C];
 
1316
                    SiJaB.matrix[Gjk][jk][dc] += Z[c][d];
 
1317
                  }
 
1318
                }
 
1319
 
 
1320
                free_block(Z);
 
1321
                dpd_free_block(Wmafe->matrix[Gid], bvirtpi[Gd], Wmafe->params->coltot[Gid^GH]);
 
1322
              }
 
1323
 
 
1324
              /* t_jKbD <-- t_ijKabC W_iDaC */
 
1325
              /* sort W1(ab,C) to W3(b,aC) */
 
1326
              for(Gab=0; Gab < nirreps; Gab++) {
 
1327
                Gc = Gab ^ Gijk ^ GX3;
 
1328
                for(ab=0; ab < FBB->params->coltot[Gab]; ab++) {
 
1329
                  A = FBB->params->colorb[Gab][ab][0];
 
1330
                  B = FBB->params->colorb[Gab][ab][1];
 
1331
                  Gb = FBB->params->ssym[B];
 
1332
                  b = B - bvir_off[Gb];
 
1333
                  for(c=0; c < avirtpi[Gc]; c++) {
 
1334
                    C = avir_off[Gc] + c;
 
1335
                    ac = WmAfE->params->colidx[A][C];
 
1336
                    W3[Gb][b][ac] = W1[Gab][ab][c];
 
1337
                  }
 
1338
                }
 
1339
              }
 
1340
 
 
1341
              jk = SiJaB.params->rowidx[J][K];
 
1342
 
 
1343
              for(Gd=0; Gd < nirreps; Gd++) {
 
1344
                Gid = Gi ^ Gd;
 
1345
                Gac = Gid ^ GH;
 
1346
                Gb = Gac ^ Gijk ^ GX3;   
 
1347
 
 
1348
                id = WmAfE->row_offset[Gid][I]; 
 
1349
                WmAfE->matrix[Gid] = dpd_block_matrix(avirtpi[Gd], WmAfE->params->coltot[Gid^GH]);
 
1350
                dpd_buf4_mat_irrep_rd_block(WmAfE, Gid, id, avirtpi[Gd]);
 
1351
                Z = block_matrix(bvirtpi[Gb], avirtpi[Gd]);
 
1352
 
 
1353
                nrows = bvirtpi[Gb];
 
1354
                ncols = avirtpi[Gd];
 
1355
                nlinks = WmAfE->params->coltot[Gid^GH];
 
1356
 
 
1357
                if(nrows && ncols && nlinks)
 
1358
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3[Gb][0], nlinks,
 
1359
                          WmAfE->matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
1360
 
 
1361
                for(b=0; b < bvirtpi[Gb]; b++) {
 
1362
                  B = bvir_off[Gb] + b;
 
1363
                  for(d=0; d < avirtpi[Gd]; d++) {
 
1364
                    DD = avir_off[Gd] + d;
 
1365
                    bd = SiJaB.params->colidx[B][DD];
 
1366
                    SiJaB.matrix[Gjk][jk][bd] += Z[b][d];
 
1367
                  }
 
1368
                }
 
1369
 
 
1370
                dpd_free_block(WmAfE->matrix[Gid], avirtpi[Gd], WmAfE->params->coltot[Gid^GH]);
 
1371
                free_block(Z);
 
1372
              }
 
1373
 
 
1374
              /* S_miab <--- +t_ijKabC W_jKmC */
 
1375
              /* S_imab <--- -t_ijKabC W_jKmC */
 
1376
 
 
1377
              jk = WmNiE->params->rowidx[J][K];
 
1378
 
 
1379
              for(Gm=0; Gm < nirreps; Gm++) {
 
1380
                Gmi = Gm ^ Gi; 
 
1381
                Gab = Gmi ^ GS; 
 
1382
                Gc = Gab ^ Gijk ^ GX3;    
 
1383
 
 
1384
                mc = WmNiE->col_offset[Gjk][Gm];
 
1385
 
 
1386
                nrows = Sijab->params->coltot[Gab];
 
1387
                ncols = boccpi[Gm];
 
1388
                nlinks = avirtpi[Gc];
 
1389
 
 
1390
                Z = dpd_block_matrix(nrows, ncols);
 
1391
 
 
1392
                if(nrows && ncols && nlinks)
 
1393
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W1[Gab][0], nlinks,
 
1394
                          &(WmNiE->matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
1395
 
 
1396
                for(m=0; m < ncols; m++) {
 
1397
                  M = bocc_off[Gm] + m;
 
1398
                  mi = Sijab->params->rowidx[M][I];
 
1399
                  im = Sijab->params->rowidx[I][M];
 
1400
                  for(ab=0; ab < nrows; ab++) {
 
1401
                    Sijab->matrix[Gmi][mi][ab] += Z[ab][m];
 
1402
                    Sijab->matrix[Gmi][im][ab] -= Z[ab][m];
 
1403
                  }
 
1404
                }
 
1405
 
 
1406
                dpd_free_block(Z, nrows, ncols);
 
1407
              }
 
1408
 
 
1409
              /* t_mKbC <-- 1/2 t_ijKabC W_ijma */
 
1410
              /* sort W(ab,C) to W(a,bC) */
 
1411
              for(Gab=0; Gab < nirreps; Gab++) {
 
1412
                Gc = Gab ^ Gijk ^ GX3; 
 
1413
                for(ab=0; ab < FBB->params->coltot[Gab]; ab++ ){
 
1414
                  A = FBB->params->colorb[Gab][ab][0];
 
1415
                  B = FBB->params->colorb[Gab][ab][1];
 
1416
                  Ga = FBB->params->rsym[A];
 
1417
                  a = A - bvir_off[Ga];
 
1418
                  for(c=0; c < avirtpi[Gc]; c++) {
 
1419
                    C = avir_off[Gc] + c;
 
1420
                    bc = SiJaB.params->colidx[B][C];
 
1421
                    W3[Ga][a][bc] = W1[Gab][ab][c];
 
1422
                  }
 
1423
                }
 
1424
              }
 
1425
 
 
1426
              ij = Wmnie->params->rowidx[I][J];
 
1427
 
 
1428
              for(Gm=0; Gm < nirreps; Gm++) {
 
1429
                Gmk = Gm ^ Gk;
 
1430
                Gbc = Gmk ^ GS;
 
1431
                Ga = Gbc ^ Gijk ^ GX3;   
 
1432
 
 
1433
                ma = Wmnie->col_offset[Gij][Gm];
 
1434
 
 
1435
                nrows = SiJaB.params->coltot[Gmk^GS];
 
1436
                ncols = boccpi[Gm];
 
1437
                nlinks = bvirtpi[Ga];
 
1438
 
 
1439
                Z = dpd_block_matrix(nrows, ncols);
 
1440
 
 
1441
                if(nrows && ncols && nlinks)
 
1442
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W3[Ga][0], nrows,
 
1443
                          &(Wmnie->matrix[Gij][ij][ma]), nlinks, 0.0, Z[0], ncols);
 
1444
 
 
1445
                for(m=0; m < boccpi[Gm]; m++) {
 
1446
                  M = bocc_off[Gm] + m;
 
1447
                  mk = SiJaB.params->rowidx[M][K];
 
1448
                  for(bc=0; bc < nrows; bc++) {
 
1449
                    SiJaB.matrix[Gmk][mk][bc] += Z[bc][m];
 
1450
                  }
 
1451
                }
 
1452
 
 
1453
                dpd_free_block(Z, nrows, ncols);
 
1454
              }
 
1455
 
 
1456
              /* S_iMbC <-- t_ijKabC W_KjMa */
 
1457
              /* sort W(ab,C) to W(a,bC) */
 
1458
              for(Gab=0; Gab < nirreps; Gab++) {
 
1459
                Gc = Gab ^ Gijk ^ GX3; 
 
1460
                for(ab=0; ab < FBB->params->coltot[Gab]; ab++ ){
 
1461
                  A = FBB->params->colorb[Gab][ab][0];
 
1462
                  B = FBB->params->colorb[Gab][ab][1];
 
1463
                  Ga = FBB->params->rsym[A];
 
1464
                  a = A - bvir_off[Ga];
 
1465
                  for(c=0; c < avirtpi[Gc]; c++) {
 
1466
                    C = avir_off[Gc] + c;
 
1467
                    bc = SiJaB.params->colidx[B][C];
 
1468
                    W3[Ga][a][bc] = W1[Gab][ab][c];
 
1469
                  }
 
1470
                }
 
1471
              }
 
1472
 
 
1473
              kj = WMnIe->params->rowidx[K][J];
 
1474
 
 
1475
              for(Gm=0; Gm < nirreps; Gm++) {
 
1476
                Gim = Gi ^ Gm; 
 
1477
                Gbc = Gim ^ GS; 
 
1478
                Ga = Gbc ^ Gijk ^ GX3;    
 
1479
 
 
1480
                ma = WMnIe->col_offset[Gjk][Gm];
 
1481
 
 
1482
                nrows = SiJaB.params->coltot[Gim^GS];
 
1483
                ncols = aoccpi[Gm];
 
1484
                nlinks = bvirtpi[Ga];
 
1485
 
 
1486
                Z = dpd_block_matrix(nrows, ncols);
 
1487
 
 
1488
                if(nrows && ncols && nlinks)
 
1489
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 1.0, W3[Ga][0], nrows,
 
1490
                          &(WMnIe->matrix[Gjk][kj][ma]), nlinks, 0.0, Z[0], ncols);
 
1491
 
 
1492
                for(m=0; m < aoccpi[Gm]; m++) {
 
1493
                  M = aocc_off[Gm] + m;
 
1494
                  im = SiJaB.params->rowidx[I][M];
 
1495
                  for(bc=0; bc < nrows; bc++) {
 
1496
                    SiJaB.matrix[Gim][im][bc] += Z[bc][m];
 
1497
                  }
 
1498
                }
 
1499
 
 
1500
                dpd_free_block(Z, nrows, ncols);
 
1501
              }
 
1502
          } /* end do_doubles */
 
1503
 
 
1504
            } /* k */
 
1505
          } /* j */
 
1506
        } /* i */
 
1507
 
 
1508
        for(Gab=0; Gab < nirreps; Gab++) {
 
1509
          Gc = Gab ^ Gijk ^ GX3;
 
1510
          dpd_free_block(W1[Gab], FBB->params->coltot[Gab], avirtpi[Gc]);
 
1511
        }
 
1512
        for(Ga=0; Ga < nirreps; Ga++) {
 
1513
          Gcb = Ga ^ Gijk ^ GX3; 
 
1514
          dpd_free_block(W2[Ga], bvirtpi[Ga], WMaFe->params->coltot[Gcb]);
 
1515
        }
 
1516
        for(Gb=0; Gb < nirreps; Gb++) {
 
1517
          Gac = Gb ^ Gijk ^ GX3;
 
1518
          dpd_free_block(W3[Gb], bvirtpi[Gb], WmAfE->params->coltot[Gac]);
 
1519
        }
 
1520
      } /* Gk */
 
1521
    } /* Gj */
 
1522
  } /* Gi */
 
1523
 
 
1524
  free(W1);
 
1525
  free(W2);
 
1526
  free(W3);
 
1527
  free_int_matrix(W_offset);
 
1528
 
 
1529
  dpd_file2_close(&fIJ);
 
1530
  dpd_file2_close(&fAB);
 
1531
  dpd_file2_close(&fij);
 
1532
  dpd_file2_close(&fab);
 
1533
 
 
1534
  if (do_singles) {
 
1535
 
 
1536
    dpd_file2_mat_wrt(SIA);
 
1537
    dpd_file2_mat_close(SIA);
 
1538
    dpd_file2_mat_wrt(Sia);
 
1539
    dpd_file2_mat_close(Sia);
 
1540
 
 
1541
    for(h=0; h < nirreps; h++) {
 
1542
      dpd_buf4_mat_irrep_close(DBB, h);
 
1543
      dpd_buf4_mat_irrep_close(DBA, h);
 
1544
    }
 
1545
  }
 
1546
 
 
1547
  if (do_doubles) {
 
1548
 
 
1549
    dpd_file2_mat_close(FME);
 
1550
    dpd_file2_mat_close(Fme);
 
1551
 
 
1552
    for(h=0; h < nirreps; h++) {
 
1553
      dpd_buf4_mat_irrep_close(WmNiE, h);
 
1554
      dpd_buf4_mat_irrep_close(Wmnie, h);
 
1555
      dpd_buf4_mat_irrep_close(WMnIe, h);
 
1556
    }
 
1557
 
 
1558
    for(h=0; h < nirreps; h++) {
 
1559
      dpd_buf4_mat_irrep_wrt(Sijab, h);
 
1560
      dpd_buf4_mat_irrep_close(Sijab, h);
 
1561
    }
 
1562
    for(h=0; h < nirreps; h++) {
 
1563
      dpd_buf4_mat_irrep_wrt(&SiJaB, h);
 
1564
      dpd_buf4_mat_irrep_close(&SiJaB, h);
 
1565
    }
 
1566
    dpd_buf4_sort(&SiJaB, CC_TMP0, qpsr, 22, 28, "CC3 SIjAb");
 
1567
    dpd_buf4_close(&SiJaB);
 
1568
    dpd_buf4_init(&SiJaB, CC_TMP0, S_irr, 22, 28, 22, 28, 0, "CC3 SIjAb");
 
1569
    dpd_buf4_axpy(&SiJaB, SIjAb, 1.0);
 
1570
    dpd_buf4_close(&SiJaB);
 
1571
  }
 
1572
 
 
1573
}
 
1574
 
 
1575
} /* extern "C" */