3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
50 #include "gmx_fatal.h"
65 #include "mtop_util.h"
87 static void tpx2system(FILE *fp,gmx_mtop_t *mtop)
90 gmx_mtop_atomloop_block_t aloop;
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)
99 fprintf(fp,"A system of %d molecules (%d atoms) was simulated.\n",
100 mtop->mols.nr,mtop->natoms-nvsite);
102 fprintf(fp,"Virtual sites were used in some of the molecules.\n");
106 static void tpx2params(FILE *fp,t_inputrec *ir)
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);
120 fprintf(fp,"A single cut-off of %g was used for Van der Waals interactions.\n",ir->rlist);
122 fprintf(fp,"Temperature coupling was done with the %s algorithm.\n",
123 etcoupl_names[ir->etc]);
126 fprintf(fp,"Pressure coupling was done with the %s algorithm.\n",
127 epcoupl_names[ir->epc]);
132 static void tpx2methods(const char *tpx,const char *tex)
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);
148 static void chk_coords(int frame,int natoms,rvec *x,matrix box,real fac,real tol)
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",
160 if ((fabs(x[i][XX]) < tol) &&
161 (fabs(x[i][YY]) < tol) &&
162 (fabs(x[i][ZZ]) < tol))
168 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
172 static void chk_vels(int frame,int natoms,rvec *v)
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",
185 static void chk_forces(int frame,int natoms,rvec *f)
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",
198 static void chk_bonds(t_idef *idef,int ePBC,rvec *x,matrix box,real tol)
200 int ftype,i,k,ai,aj,type;
201 real b0,blen,deviation,devtot;
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++];
217 b0 = idef->iparams[type].harmonic.rA;
220 b0 = idef->iparams[type].morse.b0A;
223 b0 = idef->iparams[type].cubic.b0;
226 b0 = idef->iparams[type].constr.dA;
232 pbc_dx(&pbc,x[ai],x[aj],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);
243 void chk_trj(const output_env_t oenv,const char *fn,const char *tpr,real tol)
247 t_fr_time first,last;
248 int j=-1,new_natoms,natoms;
250 real rdum,tt,old_t1,old_t2,prec;
251 gmx_bool bShowTimestep=TRUE,bOK,newline=FALSE;
254 gmx_localtop_t *top=NULL;
259 read_tpx_state(tpr,&ir,&state,NULL,&mtop);
260 top = gmx_mtop_generate_local_top(&mtop,&ir);
265 printf("Checking file %s\n",fn);
296 read_first_frame(oenv,&status,fn,&fr,TRX_READ_X | TRX_READ_V | TRX_READ_F);
300 fprintf(stderr,"\n# Atoms %d\n",fr.natoms);
302 fprintf(stderr,"Precision %g (nm)\n",1/fr.prec);
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);
311 if ( fabs((fr.time-old_t1)-(old_t1-old_t2)) >
312 0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) ) {
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);
320 chk_bonds(&top->idef,ir.ePBC,fr.x,fr.box,tol);
323 chk_coords(j,natoms,fr.x,fr.box,1e5,tol);
325 chk_vels(j,natoms,fr.v);
327 chk_forces(j,natoms,fr.f);
331 /*if (fpos && ((j<10 || j%10==0)))
332 fprintf(stderr," byte: %10lu",(unsigned long)fpos);*/
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);
344 fpos = gmx_fio_ftell(trx_get_fileio(status));
345 } while (read_next_frame(oenv,status,&fr));
347 fprintf(stderr,"\n");
351 fprintf(stderr,"\nItem #frames");
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 );
365 void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
376 gmx_bool bV,bX,bB,bFirst,bOut;
377 real r2,ekin,temp1,temp2,dist2,vdwfac2,bonlo2,bonhi2;
381 fprintf(stderr,"Checking coordinate file %s\n",fn);
382 read_tps_conf(fn,title,&top,&ePBC,&x,&v,box,TRUE);
385 fprintf(stderr,"%d atoms in file\n",atoms->nr);
387 /* check coordinates and box */
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);
396 for (i=0; (i<DIM) && !bB; i++)
397 for (j=0; (j<DIM) && !bB; j++)
398 bB=bB || (box[i][j]!=0);
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");
405 /* check velocities */
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);
420 /* check coordinates */
422 vdwfac2=sqr(vdw_fac);
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]),
441 gmx_atomprop_destroy(aps);
443 set_pbc(&pbc,ePBC,box);
446 for (i=0; (i<natom); i++) {
448 fprintf(stderr,"\r%5d",i+1);
449 for (j=i+1; (j<natom); j++) {
451 pbc_dx(&pbc,x[i],x[j],dx);
453 rvec_sub(x[i],x[j],dx);
455 dist2=sqr(atom_vdw[i]+atom_vdw[j]);
456 if ( (r2<=dist2*bonlo2) ||
457 ( (r2>=dist2*bonhi2) && (r2<=dist2*vdwfac2) ) ) {
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");
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,
470 j+1,*(atoms->atomname[j]),
471 *(atoms->resinfo[atoms->atom[j].resind].name),
472 atoms->resinfo[atoms->atom[j].resind].nr,
479 fprintf(stderr,"\rno close atoms found\n");
480 fprintf(stderr,"\r \n");
486 for (i=0; (i<natom) && (k<10); i++) {
488 for(j=0; (j<DIM) && !bOut; j++)
489 bOut = bOut || (x[i][j]<0) || (x[i][j]>box[j][j]);
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");
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");
513 fprintf(stderr,"(maybe more)\n");
515 fprintf(stderr,"no atoms found outside box\n");
516 fprintf(stderr,"\n");
521 void chk_ndx(const char *fn)
527 grps = init_index(fn,&grpname);
529 pr_blocka(debug,0,fn,grps,FALSE);
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);
541 for(i=0; (i<grps->nr); i++)
547 void chk_enx(const char *fn)
551 gmx_enxnm_t *enm=NULL;
554 real t0,old_t1,old_t2;
557 fprintf(stderr,"Checking energy file %s\n\n",fn);
559 in = open_enx(fn,"r");
560 do_enxnms(in,&nre,&enm);
561 fprintf(stderr,"%d groups in energy file",nre);
569 while (do_enx(in,fr)) {
571 if ( fabs((fr->t-old_t1)-(old_t1-old_t2)) >
572 0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) ) {
574 fprintf(stderr,"\nTimesteps at t=%g don't match (%g, %g)\n",
575 old_t1,old_t1-old_t2,fr->t-old_t1);
580 if (t0 == NOTSET) t0=fr->t;
582 fprintf(stderr,"\rframe: %6s (index %6d), t: %10.3f\n",
583 gmx_step_str(fr->step,buf),fnr,fr->t);
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");
592 free_enxnms(nre,enm);
596 int main(int argc,char *argv[])
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."
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 }
637 #define NFILE asize(fnm)
638 const char *fn1=NULL,*fn2=NULL,*tex=NULL;
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." }
668 CopyRight(stderr,argv[0]);
669 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
670 asize(desc),desc,0,NULL,&oenv);
672 fn1 = opt2fn_null("-f",NFILE,fnm);
673 fn2 = opt2fn_null("-f2",NFILE,fnm);
674 tex = opt2fn_null("-m",NFILE,fnm);
676 comp_trx(oenv,fn1,fn2,bRMSD,ftol,abstol);
678 chk_trj(oenv,fn1,opt2fn_null("-s1",NFILE,fnm),ftol);
680 fprintf(stderr,"Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
682 fn1 = opt2fn_null("-s1",NFILE,fnm);
683 fn2 = opt2fn_null("-s2",NFILE,fnm);
684 if ((fn1 && fn2) || bCompAB) {
687 gmx_fatal(FARGS,"With -ab you need to set the -s1 option");
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");
698 fn1 = opt2fn_null("-e",NFILE,fnm);
699 fn2 = opt2fn_null("-e2",NFILE,fnm);
701 comp_enx(fn1,fn2,ftol,abstol,lastener);
703 chk_enx(ftp2fn(efEDR,NFILE,fnm));
705 fprintf(stderr,"Please give me TWO energy (.edr/.ene) files!\n");
707 if (ftp2bSet(efTPS,NFILE,fnm))
708 chk_tps(ftp2fn(efTPS,NFILE,fnm), vdw_fac, bon_lo, bon_hi);
710 if (ftp2bSet(efNDX,NFILE,fnm))
711 chk_ndx(ftp2fn(efNDX,NFILE,fnm));