src/tools/gmx_editconf.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Green Red Orange Magenta Azure Cyan Skyblue
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include <ctype.h>
  43 #include "pdbio.h"
  44 #include "confio.h"
  45 #include "symtab.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "copyrite.h"
  49 #include "statutil.h"
  50 #include "string2.h"
  51 #include "strdb.h"
  52 #include "index.h"
  53 #include "vec.h"
  54 #include "typedefs.h"
  55 #include "gbutil.h"
  56 #include "strdb.h"
  57 #include "physics.h"
  58 #include "atomprop.h"
  59 #include "tpxio.h"
  60 #include "pbc.h"
  61 #include "princ.h"
  62 #include "txtdump.h"
  63 #include "viewit.h"
  64 #include "rmpbc.h"
  65 #include "gmx_ana.h"
  66
  67 typedef struct
  68 {
  69     char sanm[12];
  70     int natm;
  71     int nw;
  72     char anm[6][12];
  73 } t_simat;
  74
  75 typedef struct
  76 {
  77     char reso[12];
  78     char resn[12];
  79     int nsatm;
  80     t_simat sat[3];
  81 } t_simlist;
  82 static const char *pdbtp[epdbNR] =
  83     { "ATOM  ", "HETATM" };
  84
  85 real calc_mass(t_atoms *atoms, bool bGetMass, gmx_atomprop_t aps)
  86 {
  87     real tmass;
  88     int i;
  89
  90     tmass = 0;
  91     for (i = 0; (i < atoms->nr); i++)
  92     {
  93         if (bGetMass)
  94         {
  95             gmx_atomprop_query(aps, epropMass,
  96                                *atoms->resinfo[atoms->atom[i].resind].name,
  97                                *atoms->atomname[i], &(atoms->atom[i].m));
  98         }
  99         tmass += atoms->atom[i].m;
 100     }
 101
 102     return tmass;
 103 }
 104
 105 real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
 106                rvec max, bool bDiam)
 107 {
 108     real diam2, d;
 109     char *grpnames;
 110     int ii, i, j;
 111
 112     clear_rvec(geom_center);
 113     diam2 = 0;
 114     if (isize == 0)
 115     {
 116         clear_rvec(min);
 117         clear_rvec(max);
 118     }
 119     else
 120     {
 121         if (index)
 122             ii = index[0];
 123         else
 124             ii = 0;
 125         for (j = 0; j < DIM; j++)
 126             min[j] = max[j] = x[ii][j];
 127         for (i = 0; i < isize; i++)
 128         {
 129             if (index)
 130                 ii = index[i];
 131             else
 132                 ii = i;
 133             rvec_inc(geom_center, x[ii]);
 134             for (j = 0; j < DIM; j++)
 135             {
 136                 if (x[ii][j] < min[j])
 137                     min[j] = x[ii][j];
 138                 if (x[ii][j] > max[j])
 139                     max[j] = x[ii][j];
 140             }
 141             if (bDiam)
 142             {
 143                 if (index)
 144                     for (j = i + 1; j < isize; j++)
 145                     {
 146                         d = distance2(x[ii], x[index[j]]);
 147                         diam2 = max(d,diam2);
 148                     }
 149                 else
 150                     for (j = i + 1; j < isize; j++)
 151                     {
 152                         d = distance2(x[i], x[j]);
 153                         diam2 = max(d,diam2);
 154                     }
 155             }
 156         }
 157         svmul(1.0 / isize, geom_center, geom_center);
 158     }
 159
 160     return sqrt(diam2);
 161 }
 162
 163 void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
 164 {
 165     int i;
 166     rvec shift;
 167
 168     rvec_sub(center, geom_cent, shift);
 169
 170     printf("    shift       :%7.3f%7.3f%7.3f (nm)\n", shift[XX], shift[YY],
 171            shift[ZZ]);
 172
 173     for (i = 0; (i < natom); i++)
 174         rvec_inc(x[i], shift);
 175 }
 176
 177 void scale_conf(int natom, rvec x[], matrix box, rvec scale)
 178 {
 179     int i, j;
 180
 181     for (i = 0; i < natom; i++)
 182     {
 183         for (j = 0; j < DIM; j++)
 184             x[i][j] *= scale[j];
 185     }
 186     for (i = 0; i < DIM; i++)
 187         for (j = 0; j < DIM; j++)
 188             box[i][j] *= scale[j];
 189 }
 190
 191 void read_bfac(const char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
 192 {
 193     int i;
 194     char **bfac_lines;
 195
 196     *n_bfac = get_lines(fn, &bfac_lines);
 197     snew(*bfac_val, *n_bfac);
 198     snew(*bfac_nr, *n_bfac);
 199     fprintf(stderr, "Reading %d B-factors from %s\n", *n_bfac, fn);
 200     for (i = 0; (i < *n_bfac); i++)
 201     {
 202         /*fprintf(stderr, "Line %d: %s",i,bfac_lines[i]);*/
 203         sscanf(bfac_lines[i], "%d %lf", &(*bfac_nr)[i], &(*bfac_val)[i]);
 204         /*fprintf(stderr," nr %d val %g\n",(*bfac_nr)[i],(*bfac_val)[i]);*/
 205     }
 206
 207 }
 208
 209 void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
 210                        double *bfac, int *bfac_nr, bool peratom)
 211 {
 212     FILE *out;
 213     real bfac_min, bfac_max;
 214     int i, n;
 215     bool found;
 216
 217     bfac_max = -1e10;
 218     bfac_min = 1e10;
 219     for (i = 0; (i < n_bfac); i++)
 220     {
 221         if (bfac_nr[i] - 1 >= atoms->nres)
 222             peratom = TRUE;
 223         /*    if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
 224          gmx_fatal(FARGS,"Index of B-Factor %d is out of range: %d (%g)",
 225          i+1,bfac_nr[i],bfac[i]); */
 226         if (bfac[i] > bfac_max)
 227             bfac_max = bfac[i];
 228         if (bfac[i] < bfac_min)
 229             bfac_min = bfac[i];
 230     }
 231     while ((bfac_max > 99.99) || (bfac_min < -99.99))
 232     {
 233         fprintf(stderr,
 234                 "Range of values for B-factors too large (min %g, max %g) "
 235                     "will scale down a factor 10\n", bfac_min, bfac_max);
 236         for (i = 0; (i < n_bfac); i++)
 237             bfac[i] /= 10;
 238         bfac_max /= 10;
 239         bfac_min /= 10;
 240     }
 241     while ((fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5))
 242     {
 243         fprintf(stderr,
 244                 "Range of values for B-factors too small (min %g, max %g) "
 245                     "will scale up a factor 10\n", bfac_min, bfac_max);
 246         for (i = 0; (i < n_bfac); i++)
 247             bfac[i] *= 10;
 248         bfac_max *= 10;
 249         bfac_min *= 10;
 250     }
 251
 252     for (i = 0; (i < natoms); i++)
 253         atoms->pdbinfo[i].bfac = 0;
 254
 255     if (!peratom)
 256     {
 257         fprintf(stderr, "Will attach %d B-factors to %d residues\n", n_bfac,
 258                 nres);
 259         for (i = 0; (i < n_bfac); i++)
 260         {
 261             found = FALSE;
 262             for (n = 0; (n < natoms); n++)
 263                 if (bfac_nr[i] == atoms->resinfo[atoms->atom[n].resind].nr)
 264                 {
 265                     atoms->pdbinfo[n].bfac = bfac[i];
 266                     found = TRUE;
 267                 }
 268             if (!found)
 269             {
 270                 gmx_warning("Residue nr %d not found\n", bfac_nr[i]);
 271             }
 272         }
 273     }
 274     else
 275     {
 276         fprintf(stderr, "Will attach %d B-factors to %d atoms\n", n_bfac,
 277                 natoms);
 278         for (i = 0; (i < n_bfac); i++)
 279         {
 280             atoms->pdbinfo[bfac_nr[i] - 1].bfac = bfac[i];
 281         }
 282     }
 283 }
 284
 285 void pdb_legend(FILE *out, int natoms, int nres, t_atoms *atoms, rvec x[])
 286 {
 287     real bfac_min, bfac_max, xmin, ymin, zmin;
 288     int i;
 289     int space = ' ';
 290
 291     bfac_max = -1e10;
 292     bfac_min = 1e10;
 293     xmin = 1e10;
 294     ymin = 1e10;
 295     zmin = 1e10;
 296     for (i = 0; (i < natoms); i++)
 297     {
 298         xmin = min(xmin,x[i][XX]);
 299         ymin = min(ymin,x[i][YY]);
 300         zmin = min(zmin,x[i][ZZ]);
 301         bfac_min = min(bfac_min,atoms->pdbinfo[i].bfac);
 302         bfac_max = max(bfac_max,atoms->pdbinfo[i].bfac);
 303     }
 304     fprintf(stderr, "B-factors range from %g to %g\n", bfac_min, bfac_max);
 305     for (i = 1; (i < 12); i++)
 306     {
 307         fprintf(out,
 308                 "%-6s%5u  %-4.4s%3.3s %c%4d%c   %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 309                 "ATOM  ", natoms + 1 + i, "CA", "LEG", space, nres + 1, space,
 310                 (xmin + (i * 0.12)) * 10, ymin * 10, zmin * 10, 1.0, bfac_min
 311                     + ((i - 1.0) * (bfac_max - bfac_min) / 10));
 312     }
 313 }
 314
 315 void visualize_images(const char *fn, int ePBC, matrix box)
 316 {
 317     t_atoms atoms;
 318     rvec *img;
 319     char *c, *ala;
 320     int nat, i;
 321
 322     nat = NTRICIMG + 1;
 323     init_t_atoms(&atoms, nat, FALSE);
 324     atoms.nr = nat;
 325     snew(img,nat);
 326     /* FIXME: Constness should not be cast away */
 327     c = (char *) "C";
 328     ala = (char *) "ALA";
 329     for (i = 0; i < nat; i++)
 330     {
 331         atoms.atomname[i] = &c;
 332         atoms.atom[i].resind = i;
 333         atoms.resinfo[i].name = &ala;
 334         atoms.resinfo[i].nr = i + 1;
 335         atoms.resinfo[i].chainid = 'A' + i / NCUCVERT;
 336     }
 337     calc_triclinic_images(box, img + 1);
 338
 339     write_sto_conf(fn, "Images", &atoms, img, NULL, ePBC, box);
 340
 341     free_t_atoms(&atoms, FALSE);
 342     sfree(img);
 343 }
 344
 345 void visualize_box(FILE *out, int a0, int r0, matrix box, rvec gridsize)
 346 {
 347     int *edge;
 348     rvec *vert, shift;
 349     int nx, ny, nz, nbox, nat;
 350     int i, j, x, y, z;
 351     int rectedge[24] =
 352         { 0, 1, 1, 3, 3, 2, 0, 2, 0, 4, 1, 5, 3, 7, 2, 6, 4, 5, 5, 7, 7, 6, 6,
 353             4 };
 354
 355     a0++;
 356     r0++;
 357
 358     nx = (int) (gridsize[XX] + 0.5);
 359     ny = (int) (gridsize[YY] + 0.5);
 360     nz = (int) (gridsize[ZZ] + 0.5);
 361     nbox = nx * ny * nz;
 362     if (TRICLINIC(box))
 363     {
 364         nat = nbox * NCUCVERT;
 365         snew(vert,nat);
 366         calc_compact_unitcell_vertices(ecenterDEF, box, vert);
 367         j = 0;
 368         for (z = 0; z < nz; z++)
 369             for (y = 0; y < ny; y++)
 370                 for (x = 0; x < nx; x++)
 371                 {
 372                     for (i = 0; i < DIM; i++)
 373                         shift[i] = x * box[0][i] + y * box[1][i] + z
 374                             * box[2][i];
 375                     for (i = 0; i < NCUCVERT; i++)
 376                     {
 377                         rvec_add(vert[i], shift, vert[j]);
 378                         j++;
 379                     }
 380                 }
 381
 382         for (i = 0; i < nat; i++)
 383         {
 384             fprintf(out, pdbformat, "ATOM", a0 + i, "C", "BOX", 'K' + i
 385                 / NCUCVERT, r0 + i, 10 * vert[i][XX], 10 * vert[i][YY], 10
 386                 * vert[i][ZZ]);
 387             fprintf(out, "\n");
 388         }
 389
 390         edge = compact_unitcell_edges();
 391         for (j = 0; j < nbox; j++)
 392             for (i = 0; i < NCUCEDGE; i++)
 393                 fprintf(out, "CONECT%5d%5d\n", a0 + j * NCUCVERT + edge[2 * i],
 394                         a0 + j * NCUCVERT + edge[2 * i + 1]);
 395
 396         sfree(vert);
 397     }
 398     else
 399     {
 400         i = 0;
 401         for (z = 0; z <= 1; z++)
 402             for (y = 0; y <= 1; y++)
 403                 for (x = 0; x <= 1; x++)
 404                 {
 405                     fprintf(out, pdbformat, "ATOM", a0 + i, "C", "BOX", 'K' + i
 406                         / 8, r0 + i, x * 10 * box[XX][XX],
 407                             y * 10 * box[YY][YY], z * 10 * box[ZZ][ZZ]);
 408                     fprintf(out, "\n");
 409                     i++;
 410                 }
 411         for (i = 0; i < 24; i += 2)
 412             fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i
 413                 + 1]);
 414     }
 415 }
 416
 417 void calc_rotmatrix(rvec principal_axis, rvec targetvec, matrix rotmatrix)
 418 {
 419         rvec rotvec;
 420         real ux,uy,uz,costheta,sintheta;
 421
 422         costheta = cos_angle(principal_axis,targetvec);
 423         sintheta=sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
 424
 425         /* Determine rotation from cross product with target vector */
 426         cprod(principal_axis,targetvec,rotvec);
 427         unitv(rotvec,rotvec);
 428         printf("Aligning %g %g %g to %g %g %g : xprod  %g %g %g\n",
 429                 principal_axis[XX],principal_axis[YY],principal_axis[ZZ],targetvec[XX],targetvec[YY],targetvec[ZZ],
 430                 rotvec[XX],rotvec[YY],rotvec[ZZ]);
 431
 432         ux=rotvec[XX];
 433         uy=rotvec[YY];
 434         uz=rotvec[ZZ];
 435         rotmatrix[0][0]=ux*ux + (1.0-ux*ux)*costheta;
 436         rotmatrix[0][1]=ux*uy*(1-costheta)-uz*sintheta;
 437         rotmatrix[0][2]=ux*uz*(1-costheta)+uy*sintheta;
 438         rotmatrix[1][0]=ux*uy*(1-costheta)+uz*sintheta;
 439         rotmatrix[1][1]=uy*uy + (1.0-uy*uy)*costheta;
 440         rotmatrix[1][2]=uy*uz*(1-costheta)-ux*sintheta;
 441         rotmatrix[2][0]=ux*uz*(1-costheta)-uy*sintheta;
 442         rotmatrix[2][1]=uy*uz*(1-costheta)+ux*sintheta;
 443         rotmatrix[2][2]=uz*uz + (1.0-uz*uz)*costheta;
 444
 445         printf("Rotation matrix: \n%g %g %g\n%g %g %g\n%g %g %g\n",
 446                 rotmatrix[0][0],rotmatrix[0][1],rotmatrix[0][2],
 447                 rotmatrix[1][0],rotmatrix[1][1],rotmatrix[1][2],
 448                 rotmatrix[2][0],rotmatrix[2][1],rotmatrix[2][2]);
 449 }
 450
 451 int gmx_editconf(int argc, char *argv[])
 452 {
 453     const char
 454         *desc[] =
 455             {
 456                 "editconf converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
 457                 "or [TT].pdb[tt].",
 458                 "[PAR]",
 459                 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
 460                 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
 461                 "will center the system in the box, unless [TT]-noc[tt] is used.",
 462                 "[PAR]",
 463                 "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
 464                 "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
 465                 "[TT]dodecahedron[tt] represents a rhombic dodecahedron and",
 466                 "[TT]octahedron[tt] is a truncated octahedron.",
 467                 "The last two are special cases of a triclinic box.",
 468                 "The length of the three box vectors of the truncated octahedron is the",
 469                 "shortest distance between two opposite hexagons.",
 470                 "The volume of a dodecahedron is 0.71 and that of a truncated octahedron",
 471                 "is 0.77 of that of a cubic box with the same periodic image distance.",
 472                 "[PAR]",
 473                 "Option [TT]-box[tt] requires only",
 474                 "one value for a cubic box, dodecahedron and a truncated octahedron.",
 475                 "[PAR]",
 476                 "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the x, y",
 477                 "and z directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
 478                 "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
 479                 "to the diameter of the system (largest distance between atoms) plus twice",
 480                 "the specified distance.",
 481                 "[PAR]",
 482                 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
 483                 "a triclinic box and can not be used with option [TT]-d[tt].",
 484                 "[PAR]",
 485                 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
 486                 "can be selected for calculating the size and the geometric center,",
 487                 "otherwise the whole system is used.",
 488                 "[PAR]",
 489                 "[TT]-rotate[tt] rotates the coordinates and velocities.",
 490                 "[PAR]",
 491                 "[TT]-princ[tt] aligns the principal axes of the system along the",
 492                 "coordinate axes, this may allow you to decrease the box volume,",
 493                 "but beware that molecules can rotate significantly in a nanosecond.",
 494                 "[PAR]",
 495                 "Scaling is applied before any of the other operations are",
 496                 "performed. Boxes and coordinates can be scaled to give a certain density (option",
 497                 "[TT]-density[tt]). Note that this may be inaccurate in case a gro",
 498                 "file is given as input. A special feature of the scaling option, when the",
 499                 "factor -1 is given in one dimension, one obtains a mirror image,",
 500                 "mirrored in one of the plains, when one uses -1 in three dimensions",
 501                 "a point-mirror image is obtained.[PAR]",
 502                 "Groups are selected after all operations have been applied.[PAR]",
 503                 "Periodicity can be removed in a crude manner.",
 504                 "It is important that the box sizes at the bottom of your input file",
 505                 "are correct when the periodicity is to be removed.",
 506                 "[PAR]",
 507                 "When writing [TT].pdb[tt] files, B-factors can be",
 508                 "added with the [TT]-bf[tt] option. B-factors are read",
 509                 "from a file with with following format: first line states number of",
 510                 "entries in the file, next lines state an index",
 511                 "followed by a B-factor. The B-factors will be attached per residue",
 512                 "unless an index is larger than the number of residues or unless the",
 513                 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
 514                 "be added instead of B-factors. [TT]-legend[tt] will produce",
 515                 "a row of CA atoms with B-factors ranging from the minimum to the",
 516                 "maximum value found, effectively making a legend for viewing.",
 517                 "[PAR]",
 518                 "With the option -mead a special pdb (pqr) file for the MEAD electrostatics",
 519                 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
 520                 "is that the input file is a run input file.",
 521                 "The B-factor field is then filled with the Van der Waals radius",
 522                 "of the atoms while the occupancy field will hold the charge.",
 523                 "[PAR]",
 524                 "The option -grasp is similar, but it puts the charges in the B-factor",
 525                 "and the radius in the occupancy.",
 526                 "[PAR]",
 527                 "Option [TT]-align[tt] allows alignment",
 528                 "of the principal axis of a specified group against the given vector, ",
 529                                 "with an optional center of rotation specified by [TT]-aligncenter[tt].",
 530                 "[PAR]",
 531                 "Finally with option [TT]-label[tt] editconf can add a chain identifier",
 532                 "to a pdb file, which can be useful for analysis with e.g. rasmol.",
 533                     "[PAR]",
 534                 "To convert a truncated octrahedron file produced by a package which uses",
 535                 "a cubic box with the corners cut off (such as Gromos) use:[BR]",
 536                 "[TT]editconf -f <in> -rotate 0 45 35.264 -bt o -box <veclen> -o <out>[tt][BR]",
 537                 "where [TT]veclen[tt] is the size of the cubic box times sqrt(3)/2." };
 538     const char *bugs[] =
 539         {
 540             "For complex molecules, the periodicity removal routine may break down, ",
 541                 "in that case you can use trjconv." };
 542     static real dist = 0.0, rbox = 0.0, to_diam = 0.0;
 543     static bool bNDEF = FALSE, bRMPBC = FALSE, bCenter = FALSE, bReadVDW =
 544         FALSE, bCONECT = FALSE;
 545     static bool peratom = FALSE, bLegend = FALSE, bOrient = FALSE, bMead =
 546         FALSE, bGrasp = FALSE, bSig56 = FALSE;
 547     static rvec scale =
 548         { 1, 1, 1 }, newbox =
 549         { 0, 0, 0 }, newang =
 550         { 90, 90, 90 };
 551     static real rho = 1000.0, rvdw = 0.12;
 552     static rvec center =
 553         { 0, 0, 0 }, translation =
 554         { 0, 0, 0 }, rotangles =
 555         { 0, 0, 0 }, aligncenter =
 556                 { 0, 0, 0 }, targetvec =
 557         { 0, 0, 0 };
 558     static const char *btype[] =
 559         { NULL, "triclinic", "cubic", "dodecahedron", "octahedron", NULL },
 560         *label = "A";
 561     static rvec visbox =
 562         { 0, 0, 0 };
 563     t_pargs
 564         pa[] =
 565             {
 566                     { "-ndef", FALSE, etBOOL,
 567                         { &bNDEF }, "Choose output from default index groups" },
 568                     { "-visbox", FALSE, etRVEC,
 569                         { visbox },
 570                         "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
 571                     { "-bt", FALSE, etENUM,
 572                         { btype }, "Box type for -box and -d" },
 573                     { "-box", FALSE, etRVEC,
 574                         { newbox }, "Box vector lengths (a,b,c)" },
 575                     { "-angles", FALSE, etRVEC,
 576                         { newang }, "Angles between the box vectors (bc,ac,ab)" },
 577                     { "-d", FALSE, etREAL,
 578                         { &dist }, "Distance between the solute and the box" },
 579                     { "-c", FALSE, etBOOL,
 580                         { &bCenter },
 581                         "Center molecule in box (implied by -box and -d)" },
 582                     { "-center", FALSE, etRVEC,
 583                         { center }, "Coordinates of geometrical center" },
 584                     { "-aligncenter", FALSE, etRVEC,
 585                         { aligncenter }, "Center of rotation for alignment" },
 586                     { "-align", FALSE, etRVEC,
 587                         { targetvec },
 588                         "Align to target vector" },
 589                     { "-translate", FALSE, etRVEC,
 590                         { translation }, "Translation" },
 591                     { "-rotate", FALSE, etRVEC,
 592                         { rotangles },
 593                         "Rotation around the X, Y and Z axes in degrees" },
 594                     { "-princ", FALSE, etBOOL,
 595                         { &bOrient },
 596                         "Orient molecule(s) along their principal axes" },
 597                     { "-scale", FALSE, etRVEC,
 598                         { scale }, "Scaling factor" },
 599                     { "-density", FALSE, etREAL,
 600                         { &rho },
 601                         "Density (g/l) of the output box achieved by scaling" },
 602                     { "-pbc", FALSE, etBOOL,
 603                         { &bRMPBC },
 604                         "Remove the periodicity (make molecule whole again)" },
 605                     { "-grasp", FALSE, etBOOL,
 606                         { &bGrasp },
 607                         "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
 608                     {
 609                         "-rvdw", FALSE, etREAL,
 610                          { &rvdw },
 611                         "Default Van der Waals radius (in nm) if one can not be found in the database or if no parameters are present in the topology file" },
 612                     { "-sig56", FALSE, etREAL,
 613                         { &bSig56 },
 614                         "Use rmin/2 (minimum in the Van der Waals potential) rather than sigma/2 " },
 615                     {
 616                         "-vdwread", FALSE, etBOOL,
 617                         { &bReadVDW },
 618                         "Read the Van der Waals radii from the file vdwradii.dat rather than computing the radii based on the force field" },
 619                     { "-atom", FALSE, etBOOL,
 620                         { &peratom }, "Force B-factor attachment per atom" },
 621                     { "-legend", FALSE, etBOOL,
 622                         { &bLegend }, "Make B-factor legend" },
 623                     { "-label", FALSE, etSTR,
 624                         { &label }, "Add chain label for all residues" },
 625                     {
 626                         "-conect", FALSE, etBOOL,
 627                         { &bCONECT },
 628                         "Add CONECT records to a pdb file when written. Can only be done when a topology is present" } };
 629 #define NPA asize(pa)
 630
 631     FILE *out;
 632     const char *infile, *outfile;
 633     char title[STRLEN];
 634     int outftp, inftp, natom, i, j, n_bfac, itype, ntype;
 635     double *bfac = NULL, c6, c12;
 636     int *bfac_nr = NULL;
 637     t_topology *top = NULL;
 638     t_atoms atoms;
 639     char *grpname, *sgrpname, *agrpname;
 640     int isize, ssize, tsize, asize;
 641     atom_id *index, *sindex, *tindex, *aindex;
 642     rvec *x, *v, gc, min, max, size;
 643     int ePBC;
 644     matrix box,rotmatrix,trans;
 645         rvec princd,tmpvec;
 646     bool bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
 647     bool bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
 648     real xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
 649     gmx_atomprop_t aps;
 650     gmx_conect conect;
 651     output_env_t oenv;
 652     t_filenm fnm[] =
 653         {
 654             { efSTX, "-f", NULL, ffREAD },
 655             { efNDX, "-n", NULL, ffOPTRD },
 656             { efSTO, NULL, NULL, ffOPTWR },
 657             { efPQR, "-mead", "mead", ffOPTWR },
 658             { efDAT, "-bf", "bfact", ffOPTRD } };
 659 #define NFILE asize(fnm)
 660
 661     CopyRight(stderr, argv[0]);
 662     parse_common_args(&argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa,
 663                       asize(desc), desc, asize(bugs), bugs, &oenv);
 664
 665     bIndex = opt2bSet("-n", NFILE, fnm) || bNDEF;
 666     bMead = opt2bSet("-mead", NFILE, fnm);
 667     bSetSize = opt2parg_bSet("-box", NPA, pa);
 668     bSetAng = opt2parg_bSet("-angles", NPA, pa);
 669     bSetCenter = opt2parg_bSet("-center", NPA, pa);
 670     bDist = opt2parg_bSet("-d", NPA, pa);
 671         bAlign = opt2parg_bSet("-align", NPA, pa);
 672     /* Only automatically turn on centering without -noc */
 673     if ((bDist || bSetSize || bSetCenter) && !opt2parg_bSet("-c", NPA, pa))
 674     {
 675         bCenter = TRUE;
 676     }
 677     bScale = opt2parg_bSet("-scale", NPA, pa);
 678     bRho = opt2parg_bSet("-density", NPA, pa);
 679     bTranslate = opt2parg_bSet("-translate", NPA, pa);
 680     bRotate = opt2parg_bSet("-rotate", NPA, pa);
 681     if (bScale && bRho)
 682         fprintf(stderr, "WARNING: setting -density overrides -scale\n");
 683     bScale = bScale || bRho;
 684     bCalcGeom = bCenter || bRotate || bOrient || bScale;
 685     bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
 686
 687     infile = ftp2fn(efSTX, NFILE, fnm);
 688     if (bMead)
 689         outfile = ftp2fn(efPQR, NFILE, fnm);
 690     else
 691         outfile = ftp2fn(efSTO, NFILE, fnm);
 692     outftp = fn2ftp(outfile);
 693     inftp = fn2ftp(infile);
 694
 695     aps = gmx_atomprop_init();
 696
 697     if (bMead && bGrasp)
 698     {
 699         printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
 700         bGrasp = FALSE;
 701     }
 702     if (bGrasp && (outftp != efPDB))
 703         gmx_fatal(FARGS, "Output file should be a .pdb file"
 704         " when using the -grasp option\n");
 705         if ((bMead || bGrasp) && !((fn2ftp(infile) == efTPR) ||
 706                 (fn2ftp(infile) == efTPA) ||
 707                 (fn2ftp(infile) == efTPB)))
 708         gmx_fatal(FARGS,"Input file should be a .tp[abr] file"
 709             " when using the -mead option\n");
 710
 711         get_stx_coordnum(infile,&natom);
 712         init_t_atoms(&atoms,natom,TRUE);
 713         snew(x,natom);
 714         snew(v,natom);
 715         read_stx_conf(infile,title,&atoms,x,v,&ePBC,box);
 716         if (fn2ftp(infile) == efPDB)
 717         {
 718             get_pdb_atomnumber(&atoms,aps);
 719         }
 720         printf("Read %d atoms\n",atoms.nr);
 721
 722         /* Get the element numbers if available in a pdb file */
 723         if (fn2ftp(infile) == efPDB)
 724         get_pdb_atomnumber(&atoms,aps);
 725
 726         if (ePBC != epbcNONE)
 727         {
 728             real vol = det(box);
 729             printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
 730                 vol,100*((int)(vol*4.5)));
 731         }
 732
 733         if (bMead || bGrasp || bCONECT)
 734         top = read_top(infile,NULL);
 735
 736         if (bMead || bGrasp)
 737         {
 738             if (atoms.nr != top->atoms.nr)
 739             gmx_fatal(FARGS,"Atom numbers don't match (%d vs. %d)",atoms.nr,top->atoms.nr);
 740         snew(atoms.pdbinfo,top->atoms.nr);
 741         ntype = top->idef.atnr;
 742         for(i=0; (i<atoms.nr); i++) {
 743             /* Determine the Van der Waals radius from the force field */
 744             if (bReadVDW) {
 745                 if (!gmx_atomprop_query(aps,epropVDW,
 746                                         *top->atoms.resinfo[top->atoms.atom[i].resind].name,
 747                                         *top->atoms.atomname[i],&vdw))
 748                     vdw = rvdw;
 749             }
 750             else {
 751                 itype = top->atoms.atom[i].type;
 752                 c12   = top->idef.iparams[itype*ntype+itype].lj.c12;
 753                 c6    = top->idef.iparams[itype*ntype+itype].lj.c6;
 754                 if ((c6 != 0) && (c12 != 0)) {
 755                     real sig6;
 756                     if (bSig56)
 757                         sig6 = 2*c12/c6;
 758                     else
 759                         sig6 = c12/c6;
 760                     vdw   = 0.5*pow(sig6,1.0/6.0);
 761                 }
 762                 else
 763                     vdw = rvdw;
 764             }
 765             /* Factor of 10 for nm -> Angstroms */
 766             vdw *= 10;
 767
 768             if (bMead) {
 769                 atoms.pdbinfo[i].occup = top->atoms.atom[i].q;
 770                 atoms.pdbinfo[i].bfac  = vdw;
 771             }
 772             else {
 773                 atoms.pdbinfo[i].occup = vdw;
 774                 atoms.pdbinfo[i].bfac  = top->atoms.atom[i].q;
 775             }
 776         }
 777     }
 778     bHaveV=FALSE;
 779     for (i=0; (i<natom) && !bHaveV; i++)
 780         for (j=0; (j<DIM) && !bHaveV; j++)
 781             bHaveV=bHaveV || (v[i][j]!=0);
 782     printf("%selocities found\n",bHaveV?"V":"No v");
 783
 784     if (visbox[0] > 0) {
 785         if (bIndex)
 786             gmx_fatal(FARGS,"Sorry, can not visualize box with index groups");
 787         if (outftp != efPDB)
 788             gmx_fatal(FARGS,"Sorry, can only visualize box with a pdb file");
 789     } else if (visbox[0] == -1)
 790         visualize_images("images.pdb",ePBC,box);
 791
 792     /* remove pbc */
 793     if (bRMPBC)
 794         rm_gropbc(&atoms,x,box);
 795
 796     if (bCalcGeom) {
 797         if (bIndex) {
 798             fprintf(stderr,"\nSelect a group for determining the system size:\n");
 799             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 800                       1,&ssize,&sindex,&sgrpname);
 801         } else {
 802             ssize = atoms.nr;
 803             sindex = NULL;
 804         }
 805         diam=calc_geom(ssize,sindex,x,gc,min,max,bCalcDiam);
 806         rvec_sub(max, min, size);
 807         printf("    system size :%7.3f%7.3f%7.3f (nm)\n",
 808                size[XX], size[YY], size[ZZ]);
 809         if (bCalcDiam)
 810             printf("    diameter    :%7.3f               (nm)\n",diam);
 811         printf("    center      :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
 812         printf("    box vectors :%7.3f%7.3f%7.3f (nm)\n",
 813                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
 814         printf("    box angles  :%7.2f%7.2f%7.2f (degrees)\n",
 815                norm2(box[ZZ])==0 ? 0 :
 816         RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
 817         norm2(box[ZZ])==0 ? 0 :
 818         RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
 819         norm2(box[YY])==0 ? 0 :
 820         RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
 821         printf("    box volume  :%7.2f               (nm^3)\n",det(box));
 822     }
 823
 824     if (bRho || bOrient || bAlign)
 825         mass = calc_mass(&atoms,!fn2bTPX(infile),aps);
 826
 827     if (bOrient) {
 828         atom_id *index;
 829         char    *grpnames;
 830
 831         /* Get a group for principal component analysis */
 832         fprintf(stderr,"\nSelect group for the determining the orientation\n");
 833         get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpnames);
 834
 835         /* Orient the principal axes along the coordinate axes */
 836         orient_princ(&atoms,isize,index,natom,x,bHaveV ? v : NULL, NULL);
 837         sfree(index);
 838         sfree(grpnames);
 839     }
 840
 841     if ( bScale ) {
 842         /* scale coordinates and box */
 843         if (bRho) {
 844             /* Compute scaling constant */
 845             real vol,dens;
 846
 847             vol = det(box);
 848             dens = (mass*AMU)/(vol*NANO*NANO*NANO);
 849             fprintf(stderr,"Volume  of input %g (nm^3)\n",vol);
 850             fprintf(stderr,"Mass    of input %g (a.m.u.)\n",mass);
 851             fprintf(stderr,"Density of input %g (g/l)\n",dens);
 852             if (vol==0 || mass==0)
 853                 gmx_fatal(FARGS,"Cannot scale density with "
 854                           "zero mass (%g) or volume (%g)\n",mass,vol);
 855
 856             scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho,1.0/3.0);
 857             fprintf(stderr,"Scaling all box vectors by %g\n",scale[XX]);
 858         }
 859         scale_conf(atoms.nr,x,box,scale);
 860     }
 861
 862         if (bAlign) {
 863                 if (bIndex) {
 864             fprintf(stderr,"\nSelect a group that you want to align:\n");
 865             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 866                       1,&asize,&aindex,&agrpname);
 867         } else {
 868             asize = atoms.nr;
 869             snew(aindex,asize);
 870                         for (i=0;i<asize;i++)
 871                                 aindex[i]=i;
 872         }
 873                 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n",asize,natom,
 874                         targetvec[XX],targetvec[YY],targetvec[ZZ],
 875                         aligncenter[XX],aligncenter[YY],aligncenter[ZZ]);
 876                 /*subtract out pivot point*/
 877                 for(i=0; i<asize; i++)
 878                         rvec_dec(x[aindex[i]],aligncenter);
 879                 /*now determine transform and rotate*/
 880                 /*will this work?*/
 881                 principal_comp(asize,aindex,atoms.atom,x, trans,princd);
 882
 883                 unitv(targetvec,targetvec);
 884                 printf("Using %g %g %g as principal axis\n", trans[0][2],trans[1][2],trans[2][2]);
 885                 tmpvec[XX]=trans[0][2]; tmpvec[YY]=trans[1][2]; tmpvec[ZZ]=trans[2][2];
 886                 calc_rotmatrix(tmpvec, targetvec, rotmatrix);
 887                 /* rotmatrix finished */
 888
 889                 for (i=0;i<asize;++i)
 890                 {
 891                         mvmul(rotmatrix,x[aindex[i]],tmpvec);
 892                         copy_rvec(tmpvec,x[aindex[i]]);
 893                 }
 894
 895                 /*add pivot point back*/
 896                 for(i=0; i<asize; i++)
 897                         rvec_inc(x[aindex[i]],aligncenter);
 898                 if (!bIndex)
 899                         sfree(aindex);
 900         }
 901
 902     if (bTranslate) {
 903         if (bIndex) {
 904             fprintf(stderr,"\nSelect a group that you want to translate:\n");
 905             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 906                       1,&ssize,&sindex,&sgrpname);
 907         } else {
 908             ssize = atoms.nr;
 909             sindex = NULL;
 910         }
 911         printf("Translating %d atoms (out of %d) by %g %g %g nm\n",ssize,natom,
 912                translation[XX],translation[YY],translation[ZZ]);
 913         if (sindex) {
 914             for(i=0; i<ssize; i++)
 915                 rvec_inc(x[sindex[i]],translation);
 916         }
 917         else {
 918             for(i=0; i<natom; i++)
 919                 rvec_inc(x[i],translation);
 920         }
 921     }
 922     if (bRotate) {
 923         /* Rotate */
 924         printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles[XX],rotangles[YY],rotangles[ZZ]);
 925         for(i=0; i<DIM; i++)
 926             rotangles[i] *= DEG2RAD;
 927         rotate_conf(natom,x,v,rotangles[XX],rotangles[YY],rotangles[ZZ]);
 928     }
 929
 930     if (bCalcGeom) {
 931         /* recalc geometrical center and max and min coordinates and size */
 932         calc_geom(ssize,sindex,x,gc,min,max,FALSE);
 933         rvec_sub(max, min, size);
 934         if (bScale || bOrient || bRotate)
 935             printf("new system size : %6.3f %6.3f %6.3f\n",
 936                    size[XX],size[YY],size[ZZ]);
 937     }
 938
 939     if (bSetSize || bDist || (btype[0][0]=='t' && bSetAng)) {
 940         ePBC = epbcXYZ;
 941         if (!(bSetSize || bDist))
 942             for (i=0; i<DIM; i++)
 943                 newbox[i] = norm(box[i]);
 944         clear_mat(box);
 945         /* calculate new boxsize */
 946         switch(btype[0][0]){
 947         case 't':
 948             if (bDist)
 949                 for(i=0; i<DIM; i++)
 950                     newbox[i] = size[i]+2*dist;
 951             if (!bSetAng) {
 952                 box[XX][XX] = newbox[XX];
 953                 box[YY][YY] = newbox[YY];
 954                 box[ZZ][ZZ] = newbox[ZZ];
 955             } else {
 956                 matrix_convert(box,newbox,newang);
 957             }
 958             break;
 959         case 'c':
 960         case 'd':
 961         case 'o':
 962             if (bSetSize)
 963                 d = newbox[0];
 964             else
 965                 d = diam+2*dist;
 966             if (btype[0][0] == 'c')
 967                 for(i=0; i<DIM; i++)
 968                     box[i][i] = d;
 969             else if (btype[0][0] == 'd') {
 970                 box[XX][XX] = d;
 971                 box[YY][YY] = d;
 972                 box[ZZ][XX] = d/2;
 973                 box[ZZ][YY] = d/2;
 974                 box[ZZ][ZZ] = d*sqrt(2)/2;
 975             } else {
 976                 box[XX][XX] = d;
 977                 box[YY][XX] = d/3;
 978                 box[YY][YY] = d*sqrt(2)*2/3;
 979                 box[ZZ][XX] = -d/3;
 980                 box[ZZ][YY] = d*sqrt(2)/3;
 981                 box[ZZ][ZZ] = d*sqrt(6)/3;
 982             }
 983             break;
 984         }
 985     }
 986
 987     /* calculate new coords for geometrical center */
 988     if (!bSetCenter)
 989         calc_box_center(ecenterDEF,box,center);
 990
 991     /* center molecule on 'center' */
 992     if (bCenter)
 993         center_conf(natom,x,center,gc);
 994
 995     /* print some */
 996     if (bCalcGeom) {
 997         calc_geom(ssize,sindex,x, gc, min, max, FALSE);
 998         printf("new center      :%7.3f%7.3f%7.3f (nm)\n",gc[XX],gc[YY],gc[ZZ]);
 999     }
