← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cceom/rzero.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/cceom/rzero.cc
 
1
/*! \file
 
2
    \ingroup CCEOM
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cmath>
 
7
#include <cstring>
 
8
#include "MOInfo.h"
 
9
#include "Params.h"
 
10
#include "Local.h"
 
11
#define EXTERN
 
12
#include "globals.h"
 
13
 
 
14
namespace psi { namespace cceom {
 
15
 
 
16
extern double norm_C(dpdfile2 *CME, dpdfile2 *Cme,
 
17
            dpdbuf4 *CMNEF, dpdbuf4 *Cmnef, dpdbuf4 *CMnEf);
 
18
extern double dot_C(dpdfile2 *CME, dpdfile2 *Cme,
 
19
            dpdbuf4 *CMNEF, dpdbuf4 *Cmnef, dpdbuf4 *CMnEf);
 
20
extern double norm_C_rhf(dpdfile2 *CME, dpdbuf4 *CMnEf, dpdbuf4 *CMnfE);
 
21
extern void scm_C(dpdfile2 *CME, dpdfile2 *Cme, dpdbuf4 *CMNEF,
 
22
            dpdbuf4 *Cmnef, dpdbuf4 *CMnEf, double a);
 
23
 
 
24
/* this function determines R0, properly normalizes R, and checks orthogonality
 
25
 * with the ground state left eigenvector (1+lambda) */
 
26
 
 
27
/* for ROHF and UHF */
 
28
void rzero(int C_irr, int *converged) {
 
29
  double rzero=0.0, energy, norm, dotval;
 
30
  double dot_IA, dot_ia, dot_IJAB, dot_ijab, dot_IjAb;
 
31
  dpdfile2 RIA, Ria, RIA2, Ria2, FIA, Fia, LIA, Lia;
 
32
  dpdbuf4 RIJAB, Rijab, RIjAb, D, R2, LIJAB, Lijab, LIjAb;
 
33
  dpdbuf4 fRIJAB, fRijab, fRIjAb;
 
34
  int L_irr, i;
 
35
  int A_OCC, B_OCC, A_VIR, B_VIR;
 
36
  int AA_OCC, AA_VIR, BB_OCC, BB_VIR, AB_OCC, AB_VIR;
 
37
  char lbl[32], E_lbl[32], R1A_lbl[32], R1B_lbl[32];
 
38
  char R2AA_lbl[32], R2BB_lbl[32], R2AB_lbl[32];
 
39
  int R_index = -1;
 
40
 
 
41
  A_OCC = 0; A_VIR = 1;
 
42
  AA_OCC = 2; AA_VIR = 7;
 
43
  if (params.eom_ref <= 1) {
 
44
    B_OCC = 0; B_VIR = 1;
 
45
    BB_OCC = 2; BB_VIR = 7;
 
46
    AB_OCC = 0; AB_VIR = 5;
 
47
  }
 
48
  else if (params.eom_ref == 2) {
 
49
    B_OCC = 2; B_VIR = 3;
 
50
    BB_OCC = 12; BB_VIR = 17;
 
51
    AB_OCC = 22; AB_VIR = 28;
 
52
  }
 
53
  L_irr = eom_params.L_irr;
 
54
 
 
55
  for(i=0; i < eom_params.cs_per_irrep[C_irr]; i++) {
 
56
    if (!converged[i]) continue; /* this root did not converged */
 
57
    ++R_index;
 
58
 
 
59
    if(!strcmp(params.wfn,"EOM_CC2")) {
 
60
      sprintf(E_lbl, "EOM CC2 Energy for root %d %d", C_irr, R_index);
 
61
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) {
 
62
        fprintf(outfile,"No EOM CC2 Energy found in CC_INFO.  Not normalizing R.\n");
 
63
        return;
 
64
      }
 
65
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
66
    }
 
67
    else if(!strcmp(params.wfn,"EOM_CCSD")) {
 
68
      sprintf(E_lbl, "EOM CCSD Energy for root %d %d", C_irr, R_index);
 
69
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) {
 
70
        fprintf(outfile,"No EOM CCSD Energy found in CC_INFO.  Not normalizing R.\n");
 
71
        return;
 
72
      }
 
73
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
74
    }
 
