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,2016, 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/trxio.h"
50 #include "gromacs/mdrunutility/mdmodules.h"
51 #include "gromacs/mdtypes/inputrec.h"
52 #include "gromacs/mdtypes/md_enums.h"
53 #include "gromacs/mdtypes/pull-params.h"
54 #include "gromacs/topology/ifunc.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/trajectory/trajectoryframe.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
62 #include "gromacs/utility/stringutil.h"
64 static void cmp_int(FILE *fp
, const char *s
, int index
, int i1
, int i2
)
70 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
74 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
79 static void cmp_int64(FILE *fp
, const char *s
, gmx_int64_t i1
, gmx_int64_t i2
)
83 fprintf(fp
, "%s (", s
);
84 fprintf(fp
, "%" GMX_PRId64
, i1
);
86 fprintf(fp
, "%" GMX_PRId64
, i2
);
91 static void cmp_us(FILE *fp
, const char *s
, int index
, unsigned short i1
, unsigned short i2
)
97 fprintf(fp
, "%s[%d] (%hu - %hu)\n", s
, index
, i1
, i2
);
101 fprintf(fp
, "%s (%hu - %hu)\n", s
, i1
, i2
);
106 static void cmp_uc(FILE *fp
, const char *s
, int index
, unsigned char i1
, unsigned char i2
)
112 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
116 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
121 static gmx_bool
cmp_bool(FILE *fp
, const char *s
, int index
, gmx_bool b1
, gmx_bool b2
)
143 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
,
144 gmx::boolToString(b1
), gmx::boolToString(b2
));
148 fprintf(fp
, "%s (%s - %s)\n", s
,
149 gmx::boolToString(b1
), gmx::boolToString(b2
));
155 static void cmp_str(FILE *fp
, const char *s
, int index
,
156 const char *s1
, const char *s2
)
158 if (std::strcmp(s1
, s2
) != 0)
162 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
, s1
, s2
);
166 fprintf(fp
, "%s (%s - %s)\n", s
, s1
, s2
);
171 static gmx_bool
equal_real(real i1
, real i2
, real ftol
, real abstol
)
173 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
176 static gmx_bool
equal_float(float i1
, float i2
, float ftol
, float abstol
)
178 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
181 static gmx_bool
equal_double(double i1
, double i2
, real ftol
, real abstol
)
183 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
187 cmp_real(FILE *fp
, const char *s
, int index
, real i1
, real i2
, real ftol
, real abstol
)
189 if (!equal_real(i1
, i2
, ftol
, abstol
))
193 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
197 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
203 cmp_float(FILE *fp
, const char *s
, int index
, float i1
, float i2
, float ftol
, float abstol
)
205 if (!equal_float(i1
, i2
, ftol
, abstol
))
209 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
213 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
221 cmp_double(FILE *fp
, const char *s
, int index
, double i1
, double i2
, double ftol
, double abstol
)
223 if (!equal_double(i1
, i2
, ftol
, abstol
))
227 fprintf(fp
, "%s[%2d] (%16.9e - %16.9e)\n", s
, index
, i1
, i2
);
231 fprintf(fp
, "%s (%16.9e - %16.9e)\n", s
, i1
, i2
);
236 static void cmp_rvec(FILE *fp
, const char *s
, int index
, const rvec i1
, const rvec i2
, real ftol
, real abstol
)
238 if (!equal_real(i1
[XX
], i2
[XX
], ftol
, abstol
) ||
239 !equal_real(i1
[YY
], i2
[YY
], ftol
, abstol
) ||
240 !equal_real(i1
[ZZ
], i2
[ZZ
], ftol
, abstol
))
244 fprintf(fp
, "%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
245 s
, index
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
249 fprintf(fp
, "%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
250 s
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
255 static void cmp_ivec(FILE *fp
, const char *s
, int index
, const ivec i1
, const ivec i2
)
257 if ((i1
[XX
] != i2
[XX
]) || (i1
[YY
] != i2
[YY
]) || (i1
[ZZ
] != i2
[ZZ
]))
261 fprintf(fp
, "%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
, index
,
262 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
266 fprintf(fp
, "%s (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
,
267 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
272 static void cmp_ilist(FILE *fp
, int ftype
, const t_ilist
*il1
, const t_ilist
*il2
)
277 fprintf(fp
, "comparing ilist %s\n", interaction_function
[ftype
].name
);
278 sprintf(buf
, "%s->nr", interaction_function
[ftype
].name
);
279 cmp_int(fp
, buf
, -1, il1
->nr
, il2
->nr
);
280 sprintf(buf
, "%s->iatoms", interaction_function
[ftype
].name
);
281 if (((il1
->nr
> 0) && (!il1
->iatoms
)) ||
282 ((il2
->nr
> 0) && (!il2
->iatoms
)) ||
283 ((il1
->nr
!= il2
->nr
)))
285 fprintf(fp
, "Comparing radically different topologies - %s is different\n",
290 for (i
= 0; (i
< il1
->nr
); i
++)
292 cmp_int(fp
, buf
, i
, il1
->iatoms
[i
], il2
->iatoms
[i
]);
297 void cmp_iparm(FILE *fp
, const char *s
, t_functype ft
,
298 const t_iparams
&ip1
, const t_iparams
&ip2
, real ftol
, real abstol
)
304 for (i
= 0; i
< MAXFORCEPARAM
&& !bDiff
; i
++)
306 bDiff
= !equal_real(ip1
.generic
.buf
[i
], ip2
.generic
.buf
[i
], ftol
, abstol
);
310 fprintf(fp
, "%s1: ", s
);
311 pr_iparams(fp
, ft
, &ip1
);
312 fprintf(fp
, "%s2: ", s
);
313 pr_iparams(fp
, ft
, &ip2
);
317 void cmp_iparm_AB(FILE *fp
, const char *s
, t_functype ft
, const t_iparams
&ip1
, real ftol
, real abstol
)
319 int nrfpA
, nrfpB
, p0
, i
;
322 /* Normally the first parameter is perturbable */
324 nrfpA
= interaction_function
[ft
].nrfpA
;
325 nrfpB
= interaction_function
[ft
].nrfpB
;
330 else if (interaction_function
[ft
].flags
& IF_TABULATED
)
332 /* For tabulated interactions only the second parameter is perturbable */
337 for (i
= 0; i
< nrfpB
&& !bDiff
; i
++)
339 bDiff
= !equal_real(ip1
.generic
.buf
[p0
+i
], ip1
.generic
.buf
[nrfpA
+i
], ftol
, abstol
);
343 fprintf(fp
, "%s: ", s
);
344 pr_iparams(fp
, ft
, &ip1
);
348 static void cmp_cmap(FILE *fp
, const gmx_cmap_t
*cmap1
, const gmx_cmap_t
*cmap2
, real ftol
, real abstol
)
350 cmp_int(fp
, "cmap ngrid", -1, cmap1
->ngrid
, cmap2
->ngrid
);
351 cmp_int(fp
, "cmap grid_spacing", -1, cmap1
->grid_spacing
, cmap2
->grid_spacing
);
352 if (cmap1
->ngrid
== cmap2
->ngrid
&&
353 cmap1
->grid_spacing
== cmap2
->grid_spacing
)
357 for (g
= 0; g
< cmap1
->ngrid
; g
++)
361 fprintf(fp
, "comparing cmap %d\n", g
);
363 for (i
= 0; i
< 4*cmap1
->grid_spacing
*cmap1
->grid_spacing
; i
++)
365 cmp_real(fp
, "", i
, cmap1
->cmapdata
[g
].cmap
[i
], cmap2
->cmapdata
[g
].cmap
[i
], ftol
, abstol
);
371 static void cmp_idef(FILE *fp
, const t_idef
*id1
, const t_idef
*id2
, real ftol
, real abstol
)
374 char buf1
[64], buf2
[64];
376 fprintf(fp
, "comparing idef\n");
379 cmp_int(fp
, "idef->ntypes", -1, id1
->ntypes
, id2
->ntypes
);
380 cmp_int(fp
, "idef->atnr", -1, id1
->atnr
, id2
->atnr
);
381 for (i
= 0; (i
< std::min(id1
->ntypes
, id2
->ntypes
)); i
++)
383 sprintf(buf1
, "idef->functype[%d]", i
);
384 sprintf(buf2
, "idef->iparam[%d]", i
);
385 cmp_int(fp
, buf1
, i
, (int)id1
->functype
[i
], (int)id2
->functype
[i
]);
386 cmp_iparm(fp
, buf2
, id1
->functype
[i
],
387 id1
->iparams
[i
], id2
->iparams
[i
], ftol
, abstol
);
389 cmp_real(fp
, "fudgeQQ", -1, id1
->fudgeQQ
, id2
->fudgeQQ
, ftol
, abstol
);
390 cmp_cmap(fp
, &id1
->cmap_grid
, &id2
->cmap_grid
, ftol
, abstol
);
391 for (i
= 0; (i
< F_NRE
); i
++)
393 cmp_ilist(fp
, i
, &(id1
->il
[i
]), &(id2
->il
[i
]));
398 for (i
= 0; (i
< id1
->ntypes
); i
++)
400 cmp_iparm_AB(fp
, "idef->iparam", id1
->functype
[i
], id1
->iparams
[i
], ftol
, abstol
);
405 static void cmp_block(FILE *fp
, const t_block
*b1
, const t_block
*b2
, const char *s
)
409 fprintf(fp
, "comparing block %s\n", s
);
410 sprintf(buf
, "%s.nr", s
);
411 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
414 static void cmp_blocka(FILE *fp
, const t_blocka
*b1
, const t_blocka
*b2
, const char *s
)
418 fprintf(fp
, "comparing blocka %s\n", s
);
419 sprintf(buf
, "%s.nr", s
);
420 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
421 sprintf(buf
, "%s.nra", s
);
422 cmp_int(fp
, buf
, -1, b1
->nra
, b2
->nra
);
425 static void cmp_atom(FILE *fp
, int index
, const t_atom
*a1
, const t_atom
*a2
, real ftol
, real abstol
)
429 cmp_us(fp
, "atom.type", index
, a1
->type
, a2
->type
);
430 cmp_us(fp
, "atom.ptype", index
, a1
->ptype
, a2
->ptype
);
431 cmp_int(fp
, "atom.resind", index
, a1
->resind
, a2
->resind
);
432 cmp_int(fp
, "atom.atomnumber", index
, a1
->atomnumber
, a2
->atomnumber
);
433 cmp_real(fp
, "atom.m", index
, a1
->m
, a2
->m
, ftol
, abstol
);
434 cmp_real(fp
, "atom.q", index
, a1
->q
, a2
->q
, ftol
, abstol
);
435 cmp_us(fp
, "atom.typeB", index
, a1
->typeB
, a2
->typeB
);
436 cmp_real(fp
, "atom.mB", index
, a1
->mB
, a2
->mB
, ftol
, abstol
);
437 cmp_real(fp
, "atom.qB", index
, a1
->qB
, a2
->qB
, ftol
, abstol
);
441 cmp_us(fp
, "atom.type", index
, a1
->type
, a1
->typeB
);
442 cmp_real(fp
, "atom.m", index
, a1
->m
, a1
->mB
, ftol
, abstol
);
443 cmp_real(fp
, "atom.q", index
, a1
->q
, a1
->qB
, ftol
, abstol
);
447 static void cmp_atoms(FILE *fp
, const t_atoms
*a1
, const t_atoms
*a2
, real ftol
, real abstol
)
451 fprintf(fp
, "comparing atoms\n");
455 cmp_int(fp
, "atoms->nr", -1, a1
->nr
, a2
->nr
);
456 for (i
= 0; (i
< a1
->nr
); i
++)
458 cmp_atom(fp
, i
, &(a1
->atom
[i
]), &(a2
->atom
[i
]), ftol
, abstol
);
463 for (i
= 0; (i
< a1
->nr
); i
++)
465 cmp_atom(fp
, i
, &(a1
->atom
[i
]), NULL
, ftol
, abstol
);
470 static void cmp_top(FILE *fp
, const t_topology
*t1
, const t_topology
*t2
, real ftol
, real abstol
)
472 fprintf(fp
, "comparing top\n");
475 cmp_idef(fp
, &(t1
->idef
), &(t2
->idef
), ftol
, abstol
);
476 cmp_atoms(fp
, &(t1
->atoms
), &(t2
->atoms
), ftol
, abstol
);
477 cmp_block(fp
, &t1
->cgs
, &t2
->cgs
, "cgs");
478 cmp_block(fp
, &t1
->mols
, &t2
->mols
, "mols");
479 cmp_bool(fp
, "bIntermolecularInteractions", -1, t1
->bIntermolecularInteractions
, t2
->bIntermolecularInteractions
);
480 cmp_blocka(fp
, &t1
->excls
, &t2
->excls
, "excls");
484 cmp_idef(fp
, &(t1
->idef
), NULL
, ftol
, abstol
);
485 cmp_atoms(fp
, &(t1
->atoms
), NULL
, ftol
, abstol
);
489 static void cmp_groups(FILE *fp
, const gmx_groups_t
*g0
, const gmx_groups_t
*g1
,
490 int natoms0
, int natoms1
)
495 fprintf(fp
, "comparing groups\n");
497 for (i
= 0; i
< egcNR
; i
++)
499 sprintf(buf
, "grps[%d].nr", i
);
500 cmp_int(fp
, buf
, -1, g0
->grps
[i
].nr
, g1
->grps
[i
].nr
);
501 if (g0
->grps
[i
].nr
== g1
->grps
[i
].nr
)
503 for (j
= 0; j
< g0
->grps
[i
].nr
; j
++)
505 sprintf(buf
, "grps[%d].name[%d]", i
, j
);
507 *g0
->grpname
[g0
->grps
[i
].nm_ind
[j
]],
508 *g1
->grpname
[g1
->grps
[i
].nm_ind
[j
]]);
511 cmp_int(fp
, "ngrpnr", i
, g0
->ngrpnr
[i
], g1
->ngrpnr
[i
]);
512 if (g0
->ngrpnr
[i
] == g1
->ngrpnr
[i
] && natoms0
== natoms1
&&
513 (g0
->grpnr
[i
] != NULL
|| g1
->grpnr
[i
] != NULL
))
515 for (j
= 0; j
< natoms0
; j
++)
517 cmp_int(fp
, gtypes
[i
], j
, ggrpnr(g0
, i
, j
), ggrpnr(g1
, i
, j
));
521 /* We have compared the names in the groups lists,
522 * so we can skip the grpname list comparison.
526 static void cmp_rvecs_rmstol(FILE *fp
, const char *title
, int n
, const rvec x1
[], const rvec x2
[],
527 real ftol
, real abstol
)
532 /* For a vector you are usally not interested in a relative difference
533 * on a component that is very small compared to the other components.
534 * Therefore we do the relative comparision relative to the RMS component.
537 for (i
= 0; (i
< n
); i
++)
539 for (m
= 0; m
< DIM
; m
++)
541 rms
+= x1
[i
][m
]*x1
[i
][m
] + x2
[i
][m
]*x2
[i
][m
];
544 rms
= sqrt(rms
/(2*n
*DIM
));
546 /* Convert the relative tolerance into an absolute tolerance */
547 if (ftol
*rms
< abstol
)
552 /* And now do the actual comparision */
553 for (i
= 0; (i
< n
); i
++)
555 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], 0.0, abstol
);
559 static void cmp_rvecs(FILE *fp
, const char *title
, int n
, const rvec x1
[], const rvec x2
[],
560 gmx_bool bRMSD
, real ftol
, real abstol
)
568 for (i
= 0; (i
< n
); i
++)
570 for (m
= 0; m
< DIM
; m
++)
572 d
= x1
[i
][m
] - x2
[i
][m
];
576 fprintf(fp
, "%s RMSD %g\n", title
, std::sqrt(ssd
/n
));
580 cmp_rvecs_rmstol(fp
, title
, n
, x1
, x2
, ftol
, abstol
);
584 static void cmp_grpopts(FILE *fp
, const t_grpopts
*opt1
, const t_grpopts
*opt2
, real ftol
, real abstol
)
587 char buf1
[256], buf2
[256];
589 cmp_int(fp
, "inputrec->grpopts.ngtc", -1, opt1
->ngtc
, opt2
->ngtc
);
590 cmp_int(fp
, "inputrec->grpopts.ngacc", -1, opt1
->ngacc
, opt2
->ngacc
);
591 cmp_int(fp
, "inputrec->grpopts.ngfrz", -1, opt1
->ngfrz
, opt2
->ngfrz
);
592 cmp_int(fp
, "inputrec->grpopts.ngener", -1, opt1
->ngener
, opt2
->ngener
);
593 for (i
= 0; (i
< std::min(opt1
->ngtc
, opt2
->ngtc
)); i
++)
595 cmp_real(fp
, "inputrec->grpopts.nrdf", i
, opt1
->nrdf
[i
], opt2
->nrdf
[i
], ftol
, abstol
);
596 cmp_real(fp
, "inputrec->grpopts.ref_t", i
, opt1
->ref_t
[i
], opt2
->ref_t
[i
], ftol
, abstol
);
597 cmp_real(fp
, "inputrec->grpopts.tau_t", i
, opt1
->tau_t
[i
], opt2
->tau_t
[i
], ftol
, abstol
);
598 cmp_int(fp
, "inputrec->grpopts.annealing", i
, opt1
->annealing
[i
], opt2
->annealing
[i
]);
599 cmp_int(fp
, "inputrec->grpopts.anneal_npoints", i
,
600 opt1
->anneal_npoints
[i
], opt2
->anneal_npoints
[i
]);
601 if (opt1
->anneal_npoints
[i
] == opt2
->anneal_npoints
[i
])
603 sprintf(buf1
, "inputrec->grpopts.anneal_time[%d]", i
);
604 sprintf(buf2
, "inputrec->grpopts.anneal_temp[%d]", i
);
605 for (j
= 0; j
< opt1
->anneal_npoints
[i
]; j
++)
607 cmp_real(fp
, buf1
, j
, opt1
->anneal_time
[i
][j
], opt2
->anneal_time
[i
][j
], ftol
, abstol
);
608 cmp_real(fp
, buf2
, j
, opt1
->anneal_temp
[i
][j
], opt2
->anneal_temp
[i
][j
], ftol
, abstol
);
612 if (opt1
->ngener
== opt2
->ngener
)
614 for (i
= 0; i
< opt1
->ngener
; i
++)
616 for (j
= i
; j
< opt1
->ngener
; j
++)
618 sprintf(buf1
, "inputrec->grpopts.egp_flags[%d]", i
);
620 opt1
->egp_flags
[opt1
->ngener
*i
+j
],
621 opt2
->egp_flags
[opt1
->ngener
*i
+j
]);
625 for (i
= 0; (i
< std::min(opt1
->ngacc
, opt2
->ngacc
)); i
++)
627 cmp_rvec(fp
, "inputrec->grpopts.acc", i
, opt1
->acc
[i
], opt2
->acc
[i
], ftol
, abstol
);
629 for (i
= 0; (i
< std::min(opt1
->ngfrz
, opt2
->ngfrz
)); i
++)
631 cmp_ivec(fp
, "inputrec->grpopts.nFreeze", i
, opt1
->nFreeze
[i
], opt2
->nFreeze
[i
]);
635 static void cmp_cosines(FILE *fp
, const char *s
, const t_cosines c1
[DIM
], const t_cosines c2
[DIM
], real ftol
, real abstol
)
640 for (m
= 0; (m
< DIM
); m
++)
642 sprintf(buf
, "inputrec->%s[%d]", s
, m
);
643 cmp_int(fp
, buf
, 0, c1
->n
, c2
->n
);
644 for (i
= 0; (i
< std::min(c1
->n
, c2
->n
)); i
++)
646 cmp_real(fp
, buf
, i
, c1
->a
[i
], c2
->a
[i
], ftol
, abstol
);
647 cmp_real(fp
, buf
, i
, c1
->phi
[i
], c2
->phi
[i
], ftol
, abstol
);
651 static void cmp_pull(FILE *fp
)
653 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");
656 static void cmp_simtempvals(FILE *fp
, const t_simtemp
*simtemp1
, const t_simtemp
*simtemp2
, int n_lambda
, real ftol
, real abstol
)
659 cmp_int(fp
, "inputrec->simtempvals->eSimTempScale", -1, simtemp1
->eSimTempScale
, simtemp2
->eSimTempScale
);
660 cmp_real(fp
, "inputrec->simtempvals->simtemp_high", -1, simtemp1
->simtemp_high
, simtemp2
->simtemp_high
, ftol
, abstol
);
661 cmp_real(fp
, "inputrec->simtempvals->simtemp_low", -1, simtemp1
->simtemp_low
, simtemp2
->simtemp_low
, ftol
, abstol
);
662 for (i
= 0; i
< n_lambda
; i
++)
664 cmp_real(fp
, "inputrec->simtempvals->temperatures", -1, simtemp1
->temperatures
[i
], simtemp2
->temperatures
[i
], ftol
, abstol
);
668 static void cmp_expandedvals(FILE *fp
, const t_expanded
*expand1
, const t_expanded
*expand2
, int n_lambda
, real ftol
, real abstol
)
672 cmp_bool(fp
, "inputrec->fepvals->bInit_weights", -1, expand1
->bInit_weights
, expand2
->bInit_weights
);
673 cmp_bool(fp
, "inputrec->fepvals->bWLoneovert", -1, expand1
->bWLoneovert
, expand2
->bWLoneovert
);
675 for (i
= 0; i
< n_lambda
; i
++)
677 cmp_real(fp
, "inputrec->expandedvals->init_lambda_weights", -1,
678 expand1
->init_lambda_weights
[i
], expand2
->init_lambda_weights
[i
], ftol
, abstol
);
681 cmp_int(fp
, "inputrec->expandedvals->lambda-stats", -1, expand1
->elamstats
, expand2
->elamstats
);
682 cmp_int(fp
, "inputrec->expandedvals->lambda-mc-move", -1, expand1
->elmcmove
, expand2
->elmcmove
);
683 cmp_int(fp
, "inputrec->expandedvals->lmc-repeats", -1, expand1
->lmc_repeats
, expand2
->lmc_repeats
);
684 cmp_int(fp
, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1
->gibbsdeltalam
, expand2
->gibbsdeltalam
);
685 cmp_int(fp
, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1
->lmc_forced_nstart
, expand2
->lmc_forced_nstart
);
686 cmp_int(fp
, "inputrec->expandedvals->lambda-weights-equil", -1, expand1
->elmceq
, expand2
->elmceq
);
687 cmp_int(fp
, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1
->equil_n_at_lam
, expand2
->equil_n_at_lam
);
688 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1
->equil_samples
, expand2
->equil_samples
);
689 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1
->equil_steps
, expand2
->equil_steps
);
690 cmp_real(fp
, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1
->equil_wl_delta
, expand2
->equil_wl_delta
, ftol
, abstol
);
691 cmp_real(fp
, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1
->equil_ratio
, expand2
->equil_ratio
, ftol
, abstol
);
692 cmp_bool(fp
, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1
->bSymmetrizedTMatrix
, expand2
->bSymmetrizedTMatrix
);
693 cmp_int(fp
, "inputrec->expandedvals->nstTij", -1, expand1
->nstTij
, expand2
->nstTij
);
694 cmp_int(fp
, "inputrec->expandedvals->mininum-var-min", -1, expand1
->minvarmin
, expand2
->minvarmin
); /*default is reasonable */
695 cmp_int(fp
, "inputrec->expandedvals->weight-c-range", -1, expand1
->c_range
, expand2
->c_range
); /* default is just C=0 */
696 cmp_real(fp
, "inputrec->expandedvals->wl-scale", -1, expand1
->wl_scale
, expand2
->wl_scale
, ftol
, abstol
);
697 cmp_real(fp
, "inputrec->expandedvals->init-wl-delta", -1, expand1
->init_wl_delta
, expand2
->init_wl_delta
, ftol
, abstol
);
698 cmp_real(fp
, "inputrec->expandedvals->wl-ratio", -1, expand1
->wl_ratio
, expand2
->wl_ratio
, ftol
, abstol
);
699 cmp_int(fp
, "inputrec->expandedvals->nstexpanded", -1, expand1
->nstexpanded
, expand2
->nstexpanded
);
700 cmp_int(fp
, "inputrec->expandedvals->lmc-seed", -1, expand1
->lmc_seed
, expand2
->lmc_seed
);
701 cmp_real(fp
, "inputrec->expandedvals->mc-temperature", -1, expand1
->mc_temp
, expand2
->mc_temp
, ftol
, abstol
);
704 static void cmp_fepvals(FILE *fp
, const t_lambda
*fep1
, const t_lambda
*fep2
, real ftol
, real abstol
)
707 cmp_int(fp
, "inputrec->nstdhdl", -1, fep1
->nstdhdl
, fep2
->nstdhdl
);
708 cmp_double(fp
, "inputrec->fepvals->init_fep_state", -1, fep1
->init_fep_state
, fep2
->init_fep_state
, ftol
, abstol
);
709 cmp_double(fp
, "inputrec->fepvals->delta_lambda", -1, fep1
->delta_lambda
, fep2
->delta_lambda
, ftol
, abstol
);
710 cmp_int(fp
, "inputrec->fepvals->n_lambda", -1, fep1
->n_lambda
, fep2
->n_lambda
);
711 for (i
= 0; i
< efptNR
; i
++)
713 for (j
= 0; j
< std::min(fep1
->n_lambda
, fep2
->n_lambda
); j
++)
715 cmp_double(fp
, "inputrec->fepvals->all_lambda", -1, fep1
->all_lambda
[i
][j
], fep2
->all_lambda
[i
][j
], ftol
, abstol
);
718 cmp_int(fp
, "inputrec->fepvals->lambda_neighbors", 1, fep1
->lambda_neighbors
,
719 fep2
->lambda_neighbors
);
720 cmp_real(fp
, "inputrec->fepvals->sc_alpha", -1, fep1
->sc_alpha
, fep2
->sc_alpha
, ftol
, abstol
);
721 cmp_int(fp
, "inputrec->fepvals->sc_power", -1, fep1
->sc_power
, fep2
->sc_power
);
722 cmp_real(fp
, "inputrec->fepvals->sc_r_power", -1, fep1
->sc_r_power
, fep2
->sc_r_power
, ftol
, abstol
);
723 cmp_real(fp
, "inputrec->fepvals->sc_sigma", -1, fep1
->sc_sigma
, fep2
->sc_sigma
, ftol
, abstol
);
724 cmp_int(fp
, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1
->edHdLPrintEnergy
, fep1
->edHdLPrintEnergy
);
725 cmp_bool(fp
, "inputrec->fepvals->bScCoul", -1, fep1
->bScCoul
, fep1
->bScCoul
);
726 cmp_int(fp
, "inputrec->separate_dhdl_file", -1, fep1
->separate_dhdl_file
, fep2
->separate_dhdl_file
);
727 cmp_int(fp
, "inputrec->dhdl_derivatives", -1, fep1
->dhdl_derivatives
, fep2
->dhdl_derivatives
);
728 cmp_int(fp
, "inputrec->dh_hist_size", -1, fep1
->dh_hist_size
, fep2
->dh_hist_size
);
729 cmp_double(fp
, "inputrec->dh_hist_spacing", -1, fep1
->dh_hist_spacing
, fep2
->dh_hist_spacing
, ftol
, abstol
);
732 static void cmp_inputrec(FILE *fp
, const t_inputrec
*ir1
, const t_inputrec
*ir2
, real ftol
, real abstol
)
734 fprintf(fp
, "comparing inputrec\n");
736 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
737 * of warnings. Maybe it will change in future versions, but for the
738 * moment I've spelled them out instead. /EL 000820
739 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
740 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
741 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
743 cmp_int(fp
, "inputrec->eI", -1, ir1
->eI
, ir2
->eI
);
744 cmp_int64(fp
, "inputrec->nsteps", ir1
->nsteps
, ir2
->nsteps
);
745 cmp_int64(fp
, "inputrec->init_step", ir1
->init_step
, ir2
->init_step
);
746 cmp_int(fp
, "inputrec->simulation_part", -1, ir1
->simulation_part
, ir2
->simulation_part
);
747 cmp_int(fp
, "inputrec->ePBC", -1, ir1
->ePBC
, ir2
->ePBC
);
748 cmp_int(fp
, "inputrec->bPeriodicMols", -1, ir1
->bPeriodicMols
, ir2
->bPeriodicMols
);
749 cmp_int(fp
, "inputrec->cutoff_scheme", -1, ir1
->cutoff_scheme
, ir2
->cutoff_scheme
);
750 cmp_int(fp
, "inputrec->ns_type", -1, ir1
->ns_type
, ir2
->ns_type
);
751 cmp_int(fp
, "inputrec->nstlist", -1, ir1
->nstlist
, ir2
->nstlist
);
752 cmp_int(fp
, "inputrec->nstcomm", -1, ir1
->nstcomm
, ir2
->nstcomm
);
753 cmp_int(fp
, "inputrec->comm_mode", -1, ir1
->comm_mode
, ir2
->comm_mode
);
754 cmp_int(fp
, "inputrec->nstlog", -1, ir1
->nstlog
, ir2
->nstlog
);
755 cmp_int(fp
, "inputrec->nstxout", -1, ir1
->nstxout
, ir2
->nstxout
);
756 cmp_int(fp
, "inputrec->nstvout", -1, ir1
->nstvout
, ir2
->nstvout
);
757 cmp_int(fp
, "inputrec->nstfout", -1, ir1
->nstfout
, ir2
->nstfout
);
758 cmp_int(fp
, "inputrec->nstcalcenergy", -1, ir1
->nstcalcenergy
, ir2
->nstcalcenergy
);
759 cmp_int(fp
, "inputrec->nstenergy", -1, ir1
->nstenergy
, ir2
->nstenergy
);
760 cmp_int(fp
, "inputrec->nstxout_compressed", -1, ir1
->nstxout_compressed
, ir2
->nstxout_compressed
);
761 cmp_double(fp
, "inputrec->init_t", -1, ir1
->init_t
, ir2
->init_t
, ftol
, abstol
);
762 cmp_double(fp
, "inputrec->delta_t", -1, ir1
->delta_t
, ir2
->delta_t
, ftol
, abstol
);
763 cmp_real(fp
, "inputrec->x_compression_precision", -1, ir1
->x_compression_precision
, ir2
->x_compression_precision
, ftol
, abstol
);
764 cmp_real(fp
, "inputrec->fourierspacing", -1, ir1
->fourier_spacing
, ir2
->fourier_spacing
, ftol
, abstol
);
765 cmp_int(fp
, "inputrec->nkx", -1, ir1
->nkx
, ir2
->nkx
);
766 cmp_int(fp
, "inputrec->nky", -1, ir1
->nky
, ir2
->nky
);
767 cmp_int(fp
, "inputrec->nkz", -1, ir1
->nkz
, ir2
->nkz
);
768 cmp_int(fp
, "inputrec->pme_order", -1, ir1
->pme_order
, ir2
->pme_order
);
769 cmp_real(fp
, "inputrec->ewald_rtol", -1, ir1
->ewald_rtol
, ir2
->ewald_rtol
, ftol
, abstol
);
770 cmp_int(fp
, "inputrec->ewald_geometry", -1, ir1
->ewald_geometry
, ir2
->ewald_geometry
);
771 cmp_real(fp
, "inputrec->epsilon_surface", -1, ir1
->epsilon_surface
, ir2
->epsilon_surface
, ftol
, abstol
);
772 cmp_int(fp
, "inputrec->bContinuation", -1, ir1
->bContinuation
, ir2
->bContinuation
);
773 cmp_int(fp
, "inputrec->bShakeSOR", -1, ir1
->bShakeSOR
, ir2
->bShakeSOR
);
774 cmp_int(fp
, "inputrec->etc", -1, ir1
->etc
, ir2
->etc
);
775 cmp_int(fp
, "inputrec->bPrintNHChains", -1, ir1
->bPrintNHChains
, ir2
->bPrintNHChains
);
776 cmp_int(fp
, "inputrec->epc", -1, ir1
->epc
, ir2
->epc
);
777 cmp_int(fp
, "inputrec->epct", -1, ir1
->epct
, ir2
->epct
);
778 cmp_real(fp
, "inputrec->tau_p", -1, ir1
->tau_p
, ir2
->tau_p
, ftol
, abstol
);
779 cmp_rvec(fp
, "inputrec->ref_p(x)", -1, ir1
->ref_p
[XX
], ir2
->ref_p
[XX
], ftol
, abstol
);
780 cmp_rvec(fp
, "inputrec->ref_p(y)", -1, ir1
->ref_p
[YY
], ir2
->ref_p
[YY
], ftol
, abstol
);
781 cmp_rvec(fp
, "inputrec->ref_p(z)", -1, ir1
->ref_p
[ZZ
], ir2
->ref_p
[ZZ
], ftol
, abstol
);
782 cmp_rvec(fp
, "inputrec->compress(x)", -1, ir1
->compress
[XX
], ir2
->compress
[XX
], ftol
, abstol
);
783 cmp_rvec(fp
, "inputrec->compress(y)", -1, ir1
->compress
[YY
], ir2
->compress
[YY
], ftol
, abstol
);
784 cmp_rvec(fp
, "inputrec->compress(z)", -1, ir1
->compress
[ZZ
], ir2
->compress
[ZZ
], ftol
, abstol
);
785 cmp_int(fp
, "refcoord_scaling", -1, ir1
->refcoord_scaling
, ir2
->refcoord_scaling
);
786 cmp_rvec(fp
, "inputrec->posres_com", -1, ir1
->posres_com
, ir2
->posres_com
, ftol
, abstol
);
787 cmp_rvec(fp
, "inputrec->posres_comB", -1, ir1
->posres_comB
, ir2
->posres_comB
, ftol
, abstol
);
788 cmp_real(fp
, "inputrec->verletbuf_tol", -1, ir1
->verletbuf_tol
, ir2
->verletbuf_tol
, ftol
, abstol
);
789 cmp_real(fp
, "inputrec->rlist", -1, ir1
->rlist
, ir2
->rlist
, ftol
, abstol
);
790 cmp_real(fp
, "inputrec->rtpi", -1, ir1
->rtpi
, ir2
->rtpi
, ftol
, abstol
);
791 cmp_int(fp
, "inputrec->coulombtype", -1, ir1
->coulombtype
, ir2
->coulombtype
);
792 cmp_int(fp
, "inputrec->coulomb_modifier", -1, ir1
->coulomb_modifier
, ir2
->coulomb_modifier
);
793 cmp_real(fp
, "inputrec->rcoulomb_switch", -1, ir1
->rcoulomb_switch
, ir2
->rcoulomb_switch
, ftol
, abstol
);
794 cmp_real(fp
, "inputrec->rcoulomb", -1, ir1
->rcoulomb
, ir2
->rcoulomb
, ftol
, abstol
);
795 cmp_int(fp
, "inputrec->vdwtype", -1, ir1
->vdwtype
, ir2
->vdwtype
);
796 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
);
797 cmp_real(fp
, "inputrec->rvdw", -1, ir1
->rvdw
, ir2
->rvdw
, ftol
, abstol
);
798 cmp_real(fp
, "inputrec->epsilon_r", -1, ir1
->epsilon_r
, ir2
->epsilon_r
, ftol
, abstol
);
799 cmp_real(fp
, "inputrec->epsilon_rf", -1, ir1
->epsilon_rf
, ir2
->epsilon_rf
, ftol
, abstol
);
800 cmp_real(fp
, "inputrec->tabext", -1, ir1
->tabext
, ir2
->tabext
, ftol
, abstol
);
801 cmp_int(fp
, "inputrec->implicit_solvent", -1, ir1
->implicit_solvent
, ir2
->implicit_solvent
);
802 cmp_int(fp
, "inputrec->gb_algorithm", -1, ir1
->gb_algorithm
, ir2
->gb_algorithm
);
803 cmp_int(fp
, "inputrec->nstgbradii", -1, ir1
->nstgbradii
, ir2
->nstgbradii
);
804 cmp_real(fp
, "inputrec->rgbradii", -1, ir1
->rgbradii
, ir2
->rgbradii
, ftol
, abstol
);
805 cmp_real(fp
, "inputrec->gb_saltconc", -1, ir1
->gb_saltconc
, ir2
->gb_saltconc
, ftol
, abstol
);
806 cmp_real(fp
, "inputrec->gb_epsilon_solvent", -1, ir1
->gb_epsilon_solvent
, ir2
->gb_epsilon_solvent
, ftol
, abstol
);
807 cmp_real(fp
, "inputrec->gb_obc_alpha", -1, ir1
->gb_obc_alpha
, ir2
->gb_obc_alpha
, ftol
, abstol
);
808 cmp_real(fp
, "inputrec->gb_obc_beta", -1, ir1
->gb_obc_beta
, ir2
->gb_obc_beta
, ftol
, abstol
);
809 cmp_real(fp
, "inputrec->gb_obc_gamma", -1, ir1
->gb_obc_gamma
, ir2
->gb_obc_gamma
, ftol
, abstol
);
810 cmp_real(fp
, "inputrec->gb_dielectric_offset", -1, ir1
->gb_dielectric_offset
, ir2
->gb_dielectric_offset
, ftol
, abstol
);
811 cmp_int(fp
, "inputrec->sa_algorithm", -1, ir1
->sa_algorithm
, ir2
->sa_algorithm
);
812 cmp_real(fp
, "inputrec->sa_surface_tension", -1, ir1
->sa_surface_tension
, ir2
->sa_surface_tension
, ftol
, abstol
);
814 cmp_int(fp
, "inputrec->eDispCorr", -1, ir1
->eDispCorr
, ir2
->eDispCorr
);
815 cmp_real(fp
, "inputrec->shake_tol", -1, ir1
->shake_tol
, ir2
->shake_tol
, ftol
, abstol
);
816 cmp_int(fp
, "inputrec->efep", -1, ir1
->efep
, ir2
->efep
);
817 cmp_fepvals(fp
, ir1
->fepvals
, ir2
->fepvals
, ftol
, abstol
);
818 cmp_int(fp
, "inputrec->bSimTemp", -1, ir1
->bSimTemp
, ir2
->bSimTemp
);
819 if ((ir1
->bSimTemp
== ir2
->bSimTemp
) && (ir1
->bSimTemp
))
821 cmp_simtempvals(fp
, ir1
->simtempvals
, ir2
->simtempvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
823 cmp_int(fp
, "inputrec->bExpanded", -1, ir1
->bExpanded
, ir2
->bExpanded
);
824 if ((ir1
->bExpanded
== ir2
->bExpanded
) && (ir1
->bExpanded
))
826 cmp_expandedvals(fp
, ir1
->expandedvals
, ir2
->expandedvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
828 cmp_int(fp
, "inputrec->nwall", -1, ir1
->nwall
, ir2
->nwall
);
829 cmp_int(fp
, "inputrec->wall_type", -1, ir1
->wall_type
, ir2
->wall_type
);
830 cmp_int(fp
, "inputrec->wall_atomtype[0]", -1, ir1
->wall_atomtype
[0], ir2
->wall_atomtype
[0]);
831 cmp_int(fp
, "inputrec->wall_atomtype[1]", -1, ir1
->wall_atomtype
[1], ir2
->wall_atomtype
[1]);
832 cmp_real(fp
, "inputrec->wall_density[0]", -1, ir1
->wall_density
[0], ir2
->wall_density
[0], ftol
, abstol
);
833 cmp_real(fp
, "inputrec->wall_density[1]", -1, ir1
->wall_density
[1], ir2
->wall_density
[1], ftol
, abstol
);
834 cmp_real(fp
, "inputrec->wall_ewald_zfac", -1, ir1
->wall_ewald_zfac
, ir2
->wall_ewald_zfac
, ftol
, abstol
);
836 cmp_bool(fp
, "inputrec->bPull", -1, ir1
->bPull
, ir2
->bPull
);
837 if (ir1
->bPull
&& ir2
->bPull
)
842 cmp_int(fp
, "inputrec->eDisre", -1, ir1
->eDisre
, ir2
->eDisre
);
843 cmp_real(fp
, "inputrec->dr_fc", -1, ir1
->dr_fc
, ir2
->dr_fc
, ftol
, abstol
);
844 cmp_int(fp
, "inputrec->eDisreWeighting", -1, ir1
->eDisreWeighting
, ir2
->eDisreWeighting
);
845 cmp_int(fp
, "inputrec->bDisreMixed", -1, ir1
->bDisreMixed
, ir2
->bDisreMixed
);
846 cmp_int(fp
, "inputrec->nstdisreout", -1, ir1
->nstdisreout
, ir2
->nstdisreout
);
847 cmp_real(fp
, "inputrec->dr_tau", -1, ir1
->dr_tau
, ir2
->dr_tau
, ftol
, abstol
);
848 cmp_real(fp
, "inputrec->orires_fc", -1, ir1
->orires_fc
, ir2
->orires_fc
, ftol
, abstol
);
849 cmp_real(fp
, "inputrec->orires_tau", -1, ir1
->orires_tau
, ir2
->orires_tau
, ftol
, abstol
);
850 cmp_int(fp
, "inputrec->nstorireout", -1, ir1
->nstorireout
, ir2
->nstorireout
);
851 cmp_real(fp
, "inputrec->em_stepsize", -1, ir1
->em_stepsize
, ir2
->em_stepsize
, ftol
, abstol
);
852 cmp_real(fp
, "inputrec->em_tol", -1, ir1
->em_tol
, ir2
->em_tol
, ftol
, abstol
);
853 cmp_int(fp
, "inputrec->niter", -1, ir1
->niter
, ir2
->niter
);
854 cmp_real(fp
, "inputrec->fc_stepsize", -1, ir1
->fc_stepsize
, ir2
->fc_stepsize
, ftol
, abstol
);
855 cmp_int(fp
, "inputrec->nstcgsteep", -1, ir1
->nstcgsteep
, ir2
->nstcgsteep
);
856 cmp_int(fp
, "inputrec->nbfgscorr", 0, ir1
->nbfgscorr
, ir2
->nbfgscorr
);
857 cmp_int(fp
, "inputrec->eConstrAlg", -1, ir1
->eConstrAlg
, ir2
->eConstrAlg
);
858 cmp_int(fp
, "inputrec->nProjOrder", -1, ir1
->nProjOrder
, ir2
->nProjOrder
);
859 cmp_real(fp
, "inputrec->LincsWarnAngle", -1, ir1
->LincsWarnAngle
, ir2
->LincsWarnAngle
, ftol
, abstol
);
860 cmp_int(fp
, "inputrec->nLincsIter", -1, ir1
->nLincsIter
, ir2
->nLincsIter
);
861 cmp_real(fp
, "inputrec->bd_fric", -1, ir1
->bd_fric
, ir2
->bd_fric
, ftol
, abstol
);
862 cmp_int64(fp
, "inputrec->ld_seed", ir1
->ld_seed
, ir2
->ld_seed
);
863 cmp_real(fp
, "inputrec->cos_accel", -1, ir1
->cos_accel
, ir2
->cos_accel
, ftol
, abstol
);
864 cmp_rvec(fp
, "inputrec->deform(a)", -1, ir1
->deform
[XX
], ir2
->deform
[XX
], ftol
, abstol
);
865 cmp_rvec(fp
, "inputrec->deform(b)", -1, ir1
->deform
[YY
], ir2
->deform
[YY
], ftol
, abstol
);
866 cmp_rvec(fp
, "inputrec->deform(c)", -1, ir1
->deform
[ZZ
], ir2
->deform
[ZZ
], ftol
, abstol
);
869 cmp_int(fp
, "inputrec->userint1", -1, ir1
->userint1
, ir2
->userint1
);
870 cmp_int(fp
, "inputrec->userint2", -1, ir1
->userint2
, ir2
->userint2
);
871 cmp_int(fp
, "inputrec->userint3", -1, ir1
->userint3
, ir2
->userint3
);
872 cmp_int(fp
, "inputrec->userint4", -1, ir1
->userint4
, ir2
->userint4
);
873 cmp_real(fp
, "inputrec->userreal1", -1, ir1
->userreal1
, ir2
->userreal1
, ftol
, abstol
);
874 cmp_real(fp
, "inputrec->userreal2", -1, ir1
->userreal2
, ir2
->userreal2
, ftol
, abstol
);
875 cmp_real(fp
, "inputrec->userreal3", -1, ir1
->userreal3
, ir2
->userreal3
, ftol
, abstol
);
876 cmp_real(fp
, "inputrec->userreal4", -1, ir1
->userreal4
, ir2
->userreal4
, ftol
, abstol
);
877 cmp_grpopts(fp
, &(ir1
->opts
), &(ir2
->opts
), ftol
, abstol
);
878 cmp_cosines(fp
, "ex", ir1
->ex
, ir2
->ex
, ftol
, abstol
);
879 cmp_cosines(fp
, "et", ir1
->et
, ir2
->et
, ftol
, abstol
);
882 static void comp_pull_AB(FILE *fp
, pull_params_t
*pull
, real ftol
, real abstol
)
886 for (i
= 0; i
< pull
->ncoord
; i
++)
888 fprintf(fp
, "comparing pull coord %d\n", i
);
889 cmp_real(fp
, "pull-coord->k", -1, pull
->coord
[i
].k
, pull
->coord
[i
].kB
, ftol
, abstol
);
893 static void comp_state(const t_state
*st1
, const t_state
*st2
,
894 gmx_bool bRMSD
, real ftol
, real abstol
)
898 fprintf(stdout
, "comparing flags\n");
899 cmp_int(stdout
, "flags", -1, st1
->flags
, st2
->flags
);
900 fprintf(stdout
, "comparing box\n");
901 cmp_rvecs(stdout
, "box", DIM
, st1
->box
, st2
->box
, FALSE
, ftol
, abstol
);
902 fprintf(stdout
, "comparing box_rel\n");
903 cmp_rvecs(stdout
, "box_rel", DIM
, st1
->box_rel
, st2
->box_rel
, FALSE
, ftol
, abstol
);
904 fprintf(stdout
, "comparing boxv\n");
905 cmp_rvecs(stdout
, "boxv", DIM
, st1
->boxv
, st2
->boxv
, FALSE
, ftol
, abstol
);
906 if (st1
->flags
& (1<<estSVIR_PREV
))
908 fprintf(stdout
, "comparing shake vir_prev\n");
909 cmp_rvecs(stdout
, "svir_prev", DIM
, st1
->svir_prev
, st2
->svir_prev
, FALSE
, ftol
, abstol
);
911 if (st1
->flags
& (1<<estFVIR_PREV
))
913 fprintf(stdout
, "comparing force vir_prev\n");
914 cmp_rvecs(stdout
, "fvir_prev", DIM
, st1
->fvir_prev
, st2
->fvir_prev
, FALSE
, ftol
, abstol
);
916 if (st1
->flags
& (1<<estPRES_PREV
))
918 fprintf(stdout
, "comparing prev_pres\n");
919 cmp_rvecs(stdout
, "pres_prev", DIM
, st1
->pres_prev
, st2
->pres_prev
, FALSE
, ftol
, abstol
);
921 cmp_int(stdout
, "ngtc", -1, st1
->ngtc
, st2
->ngtc
);
922 cmp_int(stdout
, "nhchainlength", -1, st1
->nhchainlength
, st2
->nhchainlength
);
923 if (st1
->ngtc
== st2
->ngtc
&& st1
->nhchainlength
== st2
->nhchainlength
)
925 for (i
= 0; i
< st1
->ngtc
; i
++)
927 nc
= i
*st1
->nhchainlength
;
928 for (j
= 0; j
< nc
; j
++)
930 cmp_real(stdout
, "nosehoover_xi",
931 i
, st1
->nosehoover_xi
[nc
+j
], st2
->nosehoover_xi
[nc
+j
], ftol
, abstol
);
935 cmp_int(stdout
, "nnhpres", -1, st1
->nnhpres
, st2
->nnhpres
);
936 if (st1
->nnhpres
== st2
->nnhpres
&& st1
->nhchainlength
== st2
->nhchainlength
)
938 for (i
= 0; i
< st1
->nnhpres
; i
++)
940 nc
= i
*st1
->nhchainlength
;
941 for (j
= 0; j
< nc
; j
++)
943 cmp_real(stdout
, "nosehoover_xi",
944 i
, st1
->nhpres_xi
[nc
+j
], st2
->nhpres_xi
[nc
+j
], ftol
, abstol
);
949 cmp_int(stdout
, "natoms", -1, st1
->natoms
, st2
->natoms
);
950 if (st1
->natoms
== st2
->natoms
)
952 if ((st1
->flags
& (1<<estX
)) && (st2
->flags
& (1<<estX
)))
954 fprintf(stdout
, "comparing x\n");
955 cmp_rvecs(stdout
, "x", st1
->natoms
, st1
->x
, st2
->x
, bRMSD
, ftol
, abstol
);
957 if ((st1
->flags
& (1<<estV
)) && (st2
->flags
& (1<<estV
)))
959 fprintf(stdout
, "comparing v\n");
960 cmp_rvecs(stdout
, "v", st1
->natoms
, st1
->v
, st2
->v
, bRMSD
, ftol
, abstol
);
965 void comp_tpx(const char *fn1
, const char *fn2
,
966 gmx_bool bRMSD
, real ftol
, real abstol
)
969 gmx::MDModules mdModules
[2];
970 t_inputrec
*ir
[2] = { mdModules
[0].inputrec(), mdModules
[1].inputrec() };
978 for (i
= 0; i
< (fn2
? 2 : 1); i
++)
980 read_tpx_state(ff
[i
], ir
[i
], &state
[i
], &(mtop
[i
]));
984 cmp_inputrec(stdout
, ir
[0], ir
[1], ftol
, abstol
);
985 /* Convert gmx_mtop_t to t_topology.
986 * We should implement direct mtop comparison,
987 * but it might be useful to keep t_topology comparison as an option.
989 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
990 top
[1] = gmx_mtop_t_to_t_topology(&mtop
[1]);
991 cmp_top(stdout
, &top
[0], &top
[1], ftol
, abstol
);
992 cmp_groups(stdout
, &mtop
[0].groups
, &mtop
[1].groups
,
993 mtop
[0].natoms
, mtop
[1].natoms
);
994 comp_state(&state
[0], &state
[1], bRMSD
, ftol
, abstol
);
998 if (ir
[0]->efep
== efepNO
)
1000 fprintf(stdout
, "inputrec->efep = %s\n", efep_names
[ir
[0]->efep
]);
1006 comp_pull_AB(stdout
, ir
[0]->pull
, ftol
, abstol
);
1008 /* Convert gmx_mtop_t to t_topology.
1009 * We should implement direct mtop comparison,
1010 * but it might be useful to keep t_topology comparison as an option.
1012 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
1013 cmp_top(stdout
, &top
[0], NULL
, ftol
, abstol
);
1018 void comp_frame(FILE *fp
, t_trxframe
*fr1
, t_trxframe
*fr2
,
1019 gmx_bool bRMSD
, real ftol
, real abstol
)
1022 cmp_int(fp
, "not_ok", -1, fr1
->not_ok
, fr2
->not_ok
);
1023 cmp_int(fp
, "natoms", -1, fr1
->natoms
, fr2
->natoms
);
1024 if (cmp_bool(fp
, "bTitle", -1, fr1
->bTitle
, fr2
->bTitle
))
1026 cmp_str(fp
, "title", -1, fr1
->title
, fr2
->title
);
1028 if (cmp_bool(fp
, "bStep", -1, fr1
->bStep
, fr2
->bStep
))
1030 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1032 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1033 if (cmp_bool(fp
, "bTime", -1, fr1
->bTime
, fr2
->bTime
))
1035 cmp_real(fp
, "time", -1, fr1
->time
, fr2
->time
, ftol
, abstol
);
1037 if (cmp_bool(fp
, "bLambda", -1, fr1
->bLambda
, fr2
->bLambda
))
1039 cmp_real(fp
, "lambda", -1, fr1
->lambda
, fr2
->lambda
, ftol
, abstol
);
1041 if (cmp_bool(fp
, "bAtoms", -1, fr1
->bAtoms
, fr2
->bAtoms
))
1043 cmp_atoms(fp
, fr1
->atoms
, fr2
->atoms
, ftol
, abstol
);
1045 if (cmp_bool(fp
, "bPrec", -1, fr1
->bPrec
, fr2
->bPrec
))
1047 cmp_real(fp
, "prec", -1, fr1
->prec
, fr2
->prec
, ftol
, abstol
);
1049 if (cmp_bool(fp
, "bX", -1, fr1
->bX
, fr2
->bX
))
1051 cmp_rvecs(fp
, "x", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->x
, fr2
->x
, bRMSD
, ftol
, abstol
);
1053 if (cmp_bool(fp
, "bV", -1, fr1
->bV
, fr2
->bV
))
1055 cmp_rvecs(fp
, "v", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->v
, fr2
->v
, bRMSD
, ftol
, abstol
);
1057 if (cmp_bool(fp
, "bF", -1, fr1
->bF
, fr2
->bF
))
1059 cmp_rvecs(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, bRMSD
, ftol
, abstol
);
1061 if (cmp_bool(fp
, "bBox", -1, fr1
->bBox
, fr2
->bBox
))
1063 cmp_rvecs(fp
, "box", 3, fr1
->box
, fr2
->box
, FALSE
, ftol
, abstol
);
1067 void comp_trx(const gmx_output_env_t
*oenv
, const char *fn1
, const char *fn2
,
1068 gmx_bool bRMSD
, real ftol
, real abstol
)
1073 t_trxstatus
*status
[2];
1078 fprintf(stderr
, "Comparing trajectory files %s and %s\n", fn1
, fn2
);
1079 for (i
= 0; i
< 2; i
++)
1081 b
[i
] = read_first_frame(oenv
, &status
[i
], fn
[i
], &fr
[i
], TRX_READ_X
|TRX_READ_V
|TRX_READ_F
);
1088 comp_frame(stdout
, &(fr
[0]), &(fr
[1]), bRMSD
, ftol
, abstol
);
1090 for (i
= 0; i
< 2; i
++)
1092 b
[i
] = read_next_frame(oenv
, status
[i
], &fr
[i
]);
1095 while (b
[0] && b
[1]);
1097 for (i
= 0; i
< 2; i
++)
1099 if (b
[i
] && !b
[1-i
])
1101 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn
[1-i
], fn
[i
]);
1103 close_trj(status
[i
]);
1108 fprintf(stdout
, "\nBoth files read correctly\n");
1112 static real
ener_tensor_diag(int n
, int *ind1
, int *ind2
,
1115 t_energy e1
[], t_energy e2
[])
1124 d2
= tensi
[i
] - d1
*DIM
;
1126 /* Find the diagonal elements d1 and d2 */
1127 len
= std::strlen(enm1
[ind1
[i
]].name
);
1131 for (j
= 0; j
< n
; j
++)
1133 if (tensi
[j
] >= 0 &&
1134 std::strlen(enm1
[ind1
[j
]].name
) == len
&&
1135 std::strncmp(enm1
[ind1
[i
]].name
, enm1
[ind1
[j
]].name
, len
-2) == 0 &&
1136 (tensi
[j
] == d1
*DIM
+d1
|| tensi
[j
] == d2
*DIM
+d2
))
1138 prod1
*= fabs(e1
[ind1
[j
]].e
);
1139 prod2
*= fabs(e2
[ind2
[j
]].e
);
1146 return 0.5*(std::sqrt(prod1
) + std::sqrt(prod2
));
1154 static gmx_bool
enernm_equal(const char *nm1
, const char *nm2
)
1158 len1
= std::strlen(nm1
);
1159 len2
= std::strlen(nm2
);
1161 /* Remove " (bar)" at the end of a name */
1162 if (len1
> 6 && std::strcmp(nm1
+len1
-6, " (bar)") == 0)
1166 if (len2
> 6 && std::strcmp(nm2
+len2
-6, " (bar)") == 0)
1171 return (len1
== len2
&& gmx_strncasecmp(nm1
, nm2
, len1
) == 0);
1174 static void cmp_energies(FILE *fp
, int step1
, int step2
,
1175 t_energy e1
[], t_energy e2
[],
1177 real ftol
, real abstol
,
1178 int nre
, int *ind1
, int *ind2
, int maxener
)
1181 int *tensi
, len
, d1
, d2
;
1182 real ftol_i
, abstol_i
;
1184 snew(tensi
, maxener
);
1185 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1186 for (i
= 0; (i
< maxener
); i
++)
1190 len
= std::strlen(enm1
[ii
].name
);
1191 if (len
> 3 && enm1
[ii
].name
[len
-3] == '-')
1193 d1
= enm1
[ii
].name
[len
-2] - 'X';
1194 d2
= enm1
[ii
].name
[len
-1] - 'X';
1195 if (d1
>= 0 && d1
< DIM
&&
1196 d2
>= 0 && d2
< DIM
)
1198 tensi
[i
] = d1
*DIM
+ d2
;
1203 for (i
= 0; (i
< maxener
); i
++)
1205 /* Check if this is an off-diagonal tensor element */
1206 if (tensi
[i
] >= 0 && tensi
[i
] != 0 && tensi
[i
] != 4 && tensi
[i
] != 8)
1208 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1210 /* Do the relative tolerance through an absolute tolerance times
1211 * the size of diagonal components of the tensor.
1213 abstol_i
= ftol
*ener_tensor_diag(nre
, ind1
, ind2
, enm1
, tensi
, i
, e1
, e2
);
1216 fprintf(debug
, "tensor '%s' val %f diag %f\n",
1217 enm1
[i
].name
, e1
[i
].e
, abstol_i
/ftol
);
1221 /* We found a diagonal, we need to check with the minimum tolerance */
1222 abstol_i
= std::min(abstol_i
, abstol
);
1226 /* We did not find a diagonal, ignore the relative tolerance check */
1235 if (!equal_real(e1
[ind1
[i
]].e
, e2
[ind2
[i
]].e
, ftol_i
, abstol_i
))
1237 fprintf(fp
, "%-15s step %3d: %12g, step %3d: %12g\n",
1239 step1
, e1
[ind1
[i
]].e
,
1240 step2
, e2
[ind2
[i
]].e
);
1248 static void cmp_disres(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1251 char bav
[64], bt
[64], bs
[22];
1253 cmp_int(stdout
, "ndisre", -1, fr1
->ndisre
, fr2
->ndisre
);
1254 if ((fr1
->ndisre
== fr2
->ndisre
) && (fr1
->ndisre
> 0))
1256 sprintf(bav
, "step %s: disre rav", gmx_step_str(fr1
->step
, bs
));
1257 sprintf(bt
, "step %s: disre rt", gmx_step_str(fr1
->step
, bs
));
1258 for (i
= 0; (i
< fr1
->ndisre
); i
++)
1260 cmp_real(stdout
, bav
, i
, fr1
->disre_rm3tav
[i
], fr2
->disre_rm3tav
[i
], ftol
, abstol
);
1261 cmp_real(stdout
, bt
, i
, fr1
->disre_rt
[i
], fr2
->disre_rt
[i
], ftol
, abstol
);
1267 static void cmp_eblocks(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1270 char buf
[64], bs
[22];
1272 cmp_int(stdout
, "nblock", -1, fr1
->nblock
, fr2
->nblock
);
1273 if ((fr1
->nblock
== fr2
->nblock
) && (fr1
->nblock
> 0))
1275 for (j
= 0; (j
< fr1
->nblock
); j
++)
1277 t_enxblock
*b1
, *b2
; /* convenience vars */
1279 b1
= &(fr1
->block
[j
]);
1280 b2
= &(fr2
->block
[j
]);
1282 sprintf(buf
, "step %s: block[%d]", gmx_step_str(fr1
->step
, bs
), j
);
1283 cmp_int(stdout
, buf
, -1, b1
->nsub
, b2
->nsub
);
1284 cmp_int(stdout
, buf
, -1, b1
->id
, b2
->id
);
1286 if ( (b1
->nsub
== b2
->nsub
) && (b1
->id
== b2
->id
) )
1288 for (i
= 0; i
< b1
->nsub
; i
++)
1290 t_enxsubblock
*s1
, *s2
;
1295 cmp_int(stdout
, buf
, -1, (int)s1
->type
, (int)s2
->type
);
1296 cmp_int64(stdout
, buf
, s1
->nr
, s2
->nr
);
1298 if ((s1
->type
== s2
->type
) && (s1
->nr
== s2
->nr
))
1302 case xdr_datatype_float
:
1303 for (k
= 0; k
< s1
->nr
; k
++)
1305 cmp_float(stdout
, buf
, i
,
1306 s1
->fval
[k
], s2
->fval
[k
],
1310 case xdr_datatype_double
:
1311 for (k
= 0; k
< s1
->nr
; k
++)
1313 cmp_double(stdout
, buf
, i
,
1314 s1
->dval
[k
], s2
->dval
[k
],
1318 case xdr_datatype_int
:
1319 for (k
= 0; k
< s1
->nr
; k
++)
1321 cmp_int(stdout
, buf
, i
,
1322 s1
->ival
[k
], s2
->ival
[k
]);
1325 case xdr_datatype_int64
:
1326 for (k
= 0; k
< s1
->nr
; k
++)
1328 cmp_int64(stdout
, buf
,
1329 s1
->lval
[k
], s2
->lval
[k
]);
1332 case xdr_datatype_char
:
1333 for (k
= 0; k
< s1
->nr
; k
++)
1335 cmp_uc(stdout
, buf
, i
,
1336 s1
->cval
[k
], s2
->cval
[k
]);
1339 case xdr_datatype_string
:
1340 for (k
= 0; k
< s1
->nr
; k
++)
1342 cmp_str(stdout
, buf
, i
,
1343 s1
->sval
[k
], s2
->sval
[k
]);
1347 gmx_incons("Unknown data type!!");
1356 void comp_enx(const char *fn1
, const char *fn2
, real ftol
, real abstol
, const char *lastener
)
1358 int nre
, nre1
, nre2
;
1359 ener_file_t in1
, in2
;
1360 int i
, j
, maxener
, *ind1
, *ind2
, *have
;
1361 gmx_enxnm_t
*enm1
= NULL
, *enm2
= NULL
;
1362 t_enxframe
*fr1
, *fr2
;
1365 fprintf(stdout
, "comparing energy file %s and %s\n\n", fn1
, fn2
);
1367 in1
= open_enx(fn1
, "r");
1368 in2
= open_enx(fn2
, "r");
1369 do_enxnms(in1
, &nre1
, &enm1
);
1370 do_enxnms(in2
, &nre2
, &enm2
);
1373 fprintf(stdout
, "There are %d and %d terms in the energy files\n\n",
1378 fprintf(stdout
, "There are %d terms in the energy files\n\n", nre1
);
1385 for (i
= 0; i
< nre1
; i
++)
1387 for (j
= 0; j
< nre2
; j
++)
1389 if (enernm_equal(enm1
[i
].name
, enm2
[j
].name
))
1398 if (nre
== 0 || ind1
[nre
-1] != i
)
1400 cmp_str(stdout
, "enm", i
, enm1
[i
].name
, "-");
1403 for (i
= 0; i
< nre2
; i
++)
1407 cmp_str(stdout
, "enm", i
, "-", enm2
[i
].name
);
1412 for (i
= 0; i
< nre
; i
++)
1414 if ((lastener
!= NULL
) && (std::strstr(enm1
[i
].name
, lastener
) != NULL
))
1421 fprintf(stdout
, "There are %d terms to compare in the energy files\n\n",
1424 for (i
= 0; i
< maxener
; i
++)
1426 cmp_str(stdout
, "unit", i
, enm1
[ind1
[i
]].unit
, enm2
[ind2
[i
]].unit
);
1433 b1
= do_enx(in1
, fr1
);
1434 b2
= do_enx(in2
, fr2
);
1437 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn2
, fn1
);
1441 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn1
, fn2
);
1443 else if (!b1
&& !b2
)
1445 fprintf(stdout
, "\nFiles read successfully\n");
1449 cmp_real(stdout
, "t", -1, fr1
->t
, fr2
->t
, ftol
, abstol
);
1450 cmp_int(stdout
, "step", -1, fr1
->step
, fr2
->step
);
1451 /* We don't want to print the nre mismatch for every frame */
1452 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1453 if ((fr1
->nre
>= nre
) && (fr2
->nre
>= nre
))
1455 cmp_energies(stdout
, fr1
->step
, fr1
->step
, fr1
->ener
, fr2
->ener
,
1456 enm1
, ftol
, abstol
, nre
, ind1
, ind2
, maxener
);
1458 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1459 cmp_eblocks(fr1
, fr2
, ftol
, abstol
);