← Back to branch summary

~ubuntu-branches/ubuntu/trusty/psicode/trusty

« back to all changes in this revision

Viewing changes to src/lib/libdpd/T3_RHF_ic.cc

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

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

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libdpd/T3_RHF_ic.cc
 
1
/*! \file
 
2
    \ingroup DPD
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
/* T3_RHF_ic(): Computes all T3(IJK,ABC) amplitudes for a given I, J,
 
7
** K combination for input T2, F, and E intermediates.  This function
 
8
** is specifically for AAB spin cases with RHF orbitals.
 
9
**  This _ic version assumes WAbEi is available in memory
 
10
** Arguments: 
 
11
**
 
12
**   double ***W1: The target triples amplitudes in an nirreps x AB x
 
13
**   C array.  The memory for this must be allocated externally.
 
14
**
 
15
**   int nirreps: Number of irreps.
 
16
**
 
17
**   int I: Absolute index of orbital I.
 
18
**
 
19
**   int Gi: Irrep of I.
 
20
**
 
21
**   int J: Absolute index of orbital J.
 
22
**
 
23
**   int Gj: Irrep of J.
 
24
**
 
25
**   int K: Absolute index of orbital K.
 
26
**
 
27
**   int Gk: Irrep of K.
 
28
**
 
29
**   dpdbuf4 *T2: Pointer to dpd buffer for double excitation amps,
 
30
**   ordered (IJ,AB).
 
31
**
 
32
**   dpdbuf4 *F: Pointer to dpd buffer for three-virtual-index
 
33
**   intermediate, ordered (IA,BC).
 
34
**
 
35
**   dpdbuf4 *E: Pointer to dpd buffer for three-occupied-index
 
36
**   intermediate, ordered (IJ,KA).
 
37
**
 
38
**   dpdfile2 *fIJ: Pointer to the dpd file2 for the occ-occ block of
 
39
**   the Fock matrix (or other appropriate one-electron operator).
 
40
**   
 
41
**   dpdfile2 *fAB: Pointer to the dpd file2 for the vir-vir block of
 
42
**   the Fock matrix (or other appropriate one-electron operator).
 
43
**   
 
44
**   int *occpi: Number of occupied orbitals per irrep lookup array.
 
45
**
 
46
**   int *occ_off: Offset lookup for translating between absolute and
 
47
**   relative orbital indices for occupied space.
 
48
**
 
49
**   int *virtpi: Number of virtual orbitals per irrep lookup array.
 
50
**
 
51
**   int *vir_off: Offset lookup for translating between absolute and
 
52
**   relative orbital indices for virtual space.
 
53
**
 
54
**   double omega: constant to add to denominators - needed for EOM CC3
 
55
**
 
56
** TDC, July 2004
 
57
** -modified for RHF, RAK 2004
 
58
** -omega argument added, RAK 2006
 
59
*/
 
60
 
 
61
#include <cstdio>
 
62
#include <cstdlib>
 
63
#include <cstring>
 
64
#include <cmath>
 
65
#include <libqt/qt.h>
 
66
#include <libdpd/dpd.h>
 
67
#include <ccfiles.h>
 
68
#include <pthread.h>
 