75
    else if(!strcmp(params.wfn,"EOM_CC3")) {
 
76
      sprintf(E_lbl, "EOM CC3 Energy for root %d %d", C_irr, R_index);
 
77
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) {
 
78
        fprintf(outfile,"No EOM CC3 Energy found in CC_INFO.  Not normalizing R.\n");
 
79
        return;
 
80
      }
 
81
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
82
    }
 
83
 
 
84
    sprintf(R1A_lbl, "RIA %d %d", C_irr, R_index);
 
85
    sprintf(R1B_lbl, "Ria %d %d", C_irr, R_index);
 
86
    sprintf(R2AA_lbl, "RIJAB %d %d", C_irr, R_index);
 
87
    sprintf(R2BB_lbl, "Rijab %d %d", C_irr, R_index);
 
88
    sprintf(R2AB_lbl, "RIjAb %d %d", C_irr, R_index);
 
89
 
 
90
    /* Calculate <0| Hbar R |0> */
 
91
    if (C_irr == H_IRR) {
 
92
      dpd_file2_init(&FIA, CC_OEI, H_IRR, A_OCC, A_VIR, "FME");
 
93
      dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
 
94
      dot_IA = dpd_file2_dot(&FIA, &RIA);
 
95
      dpd_file2_close(&RIA);
 
96
      dpd_file2_close(&FIA);
 
97
  
 
98
      dpd_file2_init(&Fia, CC_OEI, H_IRR, B_OCC, B_VIR, "Fme");
 
99
      dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
 
100
      dot_ia = dpd_file2_dot(&Fia, &Ria);
 
101
      dpd_file2_close(&Ria);
 
102
      dpd_file2_close(&Fia);
 
103
 
 
104
      if (params.eom_ref == 1) {
 
105
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 2, 7, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
 
106
        dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
 
107
        dot_IJAB = dpd_buf4_dot(&D, &RIJAB);
 
108
        dpd_buf4_close(&RIJAB);
 
109
        dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
 
110
        dot_ijab = dpd_buf4_dot(&D, &Rijab);
 
111
        dpd_buf4_close(&Rijab);
 
112
        dpd_buf4_close(&D);
 
113
    
 
114
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 0, 5, 0, 5, 0, "D <ij|ab>");
 
115
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
116
        dot_IjAb = dpd_buf4_dot(&D, &RIjAb);
 
117
        dpd_buf4_close(&RIjAb);
 
118
        dpd_buf4_close(&D);
 
119
      }
 
120
      else if (params.eom_ref == 2) {
 
121
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 2, 7, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
 
122
        dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
 
123
        dot_IJAB = dpd_buf4_dot(&D, &RIJAB);
 
124
        dpd_buf4_close(&RIJAB);
 
125
        dpd_buf4_close(&D);
 
126
  
 
127
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 12, 17, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
 
128
        dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 12, 17, 12, 17, 0, R2BB_lbl);
 
129
        dot_ijab = dpd_buf4_dot(&D, &Rijab);
 
130
        dpd_buf4_close(&Rijab);
 
131
        dpd_buf4_close(&D);
 
132
                           
 
133
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
134
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 22, 28, 22, 28, 0, R2AB_lbl);
 
135
        dot_IjAb = dpd_buf4_dot(&D, &RIjAb);
 
136
        dpd_buf4_close(&RIjAb);
 
137
        dpd_buf4_close(&D);
 
138
      }
 
139
      rzero = (dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb)/energy;
 
140
    }
 
141
    else { /* C and H are different irreps */
 
142
      rzero = 0.0;
 
143
    }
 
144
 
 
145
    /* Now normalize so that <R|R> = 1 */
 
146
    dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
 
147
    dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
 
148
    dpd_buf4_init(&fRIJAB, CC_RAMPS, C_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, R2AA_lbl);
 
149
    dpd_buf4_init(&fRijab, CC_RAMPS, C_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, R2BB_lbl);
 
150
    dpd_buf4_init(&fRIjAb, CC_RAMPS, C_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, R2AB_lbl);
 
