← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cceom/diagSS.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/diagSS.cc
 
1
/*! \file
 
2
    \ingroup CCEOM
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdlib>
 
6
#include <cstdio>
 
7
#include <cstring>
 
8
#include <cmath>
 
9
#include <libciomr/libciomr.h>
 
10
#include <libqt/qt.h>
 
11
#include "MOInfo.h"
 
12
#include "Params.h"
 
13
#include "Local.h"
 
14
#define EXTERN
 
15
#include "globals.h"
 
16
 
 
17
namespace psi { namespace cceom {
 
18
#include <physconst.h>
 
19
 
 
20
void sigmaSS(int index, int C_irr);
 
21
void cc2_sigmaSS(int index, int C_irr);
 
22
void precondition_SS(dpdfile2 *RIA, dpdfile2 *Ria, double eval);
 
23
void schmidt_add_SS(dpdfile2 *RIA, dpdfile2 *Ria, int C_irr, int *numCs);
 
24
void precondition_SS_RHF(dpdfile2 *RIA, double eval);
 
25
void schmidt_add_SS_RHF(dpdfile2 *RIA, int C_irr, int *numCs);
 
26
void restart_SS(double **alpha, int L, int num, int C_irr);
 
27
void dgeev_eom(int L, double **G, double *evals, double **alpha);
 
28
double norm_C1(dpdfile2 *C1A, dpdfile2 *C1B);
 
29
double norm_C1_rhf(dpdfile2 *C1A);
 
30
double scm_C1(dpdfile2 *C1A, dpdfile2 *C1B, double a);
 
31
 
 
32
void diagSS(int C_irr) {
 
33
  dpdfile2 Fmi, FMI, Fae, FAE, Fme, FME;
 
34
  dpdfile2 CME, Cme, C, SIA, Sia, RIA, Ria;
 
35
  dpdbuf4 CMNEF, Cmnef, CMnEf, W;
 
36
  char lbl[32], lbl2[32];
 
37
  int lwork, info, get_right_ev = 1, get_left_ev = 0;
 
38
  int L,h,i,j,k,a,C_index,errcod,keep_going=1,numCs,iter=0;
 
39
  double norm, tval, *lambda, *lambda_old, zero=0.0;
 
40
  double **G, *work, *evals_complex, **alpha, **evectors_left;
 
41
  int nirreps, *openpi, *occpi, *virtpi, range;
 
42
  int *aoccpi, *avirtpi, *boccpi, *bvirtpi;
 
43
  int *converged, num_converged, num_roots;
 
44
  int begin_occ, begin_virt, end_virt, dim_SS = 0;
 
45
  int pf, cnt, irr_occ, irr_virt;
 
46
 
 
47
  nirreps = moinfo.nirreps;
 
48
  openpi = moinfo.openpi;
 
49
  occpi = moinfo.occpi;
 
50
  virtpi = moinfo.virtpi;
 
51
 
 
52
  if (params.eom_ref == 2) { /* UHF */
 
53
    aoccpi = moinfo.aoccpi;
 
54
    boccpi = moinfo.boccpi;
 
55
    avirtpi = moinfo.avirtpi;
 
56
    bvirtpi = moinfo.bvirtpi;
 
57
  }
 
58
 
 
59
  range = eom_params.excitation_range;
 
60
  pf = eom_params.print_singles;
 
61
  if (pf) fprintf(outfile,"\n\n");
 
62
 
 
63
  /* a bunch of tedious code to setup reasonable HOMO-LUMO guess vectors */
 
64
  C_index=0;
 
65
  if (2 * range * range * nirreps < eom_params.cs_per_irrep[C_irr])
 
66
    range = eom_params.cs_per_irrep[C_irr] / 2;
 
