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         case F_G96BONDS:
 220           b0 = idef->iparams[type].harmonic.rA;
 221           break;
 222         case F_MORSE:
 223           b0 = idef->iparams[type].morse.b0A;
 224           break;
 225         case F_CUBICBONDS:
 226           b0 = idef->iparams[type].cubic.b0;
 227           break;
 228         case F_CONSTR:
 229           b0 = idef->iparams[type].constr.dA;
 230           break;
 231         default:
 232           break;
 233         }
 234         if (b0 != 0) {
 235           pbc_dx(&pbc,x[ai],x[aj],dx);
 236           blen = norm(dx);
 237           deviation = sqr(blen-b0);
 238           if (sqrt(deviation/sqr(b0) > tol)) {
 239             fprintf(stderr,"Distance between atoms %d and %d is %.3f, should be %.3f\n",ai+1,aj+1,blen,b0);
 240           }
 241         }
 242       }
 243     }
 244 }
 245
 246 void chk_trj(const output_env_t oenv,const char *fn,const char *tpr,real tol)
 247 {
 248   t_trxframe   fr;
 249   t_count      count;
 250   t_fr_time    first,last;
 251   int          j=-1,new_natoms,natoms;
 252   gmx_off_t    fpos;
 253   real         rdum,tt,old_t1,old_t2,prec;
 254   gmx_bool         bShowTimestep=TRUE,bOK,newline=FALSE;
 255   t_trxstatus *status;
 256   gmx_mtop_t   mtop;
 257   gmx_localtop_t *top=NULL;
 258   t_state      state;
 259   t_inputrec   ir;
 260
 261   if (tpr) {
 262     read_tpx_state(tpr,&ir,&state,NULL,&mtop);
 263     top = gmx_mtop_generate_local_top(&mtop,&ir);
 264   }
 265   new_natoms = -1;
 266   natoms = -1;
 267
 268   printf("Checking file %s\n",fn);
 269
 270   j      =  0;
 271   old_t2 = -2.0;
 272   old_t1 = -1.0;
 273   fpos   = 0;
 274
 275   count.bStep = 0;
 276   count.bTime = 0;
 277   count.bLambda = 0;
 278   count.bX = 0;
 279   count.bV = 0;
 280   count.bF = 0;
 281   count.bBox = 0;
 282
 283   first.bStep = 0;
 284   first.bTime = 0;
 285   first.bLambda = 0;
 286   first.bX = 0;
 287   first.bV = 0;
 288   first.bF = 0;
 289   first.bBox = 0;
 290
 291   last.bStep = 0;
 292   last.bTime = 0;
 293   last.bLambda = 0;
 294   last.bX = 0;
 295   last.bV = 0;
 296   last.bF = 0;
 297   last.bBox = 0;
 298
 299   read_first_frame(oenv,&status,fn,&fr,TRX_READ_X | TRX_READ_V | TRX_READ_F);
 300
 301   do {
 302     if (j == 0) {
 303       fprintf(stderr,"\n# Atoms  %d\n",fr.natoms);
 304       if (fr.bPrec)
 305         fprintf(stderr,"Precision %g (nm)\n",1/fr.prec);
 306     }
 307     newline=TRUE;
 308     if ((natoms > 0) && (new_natoms != natoms)) {
 309       fprintf(stderr,"\nNumber of atoms at t=%g don't match (%d, %d)\n",
 310               old_t1,natoms,new_natoms);
 311       newline=FALSE;
 312     }
 313     if (j>=2) {
 314       if ( fabs((fr.time-old_t1)-(old_t1-old_t2)) >
 315            0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) ) {
 316         bShowTimestep=FALSE;
 317         fprintf(stderr,"%sTimesteps at t=%g don't match (%g, %g)\n",
 318                 newline?"\n":"",old_t1,old_t1-old_t2,fr.time-old_t1);
 319       }
 320     }
 321     natoms=new_natoms;
 322     if (tpr) {
 323       chk_bonds(&top->idef,ir.ePBC,fr.x,fr.box,tol);
 324     }
 325     if (fr.bX)
 326       chk_coords(j,natoms,fr.x,fr.box,1e5,tol);
 327     if (fr.bV)
 328       chk_vels(j,natoms,fr.v);
 329     if (fr.bF)
 330       chk_forces(j,natoms,fr.f);
 331
 332     old_t2=old_t1;
 333     old_t1=fr.time;
 334     /*if (fpos && ((j<10 || j%10==0)))
 335       fprintf(stderr," byte: %10lu",(unsigned long)fpos);*/
 336     j++;
 337     new_natoms=fr.natoms;
 338 #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++; }
 339     INC(fr,count,first,last,bStep);
 340     INC(fr,count,first,last,bTime);
 341     INC(fr,count,first,last,bLambda);
 342     INC(fr,count,first,last,bX);
 343     INC(fr,count,first,last,bV);
 344     INC(fr,count,first,last,bF);
 345     INC(fr,count,first,last,bBox);
 346 #undef INC
 347     fpos = gmx_fio_ftell(trx_get_fileio(status));
 348   } while (read_next_frame(oenv,status,&fr));
 349
 350   fprintf(stderr,"\n");
 351
 352   close_trj(status);
 353
 354   fprintf(stderr,"\nItem        #frames");
 355   if (bShowTimestep)
 356     fprintf(stderr," Timestep (ps)");
 357   fprintf(stderr,"\n");
 358 #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")
 359   PRINTITEM ( "Step",       bStep );
 360   PRINTITEM ( "Time",       bTime );
 361   PRINTITEM ( "Lambda",     bLambda );
 362   PRINTITEM ( "Coords",     bX );
 363   PRINTITEM ( "Velocities", bV );
 364   PRINTITEM ( "Forces",     bF );
 365   PRINTITEM ( "Box",        bBox );
 366 }
 367
 368 void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
 369 {
 370   int       natom,i,j,k;
 371   char      title[STRLEN];
 372   t_topology top;
 373   int       ePBC;
 374   t_atoms   *atoms;
 375   rvec      *x,*v;
 376   rvec      dx;
 377   matrix    box;
 378   t_pbc     pbc;
 379   gmx_bool      bV,bX,bB,bFirst,bOut;
 380   real      r2,ekin,temp1,temp2,dist2,vdwfac2,bonlo2,bonhi2;
 381   real      *atom_vdw;
 382   gmx_atomprop_t aps;
 383
 384   fprintf(stderr,"Checking coordinate file %s\n",fn);
 385   read_tps_conf(fn,title,&top,&ePBC,&x,&v,box,TRUE);
 386   atoms=&top.atoms;
 387   natom=atoms->nr;
 388   fprintf(stderr,"%d atoms in file\n",atoms->nr);
 389
 390   /* check coordinates and box */
 391   bV=FALSE;
 392   bX=FALSE;
 393   for (i=0; (i<natom) && !(bV && bX); i++)
 394     for (j=0; (j<DIM) && !(bV && bX); j++) {
 395       bV=bV || (v[i][j]!=0);
 396       bX=bX || (x[i][j]!=0);
 397     }
 398   bB=FALSE;
 399   for (i=0; (i<DIM) && !bB; i++)
 400     for (j=0; (j<DIM) && !bB; j++)
 401       bB=bB || (box[i][j]!=0);
 402
 403   fprintf(stderr,"coordinates %s\n",bX?"found":"absent");
 404   fprintf(stderr,"box         %s\n",bB?"found":"absent");
 405   fprintf(stderr,"velocities  %s\n",bV?"found":"absent");
 406   fprintf(stderr,"\n");
 407
 408   /* check velocities */
 409   if (bV) {
 410     ekin=0.0;
 411     for (i=0; (i<natom); i++)
 412       for (j=0; (j<DIM); j++)
 413         ekin+=0.5*atoms->atom[i].m*v[i][j]*v[i][j];
 414     temp1=(2.0*ekin)/(natom*DIM*BOLTZ);
 415     temp2=(2.0*ekin)/(natom*(DIM-1)*BOLTZ);
 416     fprintf(stderr,"Kinetic energy: %g (kJ/mol)\n",ekin);
 417     fprintf(stderr,"Assuming the number of degrees of freedom to be "
 418             "Natoms * %d or Natoms * %d,\n"
 419             "the velocities correspond to a temperature of the system\n"
 420             "of %g K or %g K respectively.\n\n",DIM,DIM-1,temp1,temp2);
 421   }
 422
 423   /* check coordinates */
 424   if (bX) {
 425     vdwfac2=sqr(vdw_fac);
 426     bonlo2=sqr(bon_lo);
 427     bonhi2=sqr(bon_hi);
 428
 429     fprintf(stderr,
 430             "Checking for atoms closer than %g and not between %g and %g,\n"
 431             "relative to sum of Van der Waals distance:\n",
 432             vdw_fac,bon_lo,bon_hi);
 433     snew(atom_vdw,natom);
 434     aps = gmx_atomprop_init();
 435     for (i=0; (i<natom); i++) {
 436       gmx_atomprop_query(aps,epropVDW,
 437                          *(atoms->resinfo[atoms->atom[i].resind].name),
 438                          *(atoms->atomname[i]),&(atom_vdw[i]));
 439       if (debug) fprintf(debug,"%5d %4s %4s %7g\n",i+1,
 440                          *(atoms->resinfo[atoms->atom[i].resind].name),
 441                          *(atoms->atomname[i]),
 442                          atom_vdw[i]);
 443     }
 444     gmx_atomprop_destroy(aps);
 445     if (bB)
 446       set_pbc(&pbc,ePBC,box);
 447
 448     bFirst=TRUE;
 449     for (i=0; (i<natom); i++) {
 450       if (((i+1)%10)==0)
 451         fprintf(stderr,"\r%5d",i+1);
 452       for (j=i+1; (j<natom); j++) {
 453         if (bB)
 454           pbc_dx(&pbc,x[i],x[j],dx);
 455         else
 456           rvec_sub(x[i],x[j],dx);
 457         r2=iprod(dx,dx);
 458         dist2=sqr(atom_vdw[i]+atom_vdw[j]);
 459         if ( (r2<=dist2*bonlo2) ||
 460              ( (r2>=dist2*bonhi2) && (r2<=dist2*vdwfac2) ) ) {
 461           if (bFirst) {
 462             fprintf(stderr,"\r%5s %4s %8s %5s  %5s %4s %8s %5s  %6s\n",
 463                     "atom#","name","residue","r_vdw",
 464                     "atom#","name","residue","r_vdw","distance");
 465             bFirst=FALSE;
 466           }
 467           fprintf(stderr,
 468                   "\r%5d %4s %4s%4d %-5.3g  %5d %4s %4s%4d %-5.3g  %-6.4g\n",
 469                   i+1,*(atoms->atomname[i]),
 470                   *(atoms->resinfo[atoms->atom[i].resind].name),
 471                   atoms->resinfo[atoms->atom[i].resind].nr,
 472                   atom_vdw[i],
 473                   j+1,*(atoms->atomname[j]),
 474                   *(atoms->resinfo[atoms->atom[j].resind].name),
 475                   atoms->resinfo[atoms->atom[j].resind].nr,
 476                   atom_vdw[j],
 477                   sqrt(r2) );
 478         }
 479       }
 480     }
 481     if (bFirst)
 482       fprintf(stderr,"\rno close atoms found\n");
 483     fprintf(stderr,"\r      \n");
 484
 485     if (bB) {
 486       /* check box */
 487       bFirst=TRUE;
 488       k=0;
 489       for (i=0; (i<natom) && (k<10); i++) {
 490         bOut=FALSE;
 491         for(j=0; (j<DIM) && !bOut; j++)
 492           bOut = bOut || (x[i][j]<0) || (x[i][j]>box[j][j]);
 493         if (bOut) {
 494           k++;
 495           if (bFirst) {
 496             fprintf(stderr,"Atoms outside box ( ");
 497             for (j=0; (j<DIM); j++)
 498               fprintf(stderr,"%g ",box[j][j]);
 499             fprintf(stderr,"):\n"
 500                     "(These may occur often and are normally not a problem)\n"
 501                     "%5s %4s %8s %5s  %s\n",
 502                     "atom#","name","residue","r_vdw","coordinate");
 503             bFirst=FALSE;
 504           }
 505           fprintf(stderr,
 506                   "%5d %4s %4s%4d %-5.3g",
 507                   i,*(atoms->atomname[i]),
 508                   *(atoms->resinfo[atoms->atom[i].resind].name),
 509                   atoms->resinfo[atoms->atom[i].resind].nr,atom_vdw[i]);
 510           for (j=0; (j<DIM); j++)
 511             fprintf(stderr," %6.3g",x[i][j]);
 512           fprintf(stderr,"\n");
 513         }
 514       }
 515       if (k==10)
 516         fprintf(stderr,"(maybe more)\n");
 517       if (bFirst)
 518         fprintf(stderr,"no atoms found outside box\n");
 519       fprintf(stderr,"\n");
 520     }
 521   }
 522 }
 523
 524 void chk_ndx(const char *fn)
 525 {
 526   t_blocka *grps;
 527   char **grpname;
 528   int  i,j;
 529
 530   grps = init_index(fn,&grpname);
 531   if (debug)
 532     pr_blocka(debug,0,fn,grps,FALSE);
 533   else {
 534     printf("Contents of index file %s\n",fn);
 535     printf("--------------------------------------------------\n");
 536     printf("Nr.   Group               #Entries   First    Last\n");
 537     for(i=0; (i<grps->nr); i++) {
 538       printf("%4d  %-20s%8d%8d%8d\n",i,grpname[i],
 539              grps->index[i+1]-grps->index[i],
 540              grps->a[grps->index[i]]+1,
 541              grps->a[grps->index[i+1]-1]+1);
 542     }
 543   }
 544   for(i=0; (i<grps->nr); i++)
 545     sfree(grpname[i]);
 546   sfree(grpname);
 547   done_blocka(grps);
 548 }
 549
 550 void chk_enx(const char *fn)
 551 {
 552   int        nre,fnr,ndr;
 553   ener_file_t in;
 554   gmx_enxnm_t *enm=NULL;
 555   t_enxframe *fr;
 556   gmx_bool       bShowTStep;
 557   real       t0,old_t1,old_t2;
 558   char       buf[22];
 559
 560   fprintf(stderr,"Checking energy file %s\n\n",fn);
 561
 562   in = open_enx(fn,"r");
 563   do_enxnms(in,&nre,&enm);
 564   fprintf(stderr,"%d groups in energy file",nre);
 565   snew(fr,1);
 566   old_t2=-2.0;
 567   old_t1=-1.0;
 568   fnr=0;
 569   t0=NOTSET;
 570   bShowTStep=TRUE;
 571
 572   while (do_enx(in,fr)) {
 573     if (fnr>=2) {
 574       if ( fabs((fr->t-old_t1)-(old_t1-old_t2)) >
 575            0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) ) {
 576         bShowTStep=FALSE;
 577         fprintf(stderr,"\nTimesteps at t=%g don't match (%g, %g)\n",
 578                 old_t1,old_t1-old_t2,fr->t-old_t1);
 579       }
 580     }
 581     old_t2=old_t1;
 582     old_t1=fr->t;
 583     if (t0 == NOTSET) t0=fr->t;
 584     if (fnr == 0)
 585       fprintf(stderr,"\rframe: %6s (index %6d), t: %10.3f\n",
 586               gmx_step_str(fr->step,buf),fnr,fr->t);
 587     fnr++;
 588   }
 589   fprintf(stderr,"\n\nFound %d frames",fnr);
 590   if (bShowTStep && fnr>1)
 591     fprintf(stderr," with a timestep of %g ps",(old_t1-t0)/(fnr-1));
 592   fprintf(stderr,".\n");
 593
 594   free_enxframe(fr);
 595   free_enxnms(nre,enm);
 596   sfree(fr);
 597 }
 598
 599 int cmain(int argc,char *argv[])
 600 {
 601   const char *desc[] = {
 602     "[TT]gmxcheck[tt] reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
 603     "[TT].xtc[tt]), an energy file ([TT].ene[tt] or [TT].edr[tt])",
 604     "or an index file ([TT].ndx[tt])",
 605     "and prints out useful information about them.[PAR]",
 606     "Option [TT]-c[tt] checks for presence of coordinates,",
 607     "velocities and box in the file, for close contacts (smaller than",
 608     "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
 609     "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
 610     "radii) and atoms outside the box (these may occur often and are",
 611     "no problem). If velocities are present, an estimated temperature",
 612     "will be calculated from them.[PAR]",
 613     "If an index file, is given its contents will be summarized.[PAR]",
 614     "If both a trajectory and a [TT].tpr[tt] file are given (with [TT]-s1[tt])",
 615     "the program will check whether the bond lengths defined in the tpr",
 616     "file are indeed correct in the trajectory. If not you may have",
 617     "non-matching files due to e.g. deshuffling or due to problems with",
 618     "virtual sites. With these flags, [TT]gmxcheck[tt] provides a quick check for such problems.[PAR]",
 619     "The program can compare two run input ([TT].tpr[tt], [TT].tpb[tt] or",
 620     "[TT].tpa[tt]) files",
 621     "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
 622     "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
 623     "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
 624     "For free energy simulations the A and B state topology from one",
 625     "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
 626     "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
 627     "consisting of a rough outline for a methods section for a paper."
 628   };
 629   t_filenm fnm[] = {
 630     { efTRX, "-f",  NULL, ffOPTRD },
 631     { efTRX, "-f2",  NULL, ffOPTRD },
 632     { efTPX, "-s1", "top1", ffOPTRD },
 633     { efTPX, "-s2", "top2", ffOPTRD },
 634     { efTPS, "-c",  NULL, ffOPTRD },
 635     { efEDR, "-e",  NULL, ffOPTRD },
 636     { efEDR, "-e2", "ener2", ffOPTRD },
 637     { efNDX, "-n",  NULL, ffOPTRD },
 638     { efTEX, "-m",  NULL, ffOPTWR }
 639   };
 640 #define NFILE asize(fnm)
 641   const char *fn1=NULL,*fn2=NULL,*tex=NULL;
 642
 643   output_env_t oenv;
 644   static real vdw_fac=0.8;
 645   static real bon_lo=0.4;
 646   static real bon_hi=0.7;
 647   static gmx_bool bRMSD=FALSE;
 648   static real ftol=0.001;
 649   static real abstol=0.001;
 650   static gmx_bool bCompAB=FALSE;
 651   static char *lastener=NULL;
 652   static t_pargs pa[] = {
 653     { "-vdwfac", FALSE, etREAL, {&vdw_fac},
 654       "Fraction of sum of VdW radii used as warning cutoff" },
 655     { "-bonlo",  FALSE, etREAL, {&bon_lo},
 656       "Min. fract. of sum of VdW radii for bonded atoms" },
 657     { "-bonhi",  FALSE, etREAL, {&bon_hi},
 658       "Max. fract. of sum of VdW radii for bonded atoms" },
 659     { "-rmsd",   FALSE, etBOOL, {&bRMSD},
 660       "Print RMSD for x, v and f" },
 661     { "-tol",    FALSE, etREAL, {&ftol},
 662       "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
 663     { "-abstol",    FALSE, etREAL, {&abstol},
 664       "Absolute tolerance, useful when sums are close to zero." },
 665     { "-ab",     FALSE, etBOOL, {&bCompAB},
 666       "Compare the A and B topology from one file" },
 667     { "-lastener",FALSE, etSTR,  {&lastener},
 668       "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
 669   };
 670
 671   CopyRight(stderr,argv[0]);
 672   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 673                     asize(desc),desc,0,NULL,&oenv);
 674
 675   fn1 = opt2fn_null("-f",NFILE,fnm);
 676   fn2 = opt2fn_null("-f2",NFILE,fnm);
 677   tex = opt2fn_null("-m",NFILE,fnm);
 678   if (fn1 && fn2)
 679     comp_trx(oenv,fn1,fn2,bRMSD,ftol,abstol);
 680   else if (fn1)
 681     chk_trj(oenv,fn1,opt2fn_null("-s1",NFILE,fnm),ftol);
 682   else if (fn2)
 683     fprintf(stderr,"Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
 684
 685   fn1 = opt2fn_null("-s1",NFILE,fnm);
 686   fn2 = opt2fn_null("-s2",NFILE,fnm);
 687   if ((fn1 && fn2) || bCompAB) {
 688     if (bCompAB) {
 689       if (fn1 == NULL)
 690         gmx_fatal(FARGS,"With -ab you need to set the -s1 option");
 691       fn2 = NULL;
 692     }
 693     comp_tpx(fn1,fn2,bRMSD,ftol,abstol);
 694   } else if (fn1 && tex)
 695     tpx2methods(fn1,tex);
 696   else if ((fn1 && !opt2fn_null("-f",NFILE,fnm)) || (!fn1 && fn2)) {
 697     fprintf(stderr,"Please give me TWO run input (.tpr/.tpa/.tpb) files\n"
 698             "or specify the -m flag to generate a methods.tex file\n");
 699   }
 700
 701   fn1 = opt2fn_null("-e",NFILE,fnm);
 702   fn2 = opt2fn_null("-e2",NFILE,fnm);
 703   if (fn1 && fn2)
 704     comp_enx(fn1,fn2,ftol,abstol,lastener);
 705   else if (fn1)
 706     chk_enx(ftp2fn(efEDR,NFILE,fnm));
 707   else if (fn2)
 708     fprintf(stderr,"Please give me TWO energy (.edr/.ene) files!\n");
 709
 710   if (ftp2bSet(efTPS,NFILE,fnm))
 711     chk_tps(ftp2fn(efTPS,NFILE,fnm), vdw_fac, bon_lo, bon_hi);
 712
 713   if (ftp2bSet(efNDX,NFILE,fnm))
 714     chk_ndx(ftp2fn(efNDX,NFILE,fnm));
 715
 716   thanx(stderr);
 717
 718   return 0;
 719 }