← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cceom/cc2_sigma.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cceom/cc2_sigma.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <math.h>
 
4
#include <libqt/qt.h>
 
5
#define EXTERN
 
6
#include "globals.h"
 
7
 
 
8
void cc2_sigmaSS(int i, int C_irr);
 
9
 
 
10
void cc2_sigma(int i, int C_irr) 
 
11
{
 
12
  dpdfile2 SIA;
 
13
  dpdbuf4 SIjAb;
 
14
  dpdfile2 CME;
 
15
  dpdbuf4 CMnEf;
 
16
  dpdfile2 FAE;
 
17
  dpdfile2 FMI;
 
18
  dpdfile2 FME;
 
19
  dpdbuf4 W;
 
20
  dpdbuf4 C;
 
21
  dpdfile2 S;
 
22
  dpdbuf4 WMbEj;
 
23
  dpdbuf4 WAmEf;
 
24
  dpdbuf4 WMnIe;
 
25
  dpdbuf4 WAbEi;
 
26
  dpdbuf4 WMbIj;
 
27
  dpdbuf4 Z;
 
28
  dpdbuf4 Z2;
 
29
  char lbl[32];
 
30
  char CME_lbl[32];
 
31
  char CMnEf_lbl[32];
 
32
  char SIjAb_lbl[32];
 
33
  int Gej, Gab, Gij, Gj, Gi, Ge, nrows, length, E, e, I;
 
34
  int Gam, Gef, Gim, Ga, Gm, ncols, A, a, am;
 
35
 
 
36
  if (params.eom_ref == 0) { /* RHF */
 
37
 
 
38
    cc2_sigmaSS(i,C_irr);
 
39
 
 
40
    sprintf(lbl, "%s %d", "SIA", i);
 
41
    dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
42
    dpd_file2_init(&FME, CC_OEI, H_IRR, 0, 1, "FME");
 
43
    dpd_buf4_init(&CMnEf, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "2CMnEf - CMnfE");
 
44
    dpd_dot24(&FME,&CMnEf,&SIA, 0, 0, 1.0, 1.0);
 
45
    dpd_buf4_close(&CMnEf);
 
46
    dpd_file2_close(&FME);
 
47
    dpd_file2_close(&SIA);
 
48
 
 
49
    /*
 
50
    sprintf(lbl, "%s %d", "SIA", i);
 
51
    dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
52
    sprintf(lbl, "%s %d", "CMnEf", i);
 
53
    dpd_buf4_init(&CMnEf, EOM_CMnEf, C_irr, 0, 5, 0, 5, 0, lbl);
 
54
    dpd_buf4_init(&WAmEf, CC_HBAR, H_IRR, 11, 5, 11, 5, 0, "WAmEf 2(Am,Ef) - (Am,fE)");
 
55
    dpd_contract442(&CMnEf, &WAmEf, &SIA, 0, 0, 1.0, 1.0);
 
56
    dpd_buf4_close(&WAmEf);
 
57
    dpd_buf4_close(&CMnEf);
 
58
    dpd_file2_close(&SIA);
 
59
    */
 
60
 
 
61
    dpd_buf4_init(&C, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "2CMnEf - CMnfE");
 
62
    dpd_buf4_init(&W, CC_HBAR, H_IRR, 11, 5, 11, 5, 0, "WAmEf");
 
63
    sprintf(lbl, "%s %d", "SIA", i);
 
64
    dpd_file2_init(&S, EOM_SIA, C_irr, 0, 1, lbl);
 
65
    dpd_file2_mat_init(&S);
 
66
    dpd_file2_mat_rd(&S);
 
67
    for(Gam=0; Gam < moinfo.nirreps; Gam++) {
 
68
      Gef = Gam ^ H_IRR;
 
69
      Gim = Gef ^ C_irr;
 
70
 
 
71
      dpd_buf4_mat_irrep_init(&C, Gim);
 
72
      dpd_buf4_mat_irrep_rd(&C, Gim);
 
73
      dpd_buf4_mat_irrep_shift13(&C, Gim);
 
74
 
 
75
      for(Gi=0; Gi < moinfo.nirreps; Gi++) {
 
76
        Ga = Gi ^ C_irr;
 
77
        Gm = Ga ^ Gam;
 
78
 
 
79
        W.matrix[Gam] = dpd_block_matrix(moinfo.occpi[Gm], W.params->coltot[Gef]);
 
80
 
 
81
        nrows = moinfo.occpi[Gi];
 
82
        ncols = moinfo.occpi[Gm] * W.params->coltot[Gef];
 
83
 
 
84
        for(A=0; A < moinfo.virtpi[Ga]; A++) {
 
85
          a = moinfo.vir_off[Ga] + A;
 
86
          am = W.row_offset[Gam][a];
 
87
 
 
88
          dpd_buf4_mat_irrep_rd_block(&W, Gam, am, moinfo.occpi[Gm]);
 
89
 
 
90
          if(nrows && ncols && moinfo.virtpi[Ga])
 
91
            C_DGEMV('n',nrows,ncols,1,C.shift.matrix[Gim][Gi][0],ncols,W.matrix[Gam][0], 1,
 
92
                    1, &(S.matrix[Gi][0][A]), moinfo.virtpi[Ga]);
 
93
        }
 
94
 
 
95
        dpd_free_block(W.matrix[Gam], moinfo.occpi[Gm], W.params->coltot[Gef]);
 
96
      }
 
97
 
 
98
      dpd_buf4_mat_irrep_close(&C, Gim);
 
99
    }
 
100
    dpd_file2_mat_wrt(&S);
 
101
    dpd_file2_mat_close(&S);
 
102
    dpd_file2_close(&S);
 
103
    dpd_buf4_close(&C);
 
104
    dpd_buf4_close(&W);
 
105
 
 
106
    sprintf(lbl, "%s %d", "SIA", i);
 
107
    dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
108
    sprintf(lbl, "%s %d", "CMnEf", i);
 
109
    dpd_buf4_init(&CMnEf, EOM_CMnEf, C_irr, 0, 5, 0, 5, 0, lbl);
 
110
    dpd_buf4_init(&WMnIe, CC_HBAR, H_IRR, 0, 11, 0, 11, 0, "WMnIe - 2WnMIe (Mn,eI)");
 
111
    dpd_contract442(&WMnIe, &CMnEf, &SIA, 3, 3, 1.0, 1.0);
 
112
    dpd_buf4_close(&CMnEf);
 
113
    dpd_buf4_close(&WMnIe);
 
114
    dpd_file2_close(&SIA);
 
115
 
 
116
    sprintf(CME_lbl, "%s %d", "CME", i);
 
117
    sprintf(SIjAb_lbl, "%s %d", "SIjAb", i);
 
118
 
 
119
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "WmaijDS Z(Ij,Ab)");
 
120
    dpd_buf4_init(&WMbIj, CC2_HET1, H_IRR, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
121
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
122
    dpd_contract244(&CME, &WMbIj, &Z, 0, 0, 1, 1.0, 0.0);
 
123
    dpd_file2_close(&CME);
 
124
    dpd_buf4_close(&WMbIj);
 
125
    dpd_buf4_sort(&Z, EOM_TMP, qpsr, 0, 5, "WmaijDS Z(jI,bA)");
 
126
    dpd_buf4_init(&SIjAb, EOM_SIjAb, C_irr, 0, 5, 0, 5, 0, SIjAb_lbl);
 
127
    dpd_buf4_axpy(&Z, &SIjAb,  -1.0);
 
128
    dpd_buf4_close(&Z);
 
129
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "WmaijDS Z(jI,bA)");
 
130
    dpd_buf4_axpy(&Z, &SIjAb,  -1.0);
 
131
    dpd_buf4_close(&Z);
 
132
    dpd_buf4_close(&SIjAb);
 
133
 
 
134
    sprintf(CME_lbl, "%s %d", "CME", i);
 
135
    sprintf(SIjAb_lbl, "%s %d", "SIjAb", i);
 
136
 
 
137
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "WabejDS Z(Ij,Ab)");
 
