← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

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

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

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

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libdpd/T3_AAA.c
1
 
 
2
 
/* T3_RHF_AAA(): Computes all T3(IJK,ABC) amplitudes for a given I, J,
3
 
** K combination for input C2, F, and E intermediates.  This function
4
 
** will work for AAA or BBB spin cases, with either RHF/ROHF or UHF
5
 
** orbitals.
6
 
**
7
 
** Arguments: 
8
 
**
9
 
**   double ***W1: The target triples amplitudes in an nirreps x AB x
10
 
**   C array.  The memory for this must be allocated externally.
11
 
**
12
 
**   int nirreps: Number of irreps.
13
 
**
14
 
**   int I: Absolute index of orbital I.
15
 
**
16
 
**   int Gi: Irrep of I.
17
 
**
18
 
**   int J: Absolute index of orbital J.
19
 
**
20
 
**   int Gj: Irrep of J.
21
 
**
22
 
**   int K: Absolute index of orbital K.
23
 
**
24
 
**   int Gk: Irrep of K.
25
 
**
26
 
**   dpdbuf4 *C2: Pointer to dpd buffer for double excitation amps,
27
 
**   ordered (IJ,AB).
28
 
**
29
 
**   dpdbuf4 *F: Pointer to dpd buffer for three-virtual-index
30
 
**   intermediate, ordered (IA,BC).
31
 
**
32
 
**   dpdbuf4 *E: Pointer to dpd buffer for three-occupied-index
33
 
**   intermediate, ordered (IJ,KA).
34
 
**
35
 
**   dpdfile2 *fIJ: Pointer to the dpd file2 for the occ-occ block of
36
 
**   the Fock matrix (or other appropriate one-electron operator).
37
 
**   
38
 
**   dpdfile2 *fAB: Pointer to the dpd file2 for the vir-vir block of
39
 
**   the Fock matrix (or other appropriate one-electron operator).
40
 
**   
41
 
**   int *occpi: Number of occupied orbitals per irrep lookup array.
42
 
**
43
 
**   int *occ_off: Offset lookup for translating between absolute and
44
 
**   relative orbital indices for occupied space.
45
 
**
46
 
**   int *virtpi: Number of virtual orbitals per irrep lookup array.
47
 
**
48
 
**   int *vir_off: Offset lookup for translating between absolute and
49
 
**   relative orbital indices for virtual space.
50
 
**
51
 
**   double omega: a constant to add to the final denominators -
52
 
**   needed for CC3 EOM
53
 
**
54
 
** TDC, July 2004
55
 
*/
56
 
 
57
 
#include <stdio.h>
58
 
#include <stdlib.h>
59
 
#include <string.h>
60
 
#include <math.h>
61
 
#include <libdpd/dpd.h>
62
 
#include <ccfiles.h>
63
 
 
64
 
void T3_AAA(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk, 
65
 
            dpdbuf4 *C2, dpdbuf4 *F, dpdbuf4 *E, dpdfile2 *fIJ, dpdfile2 *fAB, 
66
 
            int *occpi, int *occ_off, int *virtpi, int *vir_off, double omega)
67
 
