← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cctriples/T3_UHF_AAA.cc

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

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

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cctriples/T3_UHF_AAA.cc
 
1
/*! \file
 
2
    \ingroup DPD
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
/* T3_UHF_AAA(): Computes all connected and disconnected T3(IJK,ABC)
 
7
** amplitudes for a given I, J, K combination for input C2, F, and E
 
8
** intermediates.  This function will work for AAA or BBB spin cases,
 
9
** with either RHF/ROHF or UHF orbitals.
 
10
**
 
11
** Arguments: 
 
12
**
 
13
**   double ***W: The target connected triples amplitudes in an
 
14
**   nirreps x AB x C array.  The memory for this must be allocated
 
15
**   externally.
 
16
**
 
17
**   double ***V: The target disconnected triples amplitudes (if
 
18
**   requested) in an nirreps x AB x C array.  The memory for this
 
19
**   must be allocated externally.
 
20
**
 
21
**   int disc: Boolean: 1 == computed disconnected triples; 0 == don't
 
22
**
 
23
**   int nirreps: Number of irreps.
 
24
**
 
25
**   int I: Absolute index of orbital I.
 
26
**
 
27
**   int Gi: Irrep of I.
 
28
**
 
29
**   int J: Absolute index of orbital J.
 
30
**
 
31
**   int Gj: Irrep of J.
 
32
**
 
33
**   int K: Absolute index of orbital K.
 
34
**
 
35
**   int Gk: Irrep of K.
 
36
**
 
37
**   dpdbuf4 *C2: Pointer to dpd buffer for double excitation amps,
 
38
**   ordered (IJ,AB).
 
39
**
 
40
**   dpdbuf4 *F: Pointer to dpd buffer for three-virtual-index
 
41
**   intermediate, ordered (IA,BC).
 
42
**
 
43
**   dpdbuf4 *E: Pointer to dpd buffer for three-occupied-index
 
44
**   intermediate, ordered (IJ,KA).
 
45
**
 
46
**   dpdfile *C1: If disconnected T3's are requested, pointer to dpd
 
47
**   buffer for single-excitation amps.
 
48
**
 
49
**   dpdbuf4 *D: If disconnected T3's are requested, pointer to dpd
 
50
**   buffer for <IJ||ab> integrals.
 
51
**
 
52
**   dpdfile2 *fIA: Pointer to the dpd file2 for the occ-vir block of
 
53
**   the Fock matrix (or other appropriate one-electron operator).
 
54
**
 
55
**   dpdfile2 *fIJ: Pointer to the dpd file2 for the occ-occ block of
 
56
**   the Fock matrix (or other appropriate one-electron operator).
 
57
**   
 
58
**   dpdfile2 *fAB: Pointer to the dpd file2 for the vir-vir block of
 
59
**   the Fock matrix (or other appropriate one-electron operator).
 
60
**   
 
61
**   int *occpi: Number of occupied orbitals per irrep lookup array.
 
62
**
 
63
**   int *occ_off: Offset lookup for translating between absolute and
 
64
**   relative orbital indices for occupied space.
 
65
**
 
66
**   int *virtpi: Number of virtual orbitals per irrep lookup array.
 
67
**
 
68
**   int *vir_off: Offset lookup for translating between absolute and
 
69
**   relative orbital indices for virtual space.
 
70
**
 
71
**   double omega: a constant to add to the final denominators -
 
72
**   needed for CC3 EOM
 
73
**
 
74
** TDC, July 2004
 
75
** Modified to return disconnected triples, TDC, Feburary 2008
 
76
*/
 
77
 
 
78
#include <stdio.h>
 
79
#include <stdlib.h>
 
80
#include <string.h>
 
81
#include <math.h>
 
82
#include <libqt/qt.h>
 
83
#include <libdpd/dpd.h>
 
84
#include <ccfiles.h>
 
