src/gromacs/tools/check.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-2013, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "check.h"
  41
  42 #include <cmath>
  43 #include <cstdio>
  44 #include <cstring>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/enxio.h"
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "gromacs/fileio/tpxio.h"
  51 #include "gromacs/fileio/trxio.h"
  52 #include "gromacs/fileio/xtcio.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/mdrun/mdmodules.h"
  57 #include "gromacs/mdtypes/inputrec.h"
  58 #include "gromacs/mdtypes/md_enums.h"
  59 #include "gromacs/mdtypes/state.h"
  60 #include "gromacs/pbcutil/pbc.h"
  61 #include "gromacs/topology/atomprop.h"
  62 #include "gromacs/topology/block.h"
  63 #include "gromacs/topology/ifunc.h"
  64 #include "gromacs/topology/index.h"
  65 #include "gromacs/topology/mtop_util.h"
  66 #include "gromacs/topology/topology.h"
  67 #include "gromacs/trajectory/energyframe.h"
  68 #include "gromacs/trajectory/trajectoryframe.h"
  69 #include "gromacs/utility/arraysize.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/futil.h"
  72 #include "gromacs/utility/smalloc.h"
  73
  74 typedef struct
  75 {
  76     int bStep;
  77     int bTime;
  78     int bLambda;
  79     int bX;
  80     int bV;
  81     int bF;
  82     int bBox;
  83 } t_count;
  84
  85 typedef struct
  86 {
  87     float bStep;
  88     float bTime;
  89     float bLambda;
  90     float bX;
  91     float bV;
  92     float bF;
  93     float bBox;
  94 } t_fr_time;
  95
  96 static void comp_tpx(const char* fn1, const char* fn2, gmx_bool bRMSD, real ftol, real abstol)
  97 {
  98     const char* ff[2];
  99     t_inputrec* ir[2];
 100     t_state     state[2];
 101     gmx_mtop_t  mtop[2];
 102     int         i;
 103
 104     ff[0] = fn1;
 105     ff[1] = fn2;
 106     for (i = 0; i < (fn2 ? 2 : 1); i++)
 107     {
 108         ir[i] = new t_inputrec();
 109         read_tpx_state(ff[i], ir[i], &state[i], &(mtop[i]));
 110         gmx::MDModules().adjustInputrecBasedOnModules(ir[i]);
 111     }
 112     if (fn2)
 113     {
 114         cmp_inputrec(stdout, ir[0], ir[1], ftol, abstol);
 115         compareMtop(stdout, mtop[0], mtop[1], ftol, abstol);
 116         comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
 117     }
 118     else
 119     {
 120         if (ir[0]->efep == efepNO)
 121         {
 122             fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0]->efep]);
 123         }
 124         else
 125         {
 126             if (ir[0]->bPull)
 127             {
 128                 comp_pull_AB(stdout, ir[0]->pull, ftol, abstol);
 129             }
 130             compareMtopAB(stdout, mtop[0], ftol, abstol);
 131         }
 132     }
 133 }
 134
 135 static void comp_trx(const gmx_output_env_t* oenv, const char* fn1, const char* fn2, gmx_bool bRMSD, real ftol, real abstol)
 136 {
 137     int          i;
 138     const char*  fn[2];
 139     t_trxframe   fr[2];
 140     t_trxstatus* status[2];
 141     gmx_bool     b[2];
 142
 143     fn[0] = fn1;
 144     fn[1] = fn2;
 145     fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
 146     for (i = 0; i < 2; i++)
 147     {
 148         b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X | TRX_READ_V | TRX_READ_F);
 149     }
 150
 151     if (b[0] && b[1])
 152     {
 153         do
 154         {
 155             comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
 156
 157             for (i = 0; i < 2; i++)
 158             {
 159                 b[i] = read_next_frame(oenv, status[i], &fr[i]);
 160             }
 161         } while (b[0] && b[1]);
 162
 163         for (i = 0; i < 2; i++)
 164         {
 165             if (b[i] && !b[1 - i])
 166             {
 167                 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1 - i], fn[i]);
 168             }
 169             close_trx(status[i]);
 170         }
 171     }
 172     if (!b[0] && !b[1])
 173     {
 174         fprintf(stdout, "\nBoth files read correctly\n");
 175     }
 176 }
 177
 178 static void chk_coords(int frame, int natoms, rvec* x, matrix box, real fac, real tol)
 179 {
 180     int  i, j;
 181     int  nNul = 0;
 182     real vol  = det(box);
 183
 184     for (i = 0; (i < natoms); i++)
 185     {
 186         for (j = 0; (j < DIM); j++)
 187         {
 188             if ((vol > 0) && (fabs(x[i][j]) > fac * box[j][j]))
 189             {
 190                 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n", frame, i, x[i][j]);
 191             }
 192         }
 193         if ((fabs(x[i][XX]) < tol) && (fabs(x[i][YY]) < tol) && (fabs(x[i][ZZ]) < tol))
 194         {
 195             nNul++;
 196         }
 197     }
 198     if (nNul > 0)
 199     {
 200         printf("Warning at frame %d: there are %d particles with all coordinates zero\n", frame, nNul);
 201     }
 202 }
 203
 204 static void chk_vels(int frame, int natoms, rvec* v)
 205 {
 206     int i, j;
 207
 208     for (i = 0; (i < natoms); i++)
 209     {
 210         for (j = 0; (j < DIM); j++)
 211         {
 212             if (fabs(v[i][j]) > 500)
 213             {
 214                 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n", frame, i, v[i][j]);
 215             }
 216         }
 217     }
 218 }
 219
 220 static void chk_forces(int frame, int natoms, rvec* f)
 221 {
 222     int i, j;
 223
 224     for (i = 0; (i < natoms); i++)
 225     {
 226         for (j = 0; (j < DIM); j++)
 227         {
 228             if (fabs(f[i][j]) > 10000)
 229             {
 230                 printf("Warning at frame %d. Forces for atom %d are large (%g)\n", frame, i, f[i][j]);
 231             }
 232         }
 233     }
 234 }
 235
 236 static void chk_bonds(const InteractionDefinitions* idef, PbcType pbcType, rvec* x, matrix box, real tol)
 237 {
 238     int   ftype, k, ai, aj, type;
 239     real  b0, blen, deviation;
 240     t_pbc pbc;
 241     rvec  dx;
 242
 243     gmx::ArrayRef<const t_iparams> iparams = idef->iparams;
 244
 245     set_pbc(&pbc, pbcType, box);
 246     for (ftype = 0; (ftype < F_NRE); ftype++)
 247     {
 248         if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
 249         {
 250             for (k = 0; (k < idef->il[ftype].size());)
 251             {
 252                 type = idef->il[ftype].iatoms[k++];
 253                 ai   = idef->il[ftype].iatoms[k++];
 254                 aj   = idef->il[ftype].iatoms[k++];
 255                 b0   = 0;
 256                 switch (ftype)
 257                 {
 258                     case F_BONDS: b0 = iparams[type].harmonic.rA; break;
 259                     case F_G96BONDS: b0 = std::sqrt(iparams[type].harmonic.rA); break;
 260                     case F_MORSE: b0 = iparams[type].morse.b0A; break;
 261                     case F_CUBICBONDS: b0 = iparams[type].cubic.b0; break;
 262                     case F_CONSTR: b0 = iparams[type].constr.dA; break;
 263                     default: break;
 264                 }
 265                 if (b0 != 0)
 266                 {
 267                     pbc_dx(&pbc, x[ai], x[aj], dx);
 268                     blen      = norm(dx);
 269                     deviation = gmx::square(blen - b0);
 270                     if (std::sqrt(deviation / gmx::square(b0)) > tol)
 271                     {
 272                         fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n",
 273                                 ai + 1, aj + 1, blen, b0);
 274                     }
 275                 }
 276             }
 277         }
 278     }
 279 }
 280
 281 static void chk_trj(const gmx_output_env_t* oenv, const char* fn, const char* tpr, real tol)
 282 {
 283     t_trxframe   fr;
 284     t_count      count;
 285     t_fr_time    first, last;
 286     int          j = -1, new_natoms, natoms;
 287     real         old_t1, old_t2;
 288     gmx_bool     bShowTimestep = TRUE, newline = FALSE;
 289     t_trxstatus* status;
 290     gmx_mtop_t   mtop;
 291     t_state      state;
 292     t_inputrec   ir;
 293
 294     std::unique_ptr<gmx_localtop_t> top;
 295     if (tpr)
 296     {
 297         read_tpx_state(tpr, &ir, &state, &mtop);
 298         top = std::make_unique<gmx_localtop_t>(mtop.ffparams);
 299         gmx_mtop_generate_local_top(mtop, top.get(), ir.efep != efepNO);
 300     }
 301     new_natoms = -1;
 302     natoms     = -1;
 303
 304     printf("Checking file %s\n", fn);
 305
 306     j      = 0;
 307     old_t2 = -2.0;
 308     old_t1 = -1.0;
 309
 310     count.bStep   = 0;
 311     count.bTime   = 0;
 312     count.bLambda = 0;
 313     count.bX      = 0;
 314     count.bV      = 0;
 315     count.bF      = 0;
 316     count.bBox    = 0;
 317
 318     first.bStep   = 0;
 319     first.bTime   = 0;
 320     first.bLambda = 0;
 321     first.bX      = 0;
 322     first.bV      = 0;
 323     first.bF      = 0;
 324     first.bBox    = 0;
 325
 326     last.bStep   = 0;
 327     last.bTime   = 0;
 328     last.bLambda = 0;
 329     last.bX      = 0;
 330     last.bV      = 0;
 331     last.bF      = 0;
 332     last.bBox    = 0;
 333
 334     read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
 335
 336     do
 337     {
 338         if (j == 0)
 339         {
 340             fprintf(stderr, "\n# Atoms  %d\n", fr.natoms);
 341             if (fr.bPrec)
 342             {
 343                 fprintf(stderr, "Precision %g (nm)\n", 1 / fr.prec);
 344             }
 345         }
 346         newline = TRUE;
 347         if ((natoms > 0) && (new_natoms != natoms))
 348         {
 349             fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n", old_t1, natoms, new_natoms);
 350             newline = FALSE;
 351         }
 352         if (j >= 2)
 353         {
 354             if (std::fabs((fr.time - old_t1) - (old_t1 - old_t2))
 355                 > 0.1 * (std::fabs(fr.time - old_t1) + std::fabs(old_t1 - old_t2)))
 356             {
 357                 bShowTimestep = FALSE;
 358                 fprintf(stderr, "%sTimesteps at t=%g don't match (%g, %g)\n", newline ? "\n" : "",
 359                         old_t1, old_t1 - old_t2, fr.time - old_t1);
 360             }
 361         }
 362         natoms = new_natoms;
 363         if (tpr)
 364         {
 365             chk_bonds(&top->idef, ir.pbcType, fr.x, fr.box, tol);
 366         }
 367         if (fr.bX)
 368         {
 369             chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
 370         }
 371         if (fr.bV)
 372         {
 373             chk_vels(j, natoms, fr.v);
 374         }
 375         if (fr.bF)
 376         {
 377             chk_forces(j, natoms, fr.f);
 378         }
 379
 380         old_t2 = old_t1;
 381         old_t1 = fr.time;
 382         j++;
 383         new_natoms = fr.natoms;
 384 #define INC(s, n, f, l, item)     \
 385     if ((s).item != 0)            \
 386     {                             \
 387         if ((n).item == 0)        \
 388         {                         \
 389             first.item = fr.time; \
 390         }                         \
 391         last.item = fr.time;      \
 392         (n).item++;               \
 393     }
 394         INC(fr, count, first, last, bStep)
 395         INC(fr, count, first, last, bTime)
 396         INC(fr, count, first, last, bLambda)
 397         INC(fr, count, first, last, bX)
 398         INC(fr, count, first, last, bV)
 399         INC(fr, count, first, last, bF)
 400         INC(fr, count, first, last, bBox)
 401 #undef INC
 402     } while (read_next_frame(oenv, status, &fr));
 403
 404     fprintf(stderr, "\n");
 405
 406     close_trx(status);
 407
 408     fprintf(stderr, "\nItem        #frames");
 409     if (bShowTimestep)
 410     {
 411         fprintf(stderr, " Timestep (ps)");
 412     }
 413     fprintf(stderr, "\n");
 414 #define PRINTITEM(label, item)                                                    \
 415     fprintf(stderr, "%-10s  %6d", label, count.item);                             \
 416     if ((bShowTimestep) && (count.item > 1))                                      \
 417     {                                                                             \
 418         fprintf(stderr, "    %g\n", (last.item - first.item) / (count.item - 1)); \
 419     }                                                                             \
 420     else                                                                          \
 421         fprintf(stderr, "\n")
 422     PRINTITEM("Step", bStep);
 423     PRINTITEM("Time", bTime);
 424     PRINTITEM("Lambda", bLambda);
 425     PRINTITEM("Coords", bX);
 426     PRINTITEM("Velocities", bV);
 427     PRINTITEM("Forces", bF);
 428     PRINTITEM("Box", bBox);
 429 }
 430
 431 static void chk_tps(const char* fn, real vdw_fac, real bon_lo, real bon_hi)
 432 {
 433     int        natom, i, j, k;
 434     t_topology top;
 435     PbcType    pbcType;
 436     t_atoms*   atoms;
 437     rvec *     x, *v;
 438     rvec       dx;
 439     matrix     box;
 440     t_pbc      pbc;
 441     gmx_bool   bV, bX, bB, bFirst, bOut;
 442     real       r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
 443     real*      atom_vdw;
 444
 445     fprintf(stderr, "Checking coordinate file %s\n", fn);
 446     read_tps_conf(fn, &top, &pbcType, &x, &v, box, TRUE);
 447     atoms = &top.atoms;
 448     natom = atoms->nr;
 449     fprintf(stderr, "%d atoms in file\n", atoms->nr);
 450
 451     /* check coordinates and box */
 452     bV = FALSE;
 453     bX = FALSE;
 454     for (i = 0; (i < natom) && !(bV && bX); i++)
 455     {
 456         for (j = 0; (j < DIM) && !(bV && bX); j++)
 457         {
 458             bV = bV || (v[i][j] != 0);
 459             bX = bX || (x[i][j] != 0);
 460         }
 461     }
 462     bB = FALSE;
 463     for (i = 0; (i < DIM) && !bB; i++)
 464     {
 465         for (j = 0; (j < DIM) && !bB; j++)
 466         {
 467             bB = bB || (box[i][j] != 0);
 468         }
 469     }
 470
 471     fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
 472     fprintf(stderr, "box         %s\n", bB ? "found" : "absent");
 473     fprintf(stderr, "velocities  %s\n", bV ? "found" : "absent");
 474     fprintf(stderr, "\n");
 475
 476     /* check velocities */
 477     if (bV)
 478     {
 479         ekin = 0.0;
 480         for (i = 0; (i < natom); i++)
 481         {
 482             for (j = 0; (j < DIM); j++)
 483             {
 484                 ekin += 0.5 * atoms->atom[i].m * v[i][j] * v[i][j];
 485             }
 486         }
 487         temp1 = (2.0 * ekin) / (natom * DIM * BOLTZ);
 488         temp2 = (2.0 * ekin) / (natom * (DIM - 1) * BOLTZ);
 489         fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
 490         fprintf(stderr,
 491                 "Assuming the number of degrees of freedom to be "
 492                 "Natoms * %d or Natoms * %d,\n"
 493                 "the velocities correspond to a temperature of the system\n"
 494                 "of %g K or %g K respectively.\n\n",
 495                 DIM, DIM - 1, temp1, temp2);
 496     }
 497
 498     /* check coordinates */
 499     if (bX)
 500     {
 501         vdwfac2 = gmx::square(vdw_fac);
 502         bonlo2  = gmx::square(bon_lo);
 503         bonhi2  = gmx::square(bon_hi);
 504
 505         fprintf(stderr,
 506                 "Checking for atoms closer than %g and not between %g and %g,\n"
 507                 "relative to sum of Van der Waals distance:\n",
 508                 vdw_fac, bon_lo, bon_hi);
 509         snew(atom_vdw, natom);
 510         AtomProperties aps;
 511         for (i = 0; (i < natom); i++)
 512         {
 513             aps.setAtomProperty(epropVDW, *(atoms->resinfo[atoms->atom[i].resind].name),
 514                                 *(atoms->atomname[i]), &(atom_vdw[i]));
 515             if (debug)
 516             {
 517                 fprintf(debug, "%5d %4s %4s %7g\n", i + 1, *(atoms->resinfo[atoms->atom[i].resind].name),
 518                         *(atoms->atomname[i]), atom_vdw[i]);
 519             }
 520         }
 521         if (bB)
 522         {
 523             set_pbc(&pbc, pbcType, box);
 524         }
 525
 526         bFirst = TRUE;
 527         for (i = 0; (i < natom); i++)
 528         {
 529             if (((i + 1) % 10) == 0)
 530             {
 531                 fprintf(stderr, "\r%5d", i + 1);
 532                 fflush(stderr);
 533             }
 534             for (j = i + 1; (j < natom); j++)
 535             {
 536                 if (bB)
 537                 {
 538                     pbc_dx(&pbc, x[i], x[j], dx);
 539                 }
 540                 else
 541                 {
 542                     rvec_sub(x[i], x[j], dx);
 543                 }
 544                 r2    = iprod(dx, dx);
 545                 dist2 = gmx::square(atom_vdw[i] + atom_vdw[j]);
 546                 if ((r2 <= dist2 * bonlo2) || ((r2 >= dist2 * bonhi2) && (r2 <= dist2 * vdwfac2)))
 547                 {
 548                     if (bFirst)
 549                     {
 550                         fprintf(stderr, "\r%5s %4s %8s %5s  %5s %4s %8s %5s  %6s\n", "atom#", "name",
 551                                 "residue", "r_vdw", "atom#", "name", "residue", "r_vdw", "distance");
 552                         bFirst = FALSE;
 553                     }
 554                     fprintf(stderr, "\r%5d %4s %4s%4d %-5.3g  %5d %4s %4s%4d %-5.3g  %-6.4g\n", i + 1,
 555                             *(atoms->atomname[i]), *(atoms->resinfo[atoms->atom[i].resind].name),
 556                             atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i], j + 1,
 557                             *(atoms->atomname[j]), *(atoms->resinfo[atoms->atom[j].resind].name),
 558                             atoms->resinfo[atoms->atom[j].resind].nr, atom_vdw[j], std::sqrt(r2));
 559                 }
 560             }
 561         }
 562         if (bFirst)
 563         {
 564             fprintf(stderr, "\rno close atoms found\n");
 565         }
 566         fprintf(stderr, "\r      \n");
 567
 568         if (bB)
 569         {
 570             /* check box */
 571             bFirst = TRUE;
 572             k      = 0;
 573             for (i = 0; (i < natom) && (k < 10); i++)
 574             {
 575                 bOut = FALSE;
 576                 for (j = 0; (j < DIM) && !bOut; j++)
 577                 {
 578                     bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
 579                 }
 580                 if (bOut)
 581                 {
 582                     k++;
 583                     if (bFirst)
 584                     {
 585                         fprintf(stderr, "Atoms outside box ( ");
 586                         for (j = 0; (j < DIM); j++)
 587                         {
 588                             fprintf(stderr, "%g ", box[j][j]);
 589                         }
 590                         fprintf(stderr,
 591                                 "):\n"
 592                                 "(These may occur often and are normally not a problem)\n"
 593                                 "%5s %4s %8s %5s  %s\n",
 594                                 "atom#", "name", "residue", "r_vdw", "coordinate");
 595                         bFirst = FALSE;
 596                     }
 597                     fprintf(stderr, "%5d %4s %4s%4d %-5.3g", i, *(atoms->atomname[i]),
 598                             *(atoms->resinfo[atoms->atom[i].resind].name),
 599                             atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
 600                     for (j = 0; (j < DIM); j++)
 601                     {
 602                         fprintf(stderr, " %6.3g", x[i][j]);
 603                     }
 604                     fprintf(stderr, "\n");
 605                 }
 606             }
 607             if (k == 10)
 608             {
 609                 fprintf(stderr, "(maybe more)\n");
 610             }
 611             if (bFirst)
 612             {
 613                 fprintf(stderr, "no atoms found outside box\n");
 614             }
 615             fprintf(stderr, "\n");
 616         }
 617     }
 618 }
 619
 620 static void chk_ndx(const char* fn)
 621 {
 622     t_blocka* grps;
 623     char**    grpname;
 624     int       i;
 625
 626     grps = init_index(fn, &grpname);
 627     if (debug)
 628     {
 629         pr_blocka(debug, 0, fn, grps, FALSE);
 630     }
 631     else
 632     {
 633         printf("Contents of index file %s\n", fn);
 634         printf("--------------------------------------------------\n");
 635         printf("Nr.   Group               #Entries   First    Last\n");
 636         for (i = 0; (i < grps->nr); i++)
 637         {
 638             printf("%4d  %-20s%8d%8d%8d\n", i, grpname[i], grps->index[i + 1] - grps->index[i],
 639                    grps->a[grps->index[i]] + 1, grps->a[grps->index[i + 1] - 1] + 1);
 640         }
 641     }
 642     for (i = 0; (i < grps->nr); i++)
 643     {
 644         sfree(grpname[i]);
 645     }
 646     sfree(grpname);
 647     done_blocka(grps);
 648 }
 649
 650 static void chk_enx(const char* fn)
 651 {
 652     int          nre, fnr;
 653     ener_file_t  in;
 654     gmx_enxnm_t* enm = nullptr;
 655     t_enxframe*  fr;
 656     gmx_bool     bShowTStep;
 657     gmx_bool     timeSet;
 658     real         t0, old_t1, old_t2;
 659     char         buf[22];
 660
 661     fprintf(stderr, "Checking energy file %s\n\n", fn);
 662
 663     in = open_enx(fn, "r");
 664     do_enxnms(in, &nre, &enm);
 665     fprintf(stderr, "%d groups in energy file", nre);
 666     snew(fr, 1);
 667     old_t2     = -2.0;
 668     old_t1     = -1.0;
 669     fnr        = 0;
 670     t0         = 0;
 671     timeSet    = FALSE;
 672     bShowTStep = TRUE;
 673
 674     while (do_enx(in, fr))
 675     {
 676         if (fnr >= 2)
 677         {
 678             if (fabs((fr->t - old_t1) - (old_t1 - old_t2))
 679                 > 0.1 * (fabs(fr->t - old_t1) + std::fabs(old_t1 - old_t2)))
 680             {
 681                 bShowTStep = FALSE;
 682                 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n", old_t1,
 683                         old_t1 - old_t2, fr->t - old_t1);
 684             }
 685         }
 686         old_t2 = old_t1;
 687         old_t1 = fr->t;
 688         if (!timeSet)
 689         {
 690             t0      = fr->t;
 691             timeSet = TRUE;
 692         }
 693         if (fnr == 0)
 694         {
 695             fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n", gmx_step_str(fr->step, buf),
 696                     fnr, fr->t);
 697         }
 698         fnr++;
 699     }
 700     fprintf(stderr, "\n\nFound %d frames", fnr);
 701     if (bShowTStep && fnr > 1)
 702     {
 703         fprintf(stderr, " with a timestep of %g ps", (old_t1 - t0) / (fnr - 1));
 704     }
 705     fprintf(stderr, ".\n");
 706
 707     free_enxframe(fr);
 708     free_enxnms(nre, enm);
 709     sfree(fr);
 710 }
 711
 712 int gmx_check(int argc, char* argv[])
 713 {
 714     const char* desc[] = {
 715         "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
 716         "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
 717         "or an index file ([REF].ndx[ref])",
 718         "and prints out useful information about them.[PAR]",
 719         "Option [TT]-c[tt] checks for presence of coordinates,",
 720         "velocities and box in the file, for close contacts (smaller than",
 721         "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
 722         "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
 723         "radii) and atoms outside the box (these may occur often and are",
 724         "no problem). If velocities are present, an estimated temperature",
 725         "will be calculated from them.[PAR]",
 726         "If an index file, is given its contents will be summarized.[PAR]",
 727         "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
 728         "the program will check whether the bond lengths defined in the tpr",
 729         "file are indeed correct in the trajectory. If not you may have",
 730         "non-matching files due to e.g. deshuffling or due to problems with",
 731         "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for ",
 732         "such problems.[PAR]",
 733         "The program can compare two run input ([REF].tpr[ref])",
 734         "files",
 735         "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied. When comparing",
 736         "run input files this way, the default relative tolerance is reduced",
 737         "to 0.000001 and the absolute tolerance set to zero to find any differences",
 738         "not due to minor compiler optimization differences, although you can",
 739         "of course still set any other tolerances through the options.",
 740         "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
 741         "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
 742         "For free energy simulations the A and B state topology from one",
 743         "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]"
 744     };
 745     t_filenm fnm[] = { { efTRX, "-f", nullptr, ffOPTRD },  { efTRX, "-f2", nullptr, ffOPTRD },
 746                        { efTPR, "-s1", "top1", ffOPTRD },  { efTPR, "-s2", "top2", ffOPTRD },
 747                        { efTPS, "-c", nullptr, ffOPTRD },  { efEDR, "-e", nullptr, ffOPTRD },
 748                        { efEDR, "-e2", "ener2", ffOPTRD }, { efNDX, "-n", nullptr, ffOPTRD },
 749                        { efTEX, "-m", nullptr, ffOPTWR } };
 750 #define NFILE asize(fnm)
 751     const char *fn1 = nullptr, *fn2 = nullptr, *tex = nullptr;
 752
 753     gmx_output_env_t* oenv;
 754     static real       vdw_fac  = 0.8;
 755     static real       bon_lo   = 0.4;
 756     static real       bon_hi   = 0.7;
 757     static gmx_bool   bRMSD    = FALSE;
 758     static real       ftol     = 0.001;
 759     static real       abstol   = 0.001;
 760     static gmx_bool   bCompAB  = FALSE;
 761     static char*      lastener = nullptr;
 762     static t_pargs    pa[]     = {
 763         { "-vdwfac",
 764           FALSE,
 765           etREAL,
 766           { &vdw_fac },
 767           "Fraction of sum of VdW radii used as warning cutoff" },
 768         { "-bonlo", FALSE, etREAL, { &bon_lo }, "Min. fract. of sum of VdW radii for bonded atoms" },
 769         { "-bonhi", FALSE, etREAL, { &bon_hi }, "Max. fract. of sum of VdW radii for bonded atoms" },
 770         { "-rmsd", FALSE, etBOOL, { &bRMSD }, "Print RMSD for x, v and f" },
 771         { "-tol",
 772           FALSE,
 773           etREAL,
 774           { &ftol },
 775           "Relative tolerance for comparing real values defined as "
 776           "[MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
 777         { "-abstol",
 778           FALSE,
 779           etREAL,
 780           { &abstol },
 781           "Absolute tolerance, useful when sums are close to zero." },
 782         { "-ab", FALSE, etBOOL, { &bCompAB }, "Compare the A and B topology from one file" },
 783         { "-lastener",
 784           FALSE,
 785           etSTR,
 786           { &lastener },
 787           "Last energy term to compare (if not given all are tested). It makes sense to go up "
 788           "until the Pressure." }
 789     };
 790
 791     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 792     {
 793         return 0;
 794     }
 795
 796     fn1 = opt2fn_null("-f", NFILE, fnm);
 797     fn2 = opt2fn_null("-f2", NFILE, fnm);
 798     tex = opt2fn_null("-m", NFILE, fnm);
 799
 800     if (tex)
 801     {
 802         fprintf(stderr,
 803                 "LaTeX file writing has been removed from gmx check. "
 804                 "Please use gmx report-methods instead for it.\n");
 805     }
 806     if (fn1 && fn2)
 807     {
 808         comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
 809     }
 810     else if (fn1)
 811     {
 812         chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
 813     }
 814     else if (fn2)
 815     {
 816         fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
 817     }
 818
 819     fn1 = opt2fn_null("-s1", NFILE, fnm);
 820     fn2 = opt2fn_null("-s2", NFILE, fnm);
 821     if ((fn1 && fn2) || bCompAB)
 822     {
 823         if (bCompAB)
 824         {
 825             if (fn1 == nullptr)
 826             {
 827                 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
 828             }
 829             fn2 = nullptr;
 830         }
 831
 832         fprintf(stderr, "Note: When comparing run input files, default tolerances are reduced.\n");
 833         if (!opt2parg_bSet("-tol", asize(pa), pa))
 834         {
 835             ftol = 0.000001;
 836         }
 837         if (!opt2parg_bSet("-abstol", asize(pa), pa))
 838         {
 839             abstol = 0;
 840         }
 841         comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
 842     }
 843     else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
 844     {
 845         fprintf(stderr, "Please give me TWO run input (.tpr) files\n");
 846     }
 847
 848     fn1 = opt2fn_null("-e", NFILE, fnm);
 849     fn2 = opt2fn_null("-e2", NFILE, fnm);
 850     if (fn1 && fn2)
 851     {
 852         comp_enx(fn1, fn2, ftol, abstol, lastener);
 853     }
 854     else if (fn1)
 855     {
 856         chk_enx(ftp2fn(efEDR, NFILE, fnm));
 857     }
 858     else if (fn2)
 859     {
 860         fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
 861     }
 862
 863     if (ftp2bSet(efTPS, NFILE, fnm))
 864     {
 865         chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
 866     }
 867
 868     if (ftp2bSet(efNDX, NFILE, fnm))
 869     {
 870         chk_ndx(ftp2fn(efNDX, NFILE, fnm));
 871     }
 872
 873     return 0;
 874 }