← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/deanti_ROHF.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/ccdensity/deanti_ROHF.c
1
 
#include <stdio.h>
2
 
#include <libiwl/iwl.h>
3
 
#include <libdpd/dpd.h>
4
 
#include <psifiles.h>
5
 
#define EXTERN
6
 
#include "globals.h"
7
 
 
8
 
/* DEANTI_ROHF(): Convert the ROHF two-particle density from Dirac to
9
 
** Mulliken ordering.  The original, Fock-adjusted density (see the
10
 
** comments in fock.c) corresponds to a two-electron energy (or energy
11
 
** derivative) expression of the form:
12
 
**
13
 
** E(TWO) = 1/4 sum_pqrs Gpqrs <pq||rs>
14
 
**
15
 
** However, the derivative two-electron integrals are produced in
16
 
** Mulliken-ordered, symmetric form rather than Dirac-ordered
17
 
** antisymmetric form.  This code alters the two-particle density
18
 
** matrix ordering for the energy expression of the form:
19
 
**
20
 
** E(TWO) = 1/2 sum_pqrs Gpqrs <pq|rs>
21
 
**
22
 
** The final conversion to Mulliken ordering is taken care of in
23
 
** dump.c
24
 
**
25
 
** The second equation above may be derived via
26
 
**
27
 
** E(TWO) = 1/4 sum_pqrs Gpqrs (<pq|rs> - <pq|sr>)
28
 
**        = 1/4 sum_pqrs Gpqrs <pq|rs> - 1/4 sum_pqrs Gpqrs <pq|sr>
29
 
**        = 1/4 sum_pqrs Gpqrs <pq|rs> - 1/4 sum_pqrs Gpqsr <pq|rs>
30
 
**        = 1/4 sum_pqrs (Gpqrs - Gpqsr) <pq|rs>
31
 
**        = 1/2 sum_pqrs Gpqrs <pq|rs>
32
 
**
33
 
** Equations for the individual components are given explicitly in
34
 
** comments below.
35
 
** */
36
 
 
37
 
void deanti_ROHF(struct RHO_Params rho_params)
38
 
{
39
 
  dpdbuf4 G1, G2;
40
 
  dpdfile2 D, F;
41
 
  double one_energy = 0.0, two_energy = 0.0, total_two_energy = 0.0;
42
 
  dpdbuf4 A, B, C, DInts, E, FInts;
43
 
 
44
 
  if(!params.aobasis) {
45
 
    fprintf(outfile, "\n\tEnergies re-computed from Mulliken density:\n");
46
 
    fprintf(outfile,   "\t-------------------------------------------\n");
47
 
 
48
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 0, rho_params.DIJ_lbl);
49
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
50
 
    one_energy += dpd_file2_dot(&D, &F);
51
 
    dpd_file2_close(&F);
52
 
    dpd_file2_close(&D);
53
 
 
54
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 0, rho_params.Dij_lbl);
55
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
56
 
    one_energy += dpd_file2_dot(&D, &F);
57
 
    dpd_file2_close(&F);
58
 
    dpd_file2_close(&D);
59
 
 
60
 
    dpd_file2_init(&D, CC_OEI, 0, 1, 1, rho_params.DAB_lbl);
61
 
    dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
62
 
    one_energy += dpd_file2_dot(&D, &F);
63
 
    dpd_file2_close(&F);
64
 
    dpd_file2_close(&D);
65
 
 
66
 
    dpd_file2_init(&D, CC_OEI, 0, 1, 1, rho_params.Dab_lbl);
67
 
    dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
68
 
    one_energy += dpd_file2_dot(&D, &F);
69
 
    dpd_file2_close(&F);
70
 
    dpd_file2_close(&D);
71
 
 
72
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DIA_lbl);
73
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
74
 
    one_energy += dpd_file2_dot(&D, &F);
75
 
    dpd_file2_close(&F);
76
 
    dpd_file2_close(&D);