85
 
 
86
namespace psi { namespace cctriples {
 
87
 
 
88
    void T3_UHF_AAA(double ***W, double ***V, int disc, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk, 
 
89
                    dpdbuf4 *C2, dpdbuf4 *F, dpdbuf4 *E, dpdfile2 *C1, dpdbuf4 *D, dpdfile2 *fIA, dpdfile2 *fIJ, dpdfile2 *fAB,
 
90
                    int *occpi, int *occ_off, int *virtpi, int *vir_off, double omega)
 
91
    {
 
92
      int h;
 
93
      int i, j, k;
 
94
      int ij, ji, ik, ki, jk, kj;
 
95
      int Gij, Gji, Gik, Gki, Gjk, Gkj, Gijk;
 
96
      int Ga, Gb, Gc;
 
97
      int Gd, Gl;
 
98
      int Gid, Gjd, Gkd;
 
99
      int Gab, Gcb, Gca, Gbc, Gac;
 
100
      int Gla, Glb, Glc;
 
101
      int Gil, Gjl, Gkl;
 
102
      int a, b, c, A, B, C;
 
103
      int ab, bc, ac;
 
104
      int cd, bd, ad;
 
105
      int id, jd, kd;
 
106
      int la, lb, lc;
 
107
      int il, jl, kl;
 
108
      int nrows, ncols, nlinks;
 
109
      double dijk, denom;
 
110
      double ***W2;
 
111
      int GE, GF, GC, GX3;
 
112
      double t_ia, t_ib, t_ic, t_ja, t_jb, t_jc, t_ka, t_kb, t_kc;
 
113
      double f_ia, f_ib, f_ic, f_ja, f_jb, f_jc, f_ka, f_kb, f_kc;
 
114
      double t_jkbc, t_jkac, t_jkba, t_ikbc, t_ikac, t_ikba, t_jibc, t_jiac, t_jiba;
 
115
      double D_jkbc, D_jkac, D_jkba, D_ikbc, D_ikac, D_ikba, D_jibc, D_jiac, D_jiba;
 
116
 
 
117
      GC = C2->file.my_irrep;
 
118
      /* F and E are assumed to have same irrep */
 
119
      GF = GE =  F->file.my_irrep;
 
120
      GX3 = GC^GF;
 
121
 
 
122
      dpd_file2_mat_init(C1);
 
123
      dpd_file2_mat_rd(C1);
 
124
 
 
125
      dpd_file2_mat_init(fIJ);
 
126
      dpd_file2_mat_init(fAB);
 
127
      dpd_file2_mat_init(fIA);
 
128
      dpd_file2_mat_rd(fIJ);
 
129
      dpd_file2_mat_rd(fAB);
 
130
      dpd_file2_mat_rd(fIA);
 
131
 
 
132
      for(h=0; h < nirreps; h++) {
 
133
        dpd_buf4_mat_irrep_init(C2, h);
 
134
        dpd_buf4_mat_irrep_rd(C2, h);
 
135
 
 
136
        if(disc) {
 
137
          dpd_buf4_mat_irrep_init(D, h);
 
138
          dpd_buf4_mat_irrep_rd(D, h);
 
139
        }
 
140
 
 
141
        dpd_buf4_mat_irrep_init(E, h);
 
142
        dpd_buf4_mat_irrep_rd(E, h);
 
143
      }
 
144
 
 
145
      i = I - occ_off[Gi];
 
146
      j = J - occ_off[Gj];
 
147
      k = K - occ_off[Gk];
 
148
 
 
149
      Gij = Gji = Gi ^ Gj;
 
150
      Gik = Gki = Gi ^ Gk;
 
151
      Gjk = Gkj = Gj ^ Gk;
 
152
      Gijk = Gi ^ Gj ^ Gk;
 
153
 
 
154
      ij = C2->params->rowidx[I][J];
 
155
      ji = C2->params->rowidx[J][I];
 
156
      jk = C2->params->rowidx[J][K];
 
157
      kj = C2->params->rowidx[K][J];
 
158
      ik = C2->params->rowidx[I][K];
 
159
      ki = C2->params->rowidx[K][I];
 
160
 
 
161
      dijk = 0.0;
 
162
      if(fIJ->params->rowtot[Gi]) dijk += fIJ->matrix[Gi][i][i];
 
163
      if(fIJ->params->rowtot[Gj]) dijk += fIJ->matrix[Gj][j][j];
 
164
      if(fIJ->params->rowtot[Gk]) dijk += fIJ->matrix[Gk][k][k];
 
165
 
 
166
      W2 = (double ***) malloc(nirreps * sizeof(double **));
 
167
 
 
168
      for(Gab=0; Gab < nirreps; Gab++) {
 
169
        Gc = Gab ^ Gijk ^ GX3; /* changed */
 
170
 
 
171
        W2[Gab] = dpd_block_matrix(F->params->coltot[Gab], virtpi[Gc]);
 
172
 
 
173
        if(F->params->coltot[Gab] && virtpi[Gc]) {
 
174
          memset(W[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
175
          if(disc) memset(V[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
176
        }
 
177
      }
 
178
 
 
179
      for(Gd=0; Gd < nirreps; Gd++) {
 
180
 
 
181
        /* +t_kjcd * F_idab */
 
182
        Gid = Gi ^ Gd;
 
183
        Gab = Gid ^ GF; /* changed */
 
184
 
 
185
        Gc = Gjk ^ Gd ^ GC; /* changed */
 
186
 
 
187
        cd = C2->col_offset[Gjk][Gc];
 
188
        id = F->row_offset[Gid][I];
 
189
 
 
190
        F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
 
191
        dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
 
192
 
 
193
        nrows = F->params->coltot[Gid^GF];
 
194
        ncols = virtpi[Gc];
 
195
        nlinks = virtpi[Gd];
 
196
 
 
197
        if(nrows && ncols && nlinks)
 
198
          C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][0], nrows,
 
199
                  &(C2->matrix[Gjk][kj][cd]), nlinks, 1.0, W[Gab][0], ncols);
 
200
 
 
201
        dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
 
202
 
 
203
        /* +t_ikcd * F_jdab */
 
204
        Gjd = Gj ^ Gd;
 
205
        Gab = Gjd ^ GF; /* changed */
 
206
        Gc = Gik ^ Gd ^ GC;  /* changed */
 
207
 
 
208
        cd = C2->col_offset[Gik][Gc];
 
209
        jd = F->row_offset[Gjd][J];
 
210
 
 
211
        F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
 
212
        dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
 
213
 
 
214
        nrows = F->params->coltot[Gjd^GF];
 
215
        ncols = virtpi[Gc];
 
216
        nlinks = virtpi[Gd];
 
217
 
 
218
        if(nrows && ncols && nlinks)
 
219
          C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][0], nrows,
 
220
                  &(C2->matrix[Gik][ik][cd]), nlinks, 1.0, W[Gab][0], ncols);
 
221
 
 
222
        dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
 
223
 
 
224
        /* -t_ijcd * F_kdab */
 
225
        Gkd = Gk ^ Gd; /*changed */
 
226
        Gab = Gkd ^ GF;
 
227
        Gc = Gij ^ Gd ^ GC; 
 
228
 
 
229
        cd = C2->col_offset[Gij][Gc];
 
230
        kd = F->row_offset[Gkd][K];
 
231
 
 
232
        F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
 
233
        dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
 
234
 
 
235
        nrows = F->params->coltot[Gkd^GF];
 
236
        ncols = virtpi[Gc];
 
237
        nlinks = virtpi[Gd];
 
238
 
 
239
        if(nrows && ncols && nlinks) 
 
240
          C_DGEMM('t', 't', nrows, ncols, nlinks, -1.0, F->matrix[Gkd][0], nrows,
 
241
                  &(C2->matrix[Gij][ij][cd]), nlinks, 1.0, W[Gab][0], ncols);
 
242
 
 
243
        dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
 
244
 
 
245
      }
 
246
 
 
247
      for(Gl=0; Gl < nirreps; Gl++) {
 
248
 
 
249
        /* -t_ilab * E_jklc */
 
250
        Gil = Gi ^ Gl; /* changed */
 
251
        Gab = Gil ^ GC; 
 
252
        Gc = Gjk ^ Gl ^ GE; 
 
253
 
 
254
        lc = E->col_offset[Gjk][Gl];
 
255
        il = C2->row_offset[Gil][I];
 
256
 
 
257
        nrows = C2->params->coltot[Gil^GC];
 
258
        ncols = virtpi[Gc];
 
259
        nlinks = occpi[Gl];
 
260
 
 
261
        if(nrows && ncols && nlinks)
 
262
          C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gil][il], nrows, 
 
263
                  &(E->matrix[Gjk][jk][lc]), ncols, 1.0, W[Gab][0], ncols);
 
264
 
 
265
        /* +t_jlab * E_iklc */
 
266
        Gjl = Gj ^ Gl; /* changed */
 
267
        Gab = Gjl ^ GC;
 
268
        Gc = Gik ^ Gl ^ GE;
 
269
 
 
270
        lc = E->col_offset[Gik][Gl];
 
271
        jl = C2->row_offset[Gjl][J];
 
272
 
 
273
        nrows = C2->params->coltot[Gjl^GC];
 
274
        ncols = virtpi[Gc];
 
275
        nlinks = occpi[Gl];
 
276
 
 
277
        if(nrows && ncols && nlinks)
 
278
          C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gjl][jl], nrows, 
 
279
                  &(E->matrix[Gik][ik][lc]), ncols, 1.0, W[Gab][0], ncols);
 
