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