67
 
 
68
  if(!params.local) { /* guesses should already be built for local mode */
 
69
 
 
70
    for (cnt=0; cnt<nirreps; ++cnt) { /* cnt loops over dp's to get C_irr */
 
71
      irr_occ = dpd_dp[C_irr][cnt][0];
 
72
      irr_virt = dpd_dp[C_irr][cnt][1];
 
73
      /* C_irr = irr_occ * irr_virt */
 
74
 
 
75
      if (params.eom_ref == 0)  { /* ref = RHF; eom_ref = RHF */
 
76
        begin_occ = MAX(occpi[irr_occ]-range, 0);
 
77
        end_virt = MIN(range, virtpi[irr_virt]);
 
78
        for (i=begin_occ; i < occpi[irr_occ]; ++i)
 
79
          for (a=0; a < end_virt; ++a, ++C_index) {
 
80
            sprintf(lbl, "%s %d", "CME", C_index);
 
81
            dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
82
            dpd_file2_mat_init(&CME);
 
83
            CME.matrix[cnt][i][a] = 1.0/sqrt(2.0);
 
84
            dpd_file2_mat_wrt(&CME);
 
85
            dpd_file2_mat_close(&CME);
 
86
            dpd_file2_close(&CME);
 
87
          }
 
88
      }
 
89
      /* eom_ref = ROHF, closed shell */
 
90
      else if ((params.eom_ref < 2) && (moinfo.iopen == 0)) {
 
91
        begin_occ = MAX(occpi[irr_occ]-range, 0);
 
92
        end_virt = MIN(range, virtpi[irr_virt]);
 
93
        for (i=begin_occ; i < occpi[irr_occ]; ++i)
 
94
          for (a=0; a < end_virt; ++a, ++C_index) {
 
95
            sprintf(lbl, "%s %d", "CME", C_index);
 
96
            dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
97
            dpd_file2_mat_init(&CME);
 
98
            CME.matrix[cnt][i][a] = 1.0/sqrt(2.0);
 
99
            dpd_file2_mat_wrt(&CME);
 
100
            dpd_file2_mat_close(&CME);
 
101
            dpd_file2_close(&CME);
 
102
            sprintf(lbl, "%s %d", "Cme", C_index);
 
103
            dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
104
            dpd_file2_mat_init(&Cme);
 
105
            Cme.matrix[cnt][i][a] = 1.0/sqrt(2.0);
 
106
            dpd_file2_mat_wrt(&Cme);
 
107
            dpd_file2_mat_close(&Cme);
 
108
            dpd_file2_close(&Cme);
 
109
            ++C_index;
 
110
            sprintf(lbl, "%s %d", "CME", C_index);
 
111
            dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
112
            dpd_file2_mat_init(&CME);
 
113
            CME.matrix[cnt][i][a] = 1.0/sqrt(2.0);
 
114
            dpd_file2_mat_wrt(&CME);
 
115
            dpd_file2_mat_close(&CME);
 
116
            dpd_file2_close(&CME);
 
117
            sprintf(lbl, "%s %d", "Cme", C_index);
 
118
            dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
119
            dpd_file2_mat_init(&Cme);
 
120
            Cme.matrix[cnt][i][a] = -1.0/sqrt(2.0);
 
121
            dpd_file2_mat_wrt(&Cme);
 
122
            dpd_file2_mat_close(&Cme);
 
123
            dpd_file2_close(&Cme);
 
124
         }
 
125
      }
 
126
      else if (params.eom_ref == 1) { /* open-shell ROHF */
 
127
      /* alpha excitations */
 
128
      begin_occ = MAX(occpi[irr_occ]-range, 0);
 
129
      end_virt = MIN( virtpi[irr_virt]-openpi[irr_virt], range);
 
130
      for (i=begin_occ; i < occpi[irr_occ] ; ++i)
 
131
        for (a=0; a<end_virt; ++a, ++C_index) {
 
132
          sprintf(lbl, "%s %d", "CME", C_index);
 
133
          dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
134
          dpd_file2_mat_init(&CME);
 
135
          CME.matrix[cnt][i][a] = 1.0;
 
136
          dpd_file2_mat_wrt(&CME);
 
137
          dpd_file2_close(&CME);
 
138
          sprintf(lbl, "%s %d", "Cme", C_index);
 
139
          dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
140
          dpd_file2_mat_init(&Cme);
 
141
          dpd_file2_mat_wrt(&Cme);
 
142
          dpd_file2_close(&Cme);
 
143
        }
 
144
      /* beta excitations into open shells */
 
145
      begin_occ = MAX(occpi[irr_occ]-openpi[irr_occ]-range, 0);
 
146
      begin_virt = virtpi[irr_virt] - openpi[irr_virt];
 
147
      for (i=begin_occ; i < occpi[irr_occ]-openpi[irr_occ]; ++i)
 
148
        for (a=begin_virt; a < virtpi[irr_virt]; ++a, ++C_index) {
 
149
          sprintf(lbl, "%s %d", "Cme", C_index);
 
150
          dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
151
          dpd_file2_mat_init(&Cme);
 
152
          Cme.matrix[cnt][i][a] = 1.0;
 
153
          dpd_file2_mat_wrt(&Cme);
 
154
          dpd_file2_close(&Cme);
 
155
          sprintf(lbl, "%s %d", "CME", C_index);
 
156
          dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
157
          dpd_file2_mat_init(&CME);
 
158
          dpd_file2_mat_wrt(&CME);
 
159
          dpd_file2_close(&CME);
 
160
        }
 
161
      /* beta excitations into unoccupied orbitals */
 
162
      begin_occ = MAX(occpi[irr_occ]-openpi[irr_occ]-range, 0);
 
163
      end_virt = MIN(range - openpi[irr_virt], 0);
 
164
      for (i=begin_occ; i < occpi[irr_occ]-openpi[irr_occ]; ++i)
 
165
        for (a=0; a < end_virt; ++a, ++C_index) {
 
166
          sprintf(lbl, "%s %d", "Cme", C_index);
 
167
          dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
168
          dpd_file2_mat_init(&Cme);
 
169
          Cme.matrix[cnt][i][a] = 1.0;
 
170
          dpd_file2_mat_wrt(&Cme);
 
171
          dpd_file2_close(&Cme);
 
172
          sprintf(lbl, "%s %d", "CME", C_index);
 
173
          dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
174
          dpd_file2_mat_init(&CME);
 
175
          dpd_file2_mat_wrt(&CME);
 
176
          dpd_file2_close(&CME);
 
177
        }
 
178
      }
 
179
      else { /* UHF */
 
180
        /* alpha excitations */
 
181
        begin_occ = MAX(aoccpi[irr_occ]-range, 0);
 
182
        end_virt = MIN(avirtpi[irr_virt], range);
 
183
        for (i=begin_occ; i < aoccpi[irr_occ] ; ++i)
 
184
          for (a=0; a<end_virt; ++a, ++C_index) {
 
185
            sprintf(lbl, "%s %d", "CME", C_index);
 
186
            dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
187
            dpd_file2_mat_init(&CME);
 
188
            CME.matrix[cnt][i][a] = 1.0;
 
189
            dpd_file2_mat_wrt(&CME);
 
190
            dpd_file2_close(&CME);
 
191
            sprintf(lbl, "%s %d", "Cme", C_index);
 
192
            dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
193
            dpd_file2_mat_init(&Cme);
 
194
            dpd_file2_mat_wrt(&Cme);
 
195
            dpd_file2_close(&Cme);
 
196
          }
 
197
        begin_occ = MAX(boccpi[irr_occ]-range, 0);
 
198
        end_virt = MIN(bvirtpi[irr_virt], range);
 
199
        for (i=begin_occ; i < boccpi[irr_occ] ; ++i)
 
200
          for (a=0; a<end_virt; ++a, ++C_index) {
 
201
            sprintf(lbl, "%s %d", "CME", C_index);
 
202
            dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
203
            dpd_file2_mat_init(&CME);
 
204
            dpd_file2_mat_wrt(&CME);
 
205
            dpd_file2_close(&CME);
 
206
            sprintf(lbl, "%s %d", "Cme", C_index);
 
207
            dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
208
            dpd_file2_mat_init(&Cme);
 
209
            Cme.matrix[cnt][i][a] = 1.0;
 
210
            dpd_file2_mat_wrt(&Cme);
 
211
            dpd_file2_close(&Cme);
 
212
          }
 
213
      }
 
214
    }
 
215
 
 
216
   if (pf) fprintf(outfile,"%d initial single excitation guesses\n",C_index);
 
217
   if (C_index == 0) {
 
218
      fprintf(outfile, "No intial guesses obtained for %s state \n",
 
219
              moinfo.labels[moinfo.sym^C_irr]);
 
220
      exit(1);
 
221
    }
 
222
  }   
 