280
 
 
281
        /* +t_klab * E_jilc */
 
282
        Gkl = Gk ^ Gl; /* changed! */
 
283
        Gab = Gkl ^ GC;
 
284
        Gc = Gji ^ Gl ^ GE;
 
285
 
 
286
        lc = E->col_offset[Gji][Gl];
 
287
        kl = C2->row_offset[Gkl][K];
 
288
 
 
289
        nrows = C2->params->coltot[Gkl^GC];
 
290
        ncols = virtpi[Gc];
 
291
        nlinks = occpi[Gl];
 
292
 
 
293
        if(nrows && ncols && nlinks)
 
294
          C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gkl][kl], nrows, 
 
295
                  &(E->matrix[Gji][ji][lc]), ncols, 1.0, W[Gab][0], ncols);
 
296
 
 
297
      }
 
298
 
 
299
      for(Gd=0; Gd < nirreps; Gd++) {
 
300
        /* +t_kjbd * F_idca */
 
301
        Gid = Gi ^ Gd; /* changed */
 
302
        Gca = Gid ^ GF;
 
303
        Gb = Gjk ^ Gd ^ GC;
 
304
 
 
305
        bd = C2->col_offset[Gjk][Gb];
 
306
        id = F->row_offset[Gid][I];
 
307
 
 
308
        F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
 
309
        dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
 
310
 
 
311
        nrows = F->params->coltot[Gid^GF];
 
312
        ncols = virtpi[Gb];
 
313
        nlinks = virtpi[Gd];
 
314
 
 
315
        if(nrows && ncols && nlinks)
 
316
          C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][0], nrows,
 
