src/gromacs/tools/compare.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38
  39 #include "gmxpre.h"
  40
  41 #include <cmath>
  42 #include <cstdio>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  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"
  62
  63 static void cmp_int(FILE *fp, const char *s, int index, int i1, int i2)
  64 {
  65     if (i1 != i2)
  66     {
  67         if (index != -1)
  68         {
  69             fprintf(fp, "%s[%d] (%d - %d)\n", s, index, i1, i2);
  70         }
  71         else
  72         {
  73             fprintf(fp, "%s (%d - %d)\n", s, i1, i2);
  74         }
  75     }
  76 }
  77
  78 static void cmp_int64(FILE *fp, const char *s, gmx_int64_t i1, gmx_int64_t i2)
  79 {
  80     if (i1 != i2)
  81     {
  82         fprintf(fp, "%s (", s);
  83         fprintf(fp, "%" GMX_PRId64, i1);
  84         fprintf(fp, " - ");
  85         fprintf(fp, "%" GMX_PRId64, i2);
  86         fprintf(fp, ")\n");
  87     }
  88 }
  89
  90 static void cmp_us(FILE *fp, const char *s, int index, unsigned short i1, unsigned short i2)
  91 {
  92     if (i1 != i2)
  93     {
  94         if (index != -1)
  95         {
  96             fprintf(fp, "%s[%d] (%hu - %hu)\n", s, index, i1, i2);
  97         }
  98         else
  99         {
 100             fprintf(fp, "%s (%hu - %hu)\n", s, i1, i2);
 101         }
 102     }
 103 }
 104
 105 static void cmp_uc(FILE *fp, const char *s, int index, unsigned char i1, unsigned char i2)
 106 {
 107     if (i1 != i2)
 108     {
 109         if (index != -1)
 110         {
 111             fprintf(fp, "%s[%d] (%d - %d)\n", s, index, i1, i2);
 112         }
 113         else
 114         {
 115             fprintf(fp, "%s (%d - %d)\n", s, i1, i2);
 116         }
 117     }
 118 }
 119
 120 static gmx_bool cmp_bool(FILE *fp, const char *s, int index, gmx_bool b1, gmx_bool b2)
 121 {
 122     if (b1)
 123     {
 124         b1 = 1;
 125     }
 126     else
 127     {
 128         b1 = 0;
 129     }
 130     if (b2)
 131     {
 132         b2 = 1;
 133     }
 134     else
 135     {
 136         b2 = 0;
 137     }
 138     if (b1 != b2)
 139     {
 140         if (index != -1)
 141         {
 142             fprintf(fp, "%s[%d] (%s - %s)\n", s, index,
 143                     gmx::boolToString(b1), gmx::boolToString(b2));
 144         }
 145         else
 146         {
 147             fprintf(fp, "%s (%s - %s)\n", s,
 148                     gmx::boolToString(b1), gmx::boolToString(b2));
 149         }
 150     }
 151     return b1 && b2;
 152 }
 153
 154 static void cmp_str(FILE *fp, const char *s, int index,
 155                     const char *s1, const char *s2)
 156 {
 157     if (std::strcmp(s1, s2) != 0)
 158     {
 159         if (index != -1)
 160         {
 161             fprintf(fp, "%s[%d] (%s - %s)\n", s, index, s1, s2);
 162         }
 163         else
 164         {
 165             fprintf(fp, "%s (%s - %s)\n", s, s1, s2);
 166         }
 167     }
 168 }
 169
 170 static gmx_bool equal_real(real i1, real i2, real ftol, real abstol)
 171 {
 172     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 173 }
 174
 175 static gmx_bool equal_float(float i1, float i2, float ftol, float abstol)
 176 {
 177     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 178 }
 179
 180 static gmx_bool equal_double(double i1, double i2, real ftol, real abstol)
 181 {
 182     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 183 }
 184
 185 static void
 186 cmp_real(FILE *fp, const char *s, int index, real i1, real i2, real ftol, real abstol)
 187 {
 188     if (!equal_real(i1, i2, ftol, abstol))
 189     {
 190         if (index != -1)
 191         {
 192             fprintf(fp, "%s[%2d] (%e - %e)\n", s, index, i1, i2);
 193         }
 194         else
 195         {
 196             fprintf(fp, "%s (%e - %e)\n", s, i1, i2);
 197         }
 198     }
 199 }
 200
 201 static void
 202 cmp_float(FILE *fp, const char *s, int index, float i1, float i2, float ftol, float abstol)
 203 {
 204     if (!equal_float(i1, i2, ftol, abstol))
 205     {
 206         if (index != -1)
 207         {
 208             fprintf(fp, "%s[%2d] (%e - %e)\n", s, index, i1, i2);
 209         }
 210         else
 211         {
 212             fprintf(fp, "%s (%e - %e)\n", s, i1, i2);
 213         }
 214     }
 215 }
 216
 217
 218
 219 static void
 220 cmp_double(FILE *fp, const char *s, int index, double i1, double i2, double ftol, double abstol)
 221 {
 222     if (!equal_double(i1, i2, ftol, abstol))
 223     {
 224         if (index != -1)
 225         {
 226             fprintf(fp, "%s[%2d] (%16.9e - %16.9e)\n", s, index, i1, i2);
 227         }
 228         else
 229         {
 230             fprintf(fp, "%s (%16.9e - %16.9e)\n", s, i1, i2);
 231         }
 232     }
 233 }
 234
 235 static void cmp_rvec(FILE *fp, const char *s, int index, rvec i1, rvec i2, real ftol, real abstol)
 236 {
 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))
 240     {
 241         if (index != -1)
 242         {
 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]);
 245         }
 246         else
 247         {
 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]);
 250         }
 251     }
 252 }
 253
 254 static void cmp_ivec(FILE *fp, const char *s, int index, ivec i1, ivec i2)
 255 {
 256     if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ]))
 257     {
 258         if (index != -1)
 259         {
 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]);
 262         }
 263         else
 264         {
 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]);
 267         }
 268     }
 269 }
 270
 271 static void cmp_ilist(FILE *fp, int ftype, t_ilist *il1, t_ilist *il2)
 272 {
 273     int  i;
 274     char buf[256];
 275
 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)))
 283     {
 284         fprintf(fp, "Comparing radically different topologies - %s is different\n",
 285                 buf);
 286     }
 287     else
 288     {
 289         for (i = 0; (i < il1->nr); i++)
 290         {
 291             cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
 292         }
 293     }
 294 }
 295
 296 void cmp_iparm(FILE *fp, const char *s, t_functype ft,
 297                t_iparams ip1, t_iparams ip2, real ftol, real abstol)
 298 {
 299     int      i;
 300     gmx_bool bDiff;
 301
 302     bDiff = FALSE;
 303     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 304     {
 305         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
 306     }
 307     if (bDiff)
 308     {
 309         fprintf(fp, "%s1: ", s);
 310         pr_iparams(fp, ft, &ip1);
 311         fprintf(fp, "%s2: ", s);
 312         pr_iparams(fp, ft, &ip2);
 313     }
 314 }
 315
 316 void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft, t_iparams ip1, real ftol, real abstol)
 317 {
 318     int      nrfpA, nrfpB, p0, i;
 319     gmx_bool bDiff;
 320
 321     /* Normally the first parameter is perturbable */
 322     p0    = 0;
 323     nrfpA = interaction_function[ft].nrfpA;
 324     nrfpB = interaction_function[ft].nrfpB;
 325     if (ft == F_PDIHS)
 326     {
 327         nrfpB = 2;
 328     }
 329     else if (interaction_function[ft].flags & IF_TABULATED)
 330     {
 331         /* For tabulated interactions only the second parameter is perturbable */
 332         p0    = 1;
 333         nrfpB = 1;
 334     }
 335     bDiff = FALSE;
 336     for (i = 0; i < nrfpB && !bDiff; i++)
 337     {
 338         bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
 339     }
 340     if (bDiff)
 341     {
 342         fprintf(fp, "%s: ", s);
 343         pr_iparams(fp, ft, &ip1);
 344     }
 345 }
 346
 347 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
 348 {
 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)
 353     {
 354         int g;
 355
 356         for (g = 0; g < cmap1->ngrid; g++)
 357         {
 358             int i;
 359
 360             fprintf(fp, "comparing cmap %d\n", g);
 361
 362             for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
 363             {
 364                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
 365             }
 366         }
 367     }
 368 }
 369
 370 static void cmp_idef(FILE *fp, t_idef *id1, t_idef *id2, real ftol, real abstol)
 371 {
 372     int  i;
 373     char buf1[64], buf2[64];
 374
 375     fprintf(fp, "comparing idef\n");
 376     if (id2)
 377     {
 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++)
 381         {
 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);
 387         }
 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++)
 391         {
 392             cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
 393         }
 394     }
 395     else
 396     {
 397         for (i = 0; (i < id1->ntypes); i++)
 398         {
 399             cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
 400         }
 401     }
 402 }
 403
 404 static void cmp_block(FILE *fp, t_block *b1, t_block *b2, const char *s)
 405 {
 406     char buf[32];
 407
 408     fprintf(fp, "comparing block %s\n", s);
 409     sprintf(buf, "%s.nr", s);
 410     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 411 }
 412
 413 static void cmp_blocka(FILE *fp, t_blocka *b1, t_blocka *b2, const char *s)
 414 {
 415     char buf[32];
 416
 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);
 422 }
 423
 424 static void cmp_atom(FILE *fp, int index, t_atom *a1, t_atom *a2, real ftol, real abstol)
 425 {
 426     if (a2)
 427     {
 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);
 437     }
 438     else
 439     {
 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);
 443     }
 444 }
 445
 446 static void cmp_atoms(FILE *fp, t_atoms *a1, t_atoms *a2, real ftol, real abstol)
 447 {
 448     int i;
 449
 450     fprintf(fp, "comparing atoms\n");
 451
 452     if (a2)
 453     {
 454         cmp_int(fp, "atoms->nr", -1, a1->nr, a2->nr);
 455         for (i = 0; (i < a1->nr); i++)
 456         {
 457             cmp_atom(fp, i, &(a1->atom[i]), &(a2->atom[i]), ftol, abstol);
 458         }
 459     }
 460     else
 461     {
 462         for (i = 0; (i < a1->nr); i++)
 463         {
 464             cmp_atom(fp, i, &(a1->atom[i]), NULL, ftol, abstol);
 465         }
 466     }
 467 }
 468
 469 static void cmp_top(FILE *fp, t_topology *t1, t_topology *t2, real ftol, real abstol)
 470 {
 471     fprintf(fp, "comparing top\n");
 472     if (t2)
 473     {
 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");
 480     }
 481     else
 482     {
 483         cmp_idef(fp, &(t1->idef), NULL, ftol, abstol);
 484         cmp_atoms(fp, &(t1->atoms), NULL, ftol, abstol);
 485     }
 486 }
 487
 488 static void cmp_groups(FILE *fp, gmx_groups_t *g0, gmx_groups_t *g1,
 489                        int natoms0, int natoms1)
 490 {
 491     int  i, j;
 492     char buf[32];
 493
 494     fprintf(fp, "comparing groups\n");
 495
 496     for (i = 0; i < egcNR; i++)
 497     {
 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)
 501         {
 502             for (j = 0; j < g0->grps[i].nr; j++)
 503             {
 504                 sprintf(buf, "grps[%d].name[%d]", i, j);
 505                 cmp_str(fp, buf, -1,
 506                         *g0->grpname[g0->grps[i].nm_ind[j]],
 507                         *g1->grpname[g1->grps[i].nm_ind[j]]);
 508             }
 509         }
 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))
 513         {
 514             for (j = 0; j < natoms0; j++)
 515             {
 516                 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
 517             }
 518         }
 519     }
 520     /* We have compared the names in the groups lists,
 521      * so we can skip the grpname list comparison.
 522      */
 523 }
 524
 525 static void cmp_rvecs_rmstol(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
 526                              real ftol, real abstol)
 527 {
 528     int    i, m;
 529     double rms;
 530
 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.
 534      */
 535     rms = 0.0;
 536     for (i = 0; (i < n); i++)
 537     {
 538         for (m = 0; m < DIM; m++)
 539         {
 540             rms += x1[i][m]*x1[i][m] + x2[i][m]*x2[i][m];
 541         }
 542     }
 543     rms = sqrt(rms/(2*n*DIM));
 544
 545     /* Convert the relative tolerance into an absolute tolerance */
 546     if (ftol*rms < abstol)
 547     {
 548         abstol = ftol*rms;
 549     }
 550
 551     /* And now do the actual comparision */
 552     for (i = 0; (i < n); i++)
 553     {
 554         cmp_rvec(fp, title, i, x1[i], x2[i], 0.0, abstol);
 555     }
 556 }
 557
 558 static void cmp_rvecs(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
 559                       gmx_bool bRMSD, real ftol, real abstol)
 560 {
 561     int    i, m;
 562     double d, ssd;
 563
 564     if (bRMSD)
 565     {
 566         ssd = 0;
 567         for (i = 0; (i < n); i++)
 568         {
 569             for (m = 0; m < DIM; m++)
 570             {
 571                 d    = x1[i][m] - x2[i][m];
 572                 ssd += d*d;
 573             }
 574         }
 575         fprintf(fp, "%s RMSD %g\n", title, std::sqrt(ssd/n));
 576     }
 577     else
 578     {
 579         cmp_rvecs_rmstol(fp, title, n, x1, x2, ftol, abstol);
 580     }
 581 }
 582
 583 static void cmp_grpopts(FILE *fp, t_grpopts *opt1, t_grpopts *opt2, real ftol, real abstol)
 584 {
 585     int  i, j;
 586     char buf1[256], buf2[256];
 587
 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++)
 593     {
 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])
 601         {
 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++)
 605             {
 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);
 608             }
 609         }
 610     }
 611     if (opt1->ngener == opt2->ngener)
 612     {
 613         for (i = 0; i < opt1->ngener; i++)
 614         {
 615             for (j = i; j < opt1->ngener; j++)
 616             {
 617                 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
 618                 cmp_int(fp, buf1, j,
 619                         opt1->egp_flags[opt1->ngener*i+j],
 620                         opt2->egp_flags[opt1->ngener*i+j]);
 621             }
 622         }
 623     }
 624     for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
 625     {
 626         cmp_rvec(fp, "inputrec->grpopts.acc", i, opt1->acc[i], opt2->acc[i], ftol, abstol);
 627     }
 628     for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
 629     {
 630         cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
 631     }
 632 }
 633
 634 static void cmp_cosines(FILE *fp, const char *s, t_cosines c1[DIM], t_cosines c2[DIM], real ftol, real abstol)
 635 {
 636     int  i, m;
 637     char buf[256];
 638
 639     for (m = 0; (m < DIM); m++)
 640     {
 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++)
 644         {
 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);
 647         }
 648     }
 649 }
 650 static void cmp_pull(FILE *fp)
 651 {
 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");
 653 }
 654
 655 static void cmp_simtempvals(FILE *fp, t_simtemp *simtemp1, t_simtemp *simtemp2, int n_lambda, real ftol, real abstol)
 656 {
 657     int i;
 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++)
 662     {
 663         cmp_real(fp, "inputrec->simtempvals->temperatures", -1, simtemp1->temperatures[i], simtemp2->temperatures[i], ftol, abstol);
 664     }
 665 }
 666
 667 static void cmp_expandedvals(FILE *fp, t_expanded *expand1, t_expanded *expand2, int n_lambda, real ftol, real abstol)
 668 {
 669     int i;
 670
 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);
 673
 674     for (i = 0; i < n_lambda; i++)
 675     {
 676         cmp_real(fp, "inputrec->expandedvals->init_lambda_weights", -1,
 677                  expand1->init_lambda_weights[i], expand2->init_lambda_weights[i], ftol, abstol);
 678     }
 679
 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);
 701 }
 702
 703 static void cmp_fepvals(FILE *fp, t_lambda *fep1, t_lambda *fep2, real ftol, real abstol)
 704 {
 705     int i, j;
 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++)
 711     {
 712         for (j = 0; j < std::min(fep1->n_lambda, fep2->n_lambda); j++)
 713         {
 714             cmp_double(fp, "inputrec->fepvals->all_lambda", -1, fep1->all_lambda[i][j], fep2->all_lambda[i][j], ftol, abstol);
 715         }
 716     }
 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);
 729 }
 730
 731 static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol, real abstol)
 732 {
 733     fprintf(fp, "comparing inputrec\n");
 734
 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)
 741      */
 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);
 812
 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))
 819     {
 820         cmp_simtempvals(fp, ir1->simtempvals, ir2->simtempvals, std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
 821     }
 822     cmp_int(fp, "inputrec->bExpanded", -1, ir1->bExpanded, ir2->bExpanded);
 823     if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
 824     {
 825         cmp_expandedvals(fp, ir1->expandedvals, ir2->expandedvals, std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
 826     }
 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);
 834
 835     cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
 836     if (ir1->bPull && ir2->bPull)
 837     {
 838         cmp_pull(fp);
 839     }
 840
 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);
 866
 867
 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);
 879 }
 880
 881 static void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol)
 882 {
 883     int i;
 884
 885     for (i = 0; i < pull->ncoord; i++)
 886     {
 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);
 889     }
 890 }
 891
 892 static void comp_state(t_state *st1, t_state *st2,
 893                        gmx_bool bRMSD, real ftol, real abstol)
 894 {
 895     int i, j, nc;
 896
 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))
 906     {
 907         fprintf(stdout, "comparing shake vir_prev\n");
 908         cmp_rvecs(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, FALSE, ftol, abstol);
 909     }
 910     if (st1->flags & (1<<estFVIR_PREV))
 911     {
 912         fprintf(stdout, "comparing force vir_prev\n");
 913         cmp_rvecs(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, FALSE, ftol, abstol);
 914     }
 915     if (st1->flags & (1<<estPRES_PREV))
 916     {
 917         fprintf(stdout, "comparing prev_pres\n");
 918         cmp_rvecs(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, FALSE, ftol, abstol);
 919     }
 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)
 923     {
 924         for (i = 0; i < st1->ngtc; i++)
 925         {
 926             nc = i*st1->nhchainlength;
 927             for (j = 0; j < nc; j++)
 928             {
 929                 cmp_real(stdout, "nosehoover_xi",
 930                          i, st1->nosehoover_xi[nc+j], st2->nosehoover_xi[nc+j], ftol, abstol);
 931             }
 932         }
 933     }
 934     cmp_int(stdout, "nnhpres", -1, st1->nnhpres, st2->nnhpres);
 935     if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength)
 936     {
 937         for (i = 0; i < st1->nnhpres; i++)
 938         {
 939             nc = i*st1->nhchainlength;
 940             for (j = 0; j < nc; j++)
 941             {
 942                 cmp_real(stdout, "nosehoover_xi",
 943                          i, st1->nhpres_xi[nc+j], st2->nhpres_xi[nc+j], ftol, abstol);
 944             }
 945         }
 946     }
 947
 948     cmp_int(stdout, "natoms", -1, st1->natoms, st2->natoms);
 949     if (st1->natoms == st2->natoms)
 950     {
 951         if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX)))
 952         {
 953             fprintf(stdout, "comparing x\n");
 954             cmp_rvecs(stdout, "x", st1->natoms, st1->x, st2->x, bRMSD, ftol, abstol);
 955         }
 956         if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV)))
 957         {
 958             fprintf(stdout, "comparing v\n");
 959             cmp_rvecs(stdout, "v", st1->natoms, st1->v, st2->v, bRMSD, ftol, abstol);
 960         }
 961     }
 962 }
 963
 964 void comp_tpx(const char *fn1, const char *fn2,
 965               gmx_bool bRMSD, real ftol, real abstol)
 966 {
 967     const char  *ff[2];
 968     t_inputrec   ir[2];
 969     t_state      state[2];
 970     gmx_mtop_t   mtop[2];
 971     t_topology   top[2];
 972     int          i;
 973
 974     ff[0] = fn1;
 975     ff[1] = fn2;
 976     for (i = 0; i < (fn2 ? 2 : 1); i++)
 977     {
 978         read_tpx_state(ff[i], &(ir[i]), &state[i], &(mtop[i]));
 979     }
 980     if (fn2)
 981     {
 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.
 986          */
 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);
 993     }
 994     else
 995     {
 996         if (ir[0].efep == efepNO)
 997         {
 998             fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0].efep]);
 999         }