223
  else {
 
224
    C_index = eom_params.cs_per_irrep[0];
 
225
  } /* end if(!params.local) */
 
226
 
 
227
 
 
228
 
 
229
  /* show initial guesses */
 
230
  /*
 
231
  for (i=0; i<C_index ;++i) {
 
232
    sprintf(lbl, "%s %d", "CME", i);
 
233
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
234
    dpd_file2_print(&CME, outfile);
 
235
    dpd_file2_close(&CME);
 
236
    if (params.eom_ref > 0) { 
 
237
      sprintf(lbl, "%s %d", "Cme", i);
 
238
      if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
239
      else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
240
      dpd_file2_print(&Cme, outfile);
 
241
      dpd_file2_close(&Cme);
 
242
    }
 
243
  }
 
244
  */
 
245
  if (eom_params.skip_diagSS) {
 
246
    fprintf(outfile,"\nSkipping diagonalization of Hbar SS block.\n");
 
247
    return;
 
248
  }
 
249
 
 
250
  fflush(outfile);
 
251
 
 
252
  /* Setup residual vector file */
 
253
  dpd_file2_init(&RIA, EOM_R, C_irr, 0, 1, "RIA");
 
254
  dpd_file2_mat_init(&RIA);
 
255
  dpd_file2_mat_wrt(&RIA);
 
256
  dpd_file2_close(&RIA);
 
257
  if (params.eom_ref > 0) {
 
258
    if (params.eom_ref == 1) dpd_file2_init(&Ria, EOM_R, C_irr, 0, 1, "Ria");
 
259
    else if (params.eom_ref == 2) dpd_file2_init(&Ria, EOM_R, C_irr, 2, 3, "Ria");
 
260
    dpd_file2_mat_init(&Ria);
 
261
    dpd_file2_mat_wrt(&Ria);
 
262
    dpd_file2_close(&Ria);
 
263
  }
 
264
 
 
265
  /* arrays must be dimensioned with at least the final number of roots - even though
 
266
     num_roots may be limited until the first collapse by the number of good
 
267
     initial guesses obtained above. */
 
268
  L = num_roots = C_index;
 
269
  i = MAX(eom_params.cs_per_irrep[C_irr], C_index);
 
270
  converged = init_int_array(i);
 
271
  lambda_old = init_array(i);
 
272
 
 
273
  if (pf) fprintf(outfile,"\n");
 
