← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/diis.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/cclambda/diis.cc
 
1
/*! \file
 
2
    \ingroup CCLAMBDA
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cmath>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libpsio/psio.h>
 
10
#include <libdpd/dpd.h>
 
11
#include <libqt/qt.h>
 
12
#include <psifiles.h>
 
13
#include "MOInfo.h"
 
14
#include "Params.h"
 
15
#define EXTERN
 
16
#include "globals.h"
 
17
 
 
18
namespace psi { namespace cclambda {
 
19
 
 
20
/*
 
21
** DIIS: Direct inversion in the iterative subspace routine to
 
22
** accelerate convergence of the CCSD amplitude equations.
 
23
**
 
24
** SubstanLially improved efficiency of this routine:
 
25
** (1) Keeping at most two error vectors in core at once.
 
26
** (2) Limiting direct product (overlap) calculation to unique pairs.
 
27
** (3) Using LAPACK's linear equation solver DGESV instead of flin.
 
28
**
 
29
** These improvements have been applied only to RHF cases so far.
 
30
**
 
31
** -TDC  12/22/01
 
32
*/
 
33
 
 
34
void diis(int iter, int L_irr)
 
35
{
 
36
  int nvector=8;  /* Number of error vectors to keep */
 
37
  int h, nirreps;
 
38
  int row, col, word, p, q, i;
 
39
  int diis_cycle;
 
40
  int vector_length=0;
 
41
  int errcod, *ipiv;
 
42
  dpdfile2 L1, L1a, L1b;
 
43
  dpdbuf4 L2, L2a, L2b, L2c;
 
44
  psio_address start, end, next;
 
45
  double **error;
 
46
  double **B, *C, **vector;
 
47
  double product, determinant, maximum;
 
48
 
 
49
  nirreps = moinfo.nirreps;
 
50
 
 
51
  if(params.ref == 0) { /** RHF **/
 
52
    /* Compute the length of a single error vector */
 
53
    dpd_file2_init(&L1, CC_LAMBDA, L_irr, 0, 1, "LIA");
 
54
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
55
    for(h=0; h < nirreps; h++) {
 
56
      vector_length += L1.params->rowtot[h] * L1.params->coltot[h^L_irr];
 
57
      vector_length += L2.params->rowtot[h] * L2.params->coltot[h^L_irr];
 
58
    }
 
59
    dpd_file2_close(&L1);
 
60
    dpd_buf4_close(&L2);
 
61
 
 
62
    /* Set the diis cycle value */
 
63
    diis_cycle = (iter-1) % nvector;
 
64
 
 
65
    /* Build the current error vector and dump it to disk */
 
66
    error = dpd_block_matrix(1,vector_length);
 
67
 
 
68
    word=0;
 
69
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
70
    /*dpd_file2_print(&L1a,outfile);*/
 
71
    dpd_file2_mat_init(&L1a);
 
72
    dpd_file2_mat_rd(&L1a);
 
73
    dpd_file2_init(&L1b, CC_LAMBDA, L_irr, 0, 1, "LIA");
 
74
    /*dpd_file2_print(&L1b,outfile);*/
 
75
    dpd_file2_mat_init(&L1b);
 
76
    dpd_file2_mat_rd(&L1b);
 
77
    for(h=0; h < nirreps; h++)
 
78
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
79
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++) {
 
80
          error[0][word++] = L1a.matrix[h][row][col] - L1b.matrix[h][row][col];
 
81
    }
 
82
 
 
83
    dpd_file2_mat_close(&L1a);
 
84
    dpd_file2_close(&L1a);
 
85
    dpd_file2_mat_close(&L1b);
 
86
    dpd_file2_close(&L1b);
 
87
 
 
88
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
89
    /*dpd_buf4_print(&L2a,outfile,1);*/
 
90
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
91
    /*dpd_buf4_print(&L2b,outfile,1);*/
 
92
    for(h=0; h < nirreps; h++) {
 
93
      dpd_buf4_mat_irrep_init(&L2a, h);
 
94
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
95
      dpd_buf4_mat_irrep_init(&L2b, h);
 
96
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
97
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
98
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++) {
 
99
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
100
/* fprintf(outfile,"%15.10lf\n", error[0][word-1]); */
 
101
    }
 
102
      dpd_buf4_mat_irrep_close(&L2a, h);
 
103
      dpd_buf4_mat_irrep_close(&L2b, h);
 
104
    }
 
105
    dpd_buf4_close(&L2a);
 
106
    dpd_buf4_close(&L2b);
 
107
 
 
108
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
109
    psio_write(CC_DIIS_ERR, "DIIS Error Vectors" , (char *) error[0], 
 
110
               vector_length*sizeof(double), start, &end);
 
111
 
 
112
 
 
113
    /* Store the current amplitude vector on disk */
 
