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