274
  while ((keep_going == 1) && (iter < eom_params.max_iter_SS)) {
 
275
    if (pf) fprintf(outfile,"Iter=%-4d L=%-4d", iter+1, L);
 
276
    keep_going = 0;
 
277
    numCs = L;
 
278
    num_converged = 0;
 
279
 
 
280
    /* zero S vectors */
 
281
    for (i=0;i<L;++i) {
 
282
      if (params.eom_ref == 0) {
 
283
        sprintf(lbl, "%s %d", "SIA", i);
 
284
        dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
285
        dpd_file2_scm(&SIA, 0.0);
 
286
        dpd_file2_close(&SIA);
 
287
                                if (params.full_matrix) {
 
288
          sprintf(lbl, "%s %d", "S0", i);
 
289
                                  psio_write_entry(EOM_SIA, lbl, (char *) &zero, sizeof(double));
 
290
                          }
 
291
      }
 
292
      if (params.eom_ref > 0) {
 
293
        sprintf(lbl, "%s %d", "SIA", i);
 
294
        dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
295
        sprintf(lbl, "%s %d", "Sia", i);
 
296
        if (params.eom_ref == 1) dpd_file2_init(&Sia, EOM_Sia, C_irr, 0, 1, lbl);
 
297
        else if (params.eom_ref == 2) dpd_file2_init(&Sia, EOM_Sia, C_irr, 2, 3, lbl);
 
298
        scm_C1(&SIA, &Sia, 0.0);
 
299
        dpd_file2_close(&SIA);
 
300
        dpd_file2_close(&Sia);
 
301
      }
 
302
    }
 
303
 
 
304
    if(!strcmp(params.wfn,"EOM_CC2")) {
 
305
      for (i=0;i<L;++i)
 
306
        cc2_sigmaSS(i,C_irr);
 
307
    }
 
308
    else {
 
309
      for (i=0;i<L;++i)
 
310
        sigmaSS(i,C_irr);
 
311
    }
 
312
 
 
313
    G = block_matrix(L,L);
 
314
 
 
315
    for (i=0;i<L;++i) {
 
316
      sprintf(lbl, "%s %d", "CME", i);
 
317
      dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
318
      if (params.eom_ref > 0) {
 
319
        sprintf(lbl, "%s %d", "Cme", i);
 
320
        if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
321
        else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
322
      }
 
323
      for (j=0;j<L;++j) {
 
324
        if(params.eom_ref == 0) {
 
325
          sprintf(lbl, "%s %d", "SIA", j);
 
326
          dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
327
          tval = 2.0 * dpd_file2_dot(&CME, &SIA);
 
328
          dpd_file2_close(&SIA);
 
329
        }
 
330
        else if (params.eom_ref > 0) {
 
331
          sprintf(lbl, "%s %d", "SIA", j);
 
332
          dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
333
          sprintf(lbl, "%s %d", "Sia", j);
 
334
          if (params.eom_ref == 1) dpd_file2_init(&Sia, EOM_Sia, C_irr, 0, 1, lbl);
 
335
          else if (params.eom_ref == 2) dpd_file2_init(&Sia, EOM_Sia, C_irr, 2, 3, lbl);
 
336
          tval = dpd_file2_dot(&CME, &SIA);
 
337
          tval += dpd_file2_dot(&Cme, &Sia);
 
338
          dpd_file2_close(&SIA);
 
339
          dpd_file2_close(&Sia);
 
340
        }
 
341
        G[i][j] = tval;
 
342
      }
 
343
      dpd_file2_close(&CME);
 
344
      if (params.eom_ref > 0) dpd_file2_close(&Cme);
 
345
    }
 
346
 
 
347
    /* print_mat(G, L, L, outfile); */
 
348
 
 
349
    lambda = init_array(L);        /* holds real part of eigenvalues of G */
 
350
    alpha = block_matrix(L,L);     /* will hold eigenvectors of G */
 
351
    dgeev_eom(L, G, lambda, alpha);
 
352
    eigsort(lambda, alpha, L);
 
353
    /*  eivout(alpha, lambda, L, L, outfile); */
 
354
    free_block(G);
 
355
 
 
356
    if (pf) fprintf(outfile,
 
357
                    "  Root    EOM Energy    Delta E     Res. Norm    Conv?\n");
 
358
    fflush(outfile);
 
359
 
 
360
    dpd_file2_init(&RIA, EOM_R, C_irr, 0, 1, "RIA");
 
361
    if (params.eom_ref == 1) dpd_file2_init(&Ria, EOM_R, C_irr, 0, 1, "Ria");
 
362
    else if (params.eom_ref == 2) dpd_file2_init(&Ria, EOM_R, C_irr, 2, 3, "Ria");
 
363
 
 
364
    for (k=0;k<num_roots;++k) {
 
365
      dpd_file2_scm(&RIA, 0.0);
 
366
      if (params.eom_ref > 0) dpd_file2_scm(&Ria, 0.0);
 
367
      converged[k] = 0;
 
368
      for (i=0;i<L;++i) { 
 
369
        sprintf(lbl, "%s %d", "SIA", i);
 
370
        dpd_file2_init(&SIA, EOM_SIA, C_irr, 0, 1, lbl);
 
371
        sprintf(lbl, "%s %d", "CME", i);
 
372
        dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
373
        dpd_file2_axpbycz(&CME, &SIA, &RIA, -1.0*lambda[k]*alpha[i][k], alpha[i][k], 1.0);
 
374
        dpd_file2_close(&CME);
 
375
        dpd_file2_close(&SIA);
 
376
        if (params.eom_ref > 0) {
 
377
          sprintf(lbl, "%s %d", "Sia", i);
 
378
          if (params.eom_ref == 1) dpd_file2_init(&Sia, EOM_Sia, C_irr, 0, 1, lbl);
 
379
          else if (params.eom_ref == 2) dpd_file2_init(&Sia, EOM_Sia, C_irr, 2, 3, lbl);
 
380
          sprintf(lbl, "%s %d", "Cme", i);
 
381
          if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
382
          else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
383
          dpd_file2_axpbycz(&Cme, &Sia, &Ria, -1.0*lambda[k]*alpha[i][k], alpha[i][k], 1.0);
 
384
          dpd_file2_close(&Cme);
 
385
          dpd_file2_close(&Sia);
 
386
        }
 
387
      }
 
388
 
 
389
      if (params.eom_ref == 0) precondition_SS_RHF(&RIA, lambda[k]);
 
390
      else precondition_SS(&RIA, &Ria, lambda[k]);
 
391
 
 
392
      if (params.eom_ref == 0)
 
393
        norm = norm_C1_rhf(&RIA);
 
394
      else 
 
395
        norm = norm_C1(&RIA, &Ria);
 
396
 
 
397
      if (pf) fprintf(outfile,"%6d%15.10lf%11.2e%12.2e",k+1,lambda[k],
 
398
                      lambda[k]-lambda_old[k], norm); 
 
399
 
 
400
      if ( (norm > eom_params.residual_tol_SS) ||
 
401
           (fabs(lambda[k]-lambda_old[k]) > eom_params.eval_tol_SS) ) {
 
402
        if (pf) fprintf(outfile,"%7s\n","N");
 
403
        /*
 
404
          if (params.eom_ref == 0) precondition_SS_RHF(&RIA, lambda[k]);
 
405
          else precondition_SS(&RIA, &Ria, lambda[k]);
 
406
        */
 
407
 
 
408
        /* normalize preconditioned residual */
 
409
        if (params.eom_ref == 0) {
 
410
          norm = norm_C1_rhf(&RIA);
 
411
          dpd_file2_scm(&RIA, 1.0/norm);
 
412
        }
 
413
        else {
 
414
          norm = norm_C1(&RIA, &Ria);
 
415
          dpd_file2_scm(&RIA, 1.0/norm);
 
416
          dpd_file2_scm(&Ria, 1.0/norm);
 
417
        }
 
418
 
 
419
        if(params.eom_ref == 0) schmidt_add_SS_RHF(&RIA, C_irr, &numCs);
 
420
        else schmidt_add_SS(&RIA, &Ria, C_irr, &numCs);
 
421
      }
 
422
      else {
 
423
        if (pf) fprintf(outfile,"%7s\n","Y");
 
424
        ++num_converged;
 
425
        converged[k] = 1;
 
426
      }
 
427
 
 
428
      fflush(outfile);
 
429
    }
 
430
 
 
431
    dpd_file2_close(&RIA);
 
432
    if (params.eom_ref > 0) dpd_file2_close(&Ria);
 
433
 
 
434
    for (i=0;i<num_roots;++i) lambda_old[i] = lambda[i];
 
435
    free(lambda);
 
436
 
 
437
    if ( (iter == 2) || (L > eom_params.vectors_per_root_SS * eom_params.cs_per_irrep[C_irr])) {
 
438
      restart_SS(alpha, L, eom_params.cs_per_irrep[C_irr], C_irr);
 
439
      L = num_roots = eom_params.cs_per_irrep[C_irr];
 
440
      keep_going = 1;
 
441
      free_block(alpha);
 
442
    }
 
443
    else {
 
444
      if (numCs > L) { /* successfully added vectors */
 
445
        keep_going = 1;
 
446
        free_block(alpha);
 
447
        L = numCs;
 
448
      }
 
449
    }
 
450
    ++iter;
 
451
  }
 