114
    word=0;
 
115
 
 
116
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
117
    dpd_file2_mat_init(&L1a);
 
118
    dpd_file2_mat_rd(&L1a);
 
119
    for(h=0; h < nirreps; h++)
 
120
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
121
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
122
          error[0][word++] = L1a.matrix[h][row][col];
 
123
    dpd_file2_mat_close(&L1a);
 
124
    dpd_file2_close(&L1a);
 
125
 
 
126
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
127
    for(h=0; h < nirreps; h++) {
 
128
      dpd_buf4_mat_irrep_init(&L2a, h);
 
129
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
130
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
131
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
132
          error[0][word++] = L2a.matrix[h][row][col];
 
133
      dpd_buf4_mat_irrep_close(&L2a, h);
 
134
    }
 
135
    dpd_buf4_close(&L2a);
 
136
 
 
137
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
138
    psio_write(CC_DIIS_AMP, "DIIS Amplitude Vectors" , (char *) error[0], 
 
139
               vector_length*sizeof(double), start, &end);
 
140
 
 
141
    /* If we haven't run through enough iterations, set the correct dimensions
 
142
       for the extrapolation */
 
143
    if(!(iter >= (nvector))) {
 
144
      if(iter < 2) { /* Leave if we can't extrapolate at all */
 
145
        dpd_free_block(error, 1, vector_length);
 
146
        return; 
 
147
      }
 
148
      nvector = iter;
 
149
    }
 
150
 
 
151
    /* Build B matrix of error vector products */
 
152
    vector = dpd_block_matrix(2, vector_length);
 
153
    B = block_matrix(nvector+1,nvector+1);
 
154
    for(p=0; p < nvector; p++) {
 
155
 
 
156
      start = psio_get_address(PSIO_ZERO, p*vector_length*sizeof(double));
 
157
 
 
158
      psio_read(CC_DIIS_ERR, "DIIS Error Vectors", (char *) vector[0],
 
159
                vector_length*sizeof(double), start, &end);
 
160
 
 
161
      /*
 
162
      for(i=0; i < vector_length; i++) 
 
163
        fprintf(outfile,"E[%d][%d] = %20.15lf\n",p,i,vector[0][i]);
 
164
      */
 
165
 
 
166
      dot_arr(vector[0], vector[0], vector_length, &product);
 
167
 
 
168
      B[p][p] = product;
 
169
 
 
170
      for(q=0; q < p; q++) {
 
171
 
 
172
        start = psio_get_address(PSIO_ZERO, q*vector_length*sizeof(double));
 
173
 
 
174
        psio_read(CC_DIIS_ERR, "DIIS Error Vectors", (char *) vector[1],
 
175
                  vector_length*sizeof(double), start, &end);
 
176
 
 
177
        dot_arr(vector[1], vector[0], vector_length, &product);
 
178
 
 
179
        B[p][q] = B[q][p] = product;
 
180
      }
 
181
    }
 
182
    dpd_free_block(vector, 2, vector_length);
 
183
 
 
184
    for(p=0; p < nvector; p++) {
 
185
      B[p][nvector] = -1;
 
186
      B[nvector][p] = -1;
 
187
    }
 
188
 
 
189
    B[nvector][nvector] = 0;
 
190
 
 
191
    /* Find the maximum value in B and scale all its elements */
 
192
    maximum = fabs(B[0][0]);
 
193
    for(p=0; p < nvector; p++)
 
194
      for(q=0; q < nvector; q++)
 
195
        if(fabs(B[p][q]) > maximum) maximum = fabs(B[p][q]);
 
196
 
 
197
    for(p=0; p < nvector; p++)
 
198
      for(q=0; q < nvector; q++)
 
199
        B[p][q] /= maximum; 
 
200
 
 
201
    /*
 
202
    fprintf(outfile,"\nDIIS B:\n");
 
203
    print_mat(B,nvector,nvector,outfile);
 
204
    */
 
205
 
 
206
    /* Build the constant vector */
 
207
    C = init_array(nvector+1);
 
208
    C[nvector] = -1;
 
209
 
 
210
    /* Solve the linear equations */
 
211
    ipiv = init_int_array(nvector+1);
 
212
 
 
213
    errcod = C_DGESV(nvector+1, 1, &(B[0][0]), nvector+1, &(ipiv[0]), &(C[0]), nvector+1);
 
214
    if(errcod) {
 
215
      fprintf(outfile, "\nError in DGESV return in diis.\n");
 
216
      exit(PSI_RETURN_FAILURE);
 
217
    }
 
218
 
 
219
    /* Build a new amplitude vector from the old ones */
 
220
    vector = dpd_block_matrix(1, vector_length);
 
221
    for(p=0; p < vector_length; p++) error[0][p] = 0.0;
 
