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,2015,2016,2017,2018,2019, 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/commandline/pargs.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/enxio.h"
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/fileio/tpxio.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/fileio/xtcio.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/mdrunutility/mdmodules.h"
56 #include "gromacs/mdtypes/inputrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/topology/atomprop.h"
61 #include "gromacs/topology/block.h"
62 #include "gromacs/topology/ifunc.h"
63 #include "gromacs/topology/index.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/trajectory/energyframe.h"
67 #include "gromacs/trajectory/trajectoryframe.h"
68 #include "gromacs/utility/arraysize.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/futil.h"
71 #include "gromacs/utility/smalloc.h"
93 static void comp_tpx(const char *fn1, const char *fn2,
94 gmx_bool bRMSD, real ftol, real abstol)
104 for (i = 0; i < (fn2 ? 2 : 1); i++)
106 ir[i] = new t_inputrec();
107 read_tpx_state(ff[i], ir[i], &state[i], &(mtop[i]));
108 gmx::MDModules().adjustInputrecBasedOnModules(ir[i]);
112 cmp_inputrec(stdout, ir[0], ir[1], ftol, abstol);
113 compareMtop(stdout, mtop[0], mtop[1], ftol, abstol);
114 comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
118 if (ir[0]->efep == efepNO)
120 fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0]->efep]);
126 comp_pull_AB(stdout, ir[0]->pull, ftol, abstol);
128 compareMtopAB(stdout, mtop[0], ftol, abstol);
133 static void comp_trx(const gmx_output_env_t *oenv, const char *fn1, const char *fn2,
134 gmx_bool bRMSD, real ftol, real abstol)
139 t_trxstatus *status[2];
144 fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
145 for (i = 0; i < 2; i++)
147 b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X|TRX_READ_V|TRX_READ_F);
154 comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
156 for (i = 0; i < 2; i++)
158 b[i] = read_next_frame(oenv, status[i], &fr[i]);
161 while (b[0] && b[1]);
163 for (i = 0; i < 2; i++)
167 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1-i], fn[i]);
169 close_trx(status[i]);
174 fprintf(stdout, "\nBoth files read correctly\n");
178 static void chk_coords(int frame, int natoms, rvec *x, matrix box, real fac, real tol)
184 for (i = 0; (i < natoms); i++)
186 for (j = 0; (j < DIM); j++)
188 if ((vol > 0) && (fabs(x[i][j]) > fac*box[j][j]))
190 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
194 if ((fabs(x[i][XX]) < tol) &&
195 (fabs(x[i][YY]) < tol) &&
196 (fabs(x[i][ZZ]) < tol))
203 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
208 static void chk_vels(int frame, int natoms, rvec *v)
212 for (i = 0; (i < natoms); i++)
214 for (j = 0; (j < DIM); j++)
216 if (fabs(v[i][j]) > 500)
218 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
225 static void chk_forces(int frame, int natoms, rvec *f)
229 for (i = 0; (i < natoms); i++)
231 for (j = 0; (j < DIM); j++)
233 if (fabs(f[i][j]) > 10000)
235 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
242 static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
244 int ftype, k, ai, aj, type;
245 real b0, blen, deviation;
249 set_pbc(&pbc, ePBC, box);
250 for (ftype = 0; (ftype < F_NRE); ftype++)
252 if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
254 for (k = 0; (k < idef->il[ftype].nr); )
256 type = idef->il[ftype].iatoms[k++];
257 ai = idef->il[ftype].iatoms[k++];
258 aj = idef->il[ftype].iatoms[k++];
263 b0 = idef->iparams[type].harmonic.rA;
266 b0 = std::sqrt(idef->iparams[type].harmonic.rA);
269 b0 = idef->iparams[type].morse.b0A;
272 b0 = idef->iparams[type].cubic.b0;
275 b0 = idef->iparams[type].constr.dA;
282 pbc_dx(&pbc, x[ai], x[aj], dx);
284 deviation = gmx::square(blen-b0);
285 if (std::sqrt(deviation/gmx::square(b0)) > tol)
287 fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
295 static void chk_trj(const gmx_output_env_t *oenv, const char *fn, const char *tpr, real tol)
299 t_fr_time first, last;
300 int j = -1, new_natoms, natoms;
302 gmx_bool bShowTimestep = TRUE, newline = FALSE;
311 read_tpx_state(tpr, &ir, &state, &mtop);
312 gmx_mtop_generate_local_top(mtop, &top, ir.efep != efepNO);
317 printf("Checking file %s\n", fn);
347 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
353 fprintf(stderr, "\n# Atoms %d\n", fr.natoms);
356 fprintf(stderr, "Precision %g (nm)\n", 1/fr.prec);
360 if ((natoms > 0) && (new_natoms != natoms))
362 fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
363 old_t1, natoms, new_natoms);
368 if (std::fabs((fr.time-old_t1)-(old_t1-old_t2)) >
369 0.1*(std::fabs(fr.time-old_t1)+std::fabs(old_t1-old_t2)) )
371 bShowTimestep = FALSE;
372 fprintf(stderr, "%sTimesteps at t=%g don't match (%g, %g)\n",
373 newline ? "\n" : "", old_t1, old_t1-old_t2, fr.time-old_t1);
379 chk_bonds(&top.idef, ir.ePBC, fr.x, fr.box, tol);
383 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
387 chk_vels(j, natoms, fr.v);
391 chk_forces(j, natoms, fr.f);
397 new_natoms = fr.natoms;
398 #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++; \
400 INC(fr, count, first, last, bStep);
401 INC(fr, count, first, last, bTime);
402 INC(fr, count, first, last, bLambda);
403 INC(fr, count, first, last, bX);
404 INC(fr, count, first, last, bV);
405 INC(fr, count, first, last, bF);
406 INC(fr, count, first, last, bBox);
409 while (read_next_frame(oenv, status, &fr));
411 fprintf(stderr, "\n");
415 fprintf(stderr, "\nItem #frames");
418 fprintf(stderr, " Timestep (ps)");
420 fprintf(stderr, "\n");
421 #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")
422 PRINTITEM ( "Step", bStep );
423 PRINTITEM ( "Time", bTime );
424 PRINTITEM ( "Lambda", bLambda );
425 PRINTITEM ( "Coords", bX );
426 PRINTITEM ( "Velocities", bV );
427 PRINTITEM ( "Forces", bF );
428 PRINTITEM ( "Box", bBox );
431 static void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
441 gmx_bool bV, bX, bB, bFirst, bOut;
442 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
445 fprintf(stderr, "Checking coordinate file %s\n", fn);
446 read_tps_conf(fn, &top, &ePBC, &x, &v, box, TRUE);
449 fprintf(stderr, "%d atoms in file\n", atoms->nr);
451 /* check coordinates and box */
454 for (i = 0; (i < natom) && !(bV && bX); i++)
456 for (j = 0; (j < DIM) && !(bV && bX); j++)
458 bV = bV || (v[i][j] != 0);
459 bX = bX || (x[i][j] != 0);
463 for (i = 0; (i < DIM) && !bB; i++)
465 for (j = 0; (j < DIM) && !bB; j++)
467 bB = bB || (box[i][j] != 0);
471 fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
472 fprintf(stderr, "box %s\n", bB ? "found" : "absent");
473 fprintf(stderr, "velocities %s\n", bV ? "found" : "absent");
474 fprintf(stderr, "\n");
476 /* check velocities */
480 for (i = 0; (i < natom); i++)
482 for (j = 0; (j < DIM); j++)
484 ekin += 0.5*atoms->atom[i].m*v[i][j]*v[i][j];
487 temp1 = (2.0*ekin)/(natom*DIM*BOLTZ);
488 temp2 = (2.0*ekin)/(natom*(DIM-1)*BOLTZ);
489 fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
490 fprintf(stderr, "Assuming the number of degrees of freedom to be "
491 "Natoms * %d or Natoms * %d,\n"
492 "the velocities correspond to a temperature of the system\n"
493 "of %g K or %g K respectively.\n\n", DIM, DIM-1, temp1, temp2);
496 /* check coordinates */
499 vdwfac2 = gmx::square(vdw_fac);
500 bonlo2 = gmx::square(bon_lo);
501 bonhi2 = gmx::square(bon_hi);
504 "Checking for atoms closer than %g and not between %g and %g,\n"
505 "relative to sum of Van der Waals distance:\n",
506 vdw_fac, bon_lo, bon_hi);
507 snew(atom_vdw, natom);
509 for (i = 0; (i < natom); i++)
511 aps.setAtomProperty(epropVDW,
512 *(atoms->resinfo[atoms->atom[i].resind].name),
513 *(atoms->atomname[i]), &(atom_vdw[i]));
516 fprintf(debug, "%5d %4s %4s %7g\n", i+1,
517 *(atoms->resinfo[atoms->atom[i].resind].name),
518 *(atoms->atomname[i]),
524 set_pbc(&pbc, ePBC, box);
528 for (i = 0; (i < natom); i++)
532 fprintf(stderr, "\r%5d", i+1);
535 for (j = i+1; (j < natom); j++)
539 pbc_dx(&pbc, x[i], x[j], dx);
543 rvec_sub(x[i], x[j], dx);
546 dist2 = gmx::square(atom_vdw[i]+atom_vdw[j]);
547 if ( (r2 <= dist2*bonlo2) ||
548 ( (r2 >= dist2*bonhi2) && (r2 <= dist2*vdwfac2) ) )
552 fprintf(stderr, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
553 "atom#", "name", "residue", "r_vdw",
554 "atom#", "name", "residue", "r_vdw", "distance");
558 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
559 i+1, *(atoms->atomname[i]),
560 *(atoms->resinfo[atoms->atom[i].resind].name),
561 atoms->resinfo[atoms->atom[i].resind].nr,
563 j+1, *(atoms->atomname[j]),
564 *(atoms->resinfo[atoms->atom[j].resind].name),
565 atoms->resinfo[atoms->atom[j].resind].nr,
573 fprintf(stderr, "\rno close atoms found\n");
575 fprintf(stderr, "\r \n");
582 for (i = 0; (i < natom) && (k < 10); i++)
585 for (j = 0; (j < DIM) && !bOut; j++)
587 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
594 fprintf(stderr, "Atoms outside box ( ");
595 for (j = 0; (j < DIM); j++)
597 fprintf(stderr, "%g ", box[j][j]);
599 fprintf(stderr, "):\n"
600 "(These may occur often and are normally not a problem)\n"
601 "%5s %4s %8s %5s %s\n",
602 "atom#", "name", "residue", "r_vdw", "coordinate");
606 "%5d %4s %4s%4d %-5.3g",
607 i, *(atoms->atomname[i]),
608 *(atoms->resinfo[atoms->atom[i].resind].name),
609 atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
610 for (j = 0; (j < DIM); j++)
612 fprintf(stderr, " %6.3g", x[i][j]);
614 fprintf(stderr, "\n");
619 fprintf(stderr, "(maybe more)\n");
623 fprintf(stderr, "no atoms found outside box\n");
625 fprintf(stderr, "\n");
630 static void chk_ndx(const char *fn)
636 grps = init_index(fn, &grpname);
639 pr_blocka(debug, 0, fn, grps, FALSE);
643 printf("Contents of index file %s\n", fn);
644 printf("--------------------------------------------------\n");
645 printf("Nr. Group #Entries First Last\n");
646 for (i = 0; (i < grps->nr); i++)
648 printf("%4d %-20s%8d%8d%8d\n", i, grpname[i],
649 grps->index[i+1]-grps->index[i],
650 grps->a[grps->index[i]]+1,
651 grps->a[grps->index[i+1]-1]+1);
654 for (i = 0; (i < grps->nr); i++)
662 static void chk_enx(const char *fn)
666 gmx_enxnm_t *enm = nullptr;
670 real t0, old_t1, old_t2;
673 fprintf(stderr, "Checking energy file %s\n\n", fn);
675 in = open_enx(fn, "r");
676 do_enxnms(in, &nre, &enm);
677 fprintf(stderr, "%d groups in energy file", nre);
686 while (do_enx(in, fr))
690 if (fabs((fr->t-old_t1)-(old_t1-old_t2)) >
691 0.1*(fabs(fr->t-old_t1)+std::fabs(old_t1-old_t2)) )
694 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n",
695 old_t1, old_t1-old_t2, fr->t-old_t1);
707 fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n",
708 gmx_step_str(fr->step, buf), fnr, fr->t);
712 fprintf(stderr, "\n\nFound %d frames", fnr);
713 if (bShowTStep && fnr > 1)
715 fprintf(stderr, " with a timestep of %g ps", (old_t1-t0)/(fnr-1));
717 fprintf(stderr, ".\n");
720 free_enxnms(nre, enm);
724 int gmx_check(int argc, char *argv[])
726 const char *desc[] = {
727 "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
728 "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
729 "or an index file ([REF].ndx[ref])",
730 "and prints out useful information about them.[PAR]",
731 "Option [TT]-c[tt] checks for presence of coordinates,",
732 "velocities and box in the file, for close contacts (smaller than",
733 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
734 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
735 "radii) and atoms outside the box (these may occur often and are",
736 "no problem). If velocities are present, an estimated temperature",
737 "will be calculated from them.[PAR]",
738 "If an index file, is given its contents will be summarized.[PAR]",
739 "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
740 "the program will check whether the bond lengths defined in the tpr",
741 "file are indeed correct in the trajectory. If not you may have",
742 "non-matching files due to e.g. deshuffling or due to problems with",
743 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
744 "The program can compare two run input ([REF].tpr[ref])",
746 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied. When comparing",
747 "run input files this way, the default relative tolerance is reduced",
748 "to 0.000001 and the absolute tolerance set to zero to find any differences",
749 "not due to minor compiler optimization differences, although you can",
750 "of course still set any other tolerances through the options.",
751 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
752 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
753 "For free energy simulations the A and B state topology from one",
754 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]"
757 { efTRX, "-f", nullptr, ffOPTRD },
758 { efTRX, "-f2", nullptr, ffOPTRD },
759 { efTPR, "-s1", "top1", ffOPTRD },
760 { efTPR, "-s2", "top2", ffOPTRD },
761 { efTPS, "-c", nullptr, ffOPTRD },
762 { efEDR, "-e", nullptr, ffOPTRD },
763 { efEDR, "-e2", "ener2", ffOPTRD },
764 { efNDX, "-n", nullptr, ffOPTRD },
765 { efTEX, "-m", nullptr, ffOPTWR }
767 #define NFILE asize(fnm)
768 const char *fn1 = nullptr, *fn2 = nullptr, *tex = nullptr;
770 gmx_output_env_t *oenv;
771 static real vdw_fac = 0.8;
772 static real bon_lo = 0.4;
773 static real bon_hi = 0.7;
774 static gmx_bool bRMSD = FALSE;
775 static real ftol = 0.001;
776 static real abstol = 0.001;
777 static gmx_bool bCompAB = FALSE;
778 static char *lastener = nullptr;
779 static t_pargs pa[] = {
780 { "-vdwfac", FALSE, etREAL, {&vdw_fac},
781 "Fraction of sum of VdW radii used as warning cutoff" },
782 { "-bonlo", FALSE, etREAL, {&bon_lo},
783 "Min. fract. of sum of VdW radii for bonded atoms" },
784 { "-bonhi", FALSE, etREAL, {&bon_hi},
785 "Max. fract. of sum of VdW radii for bonded atoms" },
786 { "-rmsd", FALSE, etBOOL, {&bRMSD},
787 "Print RMSD for x, v and f" },
788 { "-tol", FALSE, etREAL, {&ftol},
789 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
790 { "-abstol", FALSE, etREAL, {&abstol},
791 "Absolute tolerance, useful when sums are close to zero." },
792 { "-ab", FALSE, etBOOL, {&bCompAB},
793 "Compare the A and B topology from one file" },
794 { "-lastener", FALSE, etSTR, {&lastener},
795 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
798 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
799 asize(desc), desc, 0, nullptr, &oenv))
804 fn1 = opt2fn_null("-f", NFILE, fnm);
805 fn2 = opt2fn_null("-f2", NFILE, fnm);
806 tex = opt2fn_null("-m", NFILE, fnm);
810 fprintf(stderr, "LaTeX file writing has been removed from gmx check. "
811 "Please use gmx report-methods instead for it.\n");
815 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
819 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
823 fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
826 fn1 = opt2fn_null("-s1", NFILE, fnm);
827 fn2 = opt2fn_null("-s2", NFILE, fnm);
828 if ((fn1 && fn2) || bCompAB)
834 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
839 fprintf(stderr, "Note: When comparing run input files, default tolerances are reduced.\n");
840 if (!opt2parg_bSet("-tol", asize(pa), pa))
844 if (!opt2parg_bSet("-abstol", asize(pa), pa))
848 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
850 else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
852 fprintf(stderr, "Please give me TWO run input (.tpr) files\n");
855 fn1 = opt2fn_null("-e", NFILE, fnm);
856 fn2 = opt2fn_null("-e2", NFILE, fnm);
859 comp_enx(fn1, fn2, ftol, abstol, lastener);
863 chk_enx(ftp2fn(efEDR, NFILE, fnm));
867 fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
870 if (ftp2bSet(efTPS, NFILE, fnm))
872 chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
875 if (ftp2bSet(efNDX, NFILE, fnm))
877 chk_ndx(ftp2fn(efNDX, NFILE, fnm));