69
 
 
70
extern "C" {
 
71
 
 
72
void T3_RHF_ic(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk, 
 
73
                dpdbuf4 *T2, dpdbuf4 *F, dpdbuf4 *E, dpdfile2 *fIJ, dpdfile2 *fAB, 
 
74
                int *occpi, int *occ_off, int *virtpi, int *vir_off, double omega)
 
75
{
 
76
  int h;
 
77
  int i, j, k;
 
78
  int ij, ji, ik, ki, jk, kj;
 
79
  int Gij, Gji, Gik, Gki, Gjk, Gkj, Gijk;
 
80
  int Ga, Gb, Gc;
 
81
  int Gd, Gl;
 
82
  int Gid, Gjd, Gkd;
 
83
  int Gab, Gcb, Gca;
 
84
  int Gla, Glb, Glc;
 
85
  int Gil, Gjl, Gkl;
 
86
  int a, b, c, A, B, C;
 
87
  int ab;
 
88
  int cd, bd, ad,dc;
 
89
  int id, jd, kd;
 
90
  int la, lb, lc;
 
91
  int il, jl, kl;
 
92
  int nrows, ncols, nlinks;
 
93
  double dijk, denom;
 
94
  double ***W2;
 
95
  int GE, GF, GC, GX3;
 
96
 
 
97
  GC = T2->file.my_irrep;
 
98
  /* F and E are assumed to have same irrep */
 
99
  GF = GE =  F->file.my_irrep;
 
100
  GX3 = GC^GF;
 
101
 
 
102
  i = I - occ_off[Gi];
 
103
  j = J - occ_off[Gj];
 
104
  k = K - occ_off[Gk];
 
105
 
 
106
  Gij = Gji = Gi ^ Gj;
 
107
  Gik = Gki = Gi ^ Gk;
 
108
  Gjk = Gkj = Gj ^ Gk;
 
109
  Gijk = Gi ^ Gj ^ Gk;
 
110
 
 
111
  ij = T2->params->rowidx[I][J];
 
112
  ji = T2->params->rowidx[J][I];
 
113
  jk = T2->params->rowidx[J][K];
 
114
  kj = T2->params->rowidx[K][J];
 
115
  ik = T2->params->rowidx[I][K];
 
116
  ki = T2->params->rowidx[K][I];
 
117
 
 
118
  dijk = 0.0;
 
119
  if(fIJ->params->rowtot[Gi]) dijk += fIJ->matrix[Gi][i][i];
 
120
  if(fIJ->params->rowtot[Gj]) dijk += fIJ->matrix[Gj][j][j];
 
121
  if(fIJ->params->rowtot[Gk]) dijk += fIJ->matrix[Gk][k][k];
 
122
 
 
123
  W2 = (double ***) malloc(nirreps * sizeof(double **));
 
124
 
 
125
  for(Gab=0; Gab < nirreps; Gab++) {
 
126
    Gc = Gab ^ Gijk ^ GX3;
 
127
 
 
128
    W2[Gab] = dpd_block_matrix(F->params->coltot[Gab], virtpi[Gc]);
 
129
 
 
130
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
131
      memset(W1[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
132
    }
 
133
  }
 
134
 
 
135
  /***** do terms that are (ab,c) and don't need sorted *****/
 
136
  for(Gd=0; Gd < nirreps; Gd++) {
 
137
 
 
138
    /* term 1F */
 
139
    /* +t_JkDc * F_IDAB */
 
140
    Gid = Gi ^ Gd;
 
141
    Gab = Gid ^ GF;
 
142
    Gc = Gjk ^ Gd ^ GC;
 
143
 
 
144
    cd = T2->col_offset[Gjk][Gc];
 
145
    id = F->row_offset[Gid][I];
 
146
 
 
147
    nrows = F->params->coltot[Gid^GF];
 
148
    ncols = virtpi[Gc];
 
149
    nlinks = virtpi[Gd];
 
150
 
 
151
    if(nrows && ncols && nlinks)
 
152
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][id], nrows,
 
153
              &(T2->matrix[Gjk][kj][cd]), nlinks, 1.0, W1[Gab][0], ncols);
 
154
 
 
155
  }
 
156
 
 
157
  for(Gl=0; Gl < nirreps; Gl++) {
 
158
 
 
159
    /* term 1E */
 
160
    /* -t_ILAB * E_JkLc */
 
161
    Gil = Gi ^ Gl;
 
162
    Gab = Gil ^ GC ; 
 
163
    Gc = Gjk ^ Gl ^ GE;
 
164
 
 
165
    lc = E->col_offset[Gjk][Gl];
 
166
    il = T2->row_offset[Gil][I];
 
167
 
 
168
    nrows = T2->params->coltot[Gil^GC];
 
169
    ncols = virtpi[Gc];
 
170
    nlinks = occpi[Gl];
 
171
 
 
172
    if(nrows && ncols && nlinks)
 
173
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, T2->matrix[Gil][il], nrows, 
 
174
              &(E->matrix[Gjk][jk][lc]), ncols, 1.0, W1[Gab][0], ncols);
 
175
 
 
176
  }
 
177
 
 
178
  /***** Do terms (ba,c) that need bac sorts *****/
 
179
 
 
180
  /* zero out W2 memory */
 
181
  for(Gab=0; Gab < nirreps; Gab++) {
 
182
    Gc = Gab ^ Gijk ^ GX3;
 
183
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
184
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
185
    }
 
186
  }
 