222
    for(p=0; p < nvector; p++) {
 
223
      /*fprintf(outfile,"C[%d] = %20.15lf\n",p,C[p]);*/
 
224
 
 
225
      start = psio_get_address(PSIO_ZERO, p*vector_length*sizeof(double));
 
226
 
 
227
      psio_read(CC_DIIS_AMP, "DIIS Amplitude Vectors", (char *) vector[0], 
 
228
                vector_length*sizeof(double), start, &end);
 
229
 
 
230
      for(q=0; q < vector_length; q++) 
 
231
        error[0][q] += C[p] * vector[0][q];
 
232
 
 
233
    }
 
234
    dpd_free_block(vector, 1, vector_length);
 
235
 
 
236
    /* Now place these elements into the DPD amplitude arrays */
 
237
    word=0;
 
238
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
239
    dpd_file2_mat_init(&L1a);
 
240
    for(h=0; h < nirreps; h++)
 
241
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
242
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
243
          L1a.matrix[h][row][col] = error[0][word++];
 
244
    dpd_file2_mat_wrt(&L1a);
 
245
    dpd_file2_mat_close(&L1a);
 
246
    dpd_file2_copy(&L1a, CC_LAMBDA, "New Lia"); /* to be removed after spin-adaptation */
 
247
    /*dpd_file2_print(&L1a,outfile);*/
 
248
    dpd_file2_close(&L1a);
 
249
 
 
250
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
251
    for(h=0; h < nirreps; h++) {
 
252
      dpd_buf4_mat_irrep_init(&L2a, h);
 
253
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
254
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
255
          L2a.matrix[h][row][col] = error[0][word++];
 
256
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
257
      dpd_buf4_mat_irrep_close(&L2a, h);
 
258
    }
 
259
    dpd_buf4_close(&L2a);
 
260
 
 
261
    /* to be removed after spin-adaptation */
 
262
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 0, 5, 1, "New LIjAb");
 
263
    dpd_buf4_copy(&L2a, CC_LAMBDA, "New LIJAB");
 
264
    dpd_buf4_copy(&L2a, CC_LAMBDA, "New Lijab");
 
265
    dpd_buf4_close(&L2a);
 
266
 
 
267
    /* Release memory and return */
 
268
    /*    free_matrix(vector, nvector); */
 
269
    free_block(B);
 
270
    free(C);
 
271
    free(ipiv);
 
272
    dpd_free_block(error, 1, vector_length);
 
273
  }
 
274
  else if(params.ref == 1) { /** ROHF **/
 
275
  
 
276
    /* Compute the length of a single error vector */
 
277
    /* RAK changed file nums here from CC_TMP0 */
 
278
    dpd_file2_init(&L1, CC_LAMBDA, L_irr, 0, 1, "LIA");
 
279
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
 
280
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
281
    for(h=0; h < nirreps; h++) {
 
282
      vector_length += 2 * L1.params->rowtot[h] * L1.params->coltot[h^L_irr];
 
283
      vector_length += 2 * L2a.params->rowtot[h] * L2a.params->coltot[h^L_irr];
 
284
      vector_length += L2b.params->rowtot[h] * L2b.params->coltot[h^L_irr];
 
285
    }
 
286
    dpd_file2_close(&L1);
 
287
    dpd_buf4_close(&L2a);
 
288
    dpd_buf4_close(&L2b);
 
289
 
 
290
    /* Set the diis cycle value */
 
291
    diis_cycle = (iter-1) % nvector;
 
292
 
 
293
    /* Build the current error vector and dump it to disk */
 
294
    error = dpd_block_matrix(1,vector_length);
 
295
    word=0;
 
296
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
297
    dpd_file2_mat_init(&L1a);
 
298
    dpd_file2_mat_rd(&L1a);
 
299
    dpd_file2_init(&L1b, CC_LAMBDA, L_irr, 0, 1, "LIA");
 
300
    dpd_file2_mat_init(&L1b);
 
301
    dpd_file2_mat_rd(&L1b);
 
302
    for(h=0; h < nirreps; h++)
 
303
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
304
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
305
          error[0][word++] = L1a.matrix[h][row][col] - L1b.matrix[h][row][col];
 
306
    dpd_file2_mat_close(&L1a);
 
307
    dpd_file2_close(&L1a);
 
308
    dpd_file2_mat_close(&L1b);
 
309
    dpd_file2_close(&L1b);
 
310
 
 
311
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New Lia");
 
312
    dpd_file2_mat_init(&L1a);
 
313
    dpd_file2_mat_rd(&L1a);
 
314
    dpd_file2_init(&L1b, CC_LAMBDA, L_irr, 0, 1, "Lia");
 
315
    dpd_file2_mat_init(&L1b);
 
316
    dpd_file2_mat_rd(&L1b);
 
317
    for(h=0; h < nirreps; h++)
 
318
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
319
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
320
          error[0][word++] = L1a.matrix[h][row][col] - L1b.matrix[h][row][col];
 
321
    dpd_file2_mat_close(&L1a);
 
322
    dpd_file2_close(&L1a);
 
323
    dpd_file2_mat_close(&L1b);
 
324
    dpd_file2_close(&L1b);
 
325
  
 
326
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
327
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
 
328
    for(h=0; h < nirreps; h++) {
 
329
      dpd_buf4_mat_irrep_init(&L2a, h);
 
330
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
331
      dpd_buf4_mat_irrep_init(&L2b, h);
 
332
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
333
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
334
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
335
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
336
      dpd_buf4_mat_irrep_close(&L2a, h);
 
337
      dpd_buf4_mat_irrep_close(&L2b, h);
 
338
    }
 
339
    dpd_buf4_close(&L2a);
 
340
    dpd_buf4_close(&L2b);
 
341
 
 
342
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New Lijab");
 
343
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "Lijab");
 
