← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/mp2/opdm.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/mp2/opdm.cc
 
1
/*! \file
 
2
    \ingroup MP2
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cmath>
 
6
#include <cstdlib>
 
7
#include <libdpd/dpd.h>
 
8
#include <psifiles.h>
 
9
#define EXTERN
 
10
#include "globals.h"
 
11
 
 
12
namespace psi{ namespace mp2{
 
13
 
 
14
void rhf_opdm(void);
 
15
void uhf_opdm(void);
 
16
void rhf_sf_opdm(void);
 
17
void uhf_sf_opdm(void);
 
18
 
 
19
void opdm(void)
 
20
{
 
21
  if(params.gradient) {
 
22
    if(params.ref == 0) rhf_opdm();
 
23
    else if(params.ref == 2) uhf_sf_opdm();
 
24
  }
 
25
  else {
 
26
    if(params.ref == 0) rhf_opdm();
 
27
    else if(params.ref == 2) uhf_opdm();
 
28
  }
 
29
}
 
30
 
 
31
void rhf_opdm(void)
 
32
{
 
33
  int h, i, j, a, b;
 
34
  int I, J, A, B;
 
35
  double trace_IJ=0.0;
 
36
  double trace_AB=0.0;
 
37
  double alpha;
 
38
  dpdfile2 DIJ, DAB, D;
 
39
  dpdbuf4 T2, T2A, T2B;
 
40
 
 
41
  if(params.gradient) alpha = 1.0;
 
42
  else alpha = 2.0;
 
43
 
 
44
  dpd_file2_init(&DIJ, CC_OEI, 0, 0, 0, "DIJ");
 
45
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
46
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 0, 5, 0, 5, 0, "2 tIjAb - tIjBa");
 
47
  dpd_contract442(&T2A, &T2B, &DIJ, 0, 0, -alpha, 0);
 
48
  dpd_buf4_close(&T2B);
 
49
  dpd_buf4_close(&T2A);
 
50
  trace_IJ = dpd_file2_trace(&DIJ);
 
51
  if(params.gradient) dpd_file2_copy(&DIJ, CC_OEI, "Dij");
 
52
  dpd_file2_close(&DIJ);
 
53
 
 
54
  dpd_file2_init(&DAB, CC_OEI, 0, 1, 1, "DAB");
 
55
  dpd_file2_scm(&DAB,0.0);
 
56
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
57
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 0, 5, 0, 5, 0, "2 tIjAb - tIjBa");
 
58
  dpd_contract442(&T2A, &T2B, &DAB, 3, 3, alpha, 0);
 
59
  dpd_buf4_close(&T2B);
 
60
  dpd_buf4_close(&T2A);
 
61
  trace_AB = dpd_file2_trace(&DAB);
 
62
  if(params.gradient) dpd_file2_copy(&DAB, CC_OEI, "Dab");
 
63
  dpd_file2_close(&DAB);
 
64
 
 
65
  if(params.gradient) {
 
66
    /* zero out Dia DIA Dai DAI */
 
67
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DAI");
 
68
    dpd_file2_scm(&D,0.0);
 
69
    dpd_file2_close(&D);
 
70
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dai");
 
71
    dpd_file2_scm(&D,0.0);
 
72
    dpd_file2_close(&D);
 
73
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DIA");
 
74
    dpd_file2_scm(&D,0.0);
 
75
    dpd_file2_close(&D);
 
76
    dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dia");
 
77
    dpd_file2_scm(&D,0.0);
 
78
  }
 
79
  
 
80
/*
 
81
  fprintf(outfile, "\n\tTrace of OPDM(2)_IJ     = %20.15f\n", trace_IJ);
 
82
  fprintf(outfile, "\n\tTrace of OPDM(2)_AB     = %20.15f\n", trace_AB);
 
83
*/
 
84
  fprintf(outfile, "\tTrace of OPDM(2)        = %20.15f\n", trace_IJ+trace_AB);
 
85
  fflush(outfile);
 
86
 
 
87
  return;
 
88
}
 
89
 
 
90
void uhf_opdm(void)
 
