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/mdtypes/md_enums.h"
52 #include "gromacs/mdtypes/pull-params.h"
53 #include "gromacs/topology/ifunc.h"
54 #include "gromacs/topology/mtop_util.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
60 #include "gromacs/utility/stringutil.h"
62 static void cmp_int(FILE *fp
, const char *s
, int index
, int i1
, int i2
)
68 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
72 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
77 static void cmp_int64(FILE *fp
, const char *s
, gmx_int64_t i1
, gmx_int64_t i2
)
81 fprintf(fp
, "%s (", s
);
82 fprintf(fp
, "%" GMX_PRId64
, i1
);
84 fprintf(fp
, "%" GMX_PRId64
, i2
);
89 static void cmp_us(FILE *fp
, const char *s
, int index
, unsigned short i1
, unsigned short i2
)
95 fprintf(fp
, "%s[%d] (%hu - %hu)\n", s
, index
, i1
, i2
);
99 fprintf(fp
, "%s (%hu - %hu)\n", s
, i1
, i2
);
104 static void cmp_uc(FILE *fp
, const char *s
, int index
, unsigned char i1
, unsigned char i2
)
110 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
114 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
119 static gmx_bool
cmp_bool(FILE *fp
, const char *s
, int index
, gmx_bool b1
, gmx_bool b2
)
141 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
,
142 gmx::boolToString(b1
), gmx::boolToString(b2
));
146 fprintf(fp
, "%s (%s - %s)\n", s
,
147 gmx::boolToString(b1
), gmx::boolToString(b2
));
153 static void cmp_str(FILE *fp
, const char *s
, int index
,
154 const char *s1
, const char *s2
)
156 if (std::strcmp(s1
, s2
) != 0)
160 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
, s1
, s2
);
164 fprintf(fp
, "%s (%s - %s)\n", s
, s1
, s2
);
169 static gmx_bool
equal_real(real i1
, real i2
, real ftol
, real abstol
)
171 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
174 static gmx_bool
equal_float(float i1
, float i2
, float ftol
, float abstol
)
176 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
179 static gmx_bool
equal_double(double i1
, double i2
, real ftol
, real abstol
)
181 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
185 cmp_real(FILE *fp
, const char *s
, int index
, real i1
, real i2
, real ftol
, real abstol
)
187 if (!equal_real(i1
, i2
, ftol
, abstol
))
191 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
195 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
201 cmp_float(FILE *fp
, const char *s
, int index
, float i1
, float i2
, float ftol
, float abstol
)
203 if (!equal_float(i1
, i2
, ftol
, abstol
))
207 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
211 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
219 cmp_double(FILE *fp
, const char *s
, int index
, double i1
, double i2
, double ftol
, double abstol
)
221 if (!equal_double(i1
, i2
, ftol
, abstol
))
225 fprintf(fp
, "%s[%2d] (%16.9e - %16.9e)\n", s
, index
, i1
, i2
);
229 fprintf(fp
, "%s (%16.9e - %16.9e)\n", s
, i1
, i2
);
234 static void cmp_rvec(FILE *fp
, const char *s
, int index
, rvec i1
, rvec i2
, real ftol
, real abstol
)
236 if (!equal_real(i1
[XX
], i2
[XX
], ftol
, abstol
) ||
237 !equal_real(i1
[YY
], i2
[YY
], ftol
, abstol
) ||
238 !equal_real(i1
[ZZ
], i2
[ZZ
], ftol
, abstol
))
242 fprintf(fp
, "%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
243 s
, index
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
247 fprintf(fp
, "%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
248 s
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
253 static void cmp_ivec(FILE *fp
, const char *s
, int index
, ivec i1
, ivec i2
)
255 if ((i1
[XX
] != i2
[XX
]) || (i1
[YY
] != i2
[YY
]) || (i1
[ZZ
] != i2
[ZZ
]))
259 fprintf(fp
, "%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
, index
,
260 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
264 fprintf(fp
, "%s (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
,
265 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
270 static void cmp_ilist(FILE *fp
, int ftype
, t_ilist
*il1
, t_ilist
*il2
)
275 fprintf(fp
, "comparing ilist %s\n", interaction_function
[ftype
].name
);
276 sprintf(buf
, "%s->nr", interaction_function
[ftype
].name
);
277 cmp_int(fp
, buf
, -1, il1
->nr
, il2
->nr
);
278 sprintf(buf
, "%s->iatoms", interaction_function
[ftype
].name
);
279 if (((il1
->nr
> 0) && (!il1
->iatoms
)) ||
280 ((il2
->nr
> 0) && (!il2
->iatoms
)) ||
281 ((il1
->nr
!= il2
->nr
)))
283 fprintf(fp
, "Comparing radically different topologies - %s is different\n",
288 for (i
= 0; (i
< il1
->nr
); i
++)
290 cmp_int(fp
, buf
, i
, il1
->iatoms
[i
], il2
->iatoms
[i
]);
295 void cmp_iparm(FILE *fp
, const char *s
, t_functype ft
,
296 t_iparams ip1
, t_iparams ip2
, real ftol
, real abstol
)
302 for (i
= 0; i
< MAXFORCEPARAM
&& !bDiff
; i
++)
304 bDiff
= !equal_real(ip1
.generic
.buf
[i
], ip2
.generic
.buf
[i
], ftol
, abstol
);
308 fprintf(fp
, "%s1: ", s
);
309 pr_iparams(fp
, ft
, &ip1
);
310 fprintf(fp
, "%s2: ", s
);
311 pr_iparams(fp
, ft
, &ip2
);
315 void cmp_iparm_AB(FILE *fp
, const char *s
, t_functype ft
, t_iparams ip1
, real ftol
, real abstol
)
317 int nrfpA
, nrfpB
, p0
, i
;
320 /* Normally the first parameter is perturbable */
322 nrfpA
= interaction_function
[ft
].nrfpA
;
323 nrfpB
= interaction_function
[ft
].nrfpB
;
328 else if (interaction_function
[ft
].flags
& IF_TABULATED
)
330 /* For tabulated interactions only the second parameter is perturbable */
335 for (i
= 0; i
< nrfpB
&& !bDiff
; i
++)
337 bDiff
= !equal_real(ip1
.generic
.buf
[p0
+i
], ip1
.generic
.buf
[nrfpA
+i
], ftol
, abstol
);
341 fprintf(fp
, "%s: ", s
);
342 pr_iparams(fp
, ft
, &ip1
);
346 static void cmp_idef(FILE *fp
, t_idef
*id1
, t_idef
*id2
, real ftol
, real abstol
)
349 char buf1
[64], buf2
[64];
351 fprintf(fp
, "comparing idef\n");
354 cmp_int(fp
, "idef->ntypes", -1, id1
->ntypes
, id2
->ntypes
);
355 cmp_int(fp
, "idef->atnr", -1, id1
->atnr
, id2
->atnr
);
356 for (i
= 0; (i
< std::min(id1
->ntypes
, id2
->ntypes
)); i
++)
358 sprintf(buf1
, "idef->functype[%d]", i
);
359 sprintf(buf2
, "idef->iparam[%d]", i
);
360 cmp_int(fp
, buf1
, i
, (int)id1
->functype
[i
], (int)id2
->functype
[i
]);
361 cmp_iparm(fp
, buf2
, id1
->functype
[i
],
362 id1
->iparams
[i
], id2
->iparams
[i
], ftol
, abstol
);
364 cmp_real(fp
, "fudgeQQ", -1, id1
->fudgeQQ
, id2
->fudgeQQ
, ftol
, abstol
);
365 for (i
= 0; (i
< F_NRE
); i
++)
367 cmp_ilist(fp
, i
, &(id1
->il
[i
]), &(id2
->il
[i
]));
372 for (i
= 0; (i
< id1
->ntypes
); i
++)
374 cmp_iparm_AB(fp
, "idef->iparam", id1
->functype
[i
], id1
->iparams
[i
], ftol
, abstol
);
379 static void cmp_block(FILE *fp
, t_block
*b1
, t_block
*b2
, const char *s
)
383 fprintf(fp
, "comparing block %s\n", s
);
384 sprintf(buf
, "%s.nr", s
);
385 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
388 static void cmp_blocka(FILE *fp
, t_blocka
*b1
, t_blocka
*b2
, const char *s
)
392 fprintf(fp
, "comparing blocka %s\n", s
);
393 sprintf(buf
, "%s.nr", s
);
394 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
395 sprintf(buf
, "%s.nra", s
);
396 cmp_int(fp
, buf
, -1, b1
->nra
, b2
->nra
);
399 static void cmp_atom(FILE *fp
, int index
, t_atom
*a1
, t_atom
*a2
, real ftol
, real abstol
)
403 cmp_us(fp
, "atom.type", index
, a1
->type
, a2
->type
);
404 cmp_us(fp
, "atom.ptype", index
, a1
->ptype
, a2
->ptype
);
405 cmp_int(fp
, "atom.resind", index
, a1
->resind
, a2
->resind
);
406 cmp_int(fp
, "atom.atomnumber", index
, a1
->atomnumber
, a2
->atomnumber
);
407 cmp_real(fp
, "atom.m", index
, a1
->m
, a2
->m
, ftol
, abstol
);
408 cmp_real(fp
, "atom.q", index
, a1
->q
, a2
->q
, ftol
, abstol
);
409 cmp_us(fp
, "atom.typeB", index
, a1
->typeB
, a2
->typeB
);
410 cmp_real(fp
, "atom.mB", index
, a1
->mB
, a2
->mB
, ftol
, abstol
);
411 cmp_real(fp
, "atom.qB", index
, a1
->qB
, a2
->qB
, ftol
, abstol
);
415 cmp_us(fp
, "atom.type", index
, a1
->type
, a1
->typeB
);
416 cmp_real(fp
, "atom.m", index
, a1
->m
, a1
->mB
, ftol
, abstol
);
417 cmp_real(fp
, "atom.q", index
, a1
->q
, a1
->qB
, ftol
, abstol
);
421 static void cmp_atoms(FILE *fp
, t_atoms
*a1
, t_atoms
*a2
, real ftol
, real abstol
)
425 fprintf(fp
, "comparing atoms\n");
429 cmp_int(fp
, "atoms->nr", -1, a1
->nr
, a2
->nr
);
430 for (i
= 0; (i
< a1
->nr
); i
++)
432 cmp_atom(fp
, i
, &(a1
->atom
[i
]), &(a2
->atom
[i
]), ftol
, abstol
);
437 for (i
= 0; (i
< a1
->nr
); i
++)
439 cmp_atom(fp
, i
, &(a1
->atom
[i
]), NULL
, ftol
, abstol
);
444 static void cmp_top(FILE *fp
, t_topology
*t1
, t_topology
*t2
, real ftol
, real abstol
)
446 fprintf(fp
, "comparing top\n");
449 cmp_idef(fp
, &(t1
->idef
), &(t2
->idef
), ftol
, abstol
);
450 cmp_atoms(fp
, &(t1
->atoms
), &(t2
->atoms
), ftol
, abstol
);
451 cmp_block(fp
, &t1
->cgs
, &t2
->cgs
, "cgs");
452 cmp_block(fp
, &t1
->mols
, &t2
->mols
, "mols");
453 cmp_bool(fp
, "bIntermolecularInteractions", -1, t1
->bIntermolecularInteractions
, t2
->bIntermolecularInteractions
);
454 cmp_blocka(fp
, &t1
->excls
, &t2
->excls
, "excls");
458 cmp_idef(fp
, &(t1
->idef
), NULL
, ftol
, abstol
);
459 cmp_atoms(fp
, &(t1
->atoms
), NULL
, ftol
, abstol
);
463 static void cmp_groups(FILE *fp
, gmx_groups_t
*g0
, gmx_groups_t
*g1
,
464 int natoms0
, int natoms1
)
469 fprintf(fp
, "comparing groups\n");
471 for (i
= 0; i
< egcNR
; i
++)
473 sprintf(buf
, "grps[%d].nr", i
);
474 cmp_int(fp
, buf
, -1, g0
->grps
[i
].nr
, g1
->grps
[i
].nr
);
475 if (g0
->grps
[i
].nr
== g1
->grps
[i
].nr
)
477 for (j
= 0; j
< g0
->grps
[i
].nr
; j
++)
479 sprintf(buf
, "grps[%d].name[%d]", i
, j
);
481 *g0
->grpname
[g0
->grps
[i
].nm_ind
[j
]],
482 *g1
->grpname
[g1
->grps
[i
].nm_ind
[j
]]);
485 cmp_int(fp
, "ngrpnr", i
, g0
->ngrpnr
[i
], g1
->ngrpnr
[i
]);
486 if (g0
->ngrpnr
[i
] == g1
->ngrpnr
[i
] && natoms0
== natoms1
&&
487 (g0
->grpnr
[i
] != NULL
|| g1
->grpnr
[i
] != NULL
))
489 for (j
= 0; j
< natoms0
; j
++)
491 cmp_int(fp
, gtypes
[i
], j
, ggrpnr(g0
, i
, j
), ggrpnr(g1
, i
, j
));
495 /* We have compared the names in the groups lists,
496 * so we can skip the grpname list comparison.
500 static void cmp_rvecs(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
501 gmx_bool bRMSD
, real ftol
, real abstol
)
509 for (i
= 0; (i
< n
); i
++)
511 for (m
= 0; m
< DIM
; m
++)
513 d
= x1
[i
][m
] - x2
[i
][m
];
517 fprintf(fp
, "%s RMSD %g\n", title
, std::sqrt(ssd
/n
));
521 for (i
= 0; (i
< n
); i
++)
523 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], ftol
, abstol
);
529 /* Similar to cmp_rvecs, but this routine scales the allowed absolute tolerance
530 * by the RMS of the force components of x1.
532 static void cmp_rvecs_rmstol(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
533 real ftol
, real abstol
)
537 double ave_x1
, rms_x1
;
539 /* It is tricky to compare real values, in particular forces that
540 * are sums of lots of terms where the final value might be close to 0.0.
541 * To get a reference magnitude we calculate the RMS value of each
542 * component in x1, and then set the allowed absolute tolerance to the
543 * relative tolerance times this RMS magnitude.
546 for (i
= 0; i
< n
; i
++)
548 for (m
= 0; m
< DIM
; m
++)
556 for (i
= 0; (i
< n
); i
++)
558 for (m
= 0; m
< DIM
; m
++)
560 d
= x1
[i
][m
] - ave_x1
;
564 rms_x1
= std::sqrt(rms_x1
/(DIM
*n
));
565 /* And now do the actual comparision with a hopefully realistic abstol. */
566 for (i
= 0; (i
< n
); i
++)
568 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], ftol
, abstol
*rms_x1
);
572 static void cmp_grpopts(FILE *fp
, t_grpopts
*opt1
, t_grpopts
*opt2
, real ftol
, real abstol
)
575 char buf1
[256], buf2
[256];
577 cmp_int(fp
, "inputrec->grpopts.ngtc", -1, opt1
->ngtc
, opt2
->ngtc
);
578 cmp_int(fp
, "inputrec->grpopts.ngacc", -1, opt1
->ngacc
, opt2
->ngacc
);
579 cmp_int(fp
, "inputrec->grpopts.ngfrz", -1, opt1
->ngfrz
, opt2
->ngfrz
);
580 cmp_int(fp
, "inputrec->grpopts.ngener", -1, opt1
->ngener
, opt2
->ngener
);
581 for (i
= 0; (i
< std::min(opt1
->ngtc
, opt2
->ngtc
)); i
++)
583 cmp_real(fp
, "inputrec->grpopts.nrdf", i
, opt1
->nrdf
[i
], opt2
->nrdf
[i
], ftol
, abstol
);
584 cmp_real(fp
, "inputrec->grpopts.ref_t", i
, opt1
->ref_t
[i
], opt2
->ref_t
[i
], ftol
, abstol
);
585 cmp_real(fp
, "inputrec->grpopts.tau_t", i
, opt1
->tau_t
[i
], opt2
->tau_t
[i
], ftol
, abstol
);
586 cmp_int(fp
, "inputrec->grpopts.annealing", i
, opt1
->annealing
[i
], opt2
->annealing
[i
]);
587 cmp_int(fp
, "inputrec->grpopts.anneal_npoints", i
,
588 opt1
->anneal_npoints
[i
], opt2
->anneal_npoints
[i
]);
589 if (opt1
->anneal_npoints
[i
] == opt2
->anneal_npoints
[i
])
591 sprintf(buf1
, "inputrec->grpopts.anneal_time[%d]", i
);
592 sprintf(buf2
, "inputrec->grpopts.anneal_temp[%d]", i
);
593 for (j
= 0; j
< opt1
->anneal_npoints
[i
]; j
++)
595 cmp_real(fp
, buf1
, j
, opt1
->anneal_time
[i
][j
], opt2
->anneal_time
[i
][j
], ftol
, abstol
);
596 cmp_real(fp
, buf2
, j
, opt1
->anneal_temp
[i
][j
], opt2
->anneal_temp
[i
][j
], ftol
, abstol
);
600 if (opt1
->ngener
== opt2
->ngener
)
602 for (i
= 0; i
< opt1
->ngener
; i
++)
604 for (j
= i
; j
< opt1
->ngener
; j
++)
606 sprintf(buf1
, "inputrec->grpopts.egp_flags[%d]", i
);
608 opt1
->egp_flags
[opt1
->ngener
*i
+j
],
609 opt2
->egp_flags
[opt1
->ngener
*i
+j
]);
613 for (i
= 0; (i
< std::min(opt1
->ngacc
, opt2
->ngacc
)); i
++)
615 cmp_rvec(fp
, "inputrec->grpopts.acc", i
, opt1
->acc
[i
], opt2
->acc
[i
], ftol
, abstol
);
617 for (i
= 0; (i
< std::min(opt1
->ngfrz
, opt2
->ngfrz
)); i
++)
619 cmp_ivec(fp
, "inputrec->grpopts.nFreeze", i
, opt1
->nFreeze
[i
], opt2
->nFreeze
[i
]);
623 static void cmp_cosines(FILE *fp
, const char *s
, t_cosines c1
[DIM
], t_cosines c2
[DIM
], real ftol
, real abstol
)
628 for (m
= 0; (m
< DIM
); m
++)
630 sprintf(buf
, "inputrec->%s[%d]", s
, m
);
631 cmp_int(fp
, buf
, 0, c1
->n
, c2
->n
);
632 for (i
= 0; (i
< std::min(c1
->n
, c2
->n
)); i
++)
634 cmp_real(fp
, buf
, i
, c1
->a
[i
], c2
->a
[i
], ftol
, abstol
);
635 cmp_real(fp
, buf
, i
, c1
->phi
[i
], c2
->phi
[i
], ftol
, abstol
);
639 static void cmp_pull(FILE *fp
)
641 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");
644 static void cmp_simtempvals(FILE *fp
, t_simtemp
*simtemp1
, t_simtemp
*simtemp2
, int n_lambda
, real ftol
, real abstol
)
647 cmp_int(fp
, "inputrec->simtempvals->eSimTempScale", -1, simtemp1
->eSimTempScale
, simtemp2
->eSimTempScale
);
648 cmp_real(fp
, "inputrec->simtempvals->simtemp_high", -1, simtemp1
->simtemp_high
, simtemp2
->simtemp_high
, ftol
, abstol
);
649 cmp_real(fp
, "inputrec->simtempvals->simtemp_low", -1, simtemp1
->simtemp_low
, simtemp2
->simtemp_low
, ftol
, abstol
);
650 for (i
= 0; i
< n_lambda
; i
++)
652 cmp_real(fp
, "inputrec->simtempvals->temperatures", -1, simtemp1
->temperatures
[i
], simtemp2
->temperatures
[i
], ftol
, abstol
);
656 static void cmp_expandedvals(FILE *fp
, t_expanded
*expand1
, t_expanded
*expand2
, int n_lambda
, real ftol
, real abstol
)
660 cmp_bool(fp
, "inputrec->fepvals->bInit_weights", -1, expand1
->bInit_weights
, expand2
->bInit_weights
);
661 cmp_bool(fp
, "inputrec->fepvals->bWLoneovert", -1, expand1
->bWLoneovert
, expand2
->bWLoneovert
);
663 for (i
= 0; i
< n_lambda
; i
++)
665 cmp_real(fp
, "inputrec->expandedvals->init_lambda_weights", -1,
666 expand1
->init_lambda_weights
[i
], expand2
->init_lambda_weights
[i
], ftol
, abstol
);
669 cmp_int(fp
, "inputrec->expandedvals->lambda-stats", -1, expand1
->elamstats
, expand2
->elamstats
);
670 cmp_int(fp
, "inputrec->expandedvals->lambda-mc-move", -1, expand1
->elmcmove
, expand2
->elmcmove
);
671 cmp_int(fp
, "inputrec->expandedvals->lmc-repeats", -1, expand1
->lmc_repeats
, expand2
->lmc_repeats
);
672 cmp_int(fp
, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1
->gibbsdeltalam
, expand2
->gibbsdeltalam
);
673 cmp_int(fp
, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1
->lmc_forced_nstart
, expand2
->lmc_forced_nstart
);
674 cmp_int(fp
, "inputrec->expandedvals->lambda-weights-equil", -1, expand1
->elmceq
, expand2
->elmceq
);
675 cmp_int(fp
, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1
->equil_n_at_lam
, expand2
->equil_n_at_lam
);
676 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1
->equil_samples
, expand2
->equil_samples
);
677 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1
->equil_steps
, expand2
->equil_steps
);
678 cmp_real(fp
, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1
->equil_wl_delta
, expand2
->equil_wl_delta
, ftol
, abstol
);
679 cmp_real(fp
, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1
->equil_ratio
, expand2
->equil_ratio
, ftol
, abstol
);
680 cmp_bool(fp
, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1
->bSymmetrizedTMatrix
, expand2
->bSymmetrizedTMatrix
);
681 cmp_int(fp
, "inputrec->expandedvals->nstTij", -1, expand1
->nstTij
, expand2
->nstTij
);
682 cmp_int(fp
, "inputrec->expandedvals->mininum-var-min", -1, expand1
->minvarmin
, expand2
->minvarmin
); /*default is reasonable */
683 cmp_int(fp
, "inputrec->expandedvals->weight-c-range", -1, expand1
->c_range
, expand2
->c_range
); /* default is just C=0 */
684 cmp_real(fp
, "inputrec->expandedvals->wl-scale", -1, expand1
->wl_scale
, expand2
->wl_scale
, ftol
, abstol
);
685 cmp_real(fp
, "inputrec->expandedvals->init-wl-delta", -1, expand1
->init_wl_delta
, expand2
->init_wl_delta
, ftol
, abstol
);
686 cmp_real(fp
, "inputrec->expandedvals->wl-ratio", -1, expand1
->wl_ratio
, expand2
->wl_ratio
, ftol
, abstol
);
687 cmp_int(fp
, "inputrec->expandedvals->nstexpanded", -1, expand1
->nstexpanded
, expand2
->nstexpanded
);
688 cmp_int(fp
, "inputrec->expandedvals->lmc-seed", -1, expand1
->lmc_seed
, expand2
->lmc_seed
);
689 cmp_real(fp
, "inputrec->expandedvals->mc-temperature", -1, expand1
->mc_temp
, expand2
->mc_temp
, ftol
, abstol
);
692 static void cmp_fepvals(FILE *fp
, t_lambda
*fep1
, t_lambda
*fep2
, real ftol
, real abstol
)
695 cmp_int(fp
, "inputrec->nstdhdl", -1, fep1
->nstdhdl
, fep2
->nstdhdl
);
696 cmp_double(fp
, "inputrec->fepvals->init_fep_state", -1, fep1
->init_fep_state
, fep2
->init_fep_state
, ftol
, abstol
);
697 cmp_double(fp
, "inputrec->fepvals->delta_lambda", -1, fep1
->delta_lambda
, fep2
->delta_lambda
, ftol
, abstol
);
698 cmp_int(fp
, "inputrec->fepvals->n_lambda", -1, fep1
->n_lambda
, fep2
->n_lambda
);
699 for (i
= 0; i
< efptNR
; i
++)
701 for (j
= 0; j
< std::min(fep1
->n_lambda
, fep2
->n_lambda
); j
++)
703 cmp_double(fp
, "inputrec->fepvals->all_lambda", -1, fep1
->all_lambda
[i
][j
], fep2
->all_lambda
[i
][j
], ftol
, abstol
);
706 cmp_int(fp
, "inputrec->fepvals->lambda_neighbors", 1, fep1
->lambda_neighbors
,
707 fep2
->lambda_neighbors
);
708 cmp_real(fp
, "inputrec->fepvals->sc_alpha", -1, fep1
->sc_alpha
, fep2
->sc_alpha
, ftol
, abstol
);
709 cmp_int(fp
, "inputrec->fepvals->sc_power", -1, fep1
->sc_power
, fep2
->sc_power
);
710 cmp_real(fp
, "inputrec->fepvals->sc_r_power", -1, fep1
->sc_r_power
, fep2
->sc_r_power
, ftol
, abstol
);
711 cmp_real(fp
, "inputrec->fepvals->sc_sigma", -1, fep1
->sc_sigma
, fep2
->sc_sigma
, ftol
, abstol
);
712 cmp_int(fp
, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1
->edHdLPrintEnergy
, fep1
->edHdLPrintEnergy
);
713 cmp_bool(fp
, "inputrec->fepvals->bScCoul", -1, fep1
->bScCoul
, fep1
->bScCoul
);
714 cmp_int(fp
, "inputrec->separate_dhdl_file", -1, fep1
->separate_dhdl_file
, fep2
->separate_dhdl_file
);
715 cmp_int(fp
, "inputrec->dhdl_derivatives", -1, fep1
->dhdl_derivatives
, fep2
->dhdl_derivatives
);
716 cmp_int(fp
, "inputrec->dh_hist_size", -1, fep1
->dh_hist_size
, fep2
->dh_hist_size
);
717 cmp_double(fp
, "inputrec->dh_hist_spacing", -1, fep1
->dh_hist_spacing
, fep2
->dh_hist_spacing
, ftol
, abstol
);
720 static void cmp_inputrec(FILE *fp
, t_inputrec
*ir1
, t_inputrec
*ir2
, real ftol
, real abstol
)
722 fprintf(fp
, "comparing inputrec\n");
724 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
725 * of warnings. Maybe it will change in future versions, but for the
726 * moment I've spelled them out instead. /EL 000820
727 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
728 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
729 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
731 cmp_int(fp
, "inputrec->eI", -1, ir1
->eI
, ir2
->eI
);
732 cmp_int64(fp
, "inputrec->nsteps", ir1
->nsteps
, ir2
->nsteps
);
733 cmp_int64(fp
, "inputrec->init_step", ir1
->init_step
, ir2
->init_step
);
734 cmp_int(fp
, "inputrec->simulation_part", -1, ir1
->simulation_part
, ir2
->simulation_part
);
735 cmp_int(fp
, "inputrec->ePBC", -1, ir1
->ePBC
, ir2
->ePBC
);
736 cmp_int(fp
, "inputrec->bPeriodicMols", -1, ir1
->bPeriodicMols
, ir2
->bPeriodicMols
);
737 cmp_int(fp
, "inputrec->cutoff_scheme", -1, ir1
->cutoff_scheme
, ir2
->cutoff_scheme
);
738 cmp_int(fp
, "inputrec->ns_type", -1, ir1
->ns_type
, ir2
->ns_type
);
739 cmp_int(fp
, "inputrec->nstlist", -1, ir1
->nstlist
, ir2
->nstlist
);
740 cmp_int(fp
, "inputrec->nstcomm", -1, ir1
->nstcomm
, ir2
->nstcomm
);
741 cmp_int(fp
, "inputrec->comm_mode", -1, ir1
->comm_mode
, ir2
->comm_mode
);
742 cmp_int(fp
, "inputrec->nstlog", -1, ir1
->nstlog
, ir2
->nstlog
);
743 cmp_int(fp
, "inputrec->nstxout", -1, ir1
->nstxout
, ir2
->nstxout
);
744 cmp_int(fp
, "inputrec->nstvout", -1, ir1
->nstvout
, ir2
->nstvout
);
745 cmp_int(fp
, "inputrec->nstfout", -1, ir1
->nstfout
, ir2
->nstfout
);
746 cmp_int(fp
, "inputrec->nstcalcenergy", -1, ir1
->nstcalcenergy
, ir2
->nstcalcenergy
);
747 cmp_int(fp
, "inputrec->nstenergy", -1, ir1
->nstenergy
, ir2
->nstenergy
);
748 cmp_int(fp
, "inputrec->nstxout_compressed", -1, ir1
->nstxout_compressed
, ir2
->nstxout_compressed
);
749 cmp_double(fp
, "inputrec->init_t", -1, ir1
->init_t
, ir2
->init_t
, ftol
, abstol
);
750 cmp_double(fp
, "inputrec->delta_t", -1, ir1
->delta_t
, ir2
->delta_t
, ftol
, abstol
);
751 cmp_real(fp
, "inputrec->x_compression_precision", -1, ir1
->x_compression_precision
, ir2
->x_compression_precision
, ftol
, abstol
);
752 cmp_real(fp
, "inputrec->fourierspacing", -1, ir1
->fourier_spacing
, ir2
->fourier_spacing
, ftol
, abstol
);
753 cmp_int(fp
, "inputrec->nkx", -1, ir1
->nkx
, ir2
->nkx
);
754 cmp_int(fp
, "inputrec->nky", -1, ir1
->nky
, ir2
->nky
);
755 cmp_int(fp
, "inputrec->nkz", -1, ir1
->nkz
, ir2
->nkz
);
756 cmp_int(fp
, "inputrec->pme_order", -1, ir1
->pme_order
, ir2
->pme_order
);
757 cmp_real(fp
, "inputrec->ewald_rtol", -1, ir1
->ewald_rtol
, ir2
->ewald_rtol
, ftol
, abstol
);
758 cmp_int(fp
, "inputrec->ewald_geometry", -1, ir1
->ewald_geometry
, ir2
->ewald_geometry
);
759 cmp_real(fp
, "inputrec->epsilon_surface", -1, ir1
->epsilon_surface
, ir2
->epsilon_surface
, ftol
, abstol
);
760 cmp_int(fp
, "inputrec->bContinuation", -1, ir1
->bContinuation
, ir2
->bContinuation
);
761 cmp_int(fp
, "inputrec->bShakeSOR", -1, ir1
->bShakeSOR
, ir2
->bShakeSOR
);
762 cmp_int(fp
, "inputrec->etc", -1, ir1
->etc
, ir2
->etc
);
763 cmp_int(fp
, "inputrec->bPrintNHChains", -1, ir1
->bPrintNHChains
, ir2
->bPrintNHChains
);
764 cmp_int(fp
, "inputrec->epc", -1, ir1
->epc
, ir2
->epc
);
765 cmp_int(fp
, "inputrec->epct", -1, ir1
->epct
, ir2
->epct
);
766 cmp_real(fp
, "inputrec->tau_p", -1, ir1
->tau_p
, ir2
->tau_p
, ftol
, abstol
);
767 cmp_rvec(fp
, "inputrec->ref_p(x)", -1, ir1
->ref_p
[XX
], ir2
->ref_p
[XX
], ftol
, abstol
);
768 cmp_rvec(fp
, "inputrec->ref_p(y)", -1, ir1
->ref_p
[YY
], ir2
->ref_p
[YY
], ftol
, abstol
);
769 cmp_rvec(fp
, "inputrec->ref_p(z)", -1, ir1
->ref_p
[ZZ
], ir2
->ref_p
[ZZ
], ftol
, abstol
);
770 cmp_rvec(fp
, "inputrec->compress(x)", -1, ir1
->compress
[XX
], ir2
->compress
[XX
], ftol
, abstol
);
771 cmp_rvec(fp
, "inputrec->compress(y)", -1, ir1
->compress
[YY
], ir2
->compress
[YY
], ftol
, abstol
);
772 cmp_rvec(fp
, "inputrec->compress(z)", -1, ir1
->compress
[ZZ
], ir2
->compress
[ZZ
], ftol
, abstol
);
773 cmp_int(fp
, "refcoord_scaling", -1, ir1
->refcoord_scaling
, ir2
->refcoord_scaling
);
774 cmp_rvec(fp
, "inputrec->posres_com", -1, ir1
->posres_com
, ir2
->posres_com
, ftol
, abstol
);
775 cmp_rvec(fp
, "inputrec->posres_comB", -1, ir1
->posres_comB
, ir2
->posres_comB
, ftol
, abstol
);
776 cmp_real(fp
, "inputrec->verletbuf_tol", -1, ir1
->verletbuf_tol
, ir2
->verletbuf_tol
, ftol
, abstol
);
777 cmp_real(fp
, "inputrec->rlist", -1, ir1
->rlist
, ir2
->rlist
, ftol
, abstol
);
778 cmp_real(fp
, "inputrec->rtpi", -1, ir1
->rtpi
, ir2
->rtpi
, ftol
, abstol
);
779 cmp_int(fp
, "inputrec->coulombtype", -1, ir1
->coulombtype
, ir2
->coulombtype
);
780 cmp_int(fp
, "inputrec->coulomb_modifier", -1, ir1
->coulomb_modifier
, ir2
->coulomb_modifier
);
781 cmp_real(fp
, "inputrec->rcoulomb_switch", -1, ir1
->rcoulomb_switch
, ir2
->rcoulomb_switch
, ftol
, abstol
);
782 cmp_real(fp
, "inputrec->rcoulomb", -1, ir1
->rcoulomb
, ir2
->rcoulomb
, ftol
, abstol
);
783 cmp_int(fp
, "inputrec->vdwtype", -1, ir1
->vdwtype
, ir2
->vdwtype
);
784 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
);
785 cmp_real(fp
, "inputrec->rvdw", -1, ir1
->rvdw
, ir2
->rvdw
, ftol
, abstol
);
786 cmp_real(fp
, "inputrec->epsilon_r", -1, ir1
->epsilon_r
, ir2
->epsilon_r
, ftol
, abstol
);
787 cmp_real(fp
, "inputrec->epsilon_rf", -1, ir1
->epsilon_rf
, ir2
->epsilon_rf
, ftol
, abstol
);
788 cmp_real(fp
, "inputrec->tabext", -1, ir1
->tabext
, ir2
->tabext
, ftol
, abstol
);
789 cmp_int(fp
, "inputrec->implicit_solvent", -1, ir1
->implicit_solvent
, ir2
->implicit_solvent
);
790 cmp_int(fp
, "inputrec->gb_algorithm", -1, ir1
->gb_algorithm
, ir2
->gb_algorithm
);
791 cmp_int(fp
, "inputrec->nstgbradii", -1, ir1
->nstgbradii
, ir2
->nstgbradii
);
792 cmp_real(fp
, "inputrec->rgbradii", -1, ir1
->rgbradii
, ir2
->rgbradii
, ftol
, abstol
);
793 cmp_real(fp
, "inputrec->gb_saltconc", -1, ir1
->gb_saltconc
, ir2
->gb_saltconc
, ftol
, abstol
);
794 cmp_real(fp
, "inputrec->gb_epsilon_solvent", -1, ir1
->gb_epsilon_solvent
, ir2
->gb_epsilon_solvent
, ftol
, abstol
);
795 cmp_real(fp
, "inputrec->gb_obc_alpha", -1, ir1
->gb_obc_alpha
, ir2
->gb_obc_alpha
, ftol
, abstol
);
796 cmp_real(fp
, "inputrec->gb_obc_beta", -1, ir1
->gb_obc_beta
, ir2
->gb_obc_beta
, ftol
, abstol
);
797 cmp_real(fp
, "inputrec->gb_obc_gamma", -1, ir1
->gb_obc_gamma
, ir2
->gb_obc_gamma
, ftol
, abstol
);
798 cmp_real(fp
, "inputrec->gb_dielectric_offset", -1, ir1
->gb_dielectric_offset
, ir2
->gb_dielectric_offset
, ftol
, abstol
);
799 cmp_int(fp
, "inputrec->sa_algorithm", -1, ir1
->sa_algorithm
, ir2
->sa_algorithm
);
800 cmp_real(fp
, "inputrec->sa_surface_tension", -1, ir1
->sa_surface_tension
, ir2
->sa_surface_tension
, ftol
, abstol
);
802 cmp_int(fp
, "inputrec->eDispCorr", -1, ir1
->eDispCorr
, ir2
->eDispCorr
);
803 cmp_real(fp
, "inputrec->shake_tol", -1, ir1
->shake_tol
, ir2
->shake_tol
, ftol
, abstol
);
804 cmp_int(fp
, "inputrec->efep", -1, ir1
->efep
, ir2
->efep
);
805 cmp_fepvals(fp
, ir1
->fepvals
, ir2
->fepvals
, ftol
, abstol
);
806 cmp_int(fp
, "inputrec->bSimTemp", -1, ir1
->bSimTemp
, ir2
->bSimTemp
);
807 if ((ir1
->bSimTemp
== ir2
->bSimTemp
) && (ir1
->bSimTemp
))
809 cmp_simtempvals(fp
, ir1
->simtempvals
, ir2
->simtempvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
811 cmp_int(fp
, "inputrec->bExpanded", -1, ir1
->bExpanded
, ir2
->bExpanded
);
812 if ((ir1
->bExpanded
== ir2
->bExpanded
) && (ir1
->bExpanded
))
814 cmp_expandedvals(fp
, ir1
->expandedvals
, ir2
->expandedvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
816 cmp_int(fp
, "inputrec->nwall", -1, ir1
->nwall
, ir2
->nwall
);
817 cmp_int(fp
, "inputrec->wall_type", -1, ir1
->wall_type
, ir2
->wall_type
);
818 cmp_int(fp
, "inputrec->wall_atomtype[0]", -1, ir1
->wall_atomtype
[0], ir2
->wall_atomtype
[0]);
819 cmp_int(fp
, "inputrec->wall_atomtype[1]", -1, ir1
->wall_atomtype
[1], ir2
->wall_atomtype
[1]);
820 cmp_real(fp
, "inputrec->wall_density[0]", -1, ir1
->wall_density
[0], ir2
->wall_density
[0], ftol
, abstol
);
821 cmp_real(fp
, "inputrec->wall_density[1]", -1, ir1
->wall_density
[1], ir2
->wall_density
[1], ftol
, abstol
);
822 cmp_real(fp
, "inputrec->wall_ewald_zfac", -1, ir1
->wall_ewald_zfac
, ir2
->wall_ewald_zfac
, ftol
, abstol
);
824 cmp_bool(fp
, "inputrec->bPull", -1, ir1
->bPull
, ir2
->bPull
);
825 if (ir1
->bPull
&& ir2
->bPull
)
830 cmp_int(fp
, "inputrec->eDisre", -1, ir1
->eDisre
, ir2
->eDisre
);
831 cmp_real(fp
, "inputrec->dr_fc", -1, ir1
->dr_fc
, ir2
->dr_fc
, ftol
, abstol
);
832 cmp_int(fp
, "inputrec->eDisreWeighting", -1, ir1
->eDisreWeighting
, ir2
->eDisreWeighting
);
833 cmp_int(fp
, "inputrec->bDisreMixed", -1, ir1
->bDisreMixed
, ir2
->bDisreMixed
);
834 cmp_int(fp
, "inputrec->nstdisreout", -1, ir1
->nstdisreout
, ir2
->nstdisreout
);
835 cmp_real(fp
, "inputrec->dr_tau", -1, ir1
->dr_tau
, ir2
->dr_tau
, ftol
, abstol
);
836 cmp_real(fp
, "inputrec->orires_fc", -1, ir1
->orires_fc
, ir2
->orires_fc
, ftol
, abstol
);
837 cmp_real(fp
, "inputrec->orires_tau", -1, ir1
->orires_tau
, ir2
->orires_tau
, ftol
, abstol
);
838 cmp_int(fp
, "inputrec->nstorireout", -1, ir1
->nstorireout
, ir2
->nstorireout
);
839 cmp_real(fp
, "inputrec->em_stepsize", -1, ir1
->em_stepsize
, ir2
->em_stepsize
, ftol
, abstol
);
840 cmp_real(fp
, "inputrec->em_tol", -1, ir1
->em_tol
, ir2
->em_tol
, ftol
, abstol
);
841 cmp_int(fp
, "inputrec->niter", -1, ir1
->niter
, ir2
->niter
);
842 cmp_real(fp
, "inputrec->fc_stepsize", -1, ir1
->fc_stepsize
, ir2
->fc_stepsize
, ftol
, abstol
);
843 cmp_int(fp
, "inputrec->nstcgsteep", -1, ir1
->nstcgsteep
, ir2
->nstcgsteep
);
844 cmp_int(fp
, "inputrec->nbfgscorr", 0, ir1
->nbfgscorr
, ir2
->nbfgscorr
);
845 cmp_int(fp
, "inputrec->eConstrAlg", -1, ir1
->eConstrAlg
, ir2
->eConstrAlg
);
846 cmp_int(fp
, "inputrec->nProjOrder", -1, ir1
->nProjOrder
, ir2
->nProjOrder
);
847 cmp_real(fp
, "inputrec->LincsWarnAngle", -1, ir1
->LincsWarnAngle
, ir2
->LincsWarnAngle
, ftol
, abstol
);
848 cmp_int(fp
, "inputrec->nLincsIter", -1, ir1
->nLincsIter
, ir2
->nLincsIter
);
849 cmp_real(fp
, "inputrec->bd_fric", -1, ir1
->bd_fric
, ir2
->bd_fric
, ftol
, abstol
);
850 cmp_int64(fp
, "inputrec->ld_seed", ir1
->ld_seed
, ir2
->ld_seed
);
851 cmp_real(fp
, "inputrec->cos_accel", -1, ir1
->cos_accel
, ir2
->cos_accel
, ftol
, abstol
);
852 cmp_rvec(fp
, "inputrec->deform(a)", -1, ir1
->deform
[XX
], ir2
->deform
[XX
], ftol
, abstol
);
853 cmp_rvec(fp
, "inputrec->deform(b)", -1, ir1
->deform
[YY
], ir2
->deform
[YY
], ftol
, abstol
);
854 cmp_rvec(fp
, "inputrec->deform(c)", -1, ir1
->deform
[ZZ
], ir2
->deform
[ZZ
], ftol
, abstol
);
857 cmp_int(fp
, "inputrec->userint1", -1, ir1
->userint1
, ir2
->userint1
);
858 cmp_int(fp
, "inputrec->userint2", -1, ir1
->userint2
, ir2
->userint2
);
859 cmp_int(fp
, "inputrec->userint3", -1, ir1
->userint3
, ir2
->userint3
);
860 cmp_int(fp
, "inputrec->userint4", -1, ir1
->userint4
, ir2
->userint4
);
861 cmp_real(fp
, "inputrec->userreal1", -1, ir1
->userreal1
, ir2
->userreal1
, ftol
, abstol
);
862 cmp_real(fp
, "inputrec->userreal2", -1, ir1
->userreal2
, ir2
->userreal2
, ftol
, abstol
);
863 cmp_real(fp
, "inputrec->userreal3", -1, ir1
->userreal3
, ir2
->userreal3
, ftol
, abstol
);
864 cmp_real(fp
, "inputrec->userreal4", -1, ir1
->userreal4
, ir2
->userreal4
, ftol
, abstol
);
865 cmp_grpopts(fp
, &(ir1
->opts
), &(ir2
->opts
), ftol
, abstol
);
866 cmp_cosines(fp
, "ex", ir1
->ex
, ir2
->ex
, ftol
, abstol
);
867 cmp_cosines(fp
, "et", ir1
->et
, ir2
->et
, ftol
, abstol
);
870 static void comp_pull_AB(FILE *fp
, pull_params_t
*pull
, real ftol
, real abstol
)
874 for (i
= 0; i
< pull
->ncoord
; i
++)
876 fprintf(fp
, "comparing pull coord %d\n", i
);
877 cmp_real(fp
, "pull-coord->k", -1, pull
->coord
[i
].k
, pull
->coord
[i
].kB
, ftol
, abstol
);
881 static void comp_state(t_state
*st1
, t_state
*st2
,
882 gmx_bool bRMSD
, real ftol
, real abstol
)
886 fprintf(stdout
, "comparing flags\n");
887 cmp_int(stdout
, "flags", -1, st1
->flags
, st2
->flags
);
888 fprintf(stdout
, "comparing box\n");
889 cmp_rvecs(stdout
, "box", DIM
, st1
->box
, st2
->box
, FALSE
, ftol
, abstol
);
890 fprintf(stdout
, "comparing box_rel\n");
891 cmp_rvecs(stdout
, "box_rel", DIM
, st1
->box_rel
, st2
->box_rel
, FALSE
, ftol
, abstol
);
892 fprintf(stdout
, "comparing boxv\n");
893 cmp_rvecs(stdout
, "boxv", DIM
, st1
->boxv
, st2
->boxv
, FALSE
, ftol
, abstol
);
894 if (st1
->flags
& (1<<estSVIR_PREV
))
896 fprintf(stdout
, "comparing shake vir_prev\n");
897 cmp_rvecs_rmstol(stdout
, "svir_prev", DIM
, st1
->svir_prev
, st2
->svir_prev
, ftol
, abstol
);
899 if (st1
->flags
& (1<<estFVIR_PREV
))
901 fprintf(stdout
, "comparing force vir_prev\n");
902 cmp_rvecs_rmstol(stdout
, "fvir_prev", DIM
, st1
->fvir_prev
, st2
->fvir_prev
, ftol
, abstol
);
904 if (st1
->flags
& (1<<estPRES_PREV
))
906 fprintf(stdout
, "comparing prev_pres\n");
907 cmp_rvecs_rmstol(stdout
, "pres_prev", DIM
, st1
->pres_prev
, st2
->pres_prev
, ftol
, abstol
);
909 cmp_int(stdout
, "ngtc", -1, st1
->ngtc
, st2
->ngtc
);
910 cmp_int(stdout
, "nhchainlength", -1, st1
->nhchainlength
, st2
->nhchainlength
);
911 if (st1
->ngtc
== st2
->ngtc
&& st1
->nhchainlength
== st2
->nhchainlength
)
913 for (i
= 0; i
< st1
->ngtc
; i
++)
915 nc
= i
*st1
->nhchainlength
;
916 for (j
= 0; j
< nc
; j
++)
918 cmp_real(stdout
, "nosehoover_xi",
919 i
, st1
->nosehoover_xi
[nc
+j
], st2
->nosehoover_xi
[nc
+j
], ftol
, abstol
);
923 cmp_int(stdout
, "nnhpres", -1, st1
->nnhpres
, st2
->nnhpres
);
924 if (st1
->nnhpres
== st2
->nnhpres
&& st1
->nhchainlength
== st2
->nhchainlength
)
926 for (i
= 0; i
< st1
->nnhpres
; i
++)
928 nc
= i
*st1
->nhchainlength
;
929 for (j
= 0; j
< nc
; j
++)
931 cmp_real(stdout
, "nosehoover_xi",
932 i
, st1
->nhpres_xi
[nc
+j
], st2
->nhpres_xi
[nc
+j
], ftol
, abstol
);
937 cmp_int(stdout
, "natoms", -1, st1
->natoms
, st2
->natoms
);
938 if (st1
->natoms
== st2
->natoms
)
940 if ((st1
->flags
& (1<<estX
)) && (st2
->flags
& (1<<estX
)))
942 fprintf(stdout
, "comparing x\n");
943 cmp_rvecs(stdout
, "x", st1
->natoms
, st1
->x
, st2
->x
, bRMSD
, ftol
, abstol
);
945 if ((st1
->flags
& (1<<estV
)) && (st2
->flags
& (1<<estV
)))
947 fprintf(stdout
, "comparing v\n");
948 cmp_rvecs(stdout
, "v", st1
->natoms
, st1
->v
, st2
->v
, bRMSD
, ftol
, abstol
);
953 void comp_tpx(const char *fn1
, const char *fn2
,
954 gmx_bool bRMSD
, real ftol
, real abstol
)
965 for (i
= 0; i
< (fn2
? 2 : 1); i
++)
967 read_tpx_state(ff
[i
], &(ir
[i
]), &state
[i
], &(mtop
[i
]));
971 cmp_inputrec(stdout
, &ir
[0], &ir
[1], ftol
, abstol
);
972 /* Convert gmx_mtop_t to t_topology.
973 * We should implement direct mtop comparison,
974 * but it might be useful to keep t_topology comparison as an option.
976 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
977 top
[1] = gmx_mtop_t_to_t_topology(&mtop
[1]);
978 cmp_top(stdout
, &top
[0], &top
[1], ftol
, abstol
);
979 cmp_groups(stdout
, &mtop
[0].groups
, &mtop
[1].groups
,
980 mtop
[0].natoms
, mtop
[1].natoms
);
981 comp_state(&state
[0], &state
[1], bRMSD
, ftol
, abstol
);
985 if (ir
[0].efep
== efepNO
)
987 fprintf(stdout
, "inputrec->efep = %s\n", efep_names
[ir
[0].efep
]);
993 comp_pull_AB(stdout
, ir
->pull
, ftol
, abstol
);
995 /* Convert gmx_mtop_t to t_topology.
996 * We should implement direct mtop comparison,
997 * but it might be useful to keep t_topology comparison as an option.
999 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
1000 cmp_top(stdout
, &top
[0], NULL
, ftol
, abstol
);
1005 void comp_frame(FILE *fp
, t_trxframe
*fr1
, t_trxframe
*fr2
,
1006 gmx_bool bRMSD
, real ftol
, real abstol
)
1009 cmp_int(fp
, "flags", -1, fr1
->flags
, fr2
->flags
);
1010 cmp_int(fp
, "not_ok", -1, fr1
->not_ok
, fr2
->not_ok
);
1011 cmp_int(fp
, "natoms", -1, fr1
->natoms
, fr2
->natoms
);
1012 cmp_real(fp
, "t0", -1, fr1
->t0
, fr2
->t0
, ftol
, abstol
);
1013 if (cmp_bool(fp
, "bTitle", -1, fr1
->bTitle
, fr2
->bTitle
))
1015 cmp_str(fp
, "title", -1, fr1
->title
, fr2
->title
);
1017 if (cmp_bool(fp
, "bStep", -1, fr1
->bStep
, fr2
->bStep
))
1019 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1021 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1022 if (cmp_bool(fp
, "bTime", -1, fr1
->bTime
, fr2
->bTime
))
1024 cmp_real(fp
, "time", -1, fr1
->time
, fr2
->time
, ftol
, abstol
);
1026 if (cmp_bool(fp
, "bLambda", -1, fr1
->bLambda
, fr2
->bLambda
))
1028 cmp_real(fp
, "lambda", -1, fr1
->lambda
, fr2
->lambda
, ftol
, abstol
);
1030 if (cmp_bool(fp
, "bAtoms", -1, fr1
->bAtoms
, fr2
->bAtoms
))
1032 cmp_atoms(fp
, fr1
->atoms
, fr2
->atoms
, ftol
, abstol
);
1034 if (cmp_bool(fp
, "bPrec", -1, fr1
->bPrec
, fr2
->bPrec
))
1036 cmp_real(fp
, "prec", -1, fr1
->prec
, fr2
->prec
, ftol
, abstol
);
1038 if (cmp_bool(fp
, "bX", -1, fr1
->bX
, fr2
->bX
))
1040 cmp_rvecs(fp
, "x", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->x
, fr2
->x
, bRMSD
, ftol
, abstol
);
1042 if (cmp_bool(fp
, "bV", -1, fr1
->bV
, fr2
->bV
))
1044 cmp_rvecs(fp
, "v", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->v
, fr2
->v
, bRMSD
, ftol
, abstol
);
1046 if (cmp_bool(fp
, "bF", -1, fr1
->bF
, fr2
->bF
))
1050 cmp_rvecs(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, bRMSD
, ftol
, abstol
);
1054 cmp_rvecs_rmstol(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, ftol
, abstol
);
1057 if (cmp_bool(fp
, "bBox", -1, fr1
->bBox
, fr2
->bBox
))
1059 cmp_rvecs(fp
, "box", 3, fr1
->box
, fr2
->box
, FALSE
, ftol
, abstol
);
1063 void comp_trx(const gmx_output_env_t
*oenv
, const char *fn1
, const char *fn2
,
1064 gmx_bool bRMSD
, real ftol
, real abstol
)
1069 t_trxstatus
*status
[2];
1074 fprintf(stderr
, "Comparing trajectory files %s and %s\n", fn1
, fn2
);
1075 for (i
= 0; i
< 2; i
++)
1077 b
[i
] = read_first_frame(oenv
, &status
[i
], fn
[i
], &fr
[i
], TRX_READ_X
|TRX_READ_V
|TRX_READ_F
);
1084 comp_frame(stdout
, &(fr
[0]), &(fr
[1]), bRMSD
, ftol
, abstol
);
1086 for (i
= 0; i
< 2; i
++)
1088 b
[i
] = read_next_frame(oenv
, status
[i
], &fr
[i
]);
1091 while (b
[0] && b
[1]);
1093 for (i
= 0; i
< 2; i
++)
1095 if (b
[i
] && !b
[1-i
])
1097 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn
[1-i
], fn
[i
]);
1099 close_trj(status
[i
]);
1104 fprintf(stdout
, "\nBoth files read correctly\n");
1108 static real
ener_tensor_diag(int n
, int *ind1
, int *ind2
,
1111 t_energy e1
[], t_energy e2
[])
1120 d2
= tensi
[i
] - d1
*DIM
;
1122 /* Find the diagonal elements d1 and d2 */
1123 len
= std::strlen(enm1
[ind1
[i
]].name
);
1127 for (j
= 0; j
< n
; j
++)
1129 if (tensi
[j
] >= 0 &&
1130 std::strlen(enm1
[ind1
[j
]].name
) == len
&&
1131 std::strncmp(enm1
[ind1
[i
]].name
, enm1
[ind1
[j
]].name
, len
-2) == 0 &&
1132 (tensi
[j
] == d1
*DIM
+d1
|| tensi
[j
] == d2
*DIM
+d2
))
1134 prod1
*= fabs(e1
[ind1
[j
]].e
);
1135 prod2
*= fabs(e2
[ind2
[j
]].e
);
1142 return 0.5*(std::sqrt(prod1
) + std::sqrt(prod2
));
1150 static gmx_bool
enernm_equal(const char *nm1
, const char *nm2
)
1154 len1
= std::strlen(nm1
);
1155 len2
= std::strlen(nm2
);
1157 /* Remove " (bar)" at the end of a name */
1158 if (len1
> 6 && std::strcmp(nm1
+len1
-6, " (bar)") == 0)
1162 if (len2
> 6 && std::strcmp(nm2
+len2
-6, " (bar)") == 0)
1167 return (len1
== len2
&& gmx_strncasecmp(nm1
, nm2
, len1
) == 0);
1170 static void cmp_energies(FILE *fp
, int step1
, int step2
,
1171 t_energy e1
[], t_energy e2
[],
1173 real ftol
, real abstol
,
1174 int nre
, int *ind1
, int *ind2
, int maxener
)
1177 int *tensi
, len
, d1
, d2
;
1178 real ftol_i
, abstol_i
;
1180 snew(tensi
, maxener
);
1181 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1182 for (i
= 0; (i
< maxener
); i
++)
1186 len
= std::strlen(enm1
[ii
].name
);
1187 if (len
> 3 && enm1
[ii
].name
[len
-3] == '-')
1189 d1
= enm1
[ii
].name
[len
-2] - 'X';
1190 d2
= enm1
[ii
].name
[len
-1] - 'X';
1191 if (d1
>= 0 && d1
< DIM
&&
1192 d2
>= 0 && d2
< DIM
)
1194 tensi
[i
] = d1
*DIM
+ d2
;
1199 for (i
= 0; (i
< maxener
); i
++)
1201 /* Check if this is an off-diagonal tensor element */
1202 if (tensi
[i
] >= 0 && tensi
[i
] != 0 && tensi
[i
] != 4 && tensi
[i
] != 8)
1204 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1206 /* Do the relative tolerance through an absolute tolerance times
1207 * the size of diagonal components of the tensor.
1209 abstol_i
= ftol
*ener_tensor_diag(nre
, ind1
, ind2
, enm1
, tensi
, i
, e1
, e2
);
1212 fprintf(debug
, "tensor '%s' val %f diag %f\n",
1213 enm1
[i
].name
, e1
[i
].e
, abstol_i
/ftol
);
1217 /* We found a diagonal, we need to check with the minimum tolerance */
1218 abstol_i
= std::min(abstol_i
, abstol
);
1222 /* We did not find a diagonal, ignore the relative tolerance check */
1231 if (!equal_real(e1
[ind1
[i
]].e
, e2
[ind2
[i
]].e
, ftol_i
, abstol_i
))
1233 fprintf(fp
, "%-15s step %3d: %12g, step %3d: %12g\n",
1235 step1
, e1
[ind1
[i
]].e
,
1236 step2
, e2
[ind2
[i
]].e
);
1244 static void cmp_disres(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1247 char bav
[64], bt
[64], bs
[22];
1249 cmp_int(stdout
, "ndisre", -1, fr1
->ndisre
, fr2
->ndisre
);
1250 if ((fr1
->ndisre
== fr2
->ndisre
) && (fr1
->ndisre
> 0))
1252 sprintf(bav
, "step %s: disre rav", gmx_step_str(fr1
->step
, bs
));
1253 sprintf(bt
, "step %s: disre rt", gmx_step_str(fr1
->step
, bs
));
1254 for (i
= 0; (i
< fr1
->ndisre
); i
++)
1256 cmp_real(stdout
, bav
, i
, fr1
->disre_rm3tav
[i
], fr2
->disre_rm3tav
[i
], ftol
, abstol
);
1257 cmp_real(stdout
, bt
, i
, fr1
->disre_rt
[i
], fr2
->disre_rt
[i
], ftol
, abstol
);
1263 static void cmp_eblocks(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1266 char buf
[64], bs
[22];
1268 cmp_int(stdout
, "nblock", -1, fr1
->nblock
, fr2
->nblock
);
1269 if ((fr1
->nblock
== fr2
->nblock
) && (fr1
->nblock
> 0))
1271 for (j
= 0; (j
< fr1
->nblock
); j
++)
1273 t_enxblock
*b1
, *b2
; /* convenience vars */
1275 b1
= &(fr1
->block
[j
]);
1276 b2
= &(fr2
->block
[j
]);
1278 sprintf(buf
, "step %s: block[%d]", gmx_step_str(fr1
->step
, bs
), j
);
1279 cmp_int(stdout
, buf
, -1, b1
->nsub
, b2
->nsub
);
1280 cmp_int(stdout
, buf
, -1, b1
->id
, b2
->id
);
1282 if ( (b1
->nsub
== b2
->nsub
) && (b1
->id
== b2
->id
) )
1284 for (i
= 0; i
< b1
->nsub
; i
++)
1286 t_enxsubblock
*s1
, *s2
;
1291 cmp_int(stdout
, buf
, -1, (int)s1
->type
, (int)s2
->type
);
1292 cmp_int64(stdout
, buf
, s1
->nr
, s2
->nr
);
1294 if ((s1
->type
== s2
->type
) && (s1
->nr
== s2
->nr
))
1298 case xdr_datatype_float
:
1299 for (k
= 0; k
< s1
->nr
; k
++)
1301 cmp_float(stdout
, buf
, i
,
1302 s1
->fval
[k
], s2
->fval
[k
],
1306 case xdr_datatype_double
:
1307 for (k
= 0; k
< s1
->nr
; k
++)
1309 cmp_double(stdout
, buf
, i
,
1310 s1
->dval
[k
], s2
->dval
[k
],
1314 case xdr_datatype_int
:
1315 for (k
= 0; k
< s1
->nr
; k
++)
1317 cmp_int(stdout
, buf
, i
,
1318 s1
->ival
[k
], s2
->ival
[k
]);
1321 case xdr_datatype_int64
:
1322 for (k
= 0; k
< s1
->nr
; k
++)
1324 cmp_int64(stdout
, buf
,
1325 s1
->lval
[k
], s2
->lval
[k
]);
1328 case xdr_datatype_char
:
1329 for (k
= 0; k
< s1
->nr
; k
++)
1331 cmp_uc(stdout
, buf
, i
,
1332 s1
->cval
[k
], s2
->cval
[k
]);
1335 case xdr_datatype_string
:
1336 for (k
= 0; k
< s1
->nr
; k
++)
1338 cmp_str(stdout
, buf
, i
,
1339 s1
->sval
[k
], s2
->sval
[k
]);
1343 gmx_incons("Unknown data type!!");
1352 void comp_enx(const char *fn1
, const char *fn2
, real ftol
, real abstol
, const char *lastener
)
1354 int nre
, nre1
, nre2
;
1355 ener_file_t in1
, in2
;
1356 int i
, j
, maxener
, *ind1
, *ind2
, *have
;
1357 gmx_enxnm_t
*enm1
= NULL
, *enm2
= NULL
;
1358 t_enxframe
*fr1
, *fr2
;
1361 fprintf(stdout
, "comparing energy file %s and %s\n\n", fn1
, fn2
);
1363 in1
= open_enx(fn1
, "r");
1364 in2
= open_enx(fn2
, "r");
1365 do_enxnms(in1
, &nre1
, &enm1
);
1366 do_enxnms(in2
, &nre2
, &enm2
);
1369 fprintf(stdout
, "There are %d and %d terms in the energy files\n\n",
1374 fprintf(stdout
, "There are %d terms in the energy files\n\n", nre1
);
1381 for (i
= 0; i
< nre1
; i
++)
1383 for (j
= 0; j
< nre2
; j
++)
1385 if (enernm_equal(enm1
[i
].name
, enm2
[j
].name
))
1394 if (nre
== 0 || ind1
[nre
-1] != i
)
1396 cmp_str(stdout
, "enm", i
, enm1
[i
].name
, "-");
1399 for (i
= 0; i
< nre2
; i
++)
1403 cmp_str(stdout
, "enm", i
, "-", enm2
[i
].name
);
1408 for (i
= 0; i
< nre
; i
++)
1410 if ((lastener
!= NULL
) && (std::strstr(enm1
[i
].name
, lastener
) != NULL
))
1417 fprintf(stdout
, "There are %d terms to compare in the energy files\n\n",
1420 for (i
= 0; i
< maxener
; i
++)
1422 cmp_str(stdout
, "unit", i
, enm1
[ind1
[i
]].unit
, enm2
[ind2
[i
]].unit
);
1429 b1
= do_enx(in1
, fr1
);
1430 b2
= do_enx(in2
, fr2
);
1433 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn2
, fn1
);
1437 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn1
, fn2
);
1439 else if (!b1
&& !b2
)
1441 fprintf(stdout
, "\nFiles read successfully\n");
1445 cmp_real(stdout
, "t", -1, fr1
->t
, fr2
->t
, ftol
, abstol
);
1446 cmp_int(stdout
, "step", -1, fr1
->step
, fr2
->step
);
1447 /* We don't want to print the nre mismatch for every frame */
1448 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1449 if ((fr1
->nre
>= nre
) && (fr2
->nre
>= nre
))
1451 cmp_energies(stdout
, fr1
->step
, fr1
->step
, fr1
->ener
, fr2
->ener
,
1452 enm1
, ftol
, abstol
, nre
, ind1
, ind2
, maxener
);
1454 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1455 cmp_eblocks(fr1
, fr2
, ftol
, abstol
);