~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/cchbar/Wabei_ABAB_UHF.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:
 
1
/*! \file
 
2
    \ingroup CCHBAR
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <libdpd/dpd.h>
 
7
#include "MOInfo.h"
 
8
#include "Params.h"
 
9
#define EXTERN
 
10
#include "globals.h"
 
11
 
 
12
namespace psi { namespace cchbar {
 
13
 
 
14
/* WAbEi_UHF(): Computes all contributions to the AbEi spin case of
 
15
** the Wabei HBAR matrix elements.  The final product is stored in
 
16
** (Ei,Ab) ordering and is referred to on disk as "WAbEi".
 
17
**
 
18
** The spin-orbital expression for the Wabei elements is:
 
19
**
 
20
** Wabei = <ab||ei> - Fme t_mi^ab + t_i^f <ab||ef>
 
21
**         - P(ab) t_m^b <am||ef>t_i^f + 1/2 tau_mn^ab <mn||ef> t_i^f
 
22
**         + 1/2 <mn||ei> tau_mn^ab - P(ab) <mb||ef> t_mi^af
 
23
**         - P(ab) t_m^a { <mb||ei> - t_ni^bf <mn||ef> }
 
24
**
 
25
** (cf. Gauss and Stanton, JCP 103, 3561-3577 (1995).)
 
26
**
 
27
** For the AbEi spin case, we evaluate these contractions with two
 
28
** target orderings, (Ab,Ei) and (Ei,Ab), depending on the term.
 
29
** After all terms have been evaluated, the (Ab,Ei) terms are sorted
 
30
** into (Ei,Ab) ordering and both groups arer added together.
 
31
**
 
32
** TDC, June 2002
 
33
*/
 
34
 
 
35
void WAbEi_UHF(void)
 
