3
#include <libipv1/ip_lib.h>
4
#include <libciomr/libciomr.h>
5
#include <libchkpt/chkpt.h>
7
#include "globaldefs.h"
10
void pitzer_arrays(int nirreps, int *frdocc, int *fruocc, int *orbspi,
11
int *first, int *last, int *fstact, int *lstact,
13
double *** construct_evects(int nirreps, int *active, int *orbspi,
14
int *first, int *last, int *fstact, int *lstact,
21
** Reads PSIF_CHKPT & input.dat and gets all sorts of useful information about
22
** the molecular orbitals (such as their reordering array, the docc
23
** array, frozen orbitals, etc.)
25
** Created by C. David Sherrill on 24 April 1998,
26
** based on the version in DETCI
29
void get_mo_info(void)
31
int h, i, j, k, tmp, cnt, irrep, errcod, errbad;
35
chkpt_init(PSIO_OPEN_OLD);
36
CalcInfo.nirreps = chkpt_rd_nirreps();
37
CalcInfo.nbfso = chkpt_rd_nmo();
38
CalcInfo.labels = chkpt_rd_irr_labs();
39
CalcInfo.orbs_per_irr = chkpt_rd_orbspi();
40
CalcInfo.enuc = chkpt_rd_enuc();
41
CalcInfo.efzc = chkpt_rd_efzc();
42
CalcInfo.docc = chkpt_rd_clsdpi();
43
CalcInfo.socc = chkpt_rd_openpi();
46
CalcInfo.frozen_docc = init_int_array(CalcInfo.nirreps);
47
CalcInfo.frozen_uocc = init_int_array(CalcInfo.nirreps);
48
CalcInfo.rstr_docc = init_int_array(CalcInfo.nirreps);
49
CalcInfo.rstr_uocc = init_int_array(CalcInfo.nirreps);
50
CalcInfo.pitz2ci = init_int_array(CalcInfo.nbfso);
51
CalcInfo.ras_opi = init_int_matrix(MAX_RAS_SPACES,CalcInfo.nirreps);
53
if (!ras_set2(CalcInfo.nirreps, CalcInfo.nbfso, 1, 0,
54
CalcInfo.orbs_per_irr, CalcInfo.docc, CalcInfo.socc,
55
CalcInfo.frozen_docc, CalcInfo.frozen_uocc,
56
CalcInfo.rstr_docc, CalcInfo.rstr_uocc,
57
CalcInfo.ras_opi, CalcInfo.pitz2ci, 1, 0))
59
fprintf(outfile, "Error in ras_set(). Aborting.\n");
64
/* Compute maximum number of orbitals per irrep including
65
** and not including fzv
67
CalcInfo.max_orbs_per_irrep = 0;
68
CalcInfo.max_pop_per_irrep = 0;
69
for (i=0; i<CalcInfo.nirreps; i++) {
70
if (CalcInfo.max_orbs_per_irrep < CalcInfo.orbs_per_irr[i])
71
CalcInfo.max_orbs_per_irrep = CalcInfo.orbs_per_irr[i];
72
if (CalcInfo.max_pop_per_irrep < (CalcInfo.orbs_per_irr[i] -
73
CalcInfo.frozen_uocc[i]))
74
CalcInfo.max_pop_per_irrep = CalcInfo.orbs_per_irr[i] -
75
CalcInfo.frozen_uocc[i];
79
/* construct the "ordering" array, which maps the other direction */
80
/* i.e. from a CI orbital to a Pitzer orbital */
81
CalcInfo.ci2pitz = init_int_array(CalcInfo.nbfso);
82
for (i=0; i<CalcInfo.nbfso; i++) {
83
j = CalcInfo.pitz2ci[i];
84
CalcInfo.ci2pitz[j] = i;
88
/* Set up an array to map absolute ci order to relative Pitzer order */
89
CalcInfo.ci2relpitz = init_int_array(CalcInfo.nbfso);
90
for (h=0,cnt=0; h<CalcInfo.nirreps; h++) {
91
for (i=0; i<CalcInfo.orbs_per_irr[h]; i++,cnt++) {
92
j = CalcInfo.pitz2ci[cnt];
93
CalcInfo.ci2relpitz[j] = i;
97
if (Params.print_lvl > 4) {
98
fprintf(outfile, "\nPitzer to CI order array = \n");
99
for (i=0; i<CalcInfo.nbfso; i++) {
100
fprintf(outfile, "%3d ", CalcInfo.pitz2ci[i]);
102
fprintf(outfile, "\n");
106
CalcInfo.nbstri = (CalcInfo.nbfso * (CalcInfo.nbfso + 1)) / 2 ;
107
check((CalcInfo.nbstri <= IOFF_MAX),
108
"(get_mo_info): IOFF_MAX may not large enough!");
110
/* transform orbsym vector to new MO order */
111
CalcInfo.orbsym = init_int_array(CalcInfo.nbfso);
113
for (i=0,cnt=0; i<CalcInfo.nirreps; i++) {
114
for (j=0; j<CalcInfo.orbs_per_irr[i]; j++,cnt++) {
115
k = CalcInfo.pitz2ci[cnt];
116
CalcInfo.orbsym[k] = i;
120
CalcInfo.num_fzv_orbs = 0;
121
for (i=0; i<CalcInfo.nirreps; i++)
122
CalcInfo.num_fzv_orbs += CalcInfo.frozen_uocc[i];
124
CalcInfo.npop = CalcInfo.nbfso - CalcInfo.num_fzv_orbs;
126
CalcInfo.num_fzc_orbs = 0;
127
CalcInfo.num_cor_orbs = 0;
129
for (i=0; i<CalcInfo.nirreps; i++) {
130
j = CalcInfo.frozen_docc[i];
131
CalcInfo.num_fzc_orbs += j;
135
for (i=0; i<CalcInfo.nirreps; i++) {
136
CalcInfo.num_cor_orbs += CalcInfo.frozen_docc[i];
140
/* construct the CalcInfo.ras_orbs array (may not be of any use now) */
142
CalcInfo.fzc_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nbfso);
143
CalcInfo.fzv_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nbfso);
144
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
145
for (j=0; j<CalcInfo.frozen_docc[irrep]; j++)
146
CalcInfo.fzc_orbs[irrep][j] = cnt++;
148
CalcInfo.ras_orbs = (int ***) malloc (MAX_RAS_SPACES * sizeof(int **));
149
for (i=0; i<MAX_RAS_SPACES; i++) {
150
CalcInfo.ras_orbs[i] = init_int_matrix(CalcInfo.nirreps,
152
for (irrep=0; irrep<CalcInfo.nirreps; irrep++) {
153
for (j=0; j<CalcInfo.ras_opi[i][irrep]; j++) {
154
CalcInfo.ras_orbs[i][irrep][j] = cnt++;
159
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
160
for (j=0; j<CalcInfo.frozen_uocc[irrep]; j++)
161
CalcInfo.fzv_orbs[irrep][j] = cnt++;
165
/* get the Pitzer arrays first, last, fstact, lstact, and active */
166
CalcInfo.first = init_int_array(CalcInfo.nirreps);
167
CalcInfo.last = init_int_array(CalcInfo.nirreps);
168
CalcInfo.fstact = init_int_array(CalcInfo.nirreps);
169
CalcInfo.lstact = init_int_array(CalcInfo.nirreps);
170
CalcInfo.active = init_int_array(CalcInfo.nirreps);
171
pitzer_arrays(CalcInfo.nirreps, CalcInfo.frozen_docc, CalcInfo.frozen_uocc,
172
CalcInfo.orbs_per_irr, CalcInfo.first, CalcInfo.last,
173
CalcInfo.fstact, CalcInfo.lstact, CalcInfo.active);
175
/* allocate memory to store the MO coefficient matrix symm blocked */
177
CalcInfo.mo_coeffs = (double ***) malloc(CalcInfo.nirreps *
179
for (irrep=0; irrep<CalcInfo.nirreps; irrep++) {
180
i = CalcInfo.orbs_per_irr[irrep];
182
CalcInfo.mo_coeffs[irrep] = block_matrix(i,i);
185
if (Params.print_lvl > 0) {
186
fprintf(outfile, "ORBITALS:");
187
fprintf(outfile, "\n FROZEN_DOCC = ");
188
for (i=0; i<CalcInfo.nirreps; i++) {
189
fprintf(outfile, "%2d ", CalcInfo.frozen_docc[i]);
191
fprintf(outfile, "\n RESTR_DOCC = ");
192
for (i=0; i<CalcInfo.nirreps; i++) {
193
fprintf(outfile, "%2d ", CalcInfo.rstr_docc[i]);
195
fprintf(outfile, "\n DOCC = ");
196
for (i=0; i<CalcInfo.nirreps; i++) {
197
fprintf(outfile, "%2d ", CalcInfo.docc[i]);
199
fprintf(outfile, "\n SOCC = ");
200
for (i=0; i<CalcInfo.nirreps; i++) {
201
fprintf(outfile, "%2d ", CalcInfo.socc[i]);
203
fprintf(outfile, "\n RESTR_UOCC = ");
204
for (i=0; i<CalcInfo.nirreps; i++) {
205
fprintf(outfile, "%2d ", CalcInfo.rstr_uocc[i]);
207
fprintf(outfile, "\n FROZEN_UOCC = ");
208
for (i=0; i<CalcInfo.nirreps; i++) {
209
fprintf(outfile, "%2d ", CalcInfo.frozen_uocc[i]);
212
for (i=0; i<MAX_RAS_SPACES; i++) {
213
fprintf(outfile, "\n RAS %d = ",i+1);
214
for (j=0; j<CalcInfo.nirreps; j++) {
215
fprintf(outfile, "%2d ", CalcInfo.ras_opi[i][j]);
218
fprintf(outfile, "\n");
220
fprintf(outfile, " MOL ORBS = %6d\n", CalcInfo.nbfso);
221
fprintf(outfile, " FROZEN CORE = %6d RESTR CORE = %6d\n",
222
CalcInfo.num_fzc_orbs, CalcInfo.num_cor_orbs);
223
fprintf(outfile, "\n");
232
** Form the first/last/active arrays for the orbitals in Pitzer order
233
** Based on code taken from TRANSQT
238
void pitzer_arrays(int nirreps, int *frdocc, int *fruocc, int *orbspi,
239
int *first, int *last, int *fstact, int *lstact,int *active)
243
int first_offset, last_offset;
246
* Construct first and last index arrays: this defines the first
247
* absolute orbital index (Pitzer ordering) and last absolute orbital
248
* index for each irrep. When there are no orbitals for an irrep, the
249
* value is -1 for first[] and -2 for last[]. Note that there must be
250
* orbitals in the first irrep (i.e. totally symmetric) for this to work.
252
for(h=0; h < nirreps; h++) {
258
last_offset = orbspi[0] - 1;
259
first[0] = first_offset;
260
last[0] = last_offset;
262
for(h=1; h < nirreps; h++) {
263
first_offset += orbspi[h-1];
264
last_offset += orbspi[h];
266
first[h] = first_offset;
267
last[h] = last_offset;
272
* Construct first and last active index arrays: this defines the first
273
* absolute orbital index (Pitzer ordering) and last absolute orbital
274
* index for each irrep, excluding frozen orbitals. When there are no
275
* orbitals for an irrep, the value is -1 for first[] and -2 for last[].
276
* Note that there must be orbitals in the first irrep (i.e. totally
277
* symmetric) for this to work.
279
for(h=0; h < nirreps; h++) {
284
first_offset = frdocc[0];
285
last_offset = orbspi[0] - fruocc[0] - 1;
286
fstact[0] = first_offset;
287
lstact[0] = last_offset;
289
for(h=1; h < nirreps; h++) {
290
first_offset += orbspi[h-1]+frdocc[h]-frdocc[h-1];
291
last_offset += orbspi[h] - fruocc[h] + fruocc[h-1];
293
fstact[h] = first_offset;
294
lstact[h] = last_offset;
298
/* Now define active[] such that frozen orbitals are taken into account */
299
for(h=0; h < nirreps; h++) {
300
active[h] = orbspi[h]-frdocc[h]-fruocc[h];
310
** Read in the molecular orbital matrix from PSIF_CHKPT and put them in CalcInfo
312
void read_cur_orbs(void)
314
int i, j, h, dim, nirreps;
317
nirreps = CalcInfo.nirreps;
319
chkpt_init(PSIO_OPEN_OLD);
320
for (h=0; h<nirreps; h++) {
321
dim = CalcInfo.orbs_per_irr[h];
322
if (dim==0) continue;
323
tmat = chkpt_rd_scf_irrep(h);
324
for (i=0; i<dim; i++)
325
for (j=0; j<dim; j++)
326
CalcInfo.mo_coeffs[h][i][j] = tmat[i][j];
338
** This function copies SCF eigenvectors from moinfo.scf_vector into
339
** an array of matrices. Columns corresponding to inactive orbitals
340
** may be deleted if desired.
342
** Code taken from TRANSQT
346
double *** construct_evects(int nirreps, int *active, int *orbspi,
347
int *first, int *last, int *fstact, int *lstact,
351
int h, row, col, p, q;
354
evects = (double ***) malloc(nirreps * sizeof(double **));
356
for (h=0; h<nirreps; h++) {
358
evects[h] = block_matrix(orbspi[h],active[h]);
360
for(p=first[h]; p <= last[h]; p++) {
362
for(q=fstact[h]; q <= lstact[h]; q++) {
364
evects[h][row][col] = CalcInfo.mo_matrix[p][q];
369
fprintf(outfile,"\n\tMolecular Orbitals for Irrep %s\n",
371
print_mat(evects[h],orbspi[h],active[h],outfile);
384
void destruct_evects(int nirreps, double ***evects)
388
for (h=0; h<nirreps; h++)
389
if (evects[h] != NULL) free_block(evects[h]);
added
removed
src/bin/detcas/get_mo_info.c