91
{
 
92
  int h, i, j, a, b;
 
93
  int I, J, A, B;
 
94
  double traceA=0.0;
 
95
  double traceB=0.0;
 
96
  dpdfile2 D;
 
97
  dpdfile2 T1A, T1B;
 
98
  dpdbuf4 T2A, T2B;
 
99
 
 
100
  if(params.semicanonical) {
 
101
    dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
 
102
    dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
 
103
    dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tIA");
 
104
    dpd_contract222(&T1A, &T1B, &D, 0, 0, -1, 0);
 
105
    dpd_file2_close(&T1A);
 
106
    dpd_file2_close(&T1B);
 
107
    
 
108
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB"); 
 
109
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB"); 
 
110
    dpd_contract442(&T2A, &T2B, &D, 0, 0, -1, 1);
 
111
    dpd_buf4_close(&T2A);
 
112
    dpd_buf4_close(&T2B);
 
113
  }
 
114
  else {        
 
115
    dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
 
116
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB"); 
 
117
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB"); 
 
118
    dpd_contract442(&T2A, &T2B, &D, 0, 0, -1, 0);
 
119
    dpd_buf4_close(&T2A);
 
120
    dpd_buf4_close(&T2B);
 
121
  }
 
122
  
 
123
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
124
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
125
  dpd_contract442(&T2A, &T2B, &D, 0, 0, -1, 1);
 
126
  traceA += dpd_file2_trace(&D);
 
127
  dpd_buf4_close(&T2A);
 
128
  dpd_buf4_close(&T2B);
 
129
  dpd_file2_close(&D);
 
130
    
 
131
  if(params.semicanonical) {
 
132
    dpd_file2_init(&D, CC_OEI, 0, 2, 2, "Dij");
 
133
    dpd_file2_init(&T1A, CC_OEI, 0, 2, 3, "tia");
 
134
    dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
 
135
    dpd_contract222(&T1A, &T1B, &D, 0, 0, -1, 0);
 
136
    dpd_file2_close(&T1A);
 
137
    dpd_file2_close(&T1B);
 
138
    
 
139
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 10, 17, 12, 17, 0, "tijab"); 
 
140
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 10, 17, 12, 17, 0, "tijab"); 
 
141
    dpd_contract442(&T2A, &T2B, &D, 0, 0, -1, 1);
 
142
    dpd_buf4_close(&T2A);
 
143
    dpd_buf4_close(&T2B);
 
144
  }
 
145
  else { 
 
146
    dpd_file2_init(&D, CC_OEI, 0, 2, 2, "Dij");
 
147
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 10, 17, 12, 17, 0, "tijab"); 
 
148
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 10, 17, 12, 17, 0, "tijab"); 
 
149
    dpd_contract442(&T2A, &T2B, &D, 0, 0, -1, 0);
 
150
    dpd_buf4_close(&T2A);
 
151
    dpd_buf4_close(&T2B);
 
152
  }
 
153
  
 
154
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
155
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
156
  dpd_contract442(&T2A, &T2B, &D, 1, 1, -1, 1);
 
157
  traceB += dpd_file2_trace(&D);
 
158
  dpd_buf4_close(&T2A);
 
159
  dpd_buf4_close(&T2B);
 
160
  dpd_file2_close(&D);
 
161
  
 
162
  if(params.semicanonical) {
 
163
    dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
 
164
    dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
 
165
    dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tIA");
 
166
    dpd_contract222(&T1A, &T1B, &D, 1, 1, 1, 0);
 
167
    dpd_file2_close(&T1A);
 
168
    dpd_file2_close(&T1B);
 
169
    
 
170
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB"); 
 
171
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB"); 
 
172
    dpd_contract442(&T2A, &T2B, &D, 3, 3, 1, 1);
 
173
    dpd_buf4_close(&T2A);
 
174
    dpd_buf4_close(&T2B);
 
175
  }
 
176
  else {
 
177
    dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
 
178
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB"); 
 
179
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB"); 
 
180
    dpd_contract442(&T2A, &T2B, &D, 3, 3, 1, 0);
 
181
    dpd_buf4_close(&T2A);
 
182
    dpd_buf4_close(&T2B);
 
183
  }
 
184
  
 
185
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
186
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
187
  dpd_contract442(&T2A, &T2B, &D, 2, 2, 1, 1);
 
188
  traceA += dpd_file2_trace(&D);
 
189
  dpd_buf4_close(&T2A);
 
190
  dpd_buf4_close(&T2B);
 
191
  dpd_file2_close(&D);
 
