3
\brief Enter brief description of file here
7
#include <libipv1/ip_lib.h>
8
#include <libciomr/libciomr.h>
9
#include <libchkpt/chkpt.h>
11
#include "globaldefs.h"
14
namespace psi { namespace detcas {
16
void pitzer_arrays(int nirreps, int *frdocc, int *fruocc, int *orbspi,
17
int *first, int *last, int *fstact, int *lstact,
19
double *** construct_evects(int nirreps, int *active, int *orbspi,
20
int *first, int *last, int *fstact, int *lstact,
22
extern void check(int a, const char *errmsg);
27
** Reads PSIF_CHKPT & input.dat and gets all sorts of useful information about
28
** the molecular orbitals (such as their reordering array, the docc
29
** array, frozen orbitals, etc.)
31
** Created by C. David Sherrill on 24 April 1998,
32
** based on the version in DETCI
35
void get_mo_info(void)
37
int h, i, j, k, tmp, cnt, irrep, errcod, errbad;
41
/* set these to NULL so we'll know which one(s) to free in cleanup */
42
CalcInfo.mo_hess = NULL;
43
CalcInfo.mo_hess_diag = NULL;
45
/* information from checkpoint file */
46
chkpt_init(PSIO_OPEN_OLD);
47
CalcInfo.nirreps = chkpt_rd_nirreps();
48
CalcInfo.nmo = chkpt_rd_nmo();
49
CalcInfo.nbfso = chkpt_rd_nmo(); /* change to nbfso after conversion */
50
CalcInfo.labels = chkpt_rd_irr_labs();
51
CalcInfo.orbs_per_irr = chkpt_rd_orbspi();
52
CalcInfo.enuc = chkpt_rd_enuc();
53
CalcInfo.efzc = chkpt_rd_efzc();
54
CalcInfo.docc = chkpt_rd_clsdpi();
55
CalcInfo.socc = chkpt_rd_openpi();
58
CalcInfo.frozen_docc = init_int_array(CalcInfo.nirreps);
59
CalcInfo.frozen_uocc = init_int_array(CalcInfo.nirreps);
60
CalcInfo.rstr_docc = init_int_array(CalcInfo.nirreps);
61
CalcInfo.rstr_uocc = init_int_array(CalcInfo.nirreps);
62
CalcInfo.pitz2ci = init_int_array(CalcInfo.nmo);
63
CalcInfo.ras_opi = init_int_matrix(MAX_RAS_SPACES,CalcInfo.nirreps);
65
if (!ras_set2(CalcInfo.nirreps, CalcInfo.nmo, 1, 1,
66
CalcInfo.orbs_per_irr, CalcInfo.docc, CalcInfo.socc,
67
CalcInfo.frozen_docc, CalcInfo.frozen_uocc,
68
CalcInfo.rstr_docc, CalcInfo.rstr_uocc,
69
CalcInfo.ras_opi, CalcInfo.pitz2ci, 1, 0))
71
fprintf(outfile, "Error in ras_set(). Aborting.\n");
76
/* Compute maximum number of orbitals per irrep including
77
** and not including fzv
79
CalcInfo.max_orbs_per_irrep = 0;
80
CalcInfo.max_pop_per_irrep = 0;
81
for (i=0; i<CalcInfo.nirreps; i++) {
82
if (CalcInfo.max_orbs_per_irrep < CalcInfo.orbs_per_irr[i])
83
CalcInfo.max_orbs_per_irrep = CalcInfo.orbs_per_irr[i];
84
if (CalcInfo.max_pop_per_irrep < (CalcInfo.orbs_per_irr[i] -
85
CalcInfo.frozen_uocc[i]))
86
CalcInfo.max_pop_per_irrep = CalcInfo.orbs_per_irr[i] -
87
CalcInfo.frozen_uocc[i];
91
/* construct the "ordering" array, which maps the other direction */
92
/* i.e. from a CI orbital to a Pitzer orbital */
93
CalcInfo.ci2pitz = init_int_array(CalcInfo.nmo);
94
for (i=0; i<CalcInfo.nmo; i++) {
95
j = CalcInfo.pitz2ci[i];
96
CalcInfo.ci2pitz[j] = i;
100
/* Set up an array to map absolute ci order to relative Pitzer order */
101
CalcInfo.ci2relpitz = init_int_array(CalcInfo.nmo);
102
for (h=0,cnt=0; h<CalcInfo.nirreps; h++) {
103
for (i=0; i<CalcInfo.orbs_per_irr[h]; i++,cnt++) {
104
j = CalcInfo.pitz2ci[cnt];
105
CalcInfo.ci2relpitz[j] = i;
109
if (Params.print_lvl > 4) {
110
fprintf(outfile, "\nPitzer to CI order array = \n");
111
for (i=0; i<CalcInfo.nmo; i++) {
112
fprintf(outfile, "%3d ", CalcInfo.pitz2ci[i]);
114
fprintf(outfile, "\n");
118
CalcInfo.nbstri = (CalcInfo.nmo * (CalcInfo.nmo + 1)) / 2 ;
119
check((CalcInfo.nbstri <= IOFF_MAX),
120
"(get_mo_info): IOFF_MAX may not large enough!");
122
/* transform orbsym vector to new MO order */
123
CalcInfo.orbsym = init_int_array(CalcInfo.nmo);
125
for (i=0,cnt=0; i<CalcInfo.nirreps; i++) {
126
for (j=0; j<CalcInfo.orbs_per_irr[i]; j++,cnt++) {
127
k = CalcInfo.pitz2ci[cnt];
128
CalcInfo.orbsym[k] = i;
132
CalcInfo.num_fzv_orbs = 0; CalcInfo.num_vir_orbs = 0;
133
for (i=0; i<CalcInfo.nirreps; i++) {
134
CalcInfo.num_fzv_orbs += CalcInfo.frozen_uocc[i];
135
CalcInfo.num_vir_orbs += CalcInfo.rstr_uocc[i];
138
CalcInfo.npop = CalcInfo.nmo - CalcInfo.num_fzv_orbs -
139
CalcInfo.num_vir_orbs;
141
CalcInfo.num_fzc_orbs = 0;
142
CalcInfo.num_cor_orbs = 0;
143
for (i=0; i<CalcInfo.nirreps; i++) {
144
CalcInfo.num_fzc_orbs += CalcInfo.frozen_docc[i];
146
for (i=0; i<CalcInfo.nirreps; i++) {
147
CalcInfo.num_cor_orbs += CalcInfo.rstr_docc[i];
150
/* construct the CalcInfo.ras_orbs array (may not be of any use now) */
152
CalcInfo.fzc_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nmo);
153
CalcInfo.cor_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nmo);
154
CalcInfo.vir_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nmo);
155
CalcInfo.fzv_orbs = init_int_matrix(CalcInfo.nirreps,CalcInfo.nmo);
158
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
159
for (j=0; j<CalcInfo.frozen_docc[irrep]; j++)
160
CalcInfo.fzc_orbs[irrep][j] = cnt++;
163
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
164
for (j=0; j<CalcInfo.rstr_docc[irrep]; j++)
165
CalcInfo.cor_orbs[irrep][j] = cnt++;
168
CalcInfo.ras_orbs = (int ***) malloc (MAX_RAS_SPACES * sizeof(int **));
169
for (i=0; i<MAX_RAS_SPACES; i++) {
170
CalcInfo.ras_orbs[i] = init_int_matrix(CalcInfo.nirreps,
172
for (irrep=0; irrep<CalcInfo.nirreps; irrep++) {
173
for (j=0; j<CalcInfo.ras_opi[i][irrep]; j++) {
174
CalcInfo.ras_orbs[i][irrep][j] = cnt++;
180
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
181
for (j=0; j<CalcInfo.rstr_uocc[irrep]; j++)
182
CalcInfo.vir_orbs[irrep][j] = cnt++;
185
for (irrep=0; irrep<CalcInfo.nirreps; irrep++)
186
for (j=0; j<CalcInfo.frozen_uocc[irrep]; j++)
187
CalcInfo.fzv_orbs[irrep][j] = cnt++;
191
/* get the Pitzer arrays first, last, fstact, lstact, and active */
192
CalcInfo.first = init_int_array(CalcInfo.nirreps);
193
CalcInfo.last = init_int_array(CalcInfo.nirreps);
194
CalcInfo.fstact = init_int_array(CalcInfo.nirreps);
195
CalcInfo.lstact = init_int_array(CalcInfo.nirreps);
196
CalcInfo.active = init_int_array(CalcInfo.nirreps);
198
/* I think I never use this... --CDS 6/12/04
199
pitzer_arrays(CalcInfo.nirreps, CalcInfo.frozen_docc, CalcInfo.frozen_uocc,
200
CalcInfo.orbs_per_irr, CalcInfo.first, CalcInfo.last,
201
CalcInfo.fstact, CalcInfo.lstact, CalcInfo.active);
204
/* allocate memory to store the MO coefficient matrix symm blocked */
206
CalcInfo.mo_coeffs = (double ***) malloc(CalcInfo.nirreps *
208
for (irrep=0; irrep<CalcInfo.nirreps; irrep++) {
209
i = CalcInfo.orbs_per_irr[irrep];
211
CalcInfo.mo_coeffs[irrep] = block_matrix(i,i);
214
if (Params.print_lvl > 0) {
215
fprintf(outfile, "ORBITALS:");
216
fprintf(outfile, "\n FROZEN_DOCC = ");
217
for (i=0; i<CalcInfo.nirreps; i++) {
218
fprintf(outfile, "%2d ", CalcInfo.frozen_docc[i]);
220
fprintf(outfile, "\n RESTR_DOCC = ");
221
for (i=0; i<CalcInfo.nirreps; i++) {
222
fprintf(outfile, "%2d ", CalcInfo.rstr_docc[i]);
224
fprintf(outfile, "\n DOCC = ");
225
for (i=0; i<CalcInfo.nirreps; i++) {
226
fprintf(outfile, "%2d ", CalcInfo.docc[i]);
228
fprintf(outfile, "\n SOCC = ");
229
for (i=0; i<CalcInfo.nirreps; i++) {
230
fprintf(outfile, "%2d ", CalcInfo.socc[i]);
232
fprintf(outfile, "\n RESTR_UOCC = ");
233
for (i=0; i<CalcInfo.nirreps; i++) {
234
fprintf(outfile, "%2d ", CalcInfo.rstr_uocc[i]);
236
fprintf(outfile, "\n FROZEN_UOCC = ");
237
for (i=0; i<CalcInfo.nirreps; i++) {
238
fprintf(outfile, "%2d ", CalcInfo.frozen_uocc[i]);
241
for (i=0; i<MAX_RAS_SPACES; i++) {
242
fprintf(outfile, "\n RAS %d = ",i+1);
243
for (j=0; j<CalcInfo.nirreps; j++) {
244
fprintf(outfile, "%2d ", CalcInfo.ras_opi[i][j]);
247
fprintf(outfile, "\n");
249
fprintf(outfile, " MOL ORBS = %6d\n", CalcInfo.nmo);
250
fprintf(outfile, " FROZEN CORE = %6d RESTR CORE = %6d\n",
251
CalcInfo.num_fzc_orbs, CalcInfo.num_cor_orbs);
252
fprintf(outfile, "\n");
261
** Form the first/last/active arrays for the orbitals in Pitzer order
262
** Based on code taken from TRANSQT
267
void pitzer_arrays(int nirreps, int *frdocc, int *fruocc, int *orbspi,
268
int *first, int *last, int *fstact, int *lstact,int *active)
272
int first_offset, last_offset;
275
* Construct first and last index arrays: this defines the first
276
* absolute orbital index (Pitzer ordering) and last absolute orbital
277
* index for each irrep. When there are no orbitals for an irrep, the
278
* value is -1 for first[] and -2 for last[]. Note that there must be
279
* orbitals in the first irrep (i.e. totally symmetric) for this to work.
281
for(h=0; h < nirreps; h++) {
287
last_offset = orbspi[0] - 1;
288
first[0] = first_offset;
289
last[0] = last_offset;
291
for(h=1; h < nirreps; h++) {
292
first_offset += orbspi[h-1];
293
last_offset += orbspi[h];
295
first[h] = first_offset;
296
last[h] = last_offset;
301
* Construct first and last active index arrays: this defines the first
302
* absolute orbital index (Pitzer ordering) and last absolute orbital
303
* index for each irrep, excluding frozen orbitals. When there are no
304
* orbitals for an irrep, the value is -1 for first[] and -2 for last[].
305
* Note that there must be orbitals in the first irrep (i.e. totally
306
* symmetric) for this to work.
308
for(h=0; h < nirreps; h++) {
313
first_offset = frdocc[0];
314
last_offset = orbspi[0] - fruocc[0] - 1;
315
fstact[0] = first_offset;
316
lstact[0] = last_offset;
318
for(h=1; h < nirreps; h++) {
319
first_offset += orbspi[h-1]+frdocc[h]-frdocc[h-1];
320
last_offset += orbspi[h] - fruocc[h] + fruocc[h-1];
322
fstact[h] = first_offset;
323
lstact[h] = last_offset;
327
/* Now define active[] such that frozen orbitals are taken into account */
328
for(h=0; h < nirreps; h++) {
329
active[h] = orbspi[h]-frdocc[h]-fruocc[h];
339
** Read in the molecular orbital matrix from PSIF_CHKPT and put them in CalcInfo
341
void read_cur_orbs(void)
343
int i, j, h, dim, nirreps;
346
nirreps = CalcInfo.nirreps;
348
chkpt_init(PSIO_OPEN_OLD);
349
for (h=0; h<nirreps; h++) {
350
dim = CalcInfo.orbs_per_irr[h];
351
if (dim==0) continue;
352
tmat = chkpt_rd_scf_irrep(h);
353
for (i=0; i<dim; i++)
354
for (j=0; j<dim; j++)
355
CalcInfo.mo_coeffs[h][i][j] = tmat[i][j];
367
** This function copies SCF eigenvectors from moinfo.scf_vector into
368
** an array of matrices. Columns corresponding to inactive orbitals
369
** may be deleted if desired.
371
** Code taken from TRANSQT
375
double *** construct_evects(int nirreps, int *active, int *orbspi,
376
int *first, int *last, int *fstact, int *lstact,
380
int h, row, col, p, q;
383
evects = (double ***) malloc(nirreps * sizeof(double **));
385
for (h=0; h<nirreps; h++) {
387
evects[h] = block_matrix(orbspi[h],active[h]);
389
for(p=first[h]; p <= last[h]; p++) {
391
for(q=fstact[h]; q <= lstact[h]; q++) {
393
evects[h][row][col] = CalcInfo.mo_matrix[p][q];
398
fprintf(outfile,"\n\tMolecular Orbitals for Irrep %s\n",
400
print_mat(evects[h],orbspi[h],active[h],outfile);
413
void destruct_evects(int nirreps, double ***evects)
417
for (h=0; h<nirreps; h++)
418
if (evects[h] != NULL) free_block(evects[h]);
422
}} // end namespace psi::detcas
added
removed
src/bin/detcas/get_mo_info.cc