77
 
 
78
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.Dia_lbl);
79
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
80
 
    one_energy += dpd_file2_dot(&D, &F);
81
 
    dpd_file2_close(&F);
82
 
    dpd_file2_close(&D);
83
 
 
84
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DAI_lbl);
85
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
86
 
    one_energy += dpd_file2_dot(&D, &F);
87
 
    dpd_file2_close(&F);
88
 
    dpd_file2_close(&D);
89
 
 
90
 
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.Dai_lbl);
91
 
    dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
92
 
    one_energy += dpd_file2_dot(&D, &F);
93
 
    dpd_file2_close(&F);
94
 
    dpd_file2_close(&D);
95
 
 
96
 
    fprintf(outfile, "\tOne-electron energy        = %20.15f\n", one_energy);
97
 
    fflush(outfile);
98
 
  }
99
 
 
100
 
  /* G(Ij,Kl) <-- 1/2 G(IJ,KL) + 1/2 G(ij,kl) + 1/2 G(Ij,Kl) + 1/2 G(iJ,kL) */
101
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 0, 0, 0, 0, "GIjKl");
102
 
  dpd_buf4_sort(&G1, CC_TMP0, qprs, 0, 0, "GjIKl");
103
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 0, 0, 0, 0, "GjIKl");
104
 
  dpd_buf4_sort(&G2, CC_TMP0, pqsr, 0, 0, "GiJkL");
105
 
  dpd_buf4_close(&G2);
106
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 0, 0, 0, 0, "GiJkL");
107
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
108
 
  dpd_buf4_close(&G2);
109
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 0, 2, 2, 0, "GIJKL");
110
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
111
 
  dpd_buf4_close(&G2);
112
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 0, 2, 2, 0, "Gijkl");
113
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
114
 
  dpd_buf4_close(&G2);
115
 
  dpd_buf4_scm(&G1, 0.5);
116
 
  
117
 
  if(!params.aobasis) {
118
 
    dpd_buf4_init(&A, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
119
 
    two_energy = dpd_buf4_dot(&A, &G1);
120
 
    dpd_buf4_close(&A);
121
 
    fprintf(outfile, "\tIJKL energy                = %20.15f\n", two_energy);
122
 
    total_two_energy += two_energy;
123
 
  }
124
 
 
125
 
  dpd_buf4_close(&G1);
126
 
 
127
 
  /* G(IJ,KA) <-- G(IJ,KA) + G(ij,ka) + G(Ij,Ka) + G(iJ,kA) */
128
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GIjKa");
129
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GiJkA");
130
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
131
 
  dpd_buf4_close(&G2);
132
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 2, 10, 0, "GIJKA");
133
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
134
 
  dpd_buf4_close(&G2);
135
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 2, 10, 0, "Gijka");
136
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
137
 
  dpd_buf4_close(&G2);
138
 
 
139
 
  if(!params.aobasis) {
140
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
141
 
    two_energy = dpd_buf4_dot(&E, &G1);
142
 
    dpd_buf4_close(&E);
143
 
    fprintf(outfile, "\tIJKA energy                = %20.15f\n", 2*two_energy);
144
 
    total_two_energy += 2*two_energy;
145
 
  }
146
 
 
147
 
  dpd_buf4_close(&G1);
148
 
 
149
 
  /* G(Ij,Ab) <-- G(Ij,Ab) - G(Ib,jA) */
150
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
151
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbjA");
152
 
  dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "GIbjA (Ij,Ab)");
153
 
  dpd_buf4_close(&G2);
154
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "GIbjA (Ij,Ab)");
155
 
  dpd_buf4_axpy(&G2, &G1, -1.0);
156
 
  dpd_buf4_close(&G2);
157
 
  dpd_buf4_close(&G1);
158
 
  
159
 
  /* G(IJ,AB) <-- G(IJ,AB) - G(IB,JA) */
