← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/detci/calc_hd_block.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/detci/calc_hd_block.cc
 
1
/*! \file
 
2
    \ingroup DETCI
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
#define EXTERN
 
7
#include <cstdio>
 
8
#include <cstdlib>
 
9
#include <libciomr/libciomr.h>
 
10
#include <libqt/qt.h>
 
11
#include "structs.h"
 
12
#include "globals.h"
 
13
 
 
14
namespace psi { namespace detci {
 
15
 
 
16
extern int calc_orb_diff(int cnt, unsigned char *I, unsigned char *J,
 
17
   int *I_alpha_diff, int *J_alpha_diff, int *sign, int *same, int extended);
 
18
extern int *ioff;
 
19
 
 
20
/* C "GLOBAL" VARIABLES FOR THIS MODULE */
 
21
extern struct stringwr **alplist;
 
22
extern struct stringwr **betlist;
 
23
 
 
24
 
 
25
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
 
26
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
 
27
 
 
28
 
 
29
/*
 
30
** calc_hd_block(): Function calculates a block of H0, the diagonal elements of
 
31
**    the Hamiltonian matrix.
 
32
**
 
33
** Parameters:
 
34
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
35
**    betlist_local = list of beta strings with replacements
 
36
**    nas     = number of alpha strings in list
 
37
**    nbs     = number of beta strings in list
 
38
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
39
**    oei     = one-electron integrals
 
40
**    tei     = two-electron integrals
 
41
**    na      = number of explicit alpha electrons 
 
42
**    nb      = number of explicit beta electrons
 
43
**    nbf     = number of orbitals in CI
 
44
**    efzc    = frozen core energy
 
45
**
 
46
*/ 
 
47
void calc_hd_block(struct stringwr *alplist_local, struct stringwr *betlist_local,
 
48
      double **H0, double *oei, double *tei, double efzc,
 
49
      int nas, int nbs, int na, int nb, int nbf)
 
50
{
 
51
   int acnt, bcnt;
 
52
   int a1, a2, b1, b2;
 
53
   int i,j, ii, iii, jj, ij, iijj, ijij;
 
54
   double value;
 
55
   struct stringwr *betlist0;
 
56
 
 
57
   betlist0 = betlist_local;
 
58
 
 
59
   for (acnt=0; acnt<nas; acnt++) {
 
60
      
 
61
      for (bcnt=0, betlist_local=betlist0; bcnt<nbs; bcnt++) {
 
62
 
 
63
         /* add frozen core energy first */
 
64
/************************************************/
 
65
         value = efzc; 
 
66
 
 
67
         /* loop over alpha occs */
 
68
         for (a1=0; a1<na; a1++) {
 
69
            i = (int) alplist_local->occs[a1];
 
70
            ii = ioff[i] + i;
 
71
            value += oei[ii];  
 
72
            /* fprintf(outfile,"oei[%d] = %lf\n",ii,oei[ii]); */ 
 
73
            iii = ioff[ii];
 
74
 
 
75
            for (a2=0; a2<a1; a2++) {
 
76
               j = (int) alplist_local->occs[a2];
 
77
               jj = ioff[j] + j;
 
78
               iijj = iii + jj;
 
79
               ij = ioff[i] + j;
 
80
               ijij = ioff[ij] + ij;
 
81
               value += tei[iijj] - tei[ijij]; 
 
82
               }
 
83
 
 
84
            for (b1=0; b1<nb; b1++) {
 
85
               j = (int) betlist_local->occs[b1];
 
86
               jj = ioff[j] + j;
 
87
               iijj = ioff[MAX0(ii,jj)] + MIN0(ii,jj);
 
88
               value += tei[iijj]; 
 
89
               }
 
90
           }
 
91
 
 
92
         for (b1=0; b1<nb; b1++) {
 
93
            i = (int) betlist_local->occs[b1];
 
94
            ii = ioff[i] + i;
 
95
            value += oei[ii]; 
 
96
            iii = ioff[ii];
 
97
 
 
98
            for (b2=0; b2<b1; b2++) {
 
99
               j = (int) betlist_local->occs[b2];
 
100
               jj = ioff[j] + j;
 
101
               iijj = iii + jj;
 
102
               ij = ioff[i] + j;
 
103
               ijij = ioff[ij] + ij;
 
104
               value += tei[iijj] - tei[ijij]; 
 
105
               }
 
106
            } 
 
107
 
 
108
         H0[acnt][bcnt] = value;
 
109
      /*   
 
110
         fprintf(outfile,"H0[%d][%d] = %lf\n",acnt,bcnt,value); 
 
111
      */ 
 
112
         betlist_local++;
 
113
         } /* end loop over bcnt */
 
114
 
 
115
      alplist_local++;
 
116
      }
 
117
 
 
118
}
 
