← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/oeprop/main.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/oeprop/main.c
 
1
#include "includes.h"
 
2
#include "prototypes.h"
 
3
#include "globals.h"
 
4
 
 
5
 
 
6
int main(int argc, char* argv[]) {
 
7
 
 
8
 int i,j,k,l,count;
 
9
 char buffer[80];       /* buffer string */
 
10
 
 
11
        /* Setting defaults */
 
12
 read_opdm = 0;
 
13
 opdm_file = 76;
 
14
 asymm_opdm = 0;
 
15
 wrtnos = 0;
 
16
 spin_prop = 0;
 
17
 print_lvl = 1;
 
18
 wrt_dipmom = 1;
 
19
 corr = 0;
 
20
 mpmax = 1;
 
21
 mp_ref = 0;
 
22
 nuc_esp = 1;
 
23
 grid = 0;
 
24
 grid3d = 0;
 
25
 nix = 10;
 
26
 niy = 10;
 
27
 niz = 10;
 
28
 grid_zmin = 0.0;
 
29
 grid_zmax = 3.0;
 
30
 edgrad_logscale = 5;
 
31
 zvec_file = 86;
 
32
 delete_zvec = 1;
 
33
 
 
34
 
 
35
        /* Initialization and printing intro */
 
36
 start_io(argc,argv);
 
37
 print_intro();
 
38
 
 
39
/*************************** Main Code *******************************/
 
40
 
 
41
        /* Reading in basic information from checkpoint file */
 
42
 
 
43
 title = chkpt_rd_label();
 
44
 natom = chkpt_rd_natom();
 
45
 natom3 = natom * 3;
 
46
 nmo = chkpt_rd_nmo();
 
47
 nbfso = chkpt_rd_nso();
 
48
 nbfao = chkpt_rd_nao();
 
49
 natri = nbfao * (nbfao+1)/2;
 
50
 nstri = nbfso * (nbfso+1)/2;
 
51
 nshell = chkpt_rd_nshell();
 
52
 nprim = chkpt_rd_nprim();
 
53
 iopen = chkpt_rd_iopen();
 
54
 nirreps = chkpt_rd_nirreps(); 
 
55
 nsym = chkpt_rd_nsymhf();
 
56
 orbspi = chkpt_rd_orbspi();
 
57
 sopi = chkpt_rd_sopi();
 
58
 clsdpi = chkpt_rd_clsdpi();    
 
59
 openpi = chkpt_rd_openpi();
 
60
 irr_labs = chkpt_rd_irr_labs();
 
61
 geom = chkpt_rd_geom();
 
62
 zvals = chkpt_rd_zvals();
 
63
 usotao = chkpt_rd_usotao();
 
64
    
 
65
 /* Parsing */
 
66
 parsing();
 
67
 
 
68
 if (strcmp(ref,"UHF") != 0) {
 
69
   scf_evec_so = chkpt_rd_scf();
 
70
   scf_evals = chkpt_rd_evals();
 
71
   scf_evec_ao = block_matrix(nbfao,nmo);
 
72
   mmult(usotao,1,scf_evec_so,0,scf_evec_ao,0,nbfao,nbfso,nmo,0);
 
73
 }
 
74
 
 
75
 /* parse the grid information */
 
76
 if (ip_exist("GRID",0))
 
77
   grid_parse();
 
78
 
 
79
 /* Computing total charge of the system */
 
80
 
 
81
 charge = 0;
 
82
 for(i=0;i<nirreps;i++)
 
83
   charge -= 2*clsdpi[i] + openpi[i];
 
84
 for(i=0;i<natom;i++)
 
85
   charge += zvals[i];
 
86
 
 
87
 /* Setting up an offset array */
 
88
 
 
89
 ioff = init_int_array(nbfao);
 
90
 ioff[0] = 0;
 
91
 for(i=1;i<nbfao;i++) {
 
92
   ioff[i] = ioff[i-1] + i;
 
93
 }
 
94
 
 
95
        /* Computing double factorials df[i] = (i-1)!! */
 
96
 df[0] = 1.0;
 
97
 df[1] = 1.0;
 
98
 df[2] = 1.0;
 
99
 for(i=3; i<MAXFACT*2; i++) {
 
100
   df[i] = (i-1)*df[i-2];
 
101
 }
 
102
                 
 
103
 /* Printing tasks and parameters */
 
104
 
 
105
 if (print_lvl >= PRINTTASKPARAMLEVEL) {
 
106
   print_tasks();
 
107
   print_params();
 
108
 }
 
109
 
 
110
 /* Reading in basis set inforamtion */
 
111
 
 
112
 read_basset_info();
 
113
 init_xyz(); 
 
114
 
 
115
 /* Obtain a density matrix */
 
116
 
 
117
 if (read_opdm)
 
118
   read_density();
 
119
 else
 
120
   compute_density();
 
121
 
 
122
/*if (strcmp(ref,"UHF") != 0) { */
 
123
  compute_overlap();
 
124
   
 
125
  /* Obtain natural orbitals */
 
126
  if (read_opdm && wrtnos) 
 
127
    get_nmo(); 
 
128
   
 
129
  chkpt_close();
 
130
   
 
131
  /* Mulliken population analysis */
 
132
  print_pop_header();
 
133
  populate();
 
134
 
 
135
  /* Compute coordinates of the MP reference point if needed */
 
136
  if (mp_ref != -1)
 
137
    compute_mp_ref_xyz();
 
138
 
 
139
  /* Moving molecule to the reference center!
 
140
   Attention, ever since coordinates of atoms and of the grid box 
 
141
   are stored in this new coordinate system. All coordinates are 
 
142
   transformed back at the moment of printing out. */
 
143
 
 
144
  for(i=0;i<natom;i++) {
 
145
    geom[i][0] -= mp_ref_xyz[0];
 
146
    geom[i][1] -= mp_ref_xyz[1];
 
147
    geom[i][2] -= mp_ref_xyz[2];
 
148
    Lm_ref_xyz[0] -= mp_ref_xyz[0];
 
149
    Lm_ref_xyz[1] -= mp_ref_xyz[1];
 
150
    Lm_ref_xyz[2] -= mp_ref_xyz[2];
 
151
  }
 
152
  if (grid) {
 
153
    grid_origin[0] -= mp_ref_xyz[0];
 
154
    grid_origin[1] -= mp_ref_xyz[1];
 
155
    grid_origin[2] -= mp_ref_xyz[2];
 
156
  }
 
157
 
 
158
        /* Computing one-electron integrals in 
 
159
           terms of Cartesian Gaussians, 
 
160
           electric first (dipole), second and third 
 
161
           moments W.R.T. origin, electrostatic potential,
 
162
           electric field and field gradients,
 
163
           electron and spin density at atomic centers,
 
164
           and various properties over a grid, if neccessary. */
 
165
 
 
166
 compute_oeprops();
 
167
 if (grid)
 
168
     compute_grid();
 
169
 
 
170
 print_mp();
 
171
 print_lm();
 
172
 if (nuc_esp)
 
173
   print_esp();
 
174
 if (grid) {
 
175
   grid_origin[0] += mp_ref_xyz[0];
 
176
   grid_origin[1] += mp_ref_xyz[1];
 
177
   grid_origin[2] += mp_ref_xyz[2];
 
178
   print_grid();
 
179
 }
 
180
 print_misc();
 
181
 
 
182
        /* Cleaning up */
 
183
 
 
184
 free_block(usotao);
 
185
 if (strcmp(ref,"UHF") != 0) {
 
186
   free_block(scf_evec_so);
 
187
   free_block(scf_evec_ao);
 
188
 }
 
189
 free(ioff);
 
190
 free(Ptot);
 
191
 free(phi);
 
192
 free(ex); free(ey); free(ez);
 
193
 free(dexx); free(deyy); free(dezz);
 
194
 free(dexy); free(dexz); free(deyz);
 
195
 if (spin_prop) {
 
196
   free(Pspin);
 
197
   free(ahfsxx);
 
198
   free(ahfsyy);
 
199
   free(ahfszz);
 
200
   free(ahfsxy);
 
201
   free(ahfsxz);
 
202
   free(ahfsyz);
 
203
 }
 
204
 free(S);
 
205
/*}  end non-UHF routines */
 
206
 stop_io();
 
207
 exit(PSI_RETURN_SUCCESS);
 
208
 
 
209
}
 
210
 
 
211
 
 
212
char *gprgid()
 
213
{
 
214
 char *prgid = "OEPROP";
 
215
   
 
216
 return(prgid);
 
217
}