← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/mocube/mocube.c

  • 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/mocube/mocube.c
1
 
/* MOCUBE - make Gaussian-compatible cube file */
2
 
 
3
 
#include <stdio.h>
4
 
#include <stdlib.h>
5
 
#include <string.h>
6
 
#include <libipv1/ip_lib.h>
7
 
#include <libchkpt/chkpt.h>
8
 
#include <libciomr/libciomr.h>
9
 
#include <libqt/qt.h>
10
 
#include "mocube.h"
11
 
#define DEFAULT_BORDER (4.0)
12
 
#define DEFAULT_NGRID (80)
13
 
#define MAX(I,J) ((I>J) ? I : J)
14
 
#define MIN(I,J) ((I<J) ? I : J)
15
 
 
16
 
char *progid;
17
 
void init_io(int argc, char *argv[]);
18
 
void get_params(void);
19
 
void setup_delta(double **scf, double **u);
20
 
void compute_mos(double *movals, double x, double y, double z,
21
 
    double **scf, double **u);
22
 
void wrt_cube(void);
23
 
void exit_io(void);
24
 
FILE *infile, *outfile, *cubfile;
25
 
char *psi_file_prefix;
26
 
 
27
 
int main(int argc, char *argv[]) {
28
 
  int i,j, nirreps, cnt_x, cnt_y, cnt_z;
29
 
  double x, y, z, *movals, **scf, **u;
30
 
  char *filename;
31
 
 
32
 
  init_io(argc, argv);
33
 
  fprintf(outfile,"\n\t*****************************************************\n");
34
 
  fprintf(outfile,  "\t*  MOCUBE: Produces a Gaussian compatible .cub file *\n");
35
 
  fprintf(outfile,  "\t*****************************************************\n");
36
 
 
37
 
  get_params();
38
 
  strcpy(cube.title,"Title Card Provided!");
39
 
  strcpy(cube.subtitle,"Alpha MO grid values generated by PSI3's mocube");
40
 
 
41
 
  /* read zvals and geometry from chkpt */
42
 
  chkpt_init(PSIO_OPEN_OLD);
43
 
  cube.zvals = chkpt_rd_zvals();
44
 
  cube.geom = chkpt_rd_geom();
45
 
  cube.natom = params.natom;
46
 
  chkpt_close();
47
 
 
48
 
  /* print_mat(cube.geom, cube.natom, 3, stdout); */
49
 
 
50
 
  /* determine output filename */
51
 
  if (ip_exist("CUBFILE",0))
52
 
    ip_string("CUBFILE", &filename, 0);
53
 
  else {
54
 
    filename = (char *) malloc(20*sizeof(char *));
55
 
    sprintf(filename , "mos.cub");
56
 
  }
57
 
  ffile(&cubfile, filename, 0);
58
 
  fprintf(outfile,"Writing to cubfile %s\n", filename);
59
 
 
60
 
  /*** determine location of grid points ***/
61
 
    /* determine border size */
62
 
  cube.border = DEFAULT_BORDER; 
63
 
  if (ip_exist("BORDER",0)) ip_data("BORDER","%lf",&(cube.border),0);
64
 
 
65
 
    /* determine start/end of grid */
66
 
  for (i=0; i<3; ++i) {
67
 
    cube.grid_start[i] = 0.0;
68
 
    cube.grid_end[i] = 0.0;
69
 
  }
70
 
  for (i=0; i<cube.natom; ++i) {
71
 
    for (j=0; j<3; ++j) {
72
 
      cube.grid_start[j] = MIN(cube.grid_start[j], cube.geom[i][j]);
73
 
      cube.grid_end[j] = MAX(cube.grid_end[j], cube.geom[i][j]);
74
 
    }
75
 
  }
76
 
    /* include molecule plus border in box */
77
 
  for (i=0; i<3; ++i) {
78
 
    cube.grid_start[i] -= cube.border;
79
 
    cube.grid_end[i] += cube.border;
80
 
  }
81
 
  fprintf(outfile,"Grid starts at: %10.5lf %10.5lf %10.5lf\n",
82
 
      cube.grid_start[0], cube.grid_start[1], cube.grid_start[2]);
83
 
  fprintf(outfile,"Grid ends at  : %10.5lf %10.5lf %10.5lf\n",
84
 
      cube.grid_end[0], cube.grid_end[1], cube.grid_end[2]);
85
 
 
86
 
    /* determine number of grid points */
87
 
  if (ip_exist("NGRID",0)) {
88
 
    ip_data("NGRID", "%d", &j, 0);
89
 
    for (i=0;i<3;++i) cube.ngrid[i] = j;
90
 
  }
91
 
  else {
92
 
    for (i=0;i<3;++i) cube.ngrid[i] = DEFAULT_NGRID;
93
 
  }
94
 
  fprintf(outfile,"Number of grid points: %d %d %d\n",
95
 
      cube.ngrid[0], cube.ngrid[1], cube.ngrid[2]);
96
 
 
97
 
    /* calculate step_size */
98
 
  for (i=0;i<3;++i)
99
 
    cube.step_size[i] = (cube.grid_end[i]-cube.grid_start[i])
100
 
      /( cube.ngrid[i] - 1 );
101
 
  fprintf(outfile,"Step sizes between grid points: %10.5lf %10.5lf %10.5lf\n",
102
 
      cube.step_size[0], cube.step_size[1], cube.step_size[2]);
103
 
 
104
 
  fflush(outfile);
105
 
 
106
 
  /* prepare for numerous compute_mos() calls */
107
 
  chkpt_init(PSIO_OPEN_OLD);
108
 
  u = chkpt_rd_usotao();
109
 
  scf = block_matrix(params.nmo, params.nmo);
110
 
  chkpt_close();
111
 
  setup_delta(scf, u); /* also determines HOMO,LUMO */
112
 
 
113
 
  /* determine which MOS to plot */
114
 
  if (ip_exist("CUBEMOS",0)) {
115
 
    ip_count("CUBEMOS", &(cube.nmo_to_plot), 0);
116
 
    cube.mos_to_plot = (int *) malloc( cube.nmo_to_plot * sizeof(int) );
117
 
 
118
 
    for (i=0; i<cube.nmo_to_plot; ++i) {
119
 
      ip_data("CUBEMOS","%d", &(cube.mos_to_plot[i]), 1, i);
120
 
      cube.mos_to_plot[i] -= 1;
121
 
    }
122
 
  }
123
 
  else { /* do HOMO and LUMO */
124
 
    cube.nmo_to_plot = 2;
125
 
    cube.mos_to_plot = (int *) malloc( cube.nmo_to_plot * sizeof(int) );
126
 
    cube.mos_to_plot[0] = params.homo;
127
 
    cube.mos_to_plot[1] = params.lumo;
128
 
  }
129
 
  fprintf(outfile,"MOs to plot: ");
130
 
  for (i=0; i<cube.nmo_to_plot; ++i)
131
 
    fprintf(outfile,"%d ", cube.mos_to_plot[i] + 1);
132
 
  fprintf(outfile,"\n");
133
 
 
134
 
  /* malloc memory for grid */
135
 
  cube.grid = (double ****) malloc(cube.ngrid[0] * sizeof(double ***));
136
 
  for (cnt_x=0; cnt_x<cube.ngrid[0]; ++cnt_x) {
137
 
    cube.grid[cnt_x] = (double ***) malloc(cube.ngrid[1] * sizeof(double **));
138
 
    for (cnt_y=0; cnt_y<cube.ngrid[1]; ++cnt_y)
139
 
      cube.grid[cnt_x][cnt_y] = block_matrix(cube.ngrid[2], cube.nmo_to_plot);
140
 
  }
141
 
 
142
 
  /* calculate MO values at grid points */
143
 
  movals = init_array(cube.nmo_to_plot);
144
 
  x = cube.grid_start[0];
145
 
  for (cnt_x=0; cnt_x<cube.ngrid[0]; ++cnt_x, x += cube.step_size[0]) {
146
 
    y = cube.grid_start[1];
147
 
    for (cnt_y=0; cnt_y<cube.ngrid[1]; ++cnt_y, y += cube.step_size[1]) {
148
 
      z = cube.grid_start[2];
149
 
      for (cnt_z=0; cnt_z<cube.ngrid[2]; ++cnt_z, z += cube.step_size[2]) {
150
 
        compute_mos(movals, x, y, z, scf, u);
151
 
        for (i=0; i<cube.nmo_to_plot; ++i)
152
 
          cube.grid[cnt_x][cnt_y][cnt_z][i] = movals[i];
153
 
      }
154
 
    }
155
 
  }
156
 
  free(movals);
157
 
 
158
 
  wrt_cube();
159
 
  free_block(u);
160
 
  free_block(scf);
161
 
  exit_io();
162
 
  return 0 ;
163
 
}
164
 
 
165
 
