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