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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/cc3_t3x.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:
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libdpd/dpd.h>
4
 
#define EXTERN
5
 
#include "globals.h"
6
 
 
7
 
void cc3_t3x(void)
8
 
{
9
 
  if(params.ref == 0) {
10
 
    int h, nirreps;
11
 
    int *occ_off, *occpi;
12
 
    int *vir_off, *virtpi;
13
 
    int Gi, Gj, Gk, Gijk;
14
 
    int Ga, Gb, Gc, Gab;
15
 
    int i, j, k, I, J, K;
16
 
    int a, b, c, A, B, C;
17
 
    int ab;
18
 
    double ***W1;
19
 
    dpdbuf4 T2, E, F, T2AA, T2AB, T2BA, EAA, EAB, EBA, FAA, FAB, FBA;
20
 
    dpdfile2 fIJ, fAB, fij, fab;
21
 
    dpdfile2 XLD;
22
 
    dpdbuf4 L2, L2AB;
23
 
    int Gij, ij, Gbc, bc, Gjk, jk;
24
 
    int nrows, ncols;
25
 
    int **W_offset, offset;
26
 
 
27
 
    nirreps = moinfo.nirreps;
28
 
    occpi = moinfo.occpi;
29
 
    occ_off = moinfo.occ_off;
30
 
    virtpi = moinfo.virtpi;
31
 
    vir_off = moinfo.vir_off;
32
 
 
33
 
    W_offset = init_int_matrix(nirreps, nirreps);
34
 
    for(Gab=0; Gab < nirreps; Gab++) {
35
 
      for(Ga=0,offset=0; Ga < nirreps; Ga++) {
36
 
        Gb = Ga ^ Gab;
37
 
        W_offset[Gab][Ga] = offset;
38
 
        offset += virtpi[Ga] * virtpi[Gb];
39
 
      }
40
 
    }
41
 
 
42
 
    dpd_file2_init(&XLD, CC3_MISC, 0, 0, 1, "CC3 XLD");
43
 
    dpd_file2_mat_init(&XLD);
44
 
 
45
 
    dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
46
 
    dpd_buf4_init(&L2AB, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
47
 
    for(h=0; h < nirreps; h++) {
48
 
      dpd_buf4_mat_irrep_init(&L2, h);
49
 
      dpd_buf4_mat_irrep_rd(&L2, h);
50
 
 
51
 
      dpd_buf4_mat_irrep_init(&L2AB, h);
52
 
      dpd_buf4_mat_irrep_rd(&L2AB, h);
53
 
    }
54
 
 
55
 
    dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
56
 
    dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
57
 
    dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
58
 
    dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
59
 
 
60
 
    dpd_buf4_init(&T2, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
61
 
    dpd_buf4_init(&F, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
62
 
    dpd_buf4_init(&E, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
63
 
 
64
 
    T2AA = T2;
65
 
    dpd_buf4_init(&T2AB, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
66
 
    dpd_buf4_init(&T2BA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tiJaB");
67
 
    FAA = F;
68
 
    dpd_buf4_init(&FAB, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaBeI (Ie,Ba)");
69
 
    dpd_buf4_init(&FBA, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAbEi (iE,bA)");
70
 
    EAA = E;
71
 
    dpd_buf4_init(&EAB, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMbIj (Ij,Mb)");
72
 
    dpd_buf4_init(&EBA, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmBiJ (iJ,mB)");
73
 
 
74
 
    /* target T3 amplitudes go in here */
75
 
    W1 = (double ***) malloc(nirreps * sizeof(double **));
76
 
 
77
 
    for(Gi=0; Gi < nirreps; Gi++) {
78
 
      for(Gj=0; Gj < nirreps; Gj++) {
79
 
        Gij = Gi ^ Gj;
80
 
        for(Gk=0; Gk < nirreps; Gk++) {
81
 
          Gijk = Gi ^ Gj ^ Gk;
82
 
          Gjk = Gj ^ Gk;
83
 
 
84
 
          for(Gab=0; Gab < nirreps; Gab++) {
85
 
            Gc = Gab ^ Gijk; /* totally symmetric */
86
 
            W1[Gab] = dpd_block_matrix(F.params->coltot[Gab], virtpi[Gc]);
87
 
          }
88
 
 
89
 
          for(i=0; i < occpi[Gi]; i++) {
90
 
            I = occ_off[Gi] + i;
91
 
            for(j=0; j < occpi[Gj]; j++) {
92
 
              J = occ_off[Gj] + j;
93
 
              for(k=0; k < occpi[Gk]; k++) {
94
 
                K = occ_off[Gk] + k;
95
 
 
96
 
                T3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, &T2, &F, &E, &fIJ, &fAB, 
97
 
                       occpi, occ_off, virtpi, vir_off, 0.0);
98
 
 
99
 
                /* X_KC <-- 1/4 t_IJKABC <IJ||AB> */
100
 
 
101
 
                Gc = Gk;    /* assumes T1 is totally symmetric */
102
 
                Gab = Gij;  /* assumes <ij||ab> is totally symmetric */
103
 
 
104
 
                ij = L2.params->rowidx[I][J];
105
 
 
106
 
                nrows = L2.params->coltot[Gij];
107
 
                ncols = virtpi[Gc];
108
 
 
109
 
                if(nrows && ncols)
110
 
                  C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, L2.matrix[Gij][ij], 1,
111
 
                          1.0, XLD.matrix[Gk][k], 1);
112
 
 
113
 
                T3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, &T2AA, &T2AB, &T2BA, 
114
 
                       &FAA, &FAB, &FBA, &EAA, &EAB, &EBA, &fIJ, &fij, &fAB, &fab,
115
 
                       occpi, occ_off, occpi, occ_off, virtpi, vir_off, virtpi, vir_off, 0.0);
116
 
 
117
 
                /* t_IA <-- t_IJkABc <Jk|Bc> */
118
 
 
119
 
                Ga = Gi;   /* assumes T1 is totally symmetric */
120
 
                Gbc = Gjk; /* assumes <jk|bc> is totally symmetric */
121
 
 
122
 
                jk = L2AB.params->rowidx[J][K];
123
 
 
124
 
                for(Gab=0; Gab < nirreps; Gab++) {
125
 
                  Gb = Ga ^ Gab;
126
 
                  Gc = Gb ^ Gbc;
127
 
 
128
 
                  ab = W_offset[Gab][Ga];
129
 
                  bc = L2AB.col_offset[Gjk][Gb];
130
 
 
131
 
                  nrows = virtpi[Ga];
132
 
                  ncols = virtpi[Gb] * virtpi[Gc];
133
 
 
134
 
                  if(nrows && ncols)
135
 
                    C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][ab], ncols, &(L2AB.matrix[Gjk][jk][bc]), 1,
136
 
                            1.0, XLD.matrix[Gi][i], 1);
137
 
 
138
 
                }
139
 
                /* t_KC <-- 1/4 t_ijKabC <ij||ab> */
140
 
 
141
 
                Gc = Gk;  /* assumes T1 is totally symmetric */
142
 
                Gab = Gij; /* assumes <ij||ab> is totally symmetric */
143
 
 
144
 
                ij = L2.params->rowidx[I][J];
145
 
 
146
 
                nrows = L2.params->coltot[Gij];
147
 
                ncols = virtpi[Gc];
148
 
 
149
 
                if(nrows && ncols)
150
 
                  C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, L2.matrix[Gij][ij], 1,
151
 
                          1.0, XLD.matrix[Gk][k], 1);
152
 
 
153
 
 
154
 
              } /* k */
155
 
            } /* j */
156
 
          } /* i */
157
 
 
158
 
          for(Gab=0; Gab < nirreps; Gab++) {
159
 
            Gc = Gab ^ Gijk; /* totally symmetric */
160
 
            dpd_free_block(W1[Gab], F.params->coltot[Gab], virtpi[Gc]);
161
 
          }
162
 
        } /* Gk */
163
 
      } /* Gj */
164
 
    } /* Gi */
165
 
 
166
 
    free(W1);
167
 
 
168
 
    dpd_buf4_close(&E);
169
 
    dpd_buf4_close(&F);
170
 
    dpd_buf4_close(&T2);
171
 
    dpd_file2_close(&fIJ);
172
 
    dpd_file2_close(&fAB);
173
 
 
174
 
    dpd_buf4_close(&EAB);
175
 
    dpd_buf4_close(&EBA);
176
 
    dpd_buf4_close(&FAB);
177
 
    dpd_buf4_close(&FBA);
178
 
    dpd_buf4_close(&T2AB);
179
 
    dpd_buf4_close(&T2BA);
180
 
    dpd_file2_close(&fij);
181
 
    dpd_file2_close(&fab);
182
 
 
183
 
    free_int_matrix(W_offset, nirreps);
184
 
 
185
 
    for(h=0; h < nirreps; h++) {
186
 
      dpd_buf4_mat_irrep_close(&L2, h);
187
 
      dpd_buf4_mat_irrep_close(&L2AB, h);
188
 
    }
189
 
    dpd_buf4_close(&L2);
190
 
    dpd_buf4_close(&L2AB);
191
 
 
192
 
    dpd_file2_mat_wrt(&XLD);
193
 
    dpd_file2_close(&XLD);
194
 
  }
195
 
}