151
 
 
152
    /* make R0 a positive number */
 
153
    /*
 
154
    if (rzero < 0.0) {
 
155
      rzero *= -1.0;
 
156
      dpd_file2_scm(&RIA,-1.0);
 
157
      dpd_file2_scm(&Ria,-1.0);
 
158
      dpd_buf4_scm(&fRIJAB,-1.0);
 
159
      dpd_buf4_scm(&fRijab,-1.0);
 
160
      dpd_buf4_scm(&fRIjAb,-1.0);
 
161
    }
 
162
    */
 
163
 
 
164
    /*
 
165
    norm = norm_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb);
 
166
    norm *= norm;
 
167
    */
 
168
    norm = dot_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb);
 
169
    norm += rzero * rzero;
 
170
    norm = sqrt(norm);
 
171
    rzero = rzero / norm;
 
172
    scm_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb, 1.0/norm);
 
173
 
 
174
    norm = dot_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb);
 
175
    norm += rzero * rzero;
 
176
    fprintf(outfile,"<R|R> = %20.16lf\n",norm);
 
177
 
 
178
/* just debugging with converged solutions - also my need a sort_C() */
 
179
/*
 
180
dpd_file2_copy(&RIA, EOM_CME, "CME 0");
 
181
dpd_file2_copy(&Ria, EOM_Cme, "Cme 0");
 
182
dpd_buf4_copy(&fRIJAB, EOM_CMNEF, "CMNEF 0");
 
183
dpd_buf4_copy(&fRijab, EOM_Cmnef, "Cmnef 0");
 
184
dpd_buf4_copy(&fRIjAb, EOM_CMnEf, "CMnEf 0");
 
185
*/
 
186
/* end debugging stuff */ 
 
187
 
 
188
 
 
189
    dpd_file2_close(&RIA);
 
190
    dpd_file2_close(&Ria);
 
191
    dpd_buf4_close(&fRIJAB);
 
192
    dpd_buf4_close(&fRijab);
 
193
    dpd_buf4_close(&fRIjAb);
 
194
 
 
195
    if(!strcmp(params.wfn,"EOM_CC2")) {
 
196
      fprintf(outfile,"EOM CC2 R0 for root %d = %15.11lf\n", R_index, rzero);
 
197
      sprintf(lbl, "EOM CC2 R0 for root %d %d", C_irr, R_index);
 
198
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
199
    }
 
200
    else if(!strcmp(params.wfn,"EOM_CCSD")) {
 
201
      fprintf(outfile,"EOM CCSD R0 for root %d = %15.11lf\n", R_index, rzero);
 
202
      sprintf(lbl, "EOM CCSD R0 for root %d %d", C_irr, R_index);
 
203
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
204
    }
 
205
    else if(!strcmp(params.wfn,"EOM_CC3")) {
 
206
      fprintf(outfile,"EOM CC3 R0 for root %d = %15.11lf\n", R_index, rzero);
 
207
      sprintf(lbl, "EOM CC3 R0 for root %d %d", C_irr, R_index);
 
208
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
209
    }
 
210
 
 
211
    if (eom_params.dot_with_L) {
 
212
      /* evaluate check <R|L> == 0 */
 
213
      if (C_irr == L_irr ) {
 
214
        dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
 
215
        dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
 
216
        dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, R2AA_lbl);
 
217
        dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, R2BB_lbl);
 
218
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, R2AB_lbl);
 
219
  
 
220
        dpd_file2_init(&LIA, CC_OEI, L_irr, A_OCC, A_VIR, "LIA");
 
221
        dpd_file2_init(&Lia, CC_OEI, L_irr, B_OCC, B_VIR, "Lia");
 
222
        dpd_buf4_init(&LIJAB, CC_LAMPS, L_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, "LIJAB");
 
223
        dpd_buf4_init(&Lijab, CC_LAMPS, L_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, "Lijab");
 
224
        dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, "LIjAb");
 
225
  
 
226
        fprintf(outfile,"\nROHF orthogonality test\n");
 
227
        fprintf(outfile,"<L0|R0>            = %15.10lf\n", eom_params.L0 * rzero);
 