192
    
 
193
  if(params.semicanonical) {
 
194
    dpd_file2_init(&D, CC_OEI, 0, 3, 3, "Dab");
 
195
    dpd_file2_init(&T1A, CC_OEI, 0, 2, 3, "tia");
 
196
    dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
 
197
    dpd_contract222(&T1A, &T1B, &D, 1, 1, 1, 0);
 
198
    dpd_file2_close(&T1A);
 
199
    dpd_file2_close(&T1B);
 
200
    
 
201
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 12, 15, 12, 17, 0, "tijab"); 
 
202
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 12, 15, 12, 17, 0, "tijab"); 
 
203
    dpd_contract442(&T2A, &T2B, &D, 3, 3, 1, 1);
 
204
    dpd_buf4_close(&T2A);
 
205
    dpd_buf4_close(&T2B);
 
206
  }
 
207
  else {
 
208
    dpd_file2_init(&D, CC_OEI, 0, 3, 3, "Dab");
 
209
    dpd_buf4_init(&T2A, CC_TAMPS, 0, 12, 15, 12, 17, 0, "tijab"); 
 
210
    dpd_buf4_init(&T2B, CC_TAMPS, 0, 12, 15, 12, 17, 0, "tijab"); 
 
211
    dpd_contract442(&T2A, &T2B, &D, 3, 3, 1, 0);
 
212
    dpd_buf4_close(&T2A);
 
213
    dpd_buf4_close(&T2B);
 
214
  }
 