char *gprgid() {
166
 
  char *prgid = "MOCUBE";
167
 
  return (prgid);
168
 
}
169
 
 
170
 
void init_io(int argc, char *argv[]) {
171
 
  extern char *gprgid();
172
 
  char *filename;
173
 
  progid = (char *) malloc(strlen(gprgid())+2);
174
 
  sprintf(progid, ":%s",gprgid());
175
 
 
176
 
  psi_start(argc-1,argv+1,0);
177
 
  ip_cwk_add(":INPUT");
178
 
  ip_cwk_add(progid);
179
 
  psio_init();
180
 
  tstart(outfile);
181
 
}
182
 
 
183
 
void exit_io(void) {
184
 
  int i;
185
 
/*  psio_close(32,1); */
186
 
  psio_done();
187
 
  tstop(outfile);
188
 
  ip_done();
189
 
}
190
 
 
191
 
void get_params(void) {
192
 
  int i, irrep_max, index_max;
193
 
 
194
 
  chkpt_init(PSIO_OPEN_OLD);
195
 
  params.nao = chkpt_rd_nao();
196
 
  params.natom = chkpt_rd_natom();
197
 
  params.nirreps = chkpt_rd_nirreps();
198
 
  params.nirreps_present = chkpt_rd_nsymhf();
199
 
  params.nprim = chkpt_rd_nprim();
200
 
  params.nso = chkpt_rd_nso();
201
 
  params.nmo = chkpt_rd_nmo();
202
 
  params.clsdpi = chkpt_rd_clsdpi(); 
203
 
  params.label = chkpt_rd_label();
204
 
  params.openpi = chkpt_rd_openpi();
205
 
  params.orbspi = chkpt_rd_orbspi();
206
 
 
207
 
  for (i=0; i<params.nirreps; ++i)
208
 
    params.nclsd += params.clsdpi[i];
209
 
 
210
 
  chkpt_close();
211
 
}
212
 
 
213
 