228
        dot_IA = dpd_file2_dot(&LIA, &RIA);
 
229
        fprintf(outfile,"<LIA|RIA>          = %15.10lf\n", dot_IA);
 
230
        dot_ia = dpd_file2_dot(&Lia, &Ria);
 
231
        fprintf(outfile,"<Lia|Ria>          = %15.10lf\n", dot_ia);
 
232
        dot_IJAB = dpd_buf4_dot(&LIJAB, &RIJAB);
 
233
        fprintf(outfile,"<LIJAB|RIJAB>      = %15.10lf\n", dot_IJAB);
 
234
        dot_ijab = dpd_buf4_dot(&Lijab, &Rijab);
 
235
        fprintf(outfile,"<Lijab|Rijab>      = %15.10lf\n", dot_ijab);
 
236
        dot_IjAb = dpd_buf4_dot(&LIjAb, &RIjAb);
 
237
        fprintf(outfile,"<LIjAb|RIjAb>      = %15.10lf\n", dot_IjAb);
 
238
        fprintf(outfile,"<L|R>              = %15.10lf\n", (eom_params.L0 * rzero)
 
239
            + dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb);
 
240
  
 
241
        dpd_file2_close(&LIA);
 
242
        dpd_file2_close(&Lia);
 
243
        dpd_buf4_close(&LIJAB);
 
244
        dpd_buf4_close(&Lijab);
 
245
        dpd_buf4_close(&LIjAb);
 
246
        dpd_file2_close(&RIA);
 
247
        dpd_file2_close(&Ria);
 
248
        dpd_buf4_close(&RIJAB);
 
249
        dpd_buf4_close(&Rijab);
 
250
        dpd_buf4_close(&RIjAb);
 
251
      }
 
252
      else {
 
253
        fprintf(outfile,"\nOverlap <R|L> zero by symmetry\n");
 
254
      }
 
255
    } /* end dot_with_L loop */
 
256
  }
 
257
  return;
 
258
}
 
259
 
 
260
/* normalizes R and produces copies of R that are ROHF-like */
 
261
/* sort_amps then produces the others sorted versions of R */
 
262
 
 
263
void rzero_rhf(int C_irr, int *converged) {
 
264
  double r1, r2, rzero=0.0, energy, norm, dotval;
 
265
  double dot_IA, dot_ia, dot_IJAB, dot_ijab, dot_IjAb;
 
266
  dpdfile2 RIA, FIA, LIA, Lia, Ria;
 
267
  dpdbuf4 RIjAb, RIjbA, RIjAb1, RIjbA1, D, R2, LIjAb, RIJAB, Rijab;
 
268
  dpdbuf4 LIJAB, Lijab;
 
269
  int L_irr,i;
 
270
  char lbl[32], E_lbl[32], R1A_lbl[32], R1B_lbl[32], *blank;
 
271
  char R2AA_lbl[32], R2BB_lbl[32], R2AB_lbl[32];
 
272
  int R_index=-1;
 
273
 
 
274
  L_irr = eom_params.L_irr;
 
275
 
 
276
  for(i=0; i < eom_params.cs_per_irrep[C_irr]; i++) {
 
277
    if (!converged[i]) continue; /* this root did not converge */
 
278
    ++R_index;
 
279
 
 
280
    if(!strcmp(params.wfn,"EOM_CC2")) {
 
281
      sprintf(E_lbl, "EOM CC2 Energy for root %d %d", C_irr, R_index);
 
282
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) { 
 
283
        fprintf(outfile,"No EOM CC2 Energy found in CC_INFO.  Not normalizing R.\n");
 
284
        return;
 
285
      }
 
286
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
287
    }
 
288
    else if(!strcmp(params.wfn,"EOM_CCSD")) {
 
289
      sprintf(E_lbl, "EOM CCSD Energy for root %d %d", C_irr, R_index);
 
290
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) { 
 
291
        fprintf(outfile,"No EOM CCSD Energy found in CC_INFO.  Not normalizing R.\n");
 
292
        return;
 
293
      }
 
294
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
295
    }
 
