← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccresponse/ccresponse.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/ccresponse/ccresponse.c
12
12
#include <libdpd/dpd.h>
13
13
#include <libchkpt/chkpt.h>
14
14
#include <physconst.h>
15
 
#include <masses.h>
16
15
#include <psifiles.h>
17
16
#include "globals.h"
18
17
 
31
30
int **cacheprep_uhf(int level, int *cachefiles);
32
31
void cachedone_uhf(int **cachelist);
33
32
void cachedone_rhf(int **cachelist);
34
 
void transmu(void);
35
 
void sortmu(void);
36
 
void mubar(void);
37
 
void transL(double sign);
38
 
void sortL(void);
39
 
void Lbar(void);
40
33
void hbar_extra(void);
 
34
void cc2_hbar_extra(void);
41
35
void sort_lamps(void);
42
 
void compute_X(char *pert, char *cart, int irrep, double omega);
43
 
double LCX(char *pert_c, char *cart_c, int irrep_c, 
44
 
           char *pert_x, char *cart_x, int irrep_x, double omega);
45
 
double HXY(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
46
 
           char *pert_y, char *cart_y, int irrep_y, double omega_y);
47
 
double LHX1Y1(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
48
 
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
49
 
double LHX2Y2(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
50
 
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
51
 
double LHX1Y2(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
52
 
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
36
void lambda_residuals(void);
53
37
 
54
38
void local_init(void);
55
39
void local_done(void);
56
40
 
 
41
void polar(void);
 
42
void optrot(void);
 
43
 
57
44
int main(int argc, char *argv[])
58
45
{
59
46
  int **cachelist, *cachefiles;
60
 
  double polar, polar_LCX, polar_HXY, polar_LHX1Y1, polar_LHX2Y2, polar_LHX1Y2;
61
 
  char **cartcomp;
62
 
  int alpha, beta, i;
63
 
  double **tensor;
64
 
  double omega_nm, omega_ev, omega_cm;
65
 
  double prefactor, TrG, M, nu, rotation, bohr2a4, m2a, hbar;
66
47
 
67
48
  init_io(argc, argv);
68
49
  init_ioff();
88
69
 
89
70
  if(params.local) local_init();
90
71
 
91
 
  hbar_extra();
92
 
  /* sort_lamps(); now provided by cclambda */
93
 
 
94
 
  cartcomp = (char **) malloc(3 * sizeof(char *));
95
 
  cartcomp[0] = strdup("X");
96
 
  cartcomp[1] = strdup("Y");
97
 
  cartcomp[2] = strdup("Z");
98
 
 
99
 
  if(!strcmp(params.prop,"POLARIZABILITY") || !strcmp(params.prop,"ALL") || 
100
 
     !strcmp(params.prop,"ROTATION")) {
101
 
 
102
 
    tensor = block_matrix(3,3);
103
 
 
104
 
    /* prepare electric dipole integrals */
105
 
    transmu();
106
 
    sortmu();
107
 
    mubar();
108
 
 
109
 
    /* Compute the electric-dipole-perturbed CC wave functions */
110
 
    for(alpha=0; alpha < 3; alpha++) {
111
 
      compute_X("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega);
112
 
      if(params.omega != 0.0) compute_X("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega);
113
 
    }
114
 
 
115
 
    for(alpha=0; alpha < 3; alpha++) {
116
 
      for(beta=0; beta < 3; beta++) {
117
 
 
118
 
        tensor[alpha][beta] = 0.0;
119
 
        polar_LCX = 0.0;
120
 
        polar_HXY = 0.0;
121
 
        polar_LHX1Y1 = 0.0;
122
 
        polar_LHX2Y2 = 0.0;
123
 
        polar_LHX1Y2 = 0.0;
124
 
 
125
 
        if(!(moinfo.mu_irreps[alpha]^moinfo.mu_irreps[beta])) {
126
 
 
127
 
          if(params.omega != 0.0) {
128
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "Mu", cartcomp[beta], 
129
 
                            moinfo.mu_irreps[beta], params.omega);
130
 
            polar_LCX += LCX("Mu", cartcomp[beta], moinfo.mu_irreps[beta], "Mu", cartcomp[alpha],
131
 
                             moinfo.mu_irreps[alpha], -params.omega);
132
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
133
 
                            "Mu", cartcomp[beta], moinfo.mu_irreps[beta], params.omega);
134
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
135
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], params.omega);
136
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
137
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], params.omega);
138
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
139
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], params.omega);
140
 
            polar_LHX1Y2 += LHX1Y2("Mu", cartcomp[beta], moinfo.mu_irreps[beta], params.omega,
141
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega);
142
 
          }
143
 
          else {
144
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "Mu", cartcomp[beta],
145
 
                            moinfo.mu_irreps[beta], 0.0);
146
 
            polar_LCX += LCX("Mu", cartcomp[beta], moinfo.mu_irreps[beta], "Mu", cartcomp[alpha],
147
 
                             moinfo.mu_irreps[alpha], 0.0);
148
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
149
 
                            "Mu", cartcomp[beta], moinfo.mu_irreps[beta], 0.0);
150
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
151
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], 0.0);
152
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
153
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], 0.0);
154
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
155
 
                                  "Mu", cartcomp[beta], moinfo.mu_irreps[beta], 0.0);
