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, gmx_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, gmx_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, gmx_bool peratom)
 211 {
 212     FILE *out;
 213     real bfac_min, bfac_max;
 214     int i, n;
 215     gmx_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 static void renum_resnr(t_atoms *atoms,int isize,const int *index,
 452                         int resnr_start)
 453 {
 454     int i,resind_prev,resind;
 455
 456     resind_prev = -1;
 457     for(i=0; i<isize; i++)
 458     {
 459         resind = atoms->atom[index == NULL ? i : index[i]].resind;
 460         if (resind != resind_prev)
 461         {
 462             atoms->resinfo[resind].nr = resnr_start;
 463             resnr_start++;
 464         }
 465         resind_prev = resind;
 466     }
 467 }
 468
 469 int gmx_editconf(int argc, char *argv[])
 470 {
 471     const char
 472         *desc[] =
 473             {
 474                 "[TT]editconf[tt] converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
 475                 "or [TT].pdb[tt].",
 476                 "[PAR]",
 477                 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
 478                 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
 479                 "will center the system in the box, unless [TT]-noc[tt] is used.",
 480                 "[PAR]",
 481                 "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
 482                 "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
 483                 "[TT]dodecahedron[tt] represents a rhombic dodecahedron and",
 484                 "[TT]octahedron[tt] is a truncated octahedron.",
 485                 "The last two are special cases of a triclinic box.",
 486                 "The length of the three box vectors of the truncated octahedron is the",
 487                 "shortest distance between two opposite hexagons.",
 488                 "The volume of a dodecahedron is 0.71 and that of a truncated octahedron",
 489                 "is 0.77 of that of a cubic box with the same periodic image distance.",
 490                 "[PAR]",
 491                 "Option [TT]-box[tt] requires only",
 492                 "one value for a cubic box, dodecahedron and a truncated octahedron.",
 493                 "[PAR]",
 494                 "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the x, y",
 495                 "and z directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
 496                 "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
 497                 "to the diameter of the system (largest distance between atoms) plus twice",
 498                 "the specified distance.",
 499                 "[PAR]",
 500                 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
 501                 "a triclinic box and can not be used with option [TT]-d[tt].",
 502                 "[PAR]",
 503                 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
 504                 "can be selected for calculating the size and the geometric center,",
 505                 "otherwise the whole system is used.",
 506                 "[PAR]",
 507                 "[TT]-rotate[tt] rotates the coordinates and velocities.",
 508                 "[PAR]",
 509                 "[TT]-princ[tt] aligns the principal axes of the system along the",
 510                 "coordinate axes, with the longest axis aligned with the x axis. ",
 511                 "This may allow you to decrease the box volume,",
 512                 "but beware that molecules can rotate significantly in a nanosecond.",
 513                 "[PAR]",
 514                 "Scaling is applied before any of the other operations are",
 515                 "performed. Boxes and coordinates can be scaled to give a certain density (option",
 516                 "[TT]-density[tt]). Note that this may be inaccurate in case a gro",
 517                 "file is given as input. A special feature of the scaling option, when the",
 518                 "factor -1 is given in one dimension, one obtains a mirror image,",
 519                 "mirrored in one of the planes. When one uses -1 in three dimensions, ",
 520                 "a point-mirror image is obtained.[PAR]",
 521                 "Groups are selected after all operations have been applied.[PAR]",
 522                 "Periodicity can be removed in a crude manner.",
 523                 "It is important that the box vectors at the bottom of your input file",
 524                 "are correct when the periodicity is to be removed.",
 525                 "[PAR]",
 526                 "When writing [TT].pdb[tt] files, B-factors can be",
 527                 "added with the [TT]-bf[tt] option. B-factors are read",
 528                 "from a file with with following format: first line states number of",
 529                 "entries in the file, next lines state an index",
 530                 "followed by a B-factor. The B-factors will be attached per residue",
 531                 "unless an index is larger than the number of residues or unless the",
 532                 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
 533                 "be added instead of B-factors. [TT]-legend[tt] will produce",
 534                 "a row of CA atoms with B-factors ranging from the minimum to the",
 535                 "maximum value found, effectively making a legend for viewing.",
 536                 "[PAR]",
 537                 "With the option [TT]-mead[tt] a special [TT].pdb[tt] ([TT].pqr[tt])",
 538                 "file for the MEAD electrostatics",
 539                 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
 540                 "is that the input file is a run input file.",
 541                 "The B-factor field is then filled with the Van der Waals radius",
 542                 "of the atoms while the occupancy field will hold the charge.",
 543                 "[PAR]",
 544                 "The option [TT]-grasp[tt] is similar, but it puts the charges in the B-factor",
 545                 "and the radius in the occupancy.",
 546                 "[PAR]",
 547                 "Option [TT]-align[tt] allows alignment",
 548                 "of the principal axis of a specified group against the given vector, ",
 549                                 "with an optional center of rotation specified by [TT]-aligncenter[tt].",
 550                 "[PAR]",
 551                 "Finally with option [TT]-label[tt], [TT]editconf[tt] can add a chain identifier",
 552                 "to a [TT].pdb[tt] file, which can be useful for analysis with e.g. rasmol.",
 553                 "[PAR]",
 554                 "To convert a truncated octrahedron file produced by a package which uses",
 555                 "a cubic box with the corners cut off (such as GROMOS), use:[BR]",
 556                 "[TT]editconf -f in -rotate 0 45 35.264 -bt o -box veclen -o out[tt][BR]",
 557                 "where [TT]veclen[tt] is the size of the cubic box times sqrt(3)/2." };
 558     const char *bugs[] =
 559         {
 560             "For complex molecules, the periodicity removal routine may break down, "
 561                 "in that case you can use [TT]trjconv[tt]." };
 562     static real dist = 0.0, rbox = 0.0, to_diam = 0.0;
 563     static gmx_bool bNDEF = FALSE, bRMPBC = FALSE, bCenter = FALSE, bReadVDW =
 564         FALSE, bCONECT = FALSE;
 565     static gmx_bool peratom = FALSE, bLegend = FALSE, bOrient = FALSE, bMead =
 566         FALSE, bGrasp = FALSE, bSig56 = FALSE;
 567     static rvec scale =
 568         { 1, 1, 1 }, newbox =
 569         { 0, 0, 0 }, newang =
 570         { 90, 90, 90 };
 571     static real rho = 1000.0, rvdw = 0.12;
 572     static rvec center =
 573         { 0, 0, 0 }, translation =
 574         { 0, 0, 0 }, rotangles =
 575         { 0, 0, 0 }, aligncenter =
 576                 { 0, 0, 0 }, targetvec =
 577         { 0, 0, 0 };
 578     static const char *btype[] =
 579         { NULL, "triclinic", "cubic", "dodecahedron", "octahedron", NULL },
 580         *label = "A";
 581     static rvec visbox =
 582         { 0, 0, 0 };
 583     static int resnr_start = -1;
 584     t_pargs
 585         pa[] =
 586             {
 587                     { "-ndef", FALSE, etBOOL,
 588                         { &bNDEF }, "Choose output from default index groups" },
 589                     { "-visbox", FALSE, etRVEC,
 590                         { visbox },
 591                         "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
 592                     { "-bt", FALSE, etENUM,
 593                         { btype }, "Box type for -box and -d" },
 594                     { "-box", FALSE, etRVEC,
 595                         { newbox }, "Box vector lengths (a,b,c)" },
 596                     { "-angles", FALSE, etRVEC,
 597                         { newang }, "Angles between the box vectors (bc,ac,ab)" },
 598                     { "-d", FALSE, etREAL,
 599                         { &dist }, "Distance between the solute and the box" },
 600                     { "-c", FALSE, etBOOL,
 601                         { &bCenter },
 602                         "Center molecule in box (implied by -box and -d)" },
 603                     { "-center", FALSE, etRVEC,
 604                         { center }, "Coordinates of geometrical center" },
 605                     { "-aligncenter", FALSE, etRVEC,
 606                         { aligncenter }, "Center of rotation for alignment" },
 607                     { "-align", FALSE, etRVEC,
 608                         { targetvec },
 609                         "Align to target vector" },
 610                     { "-translate", FALSE, etRVEC,
 611                         { translation }, "Translation" },
 612                     { "-rotate", FALSE, etRVEC,
 613                         { rotangles },
 614                         "Rotation around the X, Y and Z axes in degrees" },
 615                     { "-princ", FALSE, etBOOL,
 616                         { &bOrient },
 617                         "Orient molecule(s) along their principal axes" },
 618                     { "-scale", FALSE, etRVEC,
 619                         { scale }, "Scaling factor" },
 620                     { "-density", FALSE, etREAL,
 621                         { &rho },
 622                         "Density (g/L) of the output box achieved by scaling" },
 623                     { "-pbc", FALSE, etBOOL,
 624                         { &bRMPBC },
 625                         "Remove the periodicity (make molecule whole again)" },
 626                     { "-resnr", FALSE, etINT,
 627                         { &resnr_start },
 628                         " Renumber residues starting from resnr" },
 629                     { "-grasp", FALSE, etBOOL,
 630                         { &bGrasp },
 631                         "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
 632                     {
 633                         "-rvdw", FALSE, etREAL,
 634                          { &rvdw },
 635                         "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" },
 636                     { "-sig56", FALSE, etREAL,
 637                         { &bSig56 },
 638                         "Use rmin/2 (minimum in the Van der Waals potential) rather than sigma/2 " },
 639                     {
 640                         "-vdwread", FALSE, etBOOL,
 641                         { &bReadVDW },
 642                         "Read the Van der Waals radii from the file vdwradii.dat rather than computing the radii based on the force field" },
 643                     { "-atom", FALSE, etBOOL,
 644                         { &peratom }, "Force B-factor attachment per atom" },
 645                     { "-legend", FALSE, etBOOL,
 646                         { &bLegend }, "Make B-factor legend" },
 647                     { "-label", FALSE, etSTR,
 648                         { &label }, "Add chain label for all residues" },
 649                     {
 650                         "-conect", FALSE, etBOOL,
 651                         { &bCONECT },
 652                         "Add CONECT records to a [TT].pdb[tt] file when written. Can only be done when a topology is present" } };
 653 #define NPA asize(pa)
 654
 655     FILE *out;
 656     const char *infile, *outfile;
 657     char title[STRLEN];
 658     int outftp, inftp, natom, i, j, n_bfac, itype, ntype;
 659     double *bfac = NULL, c6, c12;
 660     int *bfac_nr = NULL;
 661     t_topology *top = NULL;
 662     t_atoms atoms;
 663     char *grpname, *sgrpname, *agrpname;
 664     int isize, ssize, tsize, asize;
 665     atom_id *index, *sindex, *tindex, *aindex;
 666     rvec *x, *v, gc, min, max, size;
 667     int ePBC;
 668     matrix box,rotmatrix,trans;
 669         rvec princd,tmpvec;
 670     gmx_bool bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
 671     gmx_bool bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
 672     real xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
 673     gmx_atomprop_t aps;
 674     gmx_conect conect;
 675     output_env_t oenv;
 676     t_filenm fnm[] =
 677         {
 678             { efSTX, "-f", NULL, ffREAD },
 679             { efNDX, "-n", NULL, ffOPTRD },
 680             { efSTO, NULL, NULL, ffOPTWR },
 681             { efPQR, "-mead", "mead", ffOPTWR },
 682             { efDAT, "-bf", "bfact", ffOPTRD } };
 683 #define NFILE asize(fnm)
 684
 685     CopyRight(stderr, argv[0]);
 686     parse_common_args(&argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa,
 687                       asize(desc), desc, asize(bugs), bugs, &oenv);
 688
 689     bIndex = opt2bSet("-n", NFILE, fnm) || bNDEF;
 690     bMead = opt2bSet("-mead", NFILE, fnm);
 691     bSetSize = opt2parg_bSet("-box", NPA, pa);
 692     bSetAng = opt2parg_bSet("-angles", NPA, pa);
 693     bSetCenter = opt2parg_bSet("-center", NPA, pa);
 694     bDist = opt2parg_bSet("-d", NPA, pa);
 695         bAlign = opt2parg_bSet("-align", NPA, pa);
 696     /* Only automatically turn on centering without -noc */
 697     if ((bDist || bSetSize || bSetCenter) && !opt2parg_bSet("-c", NPA, pa))
 698     {
 699         bCenter = TRUE;
 700     }
 701     bScale = opt2parg_bSet("-scale", NPA, pa);
 702     bRho = opt2parg_bSet("-density", NPA, pa);
 703     bTranslate = opt2parg_bSet("-translate", NPA, pa);
 704     bRotate = opt2parg_bSet("-rotate", NPA, pa);
 705     if (bScale && bRho)
 706         fprintf(stderr, "WARNING: setting -density overrides -scale\n");
 707     bScale = bScale || bRho;
 708     bCalcGeom = bCenter || bRotate || bOrient || bScale;
 709     bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
 710
 711     infile = ftp2fn(efSTX, NFILE, fnm);
 712     if (bMead)
 713         outfile = ftp2fn(efPQR, NFILE, fnm);
 714     else
 715         outfile = ftp2fn(efSTO, NFILE, fnm);
 716     outftp = fn2ftp(outfile);
 717     inftp = fn2ftp(infile);
 718
 719     aps = gmx_atomprop_init();
 720
 721     if (bMead && bGrasp)
 722     {
 723         printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
 724         bGrasp = FALSE;
 725     }
 726     if (bGrasp && (outftp != efPDB))
 727         gmx_fatal(FARGS, "Output file should be a .pdb file"
 728         " when using the -grasp option\n");
 729         if ((bMead || bGrasp) && !((fn2ftp(infile) == efTPR) ||
 730                 (fn2ftp(infile) == efTPA) ||
 731                 (fn2ftp(infile) == efTPB)))
 732         gmx_fatal(FARGS,"Input file should be a .tp[abr] file"
 733             " when using the -mead option\n");
 734
 735         get_stx_coordnum(infile,&natom);
 736         init_t_atoms(&atoms,natom,TRUE);
 737         snew(x,natom);
 738         snew(v,natom);
 739         read_stx_conf(infile,title,&atoms,x,v,&ePBC,box);
 740         if (fn2ftp(infile) == efPDB)
 741         {
 742             get_pdb_atomnumber(&atoms,aps);
 743         }
 744         printf("Read %d atoms\n",atoms.nr);
 745
 746         /* Get the element numbers if available in a pdb file */
 747         if (fn2ftp(infile) == efPDB)
 748         get_pdb_atomnumber(&atoms,aps);
 749
 750         if (ePBC != epbcNONE)
 751         {
 752             real vol = det(box);
 753             printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
 754                 vol,100*((int)(vol*4.5)));
 755         }
 756
 757         if (bMead || bGrasp || bCONECT)
 758         top = read_top(infile,NULL);
 759
 760         if (bMead || bGrasp)
 761         {
 762             if (atoms.nr != top->atoms.nr)
 763             gmx_fatal(FARGS,"Atom numbers don't match (%d vs. %d)",atoms.nr,top->atoms.nr);
 764         snew(atoms.pdbinfo,top->atoms.nr);
 765         ntype = top->idef.atnr;
 766         for(i=0; (i<atoms.nr); i++) {
 767             /* Determine the Van der Waals radius from the force field */
 768             if (bReadVDW) {
 769                 if (!gmx_atomprop_query(aps,epropVDW,
 770                                         *top->atoms.resinfo[top->atoms.atom[i].resind].name,
 771                                         *top->atoms.atomname[i],&vdw))
 772                     vdw = rvdw;
 773             }
 774             else {
 775                 itype = top->atoms.atom[i].type;
 776                 c12   = top->idef.iparams[itype*ntype+itype].lj.c12;
 777                 c6    = top->idef.iparams[itype*ntype+itype].lj.c6;
 778                 if ((c6 != 0) && (c12 != 0)) {
 779                     real sig6;
 780                     if (bSig56)
 781                         sig6 = 2*c12/c6;
 782                     else
 783                         sig6 = c12/c6;
 784                     vdw   = 0.5*pow(sig6,1.0/6.0);
 785                 }
 786                 else
 787                     vdw = rvdw;
 788             }
 789             /* Factor of 10 for nm -> Angstroms */
 790             vdw *= 10;
 791
 792             if (bMead) {
 793                 atoms.pdbinfo[i].occup = top->atoms.atom[i].q;
 794                 atoms.pdbinfo[i].bfac  = vdw;
 795             }
 796             else {
 797                 atoms.pdbinfo[i].occup = vdw;
 798                 atoms.pdbinfo[i].bfac  = top->atoms.atom[i].q;
 799             }
 800         }
 801     }
 802     bHaveV=FALSE;
 803     for (i=0; (i<natom) && !bHaveV; i++)
 804         for (j=0; (j<DIM) && !bHaveV; j++)
 805             bHaveV=bHaveV || (v[i][j]!=0);
 806     printf("%selocities found\n",bHaveV?"V":"No v");
 807
 808     if (visbox[0] > 0) {
 809         if (bIndex)
 810             gmx_fatal(FARGS,"Sorry, can not visualize box with index groups");
 811         if (outftp != efPDB)
 812             gmx_fatal(FARGS,"Sorry, can only visualize box with a pdb file");
 813     } else if (visbox[0] == -1)
 814         visualize_images("images.pdb",ePBC,box);
 815
 816     /* remove pbc */
 817     if (bRMPBC)
 818         rm_gropbc(&atoms,x,box);
 819
 820     if (bCalcGeom) {
 821         if (bIndex) {
 822             fprintf(stderr,"\nSelect a group for determining the system size:\n");
 823             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 824                       1,&ssize,&sindex,&sgrpname);
 825         } else {
 826             ssize = atoms.nr;
 827             sindex = NULL;
 828         }
 829         diam=calc_geom(ssize,sindex,x,gc,min,max,bCalcDiam);
 830         rvec_sub(max, min, size);
 831         printf("    system size :%7.3f%7.3f%7.3f (nm)\n",
 832                size[XX], size[YY], size[ZZ]);
 833         if (bCalcDiam)
 834             printf("    diameter    :%7.3f               (nm)\n",diam);
 835         printf("    center      :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
 836         printf("    box vectors :%7.3f%7.3f%7.3f (nm)\n",
 837                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
 838         printf("    box angles  :%7.2f%7.2f%7.2f (degrees)\n",
 839                norm2(box[ZZ])==0 ? 0 :
 840         RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
 841         norm2(box[ZZ])==0 ? 0 :
 842         RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
 843         norm2(box[YY])==0 ? 0 :
 844         RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
 845         printf("    box volume  :%7.2f               (nm^3)\n",det(box));
 846     }
 847
 848     if (bRho || bOrient || bAlign)
 849         mass = calc_mass(&atoms,!fn2bTPX(infile),aps);
 850
 851     if (bOrient) {
 852         atom_id *index;
 853         char    *grpnames;
 854
 855         /* Get a group for principal component analysis */
 856         fprintf(stderr,"\nSelect group for the determining the orientation\n");
 857         get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpnames);
 858
 859         /* Orient the principal axes along the coordinate axes */
 860         orient_princ(&atoms,isize,index,natom,x,bHaveV ? v : NULL, NULL);
 861         sfree(index);
 862         sfree(grpnames);
 863     }
 864
 865     if ( bScale ) {
 866         /* scale coordinates and box */
 867         if (bRho) {
 868             /* Compute scaling constant */
 869             real vol,dens;
 870
 871             vol = det(box);
 872             dens = (mass*AMU)/(vol*NANO*NANO*NANO);
 873             fprintf(stderr,"Volume  of input %g (nm^3)\n",vol);
 874             fprintf(stderr,"Mass    of input %g (a.m.u.)\n",mass);
 875             fprintf(stderr,"Density of input %g (g/l)\n",dens);
 876             if (vol==0 || mass==0)
 877                 gmx_fatal(FARGS,"Cannot scale density with "
 878                           "zero mass (%g) or volume (%g)\n",mass,vol);
 879
 880             scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho,1.0/3.0);
 881             fprintf(stderr,"Scaling all box vectors by %g\n",scale[XX]);
 882         }
 883         scale_conf(atoms.nr,x,box,scale);
 884     }
 885
 886         if (bAlign) {
 887                 if (bIndex) {
 888             fprintf(stderr,"\nSelect a group that you want to align:\n");
 889             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 890                       1,&asize,&aindex,&agrpname);
 891         } else {
 892             asize = atoms.nr;
 893             snew(aindex,asize);
 894                         for (i=0;i<asize;i++)
 895                                 aindex[i]=i;
 896         }
 897                 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n",asize,natom,
 898                         targetvec[XX],targetvec[YY],targetvec[ZZ],
 899                         aligncenter[XX],aligncenter[YY],aligncenter[ZZ]);
 900                 /*subtract out pivot point*/
 901                 for(i=0; i<asize; i++)
 902                         rvec_dec(x[aindex[i]],aligncenter);
 903                 /*now determine transform and rotate*/
 904                 /*will this work?*/
 905                 principal_comp(asize,aindex,atoms.atom,x, trans,princd);
 906
 907                 unitv(targetvec,targetvec);
 908                 printf("Using %g %g %g as principal axis\n", trans[0][2],trans[1][2],trans[2][2]);
 909                 tmpvec[XX]=trans[0][2]; tmpvec[YY]=trans[1][2]; tmpvec[ZZ]=trans[2][2];
 910                 calc_rotmatrix(tmpvec, targetvec, rotmatrix);
 911                 /* rotmatrix finished */
 912
 913                 for (i=0;i<asize;++i)
 914                 {
 915                         mvmul(rotmatrix,x[aindex[i]],tmpvec);
 916                         copy_rvec(tmpvec,x[aindex[i]]);
 917                 }
 918
 919                 /*add pivot point back*/
 920                 for(i=0; i<asize; i++)
 921                         rvec_inc(x[aindex[i]],aligncenter);
 922                 if (!bIndex)
 923                         sfree(aindex);
 924         }
 925
 926     if (bTranslate) {
 927         if (bIndex) {
 928             fprintf(stderr,"\nSelect a group that you want to translate:\n");
 929             get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
 930                       1,&ssize,&sindex,&sgrpname);
 931         } else {
 932             ssize = atoms.nr;
 933             sindex = NULL;
 934         }
 935         printf("Translating %d atoms (out of %d) by %g %g %g nm\n",ssize,natom,
 936                translation[XX],translation[YY],translation[ZZ]);
 937         if (sindex) {
 938             for(i=0; i<ssize; i++)
 939                 rvec_inc(x[sindex[i]],translation);
 940         }
 941         else {
 942             for(i=0; i<natom; i++)
 943                 rvec_inc(x[i],translation);
 944         }
 945     }
 946     if (bRotate) {
 947         /* Rotate */
 948         printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles[XX],rotangles[YY],rotangles[ZZ]);
 949         for(i=0; i<DIM; i++)
 950             rotangles[i] *= DEG2RAD;
 951         rotate_conf(natom,x,v,rotangles[XX],rotangles[YY],rotangles[ZZ]);
 952     }
 953
 954     if (bCalcGeom) {
 955         /* recalc geometrical center and max and min coordinates and size */
 956         calc_geom(ssize,sindex,x,gc,min,max,FALSE);
 957         rvec_sub(max, min, size);
 958         if (bScale || bOrient || bRotate)
 959             printf("new system size : %6.3f %6.3f %6.3f\n",
 960                    size[XX],size[YY],size[ZZ]);
 961     }
 962
 963     if (bSetSize || bDist || (btype[0][0]=='t' && bSetAng)) {
 964         ePBC = epbcXYZ;
 965         if (!(bSetSize || bDist))
 966             for (i=0; i<DIM; i++)
 967                 newbox[i] = norm(box[i]);
 968         clear_mat(box);
 969         /* calculate new boxsize */
 970         switch(btype[0][0]){
 971         case 't':
 972             if (bDist)
 973                 for(i=0; i<DIM; i++)
 974                     newbox[i] = size[i]+2*dist;
 975             if (!bSetAng) {
 976                 box[XX][XX] = newbox[XX];
 977                 box[YY][YY] = newbox[YY];
 978                 box[ZZ][ZZ] = newbox[ZZ];
 979             } else {
 980                 matrix_convert(box,newbox,newang);
 981             }
 982             break;
 983         case 'c':
 984         case 'd':
 985         case 'o':
 986             if (bSetSize)
 987                 d = newbox[0];
 988             else
 989                 d = diam+2*dist;
 990             if (btype[0][0] == 'c')
 991                 for(i=0; i<DIM; i++)
 992                     box[i][i] = d;
 993             else if (btype[0][0] == 'd') {
 994                 box[XX][XX] = d;
 995                 box[YY][YY] = d;
 996                 box[ZZ][XX] = d/2;
 997                 box[ZZ][YY] = d/2;
 998                 box[ZZ][ZZ] = d*sqrt(2)/2;
 999             } else {
1000                 box[XX][XX] = d;
1001                 box[YY][XX] = d/3;
1002                 box[YY][YY] = d*sqrt(2)*2/3;
1003                 box[ZZ][XX] = -d/3;
1004                 box[ZZ][YY] = d*sqrt(2)/3;
1005                 box[ZZ][ZZ] = d*sqrt(6)/3;
1006             }
1007             break;
1008         }
1009     }
1010
1011     /* calculate new coords for geometrical center */
1012     if (!bSetCenter)
1013         calc_box_center(ecenterDEF,box,center);
1014
1015     /* center molecule on 'center' */
1016     if (bCenter)
1017         center_conf(natom,x,center,gc);
1018
1019     /* print some */
1020     if (bCalcGeom) {
1021         calc_geom(ssize,sindex,x, gc, min, max, FALSE);
1022         printf("new center      :%7.3f%7.3f%7.3f (nm)\n",gc[XX],gc[YY],gc[ZZ]);
1023     }
1024     if (bOrient || bScale || bDist || bSetSize) {
1025         printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1026                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
1027         printf("new box angles  :%7.2f%7.2f%7.2f (degrees)\n",
1028                norm2(box[ZZ])==0 ? 0 :
1029         RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
1030         norm2(box[ZZ])==0 ? 0 :
1031         RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
1032         norm2(box[YY])==0 ? 0 :
1033         RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
1034         printf("new box volume  :%7.2f               (nm^3)\n",det(box));
1035     }
1036
1037     if (check_box(epbcXYZ,box))
1038         printf("\nWARNING: %s\n",check_box(epbcXYZ,box));
1039
1040     if (bDist && btype[0][0]=='t')
1041     {
1042         if(TRICLINIC(box))
1043         {
1044             printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
1045                 "distance from the solute to a box surface along the corresponding normal\n"
1046                 "vector might be somewhat smaller than your specified value %f.\n"
1047                 "You can check the actual value with g_mindist -pi\n",dist);
1048         }
1049         else if (!opt2parg_bSet("-bt", NPA, pa))
1050         {
1051             printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
1052                 "If the molecule rotates the actual distance will be smaller. You might want\n"
1053                 "to use a cubic box instead, or why not try a dodecahedron today?\n");
1054         }
1055     }
1056     if (bCONECT && (outftp == efPDB) && (inftp == efTPR))
1057         conect = gmx_conect_generate(top);
1058     else
1059         conect = NULL;
1060
1061     if (bIndex) {
1062         fprintf(stderr,"\nSelect a group for output:\n");
1063         get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
1064                   1,&isize,&index,&grpname);
1065
1066         if (resnr_start >= 0)
1067         {
1068             renum_resnr(&atoms,isize,index,resnr_start);
1069         }
1070
1071         if (opt2parg_bSet("-label",NPA,pa)) {
1072             for(i=0; (i<atoms.nr); i++)
1073                 atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
1074         }
1075
1076         if (opt2bSet("-bf",NFILE,fnm) || bLegend)
1077         {
1078             gmx_fatal(FARGS,"Sorry, cannot do bfactors with an index group.");
1079         }
1080
1081         if (outftp == efPDB)
1082         {
1083             out=ffopen(outfile,"w");
1084             write_pdbfile_indexed(out,title,&atoms,x,ePBC,box,' ',1,isize,index,conect,TRUE);
1085             ffclose(out);
1086         }
1087         else
1088         {
1089             write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box,isize,index);
1090         }
1091     }
1092     else
1093     {
1094         if (resnr_start >= 0)
1095         {
1096             renum_resnr(&atoms,atoms.nr,NULL,resnr_start);
1097         }
1098
1099         if ((outftp == efPDB) || (outftp == efPQR)) {
1100             out=ffopen(outfile,"w");
1101             if (bMead) {
1102                 set_pdb_wide_format(TRUE);
1103                 fprintf(out,"REMARK    "
1104                         "The B-factors in this file hold atomic radii\n");
1105                 fprintf(out,"REMARK    "
1106                         "The occupancy in this file hold atomic charges\n");
1107             }
1108             else if (bGrasp) {
1109                 fprintf(out,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
1110                 fprintf(out,"REMARK    "
1111                         "The B-factors in this file hold atomic charges\n");
1112                 fprintf(out,"REMARK    "
1113                         "The occupancy in this file hold atomic radii\n");
1114             }
1115             else if (opt2bSet("-bf",NFILE,fnm)) {
1116                 read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
1117                 set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,
1118                                   n_bfac,bfac,bfac_nr,peratom);
1119             }
1120             if (opt2parg_bSet("-label",NPA,pa)) {
1121                 for(i=0; (i<atoms.nr); i++)
1122                     atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
1123             }
1124             write_pdbfile(out,title,&atoms,x,ePBC,box,' ',-1,conect,TRUE);
1125             if (bLegend)
1126                 pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
1127             if (visbox[0] > 0)
1128                 visualize_box(out,bLegend ? atoms.nr+12 : atoms.nr,
1129                     bLegend? atoms.nres=12 : atoms.nres,box,visbox);
1130             ffclose(out);
1131         }
1132         else
1133             write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box);
1134     }
1135     gmx_atomprop_destroy(aps);
1136
1137     do_view(oenv,outfile,NULL);
1138
1139     thanx(stderr);
1140
1141     return 0;
1142 }