{
68
 
  int h;
69
 
  int i, j, k;
70
 
  int ij, ji, ik, ki, jk, kj;
71
 
  int Gij, Gji, Gik, Gki, Gjk, Gkj, Gijk;
72
 
  int Ga, Gb, Gc;
73
 
  int Gd, Gl;
74
 
  int Gid, Gjd, Gkd;
75
 
  int Gab, Gcb, Gca;
76
 
  int Gla, Glb, Glc;
77
 
  int Gil, Gjl, Gkl;
78
 
  int a, b, c, A, B, C;
79
 
  int ab;
80
 
  int cd, bd, ad;
81
 
  int id, jd, kd;
82
 
  int la, lb, lc;
83
 
  int il, jl, kl;
84
 
  int nrows, ncols, nlinks;
85
 
  double dijk, denom;
86
 
  double ***W2;
87
 
  int GE, GF, GC, GX3;
88
 
 
89
 
  GC = C2->file.my_irrep;
90
 
  /* F and E are assumed to have same irrep */
91
 
  GF = GE =  F->file.my_irrep;
92
 
  GX3 = GC^GF;
93
 
 
94
 
  dpd_file2_mat_init(fIJ);
95
 
  dpd_file2_mat_init(fAB);
96
 
  dpd_file2_mat_rd(fIJ);
97
 
  dpd_file2_mat_rd(fAB);
98
 
 
99
 
  for(h=0; h < nirreps; h++) {
100
 
    dpd_buf4_mat_irrep_init(C2, h);
101
 
    dpd_buf4_mat_irrep_rd(C2, h);
102
 
 
103
 
    dpd_buf4_mat_irrep_init(E, h);
104
 
    dpd_buf4_mat_irrep_rd(E, h);
105
 
  }
106
 
 
107
 
  i = I - occ_off[Gi];
108
 
  j = J - occ_off[Gj];
109
 
  k = K - occ_off[Gk];
110
 
 
111
 
  Gij = Gji = Gi ^ Gj;
112
 
  Gik = Gki = Gi ^ Gk;
113
 
  Gjk = Gkj = Gj ^ Gk;
114
 
  Gijk = Gi ^ Gj ^ Gk;
115
 
 
116
 
  ij = C2->params->rowidx[I][J];
117
 
  ji = C2->params->rowidx[J][I];
118
 
  jk = C2->params->rowidx[J][K];
119
 
  kj = C2->params->rowidx[K][J];
120
 
  ik = C2->params->rowidx[I][K];
121
 
  ki = C2->params->rowidx[K][I];
122
 
 
123
 
  dijk = 0.0;
124
 
  if(fIJ->params->rowtot[Gi]) dijk += fIJ->matrix[Gi][i][i];
125
 
  if(fIJ->params->rowtot[Gj]) dijk += fIJ->matrix[Gj][j][j];
126
 
  if(fIJ->params->rowtot[Gk]) dijk += fIJ->matrix[Gk][k][k];
127
 
 
128
 
  W2 = (double ***) malloc(nirreps * sizeof(double **));
129
 
 
130
 
  for(Gab=0; Gab < nirreps; Gab++) {
131
 
    Gc = Gab ^ Gijk ^ GX3; /* changed */
132
 
 
133
 
    W2[Gab] = dpd_block_matrix(F->params->coltot[Gab], virtpi[Gc]);
134
 
 
135
 
    if(F->params->coltot[Gab] && virtpi[Gc]) {
136
 
      memset(W1[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
137
 
    }
138
 
  }
139
 
 
140
 
  for(Gd=0; Gd < nirreps; Gd++) {
141
 
 
142
 
    /* +t_kjcd * F_idab */
143
 
    Gid = Gi ^ Gd;
144
 
    Gab = Gid ^ GF; /* changed */
145
 
 
146
 
    Gc = Gjk ^ Gd ^ GC; /* changed */
147
 
 
148
 
    cd = C2->col_offset[Gjk][Gc];
149
 
    id = F->row_offset[Gid][I];
150
 
 
151
 
    F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
152
 
    dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
153
 
 
154
 
    nrows = F->params->coltot[Gid^GF];
155
 
    ncols = virtpi[Gc];
156
 
    nlinks = virtpi[Gd];
157
 
 
158
 
    if(nrows && ncols && nlinks)
159
 
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][0], nrows,
160
 
              &(C2->matrix[Gjk][kj][cd]), nlinks, 1.0, W1[Gab][0], ncols);
161
 
 
162
 
    dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
163
 
 
164
 
    /* +t_ikcd * F_jdab */
165
 
    Gjd = Gj ^ Gd;
166
 
    Gab = Gjd ^ GF; /* changed */
167
 
    Gc = Gik ^ Gd ^ GC;  /* changed */
168
 
 
169
 
    cd = C2->col_offset[Gik][Gc];
170
 
    jd = F->row_offset[Gjd][J];
171
 
 
172
 
    F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
173
 
    dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
174
 
 
175
 
    nrows = F->params->coltot[Gjd^GF];
176
 
    ncols = virtpi[Gc];
177
 
    nlinks = virtpi[Gd];
178
 
 
179
 
    if(nrows && ncols && nlinks)
180
 
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][0], nrows,
181
 
              &(C2->matrix[Gik][ik][cd]), nlinks, 1.0, W1[Gab][0], ncols);
182
 
 
183
 
    dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
