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