296
    else if(!strcmp(params.wfn,"EOM_CC3")) {
 
297
      sprintf(E_lbl, "EOM CC3 Energy for root %d %d", C_irr, R_index);
 
298
      if(psio_tocscan(CC_INFO, E_lbl) == NULL) { 
 
299
        fprintf(outfile,"No EOM CC3 Energy found in CC_INFO.  Not normalizing R.\n");
 
300
        return;
 
301
      }
 
302
      psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
 
303
    }
 
304
 
 
305
    sprintf(R1A_lbl, "RIA %d %d", C_irr, R_index);
 
306
    sprintf(R1B_lbl, "Ria %d %d", C_irr, R_index);
 
307
    sprintf(R2AB_lbl, "RIjAb %d %d", C_irr, R_index);
 
308
    sprintf(R2AA_lbl, "RIJAB %d %d", C_irr, R_index);
 
309
    sprintf(R2BB_lbl, "Rijab %d %d", C_irr, R_index);
 
310
 
 
311
    /* produce RIjbA and 2RIjAb-RIjbA copies - not yet normalized */
 
312
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
313
    dpd_buf4_sort(&RIjAb, CC_TMP, pqsr, 0, 5, "RIjbA");
 
314
    sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
315
    dpd_buf4_copy(&RIjAb, CC_RAMPS, lbl);
 
316
    dpd_buf4_close(&RIjAb);
 
317
 
 
318
    sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
319
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
 
320
    dpd_buf4_scm(&RIjAb, 2.0);
 
321
    dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
 
322
    dpd_buf4_axpy(&RIjbA, &RIjAb, -1.0);
 
323
    dpd_buf4_close(&RIjbA);
 
324
    dpd_buf4_close(&RIjAb);
 
325
 
 
326
    /* calculate <0| hbar | 0> */
 
327
    if (!params.full_matrix) {
 
328
      if (C_irr == H_IRR) {
 
329
        dpd_file2_init(&FIA, CC_OEI, H_IRR, 0, 1, "FME");
 
330
        dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
331
        r1 = 2.0 * dpd_file2_dot(&FIA, &RIA);
 
332
        dpd_file2_close(&RIA);
 
333
        dpd_file2_close(&FIA);
 
334
    
 
335
        sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
336
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
 
337
        dpd_buf4_init(&D, CC_DINTS, H_IRR, 0, 5, 0, 5, 0, "D <ij|ab>");
 
338
        r2 = dpd_buf4_dot(&D, &RIjAb);
 
339
        dpd_buf4_close(&D);
 
340
        dpd_buf4_close(&RIjAb);
 
341
        rzero = (r1 + r2)/energy;
 
342
      }
 
343
      else {
 
344
        rzero = 0.0;
 
345
      }
 
346
    }
 
347
    else { /* full matrix */
 
348
      sprintf(lbl, "%s %d %d", "R0", C_irr, R_index);
 
349
      psio_read_entry(CC_RAMPS, lbl, (char *) &rzero, sizeof(double));
 
350
    }
 
351
 
 
352
    /* Now normalize R so that <R|R> = 1 */
 
353
    dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
354
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
355
    dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
 
356
 
 
357
    /*
 
358
    if (rzero < 0.0) {
 
359
      rzero *= -1.0;
 
360
      dpd_file2_scm(&RIA,-1.0);
 
361
      dpd_buf4_scm(&RIjAb,-1.0);
 
362
      dpd_buf4_scm(&RIjbA,-1.0);
 
363
    }
 
364
    */
 
365
 
 
366
    norm = norm_C_rhf(&RIA, &RIjAb, &RIjbA);
 
367
    norm *= norm;
 
368
    norm += rzero * rzero;
 
369
    norm = sqrt(norm);
 
370
    rzero = rzero / norm;
 
371
    dpd_file2_scm(&RIA, 1.0/norm);
 
372
    dpd_buf4_scm(&RIjAb, 1.0/norm);
 
373
    dpd_buf4_scm(&RIjbA, 1.0/norm);
 
374
 
 
375
/* just for debugging cc3 put normalized vector back into C as well */
 
