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, 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/trajectoryframe.h"
67 #include "gromacs/utility/arraysize.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/futil.h"
70 #include "gromacs/utility/smalloc.h"
92 static void comp_tpx(const char *fn1, const char *fn2,
93 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 /* Convert gmx_mtop_t to t_topology.
114 * We should implement direct mtop comparison,
115 * but it might be useful to keep t_topology comparison as an option.
117 top[0] = gmx_mtop_t_to_t_topology(&mtop[0], false);
118 top[1] = gmx_mtop_t_to_t_topology(&mtop[1], false);
119 cmp_top(stdout, &top[0], &top[1], ftol, abstol);
120 cmp_groups(stdout, &mtop[0].groups, &mtop[1].groups,
121 mtop[0].natoms, mtop[1].natoms);
122 comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
126 if (ir[0]->efep == efepNO)
128 fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0]->efep]);
134 comp_pull_AB(stdout, ir[0]->pull, ftol, abstol);
136 /* Convert gmx_mtop_t to t_topology.
137 * We should implement direct mtop comparison,
138 * but it might be useful to keep t_topology comparison as an option.
140 top[0] = gmx_mtop_t_to_t_topology(&mtop[0], true);
141 cmp_top(stdout, &top[0], nullptr, ftol, abstol);
146 static void comp_trx(const gmx_output_env_t *oenv, const char *fn1, const char *fn2,
147 gmx_bool bRMSD, real ftol, real abstol)
152 t_trxstatus *status[2];
157 fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
158 for (i = 0; i < 2; i++)
160 b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X|TRX_READ_V|TRX_READ_F);
167 comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
169 for (i = 0; i < 2; i++)
171 b[i] = read_next_frame(oenv, status[i], &fr[i]);
174 while (b[0] && b[1]);
176 for (i = 0; i < 2; i++)
180 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1-i], fn[i]);
182 close_trx(status[i]);
187 fprintf(stdout, "\nBoth files read correctly\n");
191 static void tpx2system(FILE *fp, const gmx_mtop_t *mtop)
193 int nmol, nvsite = 0;
194 gmx_mtop_atomloop_block_t aloop;
197 fprintf(fp, "\\subsection{Simulation system}\n");
198 aloop = gmx_mtop_atomloop_block_init(mtop);
199 while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
201 if (atom->ptype == eptVSite)
206 fprintf(fp, "A system of %d molecules (%d atoms) was simulated.\n",
207 mtop->mols.nr, mtop->natoms-nvsite);
210 fprintf(fp, "Virtual sites were used in some of the molecules.\n");
215 static void tpx2params(FILE *fp, const t_inputrec *ir)
217 fprintf(fp, "\\subsection{Simulation settings}\n");
218 fprintf(fp, "A total of %g ns were simulated with a time step of %g fs.\n",
219 ir->nsteps*ir->delta_t*0.001, 1000*ir->delta_t);
220 fprintf(fp, "Neighbor searching was performed every %d steps.\n", ir->nstlist);
221 fprintf(fp, "The %s algorithm was used for electrostatic interactions.\n",
222 EELTYPE(ir->coulombtype));
223 fprintf(fp, "with a cut-off of %g nm.\n", ir->rcoulomb);
224 if (ir->coulombtype == eelPME)
226 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);
228 if (ir->rvdw > ir->rlist)
230 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);
234 fprintf(fp, "A single cut-off of %g was used for Van der Waals interactions.\n", ir->rlist);
238 fprintf(fp, "Temperature coupling was done with the %s algorithm.\n",
239 etcoupl_names[ir->etc]);
243 fprintf(fp, "Pressure coupling was done with the %s algorithm.\n",
244 epcoupl_names[ir->epc]);
249 static void tpx2methods(const char *tpx, const char *tex)
256 read_tpx_state(tpx, &ir, &state, &mtop);
257 fp = gmx_fio_fopen(tex, "w");
258 fprintf(fp, "\\section{Methods}\n");
259 tpx2system(fp, &mtop);
264 static void chk_coords(int frame, int natoms, rvec *x, matrix box, real fac, real tol)
270 for (i = 0; (i < natoms); i++)
272 for (j = 0; (j < DIM); j++)
274 if ((vol > 0) && (fabs(x[i][j]) > fac*box[j][j]))
276 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
280 if ((fabs(x[i][XX]) < tol) &&
281 (fabs(x[i][YY]) < tol) &&
282 (fabs(x[i][ZZ]) < tol))
289 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
294 static void chk_vels(int frame, int natoms, rvec *v)
298 for (i = 0; (i < natoms); i++)
300 for (j = 0; (j < DIM); j++)
302 if (fabs(v[i][j]) > 500)
304 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
311 static void chk_forces(int frame, int natoms, rvec *f)
315 for (i = 0; (i < natoms); i++)
317 for (j = 0; (j < DIM); j++)
319 if (fabs(f[i][j]) > 10000)
321 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
328 static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
330 int ftype, k, ai, aj, type;
331 real b0, blen, deviation;
335 set_pbc(&pbc, ePBC, box);
336 for (ftype = 0; (ftype < F_NRE); ftype++)
338 if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
340 for (k = 0; (k < idef->il[ftype].nr); )
342 type = idef->il[ftype].iatoms[k++];
343 ai = idef->il[ftype].iatoms[k++];
344 aj = idef->il[ftype].iatoms[k++];
349 b0 = idef->iparams[type].harmonic.rA;
352 b0 = std::sqrt(idef->iparams[type].harmonic.rA);
355 b0 = idef->iparams[type].morse.b0A;
358 b0 = idef->iparams[type].cubic.b0;
361 b0 = idef->iparams[type].constr.dA;
368 pbc_dx(&pbc, x[ai], x[aj], dx);
370 deviation = gmx::square(blen-b0);
371 if (std::sqrt(deviation/gmx::square(b0)) > tol)
373 fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
381 static void chk_trj(const gmx_output_env_t *oenv, const char *fn, const char *tpr, real tol)
385 t_fr_time first, last;
386 int j = -1, new_natoms, natoms;
388 gmx_bool bShowTimestep = TRUE, newline = FALSE;
391 gmx_localtop_t *top = nullptr;
397 read_tpx_state(tpr, &ir, &state, &mtop);
398 top = gmx_mtop_generate_local_top(&mtop, ir.efep != efepNO);
403 printf("Checking file %s\n", fn);
433 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X | TRX_READ_V | TRX_READ_F);
439 fprintf(stderr, "\n# Atoms %d\n", fr.natoms);
442 fprintf(stderr, "Precision %g (nm)\n", 1/fr.prec);
446 if ((natoms > 0) && (new_natoms != natoms))
448 fprintf(stderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
449 old_t1, natoms, new_natoms);
454 if (fabs((fr.time-old_t1)-(old_t1-old_t2)) >
455 0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) )
457 bShowTimestep = FALSE;
458 fprintf(stderr, "%sTimesteps at t=%g don't match (%g, %g)\n",
459 newline ? "\n" : "", old_t1, old_t1-old_t2, fr.time-old_t1);
465 chk_bonds(&top->idef, ir.ePBC, fr.x, fr.box, tol);
469 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
473 chk_vels(j, natoms, fr.v);
477 chk_forces(j, natoms, fr.f);
483 new_natoms = fr.natoms;
484 #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++; \
486 INC(fr, count, first, last, bStep);
487 INC(fr, count, first, last, bTime);
488 INC(fr, count, first, last, bLambda);
489 INC(fr, count, first, last, bX);
490 INC(fr, count, first, last, bV);
491 INC(fr, count, first, last, bF);
492 INC(fr, count, first, last, bBox);
495 while (read_next_frame(oenv, status, &fr));
497 fprintf(stderr, "\n");
501 fprintf(stderr, "\nItem #frames");
504 fprintf(stderr, " Timestep (ps)");
506 fprintf(stderr, "\n");
507 #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")
508 PRINTITEM ( "Step", bStep );
509 PRINTITEM ( "Time", bTime );
510 PRINTITEM ( "Lambda", bLambda );
511 PRINTITEM ( "Coords", bX );
512 PRINTITEM ( "Velocities", bV );
513 PRINTITEM ( "Forces", bF );
514 PRINTITEM ( "Box", bBox );
517 static void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
527 gmx_bool bV, bX, bB, bFirst, bOut;
528 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
532 fprintf(stderr, "Checking coordinate file %s\n", fn);
533 read_tps_conf(fn, &top, &ePBC, &x, &v, box, TRUE);
536 fprintf(stderr, "%d atoms in file\n", atoms->nr);
538 /* check coordinates and box */
541 for (i = 0; (i < natom) && !(bV && bX); i++)
543 for (j = 0; (j < DIM) && !(bV && bX); j++)
545 bV = bV || (v[i][j] != 0);
546 bX = bX || (x[i][j] != 0);
550 for (i = 0; (i < DIM) && !bB; i++)
552 for (j = 0; (j < DIM) && !bB; j++)
554 bB = bB || (box[i][j] != 0);
558 fprintf(stderr, "coordinates %s\n", bX ? "found" : "absent");
559 fprintf(stderr, "box %s\n", bB ? "found" : "absent");
560 fprintf(stderr, "velocities %s\n", bV ? "found" : "absent");
561 fprintf(stderr, "\n");
563 /* check velocities */
567 for (i = 0; (i < natom); i++)
569 for (j = 0; (j < DIM); j++)
571 ekin += 0.5*atoms->atom[i].m*v[i][j]*v[i][j];
574 temp1 = (2.0*ekin)/(natom*DIM*BOLTZ);
575 temp2 = (2.0*ekin)/(natom*(DIM-1)*BOLTZ);
576 fprintf(stderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
577 fprintf(stderr, "Assuming the number of degrees of freedom to be "
578 "Natoms * %d or Natoms * %d,\n"
579 "the velocities correspond to a temperature of the system\n"
580 "of %g K or %g K respectively.\n\n", DIM, DIM-1, temp1, temp2);
583 /* check coordinates */
586 vdwfac2 = gmx::square(vdw_fac);
587 bonlo2 = gmx::square(bon_lo);
588 bonhi2 = gmx::square(bon_hi);
591 "Checking for atoms closer than %g and not between %g and %g,\n"
592 "relative to sum of Van der Waals distance:\n",
593 vdw_fac, bon_lo, bon_hi);
594 snew(atom_vdw, natom);
595 aps = gmx_atomprop_init();
596 for (i = 0; (i < natom); i++)
598 gmx_atomprop_query(aps, epropVDW,
599 *(atoms->resinfo[atoms->atom[i].resind].name),
600 *(atoms->atomname[i]), &(atom_vdw[i]));
603 fprintf(debug, "%5d %4s %4s %7g\n", i+1,
604 *(atoms->resinfo[atoms->atom[i].resind].name),
605 *(atoms->atomname[i]),
609 gmx_atomprop_destroy(aps);
612 set_pbc(&pbc, ePBC, box);
616 for (i = 0; (i < natom); i++)
620 fprintf(stderr, "\r%5d", i+1);
623 for (j = i+1; (j < natom); j++)
627 pbc_dx(&pbc, x[i], x[j], dx);
631 rvec_sub(x[i], x[j], dx);
634 dist2 = gmx::square(atom_vdw[i]+atom_vdw[j]);
635 if ( (r2 <= dist2*bonlo2) ||
636 ( (r2 >= dist2*bonhi2) && (r2 <= dist2*vdwfac2) ) )
640 fprintf(stderr, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
641 "atom#", "name", "residue", "r_vdw",
642 "atom#", "name", "residue", "r_vdw", "distance");
646 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
647 i+1, *(atoms->atomname[i]),
648 *(atoms->resinfo[atoms->atom[i].resind].name),
649 atoms->resinfo[atoms->atom[i].resind].nr,
651 j+1, *(atoms->atomname[j]),
652 *(atoms->resinfo[atoms->atom[j].resind].name),
653 atoms->resinfo[atoms->atom[j].resind].nr,
661 fprintf(stderr, "\rno close atoms found\n");
663 fprintf(stderr, "\r \n");
670 for (i = 0; (i < natom) && (k < 10); i++)
673 for (j = 0; (j < DIM) && !bOut; j++)
675 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
682 fprintf(stderr, "Atoms outside box ( ");
683 for (j = 0; (j < DIM); j++)
685 fprintf(stderr, "%g ", box[j][j]);
687 fprintf(stderr, "):\n"
688 "(These may occur often and are normally not a problem)\n"
689 "%5s %4s %8s %5s %s\n",
690 "atom#", "name", "residue", "r_vdw", "coordinate");
694 "%5d %4s %4s%4d %-5.3g",
695 i, *(atoms->atomname[i]),
696 *(atoms->resinfo[atoms->atom[i].resind].name),
697 atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
698 for (j = 0; (j < DIM); j++)
700 fprintf(stderr, " %6.3g", x[i][j]);
702 fprintf(stderr, "\n");
707 fprintf(stderr, "(maybe more)\n");
711 fprintf(stderr, "no atoms found outside box\n");
713 fprintf(stderr, "\n");
718 static void chk_ndx(const char *fn)
724 grps = init_index(fn, &grpname);
727 pr_blocka(debug, 0, fn, grps, FALSE);
731 printf("Contents of index file %s\n", fn);
732 printf("--------------------------------------------------\n");
733 printf("Nr. Group #Entries First Last\n");
734 for (i = 0; (i < grps->nr); i++)
736 printf("%4d %-20s%8d%8d%8d\n", i, grpname[i],
737 grps->index[i+1]-grps->index[i],
738 grps->a[grps->index[i]]+1,
739 grps->a[grps->index[i+1]-1]+1);
742 for (i = 0; (i < grps->nr); i++)
750 static void chk_enx(const char *fn)
754 gmx_enxnm_t *enm = nullptr;
758 real t0, old_t1, old_t2;
761 fprintf(stderr, "Checking energy file %s\n\n", fn);
763 in = open_enx(fn, "r");
764 do_enxnms(in, &nre, &enm);
765 fprintf(stderr, "%d groups in energy file", nre);
774 while (do_enx(in, fr))
778 if (fabs((fr->t-old_t1)-(old_t1-old_t2)) >
779 0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) )
782 fprintf(stderr, "\nTimesteps at t=%g don't match (%g, %g)\n",
783 old_t1, old_t1-old_t2, fr->t-old_t1);
795 fprintf(stderr, "\rframe: %6s (index %6d), t: %10.3f\n",
796 gmx_step_str(fr->step, buf), fnr, fr->t);
800 fprintf(stderr, "\n\nFound %d frames", fnr);
801 if (bShowTStep && fnr > 1)
803 fprintf(stderr, " with a timestep of %g ps", (old_t1-t0)/(fnr-1));
805 fprintf(stderr, ".\n");
808 free_enxnms(nre, enm);
812 int gmx_check(int argc, char *argv[])
814 const char *desc[] = {
815 "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
816 "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
817 "or an index file ([REF].ndx[ref])",
818 "and prints out useful information about them.[PAR]",
819 "Option [TT]-c[tt] checks for presence of coordinates,",
820 "velocities and box in the file, for close contacts (smaller than",
821 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
822 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
823 "radii) and atoms outside the box (these may occur often and are",
824 "no problem). If velocities are present, an estimated temperature",
825 "will be calculated from them.[PAR]",
826 "If an index file, is given its contents will be summarized.[PAR]",
827 "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
828 "the program will check whether the bond lengths defined in the tpr",
829 "file are indeed correct in the trajectory. If not you may have",
830 "non-matching files due to e.g. deshuffling or due to problems with",
831 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
832 "The program can compare two run input ([REF].tpr[ref])",
834 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied. When comparing",
835 "run input files this way, the default relative tolerance is reduced",
836 "to 0.000001 and the absolute tolerance set to zero to find any differences",
837 "not due to minor compiler optimization differences, although you can",
838 "of course still set any other tolerances through the options."
839 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
840 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
841 "For free energy simulations the A and B state topology from one",
842 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
843 "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
844 "consisting of a rough outline for a methods section for a paper."
847 { efTRX, "-f", nullptr, ffOPTRD },
848 { efTRX, "-f2", nullptr, ffOPTRD },
849 { efTPR, "-s1", "top1", ffOPTRD },
850 { efTPR, "-s2", "top2", ffOPTRD },
851 { efTPS, "-c", nullptr, ffOPTRD },
852 { efEDR, "-e", nullptr, ffOPTRD },
853 { efEDR, "-e2", "ener2", ffOPTRD },
854 { efNDX, "-n", nullptr, ffOPTRD },
855 { efTEX, "-m", nullptr, ffOPTWR }
857 #define NFILE asize(fnm)
858 const char *fn1 = nullptr, *fn2 = nullptr, *tex = nullptr;
860 gmx_output_env_t *oenv;
861 static real vdw_fac = 0.8;
862 static real bon_lo = 0.4;
863 static real bon_hi = 0.7;
864 static gmx_bool bRMSD = FALSE;
865 static real ftol = 0.001;
866 static real abstol = 0.001;
867 static gmx_bool bCompAB = FALSE;
868 static char *lastener = nullptr;
869 static t_pargs pa[] = {
870 { "-vdwfac", FALSE, etREAL, {&vdw_fac},
871 "Fraction of sum of VdW radii used as warning cutoff" },
872 { "-bonlo", FALSE, etREAL, {&bon_lo},
873 "Min. fract. of sum of VdW radii for bonded atoms" },
874 { "-bonhi", FALSE, etREAL, {&bon_hi},
875 "Max. fract. of sum of VdW radii for bonded atoms" },
876 { "-rmsd", FALSE, etBOOL, {&bRMSD},
877 "Print RMSD for x, v and f" },
878 { "-tol", FALSE, etREAL, {&ftol},
879 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
880 { "-abstol", FALSE, etREAL, {&abstol},
881 "Absolute tolerance, useful when sums are close to zero." },
882 { "-ab", FALSE, etBOOL, {&bCompAB},
883 "Compare the A and B topology from one file" },
884 { "-lastener", FALSE, etSTR, {&lastener},
885 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
888 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
889 asize(desc), desc, 0, nullptr, &oenv))
894 fn1 = opt2fn_null("-f", NFILE, fnm);
895 fn2 = opt2fn_null("-f2", NFILE, fnm);
896 tex = opt2fn_null("-m", NFILE, fnm);
899 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
903 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
907 fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
910 fn1 = opt2fn_null("-s1", NFILE, fnm);
911 fn2 = opt2fn_null("-s2", NFILE, fnm);
912 if ((fn1 && fn2) || bCompAB)
918 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
923 fprintf(stderr, "Note: When comparing run input files, default tolerances are reduced.\n");
924 if (!opt2parg_bSet("-tol", asize(pa), pa))
928 if (!opt2parg_bSet("-abstol", asize(pa), pa))
932 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
936 tpx2methods(fn1, tex);
938 else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
940 fprintf(stderr, "Please give me TWO run input (.tpr) files\n"
941 "or specify the -m flag to generate a methods.tex file\n");
944 fn1 = opt2fn_null("-e", NFILE, fnm);
945 fn2 = opt2fn_null("-e2", NFILE, fnm);
948 comp_enx(fn1, fn2, ftol, abstol, lastener);
952 chk_enx(ftp2fn(efEDR, NFILE, fnm));
956 fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
959 if (ftp2bSet(efTPS, NFILE, fnm))
961 chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
964 if (ftp2bSet(efNDX, NFILE, fnm))
966 chk_ndx(ftp2fn(efNDX, NFILE, fnm));