/* this function writes a cubfile in Gaussian .cub format */
214
 
void wrt_cube(void) {
215
 
  int i,j,k,cnt,x,y,z;
216
 
 
217
 
  fprintf(cubfile," %s\n",cube.title);
218
 
  fprintf(cubfile," %s\n",cube.subtitle);
219
 
  fprintf(cubfile,"%5d%12.6lf%12.6lf%12.6lf\n",-1*cube.natom,
220
 
      cube.grid_start[0], cube.grid_start[1], cube.grid_start[2]);
221
 
  fprintf(cubfile,"%5d%12.6lf%12.6lf%12.6lf\n", cube.ngrid[0],
222
 
      cube.step_size[0], 0.0, 0.0);
223
 
  fprintf(cubfile,"%5d%12.6lf%12.6lf%12.6lf\n", cube.ngrid[1],
224
 
      0.0, cube.step_size[1], 0.0);
225
 
  fprintf(cubfile,"%5d%12.6lf%12.6lf%12.6lf\n", cube.ngrid[2],
226
 
      0.0, 0.0, cube.step_size[2]);
227
 
 
228
 
  for (i=0; i<params.natom; ++i) {
229
 
    fprintf(cubfile,"%5d%12.6lf%12.6lf%12.6lf%12.6lf\n", (int) cube.zvals[i],
230
 
      cube.zvals[i], cube.geom[i][0], cube.geom[i][1], cube.geom[i][2]);
231
 
  }
232
 
  cnt = 0;
233
 
  fprintf(cubfile,"%5d",cube.nmo_to_plot);
234
 
  ++cnt;
235
 
  for (i=0; i<cube.nmo_to_plot; ++i) {
236
 
    fprintf(cubfile,"%5d", cube.mos_to_plot[i] + 1);
237
 
    ++cnt;
238
 
    if ( (cnt == 8) && (i<cube.nmo_to_plot) ) {
239
 
      fprintf(cubfile,"\n");
240
 
      cnt = 0;
241
 
    }
242
 
  }
243
 
  fprintf(cubfile,"\n");
244
 
 
245
 
  /* write out grid points */
246
 
  for (x=0; x<cube.ngrid[0]; ++x) {
247
 
    for (y=0; y<cube.ngrid[1]; ++y) {
248
 
      cnt = 0;
249
 
      for (z=0; z<cube.ngrid[2]; ++z) {
250
 
        for (i=0; i<cube.nmo_to_plot; ++i) {
251
 
          fprintf(cubfile," %12.5E", cube.grid[x][y][z][i]);
252
 
          ++cnt;
253
 
          if (cnt == 6) {
254
 
            fprintf(cubfile,"\n");
255
 
            cnt = 0;
256
 
          }
257
 
        }
258
 
      }
259
 
      if (cnt != 0) fprintf(cubfile,"\n");
260
 
    }
261
 
  }
262
 
 
263
 
  fclose(cubfile);
264
 
  return;
265
 
}
266
 
 
267
 
