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