← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/lib/libdpd/cc3_sigma_UHF.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

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