187
 
 
188
  for(Gd=0; Gd < nirreps; Gd++) {
 
189
    /* term 6F */
 
190
    /* + F_JdBa * t_IkDc */
 
191
    Gjd = Gj ^ Gd;
 
192
    Gab = Gjd ^ GF;
 
193
    Gc = Gik ^ Gd ^ GC;
 
194
 
 
195
    cd = T2->col_offset[Gik][Gc];
 
196
    jd = F->row_offset[Gjd][J];
 
197
 
 
198
    nrows = F->params->coltot[Gjd^GF];
 
199
    ncols = virtpi[Gc];
 
200
    nlinks = virtpi[Gd];
 
201
 
 
202
    if(nrows && ncols && nlinks)
 
203
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][jd], nrows,
 
204
              &(T2->matrix[Gik][ki][cd]), nlinks, 1.0, W2[Gab][0], ncols);
 
205
 
 
206
  }
 
207
 
 
208
  for(Gl=0; Gl < nirreps; Gl++) {
 
209
    /* term 6E */
 
210
    /* +t_JlAb * E_IkLc */
 
211
    Gjl = Gj ^ Gl;
 
212
    Gab = Gjl ^ GC;
 
213
    Gc = Gik ^ Gl ^ GE;
 
214
 
 
215
    lc = E->col_offset[Gik][Gl];
 
216
    jl = T2->row_offset[Gjl][J];
 
217
 
 
218
    nrows = T2->params->coltot[Gjl^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, T2->matrix[Gjl][jl], nrows,
 
224
              &(E->matrix[Gik][ik][lc]), ncols, 1.0, W2[Gab][0], ncols);
 
225
  }
 
226
 
 
227
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
228
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off,
 
229
              vir_off, virtpi, vir_off, F->params->colidx, bac, 1);
 
230
 
 
231
 
 
232
  /***** do (ac,b) terms that require acb sort *****/
 
233
 
 
234
  /* zero out W2 memory */
 
235
  for(Gab=0; Gab < nirreps; Gab++) {
 
236
    Gc = Gab ^ Gijk ^ GX3;
 
237
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
238
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
239
    }
 
240
  }
 
241
 
 
242
  for(Gd=0; Gd < nirreps; Gd++) {
 
243
    /* term 2F */
 
244
    /* + F_IdAc * t_JkBd */
 
245
    Gid = Gi ^ Gd;
 
246
    Gca = Gid ^ GF;
 
247
    Gb = Gjk ^ Gd ^ GC;
 
248
 
 
249
    bd = T2->col_offset[Gjk][Gb];
 
250
    id = F->row_offset[Gid][I];
 
251
 
 
252
    nrows = F->params->coltot[Gid^GF];
 
253
    ncols = virtpi[Gb];
 
254
    nlinks = virtpi[Gd];
 
255
 
 
256
    if(nrows && ncols && nlinks)
 
257
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gid][id], nrows,
 
258
              &(T2->matrix[Gjk][jk][bd]), nlinks, 1.0, W2[Gca][0], ncols);
 
259
 
 
260
  }
 
261
 
 
262
  for(Gl=0; Gl < nirreps; Gl++) {
 
263
    /* term 2E */
 
264
    /* -t_IlAc * E_KjLb */
 
265
    Gil = Gi ^ Gl;
 
266
    Gca = Gil ^ GC;
 
267
    Gb = Gjk ^ Gl ^ GE;
 
268
 
 
269
    lb = E->col_offset[Gjk][Gl];
 
270
    il = T2->row_offset[Gil][I];
 
271
 
 
272
    nrows = T2->params->coltot[Gil^GC];
 
273
    ncols = virtpi[Gb];
 
274
    nlinks = occpi[Gl];
 
275
 
 
276
    if(nrows && ncols && nlinks)
 
277
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, T2->matrix[Gil][il], nrows, 
 
278
              &(E->matrix[Gjk][kj][lb]), ncols, 1.0, W2[Gca][0], ncols);
 
279
  }
 
280
 
 
281
 
 
282
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
283
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off,
 
284
              vir_off, virtpi, vir_off, F->params->colidx, acb, 1);
 
285
 
 
286
 
 
287
  /***** do (ca,b) terms that need a bca sort *****/
 
288
 
 
289
  /* zero out W2 memory */
 
290
  for(Gab=0; Gab < nirreps; Gab++) {
 
291
    Gc = Gab ^ Gijk ^ GX3;
 
292
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
293
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
294
    }
 
295
  }
 
