src/gromacs/tools/check.cpp

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