160
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 2, 7, 0, "GIJAB");
161
 
  dpd_buf4_copy(&G1, CC_TMP0, "G(IJ,AB)");
162
 
  dpd_buf4_close(&G1);
163
 
  dpd_buf4_init(&G1, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(IJ,AB)");
164
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
165
 
  dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "GIBJA (IJ,AB)");
166
 
  dpd_buf4_close(&G2);
167
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "GIBJA (IJ,AB)");
168
 
  dpd_buf4_axpy(&G2, &G1, -1.0); 
169
 
  dpd_buf4_close(&G2);
170
 
  dpd_buf4_close(&G1);
171
 
 
172
 
  /* G(ij,ab) <-- G(ij,ab) - G(ib,ja) */
173
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 2, 7, 0, "Gijab");
174
 
  dpd_buf4_copy(&G1, CC_TMP0, "G(ij,ab)");
175
 
  dpd_buf4_close(&G1);
176
 
  dpd_buf4_init(&G1, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(ij,ab)");
177
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
178
 
  dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "Gibja (ij,ab)");
179
 
  dpd_buf4_close(&G2);
180
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "Gibja (ij,ab)");
181
 
  dpd_buf4_axpy(&G2, &G1, -1.0);
182
 
  dpd_buf4_close(&G2);
183
 
  dpd_buf4_close(&G1);
184
 
 
185
 
  /* G(IJ,AB) <-- G(IJ,AB) + G(ij,ab) + G(Ij,Ab) + G(iJ,aB) */
186
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
187
 
  dpd_buf4_sort(&G1, CC_TMP0, qpsr, 0, 5, "G(jI,bA)");
188
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(jI,bA)");
189
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
190
 
  dpd_buf4_close(&G2);
191
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(IJ,AB)");
192
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
193
 
  dpd_buf4_close(&G2);
194
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(ij,ab)");
195
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
196
 
  dpd_buf4_close(&G2);
197
 
 
198
 
  if(!params.aobasis) {
199
 
    dpd_buf4_init(&DInts, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
200
 
    two_energy = dpd_buf4_dot(&DInts, &G1);
201
 
    dpd_buf4_close(&DInts);
202
 
    fprintf(outfile, "\tIJAB energy                = %20.15f\n", two_energy);
203
 
    total_two_energy += two_energy;
204
 
  }
205
 
 
206
 
  dpd_buf4_close(&G1);
207
 
  
208
 
  /* G(IB,JA) <-- G(IB,JA) + G(ib,ja) + G(Ib,Ja) + G(iB,jA) */
209
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
210
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
211
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
212
 
  dpd_buf4_close(&G2); 
213
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbJa");
214
 
  dpd_buf4_axpy(&G2, &G1, 1.0); 
215
 
  dpd_buf4_close(&G2); 
216
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GiBjA");
217
 
  dpd_buf4_axpy(&G2, &G1, 1.0); 
218
 
  dpd_buf4_close(&G2);
219
 
 
220
 
  if(!params.aobasis) {
221
 
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
222
 
    two_energy = dpd_buf4_dot(&C, &G1);
223
 
    dpd_buf4_close(&C);
224
 
    fprintf(outfile, "\tIBJA energy                = %20.15f\n", two_energy);
225
 
    total_two_energy += two_energy;
226
 
  }
227
 
 
228
 
  dpd_buf4_close(&G1);
229
 
 
230
 
  /* G(CI,AB) <-- G(CI,AB) + G(ci,ab) + G(Ci,Ab) + G(cI,aB) */
231
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
232
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
233
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
234
 
  dpd_buf4_close(&G2); 
235
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 7, 0, "GCIAB");
236
 
  dpd_buf4_axpy(&G2, &G1, 1.0);  
237
 
  dpd_buf4_close(&G2); 
238
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 7, 0, "Gciab");
239
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
240
 
  dpd_buf4_close(&G2);