296
 
 
297
  for(Gd=0; Gd < nirreps; Gd++) {
 
298
    /* term 3F */
 
299
    /* + F_KdCa * t_JiBd */
 
300
    Gkd = Gk ^ Gd;
 
301
    Gca = Gkd ^ GF;
 
302
    Gb = Gij ^ Gd ^ GC;
 
303
 
 
304
    bd = T2->col_offset[Gij][Gb];
 
305
    kd = F->row_offset[Gkd][K];
 
306
 
 
307
    nrows = F->params->coltot[Gkd^GF];
 
308
    ncols = virtpi[Gb];
 
309
    nlinks = virtpi[Gd];
 
310
 
 
311
    if(nrows && ncols && nlinks)
 
312
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gkd][kd], nrows,
 
313
              &(T2->matrix[Gij][ji][bd]), nlinks, 1.0, W2[Gca][0], ncols);
 
314
 
 
315
  }
 
316
 
 
317
  for(Gl=0; Gl < nirreps; Gl++) {
 
318
    /* term 3E */
 
319
    /* - t_KlCa * E_IjLb */
 
320
    Gkl = Gk ^ Gl;
 
321
    Gca = Gkl ^ GC ;
 
322
    Gb = Gji ^ Gl ^ GE;
 
323
 
 
324
    lb = E->col_offset[Gji][Gl];
 
325
    kl = T2->row_offset[Gkl][K];
 
326
 
 
327
    nrows = T2->params->coltot[Gkl^GC];
 
328
    ncols = virtpi[Gb];
 
329
    nlinks = occpi[Gl];
 
330
 
 
331
    if(nrows && ncols && nlinks)
 
332
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, T2->matrix[Gkl][kl], nrows, 
 
333
              &(E->matrix[Gji][ij][lb]), ncols, 1.0, W2[Gca][0], ncols);
 
334
  }
 
335
 
 
336
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
337
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off,
 
338
              vir_off, virtpi, vir_off, F->params->colidx, bca, 1);
 
339
 
 
340
 
 
341
  /***** do (cb,a) terms that require a cba sort *****/
 
342
  /* zero out W2 memory */
 
343
  for(Gab=0; Gab < nirreps; Gab++) {
 
344
    Gc = Gab ^ Gijk ^ GX3;
 
345
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
346
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
347
    }
 
348
  }
 
349
 
 
350
  for(Gd=0; Gd < nirreps; Gd++) {
 
351
    /* term 4F */
 
352
    /* + F_KdCb * t_IjAd */
 
353
    Gkd = Gk ^ Gd;
 
354
    Gcb = Gkd ^ GF;
 
355
    Ga = Gij ^ Gd ^ GC;
 
356
 
 
357
    ad = T2->col_offset[Gij][Ga];
 
358
    kd = F->row_offset[Gkd][K];
 
359
 
 
360
    nrows = F->params->coltot[Gkd^GF];
 
361
    ncols = virtpi[Ga];
 
362
    nlinks = virtpi[Gd];
 
363
 
 
364
    if(nrows && ncols && nlinks)
 
365
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gkd][kd], nrows,
 
366
              &(T2->matrix[Gij][ij][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
 
367
 
 
368
  }
 
369
 
 
370
  for(Gl=0; Gl < nirreps; Gl++) {
 
371
    /* term 4E */
 
372
    /*  -t_KlCb * E_JiLa */
 
373
    Gkl = Gk ^ Gl;
 
374
    Gcb  = Gkl ^ GC;
 
375
    Ga = Gji ^ Gl ^ GE;
 
376
 
 
377
    la = E->col_offset[Gji][Gl];
 
378
    kl = T2->row_offset[Gkl][K];
 
379
 
 
380
    nrows = T2->params->coltot[Gkl^GC];
 
381
    ncols = virtpi[Ga];
 
382
    nlinks = occpi[Gl];
 
383
 
 
384
    if(nrows && ncols && nlinks)
 
385
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, T2->matrix[Gkl][kl], nrows, 
 
386
              &(E->matrix[Gji][ji][la]), ncols, 1.0, W2[Gcb][0], ncols);
 
387
  }
 
388
 
 
389
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
390
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
 
391
              vir_off, virtpi, vir_off, F->params->colidx, cba, 1);
 
392
 
 
393
  /***** do (bc,a) terms that require a cab sort *****/
 
394
  /* zero out W2 memory */
 