317
                  &(C2->matrix[Gjk][kj][bd]), nlinks, 1.0, W2[Gca][0], ncols);
 
318
 
 
319
        dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
 
320
 
 
321
        /* +t_ikbd * F_jdca */
 
322
        Gjd = Gj ^ Gd; 
 
323
        Gca = Gjd ^ GF ; 
 
324
        Gb = Gik ^ Gd ^ GC;      
 
325
 
 
326
        bd = C2->col_offset[Gik][Gb];
 
327
        jd = F->row_offset[Gjd][J];
 
328
 
 
329
        F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
 
330
        dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
 
331
 
 
332
        nrows = F->params->coltot[Gjd^GF];
 
333
        ncols = virtpi[Gb];
 
334
        nlinks = virtpi[Gd];
 
335
 
 
336
        if(nrows && ncols && nlinks)
 
337
          C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][0], nrows,
 
338
                  &(C2->matrix[Gik][ik][bd]), nlinks, 1.0, W2[Gca][0], ncols);
 
339
 
 
340
        dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
 
341
 
 
342
        /* -t_ijbd * F_kdca */
 
343
        Gkd = Gk ^ Gd; 
 
344
        Gca = Gkd ^ GF; 
 
345
        Gb = Gij ^ Gd ^ GC;      
 
346
 
 
347
        bd = C2->col_offset[Gij][Gb];
 
348
        kd = F->row_offset[Gkd][K];
 
349
 
 
350
        F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
 
351
        dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
 
352
 
 
353
        nrows = F->params->coltot[Gkd^GF];
 
354
        ncols = virtpi[Gb];
 
355
        nlinks = virtpi[Gd];
 
356
 
 
357
        if(nrows && ncols && nlinks)
 
358
          C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gkd][0], nrows,
 
359
                  &(C2->matrix[Gij][ij][bd]), nlinks, 1.0, W2[Gca][0], ncols);
 
360
 
 
361
        dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
 
362
      }
 
363
 
 
364
      for(Gl=0; Gl < nirreps; Gl++) {
 
365
        /* -t_ilca * E_jklb */
 
366
        Gil = Gi ^ Gl; 
 
367
        Gca = Gil ^ GC;
 
368
        Gb = Gjk ^ Gl ^ GE;     
 
369
 
 
370
        lb = E->col_offset[Gjk][Gl];
 
371
        il = C2->row_offset[Gil][I];
 
372
 
 
373
        nrows = C2->params->coltot[Gil^GC];
 
374
        ncols = virtpi[Gb];
 
375
        nlinks = occpi[Gl];
 
376
 
 
377
        if(nrows && ncols && nlinks)
 
378
          C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gil][il], nrows, 
 
379
                  &(E->matrix[Gjk][jk][lb]), ncols, 1.0, W2[Gca][0], ncols);
 
380
 
 
381
        /* +t_jlca * E_iklb */
 
382
        Gjl = Gj ^ Gl;
 
383
        Gca = Gjl ^ GC;
 
384
        Gb = Gik ^ Gl^ GE; 
 
385
 
 
386
        lb = E->col_offset[Gik][Gl];
 
387
        jl = C2->row_offset[Gjl][J];
 
388
 
 
389
        nrows = C2->params->coltot[Gjl^GC];
 
390
        ncols = virtpi[Gb];
 
391
        nlinks = occpi[Gl];
 