36
{
 
37
  dpdfile2 Fme, T1;
 
38
  dpdbuf4 F, W, T2, B, Z, Z1, Z2, D, T, E, C;
 
39
 
 
40
  /**** Term I ****/
 
41
 
 
42
  /** W(Ei,Ab) <--- <Ei|Ab> **/
 
43
  dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
 
44
  dpd_buf4_copy(&F, CC_HBAR, "WEiAb");
 
45
  dpd_buf4_close(&F);
 
46
 
 
47
  /**** Term II ****/
 
48
 
 
49
  /** W(Ei,Ab) <--- - F_ME t_Mi^Ab **/
 
50
  dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
51
  dpd_file2_init(&Fme, CC_OEI, 0, 0, 1, "FME");
 
52
  dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
 
53
  dpd_contract244(&Fme, &T2, &W, 0, 0, 0, -1, 1);
 
54
  dpd_buf4_close(&W);
 
55
  dpd_file2_close(&Fme);
 
56
  dpd_buf4_close(&T2);
 
57
 
 
58
  /**** Term III ****/
 
59
 
 
60
  /** <Ab|Ef> t_i^f **/
 
61
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
62
  dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
63
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
64
  dpd_contract424(&B, &T1, &W, 3, 1, 0, 1, 0);
 
65
  dpd_file2_close(&T1);
 
66
  dpd_buf4_close(&B);
 
67
  dpd_buf4_close(&W);
 
68
 
 
69
  /**** Term IV ****/
 
70
 
 
71
  /** WAbEi <-- - t_m^b <mA|fE> t_i^f - t_M^A <Mb|Ef> t_i^f
 
72
      Evaluate in three steps: 
 
73
          (1) Z_mAEi = - <Am|Ef> t_i^f [stored (Am,Ei)]
 
74
          (2) Z_MbEi = <Mb|Ef> t_i^f   [stored (Mb,Ei)]
 
75
          (3) WAbEi <-- t_m^b Z_mAEi - t_M^A Z_MbEi
 
76
       Store targets in:  W(Ei,Ab)  and   W'(Ab,Ei)  
 
77
  **/
 
78
 
 
79
  /** Z(Am,Ei) <-- - <Am|Ef> t_i^f **/
 
80
  dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
 
81
  dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
 
82
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
83
  dpd_contract424(&F, &T1, &Z, 3, 1, 0, -1, 0);
 
84
  dpd_file2_close(&T1);
 
85
  dpd_buf4_close(&F);
 
86
  dpd_buf4_close(&Z);
 
87
 
 
88
  /** t_m^b Z(Am,Ei) --> W(Ei,Ab) **/
 
89
  dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
 
90
  dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
 
91
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
92
  dpd_contract424(&Z, &T1, &W, 1, 0, 0, 1, 1);
 
93
  dpd_file2_close(&T1);
 
94
  dpd_buf4_close(&Z);
 
95
  dpd_buf4_close(&W);
 
96
 
 
97
  /** Z(Mb,Ei) <-- <Mb|Ef> t_i^f **/
 
98
  dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
 
99
  dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
 
100
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
101
  dpd_contract424(&F, &T1, &Z, 3, 1, 0, 1, 0);
 
102
  dpd_file2_close(&T1);
 
103
  dpd_buf4_close(&F);
 
104
  dpd_buf4_close(&Z);
 
105
 
 
106
  /** - t_M^A Z(Mb,Ei) --> W'(Ab,Ei) **/
 
107
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
108
  dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
 
109
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
110
  dpd_contract244(&T1, &Z, &W, 0, 0, 0, -1, 1);
 
111
  dpd_file2_close(&T1);
 
112
  dpd_buf4_close(&Z);
 
113
  dpd_buf4_close(&W);
 
114
 
 
115
  /**** Term V ****/
 
116
 
 
117
  /** WAbEi <-- tau_Mn^Ab <Mn|Ef> t_i^f
 
118
      Evaluate in two steps:
 
119
         (1) Z_MnEi = <Mn|Ef> t_i^f
 
120
         (2) WAbEi <-- tau_Mn^Ab Z_MnEi
 
121
      Store target in W'(Ab,Ei)
 
122
  **/
 
123
 
 
124
  /** Z(Mn,Ei) <-- <Mn|Ef> t_i^f **/
 
125
  dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z(Mn,Ei)");
 
126
  dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
127
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
128
  dpd_contract424(&D, &T1, &Z, 3, 1, 0, 1, 0);
 
129
  dpd_file2_close(&T1);
 
130
  dpd_buf4_close(&D);
 
131
 
 
132
  /** tau_Mn^Ab Z1(Mn,Ei) --> W'(Ab,Ei) **/
 
133
  dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z(Mn,Ei)");
 
134
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
135
  dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
 
136
  dpd_contract444(&T2, &Z, &W, 1, 1, 1, 1);
 
137
  dpd_buf4_close(&T2);
 
138
  dpd_buf4_close(&W);
 
139
  dpd_buf4_close(&Z);
 
140
 
 
141
  /**** Term VI ****/
 
142
 
 
143
  /** tau_Mn^Ab <Mn|Ei> --> Z(Ab,Ei) **/
 
144
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
145
  dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
 
146
  dpd_buf4_init(&E, CC_EINTS, 0, 22, 26, 22, 26, 0, "E <Ij|Ak>");
 
147
  dpd_contract444(&T2, &E, &W, 1, 1, 1, 1);
 
148
  dpd_buf4_close(&E);
 
149
  dpd_buf4_close(&T2);
 
150
  dpd_buf4_close(&W);
 
151
 
 
152
  /**** Term VII ****/
 
153
 
 
154
  /** WAbEi <-- <bM|fE> t_iM^fA - <AM||EF> t_iM^bF - <Am|Ef> t_im^bf
 
155
      Evaluate in five steps:
 
156
        (1) Sort <bM|fE> to F(bE,Mf) ordering.  (** Note that we assume a sort has already 
 
157
            been done for <AM||EF> and <Am|Ef> in the WABEI and Wabei terms. **)
 
158
        (2) Z'(bE,iA) = F(bE,Mf) T(iA,Mf)
 
159
        (3) Sort Z'(bE,iA) --> W(Ei,Ab)
 
160
        (4) Z''(AE,ib) = - F(AE,MF) T(ib,MF) - F(AE,mf) T(ib,mf)
 
161
        (5) Sort Z''(AE,ib) --> W(Ei,Ab)
 
162
 
 
163
      NB: The storage for the sorts is expensive and will eventually require out-of-core 
 
164
          codes.
 
165
  **/
 
166
 
 
167
  /** <bM|fE> --> F(bE,Mf) **/
 
168
  dpd_buf4_init(&F, CC_FINTS, 0, 25, 29, 25, 29, 0, "F <aI|bC>");
 
169
  dpd_buf4_sort(&F, CC_FINTS, psqr, 29, 24, "F <aI|bC> (aC,Ib)");
 
170
  dpd_buf4_close(&F);
 
171
 
 
172
  /** <bM|fE> t_iM^fA --> Z(bE,iA) **/
 
173
  dpd_buf4_init(&Z, CC_TMP0, 0, 29, 27, 29, 27, 0, "Z(bE,iA)");
 
174
  dpd_buf4_init(&F, CC_FINTS, 0, 29, 24, 29, 24, 0, "F <aI|bC> (aC,Ib)");
 
175
  dpd_buf4_init(&T2, CC_TAMPS, 0, 27, 24, 27, 24, 0, "tiBJa");
 
176
  dpd_contract444(&F, &T2, &Z, 0, 0, -1, 0);
 
177
  dpd_buf4_close(&T2);
 
178
  dpd_buf4_close(&F);
 
179
  dpd_buf4_close(&Z);
 
180
 
 
181
  /** Z(bE,iA) --> W(Ei,Ab) **/
 
182
  dpd_buf4_init(&Z, CC_TMP0, 0, 29, 27, 29, 27, 0, "Z(bE,iA)");
 
183
  dpd_buf4_sort_axpy(&Z, CC_HBAR, qrsp, 26, 28, "WEiAb", 1);
 
184
  dpd_buf4_close(&Z);
 
185
 
 
186
  /** Z''(AE,ib) <-- - <AM||EF> t_iM^bF **/
 
187
  dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
 
188
  dpd_buf4_init(&F, CC_FINTS, 0, 5, 20, 5, 20, 0, "F <AI||BC> (AB,IC)");
 
189
  dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 20, 30, 20, 0, "tiaJB");
 