395
  for(Gab=0; Gab < nirreps; Gab++) {
 
396
    Gc = Gab ^ Gijk ^ GX3;
 
397
    if(F->params->coltot[Gab] && virtpi[Gc]) {
 
398
      memset(W2[Gab][0], 0, F->params->coltot[Gab]*virtpi[Gc]*sizeof(double));
 
399
    }
 
400
  }
 
401
 
 
402
  for(Gd=0; Gd < nirreps; Gd++) {
 
403
    /* term 5F */
 
404
    /* + t_IkAd * F_JdBc */
 
405
    Gjd = Gj ^ Gd;
 
406
    Gcb = Gjd ^ GF;
 
407
    Ga = Gik ^ Gd ^ GC;
 
408
 
 
409
    ad = T2->col_offset[Gik][Ga];
 
410
    jd = F->row_offset[Gjd][J];
 
411
 
 
412
    nrows = F->params->coltot[Gjd^GF];
 
413
    ncols = virtpi[Ga];
 
414
    nlinks = virtpi[Gd];
 
415
 
 
416
    if(nrows && ncols && nlinks)
 
417
      C_DGEMM('t','t',nrows, ncols, nlinks, 1.0, F->matrix[Gjd][jd], nrows,
 
418
              &(T2->matrix[Gik][ik][ad]), nlinks, 1.0, W2[Gcb][0], ncols);
 
419
 
 
420
  }
 
421
 
 
422
  for(Gl=0; Gl < nirreps; Gl++) {
 
423
    /* term 5E */
 
424
    /* - t_JlBc * E_Kila */
 
425
    Gjl = Gj ^ Gl;
 
426
    Gcb  = Gjl ^ GC;
 
427
    Ga = Gik ^ Gl ^ GE;
 
428
 
 
429
    la = E->col_offset[Gik][Gl];
 
430
    jl = T2->row_offset[Gjl][J];
 
431
 
 
432
    nrows = T2->params->coltot[Gjl^GC];
 
433
    ncols = virtpi[Ga];
 
434
    nlinks = occpi[Gl];
 
435
 
 
436
    if(nrows && ncols && nlinks)
 
437
      C_DGEMM('t', 'n', nrows, ncols, nlinks, -1.0, T2->matrix[Gjl][jl], nrows, 
 
438
              &(E->matrix[Gik][ki][la]), ncols, 1.0, W2[Gcb][0], ncols);
 
439
  }
 
440
 
 
441
  dpd_3d_sort(W2, W1, nirreps, Gijk^GX3, F->params->coltot, F->params->colidx,
 
442
              F->params->colorb, F->params->rsym, F->params->ssym, vir_off, 
 
443
              vir_off, virtpi, vir_off, F->params->colidx, cab, 1);
 
444
 
 
445
  for(Gab=0; Gab < nirreps; Gab++) {
 
446
    Gc = Gab ^ Gijk ^ GX3;
 
447
 
 
448
    for(ab=0; ab < F->params->coltot[Gab]; ab++) {
 
449
      A = F->params->colorb[Gab][ab][0];
 
450
      B = F->params->colorb[Gab][ab][1];
 
451
      Ga = F->params->rsym[A];
 
452
      Gb = F->params->ssym[B];
 
453
      a = A - vir_off[Ga];
 
454
      b = B - vir_off[Gb];
 
455
 
 
456
      for(c=0; c < virtpi[Gc]; c++) {
 
457
        C = vir_off[Gc] + c;
 
458
 
 
459
        denom = dijk;
 
460
        if(fAB->params->rowtot[Ga]) denom -= fAB->matrix[Ga][a][a];
 
461
        if(fAB->params->rowtot[Gb]) denom -= fAB->matrix[Gb][b][b];
 
462
        if(fAB->params->rowtot[Gc]) denom -= fAB->matrix[Gc][c][c];
 
463
 
 
464
        W1[Gab][ab][c] /= (denom + omega);
 
465
 
 
466
      } /* c */
 
467
    } /* ab */
 
468
  } /* Gab */
 
469
 
 
470
  for(Gab=0; Gab < nirreps; Gab++) {
 
471
    Gc = Gab ^ Gijk ^ GX3;
 
472
    dpd_free_block(W2[Gab], F->params->coltot[Gab], virtpi[Gc]);
 
473
  }
 
474
  free(W2);
 
475
}
 
476
 
 
477
} /* extern "C" */