184
 
 
185
 
    /* -t_ijcd * F_kdab */
186
 
    Gkd = Gk ^ Gd; /*changed */
187
 
    Gab = Gkd ^ GF;
188
 
    Gc = Gij ^ Gd ^ GC; 
189
 
 
190
 
    cd = C2->col_offset[Gij][Gc];
191
 
    kd = F->row_offset[Gkd][K];
192
 
 
193
 
    F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
194
 
    dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
195
 
 
196
 
    nrows = F->params->coltot[Gkd^GF];
197
 
    ncols = virtpi[Gc];
198
 
    nlinks = virtpi[Gd];
199
 
 
200
 
    if(nrows && ncols && nlinks) 
201
 
      C_DGEMM('t', 't', nrows, ncols, nlinks, -1.0, F->matrix[Gkd][0], nrows,
202
 
              &(C2->matrix[Gij][ij][cd]), nlinks, 1.0, W1[Gab][0], ncols);
203
 
 
204
 
    dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
205
 
 
206
 
  }
207
 
 
208
 
  for(Gl=0; Gl < nirreps; Gl++) {
209
 
 
210
 
    /* -t_ilab * E_jklc */
211
 
    Gil = Gi ^ Gl; /* changed */
212
 
    Gab = Gil ^ GC; 
213
 
    Gc = Gjk ^ Gl ^ GE; 
214
 
 
215
 
    lc = E->col_offset[Gjk][Gl];
216
 
    il = C2->row_offset[Gil][I];
217
 
 
218
 
    nrows = C2->params->coltot[Gil^GC];
219
 
    ncols = virtpi[Gc];
220
 
    nlinks = occpi[Gl];
221
 
 
222
 
    if(nrows && ncols && nlinks)
223
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gil][il], nrows, 
224
 
              &(E->matrix[Gjk][jk][lc]), ncols, 1.0, W1[Gab][0], ncols);
225
 
 
226
 
    /* +t_jlab * E_iklc */
227
 
    Gjl = Gj ^ Gl; /* changed */
228
 
    Gab = Gjl ^ GC;
229
 
    Gc = Gik ^ Gl ^ GE;
230
 
 
231
 
    lc = E->col_offset[Gik][Gl];
232
 
    jl = C2->row_offset[Gjl][J];
233
 
 
234
 
    nrows = C2->params->coltot[Gjl^GC];
235
 
    ncols = virtpi[Gc];
236
 
    nlinks = occpi[Gl];
237
 
 
238
 
    if(nrows && ncols && nlinks)
239
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gjl][jl], nrows, 
240
 
              &(E->matrix[Gik][ik][lc]), ncols, 1.0, W1[Gab][0], ncols);
241
 
 
242
 
    /* +t_klab * E_jilc */
243
 
    Gkl = Gk ^ Gl; /* changed! */
244
 
    Gab = Gkl ^ GC;
245
 
    Gc = Gji ^ Gl ^ GE;
246
 
 
247
 
    lc = E->col_offset[Gji][Gl];
248
 
    kl = C2->row_offset[Gkl][K];
249
 
 
250
 
    nrows = C2->params->coltot[Gkl^GC];
251
 
    ncols = virtpi[Gc];
252
 
    nlinks = occpi[Gl];
253
 
 
254
 
    if(nrows && ncols && nlinks)
255
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gkl][kl], nrows, 
256
 
              &(E->matrix[Gji][ji][lc]), ncols, 1.0, W1[Gab][0], ncols);
257
 
 
258
 
  }
259
 
 
260
 
  for(Gd=0; Gd < nirreps; Gd++) {
261
 
    /* +t_kjbd * F_idca */
262
 
    Gid = Gi ^ Gd; /* changed */
263
 
    Gca = Gid ^ GF;
264
 
    Gb = Gjk ^ Gd ^ GC;
265
 
 
266
 
    bd = C2->col_offset[Gjk][Gb];
267
 
    id = F->row_offset[Gid][I];
268
 
 
269
 
    F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
270
 
    dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
271
 
 
272
 
    nrows = F->params->coltot[Gid^GF];
273
 
    ncols = virtpi[Gb];
274
 
    nlinks = virtpi[Gd];
275
 
 
276
 
    if(nrows && ncols && nlinks)
277
 
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][0], nrows,
278
 
              &(C2->matrix[Gjk][kj][bd]), nlinks, 1.0, W2[Gca][0], ncols);
