3
#include "prototypes.h"
6
static void print_2d_summary();
7
static void print_3d_summary();
8
static void print_grid_plotmtv();
9
static void print_grid_megapovray();
10
static void create_megapovray_file();
11
static void print_grid_gausscube();
12
static FILE* grid_file;
16
if (!strcmp(grid_format,"MEGAPOVPLUS"))
17
print_grid_megapovray();
18
else if (!strcmp(grid_format,"PLOTMTV"))
20
else if (!strcmp(grid_format,"GAUSSCUBE"))
21
print_grid_gausscube();
26
void print_2d_summary()
28
fprintf(outfile," --------------------------------------------------------------\n");
29
fprintf(outfile," *** Evaluating properties over a rectangular 2D grid ***\n");
30
fprintf(outfile," --------------------------------------------------------------\n\n");
31
fprintf(outfile," -Coordinates of the lower left, lower right, and upper left corners of\n");
32
fprintf(outfile," the grid rectangle (a.u.):\n");
33
fprintf(outfile," ** x");
34
fprintf(outfile," y z\n");
35
fprintf(outfile," ---- --------------------");
36
fprintf(outfile," -------------------- --------------------\n");
37
fprintf(outfile," LL %20.10lf %20.10lf %20.10lf\n",
38
grid_origin[0],grid_origin[1],grid_origin[2]);
39
fprintf(outfile," LR %20.10lf %20.10lf %20.10lf\n",
40
grid_origin[0]+grid_step_x[0]*nix,
41
grid_origin[1]+grid_step_x[1]*nix,
42
grid_origin[2]+grid_step_x[2]*nix);
43
fprintf(outfile," UL %20.10lf %20.10lf %20.10lf\n\n\n",
44
grid_origin[0]+grid_step_y[0]*niy,
45
grid_origin[1]+grid_step_y[1]*niy,
46
grid_origin[2]+grid_step_y[2]*niy);
49
void print_3d_summary()
51
fprintf(outfile," --------------------------------------------------------------\n");
52
fprintf(outfile," *** Evaluating properties over a rectangular 3D grid ***\n");
53
fprintf(outfile," --------------------------------------------------------------\n\n");
54
fprintf(outfile," -Coordinates of the lower left and upper right corners of\n");
55
fprintf(outfile," the grid box (a.u.):\n");
56
fprintf(outfile," ** x");
57
fprintf(outfile," y z\n");
58
fprintf(outfile," ---- --------------------");
59
fprintf(outfile," -------------------- --------------------\n");
60
fprintf(outfile," LL %20.10lf %20.10lf %20.10lf\n",
61
grid_origin[0],grid_origin[1],grid_origin[2]);
62
fprintf(outfile," UR %20.10lf %20.10lf %20.10lf\n",
63
grid_origin[0]+grid_step_x[0]*nix+grid_step_y[0]*niy+grid_step_z[0]*niz,
64
grid_origin[1]+grid_step_x[1]*nix+grid_step_y[1]*niy+grid_step_z[1]*niz,
65
grid_origin[2]+grid_step_x[2]*nix+grid_step_y[2]*niy+grid_step_z[2]*niz);
68
void print_grid_plotmtv()
71
double step_x, step_y, x, y;
77
grid_file = fopen("esp.dat","w");
82
grid_file = fopen("edens.dat","w");
84
grid_file = fopen("sdens.dat","w");
89
grid_file = fopen("edgrad.dat","w");
91
grid_file = fopen("sdgrad.dat","w");
96
grid_file = fopen("edlapl.dat","w");
98
grid_file = fopen("sdlapl.dat","w");
103
/* fprintf(grid_file,"%8.4lf %d %8.4lf %8.4lf %d %8.4lf %8.4lf %d %8.4lf\n",
104
grid_xyz0[0],nix,grid_xyz1[0],
105
grid_xyz0[1],niy,grid_xyz1[1],
106
grid_xyz0[2],niz,grid_xyz1[2]); */
111
fprintf(grid_file,"$DATA = CONTOUR\n");
112
fprintf(grid_file,"%% xmin = %lf xmax = %lf nx = %d\n",grid_xy0[0],grid_xy1[0],nix+1);
113
fprintf(grid_file,"%% ymin = %lf ymax = %lf ny = %d\n",grid_xy0[1],grid_xy1[1],niy+1);
114
fprintf(grid_file,"%% zmin = %lf zmax = %lf\n",grid_zmin,grid_zmax);
115
fprintf(grid_file,"%% contfill = T meshplot = T\n");
116
for(i=0;i<=niy;i++) {
118
if (grid_pts[j][i] < grid_zmin)
119
fprintf(grid_file," %lf ",grid_zmin);
121
if (grid_pts[j][i] > grid_zmax)
122
fprintf(grid_file," %lf ",grid_zmax);
124
fprintf(grid_file," %lf ",grid_pts[j][i]);
125
fprintf(grid_file,"\n");
130
fprintf(grid_file,"$DATA = VECTOR\n");
131
fprintf(grid_file,"%% xmin = %lf xmax = %lf\n",grid_xy0[0],grid_xy1[0]);
132
fprintf(grid_file,"%% ymin = %lf ymax = %lf\n",grid_xy0[1],grid_xy1[1]);
133
fprintf(grid_file,"%% zmin = %lf zmax = %lf\n",0.0,0.0);
134
fprintf(grid_file,"%% linecolor = 3\n");
135
fprintf(grid_file,"%% xlog = off vscale = 0.2\n\n");
136
step_x = (grid_xy1[0]-grid_xy0[0])/nix;
137
step_y = (grid_xy1[1]-grid_xy0[1])/niy;
138
for(i=0;i<=nix;i++) {
139
x = grid_xy0[0] + step_x*i;
140
for(j=0;j<=niy;j++) {
141
y = grid_xy0[1] + step_y*j;
142
if (fabs(grid_pts[i][j]) <= MAXDENSGRAD)
143
fprintf(grid_file,"%9.5lf %9.5lf %9.5lf %12.8lf %12.8lf %12.8lf\n",x,y,0.0,
144
grid_vecX[i][j],grid_vecY[i][j],0.0);
150
fprintf(grid_file,"$END\n");
155
void print_grid_megapovray()
158
double step_x, step_y, step_z, x, y, z;
162
/*--- Write out a data file ---*/
166
grid_file = fopen("mo.dat","w");
175
for(i=0;i<=nix;i++) {
176
fprintf(grid_file,"%15.8lf\n",grid3d_pts[0][i][j][k]);
183
/*--- Write out a command file ---*/
184
create_megapovray_file();
190
void create_megapovray_file()
193
double radius, midx, midy, midz;
194
double camera_distance, frame_width, frame_height, frame_center[3];
195
double dimx, dimy, dimz, maxdim;
196
double **grid_geom, **Rgrid;
200
dimx = grid_xyz1[0] - grid_xyz0[0];
201
dimy = grid_xyz1[1] - grid_xyz0[1];
202
dimz = grid_xyz1[2] - grid_xyz0[2];
203
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
204
maxdim = MAX(dimz,MAX(dimx,dimy));
205
frame_width = maxdim + 2.0;
206
frame_height = maxdim + 2.0;
207
frame_center[0] = grid_origin[0] + 0.5*dimx;
208
frame_center[1] = grid_origin[1] + 0.5*dimy;
209
frame_center[2] = grid_origin[2] + 0.5*dimz;
210
camera_distance = 3.0*frame_height;
212
/* compute atomic coordinates in the grid coordinate system,
213
since MegaPov simply assumes a simple XYZ grid */
214
grid_geom = block_matrix(natom,3);
215
Rgrid = block_matrix(3,3);
217
Rgrid[i][0] = grid_unit_x[i];
218
Rgrid[i][1] = grid_unit_y[i];
219
Rgrid[i][2] = grid_unit_z[i];
221
mmult(geom,0,Rgrid,0,grid_geom,0,natom,3,3,0);
224
compute_connectivity();
226
mpvfile = fopen("mo.pov","w");
227
fprintf(mpvfile,"// File: mo.pov\n");
228
fprintf(mpvfile,"// Creator: oeprop (Psi 3.2)\n");
229
fprintf(mpvfile,"// Version: MegaPov 0.5\n\n");
231
fprintf(mpvfile,"#version unofficial MegaPov 0.5;\n\n");
233
fprintf(mpvfile,"// Camera\n");
234
fprintf(mpvfile,"camera {\n");
235
fprintf(mpvfile," orthographic\n");
236
fprintf(mpvfile," location <%lf, %lf, %lf>\n",
237
camera_distance,frame_center[1],frame_center[2]);
238
fprintf(mpvfile," look_at <%lf, %lf, %lf>\n",
239
frame_center[0],frame_center[1],frame_center[2]);
240
fprintf(mpvfile," up <%lf, %lf, %lf>\n",
241
frame_center[0],frame_center[1],frame_center[2]+frame_height);
242
fprintf(mpvfile," right <%lf, %lf, %lf>\n",
243
frame_center[0],frame_center[1]+frame_width,frame_center[2]);
244
fprintf(mpvfile,"}\n\n");
246
fprintf(mpvfile,"// Light\n");
247
fprintf(mpvfile,"light_source {<%lf, %lf, %lf> color rgb <1.0, 1.0, 1.0>}\n\n",
248
camera_distance*10.0,camera_distance*(-10.0),camera_distance*(-10.0));
250
fprintf(mpvfile,"// Background\n");
251
fprintf(mpvfile,"background { color rgb <1.0, 1.0, 1.0>}\n\n");
253
fprintf(mpvfile,"union {\n");
254
fprintf(mpvfile,"// Objects\n");
255
fprintf(mpvfile," // union finish\n");
256
fprintf(mpvfile," #declare F = finish {specular 0.4 roughness 0.005 diffuse 0.8 ambient 0.2}\n");
257
fprintf(mpvfile," // transparency of atomic surfaces \n");
258
fprintf(mpvfile," #declare T = 0;\n");
259
fprintf(mpvfile," // no_shadow\n");
260
fprintf(mpvfile," // hollow\n\n");
262
fprintf(mpvfile," // Atoms\n");
263
for(i=0;i<natom;i++) {
265
fprintf(mpvfile," // Atom # %d, charge %d\n",i,z);
266
fprintf(mpvfile," object {\n");
267
fprintf(mpvfile," sphere { < %lf, %lf, %lf > ",grid_geom[i][0],grid_geom[i][1],grid_geom[i][2]);
274
fprintf(mpvfile,"%lf }\n",radius);
275
fprintf(mpvfile," texture {\n");
276
if (z <= LAST_RGB_INDEX)
277
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[z][0],atomic_rgb[z][1],atomic_rgb[z][2]);
279
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[0][0],atomic_rgb[0][1],atomic_rgb[0][2]);
280
fprintf(mpvfile," finish {F}\n");
281
fprintf(mpvfile," }\n");
282
fprintf(mpvfile," }\n");
284
fprintf(mpvfile,"\n");
286
fprintf(mpvfile," // Bonds\n");
289
if (connectivity[i][j]) {
290
fprintf(mpvfile," // Bond between atoms %d and %d\n",i,j);
291
midx = 0.5*(grid_geom[i][0] + grid_geom[j][0]);
292
midy = 0.5*(grid_geom[i][1] + grid_geom[j][1]);
293
midz = 0.5*(grid_geom[i][2] + grid_geom[j][2]);
294
fprintf(mpvfile," object {\n");
295
fprintf(mpvfile," cylinder { < %lf, %lf, %lf > < %lf, %lf, %lf > 0.15 }\n",
296
grid_geom[i][0],grid_geom[i][1],grid_geom[i][2],
298
fprintf(mpvfile," texture {\n");
300
if (z <= LAST_RGB_INDEX)
301
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[z][0],atomic_rgb[z][1],atomic_rgb[z][2]);
303
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[0][0],atomic_rgb[0][1],atomic_rgb[0][2]);
304
fprintf(mpvfile," finish {F}\n");
305
fprintf(mpvfile," }\n");
306
fprintf(mpvfile," }\n");
307
fprintf(mpvfile," object {\n");
308
fprintf(mpvfile," cylinder { < %lf, %lf, %lf > < %lf, %lf, %lf > 0.15 }\n",
309
grid_geom[j][0],grid_geom[j][1],grid_geom[j][2],
311
fprintf(mpvfile," texture {\n");
313
if (z <= LAST_RGB_INDEX)
314
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[z][0],atomic_rgb[z][1],atomic_rgb[z][2]);
316
fprintf(mpvfile," pigment {color rgbt <%lf, %lf, %lf, T>}\n",atomic_rgb[0][0],atomic_rgb[0][1],atomic_rgb[0][2]);
317
fprintf(mpvfile," finish {F}\n");
318
fprintf(mpvfile," }\n");
319
fprintf(mpvfile," }\n");
322
fprintf(mpvfile," #declare ISO = 1;\n");
323
fprintf(mpvfile," #declare NX = %d;\n",nix);
324
fprintf(mpvfile," #declare NY = %d;\n",niy);
325
fprintf(mpvfile," #declare NZ = %d;\n",niz);
326
fprintf(mpvfile," #declare RX = %lf;\n",grid_origin[0]);
327
fprintf(mpvfile," #declare RY = %lf;\n",grid_origin[1]);
328
fprintf(mpvfile," #declare RZ = %lf;\n",grid_origin[2]);
329
fprintf(mpvfile," #declare LX = %lf;\n",grid_xyz1[0] - grid_xyz0[0]);
330
fprintf(mpvfile," #declare LY = %lf;\n",grid_xyz1[1] - grid_xyz0[1]);
331
fprintf(mpvfile," #declare LZ = %lf;\n",grid_xyz1[2] - grid_xyz0[2]);
332
fprintf(mpvfile," #declare LEVSCALE = 11;\n\n");
334
fprintf(mpvfile," #ifdef (ISO)\n");
335
fprintf(mpvfile," #declare wfun1 = function {\"data_3D_3\",\n");
336
fprintf(mpvfile," <-LEVSCALE>, library \"i_dat3d\", \"mo.dat\", <NX+1,NY+1,NZ+1,0>}\n");
337
fprintf(mpvfile," #declare wfun2 = function {\"data_3D_3\",\n");
338
fprintf(mpvfile," <LEVSCALE>, library \"i_dat3d\", \"mo.dat\", <NX+1,NY+1,NZ+1,0>}\n");
339
fprintf(mpvfile," union {\n");
340
fprintf(mpvfile," isosurface {\n");
341
fprintf(mpvfile," function {wfun1}\n");
342
fprintf(mpvfile," contained_by { box {<0,0,0>,<NX,NY,NZ>} }\n");
343
fprintf(mpvfile," threshold -1\n");
344
fprintf(mpvfile," max_gradient 1.000\n");
345
fprintf(mpvfile," eval\n");
346
fprintf(mpvfile," texture { pigment { color rgbt <1.0,0.0,0.0,0.8> } }\n");
347
fprintf(mpvfile," }\n");
348
fprintf(mpvfile," isosurface {\n");
349
fprintf(mpvfile," function {wfun2}\n");
350
fprintf(mpvfile," contained_by { box {<0,0,0>,<NX,NY,NZ>} }\n");
351
fprintf(mpvfile," threshold -1\n");
352
fprintf(mpvfile," max_gradient 1.000\n");
353
fprintf(mpvfile," eval\n");
354
fprintf(mpvfile," texture { pigment { color rgbt <0.0,1.0,0.0,0.8> } }\n");
355
fprintf(mpvfile," }\n");
356
fprintf(mpvfile," translate <RX*NX/LX,RY*NY/LY,RZ*NZ/LZ>\n");
357
fprintf(mpvfile," scale <LX/NX,LY/NY,LZ/NZ>\n");
358
fprintf(mpvfile," }\n");
359
fprintf(mpvfile," #end\n\n");
361
fprintf(mpvfile," // Rotate\n");
362
fprintf(mpvfile," rotate 0*z\n");
363
fprintf(mpvfile,"}\n");
367
free_block(grid_geom);
373
void print_grid_gausscube()
378
double step_x, step_y, step_z, x, y, z;
385
grid_file = fopen("mo.cube","w");
388
grid_file = fopen("dens.cube","w");
392
/*------------------------------------------------------------------------
393
Gaussian Cube standard overhead as described in G94 programmer's manual
394
------------------------------------------------------------------------*/
395
/* Comment and subcomment */
396
fprintf(grid_file,"Gaussian Cube file created by OEPROP (Psi 3.2)\n");
397
fprintf(grid_file,"Calculation title: %s\n",title);
399
fprintf(grid_file,"%5d",-1*natom);
401
fprintf(grid_file,"%5d",natom);
402
fprintf(grid_file,"%12.6lf%12.6lf%12.6lf\n", grid_origin[0], grid_origin[1], grid_origin[2]);
403
fprintf(grid_file,"%5d%12.6lf%12.6lf%12.6lf\n", nix+1, grid_step_x[0], grid_step_x[1], grid_step_x[2]);
404
fprintf(grid_file,"%5d%12.6lf%12.6lf%12.6lf\n", niy+1, grid_step_y[0], grid_step_y[1], grid_step_y[2]);
405
fprintf(grid_file,"%5d%12.6lf%12.6lf%12.6lf\n", niz+1, grid_step_z[0], grid_step_z[1], grid_step_z[2]);
406
for(atom=0;atom<natom;atom++) {
407
fprintf(grid_file,"%5d%12.6lf%12.6lf%12.6lf%12.6lf\n", (int)zvals[atom], zvals[atom],
408
geom[atom][0], geom[atom][1], geom[atom][2]);
411
fprintf(grid_file,"%5d",num_mos_to_plot);
412
for(mo=0;mo<num_mos_to_plot;mo++)
413
fprintf(grid_file,"%5d",mos_to_plot[mo]);
414
fprintf(grid_file,"\n");
420
for(j=0;j<=niy;j++) {
423
for(mo=0;mo<num_mos_to_plot;mo++) {
424
fprintf(grid_file,"%13.5E",grid3d_pts[mo][i][j][k]);
426
if (n_per_line == 6) {
427
fprintf(grid_file,"\n");
432
fprintf(grid_file,"\n");
438
for(j=0;j<=niy;j++) {
440
for(k=0;k<=niz;k++) {
441
fprintf(grid_file,"%13.5E",grid3d_pts[0][i][j][k]);
443
if (n_per_line == 6) {
444
fprintf(grid_file,"\n");
449
fprintf(grid_file,"\n");
added
removed
src/bin/oeprop/print_grid.c
