← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/energy_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:
src/bin/ccdensity/energy_UHF.cc
 
1
/*! \file
 
2
    \ingroup CCDENSITY
 
3
    \brief Calculates the one- and two-electron CC energies using the
 
4
    coresponding one- and two-particle density matrices. 
 
5
*/
 
6
#include <stdio.h>
 
7
#include <strings.h>
 
8
#include <string.h>
 
9
#include <libdpd/dpd.h>
 
10
#include "MOInfo.h"
 
11
#include "Params.h"
 
12
#include "Frozen.h"
 
13
#define EXTERN
 
14
#include "globals.h"
 
15
 
 
16
namespace psi { namespace ccdensity {
 
17
 
 
18
    /* ENERGY_UHF(): Compute the UHF CC energy using the one- and two-particle
 
19
    ** density matrices.
 
20
    */
 
21
 
 
22
    void energy_UHF(struct RHO_Params rho_params)
 
23
    {
 
24
      dpdfile2 D, F;
 
25
      dpdbuf4 G, A, B, C, DInts, E, FInts;
 
26
      double one_energy=0.0, two_energy=0.0, total_two_energy = 0.0;
 
27
      double this_energy, test_energy;
 
28
 
 
29
      fprintf(outfile, "\n\tEnergies re-computed from CC density:\n");
 
30
      fprintf(outfile,   "\t-------------------------------------\n");
 
31
 
 
32
      dpd_file2_init(&D, CC_OEI, 0, 0, 0, rho_params.DIJ_lbl);
 
33
      dpd_file2_init(&F, CC_OEI, 0, 0, 0, "fIJ");
 
34
      this_energy = dpd_file2_dot(&D, &F);
 
35
      dpd_file2_close(&F);
 
36
      dpd_file2_close(&D);
 
37
 
 
38
      /*    fprintf(outfile, "\tDIJ = %20.15f\n", this_energy); */
 
39
      one_energy += this_energy;
 
40
 
 
41
      dpd_file2_init(&D, CC_OEI, 0, 2, 2, rho_params.Dij_lbl);
 
42
      dpd_file2_init(&F, CC_OEI, 0, 2, 2, "fij");
 
43
      this_energy = dpd_file2_dot(&D, &F);
 
44
      dpd_file2_close(&F);
 
45
      dpd_file2_close(&D);
 
46
 
 
47
      /*    fprintf(outfile, "\tDij = %20.15f\n", this_energy); */
 
48
      one_energy += this_energy;
 
49
 
 
50
      dpd_file2_init(&D, CC_OEI, 0, 1, 1, rho_params.DAB_lbl);
 
51
      dpd_file2_init(&F, CC_OEI, 0, 1, 1, "fAB");
 
52
      this_energy = dpd_file2_dot(&D, &F);
 
53
      dpd_file2_close(&F);
 
54
      dpd_file2_close(&D);
 
55
 
 
56
      /*    fprintf(outfile, "\tDAB = %20.15f\n", this_energy); */
 
57
      one_energy += this_energy;
 
58
 
 
59
      dpd_file2_init(&D, CC_OEI, 0, 3, 3, rho_params.Dab_lbl);
 
60
      dpd_file2_init(&F, CC_OEI, 0, 3, 3, "fab");
 
61
      this_energy = dpd_file2_dot(&D, &F);
 
62
      dpd_file2_close(&F);
 
63
      dpd_file2_close(&D);
 
64
 
 
65
      /*    fprintf(outfile, "\tDab = %20.15f\n", this_energy); */
 
66
      one_energy += this_energy;
 
67
 
 
68
      dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DIA_lbl);
 
69
      dpd_file2_init(&F, CC_OEI, 0, 0, 1, "fIA");
 
70
      this_energy = dpd_file2_dot(&D, &F);
 
71
      dpd_file2_close(&F);
 
72
      dpd_file2_close(&D);
 
73
 
 
74
      /*    fprintf(outfile, "\tDIA = %20.15f\n", this_energy); */
 
75
      one_energy += this_energy;
 
76
 
 
77
      dpd_file2_init(&D, CC_OEI, 0, 2, 3, rho_params.Dia_lbl);
 
78
      dpd_file2_init(&F, CC_OEI, 0, 2, 3, "fia");
 
79
      this_energy = dpd_file2_dot(&D, &F);
 
80
      dpd_file2_close(&F);
 
81
      dpd_file2_close(&D);
 
82
 
 
83
      /*    fprintf(outfile, "\tDia = %20.15f\n", this_energy); */
 
84
      one_energy += this_energy;
 
85
 
 
86
      dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DAI_lbl);
 
87
      dpd_file2_init(&F, CC_OEI, 0, 0, 1, "fIA");
 
88
      this_energy = dpd_file2_dot(&D, &F);
 
89
      dpd_file2_close(&F);
 
90
      dpd_file2_close(&D);
 
91
 
 
92
      /*    fprintf(outfile, "\tDAI = %20.15f\n", this_energy); */
 
93
      one_energy += this_energy;
 
94
 
 
95
      dpd_file2_init(&D, CC_OEI, 0, 2, 3, rho_params.Dai_lbl);
 
96
      dpd_file2_init(&F, CC_OEI, 0, 2, 3, "fia");
 
97
      this_energy = dpd_file2_dot(&D, &F);
 
98
      dpd_file2_close(&F);
 
99
      dpd_file2_close(&D);
 
100
      /*    fprintf(outfile, "\tDai = %20.15f\n", this_energy); */
 
101
      one_energy += this_energy;
 
102
 
 
103
      fprintf(outfile, "\tOne-electron energy        = %20.15f\n", one_energy);
 
104
      fflush(outfile);
 
105
      if (params.onepdm) return;
 
106
 
 
107
      total_two_energy = 0.0;
 
108
 
 
109
 
 
110
      two_energy = 0.0;
 
111
 
 
112
      dpd_buf4_init(&A, CC_AINTS, 0, 2, 2, 0, 0, 1, "A <IJ|KL>");
 
113
      dpd_buf4_init(&G, CC_GAMMA, 0, 2, 2, 2, 2, 0, "GIJKL");
 
114
      two_energy += dpd_buf4_dot(&G, &A);
 
115
      dpd_buf4_close(&G);
 
116
      dpd_buf4_close(&A);
 
117
 
 
118
      dpd_buf4_init(&A, CC_AINTS, 0, 12, 12, 10, 10, 1, "A <ij|kl>");
 
119
      dpd_buf4_init(&G, CC_GAMMA, 0, 12, 12, 12, 12, 0, "Gijkl");
 
120
      two_energy += dpd_buf4_dot(&G, &A);
 
121
      dpd_buf4_close(&G);
 
122
      dpd_buf4_close(&A);
 
123
 
 
124
      dpd_buf4_init(&A, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
 
125
      dpd_buf4_init(&G, CC_GAMMA, 0, 22, 22, 22, 22, 0, "GIjKl");
 
126
      two_energy += dpd_buf4_dot(&G, &A);
 
127
      dpd_buf4_close(&G);
 
128
      dpd_buf4_close(&A);
 
129
 
 
130
      total_two_energy += two_energy;
 
131
      fprintf(outfile, "\tIJKL energy                = %20.15f\n", two_energy);
 
132
      fflush(outfile);
 
133
 
 
134
      two_energy = 0.0;
 
135
 
 
136
      dpd_buf4_init(&E, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
 
137
      dpd_buf4_init(&G, CC_GAMMA, 0, 2, 20, 2, 20, 0, "GIJKA");
 
138
      two_energy += dpd_buf4_dot(&G, &E);
 
139
      dpd_buf4_close(&G);
 
140
      dpd_buf4_close(&E);
 
141
 
 
142
      dpd_buf4_init(&E, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
 
143
      dpd_buf4_init(&G, CC_GAMMA, 0, 12, 30, 12, 30, 0, "Gijka");
 
144
      two_energy += dpd_buf4_dot(&G, &E);
 
145
      dpd_buf4_close(&G);
 
146
      dpd_buf4_close(&E);
 
147
 
 
148
      dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
 
149
      dpd_buf4_init(&G, CC_GAMMA, 0, 22, 24, 22, 24, 0, "GIjKa");
 
150
      two_energy += dpd_buf4_dot(&G, &E);
 
151
      dpd_buf4_close(&G);
 
152
      dpd_buf4_close(&E);
 
153
 
 
154
      dpd_buf4_init(&E, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
 
155
      dpd_buf4_init(&G, CC_GAMMA, 0, 23, 27, 23, 27, 0, "GiJkA");
 
156
      two_energy += dpd_buf4_dot(&G, &E);
 
157
      dpd_buf4_close(&G);
 
158
      dpd_buf4_close(&E);
 
159
 
 
160
      two_energy *= 2;
 
161
      total_two_energy += two_energy;
 
162
      fprintf(outfile, "\tIJKA energy                = %20.15f\n", two_energy);
 
163
      fflush(outfile);
 
164
 
 
165
 
 
166
      two_energy = 0.0;
 
167
 
 
168
      dpd_buf4_init(&DInts, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
 
169
      dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "GIJAB");
 
170
      two_energy += dpd_buf4_dot(&G, &DInts);
 
171
      dpd_buf4_close(&G);
 
172
      dpd_buf4_close(&DInts);
 
173
 
 
174
      dpd_buf4_init(&DInts, CC_DINTS, 0, 12, 17, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
 
175
      dpd_buf4_init(&G, CC_GAMMA, 0, 12, 17, 12, 17, 0, "Gijab");
 
176
      two_energy += dpd_buf4_dot(&G, &DInts);
 
177
      dpd_buf4_close(&G);
 
178
      dpd_buf4_close(&DInts);
 
179
 
 
180
      dpd_buf4_init(&DInts, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
181
      dpd_buf4_init(&G, CC_GAMMA, 0, 22, 28, 22, 28, 0, "GIjAb");
 
182
      two_energy += dpd_buf4_dot(&G, &DInts);
 
183
      dpd_buf4_close(&G);
 
184
      dpd_buf4_close(&DInts);
 
185
 
 
186
      two_energy *= 2;
 
187
      total_two_energy += two_energy;
 
188
      fprintf(outfile, "\tIJAB energy                = %20.15f\n", two_energy);
 
189
      fflush(outfile);
 
190
 
 
191
      /*
 
192
      ** Compute the Gibja contribution to the two-electron energy.  By
 
193
      ** spin-case this contribution looks like:
 
194
      **
 
195
      **  E(AA) <-- sum_IBJA G(IB,JA) <JA||IB>
 
196
      **  E(BB) <-- sum_ibja G(ib,ja) <ja||ib>
 
197
      **  E(AB) <-- sum_IbJa ( G(Ib,Ja) <Ja|Ib> + G(iB,jA) <jA|iB> -
 
198
      **                         G(Ib,jA) <jA|bI> - G(iB,Ja) <Ja|Bi> )
 
199
      **
 
200
      **  See Gibja.c for the definition of G here.
 
201
      */
 
202
 
 
203
      two_energy = 0.0;
 
204
 
 
205
      dpd_buf4_init(&C, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
 
206
      dpd_buf4_init(&G, CC_GAMMA, 0, 20, 20, 20, 20, 0, "GIBJA");
 
207
      two_energy += dpd_buf4_dot(&G, &C);
 
208
      dpd_buf4_close(&G);
 
209
      dpd_buf4_close(&C);
 
210
 
 
211
      dpd_buf4_init(&C, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
 
212
      dpd_buf4_init(&G, CC_GAMMA, 0, 30, 30, 30, 30, 0, "Gibja");
 
213
      two_energy += dpd_buf4_dot(&G, &C);
 
214
      dpd_buf4_close(&G);
 
215
      dpd_buf4_close(&C);
 
216
 
 
217
      dpd_buf4_init(&C, CC_CINTS, 0, 24, 24, 24, 24, 0, "C <Ia|Jb>");
 
218
      dpd_buf4_init(&G, CC_GAMMA, 0, 24, 24, 24, 24, 0, "GIbJa");
 
219
      two_energy += dpd_buf4_dot(&G, &C);
 
220
      dpd_buf4_close(&G);
 
221
      dpd_buf4_close(&C);
 
222
 
 
223
      dpd_buf4_init(&C, CC_CINTS, 0, 27, 27, 27, 27, 0, "C <iA|jB>");
 
224
      dpd_buf4_init(&G, CC_GAMMA, 0, 27, 27, 27, 27, 0, "GiBjA");
 
225
      two_energy += dpd_buf4_dot(&G, &C);
 
226
      dpd_buf4_close(&G);
 
227
      dpd_buf4_close(&C);
 
228
 
 
229
      dpd_buf4_init(&DInts, CC_DINTS, 0, 24, 27, 24, 27, 0, "D <Ij|Ab> (Ib,jA)");
 
230
      dpd_buf4_init(&G, CC_GAMMA, 0, 24, 27, 24, 27, 0, "GIbjA");
 
231
      two_energy -= dpd_buf4_dot(&G, &DInts);
 
232
      dpd_buf4_close(&G);
 
233
      dpd_buf4_close(&DInts);
 
234
 
 
235
      dpd_buf4_init(&DInts, CC_DINTS, 0, 27, 24, 27, 24, 0, "D <iJ|aB> (iB,Ja)");
 
236
      dpd_buf4_init(&G, CC_GAMMA, 0, 27, 24, 27, 24, 0, "GiBJa");
 
237
      two_energy -= dpd_buf4_dot(&G, &DInts);
 
238
      dpd_buf4_close(&G);
 
239
      dpd_buf4_close(&DInts);
 
240
 
 
241
      total_two_energy += two_energy;
 
242
      fprintf(outfile, "\tIBJA energy                = %20.15f\n", two_energy);
 
243
      fflush(outfile);
 
244
 
 
245
 
 
246
      two_energy = 0.0;
 
247
 
 
248
      dpd_buf4_init(&FInts, CC_FINTS, 0, 21, 7, 21, 5, 1, "F <AI|BC>");
 
249
      dpd_buf4_init(&G, CC_GAMMA, 0, 21, 7, 21, 7, 0, "GCIAB");
 
250
      two_energy += dpd_buf4_dot(&G, &FInts);
 
251
      dpd_buf4_close(&G);
 
252
      dpd_buf4_close(&FInts);
 
253
 
 
254
      dpd_buf4_init(&FInts, CC_FINTS, 0, 31, 17, 31, 15, 1, "F <ai|bc>");
 
255
      dpd_buf4_init(&G, CC_GAMMA, 0, 31, 17, 31, 17, 0, "Gciab");
 
256
      two_energy += dpd_buf4_dot(&G, &FInts);
 
257
      dpd_buf4_close(&G);
 
258
      dpd_buf4_close(&FInts);
 
259
 
 
260
      dpd_buf4_init(&FInts, CC_FINTS, 0, 25, 29, 25, 29, 0, "F <aI|bC>");
 
261
      dpd_buf4_init(&G, CC_GAMMA, 0, 25, 29, 25, 29, 0, "GcIaB");
 
262
      two_energy += dpd_buf4_dot(&G, &FInts);
 
263
      dpd_buf4_close(&G);
 
264
      dpd_buf4_close(&FInts);
 
265
 
 
266
      dpd_buf4_init(&FInts, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
 
267
      dpd_buf4_init(&G, CC_GAMMA, 0, 26, 28, 26, 28, 0, "GCiAb");
 
268
      two_energy += dpd_buf4_dot(&G, &FInts); 
 
269
      dpd_buf4_close(&G);
 
270
      dpd_buf4_close(&FInts);
 
271
 
 
272
      two_energy *= 2;
 
273
      total_two_energy += two_energy;
 
274
      fprintf(outfile, "\tCIAB energy                = %20.15f\n", two_energy);
 
275
      fflush(outfile);
 
276
 
 
277
      two_energy = 0.0;
 
278
 
 
279
      dpd_buf4_init(&B, CC_BINTS, 0, 7, 7, 5, 5, 1, "B <AB|CD>");
 
280
      dpd_buf4_init(&G, CC_GAMMA, 0, 7, 7, 7, 7, 0, "GABCD");
 
281
      two_energy += dpd_buf4_dot(&G, &B);
 
282
      dpd_buf4_close(&G);
 
283
      dpd_buf4_close(&B);
 
284
 
 
285
      dpd_buf4_init(&B, CC_BINTS, 0, 17, 17, 15, 15, 1, "B <ab|cd>");
 
286
      dpd_buf4_init(&G, CC_GAMMA, 0, 17, 17, 17, 17, 0, "Gabcd");
 
287
      two_energy += dpd_buf4_dot(&G, &B);
 
288
      dpd_buf4_close(&G);
 
289
      dpd_buf4_close(&B);
 
290
 
 
291
      dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
292
      dpd_buf4_init(&G, CC_GAMMA, 0, 28, 28, 28, 28, 0, "GAbCd");
 
293
      two_energy += dpd_buf4_dot(&G, &B);
 
294
      dpd_buf4_close(&G);
 
295
      dpd_buf4_close(&B);
 
296
 
 
297
      total_two_energy += two_energy;
 
298
      fprintf(outfile, "\tABCD energy                = %20.15f\n", two_energy);
 
299
 
 
300
      fprintf(outfile, "\tTotal two-electron energy  = %20.15f\n", total_two_energy);
 
301
      if (params.ground) {
 
302
        fprintf(outfile, "\t%-7s correlation energy = %20.15f\n", !strcmp(params.wfn,"CCSD_T") ? "CCSD(T)" : params.wfn,
 
303
                one_energy + total_two_energy);
 
304
        fprintf(outfile, "\tTotal %-7s energy       = %20.15f\n", !strcmp(params.wfn,"CCSD_T") ? "CCSD(T)" : params.wfn,
 
305
                one_energy + total_two_energy + moinfo.eref);
 
306
      }
 
307
      else {
 
308
        fprintf(outfile, "\tTotal EOM CCSD correlation energy        = %20.15f\n",
 
309
                one_energy + total_two_energy);
 
310
        fprintf(outfile, "\tCCSD correlation + EOM excitation energy = %20.15f\n",
 
311
                moinfo.ecc + params.cceom_energy);
 
312
        fprintf(outfile, "\tTotal EOM CCSD energy                    = %20.15f\n",
 
313
                one_energy + total_two_energy + moinfo.eref);
 
314
      }
 
315
    }
 
316
 
 
317
  }} // namespace psi::ccdensity