279
 
 
280
 
    dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
281
 
 
282
 
    /* +t_ikbd * F_jdca */
283
 
    Gjd = Gj ^ Gd; 
284
 
    Gca = Gjd ^ GF ; 
285
 
    Gb = Gik ^ Gd ^ GC;      
286
 
 
287
 
    bd = C2->col_offset[Gik][Gb];
288
 
    jd = F->row_offset[Gjd][J];
289
 
 
290
 
    F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
291
 
    dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
292
 
 
293
 
    nrows = F->params->coltot[Gjd^GF];
294
 
    ncols = virtpi[Gb];
295
 
    nlinks = virtpi[Gd];
296
 
 
297
 
    if(nrows && ncols && nlinks)
298
 
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][0], nrows,
299
 
              &(C2->matrix[Gik][ik][bd]), nlinks, 1.0, W2[Gca][0], ncols);
300
 
 
301
 
    dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
302
 
 
303
 
    /* -t_ijbd * F_kdca */
304
 
    Gkd = Gk ^ Gd; 
305
 
    Gca = Gkd ^ GF; 
306
 
    Gb = Gij ^ Gd ^ GC;      
307
 
 
308
 
    bd = C2->col_offset[Gij][Gb];
309
 
    kd = F->row_offset[Gkd][K];
310
 
 
311
 
    F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
312
 
    dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
313
 
 
314
 
    nrows = F->params->coltot[Gkd^GF];
315
 
    ncols = virtpi[Gb];
316
 
    nlinks = virtpi[Gd];
317
 
 
318
 
    if(nrows && ncols && nlinks)
319
 
      C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gkd][0], nrows,
320
 
              &(C2->matrix[Gij][ij][bd]), nlinks, 1.0, W2[Gca][0], ncols);
321
 
 
322
 
    dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
323
 
  }
324
 
 
325
 
  for(Gl=0; Gl < nirreps; Gl++) {
326
 
    /* -t_ilca * E_jklb */
327
 
    Gil = Gi ^ Gl; 
328
 
    Gca = Gil ^ GC;
329
 
    Gb = Gjk ^ Gl ^ GE;     
330
 
 
331
 
    lb = E->col_offset[Gjk][Gl];
332
 
    il = C2->row_offset[Gil][I];
333
 
 
334
 
    nrows = C2->params->coltot[Gil^GC];
335
 
    ncols = virtpi[Gb];
336
 
    nlinks = occpi[Gl];
337
 
 
338
 
    if(nrows && ncols && nlinks)
339
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gil][il], nrows, 
340
 
              &(E->matrix[Gjk][jk][lb]), ncols, 1.0, W2[Gca][0], ncols);
341
 
 
342
 
    /* +t_jlca * E_iklb */
343
 
    Gjl = Gj ^ Gl;
344
 
    Gca = Gjl ^ GC;
345
 
    Gb = Gik ^ Gl^ GE; 
346
 
 
347
 
    lb = E->col_offset[Gik][Gl];
348
 
    jl = C2->row_offset[Gjl][J];
349
 
 
350
 
    nrows = C2->params->coltot[Gjl^GC];
351
 
    ncols = virtpi[Gb];
352
 
    nlinks = occpi[Gl];
353
 
 
354
 
    if(nrows && ncols && nlinks)
355
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gjl][jl], nrows, 
356
 
              &(E->matrix[Gik][ik][lb]), ncols, 1.0, W2[Gca][0], ncols);
357
 
 
358
 
    /* +t_klca * E_jilb */
359
 
    Gkl = Gk ^ Gl; 
360
 
    Gca = Gkl ^ GC; 
361
 
    Gb = Gji ^ Gl ^ GE;      
362
 
 
363
 
    lb = E->col_offset[Gji][Gl];
364
 
    kl = C2->row_offset[Gkl][K];
365
 
 
366
 
    nrows = C2->params->coltot[Gkl^GC];
367
 
    ncols = virtpi[Gb];
368
 
    nlinks = occpi[Gl];
369
 
 
370
 
    if(nrows && ncols && nlinks)