376
/*
 
377
dpd_file2_copy(&RIA, EOM_CME, "CME 0");
 
378
dpd_buf4_copy(&RIjAb, EOM_CMnEf, "CMnEf 0");
 
379
*/
 
380
/* end debugging stuff */ 
 
381
 
 
382
    dpd_file2_close(&RIA);
 
383
    dpd_buf4_close(&RIjAb);
 
384
    dpd_buf4_close(&RIjbA);
 
385
 
 
386
    if(!strcmp(params.wfn,"EOM_CC2")) {
 
387
      fprintf(outfile,"EOM CC2 R0 for root %d = %15.11lf\n", R_index, rzero);
 
388
      sprintf(lbl, "EOM CC2 R0 for root %d %d", C_irr, R_index);
 
389
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
390
    }
 
391
    else if(!strcmp(params.wfn,"EOM_CCSD")) {
 
392
      fprintf(outfile,"EOM CCSD R0 for root %d = %15.11lf\n", R_index, rzero);
 
393
      sprintf(lbl, "EOM CCSD R0 for root %d %d", C_irr, R_index);
 
394
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
395
    }
 
396
    else if(!strcmp(params.wfn,"EOM_CC3")) {
 
397
      fprintf(outfile,"EOM CC3 R0 for root %d = %15.11lf\n", R_index, rzero);
 
398
      sprintf(lbl, "EOM CC3 R0 for root %d %d", C_irr, R_index);
 
399
      psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
 
400
    }
 
401
 
 
402
    /* produce ROHF like quantities and 2RIjAb-RIjbA */
 
403
    dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
404
    dpd_file2_copy(&RIA, CC_RAMPS, R1B_lbl);
 
405
    dpd_file2_close(&RIA);
 
406
 
 
407
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
408
    sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
409
    dpd_buf4_copy(&RIjAb, CC_RAMPS, lbl);
 
410
    dpd_buf4_close(&RIjAb);
 
411
 
 
412
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 2, 7, 0, 5, 1, R2AB_lbl);
 
413
    dpd_buf4_copy(&RIjAb, CC_RAMPS, R2AA_lbl);
 
414
    dpd_buf4_copy(&RIjAb, CC_RAMPS, R2BB_lbl);
 
415
    dpd_buf4_close(&RIjAb);
 
416
      
 
417
    dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
 
418
    sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
419
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
 
420
    dpd_buf4_scm(&RIjAb, 2.0);
 
421
    dpd_buf4_axpy(&RIjbA, &RIjAb, -1.0);
 
422
    dpd_buf4_close(&RIjAb);
 
423
    dpd_buf4_close(&RIjbA);
 
424
  
 
425
    /* test normalization with produced ROHF-like quantities */
 
426
 
 
427
    dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
428
    norm = dpd_file2_dot_self(&RIA);
 
429
    dpd_file2_close(&RIA);
 
430
    dpd_file2_init(&Ria, CC_RAMPS, C_irr, 0, 1, R1B_lbl);
 
431
    norm += dpd_file2_dot_self(&Ria);
 
432
    dpd_file2_close(&Ria);
 
433
    dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
 
434
    norm += dpd_buf4_dot_self(&RIJAB);
 
435
    dpd_buf4_close(&RIJAB);
 
436
    dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
 
437
    norm += dpd_buf4_dot_self(&Rijab);
 
438
    dpd_buf4_close(&Rijab);
 
439
    dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
440
    norm += dpd_buf4_dot_self(&RIjAb);
 
441
    dpd_buf4_close(&RIjAb);
 
442
    norm += rzero * rzero;
 
443
 
 
444
#ifdef EOM_DEBUG
 
445
    fprintf(outfile,"Norm with produced ROHF-like quantities = %15.10lf\n", norm);
 
446
    dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
447
    dpd_file2_print(&RIA, outfile);
 
448
    dpd_file2_close(&RIA);
 
449
#endif
 
450
 
 
451
    /* check orthogonality with L */
 
