← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/cc3_Wmnie.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/bin/ccenergy/cc3_Wmnie.c
1
 
#include <libdpd/dpd.h>
2
 
#define EXTERN
3
 
#include "globals.h"
4
 
 
5
 
/* cc3_Wmnie(): Compute the Wmnie matrix from CC3 theory, which is
6
 
** given in spin-orbitals as:
7
 
** 
8
 
** Wmnie = <mn||ie> + t_i^f <mn||fe>
9
 
**
10
 
** TDC, Feb 2004
11
 
*/
12
 
 
13
 
void purge_Wmnie(void);
14
 
 
15
 
void cc3_Wmnie(void)
16
 
{
17
 
  dpdbuf4 E, D, W, Z;
18
 
  dpdfile2 tIA, tia;
19
 
 
20
 
  if(params.ref == 0) { /** RHF **/
21
 
 
22
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
23
 
    dpd_buf4_copy(&E, CC3_HET1, "CC3 WMnIe (Mn,Ie)");
24
 
    dpd_buf4_close(&E);
25
 
 
26
 
    dpd_buf4_init(&W, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMnIe (Mn,Ie)");
27
 
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
28
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
29
 
    dpd_contract244(&tIA, &D, &W, 1, 2, 1, 1, 1);
30
 
    dpd_file2_close(&tIA);
31
 
    dpd_buf4_close(&D);
32
 
    dpd_buf4_close(&W);
33
 
  }
34
 
 
35
 
  else if (params.ref == 1) { /* ROHF */
36
 
 
37
 
    /** W(M>N,IE) <--- <MN||IE> **/
38
 
    /** W(m>n,ie) <--- <mn||ie> **/
39
 
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
40
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 2, 11, "CC3 WMNIE (M>N,EI)");
41
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 2, 11, "CC3 Wmnie (m>n,ei)");
42
 
    dpd_buf4_close(&E);
43
 
 
44
 
    /** W(Mn,Ie) <--- <Mn|Ie> **/
45
 
    /** W(mN,iE) <--- <mN|iE> **/
46
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
47
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 0, 11, "CC3 WMnIe (Mn,eI)");
48
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 0, 11, "CC3 WmNiE (mN,Ei)");
49
 
    dpd_buf4_close(&E);
50
 
 
51
 
    /**** Term 2 ****/
52
 
 
53
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
54
 
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
55
 
 
56
 
    /* <M>N||EF> T(I,F) --> W(M>N,EI) */
57
 
    /* <m>n||ef> T(i,f) --> W(m>n,ei) */
58
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 11, 2, 11, 0, "CC3 WMNIE (M>N,EI)");
59
 
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 5, 0, "D <ij||ab> (i>j,ab)");
60
 
    dpd_contract424(&D, &tIA, &W, 3, 1, 0, -1, 1.0);
61
 
    dpd_buf4_close(&W);
62
 
 
63
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 11, 2, 11, 0, "CC3 Wmnie (m>n,ei)");
64
 
    dpd_contract424(&D, &tia, &W, 3, 1, 0, -1, 1.0);
65
 
    dpd_buf4_close(&D);
66
 
    dpd_buf4_close(&W);
67
 
 
68
 
    /* Z(nM,eI) = <nM|eF> T(I,F) */
69
 
    dpd_buf4_init(&Z, CC_TMP1, 0, 0, 11, 0, 11, 0, "Z(nM,eI)");
70
 
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
71
 
    dpd_contract424(&D, &tIA, &Z, 3, 1, 0, 1, 0);
72
 
    dpd_buf4_close(&D);
73
 
    /* Z(nM,eI) --> W(Mn,eI) */
74
 
    dpd_buf4_sort_axpy(&Z, CC3_HET1, qprs, 0, 11, "CC3 WMnIe (Mn,eI)", 1);
75
 
    dpd_buf4_close(&Z);
76
 
 
77
 
    /* Z(Nm,Ei) = <Nm|Ef> T(i,f) */
78
 
    dpd_buf4_init(&Z, CC_TMP1, 0, 0, 11, 0, 11, 0, "Z(Nm,Ei)");
