src/kernel/gmxcheck.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #include <stdio.h>
  37 #include <string.h>
  38 #include <ctype.h>
  39 #include "main.h"
  40 #include "macros.h"
  41 #include "math.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "copyrite.h"
  45 #include "sysstuff.h"
  46 #include "txtdump.h"
  47 #include "fatal.h"
  48 #include "gmxfio.h"
  49 #include "trnio.h"
  50 #include "xtcio.h"
  51 #include "tpbcmp.h"
  52 #include "atomprop.h"
  53 #include "vec.h"
  54 #include "pbc.h"
  55 #include "physics.h"
  56 #include "smalloc.h"
  57 #include "confio.h"
  58 #include "enxio.h"
  59
  60 typedef struct {
  61   int bStep;
  62   int bTime;
  63   int bLambda;
  64   int bX;
  65   int bV;
  66   int bF;
  67   int bBox;
  68 } t_count;
  69
  70 typedef struct {
  71   float bStep;
  72   float bTime;
  73   float bLambda;
  74   float bX;
  75   float bV;
  76   float bF;
  77   float bBox;
  78 } t_fr_time;
  79
  80 void chk_trj(char *fn)
  81 {
  82   t_trxframe   fr;
  83   t_count      count;
  84   t_fr_time    first,last;
  85   int          j=-1,new_natoms,natoms;
  86   off_t        fpos;
  87   real         rdum,t,tt,old_t1,old_t2,prec;
  88   bool         bShowTimestep=TRUE,bOK,newline=FALSE;
  89   int          status;
  90
  91   new_natoms = -1;
  92   natoms = -1;
  93   t      = 0;
  94
  95   printf("Checking file %s\n",fn);
  96
  97   j      =  0;
  98   t      = -1;
  99   old_t2 = -2.0;
 100   old_t1 = -1.0;
 101   fpos   = 0;
 102
 103   count.bStep = 0;
 104   count.bTime = 0;
 105   count.bLambda = 0;
 106   count.bX = 0;
 107   count.bV = 0;
 108   count.bF = 0;
 109   count.bBox = 0;
 110
 111   first.bStep = 0;
 112   first.bTime = 0;
 113   first.bLambda = 0;
 114   first.bX = 0;
 115   first.bV = 0;
 116   first.bF = 0;
 117   first.bBox = 0;
 118
 119   last.bStep = 0;
 120   last.bTime = 0;
 121   last.bLambda = 0;
 122   last.bX = 0;
 123   last.bV = 0;
 124   last.bF = 0;
 125   last.bBox = 0;
 126
 127   read_first_frame(&status,fn,&fr,TRX_READ_X | TRX_READ_V | TRX_READ_F);
 128
 129   do {
 130     if (j == 0) {
 131       fprintf(stderr,"\n# Atoms  %d\n",fr.natoms);
 132       if (fr.bPrec)
 133         fprintf(stderr,"Precision %g (nm)\n",1/fr.prec);
 134     }
 135     newline=TRUE;
 136     if ((natoms > 0) && (new_natoms != natoms)) {
 137       fprintf(stderr,"\nNumber of atoms at t=%g don't match (%d, %d)\n",
 138               old_t1,natoms,new_natoms);
 139       newline=FALSE;
 140     }
 141     if (j>=2) {
 142       if ( fabs((fr.time-old_t1)-(old_t1-old_t2)) >
 143            0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) ) {
 144         bShowTimestep=FALSE;
 145         fprintf(stderr,"%sTimesteps at t=%g don't match (%g, %g)\n",
 146                 newline?"\n":"",old_t1,old_t1-old_t2,fr.time-old_t1);
 147       }
 148     }
 149     natoms=new_natoms;
 150     old_t2=old_t1;
 151     old_t1=fr.time;
 152     if (fpos && (j<10 || j%10==0))
 153       fprintf(stderr," byte: %10lu",(unsigned long)fpos);
 154     j++;
 155     t=fr.time;
 156     new_natoms=fr.natoms;
 157 #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++; }
 158     INC(fr,count,first,last,bStep);
 159     INC(fr,count,first,last,bTime);
 160     INC(fr,count,first,last,bLambda);
 161     INC(fr,count,first,last,bX);
 162     INC(fr,count,first,last,bV);
 163     INC(fr,count,first,last,bF);
 164     INC(fr,count,first,last,bBox);
 165 #undef INC
 166     fpos = fio_ftell(status);
 167   } while (read_next_frame(status,&fr));
 168
 169   fprintf(stderr,"\n");
 170
 171   close_trj(status);
 172
 173   fprintf(stderr,"\nItem        #frames");
 174   if (bShowTimestep)
 175     fprintf(stderr," Timestep (ps)");
 176   fprintf(stderr,"\n");
 177 #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")
 178   PRINTITEM ( "Step",       bStep );
 179   PRINTITEM ( "Time",       bTime );
 180   PRINTITEM ( "Lambda",     bLambda );
 181   PRINTITEM ( "Coords",     bX );
 182   PRINTITEM ( "Velocities", bV );
 183   PRINTITEM ( "Forces",     bF );
 184   PRINTITEM ( "Box",        bBox );
 185 }
 186
 187 void chk_tps(char *fn, real vdw_fac, real bon_lo, real bon_hi)
 188 {
 189   int       natom,i,j,k;
 190   char      title[STRLEN];
 191   t_topology top;
 192   t_atoms   *atoms;
 193   rvec      *x,*v;
 194   rvec      dx;
 195   matrix    box;
 196   bool      bV,bX,bB,bFirst,bOut;
 197   real      r2,ekin,temp1,temp2,dist2,vdwfac2,bonlo2,bonhi2;
 198   real      *atom_vdw;
 199   void      *ap;
 200
 201   fprintf(stderr,"Checking coordinate file %s\n",fn);
 202   read_tps_conf(fn,title,&top,&x,&v,box,TRUE);
 203   atoms=&top.atoms;
 204   natom=atoms->nr;
 205   fprintf(stderr,"%d atoms in file\n",atoms->nr);
 206
 207   /* check coordinates and box */
 208   bV=FALSE;
 209   bX=FALSE;
 210   for (i=0; (i<natom) && !(bV && bX); i++)
 211     for (j=0; (j<DIM) && !(bV && bX); j++) {
 212       bV=bV || (v[i][j]!=0);
 213       bX=bX || (x[i][j]!=0);
 214     }
 215   bB=FALSE;
 216   for (i=0; (i<DIM) && !bB; i++)
 217     for (j=0; (j<DIM) && !bB; j++)
 218       bB=bB || (box[i][j]!=0);
 219
 220   fprintf(stderr,"coordinates %s\n",bX?"found":"absent");
 221   fprintf(stderr,"box         %s\n",bB?"found":"absent");
 222   fprintf(stderr,"velocities  %s\n",bV?"found":"absent");
 223   fprintf(stderr,"\n");
 224
 225   /* check velocities */
 226   if (bV) {
 227     ekin=0.0;
 228     for (i=0; (i<natom); i++)
 229       for (j=0; (j<DIM); j++)
 230         ekin+=0.5*atoms->atom[i].m*v[i][j]*v[i][j];
 231     temp1=(2.0*ekin)/(natom*DIM*BOLTZ);
 232     temp2=(2.0*ekin)/(natom*(DIM-1)*BOLTZ);
 233     fprintf(stderr,"Kinetic energy: %g (kJ/mol)\n",ekin);
 234     fprintf(stderr,"Assuming the number of degrees of freedom to be "
 235             "Natoms * %d or Natoms * %d,\n"
 236             "the velocities correspond to a temperature of the system\n"
 237             "of %g K or %g K respectively.\n\n",DIM,DIM-1,temp1,temp2);
 238   }
 239
 240   /* check coordinates */
 241   if (bX) {
 242     vdwfac2=sqr(vdw_fac);
 243     bonlo2=sqr(bon_lo);
 244     bonhi2=sqr(bon_hi);
 245
 246     fprintf(stderr,
 247             "Checking for atoms closer than %g and not between %g and %g,\n"
 248             "relative to sum of Van der Waals distance:\n",
 249             vdw_fac,bon_lo,bon_hi);
 250     snew(atom_vdw,natom);
 251     ap = get_atomprop();
 252     for (i=0; (i<natom); i++) {
 253       (void) query_atomprop(ap,epropVDW,*(atoms->resname[atoms->atom[i].resnr]),
 254                             *(atoms->atomname[i]),&(atom_vdw[i]));
 255       if (debug) fprintf(debug,"%5d %4s %4s %7g\n",i+1,
 256                          *(atoms->resname[atoms->atom[i].resnr]),
 257                          *(atoms->atomname[i]),
 258                          atom_vdw[i]);
 259     }
 260     done_atomprop(&ap);
 261     if (bB)
 262       init_pbc(box);
 263
 264     bFirst=TRUE;
 265     for (i=0; (i<natom); i++) {
 266       if (((i+1)%10)==0)
 267         fprintf(stderr,"\r%5d",i+1);
 268       for (j=i+1; (j<natom); j++) {
 269         if (bB)
 270           pbc_dx(x[i],x[j],dx);
 271         else
 272           rvec_sub(x[i],x[j],dx);
 273         r2=iprod(dx,dx);
 274         dist2=sqr(atom_vdw[i]+atom_vdw[j]);
 275         if ( (r2<=dist2*bonlo2) ||
 276              ( (r2>=dist2*bonhi2) && (r2<=dist2*vdwfac2) ) ) {
 277           if (bFirst) {
 278             fprintf(stderr,"\r%5s %4s %8s %5s  %5s %4s %8s %5s  %6s\n",
 279                     "atom#","name","residue","r_vdw",
 280                     "atom#","name","residue","r_vdw","distance");
 281             bFirst=FALSE;
 282           }
 283           fprintf(stderr,
 284                   "\r%5d %4s %4s%4d %-5.3g  %5d %4s %4s%4d %-5.3g  %-6.4g\n",
 285                   i+1,*(atoms->atomname[i]),
 286                   *(atoms->resname[atoms->atom[i].resnr]),atoms->atom[i].resnr+1,
 287                   atom_vdw[i],
 288                   j+1,*(atoms->atomname[j]),
 289                   *(atoms->resname[atoms->atom[j].resnr]),atoms->atom[j].resnr+1,
 290                   atom_vdw[j],
 291                   sqrt(r2) );
 292         }
 293       }
 294     }
 295     if (bFirst)
 296       fprintf(stderr,"\rno close atoms found\n");
 297     fprintf(stderr,"\r      \n");
 298
 299     if (bB) {
 300       /* check box */
 301       bFirst=TRUE;
 302       k=0;
 303       for (i=0; (i<natom) && (k<10); i++) {
 304         bOut=FALSE;
 305         for(j=0; (j<DIM) && !bOut; j++)
 306           bOut = bOut || (x[i][j]<0) || (x[i][j]>box[j][j]);
 307         if (bOut) {
 308           k++;
 309           if (bFirst) {
 310             fprintf(stderr,"Atoms outside box ( ");
 311             for (j=0; (j<DIM); j++)
 312               fprintf(stderr,"%g ",box[j][j]);
 313             fprintf(stderr,"):\n"
 314                     "(These may occur often and are normally not a problem)\n"
 315                     "%5s %4s %8s %5s  %s\n",
 316                     "atom#","name","residue","r_vdw","coordinate");
 317             bFirst=FALSE;
 318           }
 319           fprintf(stderr,
 320                   "%5d %4s %4s%4d %-5.3g",
 321                   i,*(atoms->atomname[i]),
 322                   *(atoms->resname[atoms->atom[i].resnr]),
 323                   atoms->atom[i].resnr,atom_vdw[i]);
 324           for (j=0; (j<DIM); j++)
 325             fprintf(stderr," %6.3g",x[i][j]);
 326           fprintf(stderr,"\n");
 327         }
 328       }
 329       if (k==10)
 330         fprintf(stderr,"(maybe more)\n");
 331       if (bFirst)
 332         fprintf(stderr,"no atoms found outside box\n");
 333       fprintf(stderr,"\n");
 334     }
 335   }
 336 }
 337
 338 void chk_enx(char *fn)
 339 {
 340   int        in,nre,fnr,ndr;
 341   char       **enm=NULL;
 342   t_enxframe *fr;
 343   bool       bShowTStep;
 344   real       t0,old_t1,old_t2;
 345
 346   fprintf(stderr,"Checking energy file %s\n\n",fn);
 347
 348   in = open_enx(fn,"r");
 349   do_enxnms(in,&nre,&enm);
 350   fprintf(stderr,"%d groups in energy file",nre);
 351   snew(fr,1);
 352   old_t2=-2.0;
 353   old_t1=-1.0;
 354   fnr=0;
 355   t0=NOTSET;
 356   bShowTStep=TRUE;
 357
 358   while (do_enx(in,fr)) {
 359     if (fnr>=2) {
 360       if ( fabs((fr->t-old_t1)-(old_t1-old_t2)) >
 361            0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) ) {
 362         bShowTStep=FALSE;
 363         fprintf(stderr,"\nTimesteps at t=%g don't match (%g, %g)\n",
 364                 old_t1,old_t1-old_t2,fr->t-old_t1);
 365       }
 366     }
 367     old_t2=old_t1;
 368     old_t1=fr->t;
 369     if (t0 == NOTSET) t0=fr->t;
 370     if (fnr == 0)
 371       fprintf(stderr,"\rframe: %6d (index %6d), t: %10.3f\n",
 372               fr->step,fnr,fr->t);
 373     fnr++;
 374   }
 375   fprintf(stderr,"\n\nFound %d frames",fnr);
 376   if (bShowTStep && fnr>1)
 377     fprintf(stderr," with a timestep of %g ps",(old_t1-t0)/(fnr-1));
 378   fprintf(stderr,".\n");
 379
 380   free_enxframe(fr);
 381   sfree(fr);
 382 }
 383
 384
 385 int main(int argc,char *argv[])
 386 {
 387   static char *desc[] = {
 388     "gmxcheck reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
 389     "[TT].xtc[tt]) or an energy file ([TT].ene[tt] or [TT].edr[tt])",
 390     "and prints out useful information about them.[PAR]",
 391     "Option [TT]-c[tt] checks for presence of coordinates,",
 392     "velocities and box in the file, for close contacts (smaller than",
 393     "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
 394     "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
 395     "radii) and atoms outside the box (these may occur often and are",
 396     "no problem). If velocities are present, an estimated temperature",
 397     "will be calculated from them.[PAR]",
 398     "The program will compare run input ([TT].tpr[tt], [TT].tpb[tt] or",
 399     "[TT].tpa[tt]) files",
 400     "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
 401     "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
 402     "option), or a pair of energy files (using the [TT]-e2[tt] option)."
 403   };
 404   t_filenm fnm[] = {
 405     { efTRX, "-f",  NULL, ffOPTRD },
 406     { efTRX, "-f2",  NULL, ffOPTRD },
 407     { efTPX, "-s1", "top1", ffOPTRD },
 408     { efTPX, "-s2", "top2", ffOPTRD },
 409     { efTPS, "-c",  NULL, ffOPTRD },
 410     { efENX, "-e",  NULL, ffOPTRD },
 411     { efENX, "-e2", "ener2", ffOPTRD }
 412   };
 413 #define NFILE asize(fnm)
 414   char *fn1=NULL,*fn2=NULL;
 415
 416   static real vdw_fac=0.8;
 417   static real bon_lo=0.4;
 418   static real bon_hi=0.7;
 419   static real ftol=0;
 420   static char *lastener=NULL;
 421   static t_pargs pa[] = {
 422     { "-vdwfac", FALSE, etREAL, {&vdw_fac},
 423       "Fraction of sum of VdW radii used as warning cutoff" },
 424     { "-bonlo",  FALSE, etREAL, {&bon_lo},
 425       "Min. fract. of sum of VdW radii for bonded atoms" },
 426     { "-bonhi",  FALSE, etREAL, {&bon_hi},
 427       "Max. fract. of sum of VdW radii for bonded atoms" },
 428     { "-tol",    FALSE, etREAL, {&ftol},
 429       "Relative tolerance for comparing real values defined as 2*(a-b)/(|a|+|b|)" },
 430     { "-lastener",FALSE, etSTR,  {&lastener},
 431       "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
 432   };
 433
 434   CopyRight(stdout,argv[0]);
 435   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 436                     asize(desc),desc,0,NULL);
 437
 438   fn1 = opt2fn_null("-f",NFILE,fnm);
 439   fn2 = opt2fn_null("-f2",NFILE,fnm);
 440   if (fn1 && fn2)
 441     comp_trx(fn1,fn2,ftol);
 442   else if (fn1)
 443     chk_trj(fn1);
 444   else if (fn2)
 445     fprintf(stderr,"Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
 446
 447   fn1 = opt2fn_null("-s1",NFILE,fnm);
 448   fn2 = opt2fn_null("-s2",NFILE,fnm);
 449   if (fn1 && fn2)
 450     comp_tpx(fn1,fn2,ftol);
 451   else if (fn1 || fn2)
 452     fprintf(stderr,"Please give me TWO run input (.tpr/.tpa/.tpb) files!\n");
 453
 454   fn1 = opt2fn_null("-e",NFILE,fnm);
 455   fn2 = opt2fn_null("-e2",NFILE,fnm);
 456   if (fn1 && fn2)
 457     comp_enx(fn1,fn2,ftol,lastener);
 458   else if (fn1)
 459     chk_enx(ftp2fn(efENX,NFILE,fnm));
 460   else if (fn2)
 461     fprintf(stderr,"Please give me TWO energy (.edr/.ene) files!\n");
 462
 463   if (ftp2bSet(efTPS,NFILE,fnm))
 464     chk_tps(ftp2fn(efTPS,NFILE,fnm), vdw_fac, bon_lo, bon_hi);
 465
 466   if (ftp2bSet(efENX,NFILE,fnm))
 467
 468   thanx(stderr);
 469
 470   return 0;
 471 }