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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/fileio/xtcio.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdrun/mdmodules.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/mdtypes/md_enums.h"
59 #include "gromacs/mdtypes/state.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/atomprop.h"
62 #include "gromacs/topology/block.h"
63 #include "gromacs/topology/ifunc.h"
64 #include "gromacs/topology/index.h"
65 #include "gromacs/topology/mtop_util.h"
66 #include "gromacs/topology/topology.h"
67 #include "gromacs/trajectory/energyframe.h"
68 #include "gromacs/trajectory/trajectoryframe.h"
69 #include "gromacs/utility/arraysize.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/futil.h"
72 #include "gromacs/utility/smalloc.h"
96 static void comp_tpx(const char* fn1, const char* fn2, gmx_bool bRMSD, real ftol, real abstol)
106 for (i = 0; i < (fn2 ? 2 : 1); i++)
108 ir[i] = new t_inputrec();
109 read_tpx_state(ff[i], ir[i], &state[i], &(mtop[i]));
110 gmx::MDModules().adjustInputrecBasedOnModules(ir[i]);
114 cmp_inputrec(stdout, ir[0], ir[1], ftol, abstol);
115 compareMtop(stdout, mtop[0], mtop[1], ftol, abstol);
116 comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
120 if (ir[0]->efep == FreeEnergyPerturbationType::No)
122 fprintf(stdout, "inputrec->efep = %s\n", enumValueToString(ir[0]->efep));
128 comp_pull_AB(stdout, *ir[0]->pull, ftol, abstol);
130 compareMtopAB(stdout, mtop[0], ftol, abstol);
137 static void comp_trx(const gmx_output_env_t* oenv, const char* fn1, const char* fn2, gmx_bool bRMSD, real ftol, real abstol)
142 t_trxstatus* status[2];
147 fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
148 for (i = 0; i < 2; i++)
150 b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X | TRX_READ_V | TRX_READ_F);
157 comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
159 for (i = 0; i < 2; i++)
161 b[i] = read_next_frame(oenv, status[i], &fr[i]);
163 } while (b[0] && b[1]);
165 for (i = 0; i < 2; i++)
167 if (b[i] && !b[1 - i])
169 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1 - i], fn[i]);
171 close_trx(status[i]);
176 fprintf(stdout, "\nBoth files read correctly\n");
180 static void chk_coords(int frame, int natoms, rvec* x, matrix box, real fac, real tol)
186 for (i = 0; (i < natoms); i++)
188 for (j = 0; (j < DIM); j++)
190 if ((vol > 0) && (fabs(x[i][j]) > fac * box[j][j]))
192 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n", frame, i, x[i][j]);
195 if ((fabs(x[i][XX]) < tol) && (fabs(x[i][YY]) < tol) && (fabs(x[i][ZZ]) < tol))
202 printf("Warning at frame %d: there are %d particles with all coordinates zero\n", frame, nNul);
206 static void chk_vels(int frame, int natoms, rvec* v)
210 for (i = 0; (i < natoms); i++)
212 for (j = 0; (j < DIM); j++)
214 if (fabs(v[i][j]) > 500)
216 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n", frame, i, v[i][j]);
222 static void chk_forces(int frame, int natoms, rvec* f)
226 for (i = 0; (i < natoms); i++)
228 for (j = 0; (j < DIM); j++)
230 if (fabs(f[i][j]) > 10000)
232 printf("Warning at frame %d. Forces for atom %d are large (%g)\n", frame, i, f[i][j]);
238 static void chk_bonds(const InteractionDefinitions* idef, PbcType pbcType, rvec* x, matrix box, real tol)
240 int ftype, k, ai, aj, type;
241 real b0, blen, deviation;
245 gmx::ArrayRef<const t_iparams> iparams = idef->iparams;
247 set_pbc(&pbc, pbcType, box);
248 for (ftype = 0; (ftype < F_NRE); ftype++)
250 if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
252 for (k = 0; (k < idef->il[ftype].size());)
254 type = idef->il[ftype].iatoms[k++];
255 ai = idef->il[ftype].iatoms[k++];
256 aj = idef->il[ftype].iatoms[k++];
260 case F_BONDS: b0 = iparams[type].harmonic.rA; break;
261 case F_G96BONDS: b0 = std::sqrt(iparams[type].harmonic.rA); break;
262 case F_MORSE: b0 = iparams[type].morse.b0A; break;
263 case F_CUBICBONDS: b0 = iparams[type].cubic.b0; break;
264 case F_CONSTR: b0 = iparams[type].constr.dA; break;
269 pbc_dx(&pbc, x[ai], x[aj], dx);
271 deviation = gmx::square(blen - b0);
272 if (std::sqrt(deviation / gmx::square(b0)) > tol)
275 "Distance between atoms %d and %d is %.3f, should be %.3f\n",
287 static void chk_trj(const gmx_output_env_t* oenv, const char* fn, const char* tpr, real tol)
291 t_fr_time first, last;
292 int j = -1, new_natoms, natoms;
294 gmx_bool bShowTimestep = TRUE, newline = FALSE;
300 std::unique_ptr<gmx_localtop_t> top;
303 read_tpx_state(tpr, &ir, &state, &mtop);
304 top = std::make_unique<gmx_localtop_t>(mtop.ffparams);
305 gmx_mtop_generate_local_top(mtop, top.get(), ir.efep != FreeEnergyPerturbationType::No);
310 printf("Checking file %s\n", fn);
340 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
346 fprintf(stderr, "\n# Atoms %d\n", fr.natoms);
349 fprintf(stderr, "Precision %g (nm)\n", 1 / fr.prec);
353 if ((natoms > 0) && (new_natoms != natoms))
355 fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n", old_t1, natoms, new_natoms);
360 if (std::fabs((fr.time - old_t1) - (old_t1 - old_t2))
361 > 0.1 * (std::fabs(fr.time - old_t1) + std::fabs(old_t1 - old_t2)))
363 bShowTimestep = FALSE;
365 "%sTimesteps at t=%g don't match (%g, %g)\n",
375 chk_bonds(&top->idef, ir.pbcType, fr.x, fr.box, tol);
379 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
383 chk_vels(j, natoms, fr.v);
387 chk_forces(j, natoms, fr.f);
393 new_natoms = fr.natoms;
394 #define INC(s, n, f, l, item) \
399 first.item = fr.time; \
401 last.item = fr.time; \
404 INC(fr, count, first, last, bStep)
405 INC(fr, count, first, last, bTime)
406 INC(fr, count, first, last, bLambda)
407 INC(fr, count, first, last, bX)
408 INC(fr, count, first, last, bV)
409 INC(fr, count, first, last, bF)
410 INC(fr, count, first, last, bBox)
412 } while (read_next_frame(oenv, status, &fr));
414 fprintf(stderr, "\n");
418 fprintf(stderr, "\nItem #frames");
421 fprintf(stderr, " Timestep (ps)");
423 fprintf(stderr, "\n");
424 #define PRINTITEM(label, item) \
425 fprintf(stderr, "%-10s %6d", label, count.item); \
426 if ((bShowTimestep) && (count.item > 1)) \
428 fprintf(stderr, " %g\n", (last.item - first.item) / (count.item - 1)); \
431 fprintf(stderr, "\n")
432 PRINTITEM("Step", bStep);
433 PRINTITEM("Time", bTime);
434 PRINTITEM("Lambda", bLambda);
435 PRINTITEM("Coords", bX);
436 PRINTITEM("Velocities", bV);
437 PRINTITEM("Forces", bF);
438 PRINTITEM("Box", bBox);
441 static void chk_tps(const char* fn, real vdw_fac, real bon_lo, real bon_hi)
451 gmx_bool bV, bX, bB, bFirst, bOut;
452 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
455 fprintf(stderr, "Checking coordinate file %s\n", fn);
456 read_tps_conf(fn, &top, &pbcType, &x, &v, box, TRUE);
459 fprintf(stderr, "%d atoms in file\n", atoms->nr);
461 /* check coordinates and box */
464 for (i = 0; (i < natom) && !(bV && bX); i++)
466 for (j = 0; (j < DIM) && !(bV && bX); j++)
468 bV = bV || (v[i][j] != 0);
469 bX = bX || (x[i][j] != 0);
473 for (i = 0; (i < DIM) && !bB; i++)
475 for (j = 0; (j < DIM) && !bB; j++)
477 bB = bB || (box[i][j] != 0);
481 fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
482 fprintf(stderr, "box %s\n", bB ? "found" : "absent");
483 fprintf(stderr, "velocities %s\n", bV ? "found" : "absent");
484 fprintf(stderr, "\n");
486 /* check velocities */
490 for (i = 0; (i < natom); i++)
492 for (j = 0; (j < DIM); j++)
494 ekin += 0.5 * atoms->atom[i].m * v[i][j] * v[i][j];
497 temp1 = (2.0 * ekin) / (natom * DIM * gmx::c_boltz);
498 temp2 = (2.0 * ekin) / (natom * (DIM - 1) * gmx::c_boltz);
499 fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
501 "Assuming the number of degrees of freedom to be "
502 "Natoms * %d or Natoms * %d,\n"
503 "the velocities correspond to a temperature of the system\n"
504 "of %g K or %g K respectively.\n\n",
511 /* check coordinates */
514 vdwfac2 = gmx::square(vdw_fac);
515 bonlo2 = gmx::square(bon_lo);
516 bonhi2 = gmx::square(bon_hi);
519 "Checking for atoms closer than %g and not between %g and %g,\n"
520 "relative to sum of Van der Waals distance:\n",
524 snew(atom_vdw, natom);
526 for (i = 0; (i < natom); i++)
528 aps.setAtomProperty(epropVDW,
529 *(atoms->resinfo[atoms->atom[i].resind].name),
530 *(atoms->atomname[i]),
537 *(atoms->resinfo[atoms->atom[i].resind].name),
538 *(atoms->atomname[i]),
544 set_pbc(&pbc, pbcType, box);
548 for (i = 0; (i < natom); i++)
550 if (((i + 1) % 10) == 0)
552 fprintf(stderr, "\r%5d", i + 1);
555 for (j = i + 1; (j < natom); j++)
559 pbc_dx(&pbc, x[i], x[j], dx);
563 rvec_sub(x[i], x[j], dx);
566 dist2 = gmx::square(atom_vdw[i] + atom_vdw[j]);
567 if ((r2 <= dist2 * bonlo2) || ((r2 >= dist2 * bonhi2) && (r2 <= dist2 * vdwfac2)))
572 "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
585 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
587 *(atoms->atomname[i]),
588 *(atoms->resinfo[atoms->atom[i].resind].name),
589 atoms->resinfo[atoms->atom[i].resind].nr,
592 *(atoms->atomname[j]),
593 *(atoms->resinfo[atoms->atom[j].resind].name),
594 atoms->resinfo[atoms->atom[j].resind].nr,
602 fprintf(stderr, "\rno close atoms found\n");
604 fprintf(stderr, "\r \n");
611 for (i = 0; (i < natom) && (k < 10); i++)
614 for (j = 0; (j < DIM) && !bOut; j++)
616 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
623 fprintf(stderr, "Atoms outside box ( ");
624 for (j = 0; (j < DIM); j++)
626 fprintf(stderr, "%g ", box[j][j]);
630 "(These may occur often and are normally not a problem)\n"
631 "%5s %4s %8s %5s %s\n",
640 "%5d %4s %4s%4d %-5.3g",
642 *(atoms->atomname[i]),
643 *(atoms->resinfo[atoms->atom[i].resind].name),
644 atoms->resinfo[atoms->atom[i].resind].nr,
646 for (j = 0; (j < DIM); j++)
648 fprintf(stderr, " %6.3g", x[i][j]);
650 fprintf(stderr, "\n");
655 fprintf(stderr, "(maybe more)\n");
659 fprintf(stderr, "no atoms found outside box\n");
661 fprintf(stderr, "\n");
666 static void chk_ndx(const char* fn)
672 grps = init_index(fn, &grpname);
675 pr_blocka(debug, 0, fn, grps, FALSE);
679 printf("Contents of index file %s\n", fn);
680 printf("--------------------------------------------------\n");
681 printf("Nr. Group #Entries First Last\n");
682 for (i = 0; (i < grps->nr); i++)
684 printf("%4d %-20s%8d%8d%8d\n",
687 grps->index[i + 1] - grps->index[i],
688 grps->a[grps->index[i]] + 1,
689 grps->a[grps->index[i + 1] - 1] + 1);
692 for (i = 0; (i < grps->nr); i++)
700 static void chk_enx(const char* fn)
704 gmx_enxnm_t* enm = nullptr;
708 real t0, old_t1, old_t2;
711 fprintf(stderr, "Checking energy file %s\n\n", fn);
713 in = open_enx(fn, "r");
714 do_enxnms(in, &nre, &enm);
715 fprintf(stderr, "%d groups in energy file", nre);
724 while (do_enx(in, fr))
728 if (fabs((fr->t - old_t1) - (old_t1 - old_t2))
729 > 0.1 * (fabs(fr->t - old_t1) + std::fabs(old_t1 - old_t2)))
732 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n", old_t1, old_t1 - old_t2, fr->t - old_t1);
744 fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n", gmx_step_str(fr->step, buf), fnr, fr->t);
748 fprintf(stderr, "\n\nFound %d frames", fnr);
749 if (bShowTStep && fnr > 1)
751 fprintf(stderr, " with a timestep of %g ps", (old_t1 - t0) / (fnr - 1));
753 fprintf(stderr, ".\n");
756 free_enxnms(nre, enm);
760 int gmx_check(int argc, char* argv[])
762 const char* desc[] = {
763 "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
764 "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
765 "or an index file ([REF].ndx[ref])",
766 "and prints out useful information about them.[PAR]",
767 "Option [TT]-c[tt] checks for presence of coordinates,",
768 "velocities and box in the file, for close contacts (smaller than",
769 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
770 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
771 "radii) and atoms outside the box (these may occur often and are",
772 "no problem). If velocities are present, an estimated temperature",
773 "will be calculated from them.[PAR]",
774 "If an index file, is given its contents will be summarized.[PAR]",
775 "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
776 "the program will check whether the bond lengths defined in the tpr",
777 "file are indeed correct in the trajectory. If not you may have",
778 "non-matching files due to e.g. deshuffling or due to problems with",
779 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for ",
780 "such problems.[PAR]",
781 "The program can compare two run input ([REF].tpr[ref])",
783 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied. When comparing",
784 "run input files this way, the default relative tolerance is reduced",
785 "to 0.000001 and the absolute tolerance set to zero to find any differences",
786 "not due to minor compiler optimization differences, although you can",
787 "of course still set any other tolerances through the options.",
788 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
789 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
790 "For free energy simulations the A and B state topology from one",
791 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]"
793 t_filenm fnm[] = { { efTRX, "-f", nullptr, ffOPTRD }, { efTRX, "-f2", nullptr, ffOPTRD },
794 { efTPR, "-s1", "top1", ffOPTRD }, { efTPR, "-s2", "top2", ffOPTRD },
795 { efTPS, "-c", nullptr, ffOPTRD }, { efEDR, "-e", nullptr, ffOPTRD },
796 { efEDR, "-e2", "ener2", ffOPTRD }, { efNDX, "-n", nullptr, ffOPTRD },
797 { efTEX, "-m", nullptr, ffOPTWR } };
798 #define NFILE asize(fnm)
799 const char *fn1 = nullptr, *fn2 = nullptr, *tex = nullptr;
801 gmx_output_env_t* oenv;
805 gmx_bool bRMSD = FALSE;
808 gmx_bool bCompAB = FALSE;
809 char* lastener = nullptr;
810 static t_pargs pa[] = {
815 "Fraction of sum of VdW radii used as warning cutoff" },
816 { "-bonlo", FALSE, etREAL, { &bon_lo }, "Min. fract. of sum of VdW radii for bonded atoms" },
817 { "-bonhi", FALSE, etREAL, { &bon_hi }, "Max. fract. of sum of VdW radii for bonded atoms" },
818 { "-rmsd", FALSE, etBOOL, { &bRMSD }, "Print RMSD for x, v and f" },
823 "Relative tolerance for comparing real values defined as "
824 "[MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
829 "Absolute tolerance, useful when sums are close to zero." },
830 { "-ab", FALSE, etBOOL, { &bCompAB }, "Compare the A and B topology from one file" },
835 "Last energy term to compare (if not given all are tested). It makes sense to go up "
836 "until the Pressure." }
839 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
844 fn1 = opt2fn_null("-f", NFILE, fnm);
845 fn2 = opt2fn_null("-f2", NFILE, fnm);
846 tex = opt2fn_null("-m", NFILE, fnm);
851 "LaTeX file writing has been removed from gmx check. "
852 "Please use gmx report-methods instead for it.\n");
856 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
860 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
864 fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
866 output_env_done(oenv);
868 fn1 = opt2fn_null("-s1", NFILE, fnm);
869 fn2 = opt2fn_null("-s2", NFILE, fnm);
870 if ((fn1 && fn2) || bCompAB)
876 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
881 fprintf(stderr, "Note: When comparing run input files, default tolerances are reduced.\n");
882 if (!opt2parg_bSet("-tol", asize(pa), pa))
886 if (!opt2parg_bSet("-abstol", asize(pa), pa))
890 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
892 else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
894 fprintf(stderr, "Please give me TWO run input (.tpr) files\n");
897 fn1 = opt2fn_null("-e", NFILE, fnm);
898 fn2 = opt2fn_null("-e2", NFILE, fnm);
901 comp_enx(fn1, fn2, ftol, abstol, lastener);
905 chk_enx(ftp2fn(efEDR, NFILE, fnm));
909 fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
912 if (ftp2bSet(efTPS, NFILE, fnm))
914 chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
917 if (ftp2bSet(efNDX, NFILE, fnm))
919 chk_ndx(ftp2fn(efNDX, NFILE, fnm));