119
/*
 
120
** calc_hd_block_ave(): Function calculates a block of H0 and the diagonal elements
 
121
** of the Hamiltonian matrix averaged over spin-coupling sets to correct any
 
122
** spin contamination of the c and sigma vectors. 
 
123
**
 
124
** Parameters:
 
125
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
126
**    betlist_local = list of beta strings with replacements
 
127
**    nas     = number of alpha strings in list
 
128
**    nbs     = number of beta strings in list
 
129
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
130
**    oei     = one-electron integrals
 
131
**    tei     = two-electron integrals
 
132
**    na      = number of explicit alpha electrons 
 
133
**    nb      = number of explicit beta electrons
 
134
**    nbf     = number of orbitals in CI
 
135
**    efzc    = frozen core energy
 
136
**
 
137
*/ 
 
138
void calc_hd_block_ave(struct stringwr *alplist_local, struct stringwr *betlist_local,
 
139
      double **H0, double *tf_oei, double *tei, double efzc,
 
140
      int nas, int nbs, int na, int nb, int nbf)
 
141
{
 
142
   int acnt, bcnt;
 
143
   int a1, a2, a3, b1, b2, b3;
 
144
   int i,j, ii, iii, jj, ij, iijj, ijij;
 
145
   double value, tval, tval2, Kave;
 
146
   struct stringwr *betlist0;
 
147
   double k_total; /* total number of K ints in energy expression */
 
148
   int k_combo; /* total combination of unique K ints over spin-coupling set */
 
149
   int *unique_occs; /* the uniquely occupied orbitals for a given determinant */
 
150
   int num_el;  /* total number of electrons explicitly treated */
 
151
   int num_unique; /* number of unique orbitals */
 
152
   betlist0 = betlist_local;
 
153
   
 
154
 
 
155
   k_total = combinations(na,2) + combinations(nb,2); 
 
156
 
 
157
   num_el = na + nb;
 
158
   unique_occs = init_int_array(num_el);
 
159
 
 
160
   for (acnt=0; acnt<nas; acnt++) {
 
161
      
 
162
      for (bcnt=0, betlist_local=betlist0; bcnt<nbs; bcnt++) {
 
163
 
 
164
         /* add frozen core energy first */
 
165
         value = efzc; 
 
166
 
 
167
         /* loop over alpha occs */
 
168
         for (a1=0; a1<na; a1++) {
 
169
            i = (int) alplist_local->occs[a1];
 
170
            ii = ioff[i] + i;
 
171
            /* h_ii bar alpha alpha */
 
172
            value += tf_oei[ii];
 
173
            /* fprintf(outfile,"tf_oei[%d] = %lf\n",ii,tf_oei[ii]); */
 
174
            iii = ioff[ii];
 
175
 
 
176
            /* loop over alpha occs */
 
177
            for (a2=0; a2<a1; a2++) {
 
178
               j = (int) alplist_local->occs[a2];
 
179
               jj = ioff[j] + j;
 
180
               iijj = iii + jj;
 
181
               /* J alpha alpha */
 
182
               value += tei[iijj]; 
 
183
               }
 
184
 
 
185
            /* loop over beta occs */
 
186
            for (b1=0; b1<nb; b1++) {
 
187
               j = (int) betlist_local->occs[b1];
 
188
               jj = ioff[j] + j;
 
189
               iijj = ioff[MAX0(ii,jj)] + MIN0(ii,jj);
 
190
               value += tei[iijj]; 
 
191
               }
 
192
           }
 
193
 
 
194
         /* loop over beta occs */
 
195
         for (b1=0; b1<nb; b1++) {
 
196
            i = (int) betlist_local->occs[b1];
 
197
            ii = ioff[i] + i;
 
198
            value += tf_oei[ii];
 
199
            /* fprintf(outfile,"tf_oei[%d] = %lf\n",ii,tf_oei[ii]); */
 
200
            iii = ioff[ii];
 
201
 
 
202
            /* loop over beta occs */
 
203
            for (b2=0; b2<b1; b2++) {
 
204
               j = (int) betlist_local->occs[b2];
 
205
               jj = ioff[j] + j;
 
206
               iijj = iii + jj;
 
207
               ij = ioff[i] + j; 
 
208
               ijij = ioff[ij] + ij;
 
209
               value += tei[iijj];  
 
210
               }
 
211
            }
 
212
 
 
213
         /* determine average K over spin-coupling set */
 
214
         num_unique = 0;
 
215
         for (a1=0; a1<na; a1++) unique_occs[num_unique++] = (int) alplist_local->occs[a1];
 
216
         /* for (j=0; j<num_unique; j++) 
 
217
            fprintf(outfile,"unique_occs[%d] = %d\n",j,unique_occs[j]); */
 
218
            for (b1=0; b1<nb; b1++) {
 
219
               j = (int) betlist_local->occs[b1];
 
220
               for (a1=0; a1<na; a1++) {
 
221
                  if (j==unique_occs[a1]) break;
 
222
                  if (a1==(na-1)) unique_occs[num_unique++] = j;
 
223
                  }
 
224
               }
 
225
         /* fprintf(outfile,"num_unique = %d\n",num_unique);
 
226
         fprintf(outfile,"num_el = %d\n",num_el);
 
227
         */
 
228
         if (num_unique>num_el) fprintf(outfile,"WARNING: The number of explicit electrons" \
 
229
                             "!= num_el\n");
 
230
               
 
231
       /*   
 
232
         for (j=0; j<na; j++) 
 
233
            fprintf(outfile,"alp_occs[%d] = %d\n",j,(int)alplist_local->occs[j]);
 
234
         for (j=0; j<nb; j++) 
 
235
            fprintf(outfile,"bet_occs[%d] = %d\n",j,(int)betlist_local->occs[j]);
 
236
         for (j=0; j<num_unique; j++) 
 
237
            fprintf(outfile,"unique_occs[%d] = %d\n",j,unique_occs[j]);
 
238
       */
 
239
 
 
240
         Kave = 0.0;
 
241
         for (a1=0; a1<num_unique; a1++) {
 
242
            i = unique_occs[a1];
 
243
            for (b1=0; b1<a1; b1++) {
 
244
               j = unique_occs[b1];
 
245
               ij = ioff[MAX0(i,j)] + MIN0(i,j);
 
246
               ijij = ioff[ij] + ij;
 
247
               Kave += tei[ijij];
 
248
               /* fprintf(outfile,"tei[%d] = %lf\n",ijij,tei[ijij]); */ 
 
249
               }
 
250
            }
 
251
         
 
252
         /* fprintf(outfile,"num_unique = %d\n",num_unique);
 
253
         fprintf(outfile,"ioff[num_unique-1] = %d\n",ioff[num_unique]);
 
254
         fprintf(outfile,"k_total = %d\n",k_total);
 
255
         */
 
256
 
 
257
         if (num_unique > 1) Kave /= ioff[num_unique-1];
 
258
         value -= 0.5 * Kave * k_total; 
 
259
         /* fprintf(outfile,"Kave = %lf\n",Kave); */
 
260
 
 
261
         if (Parameters.print_lvl > 5) {
 
262
           fprintf(outfile,"acnt = %d\t bcnt = %d\n",acnt,bcnt); 
 
263
           fprintf(outfile,"tval = %lf\n",tval);
 
264
           for(a1=0; a1<na; a1++)
 
265
             fprintf(outfile," %d",alplist_local->occs[a1]);
 
266
           fprintf(outfile," \n");
 
267
           for(b1=0; b1<nb; b1++)
 
268
             fprintf(outfile," %d",betlist_local->occs[b1]);
 
269
           fprintf(outfile," \n");
 
270
           } 
 
271
 
 
272
         H0[acnt][bcnt] = value;
 
273
         /* fprintf(outfile,"H0[%d][%d] = %lf\n",acnt,bcnt,value); */
 
274
         betlist_local++;
 
275
         } /* end loop over bcnt */
 
276
 
 
277
      alplist_local++;
 
278
      }
 
279
 
 
280
}
 