392
 
 
393
        if(nrows && ncols && nlinks)
 
394
          C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gjl][jl], nrows, 
 
395
                  &(E->matrix[Gik][ik][lb]), ncols, 1.0, W2[Gca][0], ncols);
 
396
 
 
397
        /* +t_klca * E_jilb */
 
398
        Gkl = Gk ^ Gl; 
 
399
        Gca = Gkl ^ GC; 
 
400
        Gb = Gji ^ Gl ^ GE;      
 
401
 
 
402
        lb = E->col_offset[Gji][Gl];
 
403
        kl = C2->row_offset[Gkl][K];
 
404
 
 
405
        nrows = C2->params->coltot[Gkl^GC];
 
406
        ncols = virtpi[Gb];
 
407
        nlinks = occpi[Gl];
 
408
 
 
409
        if(nrows && ncols && nlinks)
 
410
          C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gkl][kl], nrows, 
 
411
                  &(E->matrix[Gji][ji][lb]), ncols, 1.0, W2[Gca][0], ncols);
 
412
      }
 
413
 
 
414
      dpd_3d_sort(W2, W, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
415
                  F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
 
416
                  vir_off, virtpi, vir_off, F->params->colidx, bca, 1);
 
417
 
 
418
      for(Gab=0; Gab < nirreps; Gab++) {
 
419
        Gc = Gab ^ Gijk ^ GX3; /* changed */
 
420
        if(F->params->coltot[Gab] && virtpi[Gc]) {
 
421
          memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
422
        }
 
423
      }
 
424
 
 
425
      for(Gd=0; Gd < nirreps; Gd++) {
 
426
        /* -t_kjad * F_idcb */
 
427
        Gid = Gi ^ Gd;
 
428
        Gcb = Gid ^ GF;
 
429
        Ga = Gkj ^ Gd ^ GC;     
 
430
 
 
431
        ad = C2->col_offset[Gkj][Ga];
 
432
        id = F->row_offset[Gid][I];
 
433
 
 
434
        F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
 
435
        dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
 
436
 
 
437
        nrows = F->params->coltot[Gid^GF];
 
438
        ncols = virtpi[Ga];
 
439
        nlinks = virtpi[Gd];
 
440
 
 
441
        if(nrows && ncols && nlinks)
 
442
          C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gid][0], nrows,
 
443
                  &(C2->matrix[Gkj][kj][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
 
444
 
 
445
        dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
 
446
 
 
447
        /* -t_ikad * F_jdcb */
 
448
        Gjd = Gj ^ Gd;
 
449
        Gcb = Gjd ^ GF;
 
450
        Ga = Gik ^ Gd ^ GC;     
 
451
 
 
452
        ad = C2->col_offset[Gik][Ga];
 
453
        jd = F->row_offset[Gjd][J];
 
454
 
 
455
        F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
 
456
        dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
 
457
 
 
458
        nrows = F->params->coltot[Gjd^GF];
 
459
        ncols = virtpi[Ga];
 
460
        nlinks = virtpi[Gd];
 
461
 
 
462
        if(nrows && ncols && nlinks)
 
463
          C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gjd][0], nrows,
 
464
                  &(C2->matrix[Gik][ik][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
 
465
 
 
466
        dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
 
467
 
 
468
        /* +t_ijad * F_kdcb */
 
469
        Gkd = Gk ^ Gd;
 
470
        Gcb = Gkd ^ GF;
 
471
        Ga = Gij ^ Gd ^ GC;     
 
472
 
 
473
        ad = C2->col_offset[Gij][Ga];
 
474
        kd = F->row_offset[Gkd][K];
 
475
 
 
476
        F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
 
477
        dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
 
478
 
 
479
        nrows = F->params->coltot[Gkd^GF];
 
480
        ncols = virtpi[Ga];
 
481
        nlinks = virtpi[Gd];
 
482
 
 
483
        if(nrows && ncols && nlinks)
 
484
          C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gkd][0], nrows,
 
485
                  &(C2->matrix[Gij][ij][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
 
486
 
 
487
        dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
 
488
 
 
489
      }
 
490
 
 
491
      for(Gl=0; Gl < nirreps; Gl++) {
 
492
        /* +t_ilcb * E_jkla */
 
493
        Gil = Gi ^ Gl; 
 
494
        Gcb  = Gil ^ GC; 
 
495
        Ga = Gjk ^ Gl ^ GE;       
 
496
 
 
497
        la = E->col_offset[Gjk][Gl];
 
498
        il = C2->row_offset[Gil][I];
 
499
 
 
500
        nrows = C2->params->coltot[Gil^GC];
 
501
        ncols = virtpi[Ga];
 
502
        nlinks = occpi[Gl];
 
503
 
 
504
        if(nrows && ncols && nlinks)
 
505
          C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gil][il], nrows, 
 
506
                  &(E->matrix[Gjk][jk][la]), ncols, 1.0, W2[Gcb][0], ncols);
 
507
 
 
508
        /* -t_jlcb * E_ikla */
 
509
        Gjl = Gj ^ Gl; 
 
510
        Gcb = Gjl ^ GC; 
 
511
        Ga = Gik ^ Gl ^ GE;       
 
512
 
 
513
        la = E->col_offset[Gik][Gl];
 
514
        jl = C2->row_offset[Gjl][J];
 
515
 
 
516
        nrows = C2->params->coltot[Gjl^GC];
 
517
        ncols = virtpi[Ga];
 
518
        nlinks = occpi[Gl];
 
519
 
 
520
        if(nrows && ncols && nlinks)
 
521
          C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gjl][jl], nrows, 
 
522
                  &(E->matrix[Gik][ik][la]), ncols, 1.0, W2[Gcb][0], ncols);
 