79
 
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
80
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
81
 
    dpd_buf4_close(&D);
82
 
    /* Z(Nm,Ei) --> W(mN,Ei) */
83
 
    dpd_buf4_sort_axpy(&Z, CC3_HET1, qprs, 0, 11, "CC3 WmNiE (mN,Ei)", 1);
84
 
    dpd_buf4_close(&Z);
85
 
 
86
 
    /* purge (mn,ei)'s before sorting */
87
 
    purge_Wmnie();
88
 
 
89
 
    /* also put "normal" sorted versions in CC3_HET1 */
90
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 11, 2, 11, 0, "CC3 WMNIE (M>N,EI)");
91
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 2, 10, "CC3 WMNIE (M>N,IE)");
92
 
    dpd_buf4_close(&W);
93
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 11, 2, 11, 0, "CC3 Wmnie (m>n,ei)");
94
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 2, 10, "CC3 Wmnie (m>n,ie)");
95
 
    dpd_buf4_close(&W);
96
 
    dpd_buf4_init(&W, CC3_HET1, 0, 0, 11, 0, 11, 0, "CC3 WMnIe (Mn,eI)");
97
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 0, 10, "CC3 WMnIe (Mn,Ie)");
98
 
    dpd_buf4_close(&W);
99
 
    dpd_buf4_init(&W, CC3_HET1, 0, 0, 11, 0, 11, 0, "CC3 WmNiE (mN,Ei)");
100
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 0, 10, "CC3 WmNiE (mN,iE)");
101
 
    dpd_buf4_close(&W);
102
 
 
103
 
    dpd_file2_close(&tIA);
104
 
    dpd_file2_close(&tia);
105
 
 
106
 
  }
107
 
 
108
 
  else if (params.ref == 2) {
109
 
 
110
 
    /** W(M>N,IE) <--- <MN||IE> **/
111
 
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
112
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 2, 21, "CC3 WMNIE (M>N,EI)");
113
 
    dpd_buf4_close(&E);
114
 
 
115
 
    /** W(m>n,ie) <--- <mn||ie> **/
116
 
    dpd_buf4_init(&E, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
117
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 12, 31, "CC3 Wmnie (m>n,ei)");
118
 
    dpd_buf4_close(&E);
119
 
 
120
 
    /** W(Mn,Ie) <--- <Mn|Ie> **/
121
 
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
122
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 22, 25, "CC3 WMnIe (Mn,eI)");
123
 
    dpd_buf4_close(&E);
124
 
 
125
 
    /** W(mN,iE) <--- <mN|iE> **/
126
 
    dpd_buf4_init(&E, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
127
 
    dpd_buf4_sort(&E, CC3_HET1, pqsr, 23, 26, "CC3 WmNiE (mN,Ei)");
128
 
    dpd_buf4_close(&E);
129
 
 
130
 
    /**** Term 2 ****/
131
 
 
132
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
133
 
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
134
 
 
135
 
    /* <M>N||EF> T(I,F) --> W(M>N,EI) */
136
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 21, 2, 21, 0, "CC3 WMNIE (M>N,EI)");
137
 
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 5, 0, "D <IJ||AB> (I>J,AB)");
138
 
    dpd_contract424(&D, &tIA, &W, 3, 1, 0, -1, 1.0);
139
 
    dpd_buf4_close(&D);
140
 
    dpd_buf4_close(&W);
141
 
 
142
 
    /* <m>n||ef> T(i,f) --> W(m>n,ei) */
143
 
    dpd_buf4_init(&W, CC3_HET1, 0, 12, 31, 12, 31, 0, "CC3 Wmnie (m>n,ei)");
144
 
    dpd_buf4_init(&D, CC_DINTS, 0, 12, 15, 12, 15, 0, "D <ij||ab> (i>j,ab)");
145
 
    dpd_contract424(&D, &tia, &W, 3, 1, 0, -1, 1.0);
146
 
    dpd_buf4_close(&D);
147
 
    dpd_buf4_close(&W);
148
 
 
149
 
    /* Z(nM,eI) = <nM|eF> T(I,F) */