344
    for(h=0; h < nirreps; h++) {
 
345
      dpd_buf4_mat_irrep_init(&L2a, h);
 
346
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
347
      dpd_buf4_mat_irrep_init(&L2b, h);
 
348
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
349
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
350
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
351
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
352
      dpd_buf4_mat_irrep_close(&L2a, h);
 
353
      dpd_buf4_mat_irrep_close(&L2b, h);
 
354
    }
 
355
    dpd_buf4_close(&L2a);
 
356
    dpd_buf4_close(&L2b);
 
357
 
 
358
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
359
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
 
360
    for(h=0; h < nirreps; h++) {
 
361
      dpd_buf4_mat_irrep_init(&L2a, h);
 
362
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
363
      dpd_buf4_mat_irrep_init(&L2b, h);
 
364
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
365
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
366
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
367
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
368
      dpd_buf4_mat_irrep_close(&L2a, h);
 
369
      dpd_buf4_mat_irrep_close(&L2b, h);
 
370
    }
 
371
    dpd_buf4_close(&L2a);
 
372
    dpd_buf4_close(&L2b);
 
373
 
 
374
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
375
    psio_write(CC_DIIS_ERR, "DIIS Error Vectors" , (char *) error[0], 
 
376
               vector_length*sizeof(double), start, &end);
 
377
 
 
378
    /* Store the current amplitude vector on disk */
 
379
    word=0;
 
380
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
381
    dpd_file2_mat_init(&L1a);
 
382
    dpd_file2_mat_rd(&L1a);
 
383
    for(h=0; h < nirreps; h++)
 
384
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
385
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
386
          error[0][word++] = L1a.matrix[h][row][col];
 
387
    dpd_file2_mat_close(&L1a);
 
388
    dpd_file2_close(&L1a);
 
389
 
 
390
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New Lia");
 
391
    dpd_file2_mat_init(&L1a);
 
392
    dpd_file2_mat_rd(&L1a);
 
393
    for(h=0; h < nirreps; h++)
 
394
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
395
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
396
          error[0][word++] = L1a.matrix[h][row][col];
 
397
    dpd_file2_mat_close(&L1a);
 
398
    dpd_file2_close(&L1a);
 
399
  
 
400
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
401
    for(h=0; h < nirreps; h++) {
 
402
      dpd_buf4_mat_irrep_init(&L2a, h);
 
403
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
404
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
405
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
406
          error[0][word++] = L2a.matrix[h][row][col];
 
407
      dpd_buf4_mat_irrep_close(&L2a, h);
 
408
    }
 
409
    dpd_buf4_close(&L2a);
 
410
 
 
411
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New Lijab");
 
412
    for(h=0; h < nirreps; h++) {
 
413
      dpd_buf4_mat_irrep_init(&L2a, h);
 
414
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
415
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
416
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
417
          error[0][word++] = L2a.matrix[h][row][col];
 
418
      dpd_buf4_mat_irrep_close(&L2a, h);
 
419
    }
 
420
    dpd_buf4_close(&L2a);
 
421
 
 
422
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
423
    for(h=0; h < nirreps; h++) {
 
424
      dpd_buf4_mat_irrep_init(&L2a, h);
 
425
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
426
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
427
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
428
          error[0][word++] = L2a.matrix[h][row][col];
 
429
      dpd_buf4_mat_irrep_close(&L2a, h);
 
430
    }
 
431
    dpd_buf4_close(&L2a);
 
432
 
 
433
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
434
    psio_write(CC_DIIS_AMP, "DIIS Amplitude Vectors" , (char *) error[0], 
 
435
               vector_length*sizeof(double), start, &end);
 
436
 
 
437
    /* If we haven't run through enough iterations, set the correct dimensions
 
438
       for the extrapolation */
 
