src/programs/gmx/gmxcheck.c

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