156
 
            polar_LHX1Y2 += LHX1Y2("Mu", cartcomp[beta], moinfo.mu_irreps[beta], 0.0,
157
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0);
158
 
          }
159
 
 
160
 
          polar = polar_LCX + polar_HXY + polar_LHX1Y1 + polar_LHX2Y2 + polar_LHX1Y2;
161
 
          /*
162
 
          fprintf(outfile, "polar_LCX    = %20.12f\n", polar_LCX);
163
 
          fprintf(outfile, "polar_HXY    = %20.12f\n", polar_HXY);
164
 
          fprintf(outfile, "polar_LHX1Y1 = %20.12f\n", polar_LHX1Y1);
165
 
          fprintf(outfile, "polar_LHX2Y2 = %20.12f\n", polar_LHX2Y2);
166
 
          fprintf(outfile, "polar_LHX1Y2 = %20.12f\n", polar_LHX1Y2);
167
 
          */
168
 
 
169
 
          tensor[alpha][beta] = -polar;
170
 
        }
171
 
      }
172
 
    }
173
 
 
174
 
    fprintf(outfile, "\n                 CCSD Dipole Polarizability [(e^2 a0^2)/E_h]:\n");
175
 
    fprintf(outfile, "  -------------------------------------------------------------------------\n");
176
 
    if(params.omega != 0.0) 
177
 
      omega_nm = (_c*_h*1e9)/(_hartree2J*params.omega);
178
 
    omega_ev = _hartree2ev*params.omega;
179
 
    omega_cm = _hartree2wavenumbers*params.omega;
180
 
    if(params.omega != 0.0)
181
 
      fprintf(outfile,   "   Evaluated at omega = %8.6f E_h (%6.2f nm, %5.3f eV, %8.2f cm-1)\n", 
182
 
          params.omega, omega_nm, omega_ev, omega_cm);
183
 
    else
184
 
      fprintf(outfile,   "   Evaluated at omega = %8.6f E_h (Inf nm, %5.3f eV, %8.2f cm-1)\n", 
185
 
          params.omega, omega_ev, omega_cm);
186
 
    fprintf(outfile, "  -------------------------------------------------------------------------\n");
187
 
    mat_print(tensor, 3, 3, outfile);
188
 
 
189
 
    free_block(tensor);
190
 
  }
191
 
 
192
 
  /*** Optical rotation ***/
193
 
 
194
 
  /* prepare magnetic dipole integrals */
195
 
  if(!strcmp(params.prop,"ROTATION") || !strcmp(params.prop,"ALL")) {
196
 
 
197
 
    tensor = block_matrix(3,3);
198
 
 
199
 
    /* prepare the magnetic-dipole integrals */
200
 
    transL(+1.0);
201
 
    sortL();
202
 
    Lbar();
203
 
 
204
 
    /* Compute the +omega magnetic-dipole-perturbed CC wave functions */
205
 
    /* NB: The -omega electric-dipole perturbed wfns should already be available */
206
 
    for(alpha=0; alpha < 3; alpha++)
207
 
      compute_X("L", cartcomp[alpha], moinfo.l_irreps[alpha], params.omega);
208
 
 
209
 
    for(alpha=0; alpha < 3; alpha++) {
210
 
      for(beta=0; beta < 3; beta++) {
211
 
 
212
 
        polar_LCX = 0.0;
213
 
        polar_HXY = 0.0;
214
 
        polar_LHX1Y1 = 0.0;
215
 
        polar_LHX2Y2 = 0.0;
216
 
        polar_LHX1Y2 = 0.0;
217
 
 
218
 
        if(!(moinfo.mu_irreps[alpha]^moinfo.l_irreps[beta])) {
219
 
 
220
 
          if(params.omega != 0.0) {
221
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "L", cartcomp[beta], 
222
 
                            moinfo.l_irreps[beta], params.omega);
223
 
            polar_LCX += LCX("L", cartcomp[beta], moinfo.l_irreps[beta], "Mu", cartcomp[alpha],
224
 
                             moinfo.mu_irreps[alpha], -params.omega);
225
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
226
 
                            "L", cartcomp[beta], moinfo.l_irreps[beta], params.omega);
227
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
228
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], params.omega);
229
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
230
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], params.omega);
231
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega,
232
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], params.omega);
233
 
            polar_LHX1Y2 += LHX1Y2("L", cartcomp[beta], moinfo.l_irreps[beta], params.omega,
234
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], -params.omega);
235
 
          }
