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