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/mdtypes/inputrec.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/trajectory/trajectoryframe.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
61 #include "gromacs/utility/stringutil.h"
63 static void cmp_int(FILE *fp
, const char *s
, int index
, int i1
, int i2
)
69 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
73 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
78 static void cmp_int64(FILE *fp
, const char *s
, gmx_int64_t i1
, gmx_int64_t i2
)
82 fprintf(fp
, "%s (", s
);
83 fprintf(fp
, "%" GMX_PRId64
, i1
);
85 fprintf(fp
, "%" GMX_PRId64
, i2
);
90 static void cmp_us(FILE *fp
, const char *s
, int index
, unsigned short i1
, unsigned short i2
)
96 fprintf(fp
, "%s[%d] (%hu - %hu)\n", s
, index
, i1
, i2
);
100 fprintf(fp
, "%s (%hu - %hu)\n", s
, i1
, i2
);
105 static void cmp_uc(FILE *fp
, const char *s
, int index
, unsigned char i1
, unsigned char i2
)
111 fprintf(fp
, "%s[%d] (%d - %d)\n", s
, index
, i1
, i2
);
115 fprintf(fp
, "%s (%d - %d)\n", s
, i1
, i2
);
120 static gmx_bool
cmp_bool(FILE *fp
, const char *s
, int index
, gmx_bool b1
, gmx_bool b2
)
142 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
,
143 gmx::boolToString(b1
), gmx::boolToString(b2
));
147 fprintf(fp
, "%s (%s - %s)\n", s
,
148 gmx::boolToString(b1
), gmx::boolToString(b2
));
154 static void cmp_str(FILE *fp
, const char *s
, int index
,
155 const char *s1
, const char *s2
)
157 if (std::strcmp(s1
, s2
) != 0)
161 fprintf(fp
, "%s[%d] (%s - %s)\n", s
, index
, s1
, s2
);
165 fprintf(fp
, "%s (%s - %s)\n", s
, s1
, s2
);
170 static gmx_bool
equal_real(real i1
, real i2
, real ftol
, real abstol
)
172 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
175 static gmx_bool
equal_float(float i1
, float i2
, float ftol
, float abstol
)
177 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
180 static gmx_bool
equal_double(double i1
, double i2
, real ftol
, real abstol
)
182 return ( ( 2*fabs(i1
- i2
) <= (fabs(i1
) + fabs(i2
))*ftol
) || fabs(i1
-i2
) <= abstol
);
186 cmp_real(FILE *fp
, const char *s
, int index
, real i1
, real i2
, real ftol
, real abstol
)
188 if (!equal_real(i1
, i2
, ftol
, abstol
))
192 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
196 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
202 cmp_float(FILE *fp
, const char *s
, int index
, float i1
, float i2
, float ftol
, float abstol
)
204 if (!equal_float(i1
, i2
, ftol
, abstol
))
208 fprintf(fp
, "%s[%2d] (%e - %e)\n", s
, index
, i1
, i2
);
212 fprintf(fp
, "%s (%e - %e)\n", s
, i1
, i2
);
220 cmp_double(FILE *fp
, const char *s
, int index
, double i1
, double i2
, double ftol
, double abstol
)
222 if (!equal_double(i1
, i2
, ftol
, abstol
))
226 fprintf(fp
, "%s[%2d] (%16.9e - %16.9e)\n", s
, index
, i1
, i2
);
230 fprintf(fp
, "%s (%16.9e - %16.9e)\n", s
, i1
, i2
);
235 static void cmp_rvec(FILE *fp
, const char *s
, int index
, rvec i1
, rvec i2
, real ftol
, real abstol
)
237 if (!equal_real(i1
[XX
], i2
[XX
], ftol
, abstol
) ||
238 !equal_real(i1
[YY
], i2
[YY
], ftol
, abstol
) ||
239 !equal_real(i1
[ZZ
], i2
[ZZ
], ftol
, abstol
))
243 fprintf(fp
, "%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
244 s
, index
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
248 fprintf(fp
, "%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
249 s
, i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
254 static void cmp_ivec(FILE *fp
, const char *s
, int index
, ivec i1
, ivec i2
)
256 if ((i1
[XX
] != i2
[XX
]) || (i1
[YY
] != i2
[YY
]) || (i1
[ZZ
] != i2
[ZZ
]))
260 fprintf(fp
, "%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
, index
,
261 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
265 fprintf(fp
, "%s (%8d,%8d,%8d - %8d,%8d,%8d)\n", s
,
266 i1
[XX
], i1
[YY
], i1
[ZZ
], i2
[XX
], i2
[YY
], i2
[ZZ
]);
271 static void cmp_ilist(FILE *fp
, int ftype
, t_ilist
*il1
, t_ilist
*il2
)
276 fprintf(fp
, "comparing ilist %s\n", interaction_function
[ftype
].name
);
277 sprintf(buf
, "%s->nr", interaction_function
[ftype
].name
);
278 cmp_int(fp
, buf
, -1, il1
->nr
, il2
->nr
);
279 sprintf(buf
, "%s->iatoms", interaction_function
[ftype
].name
);
280 if (((il1
->nr
> 0) && (!il1
->iatoms
)) ||
281 ((il2
->nr
> 0) && (!il2
->iatoms
)) ||
282 ((il1
->nr
!= il2
->nr
)))
284 fprintf(fp
, "Comparing radically different topologies - %s is different\n",
289 for (i
= 0; (i
< il1
->nr
); i
++)
291 cmp_int(fp
, buf
, i
, il1
->iatoms
[i
], il2
->iatoms
[i
]);
296 void cmp_iparm(FILE *fp
, const char *s
, t_functype ft
,
297 t_iparams ip1
, t_iparams ip2
, real ftol
, real abstol
)
303 for (i
= 0; i
< MAXFORCEPARAM
&& !bDiff
; i
++)
305 bDiff
= !equal_real(ip1
.generic
.buf
[i
], ip2
.generic
.buf
[i
], ftol
, abstol
);
309 fprintf(fp
, "%s1: ", s
);
310 pr_iparams(fp
, ft
, &ip1
);
311 fprintf(fp
, "%s2: ", s
);
312 pr_iparams(fp
, ft
, &ip2
);
316 void cmp_iparm_AB(FILE *fp
, const char *s
, t_functype ft
, t_iparams ip1
, real ftol
, real abstol
)
318 int nrfpA
, nrfpB
, p0
, i
;
321 /* Normally the first parameter is perturbable */
323 nrfpA
= interaction_function
[ft
].nrfpA
;
324 nrfpB
= interaction_function
[ft
].nrfpB
;
329 else if (interaction_function
[ft
].flags
& IF_TABULATED
)
331 /* For tabulated interactions only the second parameter is perturbable */
336 for (i
= 0; i
< nrfpB
&& !bDiff
; i
++)
338 bDiff
= !equal_real(ip1
.generic
.buf
[p0
+i
], ip1
.generic
.buf
[nrfpA
+i
], ftol
, abstol
);
342 fprintf(fp
, "%s: ", s
);
343 pr_iparams(fp
, ft
, &ip1
);
347 static void cmp_cmap(FILE *fp
, const gmx_cmap_t
*cmap1
, const gmx_cmap_t
*cmap2
, real ftol
, real abstol
)
349 cmp_int(fp
, "cmap ngrid", -1, cmap1
->ngrid
, cmap2
->ngrid
);
350 cmp_int(fp
, "cmap grid_spacing", -1, cmap1
->grid_spacing
, cmap2
->grid_spacing
);
351 if (cmap1
->ngrid
== cmap2
->ngrid
&&
352 cmap1
->grid_spacing
== cmap2
->grid_spacing
)
356 for (g
= 0; g
< cmap1
->ngrid
; g
++)
360 fprintf(fp
, "comparing cmap %d\n", g
);
362 for (i
= 0; i
< 4*cmap1
->grid_spacing
*cmap1
->grid_spacing
; i
++)
364 cmp_real(fp
, "", i
, cmap1
->cmapdata
[g
].cmap
[i
], cmap2
->cmapdata
[g
].cmap
[i
], ftol
, abstol
);
370 static void cmp_idef(FILE *fp
, t_idef
*id1
, t_idef
*id2
, real ftol
, real abstol
)
373 char buf1
[64], buf2
[64];
375 fprintf(fp
, "comparing idef\n");
378 cmp_int(fp
, "idef->ntypes", -1, id1
->ntypes
, id2
->ntypes
);
379 cmp_int(fp
, "idef->atnr", -1, id1
->atnr
, id2
->atnr
);
380 for (i
= 0; (i
< std::min(id1
->ntypes
, id2
->ntypes
)); i
++)
382 sprintf(buf1
, "idef->functype[%d]", i
);
383 sprintf(buf2
, "idef->iparam[%d]", i
);
384 cmp_int(fp
, buf1
, i
, (int)id1
->functype
[i
], (int)id2
->functype
[i
]);
385 cmp_iparm(fp
, buf2
, id1
->functype
[i
],
386 id1
->iparams
[i
], id2
->iparams
[i
], ftol
, abstol
);
388 cmp_real(fp
, "fudgeQQ", -1, id1
->fudgeQQ
, id2
->fudgeQQ
, ftol
, abstol
);
389 cmp_cmap(fp
, &id1
->cmap_grid
, &id2
->cmap_grid
, ftol
, abstol
);
390 for (i
= 0; (i
< F_NRE
); i
++)
392 cmp_ilist(fp
, i
, &(id1
->il
[i
]), &(id2
->il
[i
]));
397 for (i
= 0; (i
< id1
->ntypes
); i
++)
399 cmp_iparm_AB(fp
, "idef->iparam", id1
->functype
[i
], id1
->iparams
[i
], ftol
, abstol
);
404 static void cmp_block(FILE *fp
, t_block
*b1
, t_block
*b2
, const char *s
)
408 fprintf(fp
, "comparing block %s\n", s
);
409 sprintf(buf
, "%s.nr", s
);
410 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
413 static void cmp_blocka(FILE *fp
, t_blocka
*b1
, t_blocka
*b2
, const char *s
)
417 fprintf(fp
, "comparing blocka %s\n", s
);
418 sprintf(buf
, "%s.nr", s
);
419 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
420 sprintf(buf
, "%s.nra", s
);
421 cmp_int(fp
, buf
, -1, b1
->nra
, b2
->nra
);
424 static void cmp_atom(FILE *fp
, int index
, t_atom
*a1
, t_atom
*a2
, real ftol
, real abstol
)
428 cmp_us(fp
, "atom.type", index
, a1
->type
, a2
->type
);
429 cmp_us(fp
, "atom.ptype", index
, a1
->ptype
, a2
->ptype
);
430 cmp_int(fp
, "atom.resind", index
, a1
->resind
, a2
->resind
);
431 cmp_int(fp
, "atom.atomnumber", index
, a1
->atomnumber
, a2
->atomnumber
);
432 cmp_real(fp
, "atom.m", index
, a1
->m
, a2
->m
, ftol
, abstol
);
433 cmp_real(fp
, "atom.q", index
, a1
->q
, a2
->q
, ftol
, abstol
);
434 cmp_us(fp
, "atom.typeB", index
, a1
->typeB
, a2
->typeB
);
435 cmp_real(fp
, "atom.mB", index
, a1
->mB
, a2
->mB
, ftol
, abstol
);
436 cmp_real(fp
, "atom.qB", index
, a1
->qB
, a2
->qB
, ftol
, abstol
);
440 cmp_us(fp
, "atom.type", index
, a1
->type
, a1
->typeB
);
441 cmp_real(fp
, "atom.m", index
, a1
->m
, a1
->mB
, ftol
, abstol
);
442 cmp_real(fp
, "atom.q", index
, a1
->q
, a1
->qB
, ftol
, abstol
);
446 static void cmp_atoms(FILE *fp
, t_atoms
*a1
, t_atoms
*a2
, real ftol
, real abstol
)
450 fprintf(fp
, "comparing atoms\n");
454 cmp_int(fp
, "atoms->nr", -1, a1
->nr
, a2
->nr
);
455 for (i
= 0; (i
< a1
->nr
); i
++)
457 cmp_atom(fp
, i
, &(a1
->atom
[i
]), &(a2
->atom
[i
]), ftol
, abstol
);
462 for (i
= 0; (i
< a1
->nr
); i
++)
464 cmp_atom(fp
, i
, &(a1
->atom
[i
]), NULL
, ftol
, abstol
);
469 static void cmp_top(FILE *fp
, t_topology
*t1
, t_topology
*t2
, real ftol
, real abstol
)
471 fprintf(fp
, "comparing top\n");
474 cmp_idef(fp
, &(t1
->idef
), &(t2
->idef
), ftol
, abstol
);
475 cmp_atoms(fp
, &(t1
->atoms
), &(t2
->atoms
), ftol
, abstol
);
476 cmp_block(fp
, &t1
->cgs
, &t2
->cgs
, "cgs");
477 cmp_block(fp
, &t1
->mols
, &t2
->mols
, "mols");
478 cmp_bool(fp
, "bIntermolecularInteractions", -1, t1
->bIntermolecularInteractions
, t2
->bIntermolecularInteractions
);
479 cmp_blocka(fp
, &t1
->excls
, &t2
->excls
, "excls");
483 cmp_idef(fp
, &(t1
->idef
), NULL
, ftol
, abstol
);
484 cmp_atoms(fp
, &(t1
->atoms
), NULL
, ftol
, abstol
);
488 static void cmp_groups(FILE *fp
, gmx_groups_t
*g0
, gmx_groups_t
*g1
,
489 int natoms0
, int natoms1
)
494 fprintf(fp
, "comparing groups\n");
496 for (i
= 0; i
< egcNR
; i
++)
498 sprintf(buf
, "grps[%d].nr", i
);
499 cmp_int(fp
, buf
, -1, g0
->grps
[i
].nr
, g1
->grps
[i
].nr
);
500 if (g0
->grps
[i
].nr
== g1
->grps
[i
].nr
)
502 for (j
= 0; j
< g0
->grps
[i
].nr
; j
++)
504 sprintf(buf
, "grps[%d].name[%d]", i
, j
);
506 *g0
->grpname
[g0
->grps
[i
].nm_ind
[j
]],
507 *g1
->grpname
[g1
->grps
[i
].nm_ind
[j
]]);
510 cmp_int(fp
, "ngrpnr", i
, g0
->ngrpnr
[i
], g1
->ngrpnr
[i
]);
511 if (g0
->ngrpnr
[i
] == g1
->ngrpnr
[i
] && natoms0
== natoms1
&&
512 (g0
->grpnr
[i
] != NULL
|| g1
->grpnr
[i
] != NULL
))
514 for (j
= 0; j
< natoms0
; j
++)
516 cmp_int(fp
, gtypes
[i
], j
, ggrpnr(g0
, i
, j
), ggrpnr(g1
, i
, j
));
520 /* We have compared the names in the groups lists,
521 * so we can skip the grpname list comparison.
525 static void cmp_rvecs_rmstol(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
526 real ftol
, real abstol
)
531 /* For a vector you are usally not interested in a relative difference
532 * on a component that is very small compared to the other components.
533 * Therefore we do the relative comparision relative to the RMS component.
536 for (i
= 0; (i
< n
); i
++)
538 for (m
= 0; m
< DIM
; m
++)
540 rms
+= x1
[i
][m
]*x1
[i
][m
] + x2
[i
][m
]*x2
[i
][m
];
543 rms
= sqrt(rms
/(2*n
*DIM
));
545 /* Convert the relative tolerance into an absolute tolerance */
546 if (ftol
*rms
< abstol
)
551 /* And now do the actual comparision */
552 for (i
= 0; (i
< n
); i
++)
554 cmp_rvec(fp
, title
, i
, x1
[i
], x2
[i
], 0.0, abstol
);
558 static void cmp_rvecs(FILE *fp
, const char *title
, int n
, rvec x1
[], rvec x2
[],
559 gmx_bool bRMSD
, real ftol
, real abstol
)
567 for (i
= 0; (i
< n
); i
++)
569 for (m
= 0; m
< DIM
; m
++)
571 d
= x1
[i
][m
] - x2
[i
][m
];
575 fprintf(fp
, "%s RMSD %g\n", title
, std::sqrt(ssd
/n
));
579 cmp_rvecs_rmstol(fp
, title
, n
, x1
, x2
, ftol
, abstol
);
583 static void cmp_grpopts(FILE *fp
, t_grpopts
*opt1
, t_grpopts
*opt2
, real ftol
, real abstol
)
586 char buf1
[256], buf2
[256];
588 cmp_int(fp
, "inputrec->grpopts.ngtc", -1, opt1
->ngtc
, opt2
->ngtc
);
589 cmp_int(fp
, "inputrec->grpopts.ngacc", -1, opt1
->ngacc
, opt2
->ngacc
);
590 cmp_int(fp
, "inputrec->grpopts.ngfrz", -1, opt1
->ngfrz
, opt2
->ngfrz
);
591 cmp_int(fp
, "inputrec->grpopts.ngener", -1, opt1
->ngener
, opt2
->ngener
);
592 for (i
= 0; (i
< std::min(opt1
->ngtc
, opt2
->ngtc
)); i
++)
594 cmp_real(fp
, "inputrec->grpopts.nrdf", i
, opt1
->nrdf
[i
], opt2
->nrdf
[i
], ftol
, abstol
);
595 cmp_real(fp
, "inputrec->grpopts.ref_t", i
, opt1
->ref_t
[i
], opt2
->ref_t
[i
], ftol
, abstol
);
596 cmp_real(fp
, "inputrec->grpopts.tau_t", i
, opt1
->tau_t
[i
], opt2
->tau_t
[i
], ftol
, abstol
);
597 cmp_int(fp
, "inputrec->grpopts.annealing", i
, opt1
->annealing
[i
], opt2
->annealing
[i
]);
598 cmp_int(fp
, "inputrec->grpopts.anneal_npoints", i
,
599 opt1
->anneal_npoints
[i
], opt2
->anneal_npoints
[i
]);
600 if (opt1
->anneal_npoints
[i
] == opt2
->anneal_npoints
[i
])
602 sprintf(buf1
, "inputrec->grpopts.anneal_time[%d]", i
);
603 sprintf(buf2
, "inputrec->grpopts.anneal_temp[%d]", i
);
604 for (j
= 0; j
< opt1
->anneal_npoints
[i
]; j
++)
606 cmp_real(fp
, buf1
, j
, opt1
->anneal_time
[i
][j
], opt2
->anneal_time
[i
][j
], ftol
, abstol
);
607 cmp_real(fp
, buf2
, j
, opt1
->anneal_temp
[i
][j
], opt2
->anneal_temp
[i
][j
], ftol
, abstol
);
611 if (opt1
->ngener
== opt2
->ngener
)
613 for (i
= 0; i
< opt1
->ngener
; i
++)
615 for (j
= i
; j
< opt1
->ngener
; j
++)
617 sprintf(buf1
, "inputrec->grpopts.egp_flags[%d]", i
);
619 opt1
->egp_flags
[opt1
->ngener
*i
+j
],
620 opt2
->egp_flags
[opt1
->ngener
*i
+j
]);
624 for (i
= 0; (i
< std::min(opt1
->ngacc
, opt2
->ngacc
)); i
++)
626 cmp_rvec(fp
, "inputrec->grpopts.acc", i
, opt1
->acc
[i
], opt2
->acc
[i
], ftol
, abstol
);
628 for (i
= 0; (i
< std::min(opt1
->ngfrz
, opt2
->ngfrz
)); i
++)
630 cmp_ivec(fp
, "inputrec->grpopts.nFreeze", i
, opt1
->nFreeze
[i
], opt2
->nFreeze
[i
]);
634 static void cmp_cosines(FILE *fp
, const char *s
, t_cosines c1
[DIM
], t_cosines c2
[DIM
], real ftol
, real abstol
)
639 for (m
= 0; (m
< DIM
); m
++)
641 sprintf(buf
, "inputrec->%s[%d]", s
, m
);
642 cmp_int(fp
, buf
, 0, c1
->n
, c2
->n
);
643 for (i
= 0; (i
< std::min(c1
->n
, c2
->n
)); i
++)
645 cmp_real(fp
, buf
, i
, c1
->a
[i
], c2
->a
[i
], ftol
, abstol
);
646 cmp_real(fp
, buf
, i
, c1
->phi
[i
], c2
->phi
[i
], ftol
, abstol
);
650 static void cmp_pull(FILE *fp
)
652 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");
655 static void cmp_simtempvals(FILE *fp
, t_simtemp
*simtemp1
, t_simtemp
*simtemp2
, int n_lambda
, real ftol
, real abstol
)
658 cmp_int(fp
, "inputrec->simtempvals->eSimTempScale", -1, simtemp1
->eSimTempScale
, simtemp2
->eSimTempScale
);
659 cmp_real(fp
, "inputrec->simtempvals->simtemp_high", -1, simtemp1
->simtemp_high
, simtemp2
->simtemp_high
, ftol
, abstol
);
660 cmp_real(fp
, "inputrec->simtempvals->simtemp_low", -1, simtemp1
->simtemp_low
, simtemp2
->simtemp_low
, ftol
, abstol
);
661 for (i
= 0; i
< n_lambda
; i
++)
663 cmp_real(fp
, "inputrec->simtempvals->temperatures", -1, simtemp1
->temperatures
[i
], simtemp2
->temperatures
[i
], ftol
, abstol
);
667 static void cmp_expandedvals(FILE *fp
, t_expanded
*expand1
, t_expanded
*expand2
, int n_lambda
, real ftol
, real abstol
)
671 cmp_bool(fp
, "inputrec->fepvals->bInit_weights", -1, expand1
->bInit_weights
, expand2
->bInit_weights
);
672 cmp_bool(fp
, "inputrec->fepvals->bWLoneovert", -1, expand1
->bWLoneovert
, expand2
->bWLoneovert
);
674 for (i
= 0; i
< n_lambda
; i
++)
676 cmp_real(fp
, "inputrec->expandedvals->init_lambda_weights", -1,
677 expand1
->init_lambda_weights
[i
], expand2
->init_lambda_weights
[i
], ftol
, abstol
);
680 cmp_int(fp
, "inputrec->expandedvals->lambda-stats", -1, expand1
->elamstats
, expand2
->elamstats
);
681 cmp_int(fp
, "inputrec->expandedvals->lambda-mc-move", -1, expand1
->elmcmove
, expand2
->elmcmove
);
682 cmp_int(fp
, "inputrec->expandedvals->lmc-repeats", -1, expand1
->lmc_repeats
, expand2
->lmc_repeats
);
683 cmp_int(fp
, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1
->gibbsdeltalam
, expand2
->gibbsdeltalam
);
684 cmp_int(fp
, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1
->lmc_forced_nstart
, expand2
->lmc_forced_nstart
);
685 cmp_int(fp
, "inputrec->expandedvals->lambda-weights-equil", -1, expand1
->elmceq
, expand2
->elmceq
);
686 cmp_int(fp
, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1
->equil_n_at_lam
, expand2
->equil_n_at_lam
);
687 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1
->equil_samples
, expand2
->equil_samples
);
688 cmp_int(fp
, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1
->equil_steps
, expand2
->equil_steps
);
689 cmp_real(fp
, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1
->equil_wl_delta
, expand2
->equil_wl_delta
, ftol
, abstol
);
690 cmp_real(fp
, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1
->equil_ratio
, expand2
->equil_ratio
, ftol
, abstol
);
691 cmp_bool(fp
, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1
->bSymmetrizedTMatrix
, expand2
->bSymmetrizedTMatrix
);
692 cmp_int(fp
, "inputrec->expandedvals->nstTij", -1, expand1
->nstTij
, expand2
->nstTij
);
693 cmp_int(fp
, "inputrec->expandedvals->mininum-var-min", -1, expand1
->minvarmin
, expand2
->minvarmin
); /*default is reasonable */
694 cmp_int(fp
, "inputrec->expandedvals->weight-c-range", -1, expand1
->c_range
, expand2
->c_range
); /* default is just C=0 */
695 cmp_real(fp
, "inputrec->expandedvals->wl-scale", -1, expand1
->wl_scale
, expand2
->wl_scale
, ftol
, abstol
);
696 cmp_real(fp
, "inputrec->expandedvals->init-wl-delta", -1, expand1
->init_wl_delta
, expand2
->init_wl_delta
, ftol
, abstol
);
697 cmp_real(fp
, "inputrec->expandedvals->wl-ratio", -1, expand1
->wl_ratio
, expand2
->wl_ratio
, ftol
, abstol
);
698 cmp_int(fp
, "inputrec->expandedvals->nstexpanded", -1, expand1
->nstexpanded
, expand2
->nstexpanded
);
699 cmp_int(fp
, "inputrec->expandedvals->lmc-seed", -1, expand1
->lmc_seed
, expand2
->lmc_seed
);
700 cmp_real(fp
, "inputrec->expandedvals->mc-temperature", -1, expand1
->mc_temp
, expand2
->mc_temp
, ftol
, abstol
);
703 static void cmp_fepvals(FILE *fp
, t_lambda
*fep1
, t_lambda
*fep2
, real ftol
, real abstol
)
706 cmp_int(fp
, "inputrec->nstdhdl", -1, fep1
->nstdhdl
, fep2
->nstdhdl
);
707 cmp_double(fp
, "inputrec->fepvals->init_fep_state", -1, fep1
->init_fep_state
, fep2
->init_fep_state
, ftol
, abstol
);
708 cmp_double(fp
, "inputrec->fepvals->delta_lambda", -1, fep1
->delta_lambda
, fep2
->delta_lambda
, ftol
, abstol
);
709 cmp_int(fp
, "inputrec->fepvals->n_lambda", -1, fep1
->n_lambda
, fep2
->n_lambda
);
710 for (i
= 0; i
< efptNR
; i
++)
712 for (j
= 0; j
< std::min(fep1
->n_lambda
, fep2
->n_lambda
); j
++)
714 cmp_double(fp
, "inputrec->fepvals->all_lambda", -1, fep1
->all_lambda
[i
][j
], fep2
->all_lambda
[i
][j
], ftol
, abstol
);
717 cmp_int(fp
, "inputrec->fepvals->lambda_neighbors", 1, fep1
->lambda_neighbors
,
718 fep2
->lambda_neighbors
);
719 cmp_real(fp
, "inputrec->fepvals->sc_alpha", -1, fep1
->sc_alpha
, fep2
->sc_alpha
, ftol
, abstol
);
720 cmp_int(fp
, "inputrec->fepvals->sc_power", -1, fep1
->sc_power
, fep2
->sc_power
);
721 cmp_real(fp
, "inputrec->fepvals->sc_r_power", -1, fep1
->sc_r_power
, fep2
->sc_r_power
, ftol
, abstol
);
722 cmp_real(fp
, "inputrec->fepvals->sc_sigma", -1, fep1
->sc_sigma
, fep2
->sc_sigma
, ftol
, abstol
);
723 cmp_int(fp
, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1
->edHdLPrintEnergy
, fep1
->edHdLPrintEnergy
);
724 cmp_bool(fp
, "inputrec->fepvals->bScCoul", -1, fep1
->bScCoul
, fep1
->bScCoul
);
725 cmp_int(fp
, "inputrec->separate_dhdl_file", -1, fep1
->separate_dhdl_file
, fep2
->separate_dhdl_file
);
726 cmp_int(fp
, "inputrec->dhdl_derivatives", -1, fep1
->dhdl_derivatives
, fep2
->dhdl_derivatives
);
727 cmp_int(fp
, "inputrec->dh_hist_size", -1, fep1
->dh_hist_size
, fep2
->dh_hist_size
);
728 cmp_double(fp
, "inputrec->dh_hist_spacing", -1, fep1
->dh_hist_spacing
, fep2
->dh_hist_spacing
, ftol
, abstol
);
731 static void cmp_inputrec(FILE *fp
, t_inputrec
*ir1
, t_inputrec
*ir2
, real ftol
, real abstol
)
733 fprintf(fp
, "comparing inputrec\n");
735 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
736 * of warnings. Maybe it will change in future versions, but for the
737 * moment I've spelled them out instead. /EL 000820
738 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
739 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
740 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
742 cmp_int(fp
, "inputrec->eI", -1, ir1
->eI
, ir2
->eI
);
743 cmp_int64(fp
, "inputrec->nsteps", ir1
->nsteps
, ir2
->nsteps
);
744 cmp_int64(fp
, "inputrec->init_step", ir1
->init_step
, ir2
->init_step
);
745 cmp_int(fp
, "inputrec->simulation_part", -1, ir1
->simulation_part
, ir2
->simulation_part
);
746 cmp_int(fp
, "inputrec->ePBC", -1, ir1
->ePBC
, ir2
->ePBC
);
747 cmp_int(fp
, "inputrec->bPeriodicMols", -1, ir1
->bPeriodicMols
, ir2
->bPeriodicMols
);
748 cmp_int(fp
, "inputrec->cutoff_scheme", -1, ir1
->cutoff_scheme
, ir2
->cutoff_scheme
);
749 cmp_int(fp
, "inputrec->ns_type", -1, ir1
->ns_type
, ir2
->ns_type
);
750 cmp_int(fp
, "inputrec->nstlist", -1, ir1
->nstlist
, ir2
->nstlist
);
751 cmp_int(fp
, "inputrec->nstcomm", -1, ir1
->nstcomm
, ir2
->nstcomm
);
752 cmp_int(fp
, "inputrec->comm_mode", -1, ir1
->comm_mode
, ir2
->comm_mode
);
753 cmp_int(fp
, "inputrec->nstlog", -1, ir1
->nstlog
, ir2
->nstlog
);
754 cmp_int(fp
, "inputrec->nstxout", -1, ir1
->nstxout
, ir2
->nstxout
);
755 cmp_int(fp
, "inputrec->nstvout", -1, ir1
->nstvout
, ir2
->nstvout
);
756 cmp_int(fp
, "inputrec->nstfout", -1, ir1
->nstfout
, ir2
->nstfout
);
757 cmp_int(fp
, "inputrec->nstcalcenergy", -1, ir1
->nstcalcenergy
, ir2
->nstcalcenergy
);
758 cmp_int(fp
, "inputrec->nstenergy", -1, ir1
->nstenergy
, ir2
->nstenergy
);
759 cmp_int(fp
, "inputrec->nstxout_compressed", -1, ir1
->nstxout_compressed
, ir2
->nstxout_compressed
);
760 cmp_double(fp
, "inputrec->init_t", -1, ir1
->init_t
, ir2
->init_t
, ftol
, abstol
);
761 cmp_double(fp
, "inputrec->delta_t", -1, ir1
->delta_t
, ir2
->delta_t
, ftol
, abstol
);
762 cmp_real(fp
, "inputrec->x_compression_precision", -1, ir1
->x_compression_precision
, ir2
->x_compression_precision
, ftol
, abstol
);
763 cmp_real(fp
, "inputrec->fourierspacing", -1, ir1
->fourier_spacing
, ir2
->fourier_spacing
, ftol
, abstol
);
764 cmp_int(fp
, "inputrec->nkx", -1, ir1
->nkx
, ir2
->nkx
);
765 cmp_int(fp
, "inputrec->nky", -1, ir1
->nky
, ir2
->nky
);
766 cmp_int(fp
, "inputrec->nkz", -1, ir1
->nkz
, ir2
->nkz
);
767 cmp_int(fp
, "inputrec->pme_order", -1, ir1
->pme_order
, ir2
->pme_order
);
768 cmp_real(fp
, "inputrec->ewald_rtol", -1, ir1
->ewald_rtol
, ir2
->ewald_rtol
, ftol
, abstol
);
769 cmp_int(fp
, "inputrec->ewald_geometry", -1, ir1
->ewald_geometry
, ir2
->ewald_geometry
);
770 cmp_real(fp
, "inputrec->epsilon_surface", -1, ir1
->epsilon_surface
, ir2
->epsilon_surface
, ftol
, abstol
);
771 cmp_int(fp
, "inputrec->bContinuation", -1, ir1
->bContinuation
, ir2
->bContinuation
);
772 cmp_int(fp
, "inputrec->bShakeSOR", -1, ir1
->bShakeSOR
, ir2
->bShakeSOR
);
773 cmp_int(fp
, "inputrec->etc", -1, ir1
->etc
, ir2
->etc
);
774 cmp_int(fp
, "inputrec->bPrintNHChains", -1, ir1
->bPrintNHChains
, ir2
->bPrintNHChains
);
775 cmp_int(fp
, "inputrec->epc", -1, ir1
->epc
, ir2
->epc
);
776 cmp_int(fp
, "inputrec->epct", -1, ir1
->epct
, ir2
->epct
);
777 cmp_real(fp
, "inputrec->tau_p", -1, ir1
->tau_p
, ir2
->tau_p
, ftol
, abstol
);
778 cmp_rvec(fp
, "inputrec->ref_p(x)", -1, ir1
->ref_p
[XX
], ir2
->ref_p
[XX
], ftol
, abstol
);
779 cmp_rvec(fp
, "inputrec->ref_p(y)", -1, ir1
->ref_p
[YY
], ir2
->ref_p
[YY
], ftol
, abstol
);
780 cmp_rvec(fp
, "inputrec->ref_p(z)", -1, ir1
->ref_p
[ZZ
], ir2
->ref_p
[ZZ
], ftol
, abstol
);
781 cmp_rvec(fp
, "inputrec->compress(x)", -1, ir1
->compress
[XX
], ir2
->compress
[XX
], ftol
, abstol
);
782 cmp_rvec(fp
, "inputrec->compress(y)", -1, ir1
->compress
[YY
], ir2
->compress
[YY
], ftol
, abstol
);
783 cmp_rvec(fp
, "inputrec->compress(z)", -1, ir1
->compress
[ZZ
], ir2
->compress
[ZZ
], ftol
, abstol
);
784 cmp_int(fp
, "refcoord_scaling", -1, ir1
->refcoord_scaling
, ir2
->refcoord_scaling
);
785 cmp_rvec(fp
, "inputrec->posres_com", -1, ir1
->posres_com
, ir2
->posres_com
, ftol
, abstol
);
786 cmp_rvec(fp
, "inputrec->posres_comB", -1, ir1
->posres_comB
, ir2
->posres_comB
, ftol
, abstol
);
787 cmp_real(fp
, "inputrec->verletbuf_tol", -1, ir1
->verletbuf_tol
, ir2
->verletbuf_tol
, ftol
, abstol
);
788 cmp_real(fp
, "inputrec->rlist", -1, ir1
->rlist
, ir2
->rlist
, ftol
, abstol
);
789 cmp_real(fp
, "inputrec->rtpi", -1, ir1
->rtpi
, ir2
->rtpi
, ftol
, abstol
);
790 cmp_int(fp
, "inputrec->coulombtype", -1, ir1
->coulombtype
, ir2
->coulombtype
);
791 cmp_int(fp
, "inputrec->coulomb_modifier", -1, ir1
->coulomb_modifier
, ir2
->coulomb_modifier
);
792 cmp_real(fp
, "inputrec->rcoulomb_switch", -1, ir1
->rcoulomb_switch
, ir2
->rcoulomb_switch
, ftol
, abstol
);
793 cmp_real(fp
, "inputrec->rcoulomb", -1, ir1
->rcoulomb
, ir2
->rcoulomb
, ftol
, abstol
);
794 cmp_int(fp
, "inputrec->vdwtype", -1, ir1
->vdwtype
, ir2
->vdwtype
);
795 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
);
796 cmp_real(fp
, "inputrec->rvdw", -1, ir1
->rvdw
, ir2
->rvdw
, ftol
, abstol
);
797 cmp_real(fp
, "inputrec->epsilon_r", -1, ir1
->epsilon_r
, ir2
->epsilon_r
, ftol
, abstol
);
798 cmp_real(fp
, "inputrec->epsilon_rf", -1, ir1
->epsilon_rf
, ir2
->epsilon_rf
, ftol
, abstol
);
799 cmp_real(fp
, "inputrec->tabext", -1, ir1
->tabext
, ir2
->tabext
, ftol
, abstol
);
800 cmp_int(fp
, "inputrec->implicit_solvent", -1, ir1
->implicit_solvent
, ir2
->implicit_solvent
);
801 cmp_int(fp
, "inputrec->gb_algorithm", -1, ir1
->gb_algorithm
, ir2
->gb_algorithm
);
802 cmp_int(fp
, "inputrec->nstgbradii", -1, ir1
->nstgbradii
, ir2
->nstgbradii
);
803 cmp_real(fp
, "inputrec->rgbradii", -1, ir1
->rgbradii
, ir2
->rgbradii
, ftol
, abstol
);
804 cmp_real(fp
, "inputrec->gb_saltconc", -1, ir1
->gb_saltconc
, ir2
->gb_saltconc
, ftol
, abstol
);
805 cmp_real(fp
, "inputrec->gb_epsilon_solvent", -1, ir1
->gb_epsilon_solvent
, ir2
->gb_epsilon_solvent
, ftol
, abstol
);
806 cmp_real(fp
, "inputrec->gb_obc_alpha", -1, ir1
->gb_obc_alpha
, ir2
->gb_obc_alpha
, ftol
, abstol
);
807 cmp_real(fp
, "inputrec->gb_obc_beta", -1, ir1
->gb_obc_beta
, ir2
->gb_obc_beta
, ftol
, abstol
);
808 cmp_real(fp
, "inputrec->gb_obc_gamma", -1, ir1
->gb_obc_gamma
, ir2
->gb_obc_gamma
, ftol
, abstol
);
809 cmp_real(fp
, "inputrec->gb_dielectric_offset", -1, ir1
->gb_dielectric_offset
, ir2
->gb_dielectric_offset
, ftol
, abstol
);
810 cmp_int(fp
, "inputrec->sa_algorithm", -1, ir1
->sa_algorithm
, ir2
->sa_algorithm
);
811 cmp_real(fp
, "inputrec->sa_surface_tension", -1, ir1
->sa_surface_tension
, ir2
->sa_surface_tension
, ftol
, abstol
);
813 cmp_int(fp
, "inputrec->eDispCorr", -1, ir1
->eDispCorr
, ir2
->eDispCorr
);
814 cmp_real(fp
, "inputrec->shake_tol", -1, ir1
->shake_tol
, ir2
->shake_tol
, ftol
, abstol
);
815 cmp_int(fp
, "inputrec->efep", -1, ir1
->efep
, ir2
->efep
);
816 cmp_fepvals(fp
, ir1
->fepvals
, ir2
->fepvals
, ftol
, abstol
);
817 cmp_int(fp
, "inputrec->bSimTemp", -1, ir1
->bSimTemp
, ir2
->bSimTemp
);
818 if ((ir1
->bSimTemp
== ir2
->bSimTemp
) && (ir1
->bSimTemp
))
820 cmp_simtempvals(fp
, ir1
->simtempvals
, ir2
->simtempvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
822 cmp_int(fp
, "inputrec->bExpanded", -1, ir1
->bExpanded
, ir2
->bExpanded
);
823 if ((ir1
->bExpanded
== ir2
->bExpanded
) && (ir1
->bExpanded
))
825 cmp_expandedvals(fp
, ir1
->expandedvals
, ir2
->expandedvals
, std::min(ir1
->fepvals
->n_lambda
, ir2
->fepvals
->n_lambda
), ftol
, abstol
);
827 cmp_int(fp
, "inputrec->nwall", -1, ir1
->nwall
, ir2
->nwall
);
828 cmp_int(fp
, "inputrec->wall_type", -1, ir1
->wall_type
, ir2
->wall_type
);
829 cmp_int(fp
, "inputrec->wall_atomtype[0]", -1, ir1
->wall_atomtype
[0], ir2
->wall_atomtype
[0]);
830 cmp_int(fp
, "inputrec->wall_atomtype[1]", -1, ir1
->wall_atomtype
[1], ir2
->wall_atomtype
[1]);
831 cmp_real(fp
, "inputrec->wall_density[0]", -1, ir1
->wall_density
[0], ir2
->wall_density
[0], ftol
, abstol
);
832 cmp_real(fp
, "inputrec->wall_density[1]", -1, ir1
->wall_density
[1], ir2
->wall_density
[1], ftol
, abstol
);
833 cmp_real(fp
, "inputrec->wall_ewald_zfac", -1, ir1
->wall_ewald_zfac
, ir2
->wall_ewald_zfac
, ftol
, abstol
);
835 cmp_bool(fp
, "inputrec->bPull", -1, ir1
->bPull
, ir2
->bPull
);
836 if (ir1
->bPull
&& ir2
->bPull
)
841 cmp_int(fp
, "inputrec->eDisre", -1, ir1
->eDisre
, ir2
->eDisre
);
842 cmp_real(fp
, "inputrec->dr_fc", -1, ir1
->dr_fc
, ir2
->dr_fc
, ftol
, abstol
);
843 cmp_int(fp
, "inputrec->eDisreWeighting", -1, ir1
->eDisreWeighting
, ir2
->eDisreWeighting
);
844 cmp_int(fp
, "inputrec->bDisreMixed", -1, ir1
->bDisreMixed
, ir2
->bDisreMixed
);
845 cmp_int(fp
, "inputrec->nstdisreout", -1, ir1
->nstdisreout
, ir2
->nstdisreout
);
846 cmp_real(fp
, "inputrec->dr_tau", -1, ir1
->dr_tau
, ir2
->dr_tau
, ftol
, abstol
);
847 cmp_real(fp
, "inputrec->orires_fc", -1, ir1
->orires_fc
, ir2
->orires_fc
, ftol
, abstol
);
848 cmp_real(fp
, "inputrec->orires_tau", -1, ir1
->orires_tau
, ir2
->orires_tau
, ftol
, abstol
);
849 cmp_int(fp
, "inputrec->nstorireout", -1, ir1
->nstorireout
, ir2
->nstorireout
);
850 cmp_real(fp
, "inputrec->em_stepsize", -1, ir1
->em_stepsize
, ir2
->em_stepsize
, ftol
, abstol
);
851 cmp_real(fp
, "inputrec->em_tol", -1, ir1
->em_tol
, ir2
->em_tol
, ftol
, abstol
);
852 cmp_int(fp
, "inputrec->niter", -1, ir1
->niter
, ir2
->niter
);
853 cmp_real(fp
, "inputrec->fc_stepsize", -1, ir1
->fc_stepsize
, ir2
->fc_stepsize
, ftol
, abstol
);
854 cmp_int(fp
, "inputrec->nstcgsteep", -1, ir1
->nstcgsteep
, ir2
->nstcgsteep
);
855 cmp_int(fp
, "inputrec->nbfgscorr", 0, ir1
->nbfgscorr
, ir2
->nbfgscorr
);
856 cmp_int(fp
, "inputrec->eConstrAlg", -1, ir1
->eConstrAlg
, ir2
->eConstrAlg
);
857 cmp_int(fp
, "inputrec->nProjOrder", -1, ir1
->nProjOrder
, ir2
->nProjOrder
);
858 cmp_real(fp
, "inputrec->LincsWarnAngle", -1, ir1
->LincsWarnAngle
, ir2
->LincsWarnAngle
, ftol
, abstol
);
859 cmp_int(fp
, "inputrec->nLincsIter", -1, ir1
->nLincsIter
, ir2
->nLincsIter
);
860 cmp_real(fp
, "inputrec->bd_fric", -1, ir1
->bd_fric
, ir2
->bd_fric
, ftol
, abstol
);
861 cmp_int64(fp
, "inputrec->ld_seed", ir1
->ld_seed
, ir2
->ld_seed
);
862 cmp_real(fp
, "inputrec->cos_accel", -1, ir1
->cos_accel
, ir2
->cos_accel
, ftol
, abstol
);
863 cmp_rvec(fp
, "inputrec->deform(a)", -1, ir1
->deform
[XX
], ir2
->deform
[XX
], ftol
, abstol
);
864 cmp_rvec(fp
, "inputrec->deform(b)", -1, ir1
->deform
[YY
], ir2
->deform
[YY
], ftol
, abstol
);
865 cmp_rvec(fp
, "inputrec->deform(c)", -1, ir1
->deform
[ZZ
], ir2
->deform
[ZZ
], ftol
, abstol
);
868 cmp_int(fp
, "inputrec->userint1", -1, ir1
->userint1
, ir2
->userint1
);
869 cmp_int(fp
, "inputrec->userint2", -1, ir1
->userint2
, ir2
->userint2
);
870 cmp_int(fp
, "inputrec->userint3", -1, ir1
->userint3
, ir2
->userint3
);
871 cmp_int(fp
, "inputrec->userint4", -1, ir1
->userint4
, ir2
->userint4
);
872 cmp_real(fp
, "inputrec->userreal1", -1, ir1
->userreal1
, ir2
->userreal1
, ftol
, abstol
);
873 cmp_real(fp
, "inputrec->userreal2", -1, ir1
->userreal2
, ir2
->userreal2
, ftol
, abstol
);
874 cmp_real(fp
, "inputrec->userreal3", -1, ir1
->userreal3
, ir2
->userreal3
, ftol
, abstol
);
875 cmp_real(fp
, "inputrec->userreal4", -1, ir1
->userreal4
, ir2
->userreal4
, ftol
, abstol
);
876 cmp_grpopts(fp
, &(ir1
->opts
), &(ir2
->opts
), ftol
, abstol
);
877 cmp_cosines(fp
, "ex", ir1
->ex
, ir2
->ex
, ftol
, abstol
);
878 cmp_cosines(fp
, "et", ir1
->et
, ir2
->et
, ftol
, abstol
);
881 static void comp_pull_AB(FILE *fp
, pull_params_t
*pull
, real ftol
, real abstol
)
885 for (i
= 0; i
< pull
->ncoord
; i
++)
887 fprintf(fp
, "comparing pull coord %d\n", i
);
888 cmp_real(fp
, "pull-coord->k", -1, pull
->coord
[i
].k
, pull
->coord
[i
].kB
, ftol
, abstol
);
892 static void comp_state(t_state
*st1
, t_state
*st2
,
893 gmx_bool bRMSD
, real ftol
, real abstol
)
897 fprintf(stdout
, "comparing flags\n");
898 cmp_int(stdout
, "flags", -1, st1
->flags
, st2
->flags
);
899 fprintf(stdout
, "comparing box\n");
900 cmp_rvecs(stdout
, "box", DIM
, st1
->box
, st2
->box
, FALSE
, ftol
, abstol
);
901 fprintf(stdout
, "comparing box_rel\n");
902 cmp_rvecs(stdout
, "box_rel", DIM
, st1
->box_rel
, st2
->box_rel
, FALSE
, ftol
, abstol
);
903 fprintf(stdout
, "comparing boxv\n");
904 cmp_rvecs(stdout
, "boxv", DIM
, st1
->boxv
, st2
->boxv
, FALSE
, ftol
, abstol
);
905 if (st1
->flags
& (1<<estSVIR_PREV
))
907 fprintf(stdout
, "comparing shake vir_prev\n");
908 cmp_rvecs(stdout
, "svir_prev", DIM
, st1
->svir_prev
, st2
->svir_prev
, FALSE
, ftol
, abstol
);
910 if (st1
->flags
& (1<<estFVIR_PREV
))
912 fprintf(stdout
, "comparing force vir_prev\n");
913 cmp_rvecs(stdout
, "fvir_prev", DIM
, st1
->fvir_prev
, st2
->fvir_prev
, FALSE
, ftol
, abstol
);
915 if (st1
->flags
& (1<<estPRES_PREV
))
917 fprintf(stdout
, "comparing prev_pres\n");
918 cmp_rvecs(stdout
, "pres_prev", DIM
, st1
->pres_prev
, st2
->pres_prev
, FALSE
, ftol
, abstol
);
920 cmp_int(stdout
, "ngtc", -1, st1
->ngtc
, st2
->ngtc
);
921 cmp_int(stdout
, "nhchainlength", -1, st1
->nhchainlength
, st2
->nhchainlength
);
922 if (st1
->ngtc
== st2
->ngtc
&& st1
->nhchainlength
== st2
->nhchainlength
)
924 for (i
= 0; i
< st1
->ngtc
; i
++)
926 nc
= i
*st1
->nhchainlength
;
927 for (j
= 0; j
< nc
; j
++)
929 cmp_real(stdout
, "nosehoover_xi",
930 i
, st1
->nosehoover_xi
[nc
+j
], st2
->nosehoover_xi
[nc
+j
], ftol
, abstol
);
934 cmp_int(stdout
, "nnhpres", -1, st1
->nnhpres
, st2
->nnhpres
);
935 if (st1
->nnhpres
== st2
->nnhpres
&& st1
->nhchainlength
== st2
->nhchainlength
)
937 for (i
= 0; i
< st1
->nnhpres
; i
++)
939 nc
= i
*st1
->nhchainlength
;
940 for (j
= 0; j
< nc
; j
++)
942 cmp_real(stdout
, "nosehoover_xi",
943 i
, st1
->nhpres_xi
[nc
+j
], st2
->nhpres_xi
[nc
+j
], ftol
, abstol
);
948 cmp_int(stdout
, "natoms", -1, st1
->natoms
, st2
->natoms
);
949 if (st1
->natoms
== st2
->natoms
)
951 if ((st1
->flags
& (1<<estX
)) && (st2
->flags
& (1<<estX
)))
953 fprintf(stdout
, "comparing x\n");
954 cmp_rvecs(stdout
, "x", st1
->natoms
, st1
->x
, st2
->x
, bRMSD
, ftol
, abstol
);
956 if ((st1
->flags
& (1<<estV
)) && (st2
->flags
& (1<<estV
)))
958 fprintf(stdout
, "comparing v\n");
959 cmp_rvecs(stdout
, "v", st1
->natoms
, st1
->v
, st2
->v
, bRMSD
, ftol
, abstol
);
964 void comp_tpx(const char *fn1
, const char *fn2
,
965 gmx_bool bRMSD
, real ftol
, real abstol
)
976 for (i
= 0; i
< (fn2
? 2 : 1); i
++)
978 read_tpx_state(ff
[i
], &(ir
[i
]), &state
[i
], &(mtop
[i
]));
982 cmp_inputrec(stdout
, &ir
[0], &ir
[1], ftol
, abstol
);
983 /* Convert gmx_mtop_t to t_topology.
984 * We should implement direct mtop comparison,
985 * but it might be useful to keep t_topology comparison as an option.
987 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
988 top
[1] = gmx_mtop_t_to_t_topology(&mtop
[1]);
989 cmp_top(stdout
, &top
[0], &top
[1], ftol
, abstol
);
990 cmp_groups(stdout
, &mtop
[0].groups
, &mtop
[1].groups
,
991 mtop
[0].natoms
, mtop
[1].natoms
);
992 comp_state(&state
[0], &state
[1], bRMSD
, ftol
, abstol
);
996 if (ir
[0].efep
== efepNO
)
998 fprintf(stdout
, "inputrec->efep = %s\n", efep_names
[ir
[0].efep
]);
1004 comp_pull_AB(stdout
, ir
->pull
, ftol
, abstol
);
1006 /* Convert gmx_mtop_t to t_topology.
1007 * We should implement direct mtop comparison,
1008 * but it might be useful to keep t_topology comparison as an option.
1010 top
[0] = gmx_mtop_t_to_t_topology(&mtop
[0]);
1011 cmp_top(stdout
, &top
[0], NULL
, ftol
, abstol
);
1016 void comp_frame(FILE *fp
, t_trxframe
*fr1
, t_trxframe
*fr2
,
1017 gmx_bool bRMSD
, real ftol
, real abstol
)
1020 cmp_int(fp
, "not_ok", -1, fr1
->not_ok
, fr2
->not_ok
);
1021 cmp_int(fp
, "natoms", -1, fr1
->natoms
, fr2
->natoms
);
1022 if (cmp_bool(fp
, "bTitle", -1, fr1
->bTitle
, fr2
->bTitle
))
1024 cmp_str(fp
, "title", -1, fr1
->title
, fr2
->title
);
1026 if (cmp_bool(fp
, "bStep", -1, fr1
->bStep
, fr2
->bStep
))
1028 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1030 cmp_int(fp
, "step", -1, fr1
->step
, fr2
->step
);
1031 if (cmp_bool(fp
, "bTime", -1, fr1
->bTime
, fr2
->bTime
))
1033 cmp_real(fp
, "time", -1, fr1
->time
, fr2
->time
, ftol
, abstol
);
1035 if (cmp_bool(fp
, "bLambda", -1, fr1
->bLambda
, fr2
->bLambda
))
1037 cmp_real(fp
, "lambda", -1, fr1
->lambda
, fr2
->lambda
, ftol
, abstol
);
1039 if (cmp_bool(fp
, "bAtoms", -1, fr1
->bAtoms
, fr2
->bAtoms
))
1041 cmp_atoms(fp
, fr1
->atoms
, fr2
->atoms
, ftol
, abstol
);
1043 if (cmp_bool(fp
, "bPrec", -1, fr1
->bPrec
, fr2
->bPrec
))
1045 cmp_real(fp
, "prec", -1, fr1
->prec
, fr2
->prec
, ftol
, abstol
);
1047 if (cmp_bool(fp
, "bX", -1, fr1
->bX
, fr2
->bX
))
1049 cmp_rvecs(fp
, "x", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->x
, fr2
->x
, bRMSD
, ftol
, abstol
);
1051 if (cmp_bool(fp
, "bV", -1, fr1
->bV
, fr2
->bV
))
1053 cmp_rvecs(fp
, "v", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->v
, fr2
->v
, bRMSD
, ftol
, abstol
);
1055 if (cmp_bool(fp
, "bF", -1, fr1
->bF
, fr2
->bF
))
1057 cmp_rvecs(fp
, "f", std::min(fr1
->natoms
, fr2
->natoms
), fr1
->f
, fr2
->f
, bRMSD
, ftol
, abstol
);
1059 if (cmp_bool(fp
, "bBox", -1, fr1
->bBox
, fr2
->bBox
))
1061 cmp_rvecs(fp
, "box", 3, fr1
->box
, fr2
->box
, FALSE
, ftol
, abstol
);
1065 void comp_trx(const gmx_output_env_t
*oenv
, const char *fn1
, const char *fn2
,
1066 gmx_bool bRMSD
, real ftol
, real abstol
)
1071 t_trxstatus
*status
[2];
1076 fprintf(stderr
, "Comparing trajectory files %s and %s\n", fn1
, fn2
);
1077 for (i
= 0; i
< 2; i
++)
1079 b
[i
] = read_first_frame(oenv
, &status
[i
], fn
[i
], &fr
[i
], TRX_READ_X
|TRX_READ_V
|TRX_READ_F
);
1086 comp_frame(stdout
, &(fr
[0]), &(fr
[1]), bRMSD
, ftol
, abstol
);
1088 for (i
= 0; i
< 2; i
++)
1090 b
[i
] = read_next_frame(oenv
, status
[i
], &fr
[i
]);
1093 while (b
[0] && b
[1]);
1095 for (i
= 0; i
< 2; i
++)
1097 if (b
[i
] && !b
[1-i
])
1099 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn
[1-i
], fn
[i
]);
1101 close_trj(status
[i
]);
1106 fprintf(stdout
, "\nBoth files read correctly\n");
1110 static real
ener_tensor_diag(int n
, int *ind1
, int *ind2
,
1113 t_energy e1
[], t_energy e2
[])
1122 d2
= tensi
[i
] - d1
*DIM
;
1124 /* Find the diagonal elements d1 and d2 */
1125 len
= std::strlen(enm1
[ind1
[i
]].name
);
1129 for (j
= 0; j
< n
; j
++)
1131 if (tensi
[j
] >= 0 &&
1132 std::strlen(enm1
[ind1
[j
]].name
) == len
&&
1133 std::strncmp(enm1
[ind1
[i
]].name
, enm1
[ind1
[j
]].name
, len
-2) == 0 &&
1134 (tensi
[j
] == d1
*DIM
+d1
|| tensi
[j
] == d2
*DIM
+d2
))
1136 prod1
*= fabs(e1
[ind1
[j
]].e
);
1137 prod2
*= fabs(e2
[ind2
[j
]].e
);
1144 return 0.5*(std::sqrt(prod1
) + std::sqrt(prod2
));
1152 static gmx_bool
enernm_equal(const char *nm1
, const char *nm2
)
1156 len1
= std::strlen(nm1
);
1157 len2
= std::strlen(nm2
);
1159 /* Remove " (bar)" at the end of a name */
1160 if (len1
> 6 && std::strcmp(nm1
+len1
-6, " (bar)") == 0)
1164 if (len2
> 6 && std::strcmp(nm2
+len2
-6, " (bar)") == 0)
1169 return (len1
== len2
&& gmx_strncasecmp(nm1
, nm2
, len1
) == 0);
1172 static void cmp_energies(FILE *fp
, int step1
, int step2
,
1173 t_energy e1
[], t_energy e2
[],
1175 real ftol
, real abstol
,
1176 int nre
, int *ind1
, int *ind2
, int maxener
)
1179 int *tensi
, len
, d1
, d2
;
1180 real ftol_i
, abstol_i
;
1182 snew(tensi
, maxener
);
1183 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1184 for (i
= 0; (i
< maxener
); i
++)
1188 len
= std::strlen(enm1
[ii
].name
);
1189 if (len
> 3 && enm1
[ii
].name
[len
-3] == '-')
1191 d1
= enm1
[ii
].name
[len
-2] - 'X';
1192 d2
= enm1
[ii
].name
[len
-1] - 'X';
1193 if (d1
>= 0 && d1
< DIM
&&
1194 d2
>= 0 && d2
< DIM
)
1196 tensi
[i
] = d1
*DIM
+ d2
;
1201 for (i
= 0; (i
< maxener
); i
++)
1203 /* Check if this is an off-diagonal tensor element */
1204 if (tensi
[i
] >= 0 && tensi
[i
] != 0 && tensi
[i
] != 4 && tensi
[i
] != 8)
1206 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1208 /* Do the relative tolerance through an absolute tolerance times
1209 * the size of diagonal components of the tensor.
1211 abstol_i
= ftol
*ener_tensor_diag(nre
, ind1
, ind2
, enm1
, tensi
, i
, e1
, e2
);
1214 fprintf(debug
, "tensor '%s' val %f diag %f\n",
1215 enm1
[i
].name
, e1
[i
].e
, abstol_i
/ftol
);
1219 /* We found a diagonal, we need to check with the minimum tolerance */
1220 abstol_i
= std::min(abstol_i
, abstol
);
1224 /* We did not find a diagonal, ignore the relative tolerance check */
1233 if (!equal_real(e1
[ind1
[i
]].e
, e2
[ind2
[i
]].e
, ftol_i
, abstol_i
))
1235 fprintf(fp
, "%-15s step %3d: %12g, step %3d: %12g\n",
1237 step1
, e1
[ind1
[i
]].e
,
1238 step2
, e2
[ind2
[i
]].e
);
1246 static void cmp_disres(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1249 char bav
[64], bt
[64], bs
[22];
1251 cmp_int(stdout
, "ndisre", -1, fr1
->ndisre
, fr2
->ndisre
);
1252 if ((fr1
->ndisre
== fr2
->ndisre
) && (fr1
->ndisre
> 0))
1254 sprintf(bav
, "step %s: disre rav", gmx_step_str(fr1
->step
, bs
));
1255 sprintf(bt
, "step %s: disre rt", gmx_step_str(fr1
->step
, bs
));
1256 for (i
= 0; (i
< fr1
->ndisre
); i
++)
1258 cmp_real(stdout
, bav
, i
, fr1
->disre_rm3tav
[i
], fr2
->disre_rm3tav
[i
], ftol
, abstol
);
1259 cmp_real(stdout
, bt
, i
, fr1
->disre_rt
[i
], fr2
->disre_rt
[i
], ftol
, abstol
);
1265 static void cmp_eblocks(t_enxframe
*fr1
, t_enxframe
*fr2
, real ftol
, real abstol
)
1268 char buf
[64], bs
[22];
1270 cmp_int(stdout
, "nblock", -1, fr1
->nblock
, fr2
->nblock
);
1271 if ((fr1
->nblock
== fr2
->nblock
) && (fr1
->nblock
> 0))
1273 for (j
= 0; (j
< fr1
->nblock
); j
++)
1275 t_enxblock
*b1
, *b2
; /* convenience vars */
1277 b1
= &(fr1
->block
[j
]);
1278 b2
= &(fr2
->block
[j
]);
1280 sprintf(buf
, "step %s: block[%d]", gmx_step_str(fr1
->step
, bs
), j
);
1281 cmp_int(stdout
, buf
, -1, b1
->nsub
, b2
->nsub
);
1282 cmp_int(stdout
, buf
, -1, b1
->id
, b2
->id
);
1284 if ( (b1
->nsub
== b2
->nsub
) && (b1
->id
== b2
->id
) )
1286 for (i
= 0; i
< b1
->nsub
; i
++)
1288 t_enxsubblock
*s1
, *s2
;
1293 cmp_int(stdout
, buf
, -1, (int)s1
->type
, (int)s2
->type
);
1294 cmp_int64(stdout
, buf
, s1
->nr
, s2
->nr
);
1296 if ((s1
->type
== s2
->type
) && (s1
->nr
== s2
->nr
))
1300 case xdr_datatype_float
:
1301 for (k
= 0; k
< s1
->nr
; k
++)
1303 cmp_float(stdout
, buf
, i
,
1304 s1
->fval
[k
], s2
->fval
[k
],
1308 case xdr_datatype_double
:
1309 for (k
= 0; k
< s1
->nr
; k
++)
1311 cmp_double(stdout
, buf
, i
,
1312 s1
->dval
[k
], s2
->dval
[k
],
1316 case xdr_datatype_int
:
1317 for (k
= 0; k
< s1
->nr
; k
++)
1319 cmp_int(stdout
, buf
, i
,
1320 s1
->ival
[k
], s2
->ival
[k
]);
1323 case xdr_datatype_int64
:
1324 for (k
= 0; k
< s1
->nr
; k
++)
1326 cmp_int64(stdout
, buf
,
1327 s1
->lval
[k
], s2
->lval
[k
]);
1330 case xdr_datatype_char
:
1331 for (k
= 0; k
< s1
->nr
; k
++)
1333 cmp_uc(stdout
, buf
, i
,
1334 s1
->cval
[k
], s2
->cval
[k
]);
1337 case xdr_datatype_string
:
1338 for (k
= 0; k
< s1
->nr
; k
++)
1340 cmp_str(stdout
, buf
, i
,
1341 s1
->sval
[k
], s2
->sval
[k
]);
1345 gmx_incons("Unknown data type!!");
1354 void comp_enx(const char *fn1
, const char *fn2
, real ftol
, real abstol
, const char *lastener
)
1356 int nre
, nre1
, nre2
;
1357 ener_file_t in1
, in2
;
1358 int i
, j
, maxener
, *ind1
, *ind2
, *have
;
1359 gmx_enxnm_t
*enm1
= NULL
, *enm2
= NULL
;
1360 t_enxframe
*fr1
, *fr2
;
1363 fprintf(stdout
, "comparing energy file %s and %s\n\n", fn1
, fn2
);
1365 in1
= open_enx(fn1
, "r");
1366 in2
= open_enx(fn2
, "r");
1367 do_enxnms(in1
, &nre1
, &enm1
);
1368 do_enxnms(in2
, &nre2
, &enm2
);
1371 fprintf(stdout
, "There are %d and %d terms in the energy files\n\n",
1376 fprintf(stdout
, "There are %d terms in the energy files\n\n", nre1
);
1383 for (i
= 0; i
< nre1
; i
++)
1385 for (j
= 0; j
< nre2
; j
++)
1387 if (enernm_equal(enm1
[i
].name
, enm2
[j
].name
))
1396 if (nre
== 0 || ind1
[nre
-1] != i
)
1398 cmp_str(stdout
, "enm", i
, enm1
[i
].name
, "-");
1401 for (i
= 0; i
< nre2
; i
++)
1405 cmp_str(stdout
, "enm", i
, "-", enm2
[i
].name
);
1410 for (i
= 0; i
< nre
; i
++)
1412 if ((lastener
!= NULL
) && (std::strstr(enm1
[i
].name
, lastener
) != NULL
))
1419 fprintf(stdout
, "There are %d terms to compare in the energy files\n\n",
1422 for (i
= 0; i
< maxener
; i
++)
1424 cmp_str(stdout
, "unit", i
, enm1
[ind1
[i
]].unit
, enm2
[ind2
[i
]].unit
);
1431 b1
= do_enx(in1
, fr1
);
1432 b2
= do_enx(in2
, fr2
);
1435 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn2
, fn1
);
1439 fprintf(stdout
, "\nEnd of file on %s but not on %s\n", fn1
, fn2
);
1441 else if (!b1
&& !b2
)
1443 fprintf(stdout
, "\nFiles read successfully\n");
1447 cmp_real(stdout
, "t", -1, fr1
->t
, fr2
->t
, ftol
, abstol
);
1448 cmp_int(stdout
, "step", -1, fr1
->step
, fr2
->step
);
1449 /* We don't want to print the nre mismatch for every frame */
1450 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1451 if ((fr1
->nre
>= nre
) && (fr2
->nre
>= nre
))
1453 cmp_energies(stdout
, fr1
->step
, fr1
->step
, fr1
->ener
, fr2
->ener
,
1454 enm1
, ftol
, abstol
, nre
, ind1
, ind2
, maxener
);
1456 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1457 cmp_eblocks(fr1
, fr2
, ftol
, abstol
);