138
    dpd_buf4_scm(&Z,0);
 
139
    dpd_buf4_init(&W, CC2_HET1, H_IRR, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
140
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
141
    /*dpd_contract244(&CME, &WAbEi, &Z, 1, 0, 0, 1.0, 0.0);*/
 
142
    dpd_file2_mat_init(&CME);
 
143
    dpd_file2_mat_rd(&CME);
 
144
    for(Gej=0; Gej < moinfo.nirreps; Gej++) {
 
145
      Gab = Gej ^ H_IRR;
 
146
      Gij = Gab ^ C_irr;
 
147
 
 
148
      dpd_buf4_mat_irrep_init(&Z, Gij);
 
149
      dpd_buf4_mat_irrep_shift13(&Z, Gij);
 
150
 
 
151
      for(Ge=0; Ge < moinfo.nirreps; Ge++) {
 
152
        Gj = Ge ^ Gej;
 
153
        Gi = Gj ^ Gij;
 
154
 
 
155
        nrows = moinfo.occpi[Gj];
 
156
        length = nrows * W.params->coltot[Gab];
 
157
        dpd_buf4_mat_irrep_init_block(&W, Gej, nrows);
 
158
 
 
159
        for(E=0; E < moinfo.virtpi[Ge]; E++) {
 
160
          e = moinfo.vir_off[Ge] + E;
 
161
          dpd_buf4_mat_irrep_rd_block(&W, Gej, W.row_offset[Gej][e], nrows);
 
162
 
 
163
          for(I=0; I < moinfo.occpi[Gi]; I++) {
 
164
            if(length)
 
165
              C_DAXPY(length, CME.matrix[Gi][I][E], W.matrix[Gej][0], 1,
 
166
                      Z.shift.matrix[Gij][Gi][I], 1);
 
167
          }
 
168
        }
 
169
 
 
170
        dpd_buf4_mat_irrep_close_block(&W, Gej, nrows);
 
171
      }
 
172
 
 
173
      dpd_buf4_mat_irrep_wrt(&Z, Gij);
 
174
      dpd_buf4_mat_irrep_close(&Z, Gij);
 
175
 
 
176
    }
 
177
    dpd_file2_mat_close(&CME);
 
178
    dpd_file2_close(&CME);
 
179
    dpd_buf4_close(&W);
 
180
 
 
181
    /*
 
182
    dpd_buf4_sort(&Z, EOM_TMP, qpsr, 0, 5, "WabejDS Z(jI,bA)");
 
183
    dpd_buf4_init(&SIjAb, EOM_SIjAb, C_irr, 0, 5, 0, 5, 0, SIjAb_lbl);
 
184
    dpd_buf4_axpy(&Z, &SIjAb, 1.0);
 
185
    dpd_buf4_close(&Z);
 
186
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "WabejDS Z(jI,bA)");
 
187
    dpd_buf4_axpy(&Z, &SIjAb, 1.0);
 
188
    dpd_buf4_close(&Z);
 
189
    dpd_buf4_close(&SIjAb);
 
190
    */
 
191
 
 
192
    dpd_buf4_sort_axpy(&Z, EOM_SIjAb, qpsr, 0, 5, SIjAb_lbl, 1);
 
193
    dpd_buf4_init(&SIjAb, EOM_SIjAb, C_irr, 0, 5, 0, 5, 0, SIjAb_lbl);
 
194
    dpd_buf4_axpy(&Z, &SIjAb, 1.0);
 
195
    dpd_buf4_close(&SIjAb);
 
196
    dpd_buf4_close(&Z);
 
197
 
 
198
    sprintf(CMnEf_lbl, "%s %d", "CMnEf", i);
 
199
    sprintf(SIjAb_lbl, "%s %d", "SIjAb", i);
 
200
 
 
201
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "FDD_Fbe Z(Ij,Ab)");
 