1000     if (bOrient || bScale || bDist || bSetSize) {
1001         printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1002                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
1003         printf("new box angles  :%7.2f%7.2f%7.2f (degrees)\n",
1004                norm2(box[ZZ])==0 ? 0 :
1005         RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
1006         norm2(box[ZZ])==0 ? 0 :
1007         RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
1008         norm2(box[YY])==0 ? 0 :
1009         RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
1010         printf("new box volume  :%7.2f               (nm^3)\n",det(box));
1011     }
1012
1013     if (check_box(epbcXYZ,box))
1014         printf("\nWARNING: %s\n",check_box(epbcXYZ,box));
1015
1016     if (bDist && btype[0][0]=='t')
1017     {
1018         if(TRICLINIC(box))
1019         {
1020             printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
1021                 "distance from the solute to a box surface along the corresponding normal\n"
1022                 "vector might be somewhat smaller than your specified value %f.\n"
1023                 "You can check the actual value with g_mindist -pi\n",dist);
1024         }
1025         else
1026         {
1027             printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
1028                 "If the molecule rotates the actual distance will be smaller. You might want\n"
1029                 "to use a cubic box instead, or why not try a dodecahedron today?\n");
1030         }
1031     }
1032     if (bCONECT && (outftp == efPDB) && (inftp == efTPR))
1033         conect = gmx_conect_generate(top);
1034     else
1035         conect = NULL;
1036
1037     if (bIndex) {
1038         fprintf(stderr,"\nSelect a group for output:\n");
1039         get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
1040                   1,&isize,&index,&grpname);
1041         if (opt2bSet("-bf",NFILE,fnm))
1042             gmx_fatal(FARGS,"combination not implemented: -bf -n  or -bf -ndef");
1043         else {
1044             if (outftp == efPDB) {
1045                 out=ffopen(outfile,"w");
1046                 write_pdbfile_indexed(out,title,&atoms,x,ePBC,box,' ',1,isize,index,conect,TRUE);
1047                 ffclose(out);
1048             }
1049             else
1050                 write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box,
1051                     isize,index);
1052         }
1053     }
1054     else {
1055         if ((outftp == efPDB) || (outftp == efPQR)) {
1056             out=ffopen(outfile,"w");
1057             if (bMead) {
1058                 set_pdb_wide_format(TRUE);
1059                 fprintf(out,"REMARK    "
1060                         "The B-factors in this file hold atomic radii\n");
1061                 fprintf(out,"REMARK    "
1062                         "The occupancy in this file hold atomic charges\n");
1063             }
1064             else if (bGrasp) {
1065                 fprintf(out,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
1066                 fprintf(out,"REMARK    "
1067                         "The B-factors in this file hold atomic charges\n");
1068                 fprintf(out,"REMARK    "
1069                         "The occupancy in this file hold atomic radii\n");
1070             }
1071             else if (opt2bSet("-bf",NFILE,fnm)) {
1072                 read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
1073                 set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,
1074                                   n_bfac,bfac,bfac_nr,peratom);
1075             }
1076             if (opt2parg_bSet("-label",NPA,pa)) {
1077                 for(i=0; (i<atoms.nr); i++)
1078                     atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
1079             }
1080             write_pdbfile(out,title,&atoms,x,ePBC,box,' ',-1,conect,TRUE);
1081             if (bLegend)
1082                 pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
1083             if (visbox[0] > 0)
1084                 visualize_box(out,bLegend ? atoms.nr+12 : atoms.nr,
1085                     bLegend? atoms.nres=12 : atoms.nres,box,visbox);
1086             ffclose(out);
1087         }
1088         else
1089             write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box);
1090     }
1091     gmx_atomprop_destroy(aps);
1092
1093     do_view(oenv,outfile,NULL);
1094
1095     thanx(stderr);
1096
1097     return 0;
1098 }