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