202
    dpd_buf4_init(&CMnEf, EOM_CMnEf, C_irr, 0, 5, 0, 5, 0, CMnEf_lbl);
 
203
    dpd_file2_init(&FAE, CC_OEI, H_IRR, 1, 1, "fAB");
 
204
    dpd_contract424(&CMnEf, &FAE, &Z, 3, 1, 0, 1.0, 0.0);
 
205
    dpd_file2_close(&FAE);
 
206
    dpd_buf4_close(&CMnEf);
 
207
 
 
208
    dpd_buf4_sort(&Z, EOM_TMP, qpsr, 0, 5, "FDD_Fbe Z(jI,bA)");
 
209
    dpd_buf4_init(&Z2, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "FDD_Fbe Z(jI,bA)");
 
210
 
 
211
    dpd_buf4_init(&SIjAb, EOM_SIjAb, C_irr, 0, 5, 0, 5, 0, SIjAb_lbl);
 
212
    dpd_buf4_axpy(&Z, &SIjAb, 1.0);
 
213
    dpd_buf4_axpy(&Z2, &SIjAb, 1.0);
 
214
    dpd_buf4_close(&Z);
 
215
    dpd_buf4_close(&Z2);
 
216
    dpd_buf4_close(&SIjAb);
 
217
 
 
218
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "FDD_Fmj Z(Ij,Ab)");
 
219
    dpd_buf4_init(&CMnEf, EOM_CMnEf, C_irr, 0, 5, 0, 5, 0, CMnEf_lbl);
 
220
    dpd_file2_init(&FMI, CC_OEI, H_IRR, 0, 0, "fIJ");
 