150
 
    dpd_buf4_init(&Z, CC_TMP1, 0, 23, 25, 23, 25, 0, "Z(nM,eI)");
151
 
    dpd_buf4_init(&D, CC_DINTS, 0, 23, 29, 23, 29, 0, "D <iJ|aB>");
152
 
    dpd_contract424(&D, &tIA, &Z, 3, 1, 0, 1, 0);
153
 
    dpd_buf4_close(&D);
154
 
    /* Z(nM,eI) --> W(Mn,eI) */
155
 
    dpd_buf4_sort_axpy(&Z, CC3_HET1, qprs, 22, 25, "CC3 WMnIe (Mn,eI)", 1);
156
 
    dpd_buf4_close(&Z);
157
 
 
158
 
    /* Z(Nm,Ei) = <Nm|Ef> T(i,f) */
159
 
    dpd_buf4_init(&Z, CC_TMP1, 0, 22, 26, 22, 26, 0, "Z(Nm,Ei)");
160
 
    dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
161
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
162
 
    dpd_buf4_close(&D);
163
 
    /* Z(Nm,Ei) --> W(mN,Ei) */
164
 
    dpd_buf4_sort_axpy(&Z, CC3_HET1, qprs, 23, 26, "CC3 WmNiE (mN,Ei)", 1);
165
 
    dpd_buf4_close(&Z);
166
 
 
167
 
    /* also put "normal" sorted versions in CC3_HET1 */
168
 
    dpd_buf4_init(&W, CC3_HET1, 0, 2, 21, 2, 21, 0, "CC3 WMNIE (M>N,EI)");
169
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 2, 20, "CC3 WMNIE (M>N,IE)");
170
 
    dpd_buf4_close(&W);
171
 
    dpd_buf4_init(&W, CC3_HET1, 0, 12, 31, 12, 31, 0, "CC3 Wmnie (m>n,ei)");
172
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 12, 30, "CC3 Wmnie (m>n,ie)");
173
 
    dpd_buf4_close(&W);
174
 
    dpd_buf4_init(&W, CC3_HET1, 0, 22, 25, 22, 25, 0, "CC3 WMnIe (Mn,eI)");
175
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 22, 24, "CC3 WMnIe (Mn,Ie)");
176
 
    dpd_buf4_close(&W);
177
 
    dpd_buf4_init(&W, CC3_HET1, 0, 23, 26, 23, 26, 0, "CC3 WmNiE (mN,Ei)");
178
 
    dpd_buf4_sort(&W, CC3_HET1, pqsr, 23, 27, "CC3 WmNiE (mN,iE)");
179
 
    dpd_buf4_close(&W);
180
 
 
181
 
    dpd_file2_close(&tIA);
182
 
    dpd_file2_close(&tia);
183
 
 
184
 
  }
185
 
}
186
 
 
187
 
 
188
 
/* Purge Wmnie matrix elements */
189
 
void purge_Wmnie(void) {
190
 
  dpdfile4 W;
191
 
  int *occpi, *virtpi;
192
 
  int h, a, b, e, f, i, j, m, n;
193
 
  int    A, B, E, F, I, J, M, N;
194
 
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
195
 
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
196
 
  int *occ_off, *vir_off;
197
 
  int *occ_sym, *vir_sym;
198
 
  int *openpi, nirreps;
199
 
    
200
 
  nirreps = moinfo.nirreps;
201
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
202
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
203
 
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
204
 
  openpi = moinfo.openpi;
205
 
 
206
 
  dpd_file4_init(&W, CC3_HET1, 0, 0, 11,"CC3 WMnIe (Mn,eI)");
207
 
  for(h=0; h < nirreps; h++) {
208
 
    dpd_file4_mat_irrep_init(&W, h);
209
 
    dpd_file4_mat_irrep_rd(&W, h);
210
 
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
211
 
      n = W.params->roworb[h][mn][1];
212
 
      nsym = W.params->qsym[n];
213
 
      N = n - occ_off[nsym];
214
 
      for(ei=0; ei<W.params->coltot[h]; ei++) {
215
 
        if (N >= (occpi[nsym] - openpi[nsym]))
216
 
          W.matrix[h][mn][ei] = 0.0;
217
 
      }
218
 
    }
219
 
    dpd_file4_mat_irrep_wrt(&W, h);
220
 
    dpd_file4_mat_irrep_close(&W, h);
221
 
  }
