7
extern double norm_C(dpdfile2 *CME, dpdfile2 *Cme,
8
dpdbuf4 *CMNEF, dpdbuf4 *Cmnef, dpdbuf4 *CMnEf);
9
extern double norm_C_rhf(dpdfile2 *CME, dpdbuf4 *CMnEf, dpdbuf4 *CMnfE);
10
extern void scm_C(dpdfile2 *CME, dpdfile2 *Cme, dpdbuf4 *CMNEF,
11
dpdbuf4 *Cmnef, dpdbuf4 *CMnEf, double a);
13
/* this function determines R0, properly normalizes R, and checks orthogonality
14
* with the ground state left eigenvector (1+lambda) */
16
/* for ROHF and UHF */
17
void rzero(int C_irr, int *converged) {
18
double rzero=0.0, energy, norm, dotval;
19
double dot_IA, dot_ia, dot_IJAB, dot_ijab, dot_IjAb;
20
dpdfile2 RIA, Ria, RIA2, Ria2, FIA, Fia, LIA, Lia;
21
dpdbuf4 RIJAB, Rijab, RIjAb, D, R2, LIJAB, Lijab, LIjAb;
22
dpdbuf4 fRIJAB, fRijab, fRIjAb;
24
int A_OCC, B_OCC, A_VIR, B_VIR;
25
int AA_OCC, AA_VIR, BB_OCC, BB_VIR, AB_OCC, AB_VIR;
26
char lbl[32], E_lbl[32], R1A_lbl[32], R1B_lbl[32];
27
char R2AA_lbl[32], R2BB_lbl[32], R2AB_lbl[32];
31
AA_OCC = 2; AA_VIR = 7;
32
if (params.eom_ref <= 1) {
34
BB_OCC = 2; BB_VIR = 7;
35
AB_OCC = 0; AB_VIR = 5;
37
else if (params.eom_ref == 2) {
39
BB_OCC = 12; BB_VIR = 17;
40
AB_OCC = 22; AB_VIR = 28;
42
L_irr = eom_params.L_irr;
44
for(i=0; i < eom_params.cs_per_irrep[C_irr]; i++) {
45
if (!converged[i]) continue; /* this root did not converged */
48
sprintf(E_lbl, "EOM CCSD Energy for root %d %d", C_irr, R_index);
49
if ( psio_tocscan(CC_INFO, E_lbl) == NULL) {
50
fprintf(outfile,"No EOM CCSD Energy found in CC_INFO. Not normalizing R.\n");
53
psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
55
sprintf(R1A_lbl, "RIA %d %d", C_irr, R_index);
56
sprintf(R1B_lbl, "Ria %d %d", C_irr, R_index);
57
sprintf(R2AA_lbl, "RIJAB %d %d", C_irr, R_index);
58
sprintf(R2BB_lbl, "Rijab %d %d", C_irr, R_index);
59
sprintf(R2AB_lbl, "RIjAb %d %d", C_irr, R_index);
61
/* Calculate <0| Hbar R |0> */
63
dpd_file2_init(&FIA, CC_OEI, H_IRR, A_OCC, A_VIR, "FME");
64
dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
65
dot_IA = dpd_file2_dot(&FIA, &RIA);
66
dpd_file2_close(&RIA);
67
dpd_file2_close(&FIA);
69
dpd_file2_init(&Fia, CC_OEI, H_IRR, B_OCC, B_VIR, "Fme");
70
dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
71
dot_ia = dpd_file2_dot(&Fia, &Ria);
72
dpd_file2_close(&Ria);
73
dpd_file2_close(&Fia);
75
if (params.eom_ref == 1) {
76
dpd_buf4_init(&D, CC_DINTS, H_IRR, 2, 7, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
77
dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
78
dot_IJAB = dpd_buf4_dot(&D, &RIJAB);
79
dpd_buf4_close(&RIJAB);
80
dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
81
dot_ijab = dpd_buf4_dot(&D, &Rijab);
82
dpd_buf4_close(&Rijab);
85
dpd_buf4_init(&D, CC_DINTS, H_IRR, 0, 5, 0, 5, 0, "D <ij|ab>");
86
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
87
dot_IjAb = dpd_buf4_dot(&D, &RIjAb);
88
dpd_buf4_close(&RIjAb);
91
else if (params.eom_ref == 2) {
92
dpd_buf4_init(&D, CC_DINTS, H_IRR, 2, 7, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
93
dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
94
dot_IJAB = dpd_buf4_dot(&D, &RIJAB);
95
dpd_buf4_close(&RIJAB);
98
dpd_buf4_init(&D, CC_DINTS, H_IRR, 12, 17, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
99
dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 12, 17, 12, 17, 0, R2BB_lbl);
100
dot_ijab = dpd_buf4_dot(&D, &Rijab);
101
dpd_buf4_close(&Rijab);
104
dpd_buf4_init(&D, CC_DINTS, H_IRR, 22, 28, 22, 28, 0, "D <Ij|Ab>");
105
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 22, 28, 22, 28, 0, R2AB_lbl);
106
dot_IjAb = dpd_buf4_dot(&D, &RIjAb);
107
dpd_buf4_close(&RIjAb);
110
rzero = (dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb)/energy;
112
else { /* C and H are different irreps */
116
/* Now normalize so that <R|R> = 1 */
117
dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
118
dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
119
dpd_buf4_init(&fRIJAB, CC_RAMPS, C_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, R2AA_lbl);
120
dpd_buf4_init(&fRijab, CC_RAMPS, C_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, R2BB_lbl);
121
dpd_buf4_init(&fRIjAb, CC_RAMPS, C_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, R2AB_lbl);
123
norm = norm_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb);
125
norm += rzero * rzero;
127
rzero = rzero / norm;
128
scm_C(&RIA, &Ria, &fRIJAB, &fRijab, &fRIjAb, 1.0/norm);
131
dpd_file2_print(&RIA, outfile);
132
dpd_buf4_print(&fRIJAB, outfile, 1);
133
dpd_buf4_print(&fRijab, outfile, 1);
134
dpd_buf4_print(&fRIjAb, outfile, 1);
137
dpd_file2_close(&RIA);
138
dpd_file2_close(&Ria);
139
dpd_buf4_close(&fRIJAB);
140
dpd_buf4_close(&fRijab);
141
dpd_buf4_close(&fRIjAb);
143
fprintf(outfile,"EOM CCSD R0 for root %d = %15.10lf\n", R_index, rzero);
144
sprintf(lbl, "EOM CCSD R0 for root %d %d", C_irr, R_index);
145
psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
147
if (eom_params.dot_with_L) {
148
/* evaluate check <R|L> == 0 */
149
if (C_irr == L_irr ) {
150
dpd_file2_init(&RIA, CC_RAMPS, C_irr, A_OCC, A_VIR, R1A_lbl);
151
dpd_file2_init(&Ria, CC_RAMPS, C_irr, B_OCC, B_VIR, R1B_lbl);
152
dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, R2AA_lbl);
153
dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, R2BB_lbl);
154
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, R2AB_lbl);
156
dpd_file2_init(&LIA, CC_OEI, L_irr, A_OCC, A_VIR, "LIA");
157
dpd_file2_init(&Lia, CC_OEI, L_irr, B_OCC, B_VIR, "Lia");
158
dpd_buf4_init(&LIJAB, CC_LAMPS, L_irr, AA_OCC, AA_VIR, AA_OCC, AA_VIR, 0, "LIJAB");
159
dpd_buf4_init(&Lijab, CC_LAMPS, L_irr, BB_OCC, BB_VIR, BB_OCC, BB_VIR, 0, "Lijab");
160
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, AB_OCC, AB_VIR, AB_OCC, AB_VIR, 0, "LIjAb");
162
fprintf(outfile,"\nROHF orthogonality test\n");
163
fprintf(outfile,"<L0|R0> = %15.10lf\n", eom_params.L0 * rzero);
164
dot_IA = dpd_file2_dot(&LIA, &RIA);
165
fprintf(outfile,"<LIA|RIA> = %15.10lf\n", dot_IA);
166
dot_ia = dpd_file2_dot(&Lia, &Ria);
167
fprintf(outfile,"<Lia|Ria> = %15.10lf\n", dot_ia);
168
dot_IJAB = dpd_buf4_dot(&LIJAB, &RIJAB);
169
fprintf(outfile,"<LIJAB|RIJAB> = %15.10lf\n", dot_IJAB);
170
dot_ijab = dpd_buf4_dot(&Lijab, &Rijab);
171
fprintf(outfile,"<Lijab|Rijab> = %15.10lf\n", dot_ijab);
172
dot_IjAb = dpd_buf4_dot(&LIjAb, &RIjAb);
173
fprintf(outfile,"<LIjAb|RIjAb> = %15.10lf\n", dot_IjAb);
174
fprintf(outfile,"<L|R> = %15.10lf\n", (eom_params.L0 * rzero)
175
+ dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb);
177
dpd_file2_close(&LIA);
178
dpd_file2_close(&Lia);
179
dpd_buf4_close(&LIJAB);
180
dpd_buf4_close(&Lijab);
181
dpd_buf4_close(&LIjAb);
182
dpd_file2_close(&RIA);
183
dpd_file2_close(&Ria);
184
dpd_buf4_close(&RIJAB);
185
dpd_buf4_close(&Rijab);
186
dpd_buf4_close(&RIjAb);
189
fprintf(outfile,"\nOverlap <R|L> zero by symmetry\n");
191
} /* end dot_with_L loop */
196
/* normalizes R and produces copies of R that are ROHF-like */
197
/* sort_amps then produces the others sorted versions of R */
199
void rzero_rhf(int C_irr, int *converged) {
200
double r1, r2, rzero=0.0, energy, norm, dotval;
201
double dot_IA, dot_ia, dot_IJAB, dot_ijab, dot_IjAb;
202
dpdfile2 RIA, FIA, LIA, Lia, Ria;
203
dpdbuf4 RIjAb, RIjbA, RIjAb1, RIjbA1, D, R2, LIjAb, RIJAB, Rijab;
204
dpdbuf4 LIJAB, Lijab;
206
char lbl[32], E_lbl[32], R1A_lbl[32], R1B_lbl[32], *blank;
207
char R2AA_lbl[32], R2BB_lbl[32], R2AB_lbl[32];
210
L_irr = eom_params.L_irr;
212
for(i=0; i < eom_params.cs_per_irrep[C_irr]; i++) {
213
if (!converged[i]) continue; /* this root did not converge */
216
sprintf(E_lbl, "EOM CCSD Energy for root %d %d", C_irr, R_index);
217
if ( psio_tocscan(CC_INFO, E_lbl) == NULL) {
218
fprintf(outfile,"No EOM CCSD Energy found in CC_INFO. Not normalizing R.\n");
221
psio_read_entry(CC_INFO, E_lbl, (char *) &(energy), sizeof(double));
223
sprintf(R1A_lbl, "RIA %d %d", C_irr, R_index);
224
sprintf(R1B_lbl, "Ria %d %d", C_irr, R_index);
225
sprintf(R2AB_lbl, "RIjAb %d %d", C_irr, R_index);
226
sprintf(R2AA_lbl, "RIJAB %d %d", C_irr, R_index);
227
sprintf(R2BB_lbl, "Rijab %d %d", C_irr, R_index);
229
/* produce RIjbA and 2RIjAb-RIjbA copies - not yet normalized */
230
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
231
dpd_buf4_sort(&RIjAb, CC_TMP, pqsr, 0, 5, "RIjbA");
232
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
233
dpd_buf4_copy(&RIjAb, CC_RAMPS, lbl);
234
dpd_buf4_close(&RIjAb);
236
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
237
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
238
dpd_buf4_scm(&RIjAb, 2.0);
239
dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
240
dpd_buf4_axpy(&RIjbA, &RIjAb, -1.0);
241
dpd_buf4_close(&RIjbA);
242
dpd_buf4_close(&RIjAb);
244
/* calculate <0| hbar | 0> */
245
if (C_irr == H_IRR) {
246
dpd_file2_init(&FIA, CC_OEI, H_IRR, 0, 1, "FME");
247
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
248
r1 = 2.0 * dpd_file2_dot(&FIA, &RIA);
249
dpd_file2_close(&RIA);
250
dpd_file2_close(&FIA);
252
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
253
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
254
dpd_buf4_init(&D, CC_DINTS, H_IRR, 0, 5, 0, 5, 0, "D <ij|ab>");
255
r2 = dpd_buf4_dot(&D, &RIjAb);
257
dpd_buf4_close(&RIjAb);
258
rzero = (r1 + r2)/energy;
264
/* Now normalize R so that <R|R> = 1 */
265
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
266
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
267
dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
269
norm = norm_C_rhf(&RIA, &RIjAb, &RIjbA);
271
norm += rzero * rzero;
273
rzero = rzero / norm;
274
dpd_file2_scm(&RIA, 1.0/norm);
275
dpd_buf4_scm(&RIjAb, 1.0/norm);
276
dpd_buf4_scm(&RIjbA, 1.0/norm);
278
dpd_file2_close(&RIA);
279
dpd_buf4_close(&RIjAb);
280
dpd_buf4_close(&RIjbA);
282
fprintf(outfile,"EOM CCSD R0 for root %d = %15.10lf\n", R_index, rzero);
283
sprintf(lbl, "EOM CCSD R0 for root %d %d", C_irr, R_index);
284
psio_write_entry(CC_INFO, lbl, (char *) &rzero, sizeof(double));
286
/* produce ROHF like quantities and 2RIjAb-RIjbA */
287
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
288
dpd_file2_copy(&RIA, CC_RAMPS, R1B_lbl);
289
dpd_file2_close(&RIA);
291
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
292
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
293
dpd_buf4_copy(&RIjAb, CC_RAMPS, lbl);
294
dpd_buf4_close(&RIjAb);
296
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 2, 7, 0, 5, 1, R2AB_lbl);
297
dpd_buf4_copy(&RIjAb, CC_RAMPS, R2AA_lbl);
298
dpd_buf4_copy(&RIjAb, CC_RAMPS, R2BB_lbl);
299
dpd_buf4_close(&RIjAb);
301
dpd_buf4_init(&RIjbA, CC_TMP, C_irr, 0, 5, 0, 5, 0, "RIjbA");
302
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
303
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
304
dpd_buf4_scm(&RIjAb, 2.0);
305
dpd_buf4_axpy(&RIjbA, &RIjAb, -1.0);
306
dpd_buf4_close(&RIjAb);
307
dpd_buf4_close(&RIjbA);
309
/* test normalization with produced ROHF-like quantities */
311
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
312
norm = dpd_file2_dot_self(&RIA);
313
dpd_file2_close(&RIA);
314
dpd_file2_init(&Ria, CC_RAMPS, C_irr, 0, 1, R1B_lbl);
315
norm += dpd_file2_dot_self(&Ria);
316
dpd_file2_close(&Ria);
317
dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
318
norm += dpd_buf4_dot_self(&RIJAB);
319
dpd_buf4_close(&RIJAB);
320
dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
321
norm += dpd_buf4_dot_self(&Rijab);
322
dpd_buf4_close(&Rijab);
323
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
324
norm += dpd_buf4_dot_self(&RIjAb);
325
dpd_buf4_close(&RIjAb);
326
norm += rzero * rzero;
329
fprintf(outfile,"Norm with produced ROHF-like quantities = %15.10lf\n", norm);
330
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
331
dpd_file2_print(&RIA, outfile);
332
dpd_file2_close(&RIA);
335
/* check orthogonality with L */
336
if (eom_params.dot_with_L) {
337
if (C_irr == L_irr) {
338
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
339
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
340
r1 = 2.0 * dpd_file2_dot(&LIA, &RIA);
341
dpd_file2_close(&RIA);
342
dpd_file2_close(&LIA);
344
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "LIjAb");
345
sprintf(lbl, "%s %d %d", "2RIjAb - RIjbA", C_irr, R_index);
346
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, lbl);
347
r2 = dpd_buf4_dot(&LIjAb, &RIjAb);
348
dpd_buf4_close(&RIjAb);
349
dpd_buf4_close(&LIjAb);
351
dotval = r1 + r2 + (eom_params.L0 * rzero);
352
fprintf(outfile,"Performing RHF orthogonality test\n");
353
fprintf(outfile,"<L0|R0> = %15.10lf\n", eom_params.L0 * rzero);
354
fprintf(outfile,"2*<LIA|RIA> = %15.10lf\n", r1);
355
fprintf(outfile,"<LIjAb|2RIjAb-RIjbA> = %15.10lf\n", r2);
356
fprintf(outfile,"<L|R> = %15.10lf\n", dotval);
359
fprintf(outfile,"<L|R> zero by symmetry\n");
361
if (C_irr == L_irr ) {
362
/* double check orthogonality rohf-like */
363
dpd_file2_init(&RIA, CC_RAMPS, C_irr, 0, 1, R1A_lbl);
364
dpd_file2_init(&Ria, CC_RAMPS, C_irr, 0, 1, R1B_lbl);
365
dpd_buf4_init(&RIJAB, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2AA_lbl);
366
dpd_buf4_init(&Rijab, CC_RAMPS, C_irr, 2, 7, 2, 7, 0, R2BB_lbl);
367
dpd_buf4_init(&RIjAb, CC_RAMPS, C_irr, 0, 5, 0, 5, 0, R2AB_lbl);
369
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
370
dpd_file2_init(&Lia, CC_OEI, L_irr, 0, 1, "Lia");
371
dpd_buf4_init(&LIJAB, CC_LAMPS, L_irr, 2, 7, 2, 7, 0, "LIJAB");
372
dpd_buf4_init(&Lijab, CC_LAMPS, L_irr, 2, 7, 2, 7, 0, "Lijab");
373
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "LIjAb");
375
dot_IA = dpd_file2_dot(&LIA, &RIA);
376
dot_ia = dpd_file2_dot(&Lia, &Ria);
377
dot_IJAB = dpd_buf4_dot(&LIJAB, &RIJAB);
378
dot_ijab = dpd_buf4_dot(&Lijab, &Rijab);
379
dot_IjAb = dpd_buf4_dot(&LIjAb, &RIjAb);
381
dpd_file2_close(&RIA);
382
dpd_file2_close(&Ria);
383
dpd_buf4_close(&RIJAB);
384
dpd_buf4_close(&Rijab);
385
dpd_buf4_close(&RIjAb);
387
dpd_file2_close(&LIA);
388
dpd_file2_close(&Lia);
389
dpd_buf4_close(&LIJAB);
390
dpd_buf4_close(&Lijab);
391
dpd_buf4_close(&LIjAb);
393
fprintf(outfile,"\nROHF-like orthogonality test\n");
394
fprintf(outfile,"<L0|R0> = %15.10lf\n", eom_params.L0 * rzero);
395
fprintf(outfile,"<LIA|RIA> = %15.10lf\n", dot_IA);
396
fprintf(outfile,"<Lia|Ria> = %15.10lf\n", dot_ia);
397
fprintf(outfile,"<LIJAB|RIJAB> = %15.10lf\n", dot_IJAB);
398
fprintf(outfile,"<Lijab|Rijab> = %15.10lf\n", dot_ijab);
399
fprintf(outfile,"<LIjAb|RIjAb> = %15.10lf\n", dot_IjAb);
400
fprintf(outfile,"<L|R> = %15.10lf\n", eom_params.L0 * rzero +
401
dot_IA + dot_ia + dot_IJAB + dot_ijab + dot_IjAb);
403
} /* end dot with L, <L|R> overlap checks */
404
} /* end loop over Cs */
added
removed
src/bin/cceom/rzero.c