371
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gkl][kl], nrows, 
372
 
              &(E->matrix[Gji][ji][lb]), ncols, 1.0, W2[Gca][0], ncols);
373
 
  }
374
 
 
375
 
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
376
 
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
377
 
              vir_off, virtpi, vir_off, F->params->colidx, bca, 1);
378
 
 
379
 
  for(Gab=0; Gab < nirreps; Gab++) {
380
 
    Gc = Gab ^ Gijk ^ GX3; /* changed */
381
 
    if(F->params->coltot[Gab] && virtpi[Gc]) {
382
 
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
383
 
    }
384
 
  }
385
 
 
386
 
  for(Gd=0; Gd < nirreps; Gd++) {
387
 
    /* -t_kjad * F_idcb */
388
 
    Gid = Gi ^ Gd;
389
 
    Gcb = Gid ^ GF;
390
 
    Ga = Gkj ^ Gd ^ GC;     
391
 
 
392
 
    ad = C2->col_offset[Gkj][Ga];
393
 
    id = F->row_offset[Gid][I];
394
 
 
395
 
    F->matrix[Gid] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gid^GF]);
396
 
    dpd_buf4_mat_irrep_rd_block(F, Gid, id, virtpi[Gd]);
397
 
 
398
 
    nrows = F->params->coltot[Gid^GF];
399
 
    ncols = virtpi[Ga];
400
 
    nlinks = virtpi[Gd];
401
 
 
402
 
    if(nrows && ncols && nlinks)
403
 
      C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gid][0], nrows,
404
 
              &(C2->matrix[Gkj][kj][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
405
 
 
406
 
    dpd_free_block(F->matrix[Gid], virtpi[Gd], F->params->coltot[Gid^GF]);
407
 
 
408
 
    /* -t_ikad * F_jdcb */
409
 
    Gjd = Gj ^ Gd;
410
 
    Gcb = Gjd ^ GF;
411
 
    Ga = Gik ^ Gd ^ GC;     
412
 
 
413
 
    ad = C2->col_offset[Gik][Ga];
414
 
    jd = F->row_offset[Gjd][J];
415
 
 
416
 
    F->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gjd^GF]);
417
 
    dpd_buf4_mat_irrep_rd_block(F, Gjd, jd, virtpi[Gd]);
418
 
 
419
 
    nrows = F->params->coltot[Gjd^GF];
420
 
    ncols = virtpi[Ga];
421
 
    nlinks = virtpi[Gd];
422
 
 
423
 
    if(nrows && ncols && nlinks)
424
 
      C_DGEMM('t','t',nrows, ncols, nlinks, -1.0, F->matrix[Gjd][0], nrows,
425
 
              &(C2->matrix[Gik][ik][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
426
 
 
427
 
    dpd_free_block(F->matrix[Gjd], virtpi[Gd], F->params->coltot[Gjd^GF]);
428
 
 
429
 
    /* +t_ijad * F_kdcb */
430
 
    Gkd = Gk ^ Gd;
431
 
    Gcb = Gkd ^ GF;
432
 
    Ga = Gij ^ Gd ^ GC;     
433
 
 
434
 
    ad = C2->col_offset[Gij][Ga];
435
 
    kd = F->row_offset[Gkd][K];
436
 
 
437
 
    F->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], F->params->coltot[Gkd^GF]);
438
 
    dpd_buf4_mat_irrep_rd_block(F, Gkd, kd, virtpi[Gd]);
439
 
 
440
 
    nrows = F->params->coltot[Gkd^GF];
441
 
    ncols = virtpi[Ga];
442
 
    nlinks = virtpi[Gd];
443
 
 
444
 
    if(nrows && ncols && nlinks)
445
 
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gkd][0], nrows,
446
 
              &(C2->matrix[Gij][ij][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
447
 
 
448
 
    dpd_free_block(F->matrix[Gkd], virtpi[Gd], F->params->coltot[Gkd^GF]);
449
 
 
450
 
  }
451
 
 
452
 
  for(Gl=0; Gl < nirreps; Gl++) {
453
 
    /* +t_ilcb * E_jkla */
454
 
    Gil = Gi ^ Gl; 
455
 
    Gcb  = Gil ^ GC; 
456
 
    Ga = Gjk ^ Gl ^ GE;       
457
 
 
458
 
    la = E->col_offset[Gjk][Gl];
459
 
    il = C2->row_offset[Gil][I];
460
 
 
461
 
    nrows = C2->params->coltot[Gil^GC];
462
 
    ncols = virtpi[Ga];
463
 
    nlinks = occpi[Gl];
464
 
 
465
 
    if(nrows && ncols && nlinks)
466
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, C2->matrix[Gil][il], nrows, 
467
 
              &(E->matrix[Gjk][jk][la]), ncols, 1.0, W2[Gcb][0], ncols);
468
 
 
469
 
    /* -t_jlcb * E_ikla */
470
 
    Gjl = Gj ^ Gl; 
471
 
    Gcb = Gjl ^ GC; 
472
 
    Ga = Gik ^ Gl ^ GE;       
473
 
 
474
 
    la = E->col_offset[Gik][Gl];
475
 
    jl = C2->row_offset[Gjl][J];
476
 
 
477
 
    nrows = C2->params->coltot[Gjl^GC];
478
 
    ncols = virtpi[Ga];
479
 
    nlinks = occpi[Gl];
480
 
 
481
 
    if(nrows && ncols && nlinks)
482
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gjl][jl], nrows, 
483
 
              &(E->matrix[Gik][ik][la]), ncols, 1.0, W2[Gcb][0], ncols);
