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