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