236
 
          else {
237
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "L", cartcomp[beta],
238
 
                            moinfo.l_irreps[beta], 0.0);
239
 
            polar_LCX += LCX("L", cartcomp[beta], moinfo.l_irreps[beta], "Mu", cartcomp[alpha],
240
 
                             moinfo.mu_irreps[alpha], 0.0);
241
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
242
 
                            "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
243
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
244
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
245
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
246
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
247
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
248
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
249
 
            polar_LHX1Y2 += LHX1Y2("L", cartcomp[beta], moinfo.l_irreps[beta], 0.0,
250
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0);
251
 
          }
252
 
 
253
 
          polar = polar_LCX + polar_HXY + polar_LHX1Y1 + polar_LHX2Y2 + polar_LHX1Y2;
254
 
          /*
255
 
          fprintf(outfile, "polar_LCX    = %20.12f\n", polar_LCX);
256
 
          fprintf(outfile, "polar_HXY    = %20.12f\n", polar_HXY);
257
 
          fprintf(outfile, "polar_LHX1Y1 = %20.12f\n", polar_LHX1Y1);
258
 
          fprintf(outfile, "polar_LHX2Y2 = %20.12f\n", polar_LHX2Y2);
259
 
          fprintf(outfile, "polar_LHX1Y2 = %20.12f\n", polar_LHX1Y2);
260
 
          */
261
 
 
262
 
          tensor[alpha][beta] = 0.5 * polar;
263
 
        }
264
 
        /*      fprintf(outfile, "%1s%1s polar = %20.12f\n", cartcomp[alpha], cartcomp[beta], polar); */
265
 
      }
266
 
    }
267
 
 
268
 
    /* prepare the complex-conjugate of the magnetic-dipole integrals */
269
 
    transL(-1.0);
270
 
    sortL();
271
 
    Lbar();
272
 
 
273
 
    /* Compute the -omega cc-magnetic-dipole-perturbed CC wave functions */
274
 
    /* NB: The +omega electric-dipole perturbed wfns should already be available */
275
 
    for(alpha=0; alpha < 3; alpha++)
276
 
      compute_X("L", cartcomp[alpha], moinfo.l_irreps[alpha], -params.omega);
277
 
 
278
 
    for(alpha=0; alpha < 3; alpha++) {
279
 
      for(beta=0; beta < 3; beta++) {
280
 
 
281
 
        polar_LCX = 0.0;
282
 
        polar_HXY = 0.0;
283
 
        polar_LHX1Y1 = 0.0;
284
 
        polar_LHX2Y2 = 0.0;
285
 
        polar_LHX1Y2 = 0.0;
286
 
 
287
 
        if(!(moinfo.mu_irreps[alpha]^moinfo.l_irreps[beta])) {
288
 
 
289
 
          if(params.omega != 0.0) {
290
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "L", cartcomp[beta], 
291
 
                            moinfo.l_irreps[beta], -params.omega);
292
 
            polar_LCX += LCX("L", cartcomp[beta], moinfo.l_irreps[beta], "Mu", cartcomp[alpha],
293
 
                             moinfo.mu_irreps[alpha], params.omega);
294
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega,
295
 
                            "L", cartcomp[beta], moinfo.l_irreps[beta], -params.omega);
296
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega,
297
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], -params.omega);
298
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega,
299
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], -params.omega);
300
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega,
301
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], -params.omega);
302
 
            polar_LHX1Y2 += LHX1Y2("L", cartcomp[beta], moinfo.l_irreps[beta], -params.omega,
303
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], params.omega);
304
 
          }
305
 
          else {
306
 
            polar_LCX = LCX("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], "L", cartcomp[beta],
307
 
                            moinfo.l_irreps[beta], 0.0);
308
 
            polar_LCX += LCX("L", cartcomp[beta], moinfo.l_irreps[beta], "Mu", cartcomp[alpha],
309
 
                             moinfo.mu_irreps[alpha], 0.0);
310
 
            polar_HXY = HXY("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
311
 
                            "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
312
 
            polar_LHX1Y1 = LHX1Y1("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
313
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
314
 
            polar_LHX2Y2 = LHX2Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
315
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
316
 
            polar_LHX1Y2 = LHX1Y2("Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0,
317
 
                                  "L", cartcomp[beta], moinfo.l_irreps[beta], 0.0);
318
 
            polar_LHX1Y2 += LHX1Y2("L", cartcomp[beta], moinfo.l_irreps[beta], 0.0,
319
 
                                   "Mu", cartcomp[alpha], moinfo.mu_irreps[alpha], 0.0);
320
 
          }
321
 
 
322
 
          polar = polar_LCX + polar_HXY + polar_LHX1Y1 + polar_LHX2Y2 + polar_LHX1Y2;
323
 
          /*
324
 
          fprintf(outfile, "polar_LCX    = %20.12f\n", polar_LCX);
325
 
          fprintf(outfile, "polar_HXY    = %20.12f\n", polar_HXY);
326
 
          fprintf(outfile, "polar_LHX1Y1 = %20.12f\n", polar_LHX1Y1);
327
 
          fprintf(outfile, "polar_LHX2Y2 = %20.12f\n", polar_LHX2Y2);
328
 
          fprintf(outfile, "polar_LHX1Y2 = %20.12f\n", polar_LHX1Y2);
329
 
          */
330
 
          tensor[alpha][beta] += 0.5 * polar;
331
 
        }
332
 
        /*      fprintf(outfile, "%1s%1s polar = %20.12f\n", cartcomp[alpha], cartcomp[beta], polar); */
333
 
      }