190
  dpd_contract444(&F, &T2, &Z, 0, 0, 1, 0);
 
191
  dpd_buf4_close(&T2);
 
192
  dpd_buf4_close(&F);
 
193
  dpd_buf4_close(&Z);
 
194
 
 
195
  /** Z''(AE,ib) <-- -<Am|Ef> t_im^bf **/
 
196
  dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
 
197
  dpd_buf4_init(&F, CC_FINTS, 0, 5, 30, 5, 30, 0, "F <Ai|Bc> (AB,ic)");
 
198
  dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 30, 30, 30, 0, "tiajb");
 
199
  dpd_contract444(&F, &T2, &Z, 0, 0, 1, 1);
 
200
  dpd_buf4_close(&T2);
 
201
  dpd_buf4_close(&F);
 
202
  dpd_buf4_close(&Z);
 
203
 
 
204
  /** Z''(AE,ib) --> W(Ei,Ab) **/
 
205
  dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
 
206
  dpd_buf4_sort_axpy(&Z, CC_HBAR, qrps, 26, 28, "WEiAb", 1);
 
207
  dpd_buf4_close(&Z);
 
208
 
 
209
  /**** Terms VIII and IX ****/
 
210
 
 
211
  /** WAbEi <-- - t_M^A { <Mb|Ei> + t_iN^bF <MN||EF> + t_in^bf <Mn|Ef> }
 
212
                + t_m^b {-<mA|iE> + t_iN^fA <mN|fE> } 
 
213
      Evaluate in three steps: 
 
214
         (1) Z_MbEi =  <Mb|Ei> + t_iN^bF <MN||EF> + tin^bf <Mn|Ef>  [stored (Mb,Ei)]
 
215
         (2) Z_mAEi = -<mA|iE> + t_iN^fA <mN|fE>                    [stored (Am,Ei)]
 
216
         (3) WAbEi <-- - t_M^A Z_MbEi + t_m^b Z_mAEi
 
217
      Store targets in     W'(Ab,Ei) and  W(Ei,AB)
 
218
  **/
 
219
 
 
220
  /** Z(Mb,Ei) <-- <Mb|Ei> **/
 
221
  dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
 
222
  dpd_buf4_copy(&D, CC_TMP0, "Z(Mb,Ei)");
 
223
  dpd_buf4_close(&D);
 