439
    if(!(iter >= (nvector))) {
 
440
      if(iter < 2) { /* Leave if we can't extrapolate at all */
 
441
        free(error);
 
442
        return; 
 
443
      }
 
444
      nvector = iter;
 
445
    }
 
446
 
 
447
    /* Now grab the full set of error vectors from the file */
 
448
    vector = init_matrix(nvector, vector_length);
 
449
    next = PSIO_ZERO;
 
450
    for(p=0; p < nvector; p++) 
 
451
      psio_read(CC_DIIS_ERR, "DIIS Error Vectors", (char *) vector[p], 
 
452
                vector_length*sizeof(double), next, &next);
 
453
 
 
454
    /* Build B matrix of error vector products */
 
455
    B = init_matrix(nvector+1,nvector+1);
 
456
 
 
457
    for(p=0; p < nvector; p++)
 
458
      for(q=0; q < nvector; q++) {
 
459
        dot_arr(vector[p], vector[q], vector_length, &product); 
 
460
        B[p][q] = product;
 
461
      }
 
462
 
 
463
    for(p=0; p < nvector; p++) {
 
464
      B[p][nvector] = -1;
 
465
      B[nvector][p] = -1;
 
466
    }
 
467
 
 
468
    B[nvector][nvector] = 0;
 
469
 
 
470
    /* Find the maximum value in B and scale all its elements */
 
471
    maximum = fabs(B[0][0]);
 
472
    for(p=0; p < nvector; p++)
 
473
      for(q=0; q < nvector; q++)
 
474
        if(fabs(B[p][q]) > maximum) maximum = fabs(B[p][q]);
 
475
 
 
476
    for(p=0; p < nvector; p++)
 
477
      for(q=0; q < nvector; q++)
 
478
        B[p][q] /= maximum; 
 
479
 
 
480
    /* Build the constant vector */
 
481
    C = init_array(nvector+1);
 
482
    C[nvector] = -1;
 
483
 
 
484
    /* Solve the linear equations */
 
485
    flin(B, C, nvector+1, 1, &determinant);
 
486
 
 
487
    /* Grab the old amplitude vectors */
 
488
    next = PSIO_ZERO;
 
489
    for(p=0; p < nvector; p++) 
 
490
      psio_read(CC_DIIS_AMP, "DIIS Amplitude Vectors", (char *) vector[p], 
 
491
                vector_length*sizeof(double), next, &next);
 
492
  
 
493
    /* Build the new amplitude vector from the old ones */
 
494
    for(q=0; q < vector_length; q++) {
 
495
      error[0][q] = 0.0;
 
496
      for(p=0; p < nvector; p++)
 
497
        error[0][q] += C[p] * vector[p][q];
 
498
    }
 
499
 
 
500
    /* Now place these elements into the DPD amplitude arrays */
 
501
    word=0;
 
502
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
503
    dpd_file2_mat_init(&L1a);
 
504
    for(h=0; h < nirreps; h++)
 
505
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
506
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
507
          L1a.matrix[h][row][col] = error[0][word++];
 
508
    dpd_file2_mat_wrt(&L1a);
 
509
    dpd_file2_mat_close(&L1a);
 
510
    dpd_file2_close(&L1a);
 
511
 
 
512
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New Lia");
 
513
    dpd_file2_mat_init(&L1a);
 
514
    for(h=0; h < nirreps; h++)
 
515
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
516
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
517
          L1a.matrix[h][row][col] = error[0][word++];
 
518
    dpd_file2_mat_wrt(&L1a);
 
519
    dpd_file2_mat_close(&L1a);
 
520
    dpd_file2_close(&L1a);
 
521
  
 
522
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
523
    for(h=0; h < nirreps; h++) {
 
524
      dpd_buf4_mat_irrep_init(&L2a, h);
 
525
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
526
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
527
          L2a.matrix[h][row][col] = error[0][word++];
 
528
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
529
      dpd_buf4_mat_irrep_close(&L2a, h);
 
530
    }
 
531
    dpd_buf4_close(&L2a);
 
532
 
 
533
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New Lijab");
 
534
    for(h=0; h < nirreps; h++) {
 
535
      dpd_buf4_mat_irrep_init(&L2a, h);
 
536
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
537
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
538
          L2a.matrix[h][row][col] = error[0][word++];
 
539
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
540
      dpd_buf4_mat_irrep_close(&L2a, h);
 
541
    }
 
542
    dpd_buf4_close(&L2a);
 
543
 
 
544
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
 
545
    for(h=0; h < nirreps; h++) {
 
546
      dpd_buf4_mat_irrep_init(&L2a, h);
 
547
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
548
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
549
          L2a.matrix[h][row][col] = error[0][word++];
 
550
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
551
      dpd_buf4_mat_irrep_close(&L2a, h);
 
552
    }
 