215
  
 
216
  dpd_buf4_init(&T2A, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
217
  dpd_buf4_init(&T2B, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
 
218
  dpd_contract442(&T2A, &T2B, &D, 3, 3, 1, 1);
 
219
  traceB += dpd_file2_trace(&D);
 
220
  dpd_buf4_close(&T2A);
 
221
  dpd_buf4_close(&T2B);
 
222
  dpd_file2_close(&D);
 
223
   
 
224
  /*
 
225
  fprintf(outfile,"\n");
 
226
  fprintf(outfile,"\tTrace of Alpha OPDM(2)  = %20.15f\n", fabs(traceA));
 
227
  fprintf(outfile,"\tTrace of Beta OPDM(2)   = %20.15f\n", fabs(traceB));
 
228
  */
 
229
    
 
230
  return;
 
231
}
 
232
 
 
233
 
 
234
/* This code isn't used anymore, and if you try to resurrect it, be sure to
 
235
note that the contract244 calls that point to the same file4 are not safe
 
236
in the dpd_cache structure.  -TDC, 11/21/08
 
237
*/
 
238
void rhf_sf_opdm(void) 
 
239
{
 
240
  dpdfile2 D;
 
241
  dpdbuf4 TA, TB;
 
242
  double trace_IJ = 0.0; 
 
243
  double trace_ij = 0.0;
 
244
  double trace_AB = 0.0; 
 
245
  double trace_ab = 0.0;
 
246
 
 
247
  dpd_buf4_init(&TA, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tIJAB");
 
248
  dpd_buf4_copy(&TA, CC_TMP0, "tIJAB");
 
249
  dpd_buf4_close(&TA);
 
250
  dpd_buf4_init(&TA, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tijab");
 
251
  dpd_buf4_copy(&TA, CC_TMP0, "tijab");
 
252
  dpd_buf4_close(&TA);
 
253
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
254
  dpd_buf4_copy(&TA, CC_TMP0, "tIjAb");
 
255
  dpd_buf4_close(&TA);
 
256
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tiJaB");
 
257
  dpd_buf4_copy(&TA, CC_TMP0, "tiJaB");
 
258
  dpd_buf4_close(&TA);
 
259
 
 
260
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
 
261
  dpd_file2_scm(&D, 0.0);
 
262
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB");
 
263
  dpd_buf4_init(&TB, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tIJAB");
 
264
  dpd_contract442(&TA, &TB, &D, 0, 0, -1.0, 0.0);
 
265
  dpd_buf4_close(&TA);
 
266
  dpd_buf4_close(&TB);
 
267
 
 
268
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
269
  dpd_buf4_init(&TB, CC_TMP0, 0, 0, 5, 0, 5, 0, "tIjAb");
 
270
  dpd_contract442(&TA, &TB, &D, 0, 0, -1.0, 1.0);
 
271
  dpd_buf4_close(&TA);
 
272
  dpd_buf4_close(&TB);
 
273
  trace_IJ += dpd_file2_trace(&D);
 
274
  dpd_file2_close(&D);
 
275
 
 
276
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "Dij");
 
277
  dpd_file2_scm(&D, 0.0);
 
278
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tijab");
 
279
  dpd_buf4_init(&TB, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tijab");
 
280
  dpd_contract442(&TA, &TB, &D, 0, 0, -1.0, 0.0);
 
281
  dpd_buf4_close(&TB);
 
282
  dpd_buf4_close(&TA);
 
283
 
 
284
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tiJaB");
 
285
  dpd_buf4_init(&TB, CC_TMP0, 0, 0, 5, 0, 5, 0, "tiJaB");
 
286
  dpd_contract442(&TA, &TB, &D, 0, 0, -1.0, 1.0);
 
287
  dpd_buf4_close(&TB);
 
288
  dpd_buf4_close(&TA);
 
289
  trace_ij += dpd_file2_trace(&D);
 
290
  dpd_file2_close(&D);
 
291
 
 
292
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
 
293
  dpd_file2_scm(&D, 0.0);
 
294
  dpd_buf4_init(&TA, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB");
 
295
  dpd_buf4_init(&TB, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tIJAB");
 
296
  dpd_contract442(&TA, &TB, &D, 3, 3, 1.0, 0.0);
 
297
  dpd_buf4_close(&TA);
 
298
  dpd_buf4_close(&TB);
 
299
 
 
300
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tiJaB");
 
301
  dpd_buf4_init(&TB, CC_TMP0, 0, 0, 5, 0, 5, 0, "tiJaB");
 
302
  dpd_contract442(&TA, &TB, &D, 3, 3, 1.0, 1.0);
 
303
  dpd_buf4_close(&TA);
 
304
  dpd_buf4_close(&TB);
 
305
  trace_AB += dpd_file2_trace(&D);
 
306
  dpd_file2_close(&D);
 
307
 
 
308
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "Dab");
 
309
  dpd_file2_scm(&D, 0.0);
 
310
  dpd_buf4_init(&TA, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tijab");
 
311
  dpd_buf4_init(&TB, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tijab");
 
312
  dpd_contract442(&TA, &TB, &D, 3, 3, 1.0, 0.0);
 
313
  dpd_buf4_close(&TB);
 
314
  dpd_buf4_close(&TA);
 
315
 
 
316
  dpd_buf4_init(&TA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
317
  dpd_buf4_init(&TB, CC_TMP0, 0, 0, 5, 0, 5, 0, "tIjAb");
 
318
  dpd_contract442(&TA, &TB, &D, 3, 3, 1.0, 1.0);
 
319
  dpd_buf4_close(&TB);
 
320
  dpd_buf4_close(&TA);
 
321
  trace_ab += dpd_file2_trace(&D);
 
322
  dpd_file2_close(&D);
 
323
 
 
324
  /* zero out Dia DIA Dai DAI */
 
325
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DAI");
 
326
  dpd_file2_scm(&D,0.0);
 
327
  dpd_file2_close(&D);
 
328
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dai");
 
329
  dpd_file2_scm(&D,0.0);
 
330
  dpd_file2_close(&D);
 
331
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DIA");
 
332
  dpd_file2_scm(&D,0.0);
 
333
  dpd_file2_close(&D);
 
334
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dia");
 
335
  dpd_file2_scm(&D,0.0);
 
336
  dpd_file2_close(&D);
 
337
  
 
338
 
 
339
  fprintf(outfile, "\n\tTrace of IJ onepdm = %20.15f\n", trace_IJ);
 
340
  fprintf(outfile, "\tTrace of ij onepdm = %20.15f\n", trace_ij);
 
341
  fprintf(outfile, "\tTrace of oo onepdm = %20.15f\n", trace_IJ+trace_ij);
 
342
  fprintf(outfile, "\tTrace of AB onepdm = %20.15f\n", trace_AB);
 
343
  fprintf(outfile, "\tTrace of ab onepdm = %20.15f\n", trace_ab);
 
344
  fprintf(outfile, "\tTrace of vv onepdm = %20.15f\n", (trace_AB+trace_ab));
 
345
  fprintf(outfile, "\tTrace of total onepdm = %20.15f\n", trace_IJ+trace_ij+trace_AB+trace_ab);
 
346
}
 
347
 
 
348
void uhf_sf_opdm(void)
 
349
{
 
350
  fprintf(outfile,"\n\tNot yet yo! -MLA\n");
 
351
  exit(PSI_RETURN_FAILURE);
 
352
}
 
353
 
 
354
}} /* End namespaces */