4
** This program reads the frozen core operator from TRANSQT
5
** output ** and diagonalizes the virtual-virtual block to make
6
** Modified Virtual Orbitals (MVO's) as described in
7
** C. W. Bauschlicher, J. Chem. Phys. 72, 880 (1980)
10
** Georgia Institute of Technology
18
#include <libipv1/ip_lib.h>
19
#include <libpsio/psio.h>
20
#include <libciomr/libciomr.h>
21
#include <libchkpt/chkpt.h>
24
#include <libiwl/iwl.h>
30
/* First definitions of globals */
31
FILE *infile, *outfile;
32
char *psi_file_prefix;
37
/* Max length of ioff array */
38
#define IOFF_MAX 32641
41
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
42
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
43
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
46
/* Function prototypes */
47
void init_io(int argc, char *argv[]);
50
void get_parameters();
51
void print_parameters();
55
void get_reorder_array(void);
56
void get_fzc_operator(void);
58
void get_canonical(void);
59
void get_mp2nos(void);
61
main(int argc, char *argv[])
74
else if (params.canonical) {
86
void init_io(int argc, char *argv[])
89
int num_extra_args = 0;
92
extra_args = (char **) malloc(argc*sizeof(char *));
94
for (i=1; i<argc; i++) {
95
if (strcmp(argv[i], "--quiet") == 0) {
99
extra_args[num_extra_args++] = argv[i];
103
psi_start(num_extra_args,extra_args,0);
104
if (params.print_lvl) tstart(outfile);
114
if (params.print_lvl) {
115
fprintf(outfile, "\n");
116
fprintf(outfile,"\t**************************************************\n");
117
fprintf(outfile,"\t* MVO: Obtain Modified Virtual Orbitals *\n");
118
fprintf(outfile,"\t* *\n");
119
fprintf(outfile,"\t* C. David Sherrill *\n");
120
fprintf(outfile,"\t* March 2001 *\n");
121
fprintf(outfile,"\t**************************************************\n");
122
fprintf(outfile, "\n");
129
ioff = (int *) malloc(IOFF_MAX * sizeof(int));
131
fprintf(stderr, "(transqt): error malloc'ing ioff array\n");
135
for(i=1; i < IOFF_MAX; i++) {
136
ioff[i] = ioff[i-1] + i;
143
if (params.print_lvl) tstop(outfile);
148
void get_parameters(void)
152
errcod = ip_string("WFN", &(params.wfn),0);
153
if(errcod == IPE_KEY_NOT_FOUND) {
154
params.wfn = (char *) malloc(sizeof(char)*5);
155
strcpy(params.wfn, "CCSD");
158
params.h_fzc_file = PSIF_OEI;
159
errcod = ip_data("FZC_FILE","%d",&(params.h_fzc_file),0);
161
params.print_mos = 0;
162
errcod = ip_boolean("PRINT_MOS", &(params.print_mos),0);
164
params.print_lvl = 1;
165
errcod = ip_data("PRINT", "%d", &(params.print_lvl),0);
167
params.oei_erase = 0;
168
errcod = ip_boolean("OEI_ERASE",&(params.oei_erase),0);
171
errcod = ip_boolean("FZC",&(params.fzc),0);
174
errcod = ip_boolean("MP2NOS",&(params.mp2nos),0);
176
if (strcmp(params.wfn, "CI")==0 || strcmp(params.wfn, "DETCAS")==0 ||
177
strcmp(params.wfn, "DETCI")==0) {
184
params.fzc_fock_coeff = 1.0;
185
errcod = ip_data("FZC_FOCK_COEFF", "%lf", &(params.fzc_fock_coeff),0);
187
params.fock_coeff = 0.0;
188
errcod = ip_data("FOCK_COEFF", "%lf", &(params.fock_coeff),0);
191
errcod = ip_boolean("IVO", &(params.ivo), 0);
193
params.canonical = 0;
194
errcod = ip_boolean("CANONICAL", &(params.canonical), 0);
200
void print_parameters(void)
203
if (params.print_lvl) {
204
fprintf(outfile,"\tInput Parameters:\n");
205
fprintf(outfile,"\t-----------------\n");
206
fprintf(outfile,"\tWavefunction = %s\n", params.wfn);
207
fprintf(outfile,"\tPrint MOs = %s\n",
208
(params.print_mos ? "Yes": "No"));
209
fprintf(outfile,"\tErase OEI file = %s\n",
210
(params.oei_erase ? "Yes": "No"));
211
fprintf(outfile,"\tFrozen core = %s\n",
212
(params.fzc ? "Yes": "No"));
213
fprintf(outfile,"\tIVO = %s\n",
214
(params.ivo ? "Yes": "No"));
215
fprintf(outfile,"\tCanonical = %s\n",
216
(params.canonical ? "Yes": "No"));
217
fprintf(outfile,"\tFrozen Core OEI file = %d\n",
219
fprintf(outfile,"\tfzc_fock_coeff = %lf\n",
220
params.fzc_fock_coeff);
221
fprintf(outfile,"\tfock_coeff = %lf\n",
223
fprintf(outfile,"\tPrint Level = %d\n", params.print_lvl);
230
void get_fzc_operator(void)
233
moinfo.fzc_operator = (double *) init_array(moinfo.fzc_op_size);
235
iwl_rdone(params.h_fzc_file, moinfo.fzc_operator, &(moinfo.efzc),
236
ioff, moinfo.nmo, 0, 0, params.oei_erase,
237
(params.print_lvl>4), outfile);
239
iwl_rdone(params.h_fzc_file, PSIF_MO_FZC, moinfo.fzc_operator,
240
moinfo.fzc_op_size, params.oei_erase,
241
(params.print_lvl>4), outfile);
247
void get_moinfo(void)
249
int i,j,k,h,errcod,size,row,col,p,q,offset,first_offset,last_offset,warned;
251
double **tmpmat, **so2ao;
252
double N; /* number of valence electrons */
254
chkpt_init(PSIO_OPEN_OLD);
255
moinfo.nmo = chkpt_rd_nmo();
256
moinfo.nso = chkpt_rd_nso();
257
moinfo.nao = chkpt_rd_nao();
258
moinfo.nirreps = chkpt_rd_nirreps();
259
moinfo.iopen = chkpt_rd_iopen();
260
moinfo.labels = chkpt_rd_irr_labs();
261
moinfo.sopi = chkpt_rd_sopi();
262
moinfo.orbspi = chkpt_rd_orbspi();
263
moinfo.clsdpi = chkpt_rd_clsdpi();
264
moinfo.openpi = chkpt_rd_openpi();
266
moinfo.fzc_op_size = (moinfo.nmo * (moinfo.nmo + 1)) / 2;
268
moinfo.sosym = init_int_array(moinfo.nso);
269
for (i=0,k=0; i<moinfo.nirreps; i++) {
270
for (j=0; j<moinfo.sopi[i]; j++,k++) {
275
moinfo.orbsym = init_int_array(moinfo.nmo);
276
for (i=0,k=0; i<moinfo.nirreps; i++) {
277
for (j=0; j<moinfo.orbspi[i]; j++,k++) {
278
moinfo.orbsym[k] = i;
282
moinfo.frdocc = init_int_array(moinfo.nirreps);
283
moinfo.fruocc = init_int_array(moinfo.nirreps);
284
errcod = ip_int_array("FROZEN_DOCC", moinfo.frdocc, moinfo.nirreps);
285
errcod = ip_int_array("FROZEN_UOCC", moinfo.fruocc, moinfo.nirreps);
288
for (i=0; i<moinfo.nirreps; i++) {
289
moinfo.frdocc[i] = 0;
295
for (i=0; i<moinfo.nirreps; i++) {
296
moinfo.nfzc += moinfo.frdocc[i];
297
moinfo.nfzv += moinfo.fruocc[i];
301
tmpi = init_int_array(moinfo.nirreps);
302
errcod = ip_int_array("DOCC", tmpi, moinfo.nirreps);
303
if (errcod == IPE_OK) {
304
for (i=0,warned=0; i<moinfo.nirreps; i++) {
305
if (tmpi[i] != moinfo.clsdpi[i] && !warned) {
306
fprintf(outfile, "\tWarning: DOCC doesn't match PSIF_CHKPT\n");
309
moinfo.clsdpi[i] = tmpi[i];
310
moinfo.ndocc += tmpi[i];
314
for (i=0; i<moinfo.nirreps; i++) {
315
moinfo.ndocc += moinfo.clsdpi[i];
320
errcod = ip_int_array("SOCC", tmpi, moinfo.nirreps);
321
if (errcod == IPE_OK) {
322
for (i=0,warned=0; i<moinfo.nirreps; i++) {
323
if (tmpi[i] != moinfo.openpi[i] && !warned) {
324
fprintf(outfile, "\tWarning: SOCC doesn't match PSIF_CHKPT\n");
327
moinfo.openpi[i] = tmpi[i];
328
moinfo.nsocc += tmpi[i];
332
moinfo.virtpi = init_int_array(moinfo.nirreps);
333
for(i=0; i < moinfo.nirreps; i++) {
334
moinfo.virtpi[i] = moinfo.orbspi[i]-moinfo.clsdpi[i]-moinfo.openpi[i];
337
if (params.print_lvl) {
338
fprintf(outfile,"\n\tCheckpoint file parameters:\n");
339
fprintf(outfile,"\t------------------\n");
340
fprintf(outfile,"\tNumber of irreps = %d\n",moinfo.nirreps);
341
fprintf(outfile,"\tNumber of SOs = %d\n",moinfo.nso);
342
fprintf(outfile,"\tNumber of MOs = %d\n\n",moinfo.nmo);
344
"\tLabel\t# SOs\t# MOs\t# FZDC\t# DOCC\t# SOCC\t# VIRT\t# FZVR\n");
346
"\t-----\t-----\t-----\t------\t------\t------\t------\t------\n");
347
for(i=0; i < moinfo.nirreps; i++) {
349
"\t %s\t %d\t %d\t %d\t %d\t %d\t %d\t %d\n",
350
moinfo.labels[i],moinfo.sopi[i],moinfo.orbspi[i],moinfo.frdocc[i],
351
moinfo.clsdpi[i],moinfo.openpi[i],moinfo.virtpi[i],
359
Construct first and last index arrays for SOs: this defines the first
360
absolute orbital index and last absolute orbital
361
index for each irrep. When there are no orbitals for an irrep, the
362
value is -1 for first[] and -2 for last[]. Note that there must be
363
basis functions in the first irrep (i.e. totally symmetric) for this to
366
moinfo.first_so = init_int_array(moinfo.nirreps);
367
moinfo.last_so = init_int_array(moinfo.nirreps);
368
for(h=0; h < moinfo.nirreps; h++) {
369
moinfo.first_so[h] = -1;
370
moinfo.last_so[h] = -2;
373
last_offset = moinfo.sopi[0] - 1;
374
moinfo.first_so[0] = first_offset;
375
moinfo.last_so[0] = last_offset;
376
for(h=1; h < moinfo.nirreps; h++) {
377
first_offset += moinfo.sopi[h-1];
378
last_offset += moinfo.sopi[h];
380
moinfo.first_so[h] = first_offset;
381
moinfo.last_so[h] = last_offset;
386
Construct first and last index arrays: this defines the first
387
absolute orbital index (Pitzer ordering) and last absolute orbital
388
index for each irrep. When there are no orbitals for an irrep, the
389
value is -1 for first[] and -2 for last[]. Note that there must be
390
orbitals in the first irrep (i.e. totally symmetric) for this to work.
392
moinfo.first = init_int_array(moinfo.nirreps);
393
moinfo.last = init_int_array(moinfo.nirreps);
394
for(h=0; h < moinfo.nirreps; h++) {
395
moinfo.first[h] = -1;
399
last_offset = moinfo.orbspi[0] - 1;
400
moinfo.first[0] = first_offset;
401
moinfo.last[0] = last_offset;
402
for(h=1; h < moinfo.nirreps; h++) {
403
first_offset += moinfo.orbspi[h-1];
404
last_offset += moinfo.orbspi[h];
405
if(moinfo.orbspi[h]) {
406
moinfo.first[h] = first_offset;
407
moinfo.last[h] = last_offset;
411
Construct first and last active index arrays: this defines the first
412
absolute orbital index (Pitzer ordering) and last absolute orbital
413
index for each irrep, excluding frozen orbitals. When there are no
414
orbitals for an irrep, the value is -1 for first[] and -2 for last[].
415
Note that there must be orbitals in the first irrep (i.e. totally
416
symmetric) for this to work.
418
moinfo.fstact = init_int_array(moinfo.nirreps);
419
moinfo.lstact = init_int_array(moinfo.nirreps);
420
for(h=0; h < moinfo.nirreps; h++) {
421
moinfo.fstact[h] = -1;
422
moinfo.lstact[h] = -2;
424
first_offset = moinfo.frdocc[0];
425
last_offset = moinfo.orbspi[0] - moinfo.fruocc[0] - 1;
426
moinfo.fstact[0] = first_offset;
427
moinfo.lstact[0] = last_offset;
428
for(h=1; h < moinfo.nirreps; h++) {
429
first_offset += moinfo.orbspi[h-1]+moinfo.frdocc[h]-moinfo.frdocc[h-1];
430
last_offset += moinfo.orbspi[h] - moinfo.fruocc[h] + moinfo.fruocc[h-1];
431
if(moinfo.orbspi[h]) {
432
moinfo.fstact[h] = first_offset;
433
moinfo.lstact[h] = last_offset;
437
/* Now define active[] such that frozen orbitals are taken into account */
438
moinfo.active = init_int_array(moinfo.nirreps);
439
for(h=0; h < moinfo.nirreps; h++) {
440
moinfo.active[h] = moinfo.orbspi[h]-moinfo.frdocc[h]-moinfo.fruocc[h];
445
/* in case IVO's have been asked for */
448
/* count the number of valence electrons */
449
for (h=0; h < moinfo.nirreps; h++) {
450
N += 2.0 * (double) (moinfo.clsdpi[h] - moinfo.frdocc[h]);
451
N += (double) moinfo.openpi[h];
454
/* use the number of valence electrons to compute the mixing ratio */
455
params.fock_coeff = (N-1)/N;
458
fprintf(outfile, "inside IVO routine now, N=%f, fock=%f\n", N,
462
if (params.fock_coeff < 0.0) params.fock_coeff = 0.0;
463
params.fzc_fock_coeff = 1.0 - params.fock_coeff;
471
void get_reorder_array(void)
474
int j, k, l, fzv_offset;
476
moinfo.order = init_int_array(moinfo.nmo);
478
/* the following will only be set nonzero if it is a CI related WFN */
479
moinfo.ras_opi = init_int_matrix(4,moinfo.nirreps);
481
if (strcmp(params.wfn, "CI") == 0 || strcmp(params.wfn, "DETCI") == 0
482
|| strcmp(params.wfn, "QDPT") == 0
483
|| strcmp(params.wfn, "OOCCD") == 0
484
|| strcmp(params.wfn, "DETCAS") == 0) {
486
moinfo.ras_opi = init_int_matrix(4,moinfo.nirreps);
488
if (!ras_set(moinfo.nirreps, moinfo.nmo, params.fzc, moinfo.orbspi,
489
moinfo.clsdpi, moinfo.openpi, moinfo.frdocc, moinfo.fruocc,
490
moinfo.ras_opi, moinfo.order, params.ras_type)) {
491
fprintf(outfile, "Error in ras_set(). Aborting.\n");
496
else { /* default (CC, MP2, other) */
497
reorder_qt(moinfo.clsdpi, moinfo.openpi, moinfo.frdocc, moinfo.fruocc,
498
moinfo.order, moinfo.orbspi, moinfo.nirreps);
501
/* construct an array to map the other direction, i.e., from correlated */
502
/* to Pitzer order */
503
moinfo.corr2pitz = init_int_array(moinfo.nmo);
504
for (i=0; i<moinfo.nmo; i++) {
506
moinfo.corr2pitz[j] = i;
514
free(moinfo.fzc_operator);
515
free_int_matrix(moinfo.ras_opi, 4);
521
int h, nirreps, nvir, ntri, offset, row, col;
522
int i, j, ij, iabs, jabs, icorr, jcorr, ijcorr, nocc;
523
int ndv, ndoccv, *ndvph, k, colv;
525
double *FCvv, *FC, **evecs, *evals, **Cvv, **Cvvp, **Cnew;
526
double *eig_unsrt, **scf_vector;
528
FC = moinfo.fzc_operator;
529
nirreps = moinfo.nirreps;
531
chkpt_init(PSIO_OPEN_OLD);
532
scf_vector = chkpt_rd_scf();
533
eig_unsrt = chkpt_rd_evals();
536
/* read in which of the doccs are to be treated as virtuals */
539
if (ip_exist("DOCC_VIRT",0)) {
540
errcod = ip_count("DOCC_VIRT",&ndv,0);
541
params.docc_virt = init_int_array(ndv);
542
for(k=0; k < ndv; k++) {
543
params.docc_virt[k] = 0;
544
errcod = ip_data("DOCC_VIRT","%d",(&(params.docc_virt[k])),1,k);
548
/* finished reading in which doccs are to be treated as virtuals */
550
/* form the vv block of FC for each irrep h */
551
for (h=0,offset=0; h < nirreps; h++) {
553
/* find out which doccs are to be treated as virtuals */
555
ndvph = init_int_array(moinfo.clsdpi[h]);
556
for (k=0;k<moinfo.clsdpi[h];k++) {
559
if (ip_exist("DOCC_VIRT",0)) {
560
for (k=0;k<(ndv-1);k=(k+2)) {
561
if ( (params.docc_virt[k]) == (h+1)) {
563
ndvph[(params.docc_virt[k+1])-1] = ndoccv;
567
/* ndoccv is the number of doccs to be treated as virts for irrep h */
568
/* ndvph is a 1-d array of length "number of doccs for irrep h". */
569
/* it is a string of integers, 0 meaning "leave this docc alone", */
570
/* 1 at position j means that docc number j is the first docc-virt, */
571
/* 2 at position j means that docc number j is the second docc-virt */
574
nvir = moinfo.virtpi[h] + ndoccv;
575
nocc = moinfo.clsdpi[h] + moinfo.openpi[h] - ndoccv;
579
if (params.print_lvl) {
580
fprintf(outfile, "Working on irrep %d (%d vir, %d occ)\n",
584
ntri = (nvir * (nvir+1))/2;
585
FCvv = init_array(ntri);
586
evals = init_array(ntri);
587
evecs = block_matrix(nvir,nvir);
588
Cvv = block_matrix(moinfo.sopi[h],nvir);
589
Cvvp = block_matrix(moinfo.sopi[h],nvir);
590
Cnew = block_matrix(moinfo.sopi[h],moinfo.orbspi[h]);
592
for (i=0,ij=0; i<nvir; i++) {
593
/* convert this MO index into a Correlated MO index */
595
iabs = i + nocc + offset;
598
for (k=0; k<(nocc + ndoccv); k++) {
599
if ( (i+1) == ndvph[k] )
603
icorr = moinfo.order[iabs];
604
for (j=0; j<=i; j++,ij++) {
606
jabs = j + nocc + offset;
609
for (k=0; k<(nocc + ndoccv); k++) {
610
if ( (j+1) == ndvph[k] )
614
jcorr = moinfo.order[jabs];
615
ijcorr = INDEX(icorr,jcorr);
616
FCvv[ij] = params.fzc_fock_coeff * FC[ijcorr] ;
617
if (i==j) FCvv[ij] += params.fock_coeff * eig_unsrt[iabs];
621
if (params.print_lvl) {
622
fprintf(outfile, "Old FCvv matrix for irrep %d\n", h);
623
print_array(FCvv, nvir, outfile);
626
/* diagonalize the vv block of FC */
627
rsp(nvir, nvir, ntri, FCvv, evals, 1, evecs, 1E-12);
629
if (params.print_lvl) {
630
fprintf(outfile, "\nDiagonalization matrix for FCvv irrep %d\n", h);
631
print_mat(evecs, nvir, nvir, outfile);
633
fprintf(outfile, "\nEigenvalues of FCvv for irrep %d\n", h);
634
for (i=0; i<nvir; i++) {
635
fprintf(outfile, "%lf ", evals[i]);
637
fprintf(outfile, "\n");
640
/* load up the vv part of SCF matrix C */
641
for (i=moinfo.first_so[h],row=0; i<= moinfo.last_so[h]; i++,row++) {
643
for (j=moinfo.first_so[h],k=0; j < (moinfo.first_so[h]+nocc+ndoccv); j++,k++) {
645
Cvv[row][colv] = scf_vector[i][j];
649
for (j=moinfo.first_so[h]+nocc+ndoccv,col=ndoccv; j<=moinfo.last_so[h]; j++,col++) {
650
Cvv[row][col] = scf_vector[i][j];
654
if (params.print_lvl) {
655
fprintf(outfile, "Old Cvv matrix for irrep %d\n", h);
656
print_mat(Cvv, moinfo.sopi[h], nvir, outfile);
659
/* multiply the FCvv eigenvectors by Cvv to get new Cvv */
660
mmult(Cvv, 0, evecs, 0, Cvvp, 0, moinfo.sopi[h], nvir, nvir, 0);
662
if (params.print_lvl) {
663
fprintf(outfile, "New Cvv matrix for irrep %d\n", h);
664
print_mat(Cvvp, moinfo.sopi[h], nvir, outfile);
667
/* Load up the new C matrix and write it out */
669
/* load up the occupied orbitals */
670
for (i=moinfo.first_so[h],row=0; i<= moinfo.last_so[h]; i++,row++) {
671
for (j=moinfo.first_so[h],col=0; j< moinfo.first_so[h]+nocc+ndoccv; j++,col++) {
672
if (ndvph[col] != 0) {
673
Cnew[row][col] = Cvvp[row][(ndvph[col]-1)];
676
Cnew[row][col] = scf_vector[i][j];
681
/* load up the new virtuals */
682
for (i=0; i<moinfo.sopi[h]; i++) {
683
for (j=ndoccv; j<nvir; j++) {
684
Cnew[i][j+nocc] = Cvvp[i][j];
688
if (params.print_lvl) {
689
fprintf(outfile, "New C matrix for irrep %d\n", h);
690
print_mat(Cnew, moinfo.sopi[h], moinfo.orbspi[h], outfile);
693
chkpt_init(PSIO_OPEN_OLD);
694
chkpt_wt_scf_irrep(Cnew, h);
697
offset += moinfo.orbspi[h];
704
free_block(scf_vector);
added
removed
src/bin/mvo/mvo.c
