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