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