221
    dpd_contract244(&FMI, &CMnEf, &Z, 0, 0, 0, 1.0, 0.0);
 
222
    dpd_file2_close(&FMI);
 
223
    dpd_buf4_close(&CMnEf);
 
224
    dpd_buf4_sort(&Z, EOM_TMP, qpsr, 0, 5, "FDD_Fmj Z(jI,bA)"); 
 
225
    dpd_buf4_init(&SIjAb, EOM_SIjAb, C_irr, 0, 5, 0, 5, 0, SIjAb_lbl);
 
226
    dpd_buf4_axpy(&Z, &SIjAb, -1.0);
 
227
    dpd_buf4_close(&Z);
 
228
    dpd_buf4_init(&Z, EOM_TMP, C_irr, 0, 5, 0, 5, 0, "FDD_Fmj Z(jI,bA)");
 
229
    dpd_buf4_axpy(&Z, &SIjAb, -1.0);
 
230
    dpd_buf4_close(&Z);
 
231
    dpd_buf4_close(&SIjAb);
 
232
  }
 
233
  else if (params.eom_ref == 1) { /* ROHF */
 
234
    printf("ROHF EOM_CC2 is not currently implemented\n");
 
235
    exit(PSI_RETURN_FAILURE);
 
236
  }
 
237
  else { /* UHF */
 
238
    printf("UHF EOM_CC2 is not currently implemented\n");
 
239
    exit(PSI_RETURN_FAILURE);
 
240
  }
 
241
}
 
242
 
 
243
void cc2_sigmaSS(int i, int C_irr) 
 
244
{
 
245
  dpdfile2 SIA;
 
246
  dpdfile2 CME;
 
247
  dpdfile2 FAE;
 
248
  dpdfile2 FMI;
 
249
  dpdbuf4 WMbEj;
 
250
  dpdfile2 Xme;
 
251
  dpdbuf4 T2;
 
252
  dpdbuf4 D;
 
253
  char lbl[32];
 
254
 
 
255
  if (params.eom_ref == 0) { /* RHF */
 
256
 
 
257
    /* sigmaSS */
 
258
 
 
259
    sprintf(lbl, "%s %d", "SIA", i);
 
260
    dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
261
    sprintf(lbl, "%s %d", "CME", i);
 
262
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
263
 
 
264
    dpd_file2_init(&FAE, CC_OEI, H_IRR, 1, 1, "FAE");
 
265
    dpd_contract222(&CME, &FAE, &SIA, 0, 0, 1.0, 0.0);
 
266
    dpd_file2_close(&FAE);
 
267
 
 
268
    dpd_file2_init(&FMI, CC_OEI, H_IRR, 0, 0, "FMI");
 
269
    dpd_contract222(&FMI, &CME, &SIA, 1, 1, -1.0, 1.0);
 
270
    dpd_file2_close(&FMI);
 
271
 
 
272
    dpd_buf4_init(&WMbEj, CC2_HET1, H_IRR, 10, 10, 10, 10, 0, "CC2 2 W(jb,ME) + W(Jb,Me)");
 
273
    dpd_contract422(&WMbEj, &CME, &SIA, 0, 0, 1.0, 1.0);
 
274
    dpd_buf4_close(&WMbEj);
 
275
 
 
276
    dpd_file2_init(&Xme, CC_OEI, C_irr, 0, 1, "XME");
 
277
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 10, 10, 10, 0, "D 2<ij|ab> - <ij|ba> (ia,jb)");
 
278
    dpd_contract422(&D, &CME, &Xme, 0, 0, 1, 0);
 
279
    dpd_buf4_close(&D);
 
280
 
 
281
    dpd_buf4_init(&T2, CC_TAMPS, 0, 10, 10, 10, 10, 0, "2 tIAjb - tIBja");
 
282
    dpd_contract422(&T2, &Xme, &SIA, 0, 0, 1, 1);
 
283
    dpd_buf4_close(&T2);
 
284
    dpd_file2_close(&Xme);
 
285
 
 
286
    dpd_file2_close(&CME);
 
287
    dpd_file2_close(&SIA);
 
288
  }
 
289
  else if (params.eom_ref == 1) { /* ROHF */
 
290
    printf("ROHF CC2-LR is not currently implemented\n");
 
291
    exit(PSI_RETURN_FAILURE);
 
292
  }
 
293
  else { /* UHF */
 
294
    printf("UHF CC2-LR is not currently implemented\n");
 
295
    exit(PSI_RETURN_FAILURE);
 
296
  }
 
297
}