← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/detci/mpn.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/detci/mpn.cc
 
1
/*
 
2
** Moeller Plesset Perturbation Series Generator
 
3
**
 
4
** Matt L. Leininger
 
5
** August 11, 1998
 
6
**
 
7
*/
 
8
 
 
9
/* #define DEBUG */
 
10
 
 
11
#include <math.h>
 
12
 
 
13
extern "C" {
 
14
   #include <stdlib.h>
 
15
   #include <stdio.h>
 
16
   #include <libciomr/libciomr.h>
 
17
   #include <libqt/qt.h>
 
18
   #include <libchkpt/chkpt.h>
 
19
   #include "structs.h"
 
20
   #include "ci_tol.h"
 
21
   #define EXTERN
 
22
   #include "globals.h"
 
23
   extern void print_vec(unsigned int nprint, int *Iacode, int *Ibcode,
 
24
      int *Iaidx, int *Ibidx, double *coeff,
 
25
      struct olsen_graph *AlphaG, struct olsen_graph *BetaG,
 
26
      struct stringwr **alplist, struct stringwr **betlist,
 
27
      FILE *outfile);
 
28
}
 
29
 
 
30
#include "civect.h"
 
31
 
 
32
void mpn_generator(CIvect &Hd, struct stringwr **alplist, 
 
33
      struct stringwr **betlist)
 