1000         else
1001         {
1002             if (ir[0].bPull)
1003             {
1004                 comp_pull_AB(stdout, ir->pull, ftol, abstol);
1005             }
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.
1009              */
1010             top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
1011             cmp_top(stdout, &top[0], NULL, ftol, abstol);
1012         }
1013     }
1014 }
1015
1016 void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
1017                 gmx_bool bRMSD, real ftol, real abstol)
1018 {
1019     fprintf(fp, "\n");
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))
1023     {
1024         cmp_str(fp, "title", -1, fr1->title, fr2->title);
1025     }
1026     if (cmp_bool(fp, "bStep", -1, fr1->bStep, fr2->bStep))
1027     {
1028         cmp_int(fp, "step", -1, fr1->step, fr2->step);
1029     }
1030     cmp_int(fp, "step", -1, fr1->step, fr2->step);
1031     if (cmp_bool(fp, "bTime", -1, fr1->bTime, fr2->bTime))
1032     {
1033         cmp_real(fp, "time", -1, fr1->time, fr2->time, ftol, abstol);
1034     }
1035     if (cmp_bool(fp, "bLambda", -1, fr1->bLambda, fr2->bLambda))
1036     {
1037         cmp_real(fp, "lambda", -1, fr1->lambda, fr2->lambda, ftol, abstol);
1038     }
1039     if (cmp_bool(fp, "bAtoms", -1, fr1->bAtoms, fr2->bAtoms))
1040     {
1041         cmp_atoms(fp, fr1->atoms, fr2->atoms, ftol, abstol);
1042     }
1043     if (cmp_bool(fp, "bPrec", -1, fr1->bPrec, fr2->bPrec))
1044     {
1045         cmp_real(fp, "prec", -1, fr1->prec, fr2->prec, ftol, abstol);
1046     }
1047     if (cmp_bool(fp, "bX", -1, fr1->bX, fr2->bX))
1048     {
1049         cmp_rvecs(fp, "x", std::min(fr1->natoms, fr2->natoms), fr1->x, fr2->x, bRMSD, ftol, abstol);
1050     }
1051     if (cmp_bool(fp, "bV", -1, fr1->bV, fr2->bV))
1052     {
1053         cmp_rvecs(fp, "v", std::min(fr1->natoms, fr2->natoms), fr1->v, fr2->v, bRMSD, ftol, abstol);
1054     }
1055     if (cmp_bool(fp, "bF", -1, fr1->bF, fr2->bF))
1056     {
1057         cmp_rvecs(fp, "f", std::min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, bRMSD, ftol, abstol);
1058     }
1059     if (cmp_bool(fp, "bBox", -1, fr1->bBox, fr2->bBox))
1060     {
1061         cmp_rvecs(fp, "box", 3, fr1->box, fr2->box, FALSE, ftol, abstol);
1062     }
1063 }
1064
1065 void comp_trx(const gmx_output_env_t *oenv, const char *fn1, const char *fn2,
1066               gmx_bool bRMSD, real ftol, real abstol)
1067 {
1068     int          i;
1069     const char  *fn[2];
1070     t_trxframe   fr[2];
1071     t_trxstatus *status[2];
1072     gmx_bool     b[2];
1073
1074     fn[0] = fn1;
1075     fn[1] = fn2;
1076     fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
1077     for (i = 0; i < 2; i++)
1078     {
1079         b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X|TRX_READ_V|TRX_READ_F);
1080     }
1081
1082     if (b[0] && b[1])
1083     {
1084         do
1085         {
1086             comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
1087
1088             for (i = 0; i < 2; i++)
1089             {
1090                 b[i] = read_next_frame(oenv, status[i], &fr[i]);
1091             }
1092         }
1093         while (b[0] && b[1]);
1094
1095         for (i = 0; i < 2; i++)
1096         {
1097             if (b[i] && !b[1-i])
1098             {
1099                 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1-i], fn[i]);
1100             }
1101             close_trj(status[i]);
1102         }
1103     }
1104     if (!b[0] && !b[1])
1105     {
1106         fprintf(stdout, "\nBoth files read correctly\n");
1107     }
1108 }
1109
1110 static real ener_tensor_diag(int n, int *ind1, int *ind2,
1111                              gmx_enxnm_t *enm1,
1112                              int *tensi, int i,
1113                              t_energy e1[], t_energy e2[])
1114 {
1115     int    d1, d2;
1116     int    j;
1117     real   prod1, prod2;
1118     int    nfound;
1119     size_t len;
1120
1121     d1 = tensi[i]/DIM;
1122     d2 = tensi[i] - d1*DIM;
1123
1124     /* Find the diagonal elements d1 and d2 */
1125     len    = std::strlen(enm1[ind1[i]].name);
1126     prod1  = 1;
1127     prod2  = 1;
1128     nfound = 0;
1129     for (j = 0; j < n; j++)
1130     {
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))
1135         {
1136             prod1 *= fabs(e1[ind1[j]].e);
1137             prod2 *= fabs(e2[ind2[j]].e);
1138             nfound++;
1139         }
1140     }
1141
1142     if (nfound == 2)
1143     {
1144         return 0.5*(std::sqrt(prod1) + std::sqrt(prod2));
1145     }
1146     else
1147     {
1148         return 0;
1149     }
1150 }
1151
1152 static gmx_bool enernm_equal(const char *nm1, const char *nm2)
1153 {
1154     int len1, len2;
1155
1156     len1 = std::strlen(nm1);
1157     len2 = std::strlen(nm2);
1158
1159     /* Remove " (bar)" at the end of a name */
1160     if (len1 > 6 && std::strcmp(nm1+len1-6, " (bar)") == 0)
1161     {
1162         len1 -= 6;
1163     }
1164     if (len2 > 6 && std::strcmp(nm2+len2-6, " (bar)") == 0)
1165     {
1166         len2 -= 6;
1167     }
1168
1169     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1170 }
1171
1172 static void cmp_energies(FILE *fp, int step1, int step2,
1173                          t_energy e1[], t_energy e2[],
1174                          gmx_enxnm_t *enm1,
1175                          real ftol, real abstol,
1176                          int nre, int *ind1, int *ind2, int maxener)
1177 {
1178     int   i, ii;
1179     int  *tensi, len, d1, d2;
1180     real  ftol_i, abstol_i;
1181
1182     snew(tensi, maxener);
1183     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1184     for (i = 0; (i < maxener); i++)
1185     {
1186         ii       = ind1[i];
1187         tensi[i] = -1;
1188         len      = std::strlen(enm1[ii].name);
1189         if (len > 3 && enm1[ii].name[len-3] == '-')
1190         {
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)
1195             {
1196                 tensi[i] = d1*DIM + d2;
1197             }
1198         }
1199     }
1200
1201     for (i = 0; (i < maxener); i++)
1202     {
1203         /* Check if this is an off-diagonal tensor element */
1204         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1205         {
1206             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1207             ftol_i = 5;
1208             /* Do the relative tolerance through an absolute tolerance times
1209              * the size of diagonal components of the tensor.
1210              */
1211             abstol_i = ftol*ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1212             if (debug)
1213             {
1214                 fprintf(debug, "tensor '%s' val %f diag %f\n",
1215                         enm1[i].name, e1[i].e, abstol_i/ftol);
1216             }
1217             if (abstol_i > 0)
1218             {
1219                 /* We found a diagonal, we need to check with the minimum tolerance */
1220                 abstol_i = std::min(abstol_i, abstol);
1221             }
1222             else
1223             {
1224                 /* We did not find a diagonal, ignore the relative tolerance check */
1225                 abstol_i = abstol;
1226             }
1227         }
1228         else
1229         {
1230             ftol_i   = ftol;
1231             abstol_i = abstol;
1232         }
1233         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1234         {
1235             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n",
1236                     enm1[ind1[i]].name,
1237                     step1, e1[ind1[i]].e,
1238                     step2, e2[ind2[i]].e);
1239         }
1240     }
1241
1242     sfree(tensi);
1243 }
1244
1245 #if 0
1246 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1247 {
1248     int  i;
1249     char bav[64], bt[64], bs[22];
1250
1251     cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1252     if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1253     {
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++)
1257         {
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);
1260         }
1261     }
1262 }
1263 #endif
1264
1265 static void cmp_eblocks(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1266 {
1267     int  i, j, k;
1268     char buf[64], bs[22];
1269
1270     cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1271     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1272     {
1273         for (j = 0; (j < fr1->nblock); j++)
1274         {
1275             t_enxblock *b1, *b2; /* convenience vars */
1276
1277             b1 = &(fr1->block[j]);
1278             b2 = &(fr2->block[j]);
1279
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);
1283
1284             if ( (b1->nsub == b2->nsub) && (b1->id == b2->id) )
1285             {
1286                 for (i = 0; i < b1->nsub; i++)
1287                 {
1288                     t_enxsubblock *s1, *s2;
1289
1290                     s1 = &(b1->sub[i]);
1291                     s2 = &(b2->sub[i]);
1292
1293                     cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
1294                     cmp_int64(stdout, buf, s1->nr, s2->nr);
1295
1296                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
1297                     {
1298                         switch (s1->type)
1299                         {
1300                             case xdr_datatype_float:
1301                                 for (k = 0; k < s1->nr; k++)
1302                                 {
1303                                     cmp_float(stdout, buf, i,
1304                                               s1->fval[k], s2->fval[k],
1305                                               ftol, abstol);
1306                                 }
1307                                 break;
1308                             case xdr_datatype_double:
1309                                 for (k = 0; k < s1->nr; k++)
1310                                 {
1311                                     cmp_double(stdout, buf, i,
1312                                                s1->dval[k], s2->dval[k],
1313                                                ftol, abstol);
1314                                 }
1315                                 break;
1316                             case xdr_datatype_int:
1317                                 for (k = 0; k < s1->nr; k++)
1318                                 {
1319                                     cmp_int(stdout, buf, i,
1320                                             s1->ival[k], s2->ival[k]);
1321                                 }
1322                                 break;
1323                             case xdr_datatype_int64:
1324                                 for (k = 0; k < s1->nr; k++)
1325                                 {
1326                                     cmp_int64(stdout, buf,
1327                                               s1->lval[k], s2->lval[k]);
1328                                 }
1329                                 break;
1330                             case xdr_datatype_char:
1331                                 for (k = 0; k < s1->nr; k++)
1332                                 {
1333                                     cmp_uc(stdout, buf, i,
1334                                            s1->cval[k], s2->cval[k]);
1335                                 }
1336                                 break;
1337                             case xdr_datatype_string:
1338                                 for (k = 0; k < s1->nr; k++)
1339                                 {
1340                                     cmp_str(stdout, buf, i,
1341                                             s1->sval[k], s2->sval[k]);
1342                                 }
1343                                 break;
1344                             default:
1345                                 gmx_incons("Unknown data type!!");
1346                         }
1347                     }
1348                 }
1349             }
1350         }
1351     }
1352 }
1353
1354 void comp_enx(const char *fn1, const char *fn2, real ftol, real abstol, const char *lastener)
1355 {
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;
1361     gmx_bool       b1, b2;
1362
1363     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1364
1365     in1 = open_enx(fn1, "r");
1366     in2 = open_enx(fn2, "r");
1367     do_enxnms(in1, &nre1, &enm1);
1368     do_enxnms(in2, &nre2, &enm2);
1369     if (nre1 != nre2)
1370     {
1371         fprintf(stdout, "There are %d and %d terms in the energy files\n\n",
1372                 nre1, nre2);
1373     }
1374     else
1375     {
1376         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1377     }
1378
1379     snew(ind1, nre1);
1380     snew(ind2, nre2);
1381     snew(have, nre2);
1382     nre = 0;
1383     for (i = 0; i < nre1; i++)
1384     {
1385         for (j = 0; j < nre2; j++)
1386         {
1387             if (enernm_equal(enm1[i].name, enm2[j].name))
1388             {
1389                 ind1[nre] = i;
1390                 ind2[nre] = j;
1391                 have[j]   = 1;
1392                 nre++;
1393                 break;
1394             }
1395         }
1396         if (nre == 0 || ind1[nre-1] != i)
1397         {
1398             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1399         }
1400     }
1401     for (i = 0; i < nre2; i++)
1402     {
1403         if (have[i] == 0)
1404         {
1405             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1406         }
1407     }
1408
1409     maxener = nre;
1410     for (i = 0; i < nre; i++)
1411     {
1412         if ((lastener != NULL) && (std::strstr(enm1[i].name, lastener) != NULL))
1413         {
1414             maxener = i+1;
1415             break;
1416         }
1417     }
1418
1419     fprintf(stdout, "There are %d terms to compare in the energy files\n\n",
1420             maxener);
1421
1422     for (i = 0; i < maxener; i++)
1423     {
1424         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1425     }
1426
1427     snew(fr1, 1);
1428     snew(fr2, 1);
1429     do
1430     {
1431         b1 = do_enx(in1, fr1);
1432         b2 = do_enx(in2, fr2);
1433         if (b1 && !b2)
1434         {
1435             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1436         }
1437         else if (!b1 && b2)
1438         {
1439             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1440         }
1441         else if (!b1 && !b2)
1442         {
1443             fprintf(stdout, "\nFiles read successfully\n");
1444         }
1445         else
1446         {
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))
1452             {
1453                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener,
1454                              enm1, ftol, abstol, nre, ind1, ind2, maxener);
1455             }
1456             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1457             cmp_eblocks(fr1, fr2, ftol, abstol);
1458         }
1459     }
1460     while (b1 && b2);
1461
1462     close_enx(in1);
1463     close_enx(in2);
1464
1465     free_enxframe(fr2);
1466     sfree(fr2);
1467     free_enxframe(fr1);
1468     sfree(fr1);
1469 }