553
    dpd_buf4_close(&L2a);
 
554
 
 
555
    /* Release memory and return */
 
556
    free_matrix(vector, nvector);
 
557
    free_matrix(B, nvector+1);
 
558
    free(C);
 
559
    dpd_free_block(error, 1, vector_length);
 
560
  } /** ROHF **/
 
561
  else if(params.ref == 2) { /** UHF **/
 
562
  
 
563
    /* Compute the length of a single error vector */
 
564
    dpd_file2_init(&L1a, CC_TMP0, L_irr, 0, 1, "LIA");
 
565
    dpd_file2_init(&L1b, CC_TMP0, L_irr, 2, 3, "Lia");
 
566
    dpd_buf4_init(&L2a, CC_TMP0, L_irr, 2, 7, 2, 7, 0, "LIJAB");
 
567
    dpd_buf4_init(&L2b, CC_TMP0, L_irr, 12, 17, 12, 17, 0, "Lijab");
 
568
    dpd_buf4_init(&L2c, CC_TMP0, L_irr, 22, 28, 22, 28, 0, "LIjAb");
 
569
    for(h=0; h < nirreps; h++) {
 
570
      vector_length += L1a.params->rowtot[h] * L1a.params->coltot[h^L_irr];
 
571
      vector_length += L1b.params->rowtot[h] * L1b.params->coltot[h^L_irr];
 
572
      vector_length += L2a.params->rowtot[h] * L2a.params->coltot[h^L_irr];
 
573
      vector_length += L2b.params->rowtot[h] * L2b.params->coltot[h^L_irr];
 
574
      vector_length += L2c.params->rowtot[h] * L2c.params->coltot[h^L_irr];
 
575
    }
 
576
    dpd_file2_close(&L1a);
 
577
    dpd_file2_close(&L1b);
 
578
    dpd_buf4_close(&L2a);
 
579
    dpd_buf4_close(&L2b);
 
580
    dpd_buf4_close(&L2c);
 
581
 
 
582
    /* Set the diis cycle value */
 
583
    diis_cycle = (iter-1) % nvector;
 
584
 
 
585
    /* Build the current error vector and dump it to disk */
 
586
    error = dpd_block_matrix(1,vector_length);
 
587
    word=0;
 
588
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
589
    dpd_file2_mat_init(&L1a);
 
590
    dpd_file2_mat_rd(&L1a);
 
591
    dpd_file2_init(&L1b, CC_LAMBDA, L_irr, 0, 1, "LIA");
 
592
    dpd_file2_mat_init(&L1b);
 
593
    dpd_file2_mat_rd(&L1b);
 
594
    for(h=0; h < nirreps; h++)
 
595
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
596
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
597
          error[0][word++] = L1a.matrix[h][row][col] - L1b.matrix[h][row][col];
 
598
    dpd_file2_mat_close(&L1a);
 
599
    dpd_file2_close(&L1a);
 
600
    dpd_file2_mat_close(&L1b);
 
601
    dpd_file2_close(&L1b);
 
602
 
 
603
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 2, 3, "New Lia");
 
604
    dpd_file2_mat_init(&L1a);
 
605
    dpd_file2_mat_rd(&L1a);
 
606
    dpd_file2_init(&L1b, CC_LAMBDA, L_irr, 2, 3, "Lia");
 
607
    dpd_file2_mat_init(&L1b);
 
608
    dpd_file2_mat_rd(&L1b);
 
609
    for(h=0; h < nirreps; h++)
 
610
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
611
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
612
          error[0][word++] = L1a.matrix[h][row][col] - L1b.matrix[h][row][col];
 
613
    dpd_file2_mat_close(&L1a);
 
614
    dpd_file2_close(&L1a);
 
615
    dpd_file2_mat_close(&L1b);
 
616
    dpd_file2_close(&L1b);
 
617
  
 
618
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
619
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
 
620
    for(h=0; h < nirreps; h++) {
 
621
      dpd_buf4_mat_irrep_init(&L2a, h);
 
622
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
623
      dpd_buf4_mat_irrep_init(&L2b, h);
 
624
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
625
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
626
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
627
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
628
      dpd_buf4_mat_irrep_close(&L2a, h);
 
629
      dpd_buf4_mat_irrep_close(&L2b, h);
 
630
    }
 
631
    dpd_buf4_close(&L2a);
 
632
    dpd_buf4_close(&L2b);
 
633
 
 
634
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
 
635
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "Lijab");
 
636
    for(h=0; h < nirreps; h++) {
 
637
      dpd_buf4_mat_irrep_init(&L2a, h);
 
638
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
639
      dpd_buf4_mat_irrep_init(&L2b, h);
 
640
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
641
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
642
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
643
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
644
      dpd_buf4_mat_irrep_close(&L2a, h);
 
645
      dpd_buf4_mat_irrep_close(&L2b, h);
 
646
    }
 