334
 
    }
335
 
 
336
 
    fprintf(outfile, "\n                      CCSD Optical Rotation Tensor:\n");
337
 
    fprintf(outfile, "  -------------------------------------------------------------------------\n");
338
 
    fprintf(outfile,   "   Evaluated at omega = %8.6f E_h (%6.2f nm, %5.3f eV, %8.2f cm-1)\n", params.omega,
339
 
            (_c*_h*1e9)/(_hartree2J*params.omega), _hartree2ev*params.omega,
340
 
            _hartree2wavenumbers*params.omega);
341
 
    fprintf(outfile, "  -------------------------------------------------------------------------\n");
342
 
    mat_print(tensor, 3, 3, outfile);
343
 
 
344
 
    /* compute the specific rotation */
345
 
    for(i=0,M=0.0; i < moinfo.natom ;i++) M += an2masses[(int) moinfo.zvals[i]]; /* amu */
346
 
    TrG = (tensor[0][0] + tensor[1][1] + tensor[2][2])/(3.0 * params.omega);
347
 
    nu = params.omega; /* hartree */
348
 
    bohr2a4 = _bohr2angstroms * _bohr2angstroms * _bohr2angstroms * _bohr2angstroms;
349
 
    m2a = _bohr2angstroms * 1.0e-10;
350
 
    hbar = _h/(2.0 * _pi);
351
 
    prefactor = 1.0e-2 * hbar/(_c * 2.0 * _pi * _me * m2a * m2a);
352
 
    prefactor *= prefactor;
353
 
    prefactor *= 288.0e-30 * _pi * _pi * _na * bohr2a4;
354
 
    rotation = prefactor * TrG * nu * nu / M;
355
 
    fprintf(outfile, "\n\t[alpha]_(%5.3f) = %20.12f deg/[dm (gm/cm^3)]\n", params.omega, rotation);
356
 
 
357
 
    free_block(tensor);
358
 
  }
359
 
 
360
 
  for(alpha=0; alpha < 3; alpha++) free(cartcomp[alpha]);
361
 
  free(cartcomp);
362
 
 
 
72
  if (!strcmp(params.wfn,"CC2")) {
 
73
    cc2_hbar_extra();
 
74
  }
 
75
  else {
 
76
    hbar_extra();
 
77
  }
 
78
 
 
79
  sort_lamps(); /* should be removed sometime - provided by cclambda */
 
80
  if(strcmp(params.wfn,"CC2")) lambda_residuals(); /* don't do this for CC2 */
 
81
 
 
82
  if(!strcmp(params.prop,"POLARIZABILITY")) polar();
 
83
  if(!strcmp(params.prop,"ROTATION")) optrot();
363
84
 
364
85
  if(params.local) local_done();
365
86
 
392
113
  psio_init();
393
114
 
394
115
  for(i=CC_MIN; i <= CC_MAX; i++) psio_open(i, 1);
 
116
 
 
117
  /* Clear out DIIS TOC Entries */
 
118
  psio_close(CC_DIIS_AMP, 0);
 
119
  psio_close(CC_DIIS_ERR, 0);
 
120
 
 
121
  psio_open(CC_DIIS_AMP, 0);
 
122
  psio_open(CC_DIIS_ERR, 0);
395
123
}
396
124
 
397
125
void title(void)
407
135
{
408
136
  int i;
409
137
 
410
 
  for(i=CC_MIN; i <= CC_MAX; i++) psio_close(i, 1);
 
138
  /* Close all dpd data files here */
 
139
  for(i=CC_MIN; i < CC_TMP; i++) psio_close(i,1);
 
140
  for(i=CC_TMP; i <= CC_TMP11; i++) psio_close(i,0);  /* get rid of TMP files */
 
141
  for(i=CC_TMP11+1; i <= CC_MAX; i++) psio_close(i,1);
411
142
 
412
143
  psio_done();
413
144
  tstop(outfile);