452
 
 
453
  if (pf) fprintf(outfile,"\nLowest eigenvalues of HBar Singles-Singles Block\n");
 
454
  if (pf) fprintf(outfile,"Root      Excitation Energy         Total Energy\n");
 
455
  if (pf) fprintf(outfile,"           (eV)     (cm^-1)            (au)\n");
 
456
  if (pf) for (i=0;i<num_roots;++i)
 
457
    if (pf)   fprintf(outfile,"%4d%12.3lf%12.2lf%20.10lf\n",i+1,
 
458
                      lambda_old[i]* _hartree2ev, lambda_old[i]* _hartree2wavenumbers,
 
459
                      lambda_old[i]+moinfo.eref+moinfo.ecc);
 
460
 
 
461
  free(lambda_old);
 
462
  free(converged);
 
463
  /* collapse solutions to one vector each */
 
464
  restart_SS(alpha, L, eom_params.cs_per_irrep[C_irr], C_irr);
 
465
  free_block(alpha);
 
466
  if (pf) fprintf(outfile,"\n");
 
467
  fflush(outfile);
 
468
 
 
469
  return;
 
470
}
 
471
 
 
472
void precondition_SS(dpdfile2 *RIA, dpdfile2 *Ria, double eval)
 
473
{
 
474
  dpdfile2 DIA, Dia;
 
475
  int h, nirreps, i, j, a, b, ij, ab, C_irr;
 
476
  double tval;
 
477
 
 
478
  C_irr = RIA->my_irrep;
 
479
  nirreps = RIA->params->nirreps;
 
480
 
 
481
  dpd_file2_mat_init(RIA);
 
482
  dpd_file2_mat_rd(RIA);
 
483
  dpd_file2_init(&DIA, EOM_D, C_irr, 0, 1, "DIA");
 
484
  dpd_file2_mat_init(&DIA);
 
485
  dpd_file2_mat_rd(&DIA);
 
486
  for(h=0; h < nirreps; h++)
 
487
    for(i=0; i < RIA->params->rowtot[h]; i++)
 
488
      for(a=0; a < RIA->params->coltot[h^C_irr]; a++) {
 
489
        tval = eval - DIA.matrix[h][i][a];
 
490
        if (fabs(tval) > 0.0001) RIA->matrix[h][i][a] /= tval;
 
491
      }
 
492
  dpd_file2_mat_wrt(RIA);
 
493
  dpd_file2_mat_close(RIA);
 
494
  dpd_file2_close(&DIA);
 
495
 
 
496
  dpd_file2_mat_init(Ria);
 
497
  dpd_file2_mat_rd(Ria);
 
498
 
 
499
  if (params.eom_ref == 1) dpd_file2_init(&Dia, EOM_D, C_irr, 0, 1, "Dia");
 
500
  else if (params.eom_ref == 2) dpd_file2_init(&Dia, EOM_D, C_irr, 2, 3, "Dia");
 
501
  dpd_file2_mat_init(&Dia);
 
502
  dpd_file2_mat_rd(&Dia);
 
503
  for(h=0; h < nirreps; h++)
 
504
    for(i=0; i < Ria->params->rowtot[h]; i++)
 
505
      for(a=0; a < Ria->params->coltot[h^C_irr]; a++) {
 
506
        tval = eval - Dia.matrix[h][i][a];
 
507
        if (fabs(tval) > 0.0001) Ria->matrix[h][i][a] /= tval;
 
508
      }
 
509
  dpd_file2_mat_wrt(Ria);
 
510
  dpd_file2_mat_close(Ria);
 
511
  dpd_file2_close(&Dia);
 
512
 
 
513
  return;
 
514
}
 