222
 
 
223
 
  dpd_file4_init(&W, CC3_HET1, 0, 2, 11, "CC3 WMNIE (M>N,EI)");
224
 
  for(h=0; h < W.params->nirreps; h++) {
225
 
    dpd_file4_mat_irrep_init(&W, h);
226
 
    dpd_file4_mat_irrep_rd(&W, h);
227
 
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
228
 
      for(ei=0; ei<W.params->coltot[h]; ei++) {
229
 
        e = W.params->colorb[h][ei][0];
230
 
        esym = W.params->rsym[e];
231
 
        E = e - vir_off[esym];
232
 
        if (E >= (virtpi[esym] - openpi[esym]))
233
 
          W.matrix[h][mn][ei] = 0.0;
234
 
      }
235
 
    }
236
 
    dpd_file4_mat_irrep_wrt(&W, h);
237
 
    dpd_file4_mat_irrep_close(&W, h);
238
 
  }
239
 
  dpd_file4_close(&W);
240
 
  dpd_file4_init(&W, CC3_HET1, 0, 2, 11,"CC3 Wmnie (m>n,ei)");
241
 
  for(h=0; h < nirreps; h++) {
242
 
    dpd_file4_mat_irrep_init(&W, h);
243
 
    dpd_file4_mat_irrep_rd(&W, h);
244
 
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
245
 
      m = W.params->roworb[h][mn][0];
246
 
      n = W.params->roworb[h][mn][1];
247
 
      msym = W.params->psym[m];
248
 
      nsym = W.params->qsym[n];
249
 
      M = m - occ_off[msym];
250
 
      N = n - occ_off[nsym];
251
 
      for(ei=0; ei<W.params->coltot[h]; ei++) {
252
 
        i = W.params->colorb[h][ei][1];
253
 
        isym = W.params->ssym[i];
254
 
        I = i - occ_off[isym];
255
 
        if ((M >= (occpi[msym] - openpi[msym])) ||
256
 
          (N >= (occpi[nsym] - openpi[nsym])) ||
257
 
          (I >= (occpi[isym] - openpi[isym])) )
258
 
          W.matrix[h][mn][ei] = 0.0;
259
 
      }
260
 
    }
261
 
    dpd_file4_mat_irrep_wrt(&W, h);
262
 
    dpd_file4_mat_irrep_close(&W, h);
263
 
  }
264
 
  dpd_file4_close(&W);
265
 
 
266
 
  dpd_file4_init(&W, CC3_HET1, 0, 0, 11,"CC3 WmNiE (mN,Ei)");
267
 
  for(h=0; h < nirreps; h++) {
268
 
    dpd_file4_mat_irrep_init(&W, h);
269
 
    dpd_file4_mat_irrep_rd(&W, h);
270
 
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
271
 
      m = W.params->roworb[h][mn][0];
272
 
      msym = W.params->psym[m];
273
 
      M = m - occ_off[msym];
274
 
      for(ei=0; ei<W.params->coltot[h]; ei++) {
275
 
        e = W.params->colorb[h][ei][0];
276
 
        i = W.params->colorb[h][ei][1];
277
 
        esym = W.params->rsym[e];
278
 
        isym = W.params->ssym[i];
279
 
        E = e - vir_off[esym];
280
 
        I = i - occ_off[isym];
281
 
        if ((M >= (occpi[msym] - openpi[msym])) ||
282
 
            (E >= (virtpi[esym] - openpi[esym])) ||
283
 
            (I >= (occpi[isym] - openpi[isym])) )
284
 
          W.matrix[h][mn][ei] = 0.0;
285
 
      }
286
 
    }
287
 
    dpd_file4_mat_irrep_wrt(&W, h);
288
 
    dpd_file4_mat_irrep_close(&W, h);
289
 
  }
290
 
  dpd_file4_close(&W);
291
 
  return;
292
 
}
293