647
    dpd_buf4_close(&L2a);
 
648
    dpd_buf4_close(&L2b);
 
649
 
 
650
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
 
651
    dpd_buf4_init(&L2b, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
 
652
    for(h=0; h < nirreps; h++) {
 
653
      dpd_buf4_mat_irrep_init(&L2a, h);
 
654
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
655
      dpd_buf4_mat_irrep_init(&L2b, h);
 
656
      dpd_buf4_mat_irrep_rd(&L2b, h);
 
657
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
658
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
659
          error[0][word++] = L2a.matrix[h][row][col] - L2b.matrix[h][row][col];
 
660
      dpd_buf4_mat_irrep_close(&L2a, h);
 
661
      dpd_buf4_mat_irrep_close(&L2b, h);
 
662
    }
 
663
    dpd_buf4_close(&L2a);
 
664
    dpd_buf4_close(&L2b);
 
665
 
 
666
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
667
    psio_write(CC_DIIS_ERR, "DIIS Error[0] Vectors" , (char *) error[0], 
 
668
               vector_length*sizeof(double), start, &end);
 
669
 
 
670
    /* Store the current amplitude vector on disk */
 
671
    word=0;
 
672
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
673
    dpd_file2_mat_init(&L1a);
 
674
    dpd_file2_mat_rd(&L1a);
 
675
    for(h=0; h < nirreps; h++)
 
676
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
677
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
678
          error[0][word++] = L1a.matrix[h][row][col];
 
679
    dpd_file2_mat_close(&L1a);
 
680
    dpd_file2_close(&L1a);
 
681
 
 
682
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 2, 3, "New Lia");
 
683
    dpd_file2_mat_init(&L1a);
 
684
    dpd_file2_mat_rd(&L1a);
 
685
    for(h=0; h < nirreps; h++)
 
686
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
687
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
688
          error[0][word++] = L1a.matrix[h][row][col];
 
689
    dpd_file2_mat_close(&L1a);
 
690
    dpd_file2_close(&L1a);
 
691
  
 
692
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
693
    for(h=0; h < nirreps; h++) {
 
694
      dpd_buf4_mat_irrep_init(&L2a, h);
 
695
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
696
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
697
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
698
          error[0][word++] = L2a.matrix[h][row][col];
 
699
      dpd_buf4_mat_irrep_close(&L2a, h);
 
700
    }
 
701
    dpd_buf4_close(&L2a);
 
702
 
 
703
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
 
704
    for(h=0; h < nirreps; h++) {
 
705
      dpd_buf4_mat_irrep_init(&L2a, h);
 
706
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
707
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
708
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
709
          error[0][word++] = L2a.matrix[h][row][col];
 
710
      dpd_buf4_mat_irrep_close(&L2a, h);
 
711
    }
 
712
    dpd_buf4_close(&L2a);
 
713
 
 
714
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
 
715
    for(h=0; h < nirreps; h++) {
 
716
      dpd_buf4_mat_irrep_init(&L2a, h);
 
717
      dpd_buf4_mat_irrep_rd(&L2a, h);
 
718
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
719
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
720
          error[0][word++] = L2a.matrix[h][row][col];
 
721
      dpd_buf4_mat_irrep_close(&L2a, h);
 
722
    }
 
723
    dpd_buf4_close(&L2a);
 
724
 
 
725
    start = psio_get_address(PSIO_ZERO, diis_cycle*vector_length*sizeof(double));
 
726
    psio_write(CC_DIIS_AMP, "DIIS Amplitude Vectors" , (char *) error[0], 
 
727
               vector_length*sizeof(double), start, &end);
 
728
 
 
729
    /* If we haven't run through enough iterations, set the correct dimensions
 
730
       for the extrapolation */
 
731
    if(!(iter >= (nvector))) {
 
732
      if(iter < 2) { /* Leave if we can't extrapolate at all */
 
733
        free(error[0]);
 
734
        return; 
 
735
      }
 
736
      nvector = iter;
 
737
    }
 
738
 
 
739
    /* Now grab the full set of error[0] vectors from the file */
 
740
    vector = init_matrix(nvector, vector_length);
 
741
    next = PSIO_ZERO;
 
742
    for(p=0; p < nvector; p++) 
 
743
      psio_read(CC_DIIS_ERR, "DIIS Error[0] Vectors", (char *) vector[p], 
 
744
                vector_length*sizeof(double), next, &next);
 
745
 
 
746
    /* Build B matrix of error[0] vector products */
 
747
    B = init_matrix(nvector+1,nvector+1);
 
748
 
 
749
    for(p=0; p < nvector; p++)
 
750
      for(q=0; q < nvector; q++) {
 
751
        dot_arr(vector[p], vector[q], vector_length, &product); 
 
752
        B[p][q] = product;
 
753
      }
 