523
 
 
524
        /* -t_klcb * E_jila */
 
525
        Gkl = Gk ^ Gl; 
 
526
        Gcb = Gkl ^ GC; 
 
527
        Ga = Gji ^ Gl ^ GE;      
 
528
 
 
529
        la = E->col_offset[Gji][Gl];
 
530
        kl = C2->row_offset[Gkl][K];
 
531
 
 
532
        nrows = C2->params->coltot[Gkl^GC];
 
533
        ncols = virtpi[Ga];
 
534
        nlinks = occpi[Gl];
 
535
 
 
536
        if(nrows && ncols && nlinks)
 
537
          C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gkl][kl], nrows, 
 
538
                  &(E->matrix[Gji][ji][la]), ncols, 1.0, W2[Gcb][0], ncols);
 
539
 
 
540
      }
 
541
 
 
542
      dpd_3d_sort(W2, W, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
543
                  F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
 
544
                  vir_off, virtpi, vir_off, F->params->colidx, cba, 1);
 
545
 
 
546
      /**** Compute disconnected T3s for given ijk ****/
 
547
      if(disc) {
 
548
        for(Gab=0; Gab < nirreps; Gab++) {
 
549
          Gc = Gab ^ Gijk;
 
550
 
 
551
          for(ab=0; ab < F->params->coltot[Gab]; ab++) {
 
552
            A = F->params->colorb[Gab][ab][0];
 
553
            Ga = F->params->rsym[A];
 
554
            a = A - vir_off[Ga];
 
555
            B = F->params->colorb[Gab][ab][1];
 
556
            Gb = F->params->ssym[B];
 
557
            b = B - vir_off[Gb];
 
558
 
 
559
            Gbc = Gb ^ Gc;
 
560
            Gac = Ga ^ Gc;
 
561
 
 
562
            for(c=0; c < virtpi[Gc]; c++) {
 
563
              C = vir_off[Gc] + c;
 
564
 
 
565
              bc = D->params->colidx[B][C];
 
566
              ac = D->params->colidx[A][C];
 
567
 
 
568
              /* +t_ia * D_jkbc + f_ia * t_jkbc */
 
569
              if(Gi == Ga && Gjk == Gbc) {
 
570
                t_ia = D_jkbc = f_ia = t_jkbc = 0.0;
 
571
 
 
572
                if(C1->params->rowtot[Gi] && C1->params->coltot[Gi]) {
 
573
                  t_ia = C1->matrix[Gi][i][a];
 
574
                  f_ia = fIA->matrix[Gi][i][a];
 
575
                }
 
576
 
 
577
                if(D->params->rowtot[Gjk] && D->params->coltot[Gjk]) {
 
578
                  D_jkbc = D->matrix[Gjk][jk][bc];
 
579
                  t_jkbc = C2->matrix[Gjk][jk][bc];
 
580
                }
 
581
 
 
582
                V[Gab][ab][c] += t_ia * D_jkbc + f_ia * t_jkbc;
 
583
              }
 
584
 
 
585
              /* -t_ib * D_jkac - f_ib * t_jkac */
 
586
              if(Gi == Gb && Gjk == Gac) {
 
587
                t_ib = D_jkac = f_ib = t_jkac = 0.0;
 
588
 
 
589
                if(C1->params->rowtot[Gi] && C1->params->coltot[Gi]) {
 
590
                  t_ib = C1->matrix[Gi][i][b];
 
591
                  f_ib = fIA->matrix[Gi][i][b];
 
592
                }
 
593
 
 
594
                if(D->params->rowtot[Gjk] && D->params->coltot[Gjk]) {
 
595
                  D_jkac = D->matrix[Gjk][jk][ac];
 
596
                  t_jkac = C2->matrix[Gjk][jk][ac];
 
597
                }
 
598
 
 
599
                V[Gab][ab][c] -= t_ib * D_jkac + f_ib * t_jkac;
 
600
              }
 
601
 
 
602
              /* +t_ic * D_jkab + f_ic * t_jkba */
 
603
              if(Gi == Gc && Gjk == Gab) {
 
604
                t_ic = D_jkba = f_ic = t_jkba = 0.0;
 
605
 
 
606
                if(C1->params->rowtot[Gi] && C1->params->coltot[Gi]) {
 
607
                  t_ic = C1->matrix[Gi][i][c];
 
608
                  f_ic = fIA->matrix[Gi][i][c];
 
609
                }
 
610
 
 
611
                if(D->params->rowtot[Gjk] && D->params->coltot[Gjk]) {
 
612
                  D_jkba = D->matrix[Gjk][jk][ab];
 
613
                  t_jkba = C2->matrix[Gjk][jk][ab];
 
614
                }
 
615
 
 
616
                V[Gab][ab][c] += t_ic * D_jkba + f_ic * t_jkba;
 
617
              }
 
618
 
 
619
              /* -t_ja * D_ikbc - f_ja * t_ikbc*/
 
620
              if(Gj == Ga && Gik == Gbc) {
 
621
                t_ja = D_ikbc = f_ja = t_ikbc = 0.0;
 
622
 
 
623
                if(C1->params->rowtot[Gj] && C1->params->coltot[Gj]) {
 
624
                  t_ja = C1->matrix[Gj][j][a];
 
625
                  f_ja = fIA->matrix[Gj][j][a];
 
626
                }
 
627
 
 
628
                if(D->params->rowtot[Gik] && D->params->coltot[Gik]) {
 
629
                  D_ikbc = D->matrix[Gik][ik][bc];
 
630
                  t_ikbc = C2->matrix[Gik][ik][bc];
 
631
                }
 
632
 
 
633
                V[Gab][ab][c] -= t_ja * D_ikbc + f_ja * t_ikbc;
 
634
              }
 
635
 
 
636
              /* +t_jb * D_ikac + f_jb * t_ikac */
 
637
              if(Gj == Gb && Gik == Gac) {
 
638
                t_jb = D_ikac = f_jb = t_ikac = 0.0;
 
639
 
 
640
                if(C1->params->rowtot[Gj] && C1->params->coltot[Gj]) {
 
641
                  t_jb = C1->matrix[Gj][j][b];
 
642
                  f_jb = fIA->matrix[Gj][j][b];
 
643
                }
 
644
 
 
645
                if(D->params->rowtot[Gik] && D->params->coltot[Gik]) {
 
646
                  D_ikac = D->matrix[Gik][ik][ac];
 
647
                  t_ikac = C2->matrix[Gik][ik][ac];
 
648
                }
 
649
 
 
650
                V[Gab][ab][c] += t_jb * D_ikac + f_jb * t_ikac;
 
651
              }
 
652
 
 
653
              /* -t_jc * D_ikba - f_jc * t_ikba */
 
654
              if(Gj == Gc && Gik == Gab) {
 
655
                t_jc = D_ikba = f_jc = t_ikba = 0.0;
 
656
 
 
657
                if(C1->params->rowtot[Gj] && C1->params->coltot[Gj]) {
 
658
                  t_jc = C1->matrix[Gj][j][c];
 
659
                  f_jc = fIA->matrix[Gj][j][c];
 
660
                }
 
661
 
 
662
                if(D->params->rowtot[Gik] && D->params->coltot[Gik]) {
 
663
                  D_ikba = D->matrix[Gik][ik][ab];
 
664
                  t_ikba = C2->matrix[Gik][ik][ab];
 
665
                }
 
666
 
 
667
                V[Gab][ab][c] -= t_jc * D_ikba + f_jc * t_ikba;
 
668
              }
 
669
 
 
670
              /* -t_ka * D_jibc - f_ka * t_jibc */
 
671
              if(Gk == Ga && Gji == Gbc) {
 
672
                t_ka = D_jibc = f_ka = t_jibc = 0.0;
 
673
 
 
674
                if(C1->params->rowtot[Gk] && C1->params->coltot[Gk]) {
 
675
                  t_ka = C1->matrix[Gk][k][a];
 
676
                  f_ka = fIA->matrix[Gk][k][a];
 
677
                }
 
678
 
 
679
                if(D->params->rowtot[Gji] && D->params->coltot[Gji]) {
 
680
                  D_jibc = D->matrix[Gji][ji][bc];
 
681
                  t_jibc = C2->matrix[Gji][ji][bc];
 
682
                }
 
683
 
 
684
                V[Gab][ab][c] -= t_ka * D_jibc + f_ka * t_jibc;
 
685
              }
 
686
 
 
687
              /* +t_kb * D_jiac + f_kb * t_jiac */
 
688
              if(Gk == Gb && Gji == Gac) {
 
689
                t_kb = D_jiac = f_kb = t_jiac = 0.0;
 
690
 
 
691
                if(C1->params->rowtot[Gk] && C1->params->coltot[Gk]) {
 
692
                  t_kb = C1->matrix[Gk][k][b];
 
693
                  f_kb = fIA->matrix[Gk][k][b];
 
694
                }
 
695
 
 
696
                if(D->params->rowtot[Gji] && D->params->coltot[Gji]) {
 
697
                  D_jiac = D->matrix[Gji][ji][ac];
 
698
                  t_jiac = C2->matrix[Gji][ji][ac];
 
699
                }
 
700
 
 
701
                V[Gab][ab][c] += t_kb * D_jiac + f_kb * t_jiac;
 
702
              }
 
703
 
 
704
              /* -t_kc * D_jiab - f_kc * t_jiba*/
 
705
              if(Gk == Gc && Gji == Gab) {
 
706
                t_kc = D_jiba = f_kc = t_jiba = 0.0;
 
707
 
 
708
                if(C1->params->rowtot[Gk] && C1->params->coltot[Gk]) {
 
709
                  t_kc = C1->matrix[Gk][k][c];
 
710
                  f_kc = fIA->matrix[Gk][k][c];
 
711
                }
 
712
 
 
713
                if(D->params->rowtot[Gji] && D->params->coltot[Gji]) {
 
714
                  D_jiba = D->matrix[Gji][ji][ab];
 
715
                  t_jiba = C2->matrix[Gji][ji][ab];
 
716
                }
 
717
 
 
718
                V[Gab][ab][c] -= t_kc * D_jiba + f_kc * t_jiba;
 
719
              }
 
720
 
 
721
            } /* c */
 
722
          } /* ab */
 
723
        } /* Gab */
 
724
 
 
725
 
 
726
      } /**** Disconnected T3 complete ****/
 
