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-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, 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.
37 /* This file is completely threadsafe - keep it that way! */
47 #include "gromacs/fileio/enxio.h"
48 #include "gromacs/fileio/tpxio.h"
49 #include "gromacs/fileio/trx.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/fileio/txtdump.h"
52 #include "gromacs/gmxlib/ifunc.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/mdtypes/pull-params.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
61 #include "gromacs/utility/stringutil.h"
63 static void cmp_int(FILE *fp
, const char *s
, int index
, int i1
, int i2
)
69 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
73 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
78 static void cmp_int64(FILE *fp
, const char *s
, gmx_int64_t i1
, gmx_int64_t i2
)
82 fprintf(fp
, "%s (", s
);
83 fprintf(fp
, "%" GMX_PRId64
, i1
);
85 fprintf(fp
, "%" GMX_PRId64
, i2
);
90 static void cmp_us(FILE *fp
, const char *s
, int index
, unsigned short i1
, unsigned short i2
)
96 fprintf(fp
, "%s[%d] (%hu - %hu)\n", s
, index
, i1
, i2
);
100 fprintf(fp
, "%s (%hu - %hu)\n", s
, i1
, i2
);
105 static void cmp_uc(FILE *fp
, const char *s
, int index
, unsigned char i1
, unsigned char i2
)
111 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
115 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
120 static gmx_bool
cmp_bool(FILE *fp
, const char *s
, int index
, gmx_bool b1
, gmx_bool b2
)
142 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
,
143 gmx::boolToString(b1
), gmx::boolToString(b2
));
147 fprintf(fp
, "%s (%s - %s)\n", s
,
148 gmx::boolToString(b1
), gmx::boolToString(b2
));
154 static void cmp_str(FILE *fp
, const char *s
, int index
,
155 const char *s1
, const char *s2
)
157 if (std::strcmp(s1
, s2
) != 0)
161 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
, s1
, s2
);
165 fprintf(fp
, "%s (%s - %s)\n", s
, s1
, s2
);
170 static gmx_bool
equal_real(real i1
, real i2
, real ftol
, real abstol
)
172 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
175 static gmx_bool
equal_float(float i1
, float i2
, float ftol
, float abstol
)
177 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
180 static gmx_bool
equal_double(double i1
, double i2
, real ftol
, real abstol
)
182 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
186 cmp_real(FILE *fp
, const char *s
, int index
, real i1
, real i2
, real ftol
, real abstol
)
188 if (!equal_real(i1
, i2
, ftol
, abstol
))
192 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
196 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
202 cmp_float(FILE *fp
, const char *s
, int index
, float i1
, float i2
, float ftol
, float abstol
)
204 if (!equal_float(i1
, i2
, ftol
, abstol
))
208 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
212 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
220 cmp_double(FILE *fp
, const char *s
, int index
, double i1
, double i2
, double ftol
, double abstol
)
222 if (!equal_double(i1
, i2
, ftol
, abstol
))
226 fprintf(fp
, "%s[%2d] (%16.9e - %16.9e)\n", s
, index
, i1
, i2
);
230 fprintf(fp
, "%s (%16.9e - %16.9e)\n", s
, i1
, i2
);
235 static void cmp_rvec(FILE *fp
, const char *s
, int index
, rvec i1
, rvec i2
, real ftol
, real abstol
)
237 if (!equal_real(i1
[XX
], i2
[XX
], ftol
, abstol
) ||
238 !equal_real(i1
[YY
], i2
[YY
], ftol
, abstol
) ||
239 !equal_real(i1
[ZZ
], i2
[ZZ
], ftol
, abstol
))
243 fprintf(fp
, "%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
244 s
, index
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
248 fprintf(fp
, "%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
249 s
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
254 static void cmp_ivec(FILE *fp
, const char *s
, int index
, ivec i1
, ivec i2
)
256 if ((i1
[XX
] != i2
[XX
]) || (i1
[YY
] != i2
[YY
]) || (i1
[ZZ
] != i2
[ZZ
]))
260 fprintf(fp
, "%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
, index
,
261 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
265 fprintf(fp
, "%s (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
,
266 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
271 static void cmp_ilist(FILE *fp
, int ftype
, t_ilist
*il1
, t_ilist
*il2
)
276 fprintf(fp
, "comparing ilist %s\n", interaction_function
[ftype
].name
);
277 sprintf(buf
, "%s->nr", interaction_function
[ftype
].name
);
278 cmp_int(fp
, buf
, -1, il1
->nr
, il2
->nr
);
279 sprintf(buf
, "%s->iatoms", interaction_function
[ftype
].name
);
280 if (((il1
->nr
> 0) && (!il1
->iatoms
)) ||
281 ((il2
->nr
> 0) && (!il2
->iatoms
)) ||
282 ((il1
->nr
!= il2
->nr
)))
284 fprintf(fp
, "Comparing radically different topologies - %s is different\n",
289 for (i
= 0; (i
< il1
->nr
); i
++)
291 cmp_int(fp
, buf
, i
, il1
->iatoms
[i
], il2
->iatoms
[i
]);
296 void cmp_iparm(FILE *fp
, const char *s
, t_functype ft
,
297 t_iparams ip1
, t_iparams ip2
, real ftol
, real abstol
)
303 for (i
= 0; i
< MAXFORCEPARAM
&& !bDiff
; i
++)
305 bDiff
= !equal_real(ip1
.generic
.buf
[i
], ip2
.generic
.buf
[i
], ftol
, abstol
);
309 fprintf(fp
, "%s1: ", s
);
310 pr_iparams(fp
, ft
, &ip1
);
311 fprintf(fp
, "%s2: ", s
);
312 pr_iparams(fp
, ft
, &ip2
);
316 void cmp_iparm_AB(FILE *fp
, const char *s
, t_functype ft
, t_iparams ip1
, real ftol
, real abstol
)
318 int nrfpA
, nrfpB
, p0
, i
;
321 /* Normally the first parameter is perturbable */
323 nrfpA
= interaction_function
[ft
].nrfpA
;
324 nrfpB
= interaction_function
[ft
].nrfpB
;
329 else if (interaction_function
[ft
].flags
& IF_TABULATED
)
331 /* For tabulated interactions only the second parameter is perturbable */
336 for (i
= 0; i
< nrfpB
&& !bDiff
; i
++)
338 bDiff
= !equal_real(ip1
.generic
.buf
[p0
+i
], ip1
.generic
.buf
[nrfpA
+i
], ftol
, abstol
);
342 fprintf(fp
, "%s: ", s
);
343 pr_iparams(fp
, ft
, &ip1
);
347 static void cmp_idef(FILE *fp
, t_idef
*id1
, t_idef
*id2
, real ftol
, real abstol
)
350 char buf1
[64], buf2
[64];
352 fprintf(fp
, "comparing idef\n");
355 cmp_int(fp
, "idef->ntypes", -1, id1
->ntypes
, id2
->ntypes
);
356 cmp_int(fp
, "idef->atnr", -1, id1
->atnr
, id2
->atnr
);
357 for (i
= 0; (i
< std::min(id1
->ntypes
, id2
->ntypes
)); i
++)
359 sprintf(buf1
, "idef->functype[%d]", i
);
360 sprintf(buf2
, "idef->iparam[%d]", i
);
361 cmp_int(fp
, buf1
, i
, (int)id1
->functype
[i
], (int)id2
->functype
[i
]);
362 cmp_iparm(fp
, buf2
, id1
->functype
[i
],
363 id1
->iparams
[i
], id2
->iparams
[i
], ftol
, abstol
);
365 cmp_real(fp
, "fudgeQQ", -1, id1
->fudgeQQ
, id2
->fudgeQQ
, ftol
, abstol
);
366 for (i
= 0; (i
< F_NRE
); i
++)
368 cmp_ilist(fp
, i
, &(id1
->il
[i
]), &(id2
->il
[i
]));
373 for (i
= 0; (i
< id1
->ntypes
); i
++)
375 cmp_iparm_AB(fp
, "idef->iparam", id1
->functype
[i
], id1
->iparams
[i
], ftol
, abstol
);
380 static void cmp_block(FILE *fp
, t_block
*b1
, t_block
*b2
, const char *s
)
384 fprintf(fp
, "comparing block %s\n", s
);
385 sprintf(buf
, "%s.nr", s
);
386 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
389 static void cmp_blocka(FILE *fp
, t_blocka
*b1
, t_blocka
*b2
, const char *s
)
393 fprintf(fp
, "comparing blocka %s\n", s
);
394 sprintf(buf
, "%s.nr", s
);
395 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
396 sprintf(buf
, "%s.nra", s
);
397 cmp_int(fp
, buf
, -1, b1
->nra
, b2
->nra
);
400 static void cmp_atom(FILE *fp
, int index
, t_atom
*a1
, t_atom
*a2
, real ftol
, real abstol
)
404 cmp_us(fp
, "atom.type", index
, a1
->type
, a2
->type
);
405 cmp_us(fp
, "atom.ptype", index
, a1
->ptype
, a2
->ptype
);
406 cmp_int(fp
, "atom.resind", index
, a1
->resind
, a2
->resind
);
407 cmp_int(fp
, "atom.atomnumber", index
, a1
->atomnumber
, a2
->atomnumber
);
408 cmp_real(fp
, "atom.m", index
, a1
->m
, a2
->m
, ftol
, abstol
);
409 cmp_real(fp
, "atom.q", index
, a1
->q
, a2
->q
, ftol
, abstol
);
410 cmp_us(fp
, "atom.typeB", index
, a1
->typeB
, a2
->typeB
);
411 cmp_real(fp
, "atom.mB", index
, a1
->mB
, a2
->mB
, ftol
, abstol
);
412 cmp_real(fp
, "atom.qB", index
, a1
->qB
, a2
->qB
, ftol
, abstol
);
416 cmp_us(fp
, "atom.type", index
, a1
->type
, a1
->typeB
);
417 cmp_real(fp
, "atom.m", index
, a1
->m
, a1
->mB
, ftol
, abstol
);
418 cmp_real(fp
, "atom.q", index
, a1
->q
, a1
->qB
, ftol
, abstol
);
422 static void cmp_atoms(FILE *fp
, t_atoms
*a1
, t_atoms
*a2
, real ftol
, real abstol
)
426 fprintf(fp
, "comparing atoms\n");
430 cmp_int(fp
, "atoms->nr", -1, a1
->nr
, a2
->nr
);
431 for (i
= 0; (i
< a1
->nr
); i
++)
433 cmp_atom(fp
, i
, &(a1
->atom
[i
]), &(a2
->atom
[i
]), ftol
, abstol
);
438 for (i
= 0; (i
< a1
->nr
); i
++)
440 cmp_atom(fp
, i
, &(a1
->atom
[i
]), NULL
, ftol
, abstol
);
445 static void cmp_top(FILE *fp
, t_topology
*t1
, t_topology
*t2
, real ftol
, real abstol
)
447 fprintf(fp
, "comparing top\n");
450 cmp_idef(fp
, &(t1
->idef
), &(t2
->idef
), ftol
, abstol
);
451 cmp_atoms(fp
, &(t1
->atoms
), &(t2
->atoms
), ftol
, abstol
);
452 cmp_block(fp
, &t1
->cgs
, &t2
->cgs
, "cgs");
453 cmp_block(fp
, &t1
->mols
, &t2
->mols
, "mols");
454 cmp_bool(fp
, "bIntermolecularInteractions", -1, t1
->bIntermolecularInteractions
, t2
->bIntermolecularInteractions
);
455 cmp_blocka(fp
, &t1
->excls
, &t2
->excls
, "excls");
459 cmp_idef(fp
, &(t1
->idef
), NULL
, ftol
, abstol
);
460 cmp_atoms(fp
, &(t1
->atoms
), NULL
, ftol
, abstol
);
464 static void cmp_groups(FILE *fp
, gmx_groups_t
*g0
, gmx_groups_t
*g1
,
465 int natoms0
, int natoms1
)
470 fprintf(fp
, "comparing groups\n");
472 for (i
= 0; i
< egcNR
; i
++)
474 sprintf(buf
, "grps[%d].nr", i
);
475 cmp_int(fp
, buf
, -1, g0
->grps
[i
].nr
, g1
->grps
[i
].nr
);
476 if (g0
->grps
[i
].nr
== g1
->grps
[i
].nr
)
478 for (j
= 0; j
< g0
->grps
[i
].nr
; j
++)
480 sprintf(buf
, "grps[%d].name[%d]", i
, j
);
482 *g0
->grpname
[g0
->grps
[i
].nm_ind
[j
]],
483 *g1
->grpname
[g1
->grps
[i
].nm_ind
[j
]]);
486 cmp_int(fp
, "ngrpnr", i
, g0
->ngrpnr
[i
], g1
->ngrpnr
[i
]);
487 if (g0
->ngrpnr
[i
] == g1
->ngrpnr
[i
] && natoms0
== natoms1
&&
488 (g0
->grpnr
[i
] != NULL
|| g1
->grpnr
[i
] != NULL
))
490 for (j
= 0; j
< natoms0
; j
++)
492 cmp_int(fp
, gtypes
[i
], j
, ggrpnr(g0
, i
, j
), ggrpnr(g1
, i
, j
));
496 /* We have compared the names in the groups lists,
497 * so we can skip the grpname list comparison.
501 static void cmp_rvecs(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
502 gmx_bool bRMSD
, real ftol
, real abstol
)
510 for (i
= 0; (i
< n
); i
++)
512 for (m
= 0; m
< DIM
; m
++)
514 d
= x1
[i
][m
] - x2
[i
][m
];
518 fprintf(fp
, "%s RMSD %g\n", title
, std::sqrt(ssd
/n
));
522 for (i
= 0; (i
< n
); i
++)
524 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], ftol
, abstol
);
530 /* Similar to cmp_rvecs, but this routine scales the allowed absolute tolerance
531 * by the RMS of the force components of x1.
533 static void cmp_rvecs_rmstol(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
534 real ftol
, real abstol
)
538 double ave_x1
, rms_x1
;
540 /* It is tricky to compare real values, in particular forces that
541 * are sums of lots of terms where the final value might be close to 0.0.
542 * To get a reference magnitude we calculate the RMS value of each
543 * component in x1, and then set the allowed absolute tolerance to the
544 * relative tolerance times this RMS magnitude.
547 for (i
= 0; i
< n
; i
++)
549 for (m
= 0; m
< DIM
; m
++)
557 for (i
= 0; (i
< n
); i
++)
559 for (m
= 0; m
< DIM
; m
++)
561 d
= x1
[i
][m
] - ave_x1
;
565 rms_x1
= std::sqrt(rms_x1
/(DIM
*n
));
566 /* And now do the actual comparision with a hopefully realistic abstol. */
567 for (i
= 0; (i
< n
); i
++)
569 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], ftol
, abstol
*rms_x1
);
573 static void cmp_grpopts(FILE *fp
, t_grpopts
*opt1
, t_grpopts
*opt2
, real ftol
, real abstol
)
576 char buf1
[256], buf2
[256];
578 cmp_int(fp
, "inputrec->grpopts.ngtc", -1, opt1
->ngtc
, opt2
->ngtc
);
579 cmp_int(fp
, "inputrec->grpopts.ngacc", -1, opt1
->ngacc
, opt2
->ngacc
);
580 cmp_int(fp
, "inputrec->grpopts.ngfrz", -1, opt1
->ngfrz
, opt2
->ngfrz
);
581 cmp_int(fp
, "inputrec->grpopts.ngener", -1, opt1
->ngener
, opt2
->ngener
);
582 for (i
= 0; (i
< std::min(opt1
->ngtc
, opt2
->ngtc
)); i
++)
584 cmp_real(fp
, "inputrec->grpopts.nrdf", i
, opt1
->nrdf
[i
], opt2
->nrdf
[i
], ftol
, abstol
);
585 cmp_real(fp
, "inputrec->grpopts.ref_t", i
, opt1
->ref_t
[i
], opt2
->ref_t
[i
], ftol
, abstol
);
586 cmp_real(fp
, "inputrec->grpopts.tau_t", i
, opt1
->tau_t
[i
], opt2
->tau_t
[i
], ftol
, abstol
);
587 cmp_int(fp
, "inputrec->grpopts.annealing", i
, opt1
->annealing
[i
], opt2
->annealing
[i
]);
588 cmp_int(fp
, "inputrec->grpopts.anneal_npoints", i
,
589 opt1
->anneal_npoints
[i
], opt2
->anneal_npoints
[i
]);
590 if (opt1
->anneal_npoints
[i
] == opt2
->anneal_npoints
[i
])
592 sprintf(buf1
, "inputrec->grpopts.anneal_time[%d]", i
);
593 sprintf(buf2
, "inputrec->grpopts.anneal_temp[%d]", i
);
594 for (j
= 0; j
< opt1
->anneal_npoints
[i
]; j
++)
596 cmp_real(fp
, buf1
, j
, opt1
->anneal_time
[i
][j
], opt2
->anneal_time
[i
][j
], ftol
, abstol
);
597 cmp_real(fp
, buf2
, j
, opt1
->anneal_temp
[i
][j
], opt2
->anneal_temp
[i
][j
], ftol
, abstol
);
601 if (opt1
->ngener
== opt2
->ngener
)
603 for (i
= 0; i
< opt1
->ngener
; i
++)
605 for (j
= i
; j
< opt1
->ngener
; j
++)
607 sprintf(buf1
, "inputrec->grpopts.egp_flags[%d]", i
);
609 opt1
->egp_flags
[opt1
->ngener
*i
+j
],
610 opt2
->egp_flags
[opt1
->ngener
*i
+j
]);
614 for (i
= 0; (i
< std::min(opt1
->ngacc
, opt2
->ngacc
)); i
++)
616 cmp_rvec(fp
, "inputrec->grpopts.acc", i
, opt1
->acc
[i
], opt2
->acc
[i
], ftol
, abstol
);
618 for (i
= 0; (i
< std::min(opt1
->ngfrz
, opt2
->ngfrz
)); i
++)
620 cmp_ivec(fp
, "inputrec->grpopts.nFreeze", i
, opt1
->nFreeze
[i
], opt2
->nFreeze
[i
]);
624 static void cmp_cosines(FILE *fp
, const char *s
, t_cosines c1
[DIM
], t_cosines c2
[DIM
], real ftol
, real abstol
)
629 for (m
= 0; (m
< DIM
); m
++)
631 sprintf(buf
, "inputrec->%s[%d]", s
, m
);
632 cmp_int(fp
, buf
, 0, c1
->n
, c2
->n
);
633 for (i
= 0; (i
< std::min(c1
->n
, c2
->n
)); i
++)
635 cmp_real(fp
, buf
, i
, c1
->a
[i
], c2
->a
[i
], ftol
, abstol
);
636 cmp_real(fp
, buf
, i
, c1
->phi
[i
], c2
->phi
[i
], ftol
, abstol
);
640 static void cmp_pull(FILE *fp
)
642 fprintf(fp
, "WARNING: Both files use COM pulling, but comparing of the pull struct is not implemented (yet). The pull parameters could be the same or different.\n");
645 static void cmp_simtempvals(FILE *fp
, t_simtemp
*simtemp1
, t_simtemp
*simtemp2
, int n_lambda
, real ftol
, real abstol
)
648 cmp_int(fp
, "inputrec->simtempvals->eSimTempScale", -1, simtemp1
->eSimTempScale
, simtemp2
->eSimTempScale
);
649 cmp_real(fp
, "inputrec->simtempvals->simtemp_high", -1, simtemp1
->simtemp_high
, simtemp2
->simtemp_high
, ftol
, abstol
);
650 cmp_real(fp
, "inputrec->simtempvals->simtemp_low", -1, simtemp1
->simtemp_low
, simtemp2
->simtemp_low
, ftol
, abstol
);
651 for (i
= 0; i
< n_lambda
; i
++)
653 cmp_real(fp
, "inputrec->simtempvals->temperatures", -1, simtemp1
->temperatures
[i
], simtemp2
->temperatures
[i
], ftol
, abstol
);
657 static void cmp_expandedvals(FILE *fp
, t_expanded
*expand1
, t_expanded
*expand2
, int n_lambda
, real ftol
, real abstol
)
661 cmp_bool(fp
, "inputrec->fepvals->bInit_weights", -1, expand1
->bInit_weights
, expand2
->bInit_weights
);
662 cmp_bool(fp
, "inputrec->fepvals->bWLoneovert", -1, expand1
->bWLoneovert
, expand2
->bWLoneovert
);
664 for (i
= 0; i
< n_lambda
; i
++)
666 cmp_real(fp
, "inputrec->expandedvals->init_lambda_weights", -1,
667 expand1
->init_lambda_weights
[i
], expand2
->init_lambda_weights
[i
], ftol
, abstol
);
670 cmp_int(fp
, "inputrec->expandedvals->lambda-stats", -1, expand1
->elamstats
, expand2
->elamstats
);
671 cmp_int(fp
, "inputrec->expandedvals->lambda-mc-move", -1, expand1
->elmcmove
, expand2
->elmcmove
);
672 cmp_int(fp
, "inputrec->expandedvals->lmc-repeats", -1, expand1
->lmc_repeats
, expand2
->lmc_repeats
);
673 cmp_int(fp
, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1
->gibbsdeltalam
, expand2
->gibbsdeltalam
);
674 cmp_int(fp
, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1
->lmc_forced_nstart
, expand2
->lmc_forced_nstart
);
675 cmp_int(fp
, "inputrec->expandedvals->lambda-weights-equil", -1, expand1
->elmceq
, expand2
->elmceq
);
676 cmp_int(fp
, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1
->equil_n_at_lam
, expand2
->equil_n_at_lam
);
677 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1
->equil_samples
, expand2
->equil_samples
);
678 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1
->equil_steps
, expand2
->equil_steps
);
679 cmp_real(fp
, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1
->equil_wl_delta
, expand2
->equil_wl_delta
, ftol
, abstol
);
680 cmp_real(fp
, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1
->equil_ratio
, expand2
->equil_ratio
, ftol
, abstol
);
681 cmp_bool(fp
, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1
->bSymmetrizedTMatrix
, expand2
->bSymmetrizedTMatrix
);
682 cmp_int(fp
, "inputrec->expandedvals->nstTij", -1, expand1
->nstTij
, expand2
->nstTij
);
683 cmp_int(fp
, "inputrec->expandedvals->mininum-var-min", -1, expand1
->minvarmin
, expand2
->minvarmin
); /*default is reasonable */
684 cmp_int(fp
, "inputrec->expandedvals->weight-c-range", -1, expand1
->c_range
, expand2
->c_range
); /* default is just C=0 */
685 cmp_real(fp
, "inputrec->expandedvals->wl-scale", -1, expand1
->wl_scale
, expand2
->wl_scale
, ftol
, abstol
);
686 cmp_real(fp
, "inputrec->expandedvals->init-wl-delta", -1, expand1
->init_wl_delta
, expand2
->init_wl_delta
, ftol
, abstol
);
687 cmp_real(fp
, "inputrec->expandedvals->wl-ratio", -1, expand1
->wl_ratio
, expand2
->wl_ratio
, ftol
, abstol
);
688 cmp_int(fp
, "inputrec->expandedvals->nstexpanded", -1, expand1
->nstexpanded
, expand2
->nstexpanded
);
689 cmp_int(fp
, "inputrec->expandedvals->lmc-seed", -1, expand1
->lmc_seed
, expand2
->lmc_seed
);
690 cmp_real(fp
, "inputrec->expandedvals->mc-temperature", -1, expand1
->mc_temp
, expand2
->mc_temp
, ftol
, abstol
);
693 static void cmp_fepvals(FILE *fp
, t_lambda
*fep1
, t_lambda
*fep2
, real ftol
, real abstol
)
696 cmp_int(fp
, "inputrec->nstdhdl", -1, fep1
->nstdhdl
, fep2
->nstdhdl
);
697 cmp_double(fp
, "inputrec->fepvals->init_fep_state", -1, fep1
->init_fep_state
, fep2
->init_fep_state
, ftol
, abstol
);
698 cmp_double(fp
, "inputrec->fepvals->delta_lambda", -1, fep1
->delta_lambda
, fep2
->delta_lambda
, ftol
, abstol
);
699 cmp_int(fp
, "inputrec->fepvals->n_lambda", -1, fep1
->n_lambda
, fep2
->n_lambda
);
700 for (i
= 0; i
< efptNR
; i
++)
702 for (j
= 0; j
< std::min(fep1
->n_lambda
, fep2
->n_lambda
); j
++)
704 cmp_double(fp
, "inputrec->fepvals->all_lambda", -1, fep1
->all_lambda
[i
][j
], fep2
->all_lambda
[i
][j
], ftol
, abstol
);
707 cmp_int(fp
, "inputrec->fepvals->lambda_neighbors", 1, fep1
->lambda_neighbors
,
708 fep2
->lambda_neighbors
);
709 cmp_real(fp
, "inputrec->fepvals->sc_alpha", -1, fep1
->sc_alpha
, fep2
->sc_alpha
, ftol
, abstol
);
710 cmp_int(fp
, "inputrec->fepvals->sc_power", -1, fep1
->sc_power
, fep2
->sc_power
);
711 cmp_real(fp
, "inputrec->fepvals->sc_r_power", -1, fep1
->sc_r_power
, fep2
->sc_r_power
, ftol
, abstol
);
712 cmp_real(fp
, "inputrec->fepvals->sc_sigma", -1, fep1
->sc_sigma
, fep2
->sc_sigma
, ftol
, abstol
);
713 cmp_int(fp
, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1
->edHdLPrintEnergy
, fep1
->edHdLPrintEnergy
);
714 cmp_bool(fp
, "inputrec->fepvals->bScCoul", -1, fep1
->bScCoul
, fep1
->bScCoul
);
715 cmp_int(fp
, "inputrec->separate_dhdl_file", -1, fep1
->separate_dhdl_file
, fep2
->separate_dhdl_file
);
716 cmp_int(fp
, "inputrec->dhdl_derivatives", -1, fep1
->dhdl_derivatives
, fep2
->dhdl_derivatives
);
717 cmp_int(fp
, "inputrec->dh_hist_size", -1, fep1
->dh_hist_size
, fep2
->dh_hist_size
);
718 cmp_double(fp
, "inputrec->dh_hist_spacing", -1, fep1
->dh_hist_spacing
, fep2
->dh_hist_spacing
, ftol
, abstol
);
721 static void cmp_inputrec(FILE *fp
, t_inputrec
*ir1
, t_inputrec
*ir2
, real ftol
, real abstol
)
723 fprintf(fp
, "comparing inputrec\n");
725 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
726 * of warnings. Maybe it will change in future versions, but for the
727 * moment I've spelled them out instead. /EL 000820
728 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
729 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
730 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
732 cmp_int(fp
, "inputrec->eI", -1, ir1
->eI
, ir2
->eI
);
733 cmp_int64(fp
, "inputrec->nsteps", ir1
->nsteps
, ir2
->nsteps
);
734 cmp_int64(fp
, "inputrec->init_step", ir1
->init_step
, ir2
->init_step
);
735 cmp_int(fp
, "inputrec->simulation_part", -1, ir1
->simulation_part
, ir2
->simulation_part
);
736 cmp_int(fp
, "inputrec->ePBC", -1, ir1
->ePBC
, ir2
->ePBC
);
737 cmp_int(fp
, "inputrec->bPeriodicMols", -1, ir1
->bPeriodicMols
, ir2
->bPeriodicMols
);
738 cmp_int(fp
, "inputrec->cutoff_scheme", -1, ir1
->cutoff_scheme
, ir2
->cutoff_scheme
);
739 cmp_int(fp
, "inputrec->ns_type", -1, ir1
->ns_type
, ir2
->ns_type
);
740 cmp_int(fp
, "inputrec->nstlist", -1, ir1
->nstlist
, ir2
->nstlist
);
741 cmp_int(fp
, "inputrec->nstcomm", -1, ir1
->nstcomm
, ir2
->nstcomm
);
742 cmp_int(fp
, "inputrec->comm_mode", -1, ir1
->comm_mode
, ir2
->comm_mode
);
743 cmp_int(fp
, "inputrec->nstlog", -1, ir1
->nstlog
, ir2
->nstlog
);
744 cmp_int(fp
, "inputrec->nstxout", -1, ir1
->nstxout
, ir2
->nstxout
);
745 cmp_int(fp
, "inputrec->nstvout", -1, ir1
->nstvout
, ir2
->nstvout
);
746 cmp_int(fp
, "inputrec->nstfout", -1, ir1
->nstfout
, ir2
->nstfout
);
747 cmp_int(fp
, "inputrec->nstcalcenergy", -1, ir1
->nstcalcenergy
, ir2
->nstcalcenergy
);
748 cmp_int(fp
, "inputrec->nstenergy", -1, ir1
->nstenergy
, ir2
->nstenergy
);
749 cmp_int(fp
, "inputrec->nstxout_compressed", -1, ir1
->nstxout_compressed
, ir2
->nstxout_compressed
);
750 cmp_double(fp
, "inputrec->init_t", -1, ir1
->init_t
, ir2
->init_t
, ftol
, abstol
);
751 cmp_double(fp
, "inputrec->delta_t", -1, ir1
->delta_t
, ir2
->delta_t
, ftol
, abstol
);
752 cmp_real(fp
, "inputrec->x_compression_precision", -1, ir1
->x_compression_precision
, ir2
->x_compression_precision
, ftol
, abstol
);
753 cmp_real(fp
, "inputrec->fourierspacing", -1, ir1
->fourier_spacing
, ir2
->fourier_spacing
, ftol
, abstol
);
754 cmp_int(fp
, "inputrec->nkx", -1, ir1
->nkx
, ir2
->nkx
);
755 cmp_int(fp
, "inputrec->nky", -1, ir1
->nky
, ir2
->nky
);
756 cmp_int(fp
, "inputrec->nkz", -1, ir1
->nkz
, ir2
->nkz
);
757 cmp_int(fp
, "inputrec->pme_order", -1, ir1
->pme_order
, ir2
->pme_order
);
758 cmp_real(fp
, "inputrec->ewald_rtol", -1, ir1
->ewald_rtol
, ir2
->ewald_rtol
, ftol
, abstol
);
759 cmp_int(fp
, "inputrec->ewald_geometry", -1, ir1
->ewald_geometry
, ir2
->ewald_geometry
);
760 cmp_real(fp
, "inputrec->epsilon_surface", -1, ir1
->epsilon_surface
, ir2
->epsilon_surface
, ftol
, abstol
);
761 cmp_int(fp
, "inputrec->bContinuation", -1, ir1
->bContinuation
, ir2
->bContinuation
);
762 cmp_int(fp
, "inputrec->bShakeSOR", -1, ir1
->bShakeSOR
, ir2
->bShakeSOR
);
763 cmp_int(fp
, "inputrec->etc", -1, ir1
->etc
, ir2
->etc
);
764 cmp_int(fp
, "inputrec->bPrintNHChains", -1, ir1
->bPrintNHChains
, ir2
->bPrintNHChains
);
765 cmp_int(fp
, "inputrec->epc", -1, ir1
->epc
, ir2
->epc
);
766 cmp_int(fp
, "inputrec->epct", -1, ir1
->epct
, ir2
->epct
);
767 cmp_real(fp
, "inputrec->tau_p", -1, ir1
->tau_p
, ir2
->tau_p
, ftol
, abstol
);
768 cmp_rvec(fp
, "inputrec->ref_p(x)", -1, ir1
->ref_p
[XX
], ir2
->ref_p
[XX
], ftol
, abstol
);
769 cmp_rvec(fp
, "inputrec->ref_p(y)", -1, ir1
->ref_p
[YY
], ir2
->ref_p
[YY
], ftol
, abstol
);
770 cmp_rvec(fp
, "inputrec->ref_p(z)", -1, ir1
->ref_p
[ZZ
], ir2
->ref_p
[ZZ
], ftol
, abstol
);
771 cmp_rvec(fp
, "inputrec->compress(x)", -1, ir1
->compress
[XX
], ir2
->compress
[XX
], ftol
, abstol
);
772 cmp_rvec(fp
, "inputrec->compress(y)", -1, ir1
->compress
[YY
], ir2
->compress
[YY
], ftol
, abstol
);
773 cmp_rvec(fp
, "inputrec->compress(z)", -1, ir1
->compress
[ZZ
], ir2
->compress
[ZZ
], ftol
, abstol
);
774 cmp_int(fp
, "refcoord_scaling", -1, ir1
->refcoord_scaling
, ir2
->refcoord_scaling
);
775 cmp_rvec(fp
, "inputrec->posres_com", -1, ir1
->posres_com
, ir2
->posres_com
, ftol
, abstol
);
776 cmp_rvec(fp
, "inputrec->posres_comB", -1, ir1
->posres_comB
, ir2
->posres_comB
, ftol
, abstol
);
777 cmp_real(fp
, "inputrec->verletbuf_tol", -1, ir1
->verletbuf_tol
, ir2
->verletbuf_tol
, ftol
, abstol
);
778 cmp_real(fp
, "inputrec->rlist", -1, ir1
->rlist
, ir2
->rlist
, ftol
, abstol
);
779 cmp_real(fp
, "inputrec->rtpi", -1, ir1
->rtpi
, ir2
->rtpi
, ftol
, abstol
);
780 cmp_int(fp
, "inputrec->coulombtype", -1, ir1
->coulombtype
, ir2
->coulombtype
);
781 cmp_int(fp
, "inputrec->coulomb_modifier", -1, ir1
->coulomb_modifier
, ir2
->coulomb_modifier
);
782 cmp_real(fp
, "inputrec->rcoulomb_switch", -1, ir1
->rcoulomb_switch
, ir2
->rcoulomb_switch
, ftol
, abstol
);
783 cmp_real(fp
, "inputrec->rcoulomb", -1, ir1
->rcoulomb
, ir2
->rcoulomb
, ftol
, abstol
);
784 cmp_int(fp
, "inputrec->vdwtype", -1, ir1
->vdwtype
, ir2
->vdwtype
);
785 cmp_int(fp
, "inputrec->vdw_modifier", -1, ir1
->vdw_modifier
, ir2
->vdw_modifier
); cmp_real(fp
, "inputrec->rvdw_switch", -1, ir1
->rvdw_switch
, ir2
->rvdw_switch
, ftol
, abstol
);
786 cmp_real(fp
, "inputrec->rvdw", -1, ir1
->rvdw
, ir2
->rvdw
, ftol
, abstol
);
787 cmp_real(fp
, "inputrec->epsilon_r", -1, ir1
->epsilon_r
, ir2
->epsilon_r
, ftol
, abstol
);
788 cmp_real(fp
, "inputrec->epsilon_rf", -1, ir1
->epsilon_rf
, ir2
->epsilon_rf
, ftol
, abstol
);
789 cmp_real(fp
, "inputrec->tabext", -1, ir1
->tabext
, ir2
->tabext
, ftol
, abstol
);
790 cmp_int(fp
, "inputrec->implicit_solvent", -1, ir1
->implicit_solvent
, ir2
->implicit_solvent
);
791 cmp_int(fp
, "inputrec->gb_algorithm", -1, ir1
->gb_algorithm
, ir2
->gb_algorithm
);
792 cmp_int(fp
, "inputrec->nstgbradii", -1, ir1
->nstgbradii
, ir2
->nstgbradii
);
793 cmp_real(fp
, "inputrec->rgbradii", -1, ir1
->rgbradii
, ir2
->rgbradii
, ftol
, abstol
);
794 cmp_real(fp
, "inputrec->gb_saltconc", -1, ir1
->gb_saltconc
, ir2
->gb_saltconc
, ftol
, abstol
);
795 cmp_real(fp
, "inputrec->gb_epsilon_solvent", -1, ir1
->gb_epsilon_solvent
, ir2
->gb_epsilon_solvent
, ftol
, abstol
);
796 cmp_real(fp
, "inputrec->gb_obc_alpha", -1, ir1
->gb_obc_alpha
, ir2
->gb_obc_alpha
, ftol
, abstol
);
797 cmp_real(fp
, "inputrec->gb_obc_beta", -1, ir1
->gb_obc_beta
, ir2
->gb_obc_beta
, ftol
, abstol
);
798 cmp_real(fp
, "inputrec->gb_obc_gamma", -1, ir1
->gb_obc_gamma
, ir2
->gb_obc_gamma
, ftol
, abstol
);
799 cmp_real(fp
, "inputrec->gb_dielectric_offset", -1, ir1
->gb_dielectric_offset
, ir2
->gb_dielectric_offset
, ftol
, abstol
);
800 cmp_int(fp
, "inputrec->sa_algorithm", -1, ir1
->sa_algorithm
, ir2
->sa_algorithm
);
801 cmp_real(fp
, "inputrec->sa_surface_tension", -1, ir1
->sa_surface_tension
, ir2
->sa_surface_tension
, ftol
, abstol
);
803 cmp_int(fp
, "inputrec->eDispCorr", -1, ir1
->eDispCorr
, ir2
->eDispCorr
);
804 cmp_real(fp
, "inputrec->shake_tol", -1, ir1
->shake_tol
, ir2
->shake_tol
, ftol
, abstol
);
805 cmp_int(fp
, "inputrec->efep", -1, ir1
->efep
, ir2
->efep
);
806 cmp_fepvals(fp
, ir1
->fepvals
, ir2
->fepvals
, ftol
, abstol
);
807 cmp_int(fp
, "inputrec->bSimTemp", -1, ir1
->bSimTemp
, ir2
->bSimTemp
);
808 if ((ir1
->bSimTemp
== ir2
->bSimTemp
) && (ir1
->bSimTemp
))
810 cmp_simtempvals(fp
, ir1
->simtempvals
, ir2
->simtempvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
812 cmp_int(fp
, "inputrec->bExpanded", -1, ir1
->bExpanded
, ir2
->bExpanded
);
813 if ((ir1
->bExpanded
== ir2
->bExpanded
) && (ir1
->bExpanded
))
815 cmp_expandedvals(fp
, ir1
->expandedvals
, ir2
->expandedvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
817 cmp_int(fp
, "inputrec->nwall", -1, ir1
->nwall
, ir2
->nwall
);
818 cmp_int(fp
, "inputrec->wall_type", -1, ir1
->wall_type
, ir2
->wall_type
);
819 cmp_int(fp
, "inputrec->wall_atomtype[0]", -1, ir1
->wall_atomtype
[0], ir2
->wall_atomtype
[0]);
820 cmp_int(fp
, "inputrec->wall_atomtype[1]", -1, ir1
->wall_atomtype
[1], ir2
->wall_atomtype
[1]);
821 cmp_real(fp
, "inputrec->wall_density[0]", -1, ir1
->wall_density
[0], ir2
->wall_density
[0], ftol
, abstol
);
822 cmp_real(fp
, "inputrec->wall_density[1]", -1, ir1
->wall_density
[1], ir2
->wall_density
[1], ftol
, abstol
);
823 cmp_real(fp
, "inputrec->wall_ewald_zfac", -1, ir1
->wall_ewald_zfac
, ir2
->wall_ewald_zfac
, ftol
, abstol
);
825 cmp_bool(fp
, "inputrec->bPull", -1, ir1
->bPull
, ir2
->bPull
);
826 if (ir1
->bPull
&& ir2
->bPull
)
831 cmp_int(fp
, "inputrec->eDisre", -1, ir1
->eDisre
, ir2
->eDisre
);
832 cmp_real(fp
, "inputrec->dr_fc", -1, ir1
->dr_fc
, ir2
->dr_fc
, ftol
, abstol
);
833 cmp_int(fp
, "inputrec->eDisreWeighting", -1, ir1
->eDisreWeighting
, ir2
->eDisreWeighting
);
834 cmp_int(fp
, "inputrec->bDisreMixed", -1, ir1
->bDisreMixed
, ir2
->bDisreMixed
);
835 cmp_int(fp
, "inputrec->nstdisreout", -1, ir1
->nstdisreout
, ir2
->nstdisreout
);
836 cmp_real(fp
, "inputrec->dr_tau", -1, ir1
->dr_tau
, ir2
->dr_tau
, ftol
, abstol
);
837 cmp_real(fp
, "inputrec->orires_fc", -1, ir1
->orires_fc
, ir2
->orires_fc
, ftol
, abstol
);
838 cmp_real(fp
, "inputrec->orires_tau", -1, ir1
->orires_tau
, ir2
->orires_tau
, ftol
, abstol
);
839 cmp_int(fp
, "inputrec->nstorireout", -1, ir1
->nstorireout
, ir2
->nstorireout
);
840 cmp_real(fp
, "inputrec->em_stepsize", -1, ir1
->em_stepsize
, ir2
->em_stepsize
, ftol
, abstol
);
841 cmp_real(fp
, "inputrec->em_tol", -1, ir1
->em_tol
, ir2
->em_tol
, ftol
, abstol
);
842 cmp_int(fp
, "inputrec->niter", -1, ir1
->niter
, ir2
->niter
);
843 cmp_real(fp
, "inputrec->fc_stepsize", -1, ir1
->fc_stepsize
, ir2
->fc_stepsize
, ftol
, abstol
);
844 cmp_int(fp
, "inputrec->nstcgsteep", -1, ir1
->nstcgsteep
, ir2
->nstcgsteep
);
845 cmp_int(fp
, "inputrec->nbfgscorr", 0, ir1
->nbfgscorr
, ir2
->nbfgscorr
);
846 cmp_int(fp
, "inputrec->eConstrAlg", -1, ir1
->eConstrAlg
, ir2
->eConstrAlg
);
847 cmp_int(fp
, "inputrec->nProjOrder", -1, ir1
->nProjOrder
, ir2
->nProjOrder
);
848 cmp_real(fp
, "inputrec->LincsWarnAngle", -1, ir1
->LincsWarnAngle
, ir2
->LincsWarnAngle
, ftol
, abstol
);
849 cmp_int(fp
, "inputrec->nLincsIter", -1, ir1
->nLincsIter
, ir2
->nLincsIter
);
850 cmp_real(fp
, "inputrec->bd_fric", -1, ir1
->bd_fric
, ir2
->bd_fric
, ftol
, abstol
);
851 cmp_int64(fp
, "inputrec->ld_seed", ir1
->ld_seed
, ir2
->ld_seed
);
852 cmp_real(fp
, "inputrec->cos_accel", -1, ir1
->cos_accel
, ir2
->cos_accel
, ftol
, abstol
);
853 cmp_rvec(fp
, "inputrec->deform(a)", -1, ir1
->deform
[XX
], ir2
->deform
[XX
], ftol
, abstol
);
854 cmp_rvec(fp
, "inputrec->deform(b)", -1, ir1
->deform
[YY
], ir2
->deform
[YY
], ftol
, abstol
);
855 cmp_rvec(fp
, "inputrec->deform(c)", -1, ir1
->deform
[ZZ
], ir2
->deform
[ZZ
], ftol
, abstol
);
858 cmp_int(fp
, "inputrec->userint1", -1, ir1
->userint1
, ir2
->userint1
);
859 cmp_int(fp
, "inputrec->userint2", -1, ir1
->userint2
, ir2
->userint2
);
860 cmp_int(fp
, "inputrec->userint3", -1, ir1
->userint3
, ir2
->userint3
);
861 cmp_int(fp
, "inputrec->userint4", -1, ir1
->userint4
, ir2
->userint4
);
862 cmp_real(fp
, "inputrec->userreal1", -1, ir1
->userreal1
, ir2
->userreal1
, ftol
, abstol
);
863 cmp_real(fp
, "inputrec->userreal2", -1, ir1
->userreal2
, ir2
->userreal2
, ftol
, abstol
);
864 cmp_real(fp
, "inputrec->userreal3", -1, ir1
->userreal3
, ir2
->userreal3
, ftol
, abstol
);
865 cmp_real(fp
, "inputrec->userreal4", -1, ir1
->userreal4
, ir2
->userreal4
, ftol
, abstol
);
866 cmp_grpopts(fp
, &(ir1
->opts
), &(ir2
->opts
), ftol
, abstol
);
867 cmp_cosines(fp
, "ex", ir1
->ex
, ir2
->ex
, ftol
, abstol
);
868 cmp_cosines(fp
, "et", ir1
->et
, ir2
->et
, ftol
, abstol
);
871 static void comp_pull_AB(FILE *fp
, pull_params_t
*pull
, real ftol
, real abstol
)
875 for (i
= 0; i
< pull
->ncoord
; i
++)
877 fprintf(fp
, "comparing pull coord %d\n", i
);
878 cmp_real(fp
, "pull-coord->k", -1, pull
->coord
[i
].k
, pull
->coord
[i
].kB
, ftol
, abstol
);
882 static void comp_state(t_state
*st1
, t_state
*st2
,
883 gmx_bool bRMSD
, real ftol
, real abstol
)
887 fprintf(stdout
, "comparing flags\n");
888 cmp_int(stdout
, "flags", -1, st1
->flags
, st2
->flags
);
889 fprintf(stdout
, "comparing box\n");
890 cmp_rvecs(stdout
, "box", DIM
, st1
->box
, st2
->box
, FALSE
, ftol
, abstol
);
891 fprintf(stdout
, "comparing box_rel\n");
892 cmp_rvecs(stdout
, "box_rel", DIM
, st1
->box_rel
, st2
->box_rel
, FALSE
, ftol
, abstol
);
893 fprintf(stdout
, "comparing boxv\n");
894 cmp_rvecs(stdout
, "boxv", DIM
, st1
->boxv
, st2
->boxv
, FALSE
, ftol
, abstol
);
895 if (st1
->flags
& (1<<estSVIR_PREV
))
897 fprintf(stdout
, "comparing shake vir_prev\n");
898 cmp_rvecs_rmstol(stdout
, "svir_prev", DIM
, st1
->svir_prev
, st2
->svir_prev
, ftol
, abstol
);
900 if (st1
->flags
& (1<<estFVIR_PREV
))
902 fprintf(stdout
, "comparing force vir_prev\n");
903 cmp_rvecs_rmstol(stdout
, "fvir_prev", DIM
, st1
->fvir_prev
, st2
->fvir_prev
, ftol
, abstol
);
905 if (st1
->flags
& (1<<estPRES_PREV
))
907 fprintf(stdout
, "comparing prev_pres\n");
908 cmp_rvecs_rmstol(stdout
, "pres_prev", DIM
, st1
->pres_prev
, st2
->pres_prev
, ftol
, abstol
);
910 cmp_int(stdout
, "ngtc", -1, st1
->ngtc
, st2
->ngtc
);
911 cmp_int(stdout
, "nhchainlength", -1, st1
->nhchainlength
, st2
->nhchainlength
);
912 if (st1
->ngtc
== st2
->ngtc
&& st1
->nhchainlength
== st2
->nhchainlength
)
914 for (i
= 0; i
< st1
->ngtc
; i
++)
916 nc
= i
*st1
->nhchainlength
;
917 for (j
= 0; j
< nc
; j
++)
919 cmp_real(stdout
, "nosehoover_xi",
920 i
, st1
->nosehoover_xi
[nc
+j
], st2
->nosehoover_xi
[nc
+j
], ftol
, abstol
);
924 cmp_int(stdout
, "nnhpres", -1, st1
->nnhpres
, st2
->nnhpres
);
925 if (st1
->nnhpres
== st2
->nnhpres
&& st1
->nhchainlength
== st2
->nhchainlength
)
927 for (i
= 0; i
< st1
->nnhpres
; i
++)
929 nc
= i
*st1
->nhchainlength
;
930 for (j
= 0; j
< nc
; j
++)
932 cmp_real(stdout
, "nosehoover_xi",
933 i
, st1
->nhpres_xi
[nc
+j
], st2
->nhpres_xi
[nc
+j
], ftol
, abstol
);
938 cmp_int(stdout
, "natoms", -1, st1
->natoms
, st2
->natoms
);
939 if (st1
->natoms
== st2
->natoms
)
941 if ((st1
->flags
& (1<<estX
)) && (st2
->flags
& (1<<estX
)))
943 fprintf(stdout
, "comparing x\n");
944 cmp_rvecs(stdout
, "x", st1
->natoms
, st1
->x
, st2
->x
, bRMSD
, ftol
, abstol
);
946 if ((st1
->flags
& (1<<estV
)) && (st2
->flags
& (1<<estV
)))
948 fprintf(stdout
, "comparing v\n");
949 cmp_rvecs(stdout
, "v", st1
->natoms
, st1
->v
, st2
->v
, bRMSD
, ftol
, abstol
);
954 void comp_tpx(const char *fn1
, const char *fn2
,
955 gmx_bool bRMSD
, real ftol
, real abstol
)
966 for (i
= 0; i
< (fn2
? 2 : 1); i
++)
968 read_tpx_state(ff
[i
], &(ir
[i
]), &state
[i
], &(mtop
[i
]));
972 cmp_inputrec(stdout
, &ir
[0], &ir
[1], ftol
, abstol
);
973 /* Convert gmx_mtop_t to t_topology.
974 * We should implement direct mtop comparison,
975 * but it might be useful to keep t_topology comparison as an option.
977 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
978 top
[1] = gmx_mtop_t_to_t_topology(&mtop
[1]);
979 cmp_top(stdout
, &top
[0], &top
[1], ftol
, abstol
);
980 cmp_groups(stdout
, &mtop
[0].groups
, &mtop
[1].groups
,
981 mtop
[0].natoms
, mtop
[1].natoms
);
982 comp_state(&state
[0], &state
[1], bRMSD
, ftol
, abstol
);
986 if (ir
[0].efep
== efepNO
)
988 fprintf(stdout
, "inputrec->efep = %s\n", efep_names
[ir
[0].efep
]);
994 comp_pull_AB(stdout
, ir
->pull
, ftol
, abstol
);
996 /* Convert gmx_mtop_t to t_topology.
997 * We should implement direct mtop comparison,
998 * but it might be useful to keep t_topology comparison as an option.
1000 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
1001 cmp_top(stdout
, &top
[0], NULL
, ftol
, abstol
);
1006 void comp_frame(FILE *fp
, t_trxframe
*fr1
, t_trxframe
*fr2
,
1007 gmx_bool bRMSD
, real ftol
, real abstol
)
1010 cmp_int(fp
, "flags", -1, fr1
->flags
, fr2
->flags
);
1011 cmp_int(fp
, "not_ok", -1, fr1
->not_ok
, fr2
->not_ok
);
1012 cmp_int(fp
, "natoms", -1, fr1
->natoms
, fr2
->natoms
);
1013 cmp_real(fp
, "t0", -1, fr1
->t0
, fr2
->t0
, ftol
, abstol
);
1014 if (cmp_bool(fp
, "bTitle", -1, fr1
->bTitle
, fr2
->bTitle
))
1016 cmp_str(fp
, "title", -1, fr1
->title
, fr2
->title
);
1018 if (cmp_bool(fp
, "bStep", -1, fr1
->bStep
, fr2
->bStep
))
1020 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1022 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1023 if (cmp_bool(fp
, "bTime", -1, fr1
->bTime
, fr2
->bTime
))
1025 cmp_real(fp
, "time", -1, fr1
->time
, fr2
->time
, ftol
, abstol
);
1027 if (cmp_bool(fp
, "bLambda", -1, fr1
->bLambda
, fr2
->bLambda
))
1029 cmp_real(fp
, "lambda", -1, fr1
->lambda
, fr2
->lambda
, ftol
, abstol
);
1031 if (cmp_bool(fp
, "bAtoms", -1, fr1
->bAtoms
, fr2
->bAtoms
))
1033 cmp_atoms(fp
, fr1
->atoms
, fr2
->atoms
, ftol
, abstol
);
1035 if (cmp_bool(fp
, "bPrec", -1, fr1
->bPrec
, fr2
->bPrec
))
1037 cmp_real(fp
, "prec", -1, fr1
->prec
, fr2
->prec
, ftol
, abstol
);
1039 if (cmp_bool(fp
, "bX", -1, fr1
->bX
, fr2
->bX
))
1041 cmp_rvecs(fp
, "x", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->x
, fr2
->x
, bRMSD
, ftol
, abstol
);
1043 if (cmp_bool(fp
, "bV", -1, fr1
->bV
, fr2
->bV
))
1045 cmp_rvecs(fp
, "v", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->v
, fr2
->v
, bRMSD
, ftol
, abstol
);
1047 if (cmp_bool(fp
, "bF", -1, fr1
->bF
, fr2
->bF
))
1051 cmp_rvecs(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, bRMSD
, ftol
, abstol
);
1055 cmp_rvecs_rmstol(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, ftol
, abstol
);
1058 if (cmp_bool(fp
, "bBox", -1, fr1
->bBox
, fr2
->bBox
))
1060 cmp_rvecs(fp
, "box", 3, fr1
->box
, fr2
->box
, FALSE
, ftol
, abstol
);
1064 void comp_trx(const gmx_output_env_t
*oenv
, const char *fn1
, const char *fn2
,
1065 gmx_bool bRMSD
, real ftol
, real abstol
)
1070 t_trxstatus
*status
[2];
1075 fprintf(stderr
, "Comparing trajectory files %s and %s\n", fn1
, fn2
);
1076 for (i
= 0; i
< 2; i
++)
1078 b
[i
] = read_first_frame(oenv
, &status
[i
], fn
[i
], &fr
[i
], TRX_READ_X
|TRX_READ_V
|TRX_READ_F
);
1085 comp_frame(stdout
, &(fr
[0]), &(fr
[1]), bRMSD
, ftol
, abstol
);
1087 for (i
= 0; i
< 2; i
++)
1089 b
[i
] = read_next_frame(oenv
, status
[i
], &fr
[i
]);
1092 while (b
[0] && b
[1]);
1094 for (i
= 0; i
< 2; i
++)
1096 if (b
[i
] && !b
[1-i
])
1098 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn
[1-i
], fn
[i
]);
1100 close_trj(status
[i
]);
1105 fprintf(stdout
, "\nBoth files read correctly\n");
1109 static real
ener_tensor_diag(int n
, int *ind1
, int *ind2
,
1112 t_energy e1
[], t_energy e2
[])
1121 d2
= tensi
[i
] - d1
*DIM
;
1123 /* Find the diagonal elements d1 and d2 */
1124 len
= std::strlen(enm1
[ind1
[i
]].name
);
1128 for (j
= 0; j
< n
; j
++)
1130 if (tensi
[j
] >= 0 &&
1131 std::strlen(enm1
[ind1
[j
]].name
) == len
&&
1132 std::strncmp(enm1
[ind1
[i
]].name
, enm1
[ind1
[j
]].name
, len
-2) == 0 &&
1133 (tensi
[j
] == d1
*DIM
+d1
|| tensi
[j
] == d2
*DIM
+d2
))
1135 prod1
*= fabs(e1
[ind1
[j
]].e
);
1136 prod2
*= fabs(e2
[ind2
[j
]].e
);
1143 return 0.5*(std::sqrt(prod1
) + std::sqrt(prod2
));
1151 static gmx_bool
enernm_equal(const char *nm1
, const char *nm2
)
1155 len1
= std::strlen(nm1
);
1156 len2
= std::strlen(nm2
);
1158 /* Remove " (bar)" at the end of a name */
1159 if (len1
> 6 && std::strcmp(nm1
+len1
-6, " (bar)") == 0)
1163 if (len2
> 6 && std::strcmp(nm2
+len2
-6, " (bar)") == 0)
1168 return (len1
== len2
&& gmx_strncasecmp(nm1
, nm2
, len1
) == 0);
1171 static void cmp_energies(FILE *fp
, int step1
, int step2
,
1172 t_energy e1
[], t_energy e2
[],
1174 real ftol
, real abstol
,
1175 int nre
, int *ind1
, int *ind2
, int maxener
)
1178 int *tensi
, len
, d1
, d2
;
1179 real ftol_i
, abstol_i
;
1181 snew(tensi
, maxener
);
1182 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1183 for (i
= 0; (i
< maxener
); i
++)
1187 len
= std::strlen(enm1
[ii
].name
);
1188 if (len
> 3 && enm1
[ii
].name
[len
-3] == '-')
1190 d1
= enm1
[ii
].name
[len
-2] - 'X';
1191 d2
= enm1
[ii
].name
[len
-1] - 'X';
1192 if (d1
>= 0 && d1
< DIM
&&
1193 d2
>= 0 && d2
< DIM
)
1195 tensi
[i
] = d1
*DIM
+ d2
;
1200 for (i
= 0; (i
< maxener
); i
++)
1202 /* Check if this is an off-diagonal tensor element */
1203 if (tensi
[i
] >= 0 && tensi
[i
] != 0 && tensi
[i
] != 4 && tensi
[i
] != 8)
1205 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1207 /* Do the relative tolerance through an absolute tolerance times
1208 * the size of diagonal components of the tensor.
1210 abstol_i
= ftol
*ener_tensor_diag(nre
, ind1
, ind2
, enm1
, tensi
, i
, e1
, e2
);
1213 fprintf(debug
, "tensor '%s' val %f diag %f\n",
1214 enm1
[i
].name
, e1
[i
].e
, abstol_i
/ftol
);
1218 /* We found a diagonal, we need to check with the minimum tolerance */
1219 abstol_i
= std::min(abstol_i
, abstol
);
1223 /* We did not find a diagonal, ignore the relative tolerance check */
1232 if (!equal_real(e1
[ind1
[i
]].e
, e2
[ind2
[i
]].e
, ftol_i
, abstol_i
))
1234 fprintf(fp
, "%-15s step %3d: %12g, step %3d: %12g\n",
1236 step1
, e1
[ind1
[i
]].e
,
1237 step2
, e2
[ind2
[i
]].e
);
1245 static void cmp_disres(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1248 char bav
[64], bt
[64], bs
[22];
1250 cmp_int(stdout
, "ndisre", -1, fr1
->ndisre
, fr2
->ndisre
);
1251 if ((fr1
->ndisre
== fr2
->ndisre
) && (fr1
->ndisre
> 0))
1253 sprintf(bav
, "step %s: disre rav", gmx_step_str(fr1
->step
, bs
));
1254 sprintf(bt
, "step %s: disre rt", gmx_step_str(fr1
->step
, bs
));
1255 for (i
= 0; (i
< fr1
->ndisre
); i
++)
1257 cmp_real(stdout
, bav
, i
, fr1
->disre_rm3tav
[i
], fr2
->disre_rm3tav
[i
], ftol
, abstol
);
1258 cmp_real(stdout
, bt
, i
, fr1
->disre_rt
[i
], fr2
->disre_rt
[i
], ftol
, abstol
);
1264 static void cmp_eblocks(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1267 char buf
[64], bs
[22];
1269 cmp_int(stdout
, "nblock", -1, fr1
->nblock
, fr2
->nblock
);
1270 if ((fr1
->nblock
== fr2
->nblock
) && (fr1
->nblock
> 0))
1272 for (j
= 0; (j
< fr1
->nblock
); j
++)
1274 t_enxblock
*b1
, *b2
; /* convenience vars */
1276 b1
= &(fr1
->block
[j
]);
1277 b2
= &(fr2
->block
[j
]);
1279 sprintf(buf
, "step %s: block[%d]", gmx_step_str(fr1
->step
, bs
), j
);
1280 cmp_int(stdout
, buf
, -1, b1
->nsub
, b2
->nsub
);
1281 cmp_int(stdout
, buf
, -1, b1
->id
, b2
->id
);
1283 if ( (b1
->nsub
== b2
->nsub
) && (b1
->id
== b2
->id
) )
1285 for (i
= 0; i
< b1
->nsub
; i
++)
1287 t_enxsubblock
*s1
, *s2
;
1292 cmp_int(stdout
, buf
, -1, (int)s1
->type
, (int)s2
->type
);
1293 cmp_int64(stdout
, buf
, s1
->nr
, s2
->nr
);
1295 if ((s1
->type
== s2
->type
) && (s1
->nr
== s2
->nr
))
1299 case xdr_datatype_float
:
1300 for (k
= 0; k
< s1
->nr
; k
++)
1302 cmp_float(stdout
, buf
, i
,
1303 s1
->fval
[k
], s2
->fval
[k
],
1307 case xdr_datatype_double
:
1308 for (k
= 0; k
< s1
->nr
; k
++)
1310 cmp_double(stdout
, buf
, i
,
1311 s1
->dval
[k
], s2
->dval
[k
],
1315 case xdr_datatype_int
:
1316 for (k
= 0; k
< s1
->nr
; k
++)
1318 cmp_int(stdout
, buf
, i
,
1319 s1
->ival
[k
], s2
->ival
[k
]);
1322 case xdr_datatype_int64
:
1323 for (k
= 0; k
< s1
->nr
; k
++)
1325 cmp_int64(stdout
, buf
,
1326 s1
->lval
[k
], s2
->lval
[k
]);
1329 case xdr_datatype_char
:
1330 for (k
= 0; k
< s1
->nr
; k
++)
1332 cmp_uc(stdout
, buf
, i
,
1333 s1
->cval
[k
], s2
->cval
[k
]);
1336 case xdr_datatype_string
:
1337 for (k
= 0; k
< s1
->nr
; k
++)
1339 cmp_str(stdout
, buf
, i
,
1340 s1
->sval
[k
], s2
->sval
[k
]);
1344 gmx_incons("Unknown data type!!");
1353 void comp_enx(const char *fn1
, const char *fn2
, real ftol
, real abstol
, const char *lastener
)
1355 int nre
, nre1
, nre2
;
1356 ener_file_t in1
, in2
;
1357 int i
, j
, maxener
, *ind1
, *ind2
, *have
;
1358 gmx_enxnm_t
*enm1
= NULL
, *enm2
= NULL
;
1359 t_enxframe
*fr1
, *fr2
;
1362 fprintf(stdout
, "comparing energy file %s and %s\n\n", fn1
, fn2
);
1364 in1
= open_enx(fn1
, "r");
1365 in2
= open_enx(fn2
, "r");
1366 do_enxnms(in1
, &nre1
, &enm1
);
1367 do_enxnms(in2
, &nre2
, &enm2
);
1370 fprintf(stdout
, "There are %d and %d terms in the energy files\n\n",
1375 fprintf(stdout
, "There are %d terms in the energy files\n\n", nre1
);
1382 for (i
= 0; i
< nre1
; i
++)
1384 for (j
= 0; j
< nre2
; j
++)
1386 if (enernm_equal(enm1
[i
].name
, enm2
[j
].name
))
1395 if (nre
== 0 || ind1
[nre
-1] != i
)
1397 cmp_str(stdout
, "enm", i
, enm1
[i
].name
, "-");
1400 for (i
= 0; i
< nre2
; i
++)
1404 cmp_str(stdout
, "enm", i
, "-", enm2
[i
].name
);
1409 for (i
= 0; i
< nre
; i
++)
1411 if ((lastener
!= NULL
) && (std::strstr(enm1
[i
].name
, lastener
) != NULL
))
1418 fprintf(stdout
, "There are %d terms to compare in the energy files\n\n",
1421 for (i
= 0; i
< maxener
; i
++)
1423 cmp_str(stdout
, "unit", i
, enm1
[ind1
[i
]].unit
, enm2
[ind2
[i
]].unit
);
1430 b1
= do_enx(in1
, fr1
);
1431 b2
= do_enx(in2
, fr2
);
1434 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn2
, fn1
);
1438 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn1
, fn2
);
1440 else if (!b1
&& !b2
)
1442 fprintf(stdout
, "\nFiles read successfully\n");
1446 cmp_real(stdout
, "t", -1, fr1
->t
, fr2
->t
, ftol
, abstol
);
1447 cmp_int(stdout
, "step", -1, fr1
->step
, fr2
->step
);
1448 /* We don't want to print the nre mismatch for every frame */
1449 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1450 if ((fr1
->nre
>= nre
) && (fr2
->nre
>= nre
))
1452 cmp_energies(stdout
, fr1
->step
, fr1
->step
, fr1
->ener
, fr2
->ener
,
1453 enm1
, ftol
, abstol
, nre
, ind1
, ind2
, maxener
);
1455 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1456 cmp_eblocks(fr1
, fr2
, ftol
, abstol
);