515
 
 
516
void precondition_SS_RHF(dpdfile2 *RIA, double eval)
 
517
{
 
518
  dpdfile2 DIA;
 
519
  int C_irr, h, nirreps, i, j, a, b, ij, ii, ab;
 
520
  double tval;
 
521
  psio_address next;
 
522
 
 
523
  /* Local correlation decs */
 
524
  int nso, nocc, nvir;
 
525
  double *T1tilde, *T1bar;
 
526
 
 
527
  nirreps = RIA->params->nirreps;
 
528
  C_irr = RIA->my_irrep;
 
529
 
 
530
  if(params.local && local.filter_singles) {
 
531
 
 
532
    nso = local.nso;
 
533
    nocc = local.nocc;
 
534
    nvir = local.nvir;
 
535
 
 
536
    local.pairdom_len = init_int_array(nocc*nocc);
 
537
    local.pairdom_nrlen = init_int_array(nocc*nocc);
 
538
    local.eps_occ = init_array(nocc);
 
539
    psio_read_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
 
540
                    nocc*nocc*sizeof(int));
 
541
    psio_read_entry(CC_INFO, "Local Pair Domain Length (Non-redundant basis)", (char *) local.pairdom_nrlen,
 
542
                    nocc*nocc*sizeof(int));
 
543
    psio_read_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
 
544
                    nocc*sizeof(double));
 
545
    local.W = (double ***) malloc(nocc * nocc * sizeof(double **));
 
546
    local.V = (double ***) malloc(nocc * nocc * sizeof(double **));
 
547
    local.eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
 
548
    next = PSIO_ZERO;
 
