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