← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/mvo/mvo.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/mvo/mvo.cc
 
1
/*! \defgroup MVO mvo: Compute modified virtual orbitals, etc */
 
2
 
 
3
/*! 
 
4
** \file
 
5
** \ingroup MVO
 
6
** \brief Compute modified virtual orbitals, etc.
 
7
**
 
8
** This program reads the frozen core operator from TRANSQT 
 
9
** output ** and diagonalizes the virtual-virtual block to make 
 
10
** Modified Virtual Orbitals (MVO's) as described in 
 
11
** C. W. Bauschlicher, J. Chem. Phys. 72, 880 (1980)
 
12
**
 
13
** C. David Sherrill
 
14
** Georgia Institute of Technology
 
15
** March 9, 2001
 
16
*/
 
17
 
 
18
#include <cstdio>
 
19
#include <cstdlib>
 
20
#include <cmath>
 
21
#include <cstring>
 
22
#include <libipv1/ip_lib.h>
 
23
#include <libpsio/psio.h>
 
24
#include <libciomr/libciomr.h>
 
25
#include <libchkpt/chkpt.h>
 
26
#include <libqt/qt.h>
 
27
#include <psifiles.h>
 
28
#include <libiwl/iwl.h>
 
29
#include "MOInfo.h"
 
30
#include "params.h"
 
31
#include "globals.h"
 
32
 
 
33
 
 
34
/* First definitions of globals */
 
35
extern "C" {
 
36
  FILE *infile, *outfile;
 
37
  char *psi_file_prefix;
 
38
}
 
39
 
 
40
namespace psi { namespace mvo {
 
41
 
 
42
int *ioff;
 
43
struct MOInfo moinfo;
 
44
struct Params params;
 
45
 
 
46
/* Max length of ioff array */
 
47
#define IOFF_MAX 32641
 
48
 
 
49
 
 
50
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
 
51
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
 
52
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
 
53
 
 
54
 
 
55
/* Function prototypes */
 
56
void init_io(int argc, char *argv[]);
 
57
void title();
 
58
void init_ioff();
 
59
void get_parameters();
 
60
void print_parameters();
 
61
void get_moinfo();
 
62
void cleanup();
 
63
void exit_io();
 
64
void get_reorder_array(void);
 
65
void get_fzc_operator(void);
 
66
void get_mvos(void);
 
67
void get_canonical(void);
 
68
void get_mp2nos(void);
 
69
/* void get_unos(void); */
 
70
}} // end namespace psi::mvo
 
71
 
 
72
using namespace psi::mvo;
 
73
 
 
74
int main(int argc, char *argv[])
 
75
{
 
76
  params.print_lvl = 1;
 
77
  init_io(argc, argv);
 
78
  title();
 
79
  get_parameters();
 
80
  get_moinfo();
 
81
  print_parameters();
 
82
  get_reorder_array();
 
83
  init_ioff();
 
84
  if (params.mp2nos) {
 
85
    get_mp2nos();
 
86
  }
 
87
  else if (params.unos) {
 
88
    /* get_unos(); */
 
89
  }
 
90
  else if (params.canonical) {
 
91
    get_canonical();
 
92
  }
 
93
  else {
 
94
    get_fzc_operator();
 
95
    get_mvos();
 
96
  }
 
97
  cleanup();
 
98
  exit_io();
 
99
  return(0);
 
100
}
 
