src/programs/gmx/gmxcheck.c

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