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-2013, 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
45 #include "gromacs/fileio/futil.h"
49 #include "gmx_fatal.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/trnio.h"
52 #include "gromacs/fileio/xtcio.h"
59 #include "gromacs/fileio/confio.h"
60 #include "gromacs/fileio/enxio.h"
61 #include "gromacs/fileio/tpxio.h"
62 #include "gromacs/fileio/trxio.h"
64 #include "mtop_util.h"
88 static void tpx2system(FILE *fp, gmx_mtop_t *mtop)
90 int i, nmol, nvsite = 0;
91 gmx_mtop_atomloop_block_t aloop;
94 fprintf(fp, "\\subsection{Simulation system}\n");
95 aloop = gmx_mtop_atomloop_block_init(mtop);
96 while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
98 if (atom->ptype == eptVSite)
103 fprintf(fp, "A system of %d molecules (%d atoms) was simulated.\n",
104 mtop->mols.nr, mtop->natoms-nvsite);
107 fprintf(fp, "Virtual sites were used in some of the molecules.\n");
112 static void tpx2params(FILE *fp, t_inputrec *ir)
114 fprintf(fp, "\\subsection{Simulation settings}\n");
115 fprintf(fp, "A total of %g ns were simulated with a time step of %g fs.\n",
116 ir->nsteps*ir->delta_t*0.001, 1000*ir->delta_t);
117 fprintf(fp, "Neighbor searching was performed every %d steps.\n", ir->nstlist);
118 fprintf(fp, "The %s algorithm was used for electrostatic interactions.\n",
119 EELTYPE(ir->coulombtype));
120 fprintf(fp, "with a cut-off of %g nm.\n", ir->rcoulomb);
121 if (ir->coulombtype == eelPME)
123 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);
125 if (ir->rvdw > ir->rlist)
127 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);
131 fprintf(fp, "A single cut-off of %g was used for Van der Waals interactions.\n", ir->rlist);
135 fprintf(fp, "Temperature coupling was done with the %s algorithm.\n",
136 etcoupl_names[ir->etc]);
140 fprintf(fp, "Pressure coupling was done with the %s algorithm.\n",
141 epcoupl_names[ir->epc]);
146 static void tpx2methods(const char *tpx, const char *tex)
154 read_tpx_state(tpx, &ir, &state, NULL, &mtop);
155 fp = gmx_fio_fopen(tex, "w");
156 fprintf(fp, "\\section{Methods}\n");
157 tpx2system(fp, &mtop);
162 static void chk_coords(int frame, int natoms, rvec *x, matrix box, real fac, real tol)
168 for (i = 0; (i < natoms); i++)
170 for (j = 0; (j < DIM); j++)
172 if ((vol > 0) && (fabs(x[i][j]) > fac*box[j][j]))
174 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
178 if ((fabs(x[i][XX]) < tol) &&
179 (fabs(x[i][YY]) < tol) &&
180 (fabs(x[i][ZZ]) < tol))
187 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
192 static void chk_vels(int frame, int natoms, rvec *v)
196 for (i = 0; (i < natoms); i++)
198 for (j = 0; (j < DIM); j++)
200 if (fabs(v[i][j]) > 500)
202 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
209 static void chk_forces(int frame, int natoms, rvec *f)
213 for (i = 0; (i < natoms); i++)
215 for (j = 0; (j < DIM); j++)
217 if (fabs(f[i][j]) > 10000)
219 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
226 static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
228 int ftype, i, k, ai, aj, type;
229 real b0, blen, deviation, devtot;
234 set_pbc(&pbc, ePBC, box);
235 for (ftype = 0; (ftype < F_NRE); ftype++)
237 if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
239 for (k = 0; (k < idef->il[ftype].nr); )
241 type = idef->il[ftype].iatoms[k++];
242 ai = idef->il[ftype].iatoms[k++];
243 aj = idef->il[ftype].iatoms[k++];
248 b0 = idef->iparams[type].harmonic.rA;
251 b0 = sqrt(idef->iparams[type].harmonic.rA);
254 b0 = idef->iparams[type].morse.b0A;
257 b0 = idef->iparams[type].cubic.b0;
260 b0 = idef->iparams[type].constr.dA;
267 pbc_dx(&pbc, x[ai], x[aj], dx);
269 deviation = sqr(blen-b0);
270 if (sqrt(deviation/sqr(b0) > tol))
272 fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
280 void chk_trj(const output_env_t oenv, const char *fn, const char *tpr, real tol)
284 t_fr_time first, last;
285 int j = -1, new_natoms, natoms;
287 real rdum, tt, old_t1, old_t2, prec;
288 gmx_bool bShowTimestep = TRUE, bOK, newline = FALSE;
291 gmx_localtop_t *top = NULL;
297 read_tpx_state(tpr, &ir, &state, NULL, &mtop);
298 top = gmx_mtop_generate_local_top(&mtop, &ir);
303 printf("Checking file %s\n", fn);
334 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
340 fprintf(stderr, "\n# Atoms %d\n", fr.natoms);
343 fprintf(stderr, "Precision %g (nm)\n", 1/fr.prec);
347 if ((natoms > 0) && (new_natoms != natoms))
349 fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
350 old_t1, natoms, new_natoms);
355 if (fabs((fr.time-old_t1)-(old_t1-old_t2)) >
356 0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) )
358 bShowTimestep = FALSE;
359 fprintf(stderr, "%sTimesteps at t=%g don't match (%g, %g)\n",
360 newline ? "\n" : "", old_t1, old_t1-old_t2, fr.time-old_t1);
366 chk_bonds(&top->idef, ir.ePBC, fr.x, fr.box, tol);
370 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
374 chk_vels(j, natoms, fr.v);
378 chk_forces(j, natoms, fr.f);
383 /*if (fpos && ((j<10 || j%10==0)))
384 fprintf(stderr," byte: %10lu",(unsigned long)fpos);*/
386 new_natoms = fr.natoms;
387 #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++; \
389 INC(fr, count, first, last, bStep);
390 INC(fr, count, first, last, bTime);
391 INC(fr, count, first, last, bLambda);
392 INC(fr, count, first, last, bX);
393 INC(fr, count, first, last, bV);
394 INC(fr, count, first, last, bF);
395 INC(fr, count, first, last, bBox);
397 fpos = gmx_fio_ftell(trx_get_fileio(status));
399 while (read_next_frame(oenv, status, &fr));
401 fprintf(stderr, "\n");
405 fprintf(stderr, "\nItem #frames");
408 fprintf(stderr, " Timestep (ps)");
410 fprintf(stderr, "\n");
411 #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")
412 PRINTITEM ( "Step", bStep );
413 PRINTITEM ( "Time", bTime );
414 PRINTITEM ( "Lambda", bLambda );
415 PRINTITEM ( "Coords", bX );
416 PRINTITEM ( "Velocities", bV );
417 PRINTITEM ( "Forces", bF );
418 PRINTITEM ( "Box", bBox );
421 void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
432 gmx_bool bV, bX, bB, bFirst, bOut;
433 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
437 fprintf(stderr, "Checking coordinate file %s\n", fn);
438 read_tps_conf(fn, title, &top, &ePBC, &x, &v, box, TRUE);
441 fprintf(stderr, "%d atoms in file\n", atoms->nr);
443 /* check coordinates and box */
446 for (i = 0; (i < natom) && !(bV && bX); i++)
448 for (j = 0; (j < DIM) && !(bV && bX); j++)
450 bV = bV || (v[i][j] != 0);
451 bX = bX || (x[i][j] != 0);
455 for (i = 0; (i < DIM) && !bB; i++)
457 for (j = 0; (j < DIM) && !bB; j++)
459 bB = bB || (box[i][j] != 0);
463 fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
464 fprintf(stderr, "box %s\n", bB ? "found" : "absent");
465 fprintf(stderr, "velocities %s\n", bV ? "found" : "absent");
466 fprintf(stderr, "\n");
468 /* check velocities */
472 for (i = 0; (i < natom); i++)
474 for (j = 0; (j < DIM); j++)
476 ekin += 0.5*atoms->atom[i].m*v[i][j]*v[i][j];
479 temp1 = (2.0*ekin)/(natom*DIM*BOLTZ);
480 temp2 = (2.0*ekin)/(natom*(DIM-1)*BOLTZ);
481 fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
482 fprintf(stderr, "Assuming the number of degrees of freedom to be "
483 "Natoms * %d or Natoms * %d,\n"
484 "the velocities correspond to a temperature of the system\n"
485 "of %g K or %g K respectively.\n\n", DIM, DIM-1, temp1, temp2);
488 /* check coordinates */
491 vdwfac2 = sqr(vdw_fac);
492 bonlo2 = sqr(bon_lo);
493 bonhi2 = sqr(bon_hi);
496 "Checking for atoms closer than %g and not between %g and %g,\n"
497 "relative to sum of Van der Waals distance:\n",
498 vdw_fac, bon_lo, bon_hi);
499 snew(atom_vdw, natom);
500 aps = gmx_atomprop_init();
501 for (i = 0; (i < natom); i++)
503 gmx_atomprop_query(aps, epropVDW,
504 *(atoms->resinfo[atoms->atom[i].resind].name),
505 *(atoms->atomname[i]), &(atom_vdw[i]));
508 fprintf(debug, "%5d %4s %4s %7g\n", i+1,
509 *(atoms->resinfo[atoms->atom[i].resind].name),
510 *(atoms->atomname[i]),
514 gmx_atomprop_destroy(aps);
517 set_pbc(&pbc, ePBC, box);
521 for (i = 0; (i < natom); i++)
525 fprintf(stderr, "\r%5d", i+1);
527 for (j = i+1; (j < natom); j++)
531 pbc_dx(&pbc, x[i], x[j], dx);
535 rvec_sub(x[i], x[j], dx);
538 dist2 = sqr(atom_vdw[i]+atom_vdw[j]);
539 if ( (r2 <= dist2*bonlo2) ||
540 ( (r2 >= dist2*bonhi2) && (r2 <= dist2*vdwfac2) ) )
544 fprintf(stderr, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
545 "atom#", "name", "residue", "r_vdw",
546 "atom#", "name", "residue", "r_vdw", "distance");
550 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
551 i+1, *(atoms->atomname[i]),
552 *(atoms->resinfo[atoms->atom[i].resind].name),
553 atoms->resinfo[atoms->atom[i].resind].nr,
555 j+1, *(atoms->atomname[j]),
556 *(atoms->resinfo[atoms->atom[j].resind].name),
557 atoms->resinfo[atoms->atom[j].resind].nr,
565 fprintf(stderr, "\rno close atoms found\n");
567 fprintf(stderr, "\r \n");
574 for (i = 0; (i < natom) && (k < 10); i++)
577 for (j = 0; (j < DIM) && !bOut; j++)
579 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
586 fprintf(stderr, "Atoms outside box ( ");
587 for (j = 0; (j < DIM); j++)
589 fprintf(stderr, "%g ", box[j][j]);
591 fprintf(stderr, "):\n"
592 "(These may occur often and are normally not a problem)\n"
593 "%5s %4s %8s %5s %s\n",
594 "atom#", "name", "residue", "r_vdw", "coordinate");
598 "%5d %4s %4s%4d %-5.3g",
599 i, *(atoms->atomname[i]),
600 *(atoms->resinfo[atoms->atom[i].resind].name),
601 atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
602 for (j = 0; (j < DIM); j++)
604 fprintf(stderr, " %6.3g", x[i][j]);
606 fprintf(stderr, "\n");
611 fprintf(stderr, "(maybe more)\n");
615 fprintf(stderr, "no atoms found outside box\n");
617 fprintf(stderr, "\n");
622 void chk_ndx(const char *fn)
628 grps = init_index(fn, &grpname);
631 pr_blocka(debug, 0, fn, grps, FALSE);
635 printf("Contents of index file %s\n", fn);
636 printf("--------------------------------------------------\n");
637 printf("Nr. Group #Entries First Last\n");
638 for (i = 0; (i < grps->nr); i++)
640 printf("%4d %-20s%8d%8d%8d\n", i, grpname[i],
641 grps->index[i+1]-grps->index[i],
642 grps->a[grps->index[i]]+1,
643 grps->a[grps->index[i+1]-1]+1);
646 for (i = 0; (i < grps->nr); i++)
654 void chk_enx(const char *fn)
658 gmx_enxnm_t *enm = NULL;
661 real t0, old_t1, old_t2;
664 fprintf(stderr, "Checking energy file %s\n\n", fn);
666 in = open_enx(fn, "r");
667 do_enxnms(in, &nre, &enm);
668 fprintf(stderr, "%d groups in energy file", nre);
676 while (do_enx(in, fr))
680 if (fabs((fr->t-old_t1)-(old_t1-old_t2)) >
681 0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) )
684 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n",
685 old_t1, old_t1-old_t2, fr->t-old_t1);
696 fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n",
697 gmx_step_str(fr->step, buf), fnr, fr->t);
701 fprintf(stderr, "\n\nFound %d frames", fnr);
702 if (bShowTStep && fnr > 1)
704 fprintf(stderr, " with a timestep of %g ps", (old_t1-t0)/(fnr-1));
706 fprintf(stderr, ".\n");
709 free_enxnms(nre, enm);
713 int gmx_gmxcheck(int argc, char *argv[])
715 const char *desc[] = {
716 "[TT]gmx check[tt] reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
717 "[TT].xtc[tt]), an energy file ([TT].ene[tt] or [TT].edr[tt])",
718 "or an index file ([TT].ndx[tt])",
719 "and prints out useful information about them.[PAR]",
720 "Option [TT]-c[tt] checks for presence of coordinates,",
721 "velocities and box in the file, for close contacts (smaller than",
722 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
723 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
724 "radii) and atoms outside the box (these may occur often and are",
725 "no problem). If velocities are present, an estimated temperature",
726 "will be calculated from them.[PAR]",
727 "If an index file, is given its contents will be summarized.[PAR]",
728 "If both a trajectory and a [TT].tpr[tt] file are given (with [TT]-s1[tt])",
729 "the program will check whether the bond lengths defined in the tpr",
730 "file are indeed correct in the trajectory. If not you may have",
731 "non-matching files due to e.g. deshuffling or due to problems with",
732 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
733 "The program can compare two run input ([TT].tpr[tt], [TT].tpb[tt] or",
734 "[TT].tpa[tt]) files",
735 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
736 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
737 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
738 "For free energy simulations the A and B state topology from one",
739 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
740 "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
741 "consisting of a rough outline for a methods section for a paper."
744 { efTRX, "-f", NULL, ffOPTRD },
745 { efTRX, "-f2", NULL, ffOPTRD },
746 { efTPX, "-s1", "top1", ffOPTRD },
747 { efTPX, "-s2", "top2", ffOPTRD },
748 { efTPS, "-c", NULL, ffOPTRD },
749 { efEDR, "-e", NULL, ffOPTRD },
750 { efEDR, "-e2", "ener2", ffOPTRD },
751 { efNDX, "-n", NULL, ffOPTRD },
752 { efTEX, "-m", NULL, ffOPTWR }
754 #define NFILE asize(fnm)
755 const char *fn1 = NULL, *fn2 = NULL, *tex = NULL;
758 static real vdw_fac = 0.8;
759 static real bon_lo = 0.4;
760 static real bon_hi = 0.7;
761 static gmx_bool bRMSD = FALSE;
762 static real ftol = 0.001;
763 static real abstol = 0.001;
764 static gmx_bool bCompAB = FALSE;
765 static char *lastener = NULL;
766 static t_pargs pa[] = {
767 { "-vdwfac", FALSE, etREAL, {&vdw_fac},
768 "Fraction of sum of VdW radii used as warning cutoff" },
769 { "-bonlo", FALSE, etREAL, {&bon_lo},
770 "Min. fract. of sum of VdW radii for bonded atoms" },
771 { "-bonhi", FALSE, etREAL, {&bon_hi},
772 "Max. fract. of sum of VdW radii for bonded atoms" },
773 { "-rmsd", FALSE, etBOOL, {&bRMSD},
774 "Print RMSD for x, v and f" },
775 { "-tol", FALSE, etREAL, {&ftol},
776 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
777 { "-abstol", FALSE, etREAL, {&abstol},
778 "Absolute tolerance, useful when sums are close to zero." },
779 { "-ab", FALSE, etBOOL, {&bCompAB},
780 "Compare the A and B topology from one file" },
781 { "-lastener", FALSE, etSTR, {&lastener},
782 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
785 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
786 asize(desc), desc, 0, NULL, &oenv))
791 fn1 = opt2fn_null("-f", NFILE, fnm);
792 fn2 = opt2fn_null("-f2", NFILE, fnm);
793 tex = opt2fn_null("-m", NFILE, fnm);
796 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
800 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
804 fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
807 fn1 = opt2fn_null("-s1", NFILE, fnm);
808 fn2 = opt2fn_null("-s2", NFILE, fnm);
809 if ((fn1 && fn2) || bCompAB)
815 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
819 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
823 tpx2methods(fn1, tex);
825 else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
827 fprintf(stderr, "Please give me TWO run input (.tpr/.tpa/.tpb) files\n"
828 "or specify the -m flag to generate a methods.tex file\n");
831 fn1 = opt2fn_null("-e", NFILE, fnm);
832 fn2 = opt2fn_null("-e2", NFILE, fnm);
835 comp_enx(fn1, fn2, ftol, abstol, lastener);
839 chk_enx(ftp2fn(efEDR, NFILE, fnm));
843 fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
846 if (ftp2bSet(efTPS, NFILE, fnm))
848 chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
851 if (ftp2bSet(efNDX, NFILE, fnm))
853 chk_ndx(ftp2fn(efNDX, NFILE, fnm));