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