549
    for(ij=0; ij < nocc*nocc; ij++) {
 
550
      local.eps_vir[ij] = init_array(local.pairdom_nrlen[ij]);
 
551
      psio_read(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
 
552
                local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
553
    }
 
554
    next = PSIO_ZERO;
 
555
    for(ij=0; ij < nocc*nocc; ij++) {
 
556
      local.V[ij] = block_matrix(nvir,local.pairdom_len[ij]);
 
557
      psio_read(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
 
558
                nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
 
559
    }
 
560
    next = PSIO_ZERO;
 
561
    for(ij=0; ij < nocc*nocc; ij++) {
 
562
      local.W[ij] = block_matrix(local.pairdom_len[ij],local.pairdom_nrlen[ij]);
 
563
      psio_read(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
 
564
                local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
565
    }
 
566
 
 
567
    dpd_file2_mat_init(RIA);
 
568
    dpd_file2_mat_rd(RIA);
 
569
 
 
570
    for(i=0; i < nocc; i++) {
 
571
      ii = i * nocc +i;
 
572
 
 
573
      if(!local.pairdom_len[ii]) {
 
574
        fprintf(outfile, "\n\tlocal_filter_T1: Pair ii = [%d] is zero-length, which makes no sense.\n",ii);
 
575
        exit(2);
 
576
      }
 
577
 
 
578
      T1tilde = init_array(local.pairdom_len[ii]);
 
579
      T1bar = init_array(local.pairdom_nrlen[ii]);
 
580
 
 
581
      /* Transform the virtuals to the redundant projected virtual basis */
 
582
      C_DGEMV('t', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
583
              &(RIA->matrix[0][i][0]), 1, 0.0, &(T1tilde[0]), 1);
 
584
 
 
585
      /* Transform the virtuals to the non-redundant virtual basis */
 
586
      C_DGEMV('t', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii], 
 
587
              &(T1tilde[0]), 1, 0.0, &(T1bar[0]), 1);
 
588
 
 
589
      for(a=0; a < local.pairdom_nrlen[ii]; a++) {
 
590
        tval = eval + local.eps_occ[i] - local.eps_vir[ii][a];
 
591
        if(fabs(tval) > 0.0001) T1bar[a] /= tval;
 
592
        /* else T1bar[a] = 0.0; */
 
593
      }
 
594
 
 
595
      /* Transform the new T1's to the redundant projected virtual basis */
 
596
      C_DGEMV('n', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii],
 
597
              &(T1bar[0]), 1, 0.0, &(T1tilde[0]), 1);
 
598
 
 
599
      /* Transform the new T1's to the MO basis */
 
600
      C_DGEMV('n', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
601
              &(T1tilde[0]), 1, 0.0, &(RIA->matrix[0][i][0]), 1);
 
602
 
 
603
      free(T1bar);
 
604
      free(T1tilde);
 
605
    }
 
606
 
 
607
    dpd_file2_mat_wrt(RIA);
 
608
    dpd_file2_mat_close(RIA);
 
609
 
 
610
    /* Free Local Memory */
 
611
    for(i=0; i < nocc*nocc; i++) {
 
612
      free_block(local.W[i]);
 
613
      free_block(local.V[i]);
 
614
      free(local.eps_vir[i]);
 
615
    }
 
616
    free(local.W);
 
617
    free(local.V);
 
618
    free(local.eps_vir);
 
619
 
 
620
    free(local.eps_occ);
 
621
    free(local.pairdom_len);
 
622
    free(local.pairdom_nrlen);
 
623
  }
 
624
  else {
 
625
 
 
626
    dpd_file2_mat_init(RIA);
 
627
    dpd_file2_mat_rd(RIA);
 
628
    dpd_file2_init(&DIA, EOM_D, C_irr, 0, 1, "DIA");
 
629
    dpd_file2_mat_init(&DIA);
 
630
    dpd_file2_mat_rd(&DIA);
 
631
    for(h=0; h < nirreps; h++)
 
632
      for(i=0; i < RIA->params->rowtot[h]; i++)
 
633
        for(a=0; a < RIA->params->coltot[h^C_irr]; a++) {
 
634
          tval = eval - DIA.matrix[h][i][a];
 
635
          if (fabs(tval) > 0.0001) RIA->matrix[h][i][a] /= tval;
 
636
        }
 
637
    dpd_file2_mat_wrt(RIA);
 
638
    dpd_file2_mat_close(RIA);
 
639
    dpd_file2_close(&DIA);
 
640
  }
 
641
 
 
642
  return;
 
643
}
 
644
 
 
645
void schmidt_add_SS(dpdfile2 *RIA, dpdfile2 *Ria, int C_irr, int *numCs)
 
646
{
 
647
  double dotval, norm;
 
648
  int i, I;
 
649
  dpdfile2 Cme, CME;
 
650
  char CME_lbl[32], Cme_lbl[32], CMNEF_lbl[32], Cmnef_lbl[32], CMnEf_lbl[32];
 
651
 
 
652
  for (i=0; i<*numCs; i++) {
 
653
    sprintf(CME_lbl, "%s %d", "CME", i);
 
654
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
655
    sprintf(Cme_lbl, "%s %d", "Cme", i);
 
656
    if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
657
    else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
658
    dotval  = dpd_file2_dot(RIA, &CME);
 
659
    dotval += dpd_file2_dot(Ria, &Cme);
 
660
    dpd_file2_axpy(&CME, RIA, -1.0*dotval, 0);
 
661
    dpd_file2_axpy(&Cme, Ria, -1.0*dotval, 0);
 
662
    dpd_file2_close(&CME);
 
663
    dpd_file2_close(&Cme);
 
664
  }
 
665
 
 
666
  norm = norm_C1(RIA, Ria);
 
667
 
 
668
  if (norm < eom_params.schmidt_add_residual_tol) {
 
669
    return;
 
670
  }
 
671
  else {
 
672
    scm_C1(RIA, Ria, 1.0/norm);
 
673
    sprintf(CME_lbl, "%s %d", "CME", *numCs);
 
674
    sprintf(Cme_lbl, "%s %d", "Cme", *numCs);
 
675
    dpd_file2_copy(RIA, EOM_CME, CME_lbl);
 
676
    dpd_file2_copy(Ria, EOM_Cme, Cme_lbl);
 
677
    ++(*numCs);
 
678
  }
 
679
  return;
 
680
}
 
