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