← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/detci/get_mo_info.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/detci/get_mo_info.cc
 
1
/*! \file
 
2
    \ingroup DETCI
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <libipv1/ip_lib.h>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libchkpt/chkpt.h>
 
10
#include <libqt/qt.h>
 
11
#include "structs.h"
 
12
#define EXTERN
 
13
#include "globals.h"
 
14
 
 
15
namespace psi { namespace detci {
 
16
 
 
17
/*
 
18
** GET_MO_INFO
 
19
** 
 
20
** Reads PSIF_CHKPT & input.dat and gets all sorts of useful information about 
 
21
** the molecular orbitals (such as their reordering array, the docc
 
22
** array, frozen orbitals, etc.)
 
23
**
 
24
** Created by C. David Sherrill on 17 November 1994
 
25
**
 
26
** Updated
 
27
** CDS  1/18/95 to read SCF eigenvalues also (for MP2 guess vector)
 
28
** CDS  1/ 5/97 to use nifty new ras_set() function (which transqt has been
 
29
**              using for some time).
 
30
**
 
31
*/
 
32
void get_mo_info(void)
 
33
{
 
34
   int i, j, k, tmp, cnt, irrep, errcod, errbad;
 
35
   int size;
 
36
   double *eig_unsrt;
 
37
   int parsed_ras1=0, parsed_ras2=0, do_ras4;
 
38
   int *rstr_docc, *rstr_uocc;
 
39
 
 
40
   CalcInfo.maxKlist = 0.0; 
 
41
 
 
42
   chkpt_init(PSIO_OPEN_OLD);
 
43
   CalcInfo.nirreps = chkpt_rd_nirreps();
 
44
   CalcInfo.nso = chkpt_rd_nmo();
 
45
   CalcInfo.nmo = chkpt_rd_nmo();
 
46
   CalcInfo.iopen = chkpt_rd_iopen();
 
47
   CalcInfo.labels = chkpt_rd_irr_labs();
 
48
   CalcInfo.orbs_per_irr = chkpt_rd_orbspi();
 
49
   CalcInfo.so_per_irr = chkpt_rd_sopi();
 
50
   CalcInfo.docc = chkpt_rd_clsdpi();
 
51
   CalcInfo.socc = chkpt_rd_openpi();
 
52
   CalcInfo.enuc = chkpt_rd_enuc();
 
53
   CalcInfo.escf = chkpt_rd_escf();
 
54
   CalcInfo.efzc = chkpt_rd_efzc();
 
55
   eig_unsrt = chkpt_rd_evals(); 
 
56
   chkpt_close(); 
 
57
 
 
58
   if (CalcInfo.iopen && Parameters.opentype == PARM_OPENTYPE_NONE) {
 
59
      fprintf(outfile, "Warning: iopen=1,opentype=none. Making iopen=0\n");
 
60
      CalcInfo.iopen = 0;
 
61
      }
 
62
   else if (!CalcInfo.iopen && (Parameters.opentype == PARM_OPENTYPE_HIGHSPIN
 
63
      || Parameters.opentype == PARM_OPENTYPE_SINGLET)) {
 
64
      fprintf(outfile,"Warning: iopen=0,opentype!=closed. Making iopen=1\n");
 
65
      CalcInfo.iopen = 1;
 
66
      }
 
67
   if (Parameters.ref_sym >= CalcInfo.nirreps) {
 
68
      fprintf(outfile,"Warning: ref_sym >= nirreps.  Setting ref_sym=0\n");
 
69
      Parameters.ref_sym = 0;
 
70
      }
 
71
 
 
72
   CalcInfo.frozen_docc = init_int_array(CalcInfo.nirreps);
 
73
   CalcInfo.frozen_uocc = init_int_array(CalcInfo.nirreps);
 
74
   rstr_docc = init_int_array(CalcInfo.nirreps);
 
75
   rstr_uocc = init_int_array(CalcInfo.nirreps);
 
76
   CalcInfo.explicit_core = init_int_array(CalcInfo.nirreps);
 
77
   CalcInfo.explicit_vir  = init_int_array(CalcInfo.nirreps);
 
78
   CalcInfo.reorder = init_int_array(CalcInfo.nmo);
 
79
   CalcInfo.ras_opi = init_int_matrix(4,CalcInfo.nirreps);
 
80
 
 
81
   if (!ras_set2(CalcInfo.nirreps, CalcInfo.nmo, 1, (Parameters.fzc) ?  1:0,
 
82
                CalcInfo.orbs_per_irr, CalcInfo.docc, CalcInfo.socc,
 
83
                CalcInfo.frozen_docc, CalcInfo.frozen_uocc,
 
84
                rstr_docc, rstr_uocc,
 
85
                CalcInfo.ras_opi, CalcInfo.reorder, 1, 0))
 
86
   {
 
87
     fprintf(outfile, "Error in ras_set().  Aborting.\n");
 
88
     exit(1);
 
89
   }
 
90
 
 
91
   if (1) { /* for now, always treat restricted as frozen */
 
92
     for (i=0; i<CalcInfo.nirreps; i++) {
 
93
       CalcInfo.frozen_docc[i] += rstr_docc[i];
 
94
       CalcInfo.frozen_uocc[i] += rstr_uocc[i];
 
95
     } 
 
96
   }
 
97
   else { /* for future use */
 
98
     for (i=0; i<CalcInfo.nirreps; i++) {
 
99
       CalcInfo.explicit_core[i] = rstr_docc[i];
 
100
       CalcInfo.explicit_vir[i]  = rstr_uocc[i];
 
101
     }
 
102
   }
 
103
 
 
104
   free(rstr_docc);  free(rstr_uocc);
 
105
 
 
106
 
 
107
   /* Compute maximum number of orbitals per irrep including
 
108
   ** and not including fzv
 
109
   */
 
110
  CalcInfo.max_orbs_per_irrep = 0;
 
111
  CalcInfo.max_pop_per_irrep = 0;
 
112
  for (i=0; i<CalcInfo.nirreps; i++) {
 
113
     if (CalcInfo.max_orbs_per_irrep < CalcInfo.orbs_per_irr[i])
 
114
       CalcInfo.max_orbs_per_irrep = CalcInfo.orbs_per_irr[i];
 
115
     if (CalcInfo.max_pop_per_irrep < (CalcInfo.orbs_per_irr[i] - 
 
116
                                   CalcInfo.frozen_uocc[i]))
 
117
       CalcInfo.max_pop_per_irrep = CalcInfo.orbs_per_irr[i] -
 
118
                                    CalcInfo.frozen_uocc[i];      
 
119
     }
 
120
 
 
121
 
 
122
   /* construct the "ordering" array, which maps the other direction */
 
123
   /* i.e. from a CI orbital to a Pitzer orbital                     */
 
124
   CalcInfo.order = init_int_array(CalcInfo.nmo);
 
125
   for (i=0; i<CalcInfo.nmo; i++) {
 
126
      j = CalcInfo.reorder[i];
 
127
      CalcInfo.order[j] = i;
 
128
      }
 
129
 
 
130
 
 
131
   if (Parameters.print_lvl > 4) {
 
132
      fprintf(outfile, "\nReordering array = \n");
 
133
      for (i=0; i<CalcInfo.nmo; i++) {
 
134
         fprintf(outfile, "%3d ", CalcInfo.reorder[i]);
 
135
         }
 
136
      fprintf(outfile, "\n");
 
137
      }
 
138
 
 
139
   CalcInfo.nmotri = (CalcInfo.nmo * (CalcInfo.nmo + 1)) / 2 ;
 
140
 
 
141
   /* transform orbsym vector to new MO order */
 
142
   CalcInfo.orbsym = init_int_array(CalcInfo.nmo);
 
143
   CalcInfo.scfeigval = init_array(CalcInfo.nmo);
 
144
   if(Parameters.zaptn) {
 
145
     CalcInfo.scfeigvala = init_array(CalcInfo.nmo);
 
146
     CalcInfo.scfeigvalb = init_array(CalcInfo.nmo);
 
147
     }
 
148
 
 
149
   for (i=0,cnt=0; i<CalcInfo.nirreps; i++) {
 
150
      for (j=0; j<CalcInfo.orbs_per_irr[i]; j++,cnt++) {
 
151
         k = CalcInfo.reorder[cnt];
 
152
         CalcInfo.orbsym[k] = i;
 
153
         }
 
154
      }
 
155
 
 
156
   for (i=0; i<CalcInfo.nmo; i++) {
 
157
      j = CalcInfo.reorder[i];
 
158
      CalcInfo.scfeigval[j] = eig_unsrt[i];
 
159
      if(Parameters.zaptn) {
 
160
        CalcInfo.scfeigvala[j] = eig_unsrt[i];
 
161
        CalcInfo.scfeigvalb[j] = eig_unsrt[i];
 
162
        }
 
163
      }
 
164
   free(eig_unsrt);
 
165
 
 
166
   /* calculate number of electrons */
 
167
   CalcInfo.num_alp = CalcInfo.num_bet = CalcInfo.spab = 0;
 
168
   if (Parameters.opentype == PARM_OPENTYPE_NONE ||
 
169
       Parameters.opentype == PARM_OPENTYPE_HIGHSPIN) {
 
170
      for (i=0; i<CalcInfo.nirreps; i++) {
 
171
         CalcInfo.num_alp += CalcInfo.docc[i] + CalcInfo.socc[i];
 
172
         CalcInfo.num_bet += CalcInfo.docc[i];
 
173
         }
 
174
      }
 
175
   else if (Parameters.opentype == PARM_OPENTYPE_SINGLET) {
 
176
      for (i=0; i<CalcInfo.nirreps; i++) { /* closed-shell part */
 
177
         CalcInfo.spab += CalcInfo.socc[i];      
 
178
         CalcInfo.num_alp += CalcInfo.docc[i];
 
179
         CalcInfo.num_bet += CalcInfo.docc[i];
 
180
         }
 
181
      if (CalcInfo.spab % 2) { 
 
182
         fprintf(outfile,"For opentype=singlet must have even number ");
 
183
         fprintf(outfile,"of socc electrons!\n");
 
184
         exit(1);
 
185
         }
 
186
      CalcInfo.spab /= 2;
 
187
      tmp = 0;
 
188
      for (i=0; i<CalcInfo.nirreps; i++) {
 
189
         j = CalcInfo.socc[i];
 
190
         k = 0;
 
191
         while (k < j) {
 
192
            if (tmp < CalcInfo.spab) {
 
193
               CalcInfo.num_alp++;
 
194
               tmp++; 
 
195
               k++;
 
196
               }
 
197
            else { 
 
198
               CalcInfo.num_bet++;
 
199
               tmp++;
 
200
               k++;
 
201
               } 
 
202
            }
 
203
         }
 
204
      }
 
205
   else {
 
206
      fprintf(outfile, "(get_mo_info): Can't handle opentype = %d\n",
 
207
         Parameters.opentype);
 
208
      exit(1);
 
209
      }
 
210
 
 
211
   /* at this stage I've already overwritten CalcInfo.frozen_docc,
 
212
      CalcInfo.rstr_docc, etc, to be their internal DETCI meaning
 
213
      and not necessarily the user input.  Internally frozen_docc
 
214
      and frozen_uocc refer to dropped core/virt that are never
 
215
      allowed to change occupation
 
216
   */
 
217
   CalcInfo.num_fzv_orbs = 0;
 
218
   CalcInfo.num_vir_orbs = 0;
 
219
   CalcInfo.num_fzc_orbs = 0;
 
220
   CalcInfo.num_cor_orbs = 0;
 
221
   for (i=0; i<CalcInfo.nirreps; i++) {
 
222
      CalcInfo.num_fzv_orbs += CalcInfo.frozen_uocc[i];  
 
223
      CalcInfo.num_vir_orbs += CalcInfo.explicit_vir[i];
 
224
      CalcInfo.num_fzc_orbs += CalcInfo.frozen_docc[i];  
 
225
      CalcInfo.num_cor_orbs += CalcInfo.explicit_core[i];
 
226
   }
 
227
 
 
228
   /* calculate number of orbitals active in CI */
 
229
   /* maybe this changes later for cor orbs, depends on where we go w/ it */
 
230
   CalcInfo.num_ci_orbs = CalcInfo.nmo - CalcInfo.num_fzc_orbs - 
 
231
     CalcInfo.num_fzv_orbs;
 
232
 
 
233
   if ((CalcInfo.num_ci_orbs * (CalcInfo.num_ci_orbs + 1)) / 2 > IOFF_MAX) {
 
234
      fprintf(outfile, "Error: IOFF_MAX not large enough!\n");
 
235
      exit(1);
 
236
   }
 
237
 
 
238
   CalcInfo.num_alp_expl = CalcInfo.num_alp - CalcInfo.num_fzc_orbs;
 
239
   CalcInfo.num_bet_expl = CalcInfo.num_bet - CalcInfo.num_fzc_orbs;
 
240
 
 
241
   /* construct the CalcInfo.ras_orbs array (may not be of any use now) */
 
242
   cnt = 0;
 
243
   for (i=0; i<4; i++) {
 
244
     CalcInfo.ras_orbs[i] = init_int_matrix(CalcInfo.nirreps,
 
245
       CalcInfo.num_ci_orbs);
 
246
     for (irrep=0; irrep<CalcInfo.nirreps; irrep++) {
 
247
       for (j=0; j<CalcInfo.ras_opi[i][irrep]; j++) {
 
248
         CalcInfo.ras_orbs[i][irrep][j] = cnt++;
 
249
       }
 
250
     }
 
251
   }
 
252
}
 
253
 
 
254
}} // namespace psi::detci
 
255