src/gromacs/gmxana/gmx_editconf.c

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