src/gromacs/gmxana/gmx_editconf.cpp

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