452
    if (eom_params.dot_with_L) {
 
453
      if (C_irr == L_irr) {
 
454
        dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
 
455
        dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
456
        r1 = 2.0 * dpd_file2_dot(&LIA, &RIA);
 
457
        dpd_file2_close(&RIA);
 
458
        dpd_file2_close(&LIA);
 
459
  
 
460
        dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
461
        sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
 
462
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
 
463
        r2 = dpd_buf4_dot(&LIjAb, &RIjAb);
 
464
        dpd_buf4_close(&RIjAb);
 
465
        dpd_buf4_close(&LIjAb);
 
466
  
 
467
        dotval = r1 + r2 + (eom_params.L0 * rzero);
 
468
        fprintf(outfile,"Performing RHF orthogonality test\n");
 
469
        fprintf(outfile,"<L0|R0>              = %15.10lf\n", eom_params.L0 * rzero);
 
470
        fprintf(outfile,"2*<LIA|RIA>          = %15.10lf\n", r1);
 
471
        fprintf(outfile,"<LIjAb|2RIjAb-RIjbA> = %15.10lf\n", r2);
 
472
        fprintf(outfile,"<L|R>                = %15.10lf\n", dotval);
 
473
      }
 
474
      else {
 
475
        fprintf(outfile,"<L|R> zero by symmetry\n");
 
476
      }
 
477
      if (C_irr == L_irr ) {
 
478
       /* double check orthogonality rohf-like */
 
479
        dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
 
480
        dpd_file2_init(&Ria, CC_RAMPS, C_irr, 0, 1, R1B_lbl);
 
481
        dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
 
482
        dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
 
483
        dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
 
484
  
 
485
        dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
 
486
        dpd_file2_init(&Lia, CC_OEI, L_irr, 0, 1, "Lia");
 
487
        dpd_buf4_init(&LIJAB, CC_LAMPS, L_irr, 2, 7, 2, 7, 0, "LIJAB");
 
488
        dpd_buf4_init(&Lijab, CC_LAMPS, L_irr, 2, 7, 2, 7, 0, "Lijab");
 
489
        dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
490
  
 
491
        dot_IA = dpd_file2_dot(&LIA, &RIA);
 
492
        dot_ia = dpd_file2_dot(&Lia, &Ria);
 
493
        dot_IJAB = dpd_buf4_dot(&LIJAB, &RIJAB);
 
494
        dot_ijab = dpd_buf4_dot(&Lijab, &Rijab);
 
495
        dot_IjAb = dpd_buf4_dot(&LIjAb, &RIjAb);
 
496
  
 
497
        dpd_file2_close(&RIA);
 
498
        dpd_file2_close(&Ria);
 
499
        dpd_buf4_close(&RIJAB);
 
500
        dpd_buf4_close(&Rijab);
 
501
        dpd_buf4_close(&RIjAb);
 
502
  
 
503
        dpd_file2_close(&LIA);
 
504
        dpd_file2_close(&Lia);
 
505
        dpd_buf4_close(&LIJAB);
 
506
        dpd_buf4_close(&Lijab);
 
507
        dpd_buf4_close(&LIjAb);
 
508
 
 
509
        fprintf(outfile,"\nROHF-like orthogonality test\n");
 
510
        fprintf(outfile,"<L0|R0>              = %15.10lf\n", eom_params.L0 * rzero);
 
511
        fprintf(outfile,"<LIA|RIA>            = %15.10lf\n", dot_IA);
 
512
        fprintf(outfile,"<Lia|Ria>            = %15.10lf\n", dot_ia);
 
513
        fprintf(outfile,"<LIJAB|RIJAB>        = %15.10lf\n", dot_IJAB);
 
514
        fprintf(outfile,"<Lijab|Rijab>        = %15.10lf\n", dot_ijab);
 
515
        fprintf(outfile,"<LIjAb|RIjAb>        = %15.10lf\n", dot_IjAb);
 
516
        fprintf(outfile,"<L|R>                = %15.10lf\n", eom_params.L0 * rzero +
 
517
        dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb);
 
518
      }
 
519
    } /* end dot with L, <L|R> overlap checks */
 
520
  } /* end loop over Cs */
 
521
  return;
 
522
}
 
523
 
 
524
}} // namespace psi::cceom