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