224
 
 
225
  /** <MN||EF> t_iN^bF --> Z(ME,ib) **/
 
226
  dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
 
227
  dpd_buf4_init(&D, CC_DINTS, 0, 20, 20, 20, 20, 0, "D <IJ||AB> (IA,JB)");
 
228
  dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 20, 30, 20, 0, "tiaJB");
 
229
  dpd_contract444(&D, &T2, &Z, 0, 0, 1, 0);
 
230
  dpd_buf4_close(&T2);
 
231
  dpd_buf4_close(&D);
 
232
  dpd_buf4_close(&Z);
 
233
 
 
234
  /** <Mn|Ef> t_in^bf --> Z(ME,ib) **/
 
235
  dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
 
236
  dpd_buf4_init(&D, CC_DINTS, 0, 20, 30, 20, 30, 0, "D <Ij|Ab> (IA,jb)");
 
237
  dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 30, 30, 30, 0, "tiajb");
 
238
  dpd_contract444(&D, &T2, &Z, 0, 0, 1, 1);
 
239
  dpd_buf4_close(&T2);
 
240
  dpd_buf4_close(&D);
 
241
  dpd_buf4_close(&Z);
 
242
 
 
243
  /** Z(ME,ib) --> Z(Mb,Ei) **/
 
244
  dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
 
245
  dpd_buf4_sort_axpy(&Z, CC_TMP0, psqr, 24, 26, "Z(Mb,Ei)", 1);
 
246
  dpd_buf4_close(&Z);
 
247
 
 
248
  /** W'(Ab,Ei) <-- - t_M^A Z(Mb,Ei) **/
 
249
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
250
  dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
 
251
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
252
  dpd_contract244(&T1, &Z, &W, 0, 0, 0, -1, 1);
 
253
  dpd_file2_close(&T1);
 
254
  dpd_buf4_close(&Z);
 
255
  dpd_buf4_close(&W);
 
256
 
 
257
  /** Z(Am,Ei) <-- - <mA|iE> **/
 
258
  dpd_buf4_init(&C, CC_CINTS, 0, 26, 26, 26, 26, 0, "C <Ai|Bj>");
 
259
  dpd_buf4_copy(&C, CC_TMP0, "Z(Am,Ei)");
 
260
  dpd_buf4_close(&C);
 
261
  dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
 
262
  dpd_buf4_scm(&Z, -1);
 
263
  dpd_buf4_close(&Z);
 
264
 
 
265
  /** Z(mE,iA) <-- t_iN^fA <mN|fE> **/
 
266
  dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z(mE,iA)");
 
267
  dpd_buf4_init(&D, CC_DINTS, 0, 27, 24, 27, 24, 0, "D <iJ|aB> (iB,Ja)");
 
268
  dpd_buf4_init(&T2, CC_TAMPS, 0, 27, 24, 27, 24, 0, "tiBJa");
 
269
  dpd_contract444(&D, &T2, &Z, 0, 0, 1, 0);
 
270
  dpd_buf4_close(&T2);
 
271
  dpd_buf4_close(&D);
 
272
  dpd_buf4_close(&Z);
 
273
 
 
274
  /** Z(mE,iA) --> Z(Am,Ei) **/
 
275
  dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z(mE,iA)");
 
276
  dpd_buf4_sort_axpy(&Z, CC_TMP0, spqr, 26, 26, "Z(Am,Ei)", 1);
 
277
  dpd_buf4_close(&Z);
 
278
 
 
279
  /** W(Ei,AB) <-- t_m^b Z_mAEi **/
 
280
  dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
 
281
  dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
 
282
  dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
 
283
  dpd_contract424(&Z, &T1, &W, 1, 0, 0, 1, 1);
 
284
  dpd_file2_close(&T1);
 
285
  dpd_buf4_close(&Z);
 
286
  dpd_buf4_close(&W);
 
287
 
 
288
  /**** Combine accumulated W'(Ab,Ei) and W(Ei,Ab) terms into WEiAb ****/
 
289
  dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
 
290
  dpd_buf4_sort_axpy(&W, CC_HBAR, rspq, 26, 28, "WEiAb", 1);
 
291
  dpd_buf4_close(&W);
 
292
}
 
293
 
 
294
}} // namespace psi::cchbar