101
 
 
102
namespace psi { namespace mvo {
 
103
 
 
104
void init_io(int argc, char *argv[])
 
105
{
 
106
  int i;
 
107
  int num_extra_args = 0;
 
108
  char **extra_args;
 
109
 
 
110
  extra_args = (char **) malloc(argc*sizeof(char *));
 
111
 
 
112
   for (i=1; i<argc; i++) {
 
113
     if (strcmp(argv[i], "--quiet") == 0) {
 
114
       params.print_lvl = 0;
 
115
     }
 
116
     else {
 
117
       extra_args[num_extra_args++] = argv[i];
 
118
     }
 
119
   }
 
120
 
 
121
  psi_start(&infile,&outfile,&psi_file_prefix,num_extra_args,extra_args,0); 
 
122
  if (params.print_lvl) tstart(outfile);
 
123
  ip_cwk_add(":MVO");
 
124
 
 
125
  psio_init(); psio_ipv1_config();
 
126
  free(extra_args);
 
127
}
 
128
 
 
129
 
 
130
void title(void)
 
131
{
 
132
  if (params.print_lvl) {
 
133
    fprintf(outfile, "\n");
 
134
    fprintf(outfile,"\t**************************************************\n");
 
135
    fprintf(outfile,"\t*     MVO: Obtain Modified Virtual Orbitals      *\n");
 
136
    fprintf(outfile,"\t*                                                *\n");
 
137
    fprintf(outfile,"\t*               C. David Sherrill                *\n");
 
138
    fprintf(outfile,"\t*                  March  2001                   *\n");
 
139
    fprintf(outfile,"\t**************************************************\n");
 
140
    fprintf(outfile, "\n");
 
141
  }
 
142
}
 
143
 
 
144
void init_ioff(void)
 
145
{
 
146
  int i;
 
147
  ioff = (int *) malloc(IOFF_MAX * sizeof(int));
 
148
  if(ioff == NULL) {
 
149
    fprintf(stderr, "(transqt): error malloc'ing ioff array\n");
 
150
    exit(0);
 
151
  }
 
152
  ioff[0] = 0;
 
153
  for(i=1; i < IOFF_MAX; i++) {
 
154
    ioff[i] = ioff[i-1] + i;
 
155
  }
 
156
}
 
157
 
 
158
void exit_io(void)
 
159
{
 
160
  psio_done();
 
161
  if (params.print_lvl) tstop(outfile);
 
162
  psi_stop(infile,outfile,psi_file_prefix);
 
163
}
 
164
 
 
165
 
 
166
void get_parameters(void)
 
167
{
 
168
  int errcod;
 
169
 
 
170
  errcod = ip_string("WFN", &(params.wfn),0);
 
171
  if(errcod == IPE_KEY_NOT_FOUND) {
 
172
    params.wfn = (char *) malloc(sizeof(char)*5);
 
173
    strcpy(params.wfn, "CCSD");
 
174
  }
 
175
 
 
176
  params.h_fzc_file = PSIF_OEI;
 
177
  errcod = ip_data("FZC_FILE","%d",&(params.h_fzc_file),0);
 
178
 
 
179
  params.print_mos = 0;
 
180
  errcod = ip_boolean("PRINT_MOS", &(params.print_mos),0);
 
181
 
 
182
  params.print_lvl = 1;
 
183
  errcod = ip_data("PRINT", "%d", &(params.print_lvl),0);
 
184
 
 
185
  params.oei_erase = 0;
 
186
  errcod = ip_boolean("OEI_ERASE",&(params.oei_erase),0); 
 
187
 
 
188
  params.fzc = 1;
 
189
  errcod = ip_boolean("FZC",&(params.fzc),0);
 
190
 
 
191
  /* remove restricted docc from RAS 1 ? */
 
192
  params.del_restr_docc = 1;
 
193
  errcod = ip_boolean("DELETE_RESTR_DOCC",&(params.del_restr_docc),0);
 
194
 
 
195
  params.mp2nos = 0;
 
196
  errcod = ip_boolean("MP2NOS",&(params.mp2nos),0);
 
197
 
 
198
  params.unos = 0;
 
199
  errcod = ip_boolean("UNOS",&(params.unos),0);
 
200
 
 
201
  if (strcmp(params.wfn, "CI")==0 || strcmp(params.wfn, "DETCAS")==0 ||
 
202
      strcmp(params.wfn, "DETCI")==0) {
 
203
    params.ras_type = 1;
 
204
  }
 
205
  else {
 
206
    params.ras_type = 0;
 
207
  }
 
208
 
 
209
  params.fzc_fock_coeff = 1.0;
 
210
  errcod = ip_data("FZC_FOCK_COEFF", "%lf", &(params.fzc_fock_coeff),0);
 
211
 
 
212
  params.fock_coeff = 0.0;
 
213
  errcod = ip_data("FOCK_COEFF", "%lf", &(params.fock_coeff),0);
 
214
 
 
215
  params.ivo = 0;
 
216
  errcod = ip_boolean("IVO", &(params.ivo), 0);
 
217
 
 
218
  params.canonical = 0;
 
219
  errcod = ip_boolean("CANONICAL", &(params.canonical), 0);
 
220
  return;
 
221
 
 
222
}
 
223
 
 
224
 
 
225
void print_parameters(void)
 
226
{
 
227
 
 
228
  if (params.print_lvl) {
 
229
    fprintf(outfile,"\n");
 
230
    fprintf(outfile,"\tInput Parameters:\n");
 
231
    fprintf(outfile,"\t-----------------\n");
 
232
    fprintf(outfile,"\tWavefunction           =  %s\n", params.wfn);
 
233
    fprintf(outfile,"\tPrint MOs              =  %s\n", 
 
234
                                (params.print_mos ? "Yes": "No"));
 
235
    fprintf(outfile,"\tErase OEI file         =  %s\n", 
 
236
                                (params.oei_erase ? "Yes": "No"));
 
237
    fprintf(outfile,"\tFrozen core            =  %s\n", 
 
238
                                (params.fzc ? "Yes": "No"));
 
239
    fprintf(outfile,"\tDelete Restricted Docc =  %s\n",
 
240
                                (params.del_restr_docc ? "Yes" : "No"));
 
241
    fprintf(outfile,"\tIVO                    =  %s\n", 
 
242
                                (params.ivo ? "Yes": "No"));
 
243
    fprintf(outfile,"\tCanonical              =  %s\n",
 
244
                                  (params.canonical ? "Yes": "No"));
 
245
    fprintf(outfile,"\tFrozen Core OEI file   =  %d\n", 
 
246
                                params.h_fzc_file);
 
247
    fprintf(outfile,"\tfzc_fock_coeff         =  %lf\n", 
 
248
                                params.fzc_fock_coeff);
 
249
    fprintf(outfile,"\tfock_coeff             =  %lf\n", 
 
250
                                params.fock_coeff);
 
251
    fprintf(outfile,"\tPrint Level            =  %d\n", params.print_lvl);
 
252
    fflush(outfile);
 
253
  }
 
254
  
 
255
  return;
 
256
}
 
257
 
 
258
void get_fzc_operator(void)
 
259
{
 
260
 
 
261
  moinfo.fzc_operator = (double *) init_array(moinfo.fzc_op_size);
 
262
  /*
 
263
  iwl_rdone(params.h_fzc_file, moinfo.fzc_operator, &(moinfo.efzc),
 
264
            ioff, moinfo.nmo, 0, 0, params.oei_erase,
 
265
            (params.print_lvl>4), outfile);
 
266
  */
 
267
  iwl_rdone(params.h_fzc_file, PSIF_MO_FZC, moinfo.fzc_operator, 
 
268
            moinfo.fzc_op_size, params.oei_erase,
 
269
            (params.print_lvl>4), outfile);
 
270
  
 
271
}
 
272
 
 
273
 
 
274
 
 
275
void get_moinfo(void)
 
276
{
 
277
  int i,j,k,h,errcod,size,row,col,p,q,offset,first_offset,last_offset,warned;
 
278
  int *tmpi, nopen;
 
279
  double **tmpmat, **so2ao;
 
280
  double N; /* number of valence electrons */
 
281
 
 
282
  chkpt_init(PSIO_OPEN_OLD);
 
283
  moinfo.nmo = chkpt_rd_nmo();
 
284
  moinfo.nso = chkpt_rd_nso();
 
285
  moinfo.nao = chkpt_rd_nao();
 
286
  moinfo.nirreps = chkpt_rd_nirreps();
 
287
  moinfo.iopen = chkpt_rd_iopen();
 
288
  moinfo.labels = chkpt_rd_irr_labs();
 
289
  moinfo.sopi   = chkpt_rd_sopi();
 
290
  moinfo.orbspi = chkpt_rd_orbspi();
 
291
  moinfo.clsdpi = chkpt_rd_clsdpi();
 
292
  moinfo.openpi = chkpt_rd_openpi();
 
293
 
 
294
  moinfo.fzc_op_size = (moinfo.nmo * (moinfo.nmo + 1)) / 2;
 
295
 
 
296
  moinfo.sosym = init_int_array(moinfo.nso);
 
297
  for (i=0,k=0; i<moinfo.nirreps; i++) {
 
298
    for (j=0; j<moinfo.sopi[i]; j++,k++) {
 
299
      moinfo.sosym[k] = i;
 
300
    }
 
301
  }
 
302
  
 
303
  moinfo.orbsym = init_int_array(moinfo.nmo);
 
304
  for (i=0,k=0; i<moinfo.nirreps; i++) {
 
305
    for (j=0; j<moinfo.orbspi[i]; j++,k++) {
 
306
      moinfo.orbsym[k] = i;
 
307
    }
 
308
  }
 
309
 
 
310
  moinfo.frdocc = init_int_array(moinfo.nirreps);
 
311
  moinfo.fruocc = init_int_array(moinfo.nirreps);
 
312
  errcod = ip_int_array("FROZEN_DOCC", moinfo.frdocc, moinfo.nirreps);
 
313
  errcod = ip_int_array("FROZEN_UOCC", moinfo.fruocc, moinfo.nirreps);
 
314
 
 
315
  moinfo.rstrdocc = init_int_array(moinfo.nirreps);
 
316
  moinfo.rstruocc = init_int_array(moinfo.nirreps);
 
317
  errcod = ip_int_array("RESTRICTED_DOCC",moinfo.rstrdocc,moinfo.nirreps);
 
318
  errcod = ip_int_array("RESTRICTED_UOCC",moinfo.rstruocc,moinfo.nirreps);
 
319
                                                                                
 
320
  /*
 
321
  if (params.treat_cor_as_fzc) {
 
322
    for (i=0; i<moinfo.nirreps; i++) {
 
323
      moinfo.frdocc[i] += moinfo.rstrdocc[i];
 
324
      moinfo.rstrdocc[i] = 0;
 
325
    }
 
326
  }
 
327
  */
 
328
 
 
329
  if (!params.fzc) {
 
330
    for (i=0; i<moinfo.nirreps; i++) {
 
331
      moinfo.frdocc[i] = 0;
 
332
    }
 
333
  }
 
334
 
 
335
  moinfo.nfzc = 0;
 
336
  moinfo.nfzv = 0;
 
337
  for (i=0; i<moinfo.nirreps; i++) {
 
338
    moinfo.nfzc += moinfo.frdocc[i];
 
339
    moinfo.nfzv += moinfo.fruocc[i];
 
340
  }
 
341
 
 
342
  moinfo.ndocc = 0;
 
343
  tmpi = init_int_array(moinfo.nirreps);
 
344
  errcod = ip_int_array("DOCC", tmpi, moinfo.nirreps);
 
345
  if (errcod == IPE_OK) {
 
346
    for (i=0,warned=0; i<moinfo.nirreps; i++) {
 
347
      if (tmpi[i] != moinfo.clsdpi[i] && !warned) {
 
348
        fprintf(outfile, "\tWarning: DOCC doesn't match PSIF_CHKPT\n");
 
349
        warned = 1;
 
350
      }
 
351
      moinfo.clsdpi[i] = tmpi[i];
 
352
      moinfo.ndocc += tmpi[i];
 
353
    }
 
354
  }
 
355
  else {
 
356
    for (i=0; i<moinfo.nirreps; i++) {
 
357
      moinfo.ndocc += moinfo.clsdpi[i];
 
358
    }
 
359
  }
 
360
 
 
361
  moinfo.nsocc = 0;
 
362
  errcod = ip_int_array("SOCC", tmpi, moinfo.nirreps);
 
363
  if (errcod == IPE_OK) {
 
364
    for (i=0,warned=0; i<moinfo.nirreps; i++) {
 
365
      if (tmpi[i] != moinfo.openpi[i] && !warned) {
 
366
        fprintf(outfile, "\tWarning: SOCC doesn't match PSIF_CHKPT\n");
 
367
        warned = 1;
 
368
      }
 
369
      moinfo.openpi[i] = tmpi[i];
 
370
      moinfo.nsocc += tmpi[i];
 
371
    }
 
372
  }
 
373
 
 
374
  moinfo.virtpi = init_int_array(moinfo.nirreps);
 
375
  for(i=0; i < moinfo.nirreps; i++) {
 
376
    moinfo.virtpi[i] = moinfo.orbspi[i]-moinfo.clsdpi[i]-moinfo.openpi[i];
 
377
  }
 
378
 
 
379
  if (params.print_lvl) {
 
380
    fprintf(outfile,"\n\tCheckpoint file parameters:\n");
 
381
    fprintf(outfile,"\t------------------\n");
 
382
    fprintf(outfile,"\tNumber of irreps = %d\n",moinfo.nirreps);
 
383
    fprintf(outfile,"\tNumber of SOs    = %d\n",moinfo.nso);
 
384
    fprintf(outfile,"\tNumber of MOs    = %d\n\n",moinfo.nmo);
 
385
    fprintf(outfile,
 
386
        "\tLabel\t# SOs\t# MOs\t# FZDC\t# DOCC\t# SOCC\t# VIRT\t# FZVR\n");
 
387
    fprintf(outfile,
 
388
        "\t-----\t-----\t-----\t------\t------\t------\t------\t------\n");
 
389
    for(i=0; i < moinfo.nirreps; i++) {
 
390
        fprintf(outfile,
 
391
           "\t %s\t   %d\t   %d\t    %d\t    %d\t    %d\t    %d\t    %d\n",
 
392
           moinfo.labels[i],moinfo.sopi[i],moinfo.orbspi[i],moinfo.frdocc[i],
 
393
           moinfo.clsdpi[i],moinfo.openpi[i],moinfo.virtpi[i],
 
394
           moinfo.fruocc[i]);
 
395
    }
 
396
    fflush(outfile);
 
397
  }
 
398
 
 
399
 
 
400
  /*
 
401
     Construct first and last index arrays for SOs: this defines the first
 
402
     absolute orbital index and last absolute orbital
 
403
     index for each irrep.  When there are no orbitals for an irrep, the
 
404
     value is -1 for first[] and -2 for last[].  Note that there must be
 
405
     basis functions in the first irrep (i.e. totally symmetric) for this to 
 
406
     work.
 
407
  */
 
408
  moinfo.first_so = init_int_array(moinfo.nirreps);
 
409
  moinfo.last_so = init_int_array(moinfo.nirreps);
 
410
  for(h=0; h < moinfo.nirreps; h++) {
 
411
    moinfo.first_so[h] = -1;
 
412
    moinfo.last_so[h] = -2;
 
413
  }
 
414
  first_offset = 0;
 
415
  last_offset = moinfo.sopi[0] - 1; 
 
416
  moinfo.first_so[0] = first_offset;
 
417
  moinfo.last_so[0] = last_offset;
 
418
  for(h=1; h < moinfo.nirreps; h++) {
 
419
    first_offset += moinfo.sopi[h-1];
 
420
    last_offset += moinfo.sopi[h];
 
421
    if(moinfo.sopi[h]) {
 
422
        moinfo.first_so[h] = first_offset;
 
423
        moinfo.last_so[h] = last_offset;
 
424
    }
 
425
  }
 
426
  
 
427
  /*
 
428
     Construct first and last index arrays: this defines the first
 
429
     absolute orbital index (Pitzer ordering) and last absolute orbital
 
430
     index for each irrep.  When there are no orbitals for an irrep, the
 
431
     value is -1 for first[] and -2 for last[].  Note that there must be
 
432
     orbitals in the first irrep (i.e. totally symmetric) for this to work.
 
433
  */
 
434
  moinfo.first = init_int_array(moinfo.nirreps);
 
435
  moinfo.last = init_int_array(moinfo.nirreps);
 
436
  for(h=0; h < moinfo.nirreps; h++) {
 
437
    moinfo.first[h] = -1;
 
438
    moinfo.last[h] = -2;
 
439
  }
 
440
  first_offset = 0;
 
441
  last_offset = moinfo.orbspi[0] - 1; 
 
442
  moinfo.first[0] = first_offset;
 
443
  moinfo.last[0] = last_offset;
 
444
  for(h=1; h < moinfo.nirreps; h++) {
 
445
    first_offset += moinfo.orbspi[h-1];
 
446
    last_offset += moinfo.orbspi[h];
 
447
    if(moinfo.orbspi[h]) {
 
448
      moinfo.first[h] = first_offset;
 
449
      moinfo.last[h] = last_offset;
 
450
    }
 
451
  }
 
452
  /*
 
453
     Construct first and last active index arrays: this defines the first
 
454
     absolute orbital index (Pitzer ordering) and last absolute orbital
 
455
     index for each irrep, excluding frozen orbitals.  When there are no
 
456
     orbitals for an irrep, the value is -1 for first[] and -2 for last[].
 
457
     Note that there must be orbitals in the first irrep (i.e. totally
 
458
     symmetric) for this to work.  
 
459
  */
 
460
  moinfo.fstact = init_int_array(moinfo.nirreps);
 
461
  moinfo.lstact = init_int_array(moinfo.nirreps);
 
462
  for(h=0; h < moinfo.nirreps; h++) {
 
463
    moinfo.fstact[h] = -1;
 
464
    moinfo.lstact[h] = -2;
 
465
  }
 
466
  first_offset = moinfo.frdocc[0];
 
467
  last_offset = moinfo.orbspi[0] - moinfo.fruocc[0] - 1; 
 
468
  moinfo.fstact[0] = first_offset;
 
469
  moinfo.lstact[0] = last_offset;
 
470
  for(h=1; h < moinfo.nirreps; h++) {
 
471
    first_offset += moinfo.orbspi[h-1]+moinfo.frdocc[h]-moinfo.frdocc[h-1];
 
472
    last_offset += moinfo.orbspi[h] - moinfo.fruocc[h] + moinfo.fruocc[h-1];
 
473
    if(moinfo.orbspi[h]) {
 
474
      moinfo.fstact[h] = first_offset;
 
475
      moinfo.lstact[h] = last_offset;
 
476
    }
 
477
  }
 
478
 
 
479
  /* Now define active[] such that frozen orbitals are taken into account */
 
480
  moinfo.active = init_int_array(moinfo.nirreps);
 
481
  for(h=0; h < moinfo.nirreps; h++) {
 
482
    moinfo.active[h] = moinfo.orbspi[h]-moinfo.frdocc[h]-moinfo.fruocc[h];
 
483
  }
 
484
 
 
485
  chkpt_close();
 
486
 
 
487
  /* in case IVO's have been asked for */
 
488
  if (params.ivo) {
 
489
 
 
490
    /* count the number of valence electrons */
 
491
    for (h=0; h < moinfo.nirreps; h++) {
 
492
      N += 2.0 * (double) (moinfo.clsdpi[h] - moinfo.frdocc[h]);
 
493
      N += (double) moinfo.openpi[h];
 
494
    }  
 
495
 
 
496
    /* use the number of valence electrons to compute the mixing ratio */
 
497
    params.fock_coeff = (N-1)/N;
 
498
 
 
499
    /*
 
500
    fprintf(outfile, "inside IVO routine now, N=%f, fock=%f\n", N,
 
501
      params.fock_coeff);
 
502
    */
 
503
 
 
504
    if (params.fock_coeff < 0.0) params.fock_coeff = 0.0;
 
505
    params.fzc_fock_coeff = 1.0 - params.fock_coeff;
 
506
  }
 
507
 
 
508
}
 
509
 
 
510
 
 
511
 
 
512
 
 
513
void get_reorder_array(void)
 
514
{
 
515
  int i, errcod;
 
516
  int j, k, l, fzv_offset;
 
517
 
 
518
  /* the following will only be set nonzero if it is a CI related WFN */
 
519
  moinfo.ras_opi = init_int_matrix(MAX_RAS_SPACES,moinfo.nirreps);
 
520
  moinfo.order = init_int_array(moinfo.nmo);
 
521
 
 
522
  if (strcmp(params.wfn, "CI") == 0 || strcmp(params.wfn, "DETCI") == 0
 
523
       || strcmp(params.wfn, "QDPT") == 0 
 
524
       || strcmp(params.wfn, "OOCCD") == 0 
 
525
       || strcmp(params.wfn, "DETCAS") == 0) {
 
526
    
 
527
    moinfo.ras_opi = init_int_matrix(MAX_RAS_SPACES,moinfo.nirreps); 
 
528
    
 
529
    if (!ras_set2(moinfo.nirreps, moinfo.nmo, params.fzc, 
 
530
         params.del_restr_docc, moinfo.orbspi, moinfo.clsdpi, moinfo.openpi,
 
531
         moinfo.frdocc, moinfo.fruocc, moinfo.rstrdocc, moinfo.rstruocc,
 
532
         moinfo.ras_opi, moinfo.order, params.ras_type, 0)) {
 
533
      fprintf(outfile, "Error in ras_set2().  Aborting.\n");
 
534
      exit(1);
 
535
    } 
 
536
  } 
 
537
  else { /* default (CC, MP2, other) */
 
538
    reorder_qt(moinfo.clsdpi, moinfo.openpi, moinfo.frdocc, moinfo.fruocc,
 
539
               moinfo.order, moinfo.orbspi, moinfo.nirreps);
 
540
  }
 
541
  
 
542
  /* construct an array to map the other direction, i.e., from correlated */
 
543
  /* to Pitzer order */
 
544
  moinfo.corr2pitz = init_int_array(moinfo.nmo);
 
545
  for (i=0; i<moinfo.nmo; i++) {
 
546
    j = moinfo.order[i];
 
547
    moinfo.corr2pitz[j] = i;
 
548
  }
 
549
 
 
550
}
 
551
 
 
552
 
 
553
void cleanup(void)
 
554
{
 
555
  free(moinfo.fzc_operator);
 
556
  free_int_matrix(moinfo.ras_opi);
 
557
}
 
558
 
 
559
void get_mvos(void)
 
560
{
 
561
 
 
562
  int h, nirreps, nvir, ntri, offset, row, col;
 
563
  int i, j, ij, iabs, jabs, icorr, jcorr, ijcorr, nocc;
 
564
  int ndv, ndoccv, *ndvph, k, colv;
 
565
  int errcod;
 
566
  double *FCvv, *FC, **evecs, *evals, **Cvv, **Cvvp, **Cnew;
 
567
  double *eig_unsrt, **scf_vector;
 
568
 
 
569
  FC = moinfo.fzc_operator;
 
570
  nirreps = moinfo.nirreps;
 
571
 
 
572
  chkpt_init(PSIO_OPEN_OLD);
 
573
  scf_vector = chkpt_rd_scf();
 
574
  eig_unsrt = chkpt_rd_evals();
 
575
  chkpt_close();
 
576
 
 
577
  /* read in which of the doccs are to be treated as virtuals */
 
578
  ndv = 0;
 
579
 
 
580
  if (ip_exist("DOCC_VIRT",0)) {
 
581
     errcod = ip_count("DOCC_VIRT",&ndv,0);
 
582
     params.docc_virt = init_int_array(ndv);
 
583
     for(k=0; k < ndv; k++) {
 
584
       params.docc_virt[k] = 0;
 
585
       errcod = ip_data("DOCC_VIRT","%d",(&(params.docc_virt[k])),1,k);
 
586
     }
 
587
  }
 
588
 
 
589
  /* finished reading in which doccs are to be treated as virtuals */
 
590
 
 
591
  /* form the vv block of FC for each irrep h */
 
592
  for (h=0,offset=0; h < nirreps; h++) {
 
593
 
 
594
  /* find out which doccs are to be treated as virtuals */
 
595
    ndoccv = 0;
 
596
    ndvph = init_int_array(moinfo.clsdpi[h]);
 
597
    for (k=0;k<moinfo.clsdpi[h];k++) {
 
598
      ndvph[k] = 0;
 
599
      }
 
600
    if (ip_exist("DOCC_VIRT",0)) {
 
601
      for (k=0;k<(ndv-1);k=(k+2)) {
 
602
        if ( (params.docc_virt[k]) == (h+1)) {
 
603
          ndoccv++;
 
604
          ndvph[(params.docc_virt[k+1])-1] = ndoccv;
 
605
        }
 
606
      }
 
607
    }
 
608
  /* ndoccv is the number of doccs to be treated as virts for irrep h */
 
609
  /* ndvph is a 1-d array of length "number of doccs for irrep h".    */
 
610
  /* it is a string of integers, 0 meaning "leave this docc alone",   */
 
611
  /* 1 at position j means that docc number j is the first docc-virt, */
 
612
  /* 2 at position j means that docc number j is the second docc-virt */
 
613
  /* etc.                                                             */
 
614
 
 
615
    nvir = moinfo.virtpi[h] + ndoccv;
 
616
    nocc = moinfo.clsdpi[h] + moinfo.openpi[h] - ndoccv;
 
617
 
 
618
    if (nvir < 1) break;
 
619
 
 
620
    if (params.print_lvl) {
 
621
      fprintf(outfile, "Working on irrep %d (%d vir, %d occ)\n", 
 
622
              h, nvir, nocc);
 
623
    }
 
624
 
 
625
    ntri = (nvir * (nvir+1))/2;
 
626
    FCvv = init_array(ntri);
 
627
    evals = init_array(ntri);
 
628
    evecs = block_matrix(nvir,nvir);
 
629
    Cvv = block_matrix(moinfo.sopi[h],nvir);
 
630
    Cvvp = block_matrix(moinfo.sopi[h],nvir);
 
631
    Cnew = block_matrix(moinfo.sopi[h],moinfo.orbspi[h]);
 
632
 
 
633
    for (i=0,ij=0; i<nvir; i++) {
 
634
      /* convert this MO index into a Correlated MO index */
 
635
      if (i>=ndoccv) {
 
636
        iabs = i + nocc + offset;
 
637
      }
 
638
      else {
 
639
        for (k=0; k<(nocc + ndoccv); k++) {
 
640
          if ( (i+1) == ndvph[k] )
 
641
             iabs = k + offset;
 
642
        }
 
643
      }
 
644
      icorr = moinfo.order[iabs];
 
645
      for (j=0; j<=i; j++,ij++) {
 
646
        if (j>=ndoccv) {
 
647
          jabs = j + nocc + offset;
 
648
        }
 
649
        else {
 
650
          for (k=0; k<(nocc + ndoccv); k++) {
 
651
            if ( (j+1) == ndvph[k] )
 
652
               jabs = k + offset; 
 
653
          }
 
654
        }
 
655
        jcorr = moinfo.order[jabs];  
 
656
        ijcorr = INDEX(icorr,jcorr);
 
657
        FCvv[ij] = params.fzc_fock_coeff * FC[ijcorr] ;
 
658
        if (i==j) FCvv[ij] += params.fock_coeff * eig_unsrt[iabs];
 
659
      }
 
660
    }
 
661
 
 
662
    if (params.print_lvl) {
 
663
      fprintf(outfile, "Old FCvv matrix for irrep %d\n", h);
 
664
      print_array(FCvv, nvir, outfile);
 
665
    }
 
666
 
 
667
    /* diagonalize the vv block of FC */
 
668
    rsp(nvir, nvir, ntri, FCvv, evals, 1, evecs, 1E-12);
 
669
 
 
670
    if (params.print_lvl) {
 
671
      fprintf(outfile, "\nDiagonalization matrix for FCvv irrep %d\n", h);
 
672
      print_mat(evecs, nvir, nvir, outfile);
 
673
 
 
674
      fprintf(outfile, "\nEigenvalues of FCvv for irrep %d\n", h);
 
675
      for (i=0; i<nvir; i++) {
 
676
        fprintf(outfile, "%lf ", evals[i]);
 
677
      }
 
678
      fprintf(outfile, "\n");
 
679
    }
 
680
 
 
681
    /* load up the vv part of SCF matrix C */
 
682
    for (i=moinfo.first_so[h],row=0; i<= moinfo.last_so[h]; i++,row++) {
 
683
      colv = 0;
 
684
      for (j=moinfo.first_so[h],k=0; j < (moinfo.first_so[h]+nocc+ndoccv); j++,k++) {
 
685
        if (ndvph[k] != 0) {
 
686
           Cvv[row][colv] = scf_vector[i][j];
 
687
           colv += 1;
 
688
        }
 
689
      }
 
690
      for (j=moinfo.first_so[h]+nocc+ndoccv,col=ndoccv; j<=moinfo.last_so[h]; j++,col++) {
 
691
        Cvv[row][col] = scf_vector[i][j];
 
692
      }
 
693
    } 
 
694
 
 
695
    if (params.print_lvl) {
 
696
      fprintf(outfile, "Old Cvv matrix for irrep %d\n", h);
 
697
      print_mat(Cvv, moinfo.sopi[h], nvir, outfile);
 
698
    }
 
699
     
 
700
    /* multiply the FCvv eigenvectors by Cvv to get new Cvv */
 
701
    mmult(Cvv, 0, evecs, 0, Cvvp, 0, moinfo.sopi[h], nvir, nvir, 0);
 
702
 
 
703
    if (params.print_lvl) {
 
704
      fprintf(outfile, "New Cvv matrix for irrep %d\n", h);
 
705
      print_mat(Cvvp, moinfo.sopi[h], nvir, outfile);
 
706
    }
 
707
 
 
708
    /* Load up the new C matrix and write it out */
 
709
 
 
710
    /* load up the occupied orbitals */
 
711
    for (i=moinfo.first_so[h],row=0; i<= moinfo.last_so[h]; i++,row++) {
 
712
      for (j=moinfo.first_so[h],col=0; j< moinfo.first_so[h]+nocc+ndoccv; j++,col++) {
 
713
        if (ndvph[col] != 0) {
 
714
          Cnew[row][col] = Cvvp[row][(ndvph[col]-1)];
 
715
        }
 
716
        else {
 
717
          Cnew[row][col] = scf_vector[i][j];
 
718
        }
 
719
      }
 
720
    }
 
721
 
 
722
    /* load up the new virtuals */
 
723
    for (i=0; i<moinfo.sopi[h]; i++) {
 
724
      for (j=ndoccv; j<nvir; j++) {
 
725
        Cnew[i][j+nocc] = Cvvp[i][j];
 
726
      }
 
727
    }
 
728
 
 
729
    if (params.print_lvl) {
 
730
      fprintf(outfile, "New C matrix for irrep %d\n", h);
 
731
      print_mat(Cnew, moinfo.sopi[h], moinfo.orbspi[h], outfile);
 
732
    }
 
733
 
 
734
    chkpt_init(PSIO_OPEN_OLD);
 
735
    chkpt_wt_scf_irrep(Cnew, h);
 
736
    chkpt_close();
 
737
 
 
738
    offset += moinfo.orbspi[h];
 
739
    free(FCvv);
 
740
    free(evals);
 
741
    free_block(evecs);
 
742
  }
 
743
 
 
744
  free(eig_unsrt);
 
745
  free_block(scf_vector);
 
746
}
 
747
 
 
748
}} // end namespace psi::mvo
 
749