681
 
 
682
void schmidt_add_SS_RHF(dpdfile2 *RIA, int C_irr, int *numCs)
 
683
{
 
684
  double dotval, norm;
 
685
  int i, I;
 
686
  dpdfile2 CME;
 
687
  char CME_lbl[32], Cme_lbl[32];
 
688
 
 
689
  for (i=0; i<*numCs; i++) {
 
690
    sprintf(CME_lbl, "%s %d", "CME", i);
 
691
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
692
    dotval  = 2.0 * dpd_file2_dot(RIA, &CME);
 
693
    dpd_file2_axpy(&CME, RIA, -1.0*dotval, 0);
 
694
    dpd_file2_close(&CME);
 
695
  }
 
696
 
 
697
  norm = norm_C1_rhf(RIA);
 
698
 
 
699
  if (norm < eom_params.schmidt_add_residual_tol) {
 
700
    return;
 
701
  }
 
702
  else {
 
703
    dpd_file2_scm(RIA, 1.0/norm);
 
704
    sprintf(CME_lbl, "%s %d", "CME", *numCs);
 
705
    dpd_file2_copy(RIA, EOM_CME, CME_lbl);
 
706
    ++(*numCs);
 
707
  }
 
708
  return;
 
709
}
 
710
 
 
711
void restart_SS(double **alpha, int L, int num, int C_irr) {
 
712
  int i,j,I;
 
713
  char lbl[20];
 
714
  dpdfile2 C1, CME, Cme, CME2, Cme2;
 
715
  dpdbuf4 C2, CMNEF, Cmnef, CMnEf;
 
716
  double norm, dotval;
 
717
 
 
718
  for (I=1;I<num;++I) {
 
719
    for (i=0; i<I; i++) {
 
720
      dotval = 0.0;
 
721
      for (j=0;j<L;++j) {
 
722
        dotval += alpha[j][i] * alpha[j][I];
 
723
      }
 
724
      for (j=0; j<L; j++) alpha[j][I] -= dotval * alpha[j][i];
 
725
    }
 
726
    dotval = 0.0;
 
727
    for (j=0;j<L;++j) dotval += alpha[j][I] * alpha[j][I];
 
728
    norm = sqrt(dotval);
 
729
    for (j=0;j<L;++j) alpha[j][I] = alpha[j][I]/norm;
 
730
  }
 
731
 
 
732
 
 
733
  for (i=0; i<num; ++i) {
 
734
    sprintf(lbl, "%s %d", "CME", L+i);
 
735
    dpd_file2_init(&C1, EOM_CME, C_irr, 0, 1, lbl);
 
736
    dpd_file2_scm(&C1, 0.0);
 
737
    for (j=0;j<L;++j) {
 
738
      sprintf(lbl, "%s %d", "CME", j);
 
739
      dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
740
      dpd_file2_axpy(&CME, &C1, alpha[j][i], 0);
 
741
      dpd_file2_close(&CME);
 
742
    }
 
743
    dpd_file2_close(&C1);
 
744
 
 
745
    if (params.eom_ref > 0) {
 
746
      sprintf(lbl, "%s %d", "Cme", L+i);
 
747
      if (params.eom_ref == 1) dpd_file2_init(&C1, EOM_Cme, C_irr, 0, 1, lbl);
 
748
      else if (params.eom_ref == 2) dpd_file2_init(&C1, EOM_Cme, C_irr, 2, 3, lbl);
 
749
      dpd_file2_scm(&C1, 0.0);
 
750
      for (j=0;j<L;++j) {
 
751
        sprintf(lbl, "%s %d", "Cme", j);
 
752
        if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
753
        else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
754
        dpd_file2_axpy(&Cme, &C1, alpha[j][i],0);
 
755
        dpd_file2_close(&Cme);
 
756
      }
 
757
      dpd_file2_close(&C1);
 
758
    }
 
759
  }
 
760
 
 
761
  for (i=0; i<num; ++i) {
 
762
    sprintf(lbl, "%s %d", "CME", L+i);
 
763
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, lbl);
 
764
    sprintf(lbl, "%s %d", "CME", i);
 
765
    dpd_file2_copy(&CME, EOM_CME, lbl);
 
766
    dpd_file2_close(&CME);
 
767
    if (params.eom_ref > 0) {
 
768
      sprintf(lbl, "%s %d", "Cme", L+i);
 
769
      if (params.eom_ref == 1) dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, lbl);
 
770
      else if (params.eom_ref == 2) dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, lbl);
 
771
      sprintf(lbl, "%s %d", "Cme", i);
 
772
      dpd_file2_copy(&Cme, EOM_Cme, lbl);
 
773
      dpd_file2_close(&Cme);
 
774
    }
 
775
  }
 
776
  return;
 
777
}
 
778
 
 
779
 
 
780
}} // namespace psi::cceom