34
{
 
35
  double *mpk_energy, *mp2k_energy, *oei, *tei, *buffer1, *buffer2;
 
36
  double tval, Empn = 0.0, **wfn_overlap, Empn2 = 0.0, **cvec_coeff;
 
37
  double Empn2a = 0.0;
 
38
  double *cvec_norm, norm, *tmp_coeff, tmp_norm, max_overlap = 1.0;
 
39
  int i, j, k, order, did_vec=0;
 
40
  int kvec_offset; /* offset if c_0 is not stored on disk */
 
41
 
 
42
  CIvect Cvec;
 
43
  CIvect Cvec2;
 
44
  CIvect Sigma;
 
45
 
 
46
  Cvec.set(CIblks.vectlen,CIblks.num_blocks,Parameters.icore,Parameters.Ms0,
 
47
     CIblks.Ia_code, CIblks.Ib_code, CIblks.Ia_size, CIblks.Ib_size,
 
48
     CIblks.offset, CIblks.num_alp_codes, CIblks.num_bet_codes,
 
49
     CalcInfo.nirreps, AlphaG->subgr_per_irrep, Parameters.maxnvect,
 
50
     Parameters.num_c_tmp_units, Parameters.first_c_tmp_unit,
 
51
     CIblks.first_iablk, CIblks.last_iablk, CIblks.decode);
 
52
  Sigma.set(CIblks.vectlen,CIblks.num_blocks,Parameters.icore,Parameters.Ms0,
 
53
     CIblks.Ia_code, CIblks.Ib_code, CIblks.Ia_size, CIblks.Ib_size,
 
54
     CIblks.offset, CIblks.num_alp_codes, CIblks.num_bet_codes,
 
55
     CalcInfo.nirreps, AlphaG->subgr_per_irrep, 1,
 
56
     Parameters.num_s_tmp_units, Parameters.first_s_tmp_unit,
 
57
     CIblks.first_iablk, CIblks.last_iablk, CIblks.decode);
 
58
  Cvec2.set(CIblks.vectlen,CIblks.num_blocks,Parameters.icore,Parameters.Ms0,
 
59
     CIblks.Ia_code, CIblks.Ib_code, CIblks.Ia_size, CIblks.Ib_size,
 
60
     CIblks.offset, CIblks.num_alp_codes, CIblks.num_bet_codes,
 
61
     CalcInfo.nirreps, AlphaG->subgr_per_irrep, Parameters.maxnvect,
 
62
     Parameters.num_c_tmp_units, Parameters.first_c_tmp_unit,
 
63
     CIblks.first_iablk, CIblks.last_iablk, CIblks.decode);
 
64
 
 
65
  /* set up the vector pointers/info */
 
66
  if (Cvec.read_new_first_buf() == -1) Cvec.write_new_first_buf();
 
67
  if (Sigma.read_new_first_buf() == -1) Sigma.write_new_first_buf();
 
68
 
 
69
 
 
70
  Cvec.h0block_buf_init();
 
71
  buffer1 = *(Hd.blockptr(0));
 
72
  buffer2 = Hd.buf_malloc();
 
73
  Hd.buf_unlock(); 
 
74
  
 
75
  wfn_overlap = init_matrix(Parameters.maxnvect+1, Parameters.maxnvect+1);
 
76
  cvec_coeff = init_matrix(Parameters.maxnvect+1, Parameters.maxnvect+1);
 
77
  tmp_coeff = init_array(Parameters.maxnvect+1);
 
78
  cvec_norm = init_array(Parameters.maxnvect+1);
 
79
  mpk_energy = init_array(Parameters.maxnvect+1);
 
80
  mp2k_energy = init_array(2*(Parameters.maxnvect+1)+1);
 
81
  mpk_energy[0] = mp2k_energy[0] = CalcInfo.e0;
 
82
  mpk_energy[1] = CalcInfo.e1 = mp2k_energy[1] =  
 
83
       CalcInfo.escf - CalcInfo.e0 - CalcInfo.enuc;
 
84
  if (CalcInfo.iopen) kvec_offset = 0;
 
85
  else kvec_offset = 0;
 
86
 
 
87
  for (i=0; i<=Parameters.maxnvect; i++) cvec_coeff[i][i] = 1.0; 
 
88
  cvec_norm[0] = 1.0;
 
89
  wfn_overlap[0][0] = 1.0;
 
90
 
 
91
  /* oei = CalcInfo.tf_onel_ints; */
 
92
  if (Parameters.fci) oei = CalcInfo.tf_onel_ints;
 
93
  else oei = CalcInfo.gmat[0];
 
94
  tei = CalcInfo.twoel_ints;
 
95
 
 
96
  fprintf(outfile,"   CalcInfo.escf = %25.15f\n", CalcInfo.escf);
 
97
  fprintf(outfile,"   CalcInfo.e0   = %25.15f\n", CalcInfo.e0);
 
98
  fprintf(outfile,"   CalcInfo.enuc = %25.15f\n", CalcInfo.enuc);
 
99
  fprintf(outfile,"   CalcInfo.e1   = %25.15f\n", CalcInfo.e1);
 
100
  fprintf(outfile,"   n         Corr. Energy \t\t E(MPn) \t\t"
 
101
                  "   n         Corr. Energy \t\t E(MPn)\n\n");
 
102
  fprintf(outfile,"   0  %25.15f %25.15f\n", 0.0000000000,
 
103
                  CalcInfo.e0+CalcInfo.enuc);
 
104
  fprintf(outfile,"   1  %25.15f %25.15f\n",mpk_energy[1]+CalcInfo.enuc, 
 
105
          mpk_energy[0]+mpk_energy[1]+CalcInfo.enuc);
 
106
  Empn = mpk_energy[0]+mpk_energy[1]+CalcInfo.enuc;
 
107
  Empn2 = Empn;
 
108
  fflush(outfile);
 
109
 
 
110
  Cvec.buf_lock(buffer1);
 
111
 
 
112
  Cvec.h0block_buf_init();
 
113
  tval = 1.0;
 
114
  Cvec.init_vals(0, 1, &(CalcInfo.ref_alp_list), &(CalcInfo.ref_alp_rel),
 
115
          &(CalcInfo.ref_bet_list), &(CalcInfo.ref_bet_rel), H0block.blknum, 
 
116
          &tval);
 
117
  if (Parameters.print_lvl > 5) {
 
118
    fprintf(outfile,"Zeroth-order wavefunction\n");
 
119
    Cvec.print(outfile); 
 
120
    }
 
121
 
 
122
  Sigma.buf_lock(buffer2);
 
123
  Cvec.read(0, 0);  /* Set Cvec up correctly ? */
 
124
  /* fprintf(outfile,"Cvec zero_blocks:\n");
 
125
  Cvec.print_zero_blocks(); */
 
126
  Sigma.set_zero_blocks_all();
 
127
  /* fprintf(outfile,"Sigma zero_blocks after set_zero_blocks_all.\n");
 
128
  Sigma.print_zero_blocks(); */
 
129
  sigma(alplist, betlist, Cvec, Sigma, oei, tei, Parameters.fci, 0); 
 
130
  if (Parameters.print_lvl > 5) {
 
131
    fprintf(outfile,"Sigma vector for 0 C vector\n");
 
132
    Sigma.read(0,0);
 
133
    Sigma.print(outfile);
 
134
    }
 
135
  /* fprintf(outfile,"Sigma zero_blocks after sigma call.\n");
 
136
  Sigma.print_zero_blocks(); */
 
137
 
 
138
  Cvec.read(0,0);  /* Set Cvec up correctly ? */
 
139
  Sigma.read(0,0); /* Set Sigma up correctly ? */
 
140
  tval = Cvec * Sigma;
 
141
 
 
142
 /*
 
143
  fprintf(outfile," CalcInfo.enuc = %25.15f\n", CalcInfo.enuc);
 
144
  fprintf(outfile," <psi0|Hc|psi0> = %25.15f\n", tval);
 
145
  fprintf(outfile," CalcInfo.efzc = %25.15f\n", CalcInfo.efzc);
 
146
  fprintf(outfile," mpk_energy[0] = %25.15f\n", mpk_energy[0]);
 
147
 */ 
 
148
  tval += CalcInfo.efzc - mpk_energy[0]; 
 
149
 
 
150
  fprintf(outfile,"   1  %25.15f %25.15f\n", tval+CalcInfo.enuc, 
 
151
   tval+mpk_energy[0]+CalcInfo.enuc);
 
152
  if (tval - mpk_energy[1] > ZERO) 
 
153
    fprintf(outfile, "First-order energies do not agree!\n");
 
154
  fflush(outfile);
 
155
 
 
156
  Cvec.copy_zero_blocks(Sigma); /* Probably don't need this anymore */ 
 
157
  Cvec.copy(Sigma, (1-kvec_offset), 0);
 
158
  if (Parameters.print_lvl > 5) {
 
159
    fprintf(outfile, "Cvec copying Sigma.\n");
 
160
    Cvec.print(outfile);
 
161
    }
 
162
 
 
163
  Sigma.buf_unlock();
 
164
  Hd.buf_lock(buffer2); 
 
165
 
 
166
  tval = Cvec.dcalc2(1, CalcInfo.e0, Hd, 1, alplist, betlist);
 
167
  Hd.buf_unlock();
 
168
 
 
169
  tval = 0.0;
 
170
  if (Parameters.print_lvl > 5) {
 
171
    fprintf(outfile, "Cvec after dcalc2.\n");
 
172
    Cvec.print(outfile);
 
173
    }
 
174
  Cvec.read((1-kvec_offset),0);
 
175
  Cvec.set_vals((1-kvec_offset), 1, &(CalcInfo.ref_alp_list), 
 
176
          &(CalcInfo.ref_alp_rel), &(CalcInfo.ref_bet_list), 
 
177
          &(CalcInfo.ref_bet_rel), H0block.blknum, &tval);
 
178
  Sigma.buf_lock(buffer2);
 
179
 
 
180
  if (Parameters.print_lvl > 5) {
 
181
    fprintf(outfile, "Cvec after set_vals.\n");
 
182
    Cvec.print(outfile);
 
183
    }
 
184
 
 
185
  /* Here buffer1 = Cvec and buffer2 = Sigma */
 
186
  k=1;
 
187
  while (k<Parameters.maxnvect) {
 
188
 
 
189
     /* Form Sigma */
 
190
     Sigma.read(0, 0);
 
191
     Cvec.read((k-kvec_offset), 0);  /* Set Cvec up correctly ? */
 
192
     sigma(alplist, betlist, Cvec, Sigma, oei, tei, Parameters.fci, 0);
 
193
     /* fprintf(outfile,"Sigma zero_blocks after sigma call.\n");
 
194
     Sigma.print_zero_blocks(); */
 
195
 
 
196
     /* Compute k+1, 2k, and 2k+1 th order energies from kth order wavefunction */ 
 
197
     Sigma.read(0,0); /* S_k is located in first Sigma space */
 
198
 
 
199
     #ifdef DEBUG
 
200
       tval = Sigma * Sigma;
 
201
       fprintf(outfile," Sigma * Sigma = %20.10f\n", tval);
 
202
       tval = sqrt(tval);
 
203
       fprintf(outfile," norm of Sigma = %20.10f\n", tval);
 
204
     #endif
 
205
 
 
206
     if (CalcInfo.iopen) { 
 
207
       Cvec.read(0,0); /* E_k = C_0 x S_k */
 
208
       tval = Cvec * Sigma;
 
209
       }
 
210
     else Sigma.extract_vals(0, 1, &(CalcInfo.ref_alp_list),
 
211
          &(CalcInfo.ref_alp_rel), &(CalcInfo.ref_bet_list),
 
212
          &(CalcInfo.ref_bet_rel), H0block.blknum, &tval);
 
213
     mpk_energy[k+1] = tval; 
 
214
     Empn += tval;
 
215
     fprintf(outfile,"  %2d  %25.15f %25.15f", k+1, mpk_energy[k+1], Empn); 
 
216
 
 
217
     Sigma.buf_unlock();
 
218
     if (Parameters.wigner) {
 
219
       if (Parameters.zero_blocks) Cvec2.copy_offset_filenumber(Cvec);
 
220
       Cvec.wigner_E2k_formula(Hd, Sigma, Cvec2, alplist, betlist, buffer1,
 
221
          buffer2, k, mp2k_energy, wfn_overlap, cvec_coeff, cvec_norm, kvec_offset);
 
222
       Empn2 += mp2k_energy[2*k];
 
223
       fprintf(outfile,"\t %2d %25.15f %25.15f\n",2*k,mp2k_energy[2*k],Empn2);
 
224
 
 
225
       /* 25 November 2003 - JMT
 
226
        * Moified to save MP(2n-2) energy */
 
227
       Empn2a = Empn2;
 
228
 
 
229
       Empn2 += mp2k_energy[2*k+1];
 
230
       fprintf(outfile,"\t\t\t\t\t\t\t\t"
 
231
               " %2d %25.15f %25.15f\n", 2*k+1, mp2k_energy[2*k+1], Empn2);
 
232
       }
 
233
     else
 
234
       fprintf(outfile, "\n");
 
235
 
 
236
     fflush(outfile);
 
237
 
 
238
     if (k+1 == Parameters.maxnvect) break;
 
239
 
 
240
     /* Schmidt orthogonalize vector space */
 
241
     if (Parameters.mpn_schmidt && k>1) {
 
242
       Cvec2.buf_lock(buffer2);
 
243
       if (Cvec.schmidt_add2(Cvec2,0,k-2,k-1,k-1,cvec_coeff[k-1],
 
244
           (&cvec_norm[k-1]),&max_overlap)) did_vec = 1;
 
245
       else {
 
246
           fprintf(outfile," %d vector norm = %20.15f < %20.15f\n",
 
247
                   k-1, cvec_norm[k-1], MPn_NORM_TOL);
 
248
           exit(0);
 
249
           } 
 
250
       while (max_overlap > MPn_ZERO) {
 
251
         fprintf(outfile,"Second Schmidt-Orthogonalization performed.\n");
 
252
         Cvec.read(k-1,0);
 
253
         tval = Cvec * Cvec;
 
254
         tval = 1.0/sqrt(tval);
 
255
         fprintf(outfile,"Norm constant for %d S.O. Cvec = %20.15f\n",
 
256
                 k-1, tval);
 
257
         Cvec2.read(k-1,0);
 
258
         tval = Cvec2 * Cvec2;
 
259
         tval = 1.0/sqrt(tval);
 
260
         fprintf(outfile,"Norm constant for %d S.O. Cvec2 = %20.15f\n",
 
261
                 k-1, tval);
 
262
         if (Cvec.schmidt_add2(Cvec2,0,k-2,k-1,k-1,tmp_coeff,
 
263
            &tmp_norm,&max_overlap)) did_vec = 1;
 
264
         else {
 
265
           fprintf(outfile," %d vector norm = %20.15f < %20.15f\n",
 
266
                  k-1, cvec_norm[k-1], MPn_NORM_TOL);
 
267
           exit(0);
 
268
           }
 
269
         for (i=0; i<k-1; i++) { 
 
270
            fprintf(outfile, "second coeff[%d] = %20.15f\n",i,tmp_coeff[i]);
 
271
            cvec_coeff[k-1][i] += (tmp_coeff[i]/cvec_norm[k-1]);
 
272
            }
 
273
         fprintf(outfile, "max_overlap = %20.15f\n", max_overlap);
 
274
         fprintf(outfile, "cvec_norm = %20.15f\n",cvec_norm[k-1]);
 
275
         fprintf(outfile, "tmp_norm = %20.15f\n",tmp_norm);
 
276
         cvec_norm[k-1] *= tmp_norm;
 
277
         fprintf(outfile, "cvec_norm new = %20.15f\n",cvec_norm[k-1]);
 
278
         }
 
279
       Cvec2.buf_unlock();
 
280
       }
 
281
 
 
282
     /* Construct k+1th order wavefunction */
 
283
     if (Parameters.zero_blocks) Cvec2.copy_offset_filenumber(Cvec);
 
284
     Cvec.construct_kth_order_wf(Hd, Sigma, Cvec2, alplist, betlist, buffer1,
 
285
        buffer2, k+1, mpk_energy, cvec_coeff, cvec_norm);
 
286
     if (Parameters.zero_blocks) Cvec2.copy_offset_filenumber(Cvec);
 
287
 
 
288
  /*   fprintf(outfile, "Cvec %d = \n", k+1);
 
289
     Cvec.print(outfile);
 
290
  */
 
291
 
 
292
     if (Parameters.mpn_schmidt) {
 
293
       fprintf(outfile, "cvec_coeff = \n");
 
294
       print_mat(cvec_coeff, k-1, k-1, outfile);
 
295
       }
 
296
 
 
297
     tval = 0.0;
 
298
     Cvec.set_vals(k+1, 1, &(CalcInfo.ref_alp_list), &(CalcInfo.ref_alp_rel),
 
299
             &(CalcInfo.ref_bet_list), &(CalcInfo.ref_bet_rel), H0block.blknum,
 
300
             &tval);
 
301
     Sigma.buf_lock(buffer2);
 
302
     k++;
 
303
     Cvec.copy_zero_blocks(Sigma); 
 
304
     }
 
305
 
 
306
   /* 22 Nov 2003 - JMT 
 
307
    * Save the MPn or MP(2n-1) energy
 
308
    */
 
309
   chkpt_init(PSIO_OPEN_OLD);
 
310
   if (Parameters.save_mpn2 == 1 && Parameters.wigner) {
 
311
     chkpt_wt_etot(Empn2);
 
312
     fprintf(outfile, "\nMP%d energy saved\n", (Parameters.maxnvect * 2) - 1);
 
313
   }
 
314
   else if (Parameters.save_mpn2 == 2 && Parameters.wigner) {
 
315
     chkpt_wt_etot(Empn2a);
 
316
     fprintf(outfile, "\nMP%d energy saved\n", (Parameters.maxnvect * 2) - 2);
 
317
   }
 
318
   else {
 
319
     chkpt_wt_etot(Empn);
 
320
     fprintf(outfile, "\nMP%d energy saved\n", Parameters.maxnvect);
 
321
   }
 
322
   chkpt_close();
 
323
 
 
324
   fprintf(outfile,"\n");
 
325
}
 
326