281
 
 
282
/*
 
283
** calc_hd_block_orbenergy(): Function calculates a block of H0 and the diagonal elements
 
284
** of the Hamiltonian matrix as the sum of orbital energies. 
 
285
**
 
286
** Parameters:
 
287
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
288
**    betlist_local = list of beta strings with replacements
 
289
**    nas     = number of alpha strings in list
 
290
**    nbs     = number of beta strings in list
 
291
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
292
**    oei     = one-electron integrals
 
293
**    tei     = two-electron integrals
 
294
**    na      = number of explicit alpha electrons 
 
295
**    nb      = number of explicit beta electrons
 
296
**    nbf     = number of orbitals in CI
 
297
**    efzc    = frozen core energy
 
298
**
 
299
*/ 
 
300
void calc_hd_block_orbenergy(struct stringwr *alplist_local, 
 
301
      struct stringwr *betlist_local, double **H0, double *oei, 
 
302
      double *tei, double efzc, int nas, int nbs, int na, int nb, int nbf)
 
303
{
 
304
   int acnt, bcnt;
 
305
   int a1, b1, i,j; 
 
306
   double value, tval;
 
307
   struct stringwr *betlist0, *alplist0;
 
308
   double *orb_e_diff_alp, *orb_e_diff_bet;
 
309
   double sum_orb_energies = 0.0;
 
310
 
 
311
   betlist0 = betlist_local;
 
312
   alplist0 = alplist_local; 
 
313
 
 
314
   orb_e_diff_alp = init_array(nas);
 
315
   orb_e_diff_bet = init_array(nbs);
 
316
  /* if (Parameters.Ms0) orb_e_diff_bet = &orb_e_diff_alp;
 
317
   else orb_e_diff_bet = init_array(CalcInfo.num_bet_str);
 
318
  */
 
319
 
 
320
   for (acnt=0; acnt<nas; acnt++) {
 
321
      orb_e_diff_alp[acnt] = 0.0;
 
322
      for (a1=0; a1<na; a1++) {
 
323
         i = (int) alplist_local->occs[a1];
 
324
         i += CalcInfo.num_fzc_orbs;
 
325
         if(Parameters.zaptn) 
 
326
           orb_e_diff_alp[acnt] += CalcInfo.scfeigvala[i];
 
327
         else
 
328
           orb_e_diff_alp[acnt] += CalcInfo.scfeigval[i];
 
329
         }
 
330
      alplist_local++;
 
331
      }
 
332
 
 
333
   for (bcnt=0; bcnt<nbs; bcnt++) {
 
334
      orb_e_diff_bet[bcnt] = 0.0;
 
335
      for (b1=0; b1<nb; b1++) {
 
336
         j = (int) betlist_local->occs[b1];
 
337
         j += CalcInfo.num_fzc_orbs;
 
338
         if(Parameters.zaptn) 
 
339
           orb_e_diff_bet[bcnt] += CalcInfo.scfeigvalb[j];
 
340
         else
 
341
           orb_e_diff_bet[bcnt] += CalcInfo.scfeigval[j];
 
342
         }
 
343
      betlist_local++;
 
344
      }
 
345
 
 
346
   alplist_local = alplist0;
 
347
   betlist_local = betlist0;
 
348
 
 
349
   for (acnt=0; acnt<nas; acnt++) {
 
350
         tval = efzc + orb_e_diff_alp[acnt]; 
 
351
      for (bcnt=0; bcnt<nbs; bcnt++) {
 
352
         value = orb_e_diff_bet[bcnt] + tval; 
 
353
         H0[acnt][bcnt] = value;
 
354
        /* 
 
355
         fprintf(outfile,"H0[%d][%d] = %lf\n",acnt,bcnt,value); 
 
356
        */ 
 
357
         betlist_local++;
 
358
         } /* end loop over bcnt */
 
359
      alplist_local++;
 
360
      }
 
361
 
 
362
/* Free up memory */
 
363
free(orb_e_diff_alp);
 
364
free(orb_e_diff_bet);
 
365
 
 
366
 
 
367
}
 
