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