727
 
 
728
      for(Gab=0; Gab < nirreps; Gab++) {
 
729
        Gc = Gab ^ Gijk ^ GX3; /* assumes totally symmetric! */
 
730
        for(ab=0; ab < F->params->coltot[Gab]; ab++) {
 
731
          A = F->params->colorb[Gab][ab][0];
 
732
          Ga = F->params->rsym[A];
 
733
          a = A - vir_off[Ga];
 
734
          B = F->params->colorb[Gab][ab][1];
 
735
          Gb = F->params->ssym[B];
 
736
          b = B - vir_off[Gb];
 
737
 
 
738
          for(c=0; c < virtpi[Gc]; c++) {
 
739
            denom = dijk;
 
740
            if(fAB->params->rowtot[Ga]) denom -= fAB->matrix[Ga][a][a];
 
741
            if(fAB->params->rowtot[Gb]) denom -= fAB->matrix[Gb][b][b];
 
742
            if(fAB->params->rowtot[Gc]) denom -= fAB->matrix[Gc][c][c];
 
743
 
 
744
            W[Gab][ab][c] /= (omega + denom);
 
745
            if(disc) V[Gab][ab][c] /= (omega + denom);
 
746
 
 
747
          } /* c */
 
748
        } /* ab */
 
749
      } /* Gab */
 
750
 
 
751
      for(Gab=0; Gab < nirreps; Gab++) {
 
752
        Gc = Gab ^ Gijk ^ GX3; /* changed */
 
753
        dpd_free_block(W2[Gab], F->params->coltot[Gab], virtpi[Gc]);
 
754
      }
 
755
      free(W2);
 
756
 
 
757
      dpd_file2_mat_close(fIJ);
 
758
      dpd_file2_mat_close(fAB);
 
759
      if(disc) {
 
760
        dpd_file2_mat_close(fIA);
 
761
        dpd_file2_mat_close(C1);
 
762
      }
 
763
 
 
764
      for(h=0; h < nirreps; h++) {
 
765
        dpd_buf4_mat_irrep_close(C2, h);
 
766
        if(disc) dpd_buf4_mat_irrep_close(D, h);
 
767
        dpd_buf4_mat_irrep_close(E, h);
 
768
      }
 
769
    }
 
770
 
 
771
  }}