368
 
 
369
/*
 
370
** calc_hd_block_evangelisti(): Function calculates a block of H0 and the diagonal elements
 
371
** of the Hamiltonian matrix averaged over spin-coupling sets to correct any
 
372
** spin contamination of the c and sigma vectors. 
 
373
**
 
374
** Parameters:
 
375
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
376
**    betlist_local = list of beta strings with replacements
 
377
**    nas     = number of alpha strings in list
 
378
**    nbs     = number of beta strings in list
 
379
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
380
**    oei     = one-electron integrals
 
381
**    tei     = two-electron integrals
 
382
**    na      = number of explicit alpha electrons 
 
383
**    nb      = number of explicit beta electrons
 
384
**    nbf     = number of orbitals in CI
 
385
**    efzc    = frozen core energy
 
386
**
 
387
*/ 
 
388
void calc_hd_block_evangelisti(struct stringwr *alplist_local, struct stringwr *betlist_local,
 
389
      double **H0, double *tf_oei, double *tei, double efzc,
 
390
      int nas, int nbs, int na, int nb, int nbf)
 
391
{
 
392
   int acnt, bcnt;
 
393
   int a1, b1, i,j; 
 
394
   double value, tval;
 
395
   struct stringwr *betlist0, *alplist0;
 
396
   double *orb_e_diff_alp, *orb_e_diff_bet;
 
397
   int num_alp_diff, num_bet_diff;
 
398
   int **orb_diff, *jnk;
 
399
   int sign;
 
400
 
 
401
   betlist0 = betlist_local;
 
402
   alplist0 = alplist_local; 
 
403
 
 
404
   orb_diff = init_int_matrix(2,na);
 
405
   jnk = init_int_array(na);
 
406
   orb_e_diff_alp = init_array(nas);
 
407
   orb_e_diff_bet = init_array(nbs);
 
408
 
 
409
   for (acnt=0; acnt<nas; acnt++) {
 
410
      orb_e_diff_alp[acnt] = 0.0;
 
411
      num_alp_diff = calc_orb_diff(na,
 
412
                     alplist[CalcInfo.ref_alp_list][CalcInfo.ref_alp_rel].occs, 
 
413
                     alplist_local->occs, orb_diff[0], orb_diff[1], &sign,
 
414
                     jnk, 1);
 
415
      for (a1=0; a1<num_alp_diff; a1++) {
 
416
         i = orb_diff[0][a1]; 
 
417
         j = orb_diff[1][a1]; 
 
418
         i += CalcInfo.num_fzc_orbs;
 
419
         j += CalcInfo.num_fzc_orbs;
 
420
         orb_e_diff_alp[acnt] += CalcInfo.scfeigval[j] 
 
421
                                 - CalcInfo.scfeigval[i]; 
 
422
         }
 
423
      alplist_local++;
 
424
      }
 
425
 
 
426
   for (bcnt=0; bcnt<nbs; bcnt++) {
 
427
      orb_e_diff_bet[bcnt] = 0.0;
 
428
      num_bet_diff = calc_orb_diff(nb, 
 
429
                     betlist[CalcInfo.ref_bet_list][CalcInfo.ref_bet_rel].occs,
 
430
                     betlist_local->occs, orb_diff[0], orb_diff[1], &sign, 
 
431
                     jnk, 1);
 
432
      for (b1=0; b1<num_bet_diff; b1++) {
 
433
         i = orb_diff[0][b1];
 
434
         j = orb_diff[1][b1];  
 
435
         i += CalcInfo.num_fzc_orbs;
 
436
         j += CalcInfo.num_fzc_orbs;
 
437
         orb_e_diff_bet[bcnt] += CalcInfo.scfeigval[j]
 
438
                                 - CalcInfo.scfeigval[i];
 
439
         }
 
440
      betlist_local++;
 
441
      } 
 
442
 
 
443
   alplist_local = alplist0;
 
444
   betlist_local = betlist0;
 
445
 
 
446
   for (acnt=0; acnt<nas; acnt++) {
 
447
         /* add frozen core energy first */
 
448
         tval = CalcInfo.escf - CalcInfo.enuc; 
 
449
         tval += orb_e_diff_alp[acnt]; 
 
450
      for (bcnt=0; bcnt<nbs; bcnt++) {
 
451
         value = 0.0;
 
452
         value = orb_e_diff_bet[bcnt] + tval; 
 
453
         H0[acnt][bcnt] = value;
 
454
         /* fprintf(outfile,"H0[%d][%d] = %lf\n",acnt,bcnt,value); */ 
 
455
         betlist_local++;
 
456
         } /* end loop over bcnt */
 
457
      alplist_local++;
 
458
      }
 
459
 
 
460
/* Free memory */
 
461
/*
 
462
free(jnk);
 
463
free(orb_e_diff_alp);
 
464
free(orb_e_diff_bet);
 
465
free(orb_diff);
 
466
*/
 
467
 
 
468
 
 
469
}
 