void setup_delta(double **scf, double **u)
268
 
{
269
 
  int nmo, nao, i, I, j, errcod;
270
 
  int nirreps, nfzc, nfzv;
271
 
  int *order, *clsdpi, *openpi, *orbspi, *fruocc, *frdocc;
272
 
  double **scf_pitzer, *evals_pitzer, *evals, tval;
273
 
 
274
 
  nmo = params.nmo;
275
 
  nao = params.nao;
276
 
  nirreps = params.nirreps;
277
 
  orbspi = params.orbspi;
278
 
  openpi = params.openpi;
279
 
  clsdpi = params.clsdpi;
280
 
 
281
 
  /* read SCF eigenvectors in Pitzer ordering */
282
 
  chkpt_init(PSIO_OPEN_OLD);
283
 
  scf_pitzer = chkpt_rd_scf();
284
 
  /* u = chkpt_rd_usotao(); */
285
 
  evals_pitzer = chkpt_rd_evals();
286
 
  chkpt_close();
287
 
 
288
 
  frdocc = init_int_array(nirreps);
289
 
  fruocc = init_int_array(nirreps);
290
 
  errcod = ip_int_array("FROZEN_DOCC", frdocc, nirreps);
291
 
  errcod = ip_int_array("FROZEN_UOCC", fruocc, nirreps);
292
 
 
293
 
  nfzc = nfzv = 0;
294
 
  for(i=0; i < nirreps; i++) { 
295
 
    nfzc += frdocc[i];
296
 
    nfzv += fruocc[i];
297
 
  }
298
 
 
299
 
  /*** Get the Pitzer -> QT reordering array ***/
300
 
  order = init_int_array(nmo);
301
 
  reorder_qt(clsdpi, openpi, frdocc, fruocc, order, orbspi, nirreps);
302
 
 
303
 
  /*** Arrange the SCF eigenvectors into QT ordering ***/
304
 
 /* scf = block_matrix(nmo, nmo); */
305
 
  evals = (double *) malloc(nmo * sizeof(double));
306
 
  for(i=0; i < nmo; i++) {
307
 
      I = order[i];  /* Pitzer --> QT */
308
 
      for(j=0; j < nmo; j++) scf[j][I] = scf_pitzer[j][i];
309
 
      evals[I] = evals_pitzer[i];
310
 
    }
311
 
 
312
 
  /* determine HOMO and LUMO */
313
 
  tval = -1.0E6;
314
 
  for (i=0; i<params.nclsd; ++i) {
315
 
    if (evals[i] > tval) {
316
 
        tval = evals[i];
317
 
        params.homo = i;
318
 
    }
319
 
  }
320
 
 
321
 
  tval = 1.0E6;
322
 
  for (i=params.nclsd; i<params.nmo; ++i) {
323
 
    if (evals[i] < tval) {
324
 
        tval = evals[i];
325
 
        params.lumo = i ;
326
 
    }
327
 
  }
328
 
  fprintf(outfile,"orbitals in QT ordering, HOMO: %d LUMO: %d\n",
329
 
      params.homo+1, params.lumo+1);
330
 
 
331
 
  free(evals);
332
 
  free(evals_pitzer);
333
 
 
334
 
  free(order);
335
 
  /* free(clsdpi);
336
 
  free(openpi);
337
 
  free(orbspi); */
338
 
  free(fruocc);
339
 
  free(frdocc);
340
 
  free_block(scf_pitzer);
341
 
 
342
 
  return;
343
 
}
344