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