src/gromacs/tools/check.c

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