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