2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2013, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/math/vec.h"
46 #include "gromacs/math/units.h"
47 #include "gromacs/topology/index.h"
49 #include "gromacs/utility/futil.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/trnio.h"
52 #include "gromacs/fileio/xtcio.h"
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/enxio.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
58 #include "gromacs/commandline/pargs.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/topology/atomprop.h"
61 #include "gromacs/topology/block.h"
62 #include "gromacs/topology/mtop_util.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.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;
286 real rdum, tt, old_t1, old_t2, prec;
287 gmx_bool bShowTimestep = TRUE, bOK, newline = FALSE;
290 gmx_localtop_t *top = NULL;
296 read_tpx_state(tpr, &ir, &state, NULL, &mtop);
297 top = gmx_mtop_generate_local_top(&mtop, &ir);
302 printf("Checking file %s\n", fn);
332 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
338 fprintf(stderr, "\n# Atoms %d\n", fr.natoms);
341 fprintf(stderr, "Precision %g (nm)\n", 1/fr.prec);
345 if ((natoms > 0) && (new_natoms != natoms))
347 fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
348 old_t1, natoms, new_natoms);
353 if (fabs((fr.time-old_t1)-(old_t1-old_t2)) >
354 0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) )
356 bShowTimestep = FALSE;
357 fprintf(stderr, "%sTimesteps at t=%g don't match (%g, %g)\n",
358 newline ? "\n" : "", old_t1, old_t1-old_t2, fr.time-old_t1);
364 chk_bonds(&top->idef, ir.ePBC, fr.x, fr.box, tol);
368 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
372 chk_vels(j, natoms, fr.v);
376 chk_forces(j, natoms, fr.f);
382 new_natoms = fr.natoms;
383 #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++; \
385 INC(fr, count, first, last, bStep);
386 INC(fr, count, first, last, bTime);
387 INC(fr, count, first, last, bLambda);
388 INC(fr, count, first, last, bX);
389 INC(fr, count, first, last, bV);
390 INC(fr, count, first, last, bF);
391 INC(fr, count, first, last, bBox);
394 while (read_next_frame(oenv, status, &fr));
396 fprintf(stderr, "\n");
400 fprintf(stderr, "\nItem #frames");
403 fprintf(stderr, " Timestep (ps)");
405 fprintf(stderr, "\n");
406 #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")
407 PRINTITEM ( "Step", bStep );
408 PRINTITEM ( "Time", bTime );
409 PRINTITEM ( "Lambda", bLambda );
410 PRINTITEM ( "Coords", bX );
411 PRINTITEM ( "Velocities", bV );
412 PRINTITEM ( "Forces", bF );
413 PRINTITEM ( "Box", bBox );
416 void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
427 gmx_bool bV, bX, bB, bFirst, bOut;
428 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
432 fprintf(stderr, "Checking coordinate file %s\n", fn);
433 read_tps_conf(fn, title, &top, &ePBC, &x, &v, box, TRUE);
436 fprintf(stderr, "%d atoms in file\n", atoms->nr);
438 /* check coordinates and box */
441 for (i = 0; (i < natom) && !(bV && bX); i++)
443 for (j = 0; (j < DIM) && !(bV && bX); j++)
445 bV = bV || (v[i][j] != 0);
446 bX = bX || (x[i][j] != 0);
450 for (i = 0; (i < DIM) && !bB; i++)
452 for (j = 0; (j < DIM) && !bB; j++)
454 bB = bB || (box[i][j] != 0);
458 fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
459 fprintf(stderr, "box %s\n", bB ? "found" : "absent");
460 fprintf(stderr, "velocities %s\n", bV ? "found" : "absent");
461 fprintf(stderr, "\n");
463 /* check velocities */
467 for (i = 0; (i < natom); i++)
469 for (j = 0; (j < DIM); j++)
471 ekin += 0.5*atoms->atom[i].m*v[i][j]*v[i][j];
474 temp1 = (2.0*ekin)/(natom*DIM*BOLTZ);
475 temp2 = (2.0*ekin)/(natom*(DIM-1)*BOLTZ);
476 fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
477 fprintf(stderr, "Assuming the number of degrees of freedom to be "
478 "Natoms * %d or Natoms * %d,\n"
479 "the velocities correspond to a temperature of the system\n"
480 "of %g K or %g K respectively.\n\n", DIM, DIM-1, temp1, temp2);
483 /* check coordinates */
486 vdwfac2 = sqr(vdw_fac);
487 bonlo2 = sqr(bon_lo);
488 bonhi2 = sqr(bon_hi);
491 "Checking for atoms closer than %g and not between %g and %g,\n"
492 "relative to sum of Van der Waals distance:\n",
493 vdw_fac, bon_lo, bon_hi);
494 snew(atom_vdw, natom);
495 aps = gmx_atomprop_init();
496 for (i = 0; (i < natom); i++)
498 gmx_atomprop_query(aps, epropVDW,
499 *(atoms->resinfo[atoms->atom[i].resind].name),
500 *(atoms->atomname[i]), &(atom_vdw[i]));
503 fprintf(debug, "%5d %4s %4s %7g\n", i+1,
504 *(atoms->resinfo[atoms->atom[i].resind].name),
505 *(atoms->atomname[i]),
509 gmx_atomprop_destroy(aps);
512 set_pbc(&pbc, ePBC, box);
516 for (i = 0; (i < natom); i++)
520 fprintf(stderr, "\r%5d", i+1);
522 for (j = i+1; (j < natom); j++)
526 pbc_dx(&pbc, x[i], x[j], dx);
530 rvec_sub(x[i], x[j], dx);
533 dist2 = sqr(atom_vdw[i]+atom_vdw[j]);
534 if ( (r2 <= dist2*bonlo2) ||
535 ( (r2 >= dist2*bonhi2) && (r2 <= dist2*vdwfac2) ) )
539 fprintf(stderr, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
540 "atom#", "name", "residue", "r_vdw",
541 "atom#", "name", "residue", "r_vdw", "distance");
545 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
546 i+1, *(atoms->atomname[i]),
547 *(atoms->resinfo[atoms->atom[i].resind].name),
548 atoms->resinfo[atoms->atom[i].resind].nr,
550 j+1, *(atoms->atomname[j]),
551 *(atoms->resinfo[atoms->atom[j].resind].name),
552 atoms->resinfo[atoms->atom[j].resind].nr,
560 fprintf(stderr, "\rno close atoms found\n");
562 fprintf(stderr, "\r \n");
569 for (i = 0; (i < natom) && (k < 10); i++)
572 for (j = 0; (j < DIM) && !bOut; j++)
574 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
581 fprintf(stderr, "Atoms outside box ( ");
582 for (j = 0; (j < DIM); j++)
584 fprintf(stderr, "%g ", box[j][j]);
586 fprintf(stderr, "):\n"
587 "(These may occur often and are normally not a problem)\n"
588 "%5s %4s %8s %5s %s\n",
589 "atom#", "name", "residue", "r_vdw", "coordinate");
593 "%5d %4s %4s%4d %-5.3g",
594 i, *(atoms->atomname[i]),
595 *(atoms->resinfo[atoms->atom[i].resind].name),
596 atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
597 for (j = 0; (j < DIM); j++)
599 fprintf(stderr, " %6.3g", x[i][j]);
601 fprintf(stderr, "\n");
606 fprintf(stderr, "(maybe more)\n");
610 fprintf(stderr, "no atoms found outside box\n");
612 fprintf(stderr, "\n");
617 void chk_ndx(const char *fn)
623 grps = init_index(fn, &grpname);
626 pr_blocka(debug, 0, fn, grps, FALSE);
630 printf("Contents of index file %s\n", fn);
631 printf("--------------------------------------------------\n");
632 printf("Nr. Group #Entries First Last\n");
633 for (i = 0; (i < grps->nr); i++)
635 printf("%4d %-20s%8d%8d%8d\n", i, grpname[i],
636 grps->index[i+1]-grps->index[i],
637 grps->a[grps->index[i]]+1,
638 grps->a[grps->index[i+1]-1]+1);
641 for (i = 0; (i < grps->nr); i++)
649 void chk_enx(const char *fn)
653 gmx_enxnm_t *enm = NULL;
656 real t0, old_t1, old_t2;
659 fprintf(stderr, "Checking energy file %s\n\n", fn);
661 in = open_enx(fn, "r");
662 do_enxnms(in, &nre, &enm);
663 fprintf(stderr, "%d groups in energy file", nre);
671 while (do_enx(in, fr))
675 if (fabs((fr->t-old_t1)-(old_t1-old_t2)) >
676 0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) )
679 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n",
680 old_t1, old_t1-old_t2, fr->t-old_t1);
691 fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n",
692 gmx_step_str(fr->step, buf), fnr, fr->t);
696 fprintf(stderr, "\n\nFound %d frames", fnr);
697 if (bShowTStep && fnr > 1)
699 fprintf(stderr, " with a timestep of %g ps", (old_t1-t0)/(fnr-1));
701 fprintf(stderr, ".\n");
704 free_enxnms(nre, enm);
708 int gmx_check(int argc, char *argv[])
710 const char *desc[] = {
711 "[THISMODULE] reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
712 "[TT].xtc[tt]), an energy file ([TT].ene[tt] or [TT].edr[tt])",
713 "or an index file ([TT].ndx[tt])",
714 "and prints out useful information about them.[PAR]",
715 "Option [TT]-c[tt] checks for presence of coordinates,",
716 "velocities and box in the file, for close contacts (smaller than",
717 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
718 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
719 "radii) and atoms outside the box (these may occur often and are",
720 "no problem). If velocities are present, an estimated temperature",
721 "will be calculated from them.[PAR]",
722 "If an index file, is given its contents will be summarized.[PAR]",
723 "If both a trajectory and a [TT].tpr[tt] file are given (with [TT]-s1[tt])",
724 "the program will check whether the bond lengths defined in the tpr",
725 "file are indeed correct in the trajectory. If not you may have",
726 "non-matching files due to e.g. deshuffling or due to problems with",
727 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
728 "The program can compare two run input ([TT].tpr[tt], [TT].tpb[tt] or",
729 "[TT].tpa[tt]) files",
730 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
731 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
732 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
733 "For free energy simulations the A and B state topology from one",
734 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
735 "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
736 "consisting of a rough outline for a methods section for a paper."
739 { efTRX, "-f", NULL, ffOPTRD },
740 { efTRX, "-f2", NULL, ffOPTRD },
741 { efTPX, "-s1", "top1", ffOPTRD },
742 { efTPX, "-s2", "top2", ffOPTRD },
743 { efTPS, "-c", NULL, ffOPTRD },
744 { efEDR, "-e", NULL, ffOPTRD },
745 { efEDR, "-e2", "ener2", ffOPTRD },
746 { efNDX, "-n", NULL, ffOPTRD },
747 { efTEX, "-m", NULL, ffOPTWR }
749 #define NFILE asize(fnm)
750 const char *fn1 = NULL, *fn2 = NULL, *tex = NULL;
753 static real vdw_fac = 0.8;
754 static real bon_lo = 0.4;
755 static real bon_hi = 0.7;
756 static gmx_bool bRMSD = FALSE;
757 static real ftol = 0.001;
758 static real abstol = 0.001;
759 static gmx_bool bCompAB = FALSE;
760 static char *lastener = NULL;
761 static t_pargs pa[] = {
762 { "-vdwfac", FALSE, etREAL, {&vdw_fac},
763 "Fraction of sum of VdW radii used as warning cutoff" },
764 { "-bonlo", FALSE, etREAL, {&bon_lo},
765 "Min. fract. of sum of VdW radii for bonded atoms" },
766 { "-bonhi", FALSE, etREAL, {&bon_hi},
767 "Max. fract. of sum of VdW radii for bonded atoms" },
768 { "-rmsd", FALSE, etBOOL, {&bRMSD},
769 "Print RMSD for x, v and f" },
770 { "-tol", FALSE, etREAL, {&ftol},
771 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
772 { "-abstol", FALSE, etREAL, {&abstol},
773 "Absolute tolerance, useful when sums are close to zero." },
774 { "-ab", FALSE, etBOOL, {&bCompAB},
775 "Compare the A and B topology from one file" },
776 { "-lastener", FALSE, etSTR, {&lastener},
777 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
780 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
781 asize(desc), desc, 0, NULL, &oenv))
786 fn1 = opt2fn_null("-f", NFILE, fnm);
787 fn2 = opt2fn_null("-f2", NFILE, fnm);
788 tex = opt2fn_null("-m", NFILE, fnm);
791 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
795 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
799 fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
802 fn1 = opt2fn_null("-s1", NFILE, fnm);
803 fn2 = opt2fn_null("-s2", NFILE, fnm);
804 if ((fn1 && fn2) || bCompAB)
810 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
814 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
818 tpx2methods(fn1, tex);
820 else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
822 fprintf(stderr, "Please give me TWO run input (.tpr/.tpa/.tpb) files\n"
823 "or specify the -m flag to generate a methods.tex file\n");
826 fn1 = opt2fn_null("-e", NFILE, fnm);
827 fn2 = opt2fn_null("-e2", NFILE, fnm);
830 comp_enx(fn1, fn2, ftol, abstol, lastener);
834 chk_enx(ftp2fn(efEDR, NFILE, fnm));
838 fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
841 if (ftp2bSet(efTPS, NFILE, fnm))
843 chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
846 if (ftp2bSet(efNDX, NFILE, fnm))
848 chk_ndx(ftp2fn(efNDX, NFILE, fnm));