484
 
 
485
 
    /* -t_klcb * E_jila */
486
 
    Gkl = Gk ^ Gl; 
487
 
    Gcb = Gkl ^ GC; 
488
 
    Ga = Gji ^ Gl ^ GE;      
489
 
 
490
 
    la = E->col_offset[Gji][Gl];
491
 
    kl = C2->row_offset[Gkl][K];
492
 
 
493
 
    nrows = C2->params->coltot[Gkl^GC];
494
 
    ncols = virtpi[Ga];
495
 
    nlinks = occpi[Gl];
496
 
 
497
 
    if(nrows && ncols && nlinks)
498
 
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, C2->matrix[Gkl][kl], nrows, 
499
 
              &(E->matrix[Gji][ji][la]), ncols, 1.0, W2[Gcb][0], ncols);
500
 
 
501
 
  }
502
 
 
503
 
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
504
 
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
505
 
              vir_off, virtpi, vir_off, F->params->colidx, cba, 1);
506
 
 
507
 
  for(Gab=0; Gab < nirreps; Gab++) {
508
 
    Gc = Gab ^ Gijk ^ GX3; /* assumes totally symmetric! */
509
 
 
510
 
    for(ab=0; ab < F->params->coltot[Gab]; ab++) {
511
 
      A = F->params->colorb[Gab][ab][0];
512
 
      B = F->params->colorb[Gab][ab][1];
513
 
      Ga = F->params->rsym[A];
514
 
      Gb = F->params->ssym[B];
515
 
      a = A - vir_off[Ga];
516
 
      b = B - vir_off[Gb];
517
 
 
518
 
      for(c=0; c < virtpi[Gc]; c++) {
519
 
        C = vir_off[Gc] + c;
520
 
 
521
 
        denom = dijk;
522
 
        if(fAB->params->rowtot[Ga]) denom -= fAB->matrix[Ga][a][a];
523
 
        if(fAB->params->rowtot[Gb]) denom -= fAB->matrix[Gb][b][b];
524
 
        if(fAB->params->rowtot[Gc]) denom -= fAB->matrix[Gc][c][c];
525
 
 
526
 
        W1[Gab][ab][c] /= (omega + denom);
527
 
 
528
 
      } /* c */
529
 
    } /* ab */
530
 
  } /* Gab */
531
 
 
532
 
  for(Gab=0; Gab < nirreps; Gab++) {
533
 
    Gc = Gab ^ Gijk ^ GX3; /* changed */
534
 
    dpd_free_block(W2[Gab], F->params->coltot[Gab], virtpi[Gc]);
535
 
  }
536
 
  free(W2);
537
 
 
538
 
  dpd_file2_mat_close(fIJ);
539
 
  dpd_file2_mat_close(fAB);
540
 
 
541
 
  for(h=0; h < nirreps; h++) {
542
 
    dpd_buf4_mat_irrep_close(C2, h);
543
 
    dpd_buf4_mat_irrep_close(E, h);
544
 
  }
545
 
}