241
 
 
242
 
  if(!params.aobasis) {
243
 
    dpd_buf4_init(&FInts, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
244
 
    dpd_buf4_sort(&FInts, CC_TMP0, qprs, 11, 5, "F(CI,BA)");
245
 
    dpd_buf4_close(&FInts);
246
 
    dpd_buf4_init(&FInts, CC_TMP0, 0, 11, 5, 11, 5, 0, "F(CI,BA)");
247
 
    dpd_buf4_sort(&FInts, CC_TMP1, pqsr, 11, 5, "F(CI,AB)");
248
 
    dpd_buf4_close(&FInts);
249
 
    dpd_buf4_init(&FInts, CC_TMP1, 0, 11, 5, 11, 5, 0, "F(CI,AB)");
250
 
    two_energy = dpd_buf4_dot(&FInts, &G1);
251
 
    dpd_buf4_close(&FInts);
252
 
    fprintf(outfile, "\tCIAB energy                = %20.15f\n", 2*two_energy);
253
 
    total_two_energy += 2*two_energy;
254
 
  }
255
 
 
256
 
  dpd_buf4_close(&G1);
257
 
 
258
 
  /* G(Ab,Cd) <-- 1/2 G(AB,CD) + 1/2 G(ab,cd) + G(Ab,Cd) */
259
 
  dpd_buf4_init(&G1, CC_GAMMA, 0, 5, 5, 5, 5, 0, "GAbCd");
260
 
  dpd_buf4_sort(&G1, CC_TMP0, qprs, 5, 5, "G(bA,Cd)");
261
 
  dpd_buf4_init(&G2, CC_TMP0, 0, 5, 5, 5, 5, 0, "G(bA,Cd)");
262
 
  dpd_buf4_sort(&G2, CC_TMP1, pqsr, 5, 5, "G(aB,cD)");
263
 
  dpd_buf4_close(&G2);
264
 
  dpd_buf4_init(&G2, CC_TMP1, 0, 5, 5, 5, 5, 0, "G(aB,cD)");
265
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
266
 
  dpd_buf4_close(&G2);
267
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 5, 5, 7, 7, 0, "GABCD");
268
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
269
 
  dpd_buf4_close(&G2);
270
 
  dpd_buf4_init(&G2, CC_GAMMA, 0, 5, 5, 7, 7, 0, "Gabcd");
271
 
  dpd_buf4_axpy(&G2, &G1, 1.0);
272
 
  dpd_buf4_close(&G2);
273
 
  dpd_buf4_scm(&G1, 0.5);
274
 
 
275
 
  if(!params.aobasis) {  
276
 
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
277
 
    two_energy = dpd_buf4_dot(&B, &G1);
278
 
    dpd_buf4_close(&B);
279
 
    fprintf(outfile, "\tABCD energy                = %20.15f\n", two_energy);
280
 
    total_two_energy += two_energy;
281
 
  }
282
 
 
283
 
  dpd_buf4_close(&G1);
284
 
 
285
 
  if(!params.aobasis) {
286
 
    fprintf(outfile, "\tTotal two-electron energy  = %20.15f\n", total_two_energy);
287
 
    if (params.ground) {
288
 
      fprintf(outfile, "\tCCSD correlation energy    = %20.15f\n",
289
 
              one_energy + total_two_energy);
290
 
      fprintf(outfile, "\tTotal CCSD energy          = %20.15f\n",
291
 
              one_energy + total_two_energy + moinfo.eref);
292
 
    }
293
 
    else {
294
 
      fprintf(outfile, "\tTotal EOM CCSD correlation energy        = %20.15f\n",
295
 
          one_energy + total_two_energy);
296
 
      fprintf(outfile, "\tCCSD correlation + EOM excitation energy = %20.15f\n",
297
 
          moinfo.ecc + params.cceom_energy);
298
 
      fprintf(outfile, "\tTotal EOM CCSD energy                    = %20.15f\n",
299
 
          one_energy + total_two_energy + moinfo.eref);
300
 
    }
301
 
  }
302
 
}