470
 
 
471
 
 
472
/*
 
473
** calc_hd_block_mll(): Function calculates a block of H0 and the diagonal elements
 
474
** of the Hamiltonian matrix as the sum of orbital energies. 
 
475
**
 
476
** Parameters:
 
477
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
478
**    betlist_local = list of beta strings with replacements
 
479
**    nas     = number of alpha strings in list
 
480
**    nbs     = number of beta strings in list
 
481
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
482
**    oei     = one-electron integrals
 
483
**    tei     = two-electron integrals
 
484
**    na      = number of explicit alpha electrons 
 
485
**    nb      = number of explicit beta electrons
 
486
**    nbf     = number of orbitals in CI
 
487
**    efzc    = frozen core energy
 
488
**
 
489
*/ 
 
490
void calc_hd_block_mll(struct stringwr *alplist_local, 
 
491
      struct stringwr *betlist_local, double **H0, double *oei, 
 
492
      double *tei, double efzc, int nas, int nbs, int na, int nb, int nbf)
 
493
{
 
494
   int acnt, bcnt;
 
495
   int a1, b1, i,j, i_offset, j_offset, ii, jj; 
 
496
   double value, tval;
 
497
   struct stringwr *betlist0, *alplist0;
 
498
   double *orb_e_diff_alp, *orb_e_diff_bet;
 
499
   double *oei_alp, *oei_bet, *eigval;
 
500
 
 
501
   betlist0 = betlist_local;
 
502
   alplist0 = alplist_local; 
 
503
 
 
504
   oei_alp = init_array(nas);
 
505
   oei_bet = init_array(nbs);
 
506
   orb_e_diff_alp = init_array(nas);
 
507
   orb_e_diff_bet = init_array(nbs);
 
508
  /* if (Parameters.Ms0) orb_e_diff_bet = &orb_e_diff_alp;
 
509
   else orb_e_diff_bet = init_array(nbs);
 
510
  */
 
511
 
 
512
   for (acnt=0; acnt<nas; acnt++) {
 
513
      orb_e_diff_alp[acnt] = oei_alp[acnt] = 0.0;
 
514
      for (a1=0; a1<na; a1++) {
 
515
         i = (int) alplist_local->occs[a1];
 
516
         ii = ioff[i] + i;
 
517
         i_offset = i + CalcInfo.num_fzc_orbs;
 
518
         oei_alp[acnt] += oei[ii]; 
 
519
         orb_e_diff_alp[acnt] += CalcInfo.scfeigval[i_offset] - oei[ii];
 
520
         }
 
521
      alplist_local++;
 
522
      }
 
523
 
 
524
   for (bcnt=0; bcnt<nbs; bcnt++) {
 
525
      orb_e_diff_bet[bcnt] = oei_bet[bcnt] = 0.0;
 
526
      for (b1=0; b1<nb; b1++) {
 
527
         j = (int) betlist_local->occs[b1];
 
528
         jj = ioff[j] + j;
 
529
         j_offset = j + CalcInfo.num_fzc_orbs;
 
530
         oei_bet[bcnt] += oei[jj];
 
531
         orb_e_diff_bet[bcnt] += CalcInfo.scfeigval[j_offset] - oei[jj];
 
532
         }
 
533
      betlist_local++;
 
534
      }
 
535
 
 
536
   alplist_local = alplist0;
 
537
   betlist_local = betlist0;
 
538
 
 
539
   for (acnt=0; acnt<nas; acnt++) {
 
540
         tval = efzc + 0.5 * orb_e_diff_alp[acnt] + oei_alp[acnt]; 
 
541
      for (bcnt=0; bcnt<nbs; bcnt++) {
 
542
         value = 0.5 * orb_e_diff_bet[bcnt] + oei_bet[bcnt] + tval; 
 
543
         H0[acnt][bcnt] = value;
 
544
         betlist_local++;
 
545
         } /* end loop over bcnt */
 
546
      alplist_local++;
 
547
      }
 
548
 
 
549
 free(oei_alp);
 
550
 free(oei_bet);
 
551
 free(orb_e_diff_alp);
 
552
 free(orb_e_diff_bet);
 
553
}
 
