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, 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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/fileio/enxio.h"
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/tpxio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xtcio.h"
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdrunutility/mdmodules.h"
54 #include "gromacs/mdtypes/inputrec.h"
55 #include "gromacs/mdtypes/md_enums.h"
56 #include "gromacs/mdtypes/state.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/topology/atomprop.h"
59 #include "gromacs/topology/block.h"
60 #include "gromacs/topology/ifunc.h"
61 #include "gromacs/topology/index.h"
62 #include "gromacs/topology/mtop_util.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/trajectory/trajectoryframe.h"
65 #include "gromacs/utility/arraysize.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/futil.h"
68 #include "gromacs/utility/smalloc.h"
90 static void comp_tpx(const char *fn1
, const char *fn2
,
91 gmx_bool bRMSD
, real ftol
, real abstol
)
102 for (i
= 0; i
< (fn2
? 2 : 1); i
++)
104 ir
[i
] = new t_inputrec();
105 read_tpx_state(ff
[i
], ir
[i
], &state
[i
], &(mtop
[i
]));
106 gmx::MDModules().adjustInputrecBasedOnModules(ir
[i
]);
110 cmp_inputrec(stdout
, ir
[0], ir
[1], ftol
, abstol
);
111 /* Convert gmx_mtop_t to t_topology.
112 * We should implement direct mtop comparison,
113 * but it might be useful to keep t_topology comparison as an option.
115 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0], false);
116 top
[1] = gmx_mtop_t_to_t_topology(&mtop
[1], false);
117 cmp_top(stdout
, &top
[0], &top
[1], ftol
, abstol
);
118 cmp_groups(stdout
, &mtop
[0].groups
, &mtop
[1].groups
,
119 mtop
[0].natoms
, mtop
[1].natoms
);
120 comp_state(&state
[0], &state
[1], bRMSD
, ftol
, abstol
);
124 if (ir
[0]->efep
== efepNO
)
126 fprintf(stdout
, "inputrec->efep = %s\n", efep_names
[ir
[0]->efep
]);
132 comp_pull_AB(stdout
, ir
[0]->pull
, ftol
, abstol
);
134 /* Convert gmx_mtop_t to t_topology.
135 * We should implement direct mtop comparison,
136 * but it might be useful to keep t_topology comparison as an option.
138 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0], true);
139 cmp_top(stdout
, &top
[0], nullptr, ftol
, abstol
);
144 static void comp_trx(const gmx_output_env_t
*oenv
, const char *fn1
, const char *fn2
,
145 gmx_bool bRMSD
, real ftol
, real abstol
)
150 t_trxstatus
*status
[2];
155 fprintf(stderr
, "Comparing trajectory files %s and %s\n", fn1
, fn2
);
156 for (i
= 0; i
< 2; i
++)
158 b
[i
] = read_first_frame(oenv
, &status
[i
], fn
[i
], &fr
[i
], TRX_READ_X
|TRX_READ_V
|TRX_READ_F
);
165 comp_frame(stdout
, &(fr
[0]), &(fr
[1]), bRMSD
, ftol
, abstol
);
167 for (i
= 0; i
< 2; i
++)
169 b
[i
] = read_next_frame(oenv
, status
[i
], &fr
[i
]);
172 while (b
[0] && b
[1]);
174 for (i
= 0; i
< 2; i
++)
178 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn
[1-i
], fn
[i
]);
180 close_trx(status
[i
]);
185 fprintf(stdout
, "\nBoth files read correctly\n");
189 static void tpx2system(FILE *fp
, const gmx_mtop_t
*mtop
)
191 int nmol
, nvsite
= 0;
192 gmx_mtop_atomloop_block_t aloop
;
195 fprintf(fp
, "\\subsection{Simulation system}\n");
196 aloop
= gmx_mtop_atomloop_block_init(mtop
);
197 while (gmx_mtop_atomloop_block_next(aloop
, &atom
, &nmol
))
199 if (atom
->ptype
== eptVSite
)
204 fprintf(fp
, "A system of %d molecules (%d atoms) was simulated.\n",
205 mtop
->mols
.nr
, mtop
->natoms
-nvsite
);
208 fprintf(fp
, "Virtual sites were used in some of the molecules.\n");
213 static void tpx2params(FILE *fp
, const t_inputrec
*ir
)
215 fprintf(fp
, "\\subsection{Simulation settings}\n");
216 fprintf(fp
, "A total of %g ns were simulated with a time step of %g fs.\n",
217 ir
->nsteps
*ir
->delta_t
*0.001, 1000*ir
->delta_t
);
218 fprintf(fp
, "Neighbor searching was performed every %d steps.\n", ir
->nstlist
);
219 fprintf(fp
, "The %s algorithm was used for electrostatic interactions.\n",
220 EELTYPE(ir
->coulombtype
));
221 fprintf(fp
, "with a cut-off of %g nm.\n", ir
->rcoulomb
);
222 if (ir
->coulombtype
== eelPME
)
224 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
);
226 if (ir
->rvdw
> ir
->rlist
)
228 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
);
232 fprintf(fp
, "A single cut-off of %g was used for Van der Waals interactions.\n", ir
->rlist
);
236 fprintf(fp
, "Temperature coupling was done with the %s algorithm.\n",
237 etcoupl_names
[ir
->etc
]);
241 fprintf(fp
, "Pressure coupling was done with the %s algorithm.\n",
242 epcoupl_names
[ir
->epc
]);
247 static void tpx2methods(const char *tpx
, const char *tex
)
254 read_tpx_state(tpx
, &ir
, &state
, &mtop
);
255 fp
= gmx_fio_fopen(tex
, "w");
256 fprintf(fp
, "\\section{Methods}\n");
257 tpx2system(fp
, &mtop
);
262 static void chk_coords(int frame
, int natoms
, rvec
*x
, matrix box
, real fac
, real tol
)
268 for (i
= 0; (i
< natoms
); i
++)
270 for (j
= 0; (j
< DIM
); j
++)
272 if ((vol
> 0) && (fabs(x
[i
][j
]) > fac
*box
[j
][j
]))
274 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
278 if ((fabs(x
[i
][XX
]) < tol
) &&
279 (fabs(x
[i
][YY
]) < tol
) &&
280 (fabs(x
[i
][ZZ
]) < tol
))
287 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
292 static void chk_vels(int frame
, int natoms
, rvec
*v
)
296 for (i
= 0; (i
< natoms
); i
++)
298 for (j
= 0; (j
< DIM
); j
++)
300 if (fabs(v
[i
][j
]) > 500)
302 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
309 static void chk_forces(int frame
, int natoms
, rvec
*f
)
313 for (i
= 0; (i
< natoms
); i
++)
315 for (j
= 0; (j
< DIM
); j
++)
317 if (fabs(f
[i
][j
]) > 10000)
319 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
326 static void chk_bonds(t_idef
*idef
, int ePBC
, rvec
*x
, matrix box
, real tol
)
328 int ftype
, k
, ai
, aj
, type
;
329 real b0
, blen
, deviation
;
333 set_pbc(&pbc
, ePBC
, box
);
334 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
336 if ((interaction_function
[ftype
].flags
& IF_CHEMBOND
) == IF_CHEMBOND
)
338 for (k
= 0; (k
< idef
->il
[ftype
].nr
); )
340 type
= idef
->il
[ftype
].iatoms
[k
++];
341 ai
= idef
->il
[ftype
].iatoms
[k
++];
342 aj
= idef
->il
[ftype
].iatoms
[k
++];
347 b0
= idef
->iparams
[type
].harmonic
.rA
;
350 b0
= std::sqrt(idef
->iparams
[type
].harmonic
.rA
);
353 b0
= idef
->iparams
[type
].morse
.b0A
;
356 b0
= idef
->iparams
[type
].cubic
.b0
;
359 b0
= idef
->iparams
[type
].constr
.dA
;
366 pbc_dx(&pbc
, x
[ai
], x
[aj
], dx
);
368 deviation
= gmx::square(blen
-b0
);
369 if (std::sqrt(deviation
/gmx::square(b0
)) > tol
)
371 fprintf(stderr
, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai
+1, aj
+1, blen
, b0
);
379 void chk_trj(const gmx_output_env_t
*oenv
, const char *fn
, const char *tpr
, real tol
)
383 t_fr_time first
, last
;
384 int j
= -1, new_natoms
, natoms
;
386 gmx_bool bShowTimestep
= TRUE
, newline
= FALSE
;
389 gmx_localtop_t
*top
= nullptr;
395 read_tpx_state(tpr
, &ir
, &state
, &mtop
);
396 top
= gmx_mtop_generate_local_top(&mtop
, ir
.efep
!= efepNO
);
401 printf("Checking file %s\n", fn
);
431 read_first_frame(oenv
, &status
, fn
, &fr
, TRX_READ_X
| TRX_READ_V
| TRX_READ_F
);
437 fprintf(stderr
, "\n# Atoms %d\n", fr
.natoms
);
440 fprintf(stderr
, "Precision %g (nm)\n", 1/fr
.prec
);
444 if ((natoms
> 0) && (new_natoms
!= natoms
))
446 fprintf(stderr
, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
447 old_t1
, natoms
, new_natoms
);
452 if (fabs((fr
.time
-old_t1
)-(old_t1
-old_t2
)) >
453 0.1*(fabs(fr
.time
-old_t1
)+fabs(old_t1
-old_t2
)) )
455 bShowTimestep
= FALSE
;
456 fprintf(stderr
, "%sTimesteps at t=%g don't match (%g, %g)\n",
457 newline
? "\n" : "", old_t1
, old_t1
-old_t2
, fr
.time
-old_t1
);
463 chk_bonds(&top
->idef
, ir
.ePBC
, fr
.x
, fr
.box
, tol
);
467 chk_coords(j
, natoms
, fr
.x
, fr
.box
, 1e5
, tol
);
471 chk_vels(j
, natoms
, fr
.v
);
475 chk_forces(j
, natoms
, fr
.f
);
481 new_natoms
= fr
.natoms
;
482 #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++; \
484 INC(fr
, count
, first
, last
, bStep
);
485 INC(fr
, count
, first
, last
, bTime
);
486 INC(fr
, count
, first
, last
, bLambda
);
487 INC(fr
, count
, first
, last
, bX
);
488 INC(fr
, count
, first
, last
, bV
);
489 INC(fr
, count
, first
, last
, bF
);
490 INC(fr
, count
, first
, last
, bBox
);
493 while (read_next_frame(oenv
, status
, &fr
));
495 fprintf(stderr
, "\n");
499 fprintf(stderr
, "\nItem #frames");
502 fprintf(stderr
, " Timestep (ps)");
504 fprintf(stderr
, "\n");
505 #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")
506 PRINTITEM ( "Step", bStep
);
507 PRINTITEM ( "Time", bTime
);
508 PRINTITEM ( "Lambda", bLambda
);
509 PRINTITEM ( "Coords", bX
);
510 PRINTITEM ( "Velocities", bV
);
511 PRINTITEM ( "Forces", bF
);
512 PRINTITEM ( "Box", bBox
);
515 void chk_tps(const char *fn
, real vdw_fac
, real bon_lo
, real bon_hi
)
525 gmx_bool bV
, bX
, bB
, bFirst
, bOut
;
526 real r2
, ekin
, temp1
, temp2
, dist2
, vdwfac2
, bonlo2
, bonhi2
;
530 fprintf(stderr
, "Checking coordinate file %s\n", fn
);
531 read_tps_conf(fn
, &top
, &ePBC
, &x
, &v
, box
, TRUE
);
534 fprintf(stderr
, "%d atoms in file\n", atoms
->nr
);
536 /* check coordinates and box */
539 for (i
= 0; (i
< natom
) && !(bV
&& bX
); i
++)
541 for (j
= 0; (j
< DIM
) && !(bV
&& bX
); j
++)
543 bV
= bV
|| (v
[i
][j
] != 0);
544 bX
= bX
|| (x
[i
][j
] != 0);
548 for (i
= 0; (i
< DIM
) && !bB
; i
++)
550 for (j
= 0; (j
< DIM
) && !bB
; j
++)
552 bB
= bB
|| (box
[i
][j
] != 0);
556 fprintf(stderr
, "coordinates %s\n", bX
? "found" : "absent");
557 fprintf(stderr
, "box %s\n", bB
? "found" : "absent");
558 fprintf(stderr
, "velocities %s\n", bV
? "found" : "absent");
559 fprintf(stderr
, "\n");
561 /* check velocities */
565 for (i
= 0; (i
< natom
); i
++)
567 for (j
= 0; (j
< DIM
); j
++)
569 ekin
+= 0.5*atoms
->atom
[i
].m
*v
[i
][j
]*v
[i
][j
];
572 temp1
= (2.0*ekin
)/(natom
*DIM
*BOLTZ
);
573 temp2
= (2.0*ekin
)/(natom
*(DIM
-1)*BOLTZ
);
574 fprintf(stderr
, "Kinetic energy: %g (kJ/mol)\n", ekin
);
575 fprintf(stderr
, "Assuming the number of degrees of freedom to be "
576 "Natoms * %d or Natoms * %d,\n"
577 "the velocities correspond to a temperature of the system\n"
578 "of %g K or %g K respectively.\n\n", DIM
, DIM
-1, temp1
, temp2
);
581 /* check coordinates */
584 vdwfac2
= gmx::square(vdw_fac
);
585 bonlo2
= gmx::square(bon_lo
);
586 bonhi2
= gmx::square(bon_hi
);
589 "Checking for atoms closer than %g and not between %g and %g,\n"
590 "relative to sum of Van der Waals distance:\n",
591 vdw_fac
, bon_lo
, bon_hi
);
592 snew(atom_vdw
, natom
);
593 aps
= gmx_atomprop_init();
594 for (i
= 0; (i
< natom
); i
++)
596 gmx_atomprop_query(aps
, epropVDW
,
597 *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
),
598 *(atoms
->atomname
[i
]), &(atom_vdw
[i
]));
601 fprintf(debug
, "%5d %4s %4s %7g\n", i
+1,
602 *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
),
603 *(atoms
->atomname
[i
]),
607 gmx_atomprop_destroy(aps
);
610 set_pbc(&pbc
, ePBC
, box
);
614 for (i
= 0; (i
< natom
); i
++)
618 fprintf(stderr
, "\r%5d", i
+1);
621 for (j
= i
+1; (j
< natom
); j
++)
625 pbc_dx(&pbc
, x
[i
], x
[j
], dx
);
629 rvec_sub(x
[i
], x
[j
], dx
);
632 dist2
= gmx::square(atom_vdw
[i
]+atom_vdw
[j
]);
633 if ( (r2
<= dist2
*bonlo2
) ||
634 ( (r2
>= dist2
*bonhi2
) && (r2
<= dist2
*vdwfac2
) ) )
638 fprintf(stderr
, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
639 "atom#", "name", "residue", "r_vdw",
640 "atom#", "name", "residue", "r_vdw", "distance");
644 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
645 i
+1, *(atoms
->atomname
[i
]),
646 *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
),
647 atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
,
649 j
+1, *(atoms
->atomname
[j
]),
650 *(atoms
->resinfo
[atoms
->atom
[j
].resind
].name
),
651 atoms
->resinfo
[atoms
->atom
[j
].resind
].nr
,
659 fprintf(stderr
, "\rno close atoms found\n");
661 fprintf(stderr
, "\r \n");
668 for (i
= 0; (i
< natom
) && (k
< 10); i
++)
671 for (j
= 0; (j
< DIM
) && !bOut
; j
++)
673 bOut
= bOut
|| (x
[i
][j
] < 0) || (x
[i
][j
] > box
[j
][j
]);
680 fprintf(stderr
, "Atoms outside box ( ");
681 for (j
= 0; (j
< DIM
); j
++)
683 fprintf(stderr
, "%g ", box
[j
][j
]);
685 fprintf(stderr
, "):\n"
686 "(These may occur often and are normally not a problem)\n"
687 "%5s %4s %8s %5s %s\n",
688 "atom#", "name", "residue", "r_vdw", "coordinate");
692 "%5d %4s %4s%4d %-5.3g",
693 i
, *(atoms
->atomname
[i
]),
694 *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
),
695 atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
, atom_vdw
[i
]);
696 for (j
= 0; (j
< DIM
); j
++)
698 fprintf(stderr
, " %6.3g", x
[i
][j
]);
700 fprintf(stderr
, "\n");
705 fprintf(stderr
, "(maybe more)\n");
709 fprintf(stderr
, "no atoms found outside box\n");
711 fprintf(stderr
, "\n");
716 void chk_ndx(const char *fn
)
722 grps
= init_index(fn
, &grpname
);
725 pr_blocka(debug
, 0, fn
, grps
, FALSE
);
729 printf("Contents of index file %s\n", fn
);
730 printf("--------------------------------------------------\n");
731 printf("Nr. Group #Entries First Last\n");
732 for (i
= 0; (i
< grps
->nr
); i
++)
734 printf("%4d %-20s%8d%8d%8d\n", i
, grpname
[i
],
735 grps
->index
[i
+1]-grps
->index
[i
],
736 grps
->a
[grps
->index
[i
]]+1,
737 grps
->a
[grps
->index
[i
+1]-1]+1);
740 for (i
= 0; (i
< grps
->nr
); i
++)
748 void chk_enx(const char *fn
)
752 gmx_enxnm_t
*enm
= nullptr;
756 real t0
, old_t1
, old_t2
;
759 fprintf(stderr
, "Checking energy file %s\n\n", fn
);
761 in
= open_enx(fn
, "r");
762 do_enxnms(in
, &nre
, &enm
);
763 fprintf(stderr
, "%d groups in energy file", nre
);
772 while (do_enx(in
, fr
))
776 if (fabs((fr
->t
-old_t1
)-(old_t1
-old_t2
)) >
777 0.1*(fabs(fr
->t
-old_t1
)+fabs(old_t1
-old_t2
)) )
780 fprintf(stderr
, "\nTimesteps at t=%g don't match (%g, %g)\n",
781 old_t1
, old_t1
-old_t2
, fr
->t
-old_t1
);
793 fprintf(stderr
, "\rframe: %6s (index %6d), t: %10.3f\n",
794 gmx_step_str(fr
->step
, buf
), fnr
, fr
->t
);
798 fprintf(stderr
, "\n\nFound %d frames", fnr
);
799 if (bShowTStep
&& fnr
> 1)
801 fprintf(stderr
, " with a timestep of %g ps", (old_t1
-t0
)/(fnr
-1));
803 fprintf(stderr
, ".\n");
806 free_enxnms(nre
, enm
);
810 int gmx_check(int argc
, char *argv
[])
812 const char *desc
[] = {
813 "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
814 "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
815 "or an index file ([REF].ndx[ref])",
816 "and prints out useful information about them.[PAR]",
817 "Option [TT]-c[tt] checks for presence of coordinates,",
818 "velocities and box in the file, for close contacts (smaller than",
819 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
820 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
821 "radii) and atoms outside the box (these may occur often and are",
822 "no problem). If velocities are present, an estimated temperature",
823 "will be calculated from them.[PAR]",
824 "If an index file, is given its contents will be summarized.[PAR]",
825 "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
826 "the program will check whether the bond lengths defined in the tpr",
827 "file are indeed correct in the trajectory. If not you may have",
828 "non-matching files due to e.g. deshuffling or due to problems with",
829 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
830 "The program can compare two run input ([REF].tpr[ref])",
832 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
833 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
834 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
835 "For free energy simulations the A and B state topology from one",
836 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
837 "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
838 "consisting of a rough outline for a methods section for a paper."
841 { efTRX
, "-f", nullptr, ffOPTRD
},
842 { efTRX
, "-f2", nullptr, ffOPTRD
},
843 { efTPR
, "-s1", "top1", ffOPTRD
},
844 { efTPR
, "-s2", "top2", ffOPTRD
},
845 { efTPS
, "-c", nullptr, ffOPTRD
},
846 { efEDR
, "-e", nullptr, ffOPTRD
},
847 { efEDR
, "-e2", "ener2", ffOPTRD
},
848 { efNDX
, "-n", nullptr, ffOPTRD
},
849 { efTEX
, "-m", nullptr, ffOPTWR
}
851 #define NFILE asize(fnm)
852 const char *fn1
= nullptr, *fn2
= nullptr, *tex
= nullptr;
854 gmx_output_env_t
*oenv
;
855 static real vdw_fac
= 0.8;
856 static real bon_lo
= 0.4;
857 static real bon_hi
= 0.7;
858 static gmx_bool bRMSD
= FALSE
;
859 static real ftol
= 0.001;
860 static real abstol
= 0.001;
861 static gmx_bool bCompAB
= FALSE
;
862 static char *lastener
= nullptr;
863 static t_pargs pa
[] = {
864 { "-vdwfac", FALSE
, etREAL
, {&vdw_fac
},
865 "Fraction of sum of VdW radii used as warning cutoff" },
866 { "-bonlo", FALSE
, etREAL
, {&bon_lo
},
867 "Min. fract. of sum of VdW radii for bonded atoms" },
868 { "-bonhi", FALSE
, etREAL
, {&bon_hi
},
869 "Max. fract. of sum of VdW radii for bonded atoms" },
870 { "-rmsd", FALSE
, etBOOL
, {&bRMSD
},
871 "Print RMSD for x, v and f" },
872 { "-tol", FALSE
, etREAL
, {&ftol
},
873 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
874 { "-abstol", FALSE
, etREAL
, {&abstol
},
875 "Absolute tolerance, useful when sums are close to zero." },
876 { "-ab", FALSE
, etBOOL
, {&bCompAB
},
877 "Compare the A and B topology from one file" },
878 { "-lastener", FALSE
, etSTR
, {&lastener
},
879 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
882 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
883 asize(desc
), desc
, 0, nullptr, &oenv
))
888 fn1
= opt2fn_null("-f", NFILE
, fnm
);
889 fn2
= opt2fn_null("-f2", NFILE
, fnm
);
890 tex
= opt2fn_null("-m", NFILE
, fnm
);
893 comp_trx(oenv
, fn1
, fn2
, bRMSD
, ftol
, abstol
);
897 chk_trj(oenv
, fn1
, opt2fn_null("-s1", NFILE
, fnm
), ftol
);
901 fprintf(stderr
, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
904 fn1
= opt2fn_null("-s1", NFILE
, fnm
);
905 fn2
= opt2fn_null("-s2", NFILE
, fnm
);
906 if ((fn1
&& fn2
) || bCompAB
)
912 gmx_fatal(FARGS
, "With -ab you need to set the -s1 option");
916 comp_tpx(fn1
, fn2
, bRMSD
, ftol
, abstol
);
920 tpx2methods(fn1
, tex
);
922 else if ((fn1
&& !opt2fn_null("-f", NFILE
, fnm
)) || (!fn1
&& fn2
))
924 fprintf(stderr
, "Please give me TWO run input (.tpr) files\n"
925 "or specify the -m flag to generate a methods.tex file\n");
928 fn1
= opt2fn_null("-e", NFILE
, fnm
);
929 fn2
= opt2fn_null("-e2", NFILE
, fnm
);
932 comp_enx(fn1
, fn2
, ftol
, abstol
, lastener
);
936 chk_enx(ftp2fn(efEDR
, NFILE
, fnm
));
940 fprintf(stderr
, "Please give me TWO energy (.edr/.ene) files!\n");
943 if (ftp2bSet(efTPS
, NFILE
, fnm
))
945 chk_tps(ftp2fn(efTPS
, NFILE
, fnm
), vdw_fac
, bon_lo
, bon_hi
);
948 if (ftp2bSet(efNDX
, NFILE
, fnm
))
950 chk_ndx(ftp2fn(efNDX
, NFILE
, fnm
));