src/programs/gmxcheck/gmxcheck.c

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