554
/*
 
555
** calc_hd_block_z_ave(): Function calculates a block of H0 and the diagonal elements
 
556
** of the Hamiltonian matrix averaged over spin-coupling sets to correct any
 
557
** spin contamination of the c and sigma vectors.
 
558
**
 
559
** Parameters:
 
560
**    alplist_local = list of alpha strings with replacements (used to get occs)
 
561
**    betlist_local = list of beta strings with replacements
 
562
**    nas     = number of alpha strings in list
 
563
**    nbs     = number of beta strings in list
 
564
**    H0      = matrix to hold results (stored as H0[alpidx][betidx])
 
565
**    oei     = one-electron integrals
 
566
**    tei     = two-electron integrals
 
567
**    na      = number of explicit alpha electrons
 
568
**    nb      = number of explicit beta electrons
 
569
**    nbf     = number of orbitals in CI
 
570
**    efzc    = frozen core energy
 
571
**
 
572
*/
 
573
void calc_hd_block_z_ave(struct stringwr *alplist_local, struct stringwr *betlist_local, 
 
574
double **H0, double pert_param, double *tei, double efzc, int nas, int nbs, int na, 
 
575
int nb, int nbf)
 
576
{
 
577
   int acnt, bcnt;
 
578
   int a1, a2, a3, b1, b2, b3;
 
579
   int i,j, ii, iii, jj, ij, iijj, ijij;
 
580
   double value, tval, tval2, Kave;
 
581
   struct stringwr *betlist0;
 
582
   double k_total; /* total number of K ints in energy expression */
 
583
   int k_combo; /* total combination of unique K ints over spin-coupling set */
 
584
   int *unique_occs; /* the uniquely occupied orbitals for a given determinant */
 
585
   int num_el;  /* total number of electrons explicitly treated */
 
586
   int num_unique; /* number of unique orbitals */
 
587
   betlist0 = betlist_local;
 
588
 
 
589
 
 
590
   k_total = combinations(na,2) + combinations(nb,2);
 
591
   num_el = na + nb;
 
592
   unique_occs = init_int_array(num_el);
 
593
 
 
594
   for (acnt=0; acnt<nas; acnt++) {
 
595
 
 
596
      for (bcnt=0, betlist_local=betlist0; bcnt<nbs; bcnt++) {
 
597
 
 
598
         /* add frozen core energy first */
 
599
         value = efzc;
 
600
 
 
601
         /* loop over alpha occs */
 
602
         for (a1=0; a1<na; a1++) {
 
603
            i = (int) alplist_local->occs[a1];
 
604
            value += CalcInfo.scfeigval[i+CalcInfo.num_fzc_orbs];
 
605
            ii = ioff[i] + i;
 
606
            /* h_ii bar alpha alpha */
 
607
            iii = ioff[ii];
 
608
 
 
609
            /* loop over alpha occs */
 
610
            for (a2=0; a2<a1; a2++) {
 
611
               j = (int) alplist_local->occs[a2];
 
612
               jj = ioff[j] + j;
 
613
               iijj = iii + jj;
 
614
               /* J alpha alpha */
 
615
               value -= pert_param * tei[iijj];
 
616
               }
 
617
 
 
618
            /* loop over beta occs */
 
619
            for (b1=0; b1<nb; b1++) {
 
620
               j = (int) betlist_local->occs[b1];
 
621
               jj = ioff[j] + j;
 
622
               iijj = ioff[MAX0(ii,jj)] + MIN0(ii,jj);
 
623
               value -= pert_param * tei[iijj];
 
624
               }
 
625
           }
 
626
 
 
627
         /* loop over beta occs */
 
628
         for (b1=0; b1<nb; b1++) {
 
629
            i = (int) betlist_local->occs[b1];
 
630
            value += CalcInfo.scfeigval[i+CalcInfo.num_fzc_orbs];
 
631
            ii = ioff[i] + i;
 
632
            iii = ioff[ii];
 
633
 
 
634
            /* loop over beta occs */
 
635
            for (b2=0; b2<b1; b2++) {
 
636
               j = (int) betlist_local->occs[b2];
 
637
               jj = ioff[j] + j;
 
638
               iijj = iii + jj;
 
639
               ij = ioff[i] + j;
 
640
               ijij = ioff[ij] + ij;
 
641
               value -= pert_param * tei[iijj];
 
642
               }
 
643
            }
 
644
 
 
645
         /* determine average K over spin-coupling set */
 
646
         num_unique = 0;
 
647
         for (a1=0; a1<na; a1++) unique_occs[num_unique++] = (int) alplist_local->occs[a1];
 
648
         /* for (j=0; j<num_unique; j++)
 
649
            fprintf(outfile,"unique_occs[%d] = %d\n",j,unique_occs[j]); */
 
650
            for (b1=0; b1<nb; b1++) {
 
651
               j = (int) betlist_local->occs[b1];
 
652
               for (a1=0; a1<na; a1++) {
 
653
                  if (j==unique_occs[a1]) break;
 
654
                  if (a1==(na-1)) unique_occs[num_unique++] = j;
 
655
                  }
 
656
               }
 
657
         /* fprintf(outfile,"num_unique = %d\n",num_unique);
 
658
         fprintf(outfile,"num_el = %d\n",num_el);
 
659
         */
 
660
         if (num_unique>num_el) fprintf(outfile,"WARNING: The number of explicit electrons" \
 
661
                             "!= num_el\n");
 
662
 
 
663
       /*
 
664
         for (j=0; j<na; j++)
 
665
            fprintf(outfile,"alp_occs[%d] = %d\n",j,(int)alplist_local->occs[j]);
 
666
         for (j=0; j<nb; j++)
 
667
            fprintf(outfile,"bet_occs[%d] = %d\n",j,(int)betlist_local->occs[j]);
 
668
         for (j=0; j<num_unique; j++)
 
669
            fprintf(outfile,"unique_occs[%d] = %d\n",j,unique_occs[j]);
 
670
       */
 
671
 
 
672
         Kave = 0.0;
 
673
         for (a1=0; a1<num_unique; a1++) {
 
674
            i = unique_occs[a1];
 
675
            for (b1=0; b1<a1; b1++) {
 
676
               j = unique_occs[b1];
 
677
               ij = ioff[MAX0(i,j)] + MIN0(i,j);
 
678
               ijij = ioff[ij] + ij;
 
679
               Kave += tei[ijij];
 
680
               /* fprintf(outfile,"tei[%d] = %lf\n",ijij,tei[ijij]); */
 
681
               }
 
682
            }
 
683
 
 
684
         /* fprintf(outfile,"num_unique = %d\n",num_unique);
 
685
         fprintf(outfile,"ioff[num_unique-1] = %d\n",ioff[num_unique]);
 
686
         fprintf(outfile,"k_total = %d\n",k_total);
 
687
         */
 
688
         if (num_unique > 1) Kave /= ioff[num_unique-1];
 
689
         value += 0.5 * Kave * k_total * pert_param;
 
690
         /* fprintf(outfile,"Kave = %lf\n",Kave); */
 
691
 
 
692
         if (Parameters.print_lvl > 5) {
 
693
           fprintf(outfile,"acnt = %d\t bcnt = %d\n",acnt,bcnt);
 
694
           fprintf(outfile,"tval = %lf\n",tval);
 
695
           for(a1=0; a1<na; a1++)
 
696
             fprintf(outfile," %d",alplist_local->occs[a1]);
 
697
           fprintf(outfile," \n");
 
698
           for(b1=0; b1<nb; b1++)
 
699
             fprintf(outfile," %d",betlist_local->occs[b1]);
 
700
           fprintf(outfile," \n");
 
701
           }
 
702
         H0[acnt][bcnt] = value;
 
703
       /*
 
704
         fprintf(outfile,"H0[%d][%d] = %lf\n",acnt,bcnt,value);
 
705
       */
 
706
         betlist_local++;
 
707
         } /* end loop over bcnt */
 
708
 
 
709
      alplist_local++;
 
710
      }
 
711
 
 
712
}
 
713
 
 
714
}} // namespace psi::detci
 
715