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