754
 
 
755
    for(p=0; p < nvector; p++) {
 
756
      B[p][nvector] = -1;
 
757
      B[nvector][p] = -1;
 
758
    }
 
759
 
 
760
    B[nvector][nvector] = 0;
 
761
 
 
762
    /* Find the maximum value in B and scale all its elements */
 
763
    maximum = fabs(B[0][0]);
 
764
    for(p=0; p < nvector; p++)
 
765
      for(q=0; q < nvector; q++)
 
766
        if(fabs(B[p][q]) > maximum) maximum = fabs(B[p][q]);
 
767
 
 
768
    for(p=0; p < nvector; p++)
 
769
      for(q=0; q < nvector; q++)
 
770
        B[p][q] /= maximum; 
 
771
 
 
772
    /* Build the constant vector */
 
773
    C = init_array(nvector+1);
 
774
    C[nvector] = -1;
 
775
 
 
776
    /* Solve the linear equations */
 
777
    flin(B, C, nvector+1, 1, &determinant);
 
778
 
 
779
    /* Grab the old amplitude vectors */
 
780
    next = PSIO_ZERO;
 
781
    for(p=0; p < nvector; p++) 
 
782
      psio_read(CC_DIIS_AMP, "DIIS Amplitude Vectors", (char *) vector[p], 
 
783
                vector_length*sizeof(double), next, &next);
 
784
  
 
785
    /* Build the new amplitude vector from the old ones */
 
786
    for(q=0; q < vector_length; q++) {
 
787
      error[0][q] = 0.0;
 
788
      for(p=0; p < nvector; p++)
 
789
        error[0][q] += C[p] * vector[p][q];
 
790
    }
 
791
 
 
792
    /* Now place these elements into the DPD amplitude arrays */
 
793
    word=0;
 
794
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 0, 1, "New LIA");
 
795
    dpd_file2_mat_init(&L1a);
 
796
    for(h=0; h < nirreps; h++)
 
797
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
798
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
799
          L1a.matrix[h][row][col] = error[0][word++];
 
800
    dpd_file2_mat_wrt(&L1a);
 
801
    dpd_file2_mat_close(&L1a);
 
802
    dpd_file2_close(&L1a);
 
803
 
 
804
    dpd_file2_init(&L1a, CC_LAMBDA, L_irr, 2, 3, "New Lia");
 
805
    dpd_file2_mat_init(&L1a);
 
806
    for(h=0; h < nirreps; h++)
 
807
      for(row=0; row < L1a.params->rowtot[h]; row++)
 
808
        for(col=0; col < L1a.params->coltot[h^L_irr]; col++)
 
809
          L1a.matrix[h][row][col] = error[0][word++];
 
810
    dpd_file2_mat_wrt(&L1a);
 
811
    dpd_file2_mat_close(&L1a);
 
812
    dpd_file2_close(&L1a);
 
813
  
 
814
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
 
815
    for(h=0; h < nirreps; h++) {
 
816
      dpd_buf4_mat_irrep_init(&L2a, h);
 
817
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
818
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
819
          L2a.matrix[h][row][col] = error[0][word++];
 
820
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
821
      dpd_buf4_mat_irrep_close(&L2a, h);
 
822
    }
 
823
    dpd_buf4_close(&L2a);
 
824
 
 
825
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
 
826
    for(h=0; h < nirreps; h++) {
 
827
      dpd_buf4_mat_irrep_init(&L2a, h);
 
828
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
829
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
830
          L2a.matrix[h][row][col] = error[0][word++];
 
831
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
832
      dpd_buf4_mat_irrep_close(&L2a, h);
 
833
    }
 
834
    dpd_buf4_close(&L2a);
 
835
 
 
836
    dpd_buf4_init(&L2a, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
 
837
    for(h=0; h < nirreps; h++) {
 
838
      dpd_buf4_mat_irrep_init(&L2a, h);
 
839
      for(row=0; row < L2a.params->rowtot[h]; row++)
 
840
        for(col=0; col < L2a.params->coltot[h^L_irr]; col++)
 
841
          L2a.matrix[h][row][col] = error[0][word++];
 
842
      dpd_buf4_mat_irrep_wrt(&L2a, h);
 
843
      dpd_buf4_mat_irrep_close(&L2a, h);
 
844
    }
 
845
    dpd_buf4_close(&L2a);
 
846
 
 
847
    /* Release memory and return */
 
848
    free_matrix(vector, nvector);
 
849
    free_matrix(B, nvector+1);
 
850
    free(C);
 
851
    dpd_free_block(error, 1, vector_length);
 
852
  } /** UHF **/
 
853
 
 
854
  return;
 
855
}
 
856
 
 
857
}} // namespace psi::cclambda