src/kernel/tpbcmp.c

   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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <string.h>
  46 #include "main.h"
  47 #include "macros.h"
  48 #include "smalloc.h"
  49 #include "futil.h"
  50 #include "statutil.h"
  51 #include "sysstuff.h"
  52 #include "txtdump.h"
  53 #include "gmx_fatal.h"
  54 #include "names.h"
  55 #include "tpxio.h"
  56 #include "enxio.h"
  57 #include "mtop_util.h"
  58 #include "string2.h"
  59 #include "tpbcmp.h"
  60
  61 static void cmp_int(FILE *fp, const char *s, int index, int i1, int i2)
  62 {
  63     if (i1 != i2)
  64     {
  65         if (index != -1)
  66         {
  67             fprintf(fp, "%s[%d] (%d - %d)\n", s, index, i1, i2);
  68         }
  69         else
  70         {
  71             fprintf(fp, "%s (%d - %d)\n", s, i1, i2);
  72         }
  73     }
  74 }
  75
  76 static void cmp_gmx_large_int(FILE *fp, const char *s, gmx_large_int_t i1, gmx_large_int_t i2)
  77 {
  78     if (i1 != i2)
  79     {
  80         fprintf(fp, "%s (", s);
  81         fprintf(fp, gmx_large_int_pfmt, i1);
  82         fprintf(fp, " - ");
  83         fprintf(fp, gmx_large_int_pfmt, i2);
  84         fprintf(fp, ")\n");
  85     }
  86 }
  87
  88 static void cmp_us(FILE *fp, const char *s, int index, unsigned short i1, unsigned short i2)
  89 {
  90     if (i1 != i2)
  91     {
  92         if (index != -1)
  93         {
  94             fprintf(fp, "%s[%d] (%d - %d)\n", s, index, i1, i2);
  95         }
  96         else
  97         {
  98             fprintf(fp, "%s (%d - %d)\n", s, i1, i2);
  99         }
 100     }
 101 }
 102
 103 static void cmp_uc(FILE *fp, const char *s, int index, unsigned char i1, unsigned char i2)
 104 {
 105     if (i1 != i2)
 106     {
 107         if (index != -1)
 108         {
 109             fprintf(fp, "%s[%d] (%d - %d)\n", s, index, i1, i2);
 110         }
 111         else
 112         {
 113             fprintf(fp, "%s (%d - %d)\n", s, i1, i2);
 114         }
 115     }
 116 }
 117
 118 static gmx_bool cmp_bool(FILE *fp, const char *s, int index, gmx_bool b1, gmx_bool b2)
 119 {
 120     if (b1)
 121     {
 122         b1 = 1;
 123     }
 124     else
 125     {
 126         b1 = 0;
 127     }
 128     if (b2)
 129     {
 130         b2 = 1;
 131     }
 132     else
 133     {
 134         b2 = 0;
 135     }
 136     if (b1 != b2)
 137     {
 138         if (index != -1)
 139         {
 140             fprintf(fp, "%s[%d] (%s - %s)\n", s, index,
 141                     bool_names[b1], bool_names[b2]);
 142         }
 143         else
 144         {
 145             fprintf(fp, "%s (%s - %s)\n", s,
 146                     bool_names[b1], bool_names[b2]);
 147         }
 148     }
 149     return b1 && b2;
 150 }
 151
 152 static void cmp_str(FILE *fp, const char *s, int index,
 153                     const char *s1, const char *s2)
 154 {
 155     if (strcmp(s1, s2) != 0)
 156     {
 157         if (index != -1)
 158         {
 159             fprintf(fp, "%s[%d] (%s - %s)\n", s, index, s1, s2);
 160         }
 161         else
 162         {
 163             fprintf(fp, "%s (%s - %s)\n", s, s1, s2);
 164         }
 165     }
 166 }
 167
 168 static gmx_bool equal_real(real i1, real i2, real ftol, real abstol)
 169 {
 170     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 171 }
 172
 173 static gmx_bool equal_float(float i1, float i2, float ftol, float abstol)
 174 {
 175     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 176 }
 177
 178 static gmx_bool equal_double(double i1, double i2, real ftol, real abstol)
 179 {
 180     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2) <= abstol );
 181 }
 182
 183 static void
 184 cmp_real(FILE *fp, const char *s, int index, real i1, real i2, real ftol, real abstol)
 185 {
 186     if (!equal_real(i1, i2, ftol, abstol))
 187     {
 188         if (index != -1)
 189         {
 190             fprintf(fp, "%s[%2d] (%e - %e)\n", s, index, i1, i2);
 191         }
 192         else
 193         {
 194             fprintf(fp, "%s (%e - %e)\n", s, i1, i2);
 195         }
 196     }
 197 }
 198
 199 static void
 200 cmp_float(FILE *fp, const char *s, int index, float i1, float i2, float ftol, float abstol)
 201 {
 202     if (!equal_float(i1, i2, ftol, abstol))
 203     {
 204         if (index != -1)
 205         {
 206             fprintf(fp, "%s[%2d] (%e - %e)\n", s, index, i1, i2);
 207         }
 208         else
 209         {
 210             fprintf(fp, "%s (%e - %e)\n", s, i1, i2);
 211         }
 212     }
 213 }
 214
 215
 216
 217 static void
 218 cmp_double(FILE *fp, const char *s, int index, double i1, double i2, double ftol, double abstol)
 219 {
 220     if (!equal_double(i1, i2, ftol, abstol))
 221     {
 222         if (index != -1)
 223         {
 224             fprintf(fp, "%s[%2d] (%16.9e - %16.9e)\n", s, index, i1, i2);
 225         }
 226         else
 227         {
 228             fprintf(fp, "%s (%16.9e - %16.9e)\n", s, i1, i2);
 229         }
 230     }
 231 }
 232
 233 static void cmp_rvec(FILE *fp, const char *s, int index, rvec i1, rvec i2, real ftol, real abstol)
 234 {
 235     if (!equal_real(i1[XX], i2[XX], ftol, abstol) ||
 236         !equal_real(i1[YY], i2[YY], ftol, abstol) ||
 237         !equal_real(i1[ZZ], i2[ZZ], ftol, abstol))
 238     {
 239         if (index != -1)
 240         {
 241             fprintf(fp, "%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
 242                     s, index, i1[XX], i1[YY], i1[ZZ], i2[XX], i2[YY], i2[ZZ]);
 243         }
 244         else
 245         {
 246             fprintf(fp, "%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
 247                     s, i1[XX], i1[YY], i1[ZZ], i2[XX], i2[YY], i2[ZZ]);
 248         }
 249     }
 250 }
 251
 252 static void cmp_ivec(FILE *fp, const char *s, int index, ivec i1, ivec i2)
 253 {
 254     if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ]))
 255     {
 256         if (index != -1)
 257         {
 258             fprintf(fp, "%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n", s, index,
 259                     i1[XX], i1[YY], i1[ZZ], i2[XX], i2[YY], i2[ZZ]);
 260         }
 261         else
 262         {
 263             fprintf(fp, "%s (%8d,%8d,%8d - %8d,%8d,%8d)\n", s,
 264                     i1[XX], i1[YY], i1[ZZ], i2[XX], i2[YY], i2[ZZ]);
 265         }
 266     }
 267 }
 268
 269 static void cmp_ilist(FILE *fp, int ftype, t_ilist *il1, t_ilist *il2)
 270 {
 271     int  i;
 272     char buf[256];
 273
 274     fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
 275     sprintf(buf, "%s->nr", interaction_function[ftype].name);
 276     cmp_int(fp, buf, -1, il1->nr, il2->nr);
 277     sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
 278     if (((il1->nr > 0) && (!il1->iatoms)) ||
 279         ((il2->nr > 0) && (!il2->iatoms)) ||
 280         ((il1->nr != il2->nr)))
 281     {
 282         fprintf(fp, "Comparing radically different topologies - %s is different\n",
 283                 buf);
 284     }
 285     else
 286     {
 287         for (i = 0; (i < il1->nr); i++)
 288         {
 289             cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
 290         }
 291     }
 292 }
 293
 294 void cmp_iparm(FILE *fp, const char *s, t_functype ft,
 295                t_iparams ip1, t_iparams ip2, real ftol, real abstol)
 296 {
 297     int      i;
 298     gmx_bool bDiff;
 299
 300     bDiff = FALSE;
 301     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 302     {
 303         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
 304     }
 305     if (bDiff)
 306     {
 307         fprintf(fp, "%s1: ", s);
 308         pr_iparams(fp, ft, &ip1);
 309         fprintf(fp, "%s2: ", s);
 310         pr_iparams(fp, ft, &ip2);
 311     }
 312 }
 313
 314 void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft, t_iparams ip1, real ftol, real abstol)
 315 {
 316     int      nrfpA, nrfpB, p0, i;
 317     gmx_bool bDiff;
 318
 319     /* Normally the first parameter is perturbable */
 320     p0    = 0;
 321     nrfpA = interaction_function[ft].nrfpA;
 322     nrfpB = interaction_function[ft].nrfpB;
 323     if (ft == F_PDIHS)
 324     {
 325         nrfpB = 2;
 326     }
 327     else if (interaction_function[ft].flags & IF_TABULATED)
 328     {
 329         /* For tabulated interactions only the second parameter is perturbable */
 330         p0    = 1;
 331         nrfpB = 1;
 332     }
 333     bDiff = FALSE;
 334     for (i = 0; i < nrfpB && !bDiff; i++)
 335     {
 336         bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
 337     }
 338     if (bDiff)
 339     {
 340         fprintf(fp, "%s: ", s);
 341         pr_iparams(fp, ft, &ip1);
 342     }
 343 }
 344
 345 static void cmp_idef(FILE *fp, t_idef *id1, t_idef *id2, real ftol, real abstol)
 346 {
 347     int  i;
 348     char buf1[64], buf2[64];
 349
 350     fprintf(fp, "comparing idef\n");
 351     if (id2)
 352     {
 353         cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
 354         cmp_int(fp, "idef->atnr",  -1, id1->atnr, id2->atnr);
 355         for (i = 0; (i < id1->ntypes); i++)
 356         {
 357             sprintf(buf1, "idef->functype[%d]", i);
 358             sprintf(buf2, "idef->iparam[%d]", i);
 359             cmp_int(fp, buf1, i, (int)id1->functype[i], (int)id2->functype[i]);
 360             cmp_iparm(fp, buf2, id1->functype[i],
 361                       id1->iparams[i], id2->iparams[i], ftol, abstol);
 362         }
 363         cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
 364         for (i = 0; (i < F_NRE); i++)
 365         {
 366             cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
 367         }
 368     }
 369     else
 370     {
 371         for (i = 0; (i < id1->ntypes); i++)
 372         {
 373             cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
 374         }
 375     }
 376 }
 377
 378 static void cmp_block(FILE *fp, t_block *b1, t_block *b2, const char *s)
 379 {
 380     int  i, j, k;
 381     char buf[32];
 382
 383     fprintf(fp, "comparing block %s\n", s);
 384     sprintf(buf, "%s.nr", s);
 385     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 386 }
 387
 388 static void cmp_blocka(FILE *fp, t_blocka *b1, t_blocka *b2, const char *s)
 389 {
 390     int  i, j, k;
 391     char buf[32];
 392
 393     fprintf(fp, "comparing blocka %s\n", s);
 394     sprintf(buf, "%s.nr", s);
 395     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 396     sprintf(buf, "%s.nra", s);
 397     cmp_int(fp, buf, -1, b1->nra, b2->nra);
 398 }
 399
 400 static void cmp_atom(FILE *fp, int index, t_atom *a1, t_atom *a2, real ftol, real abstol)
 401 {
 402     int  i;
 403     char buf[256];
 404
 405     if (a2)
 406     {
 407         cmp_us(fp, "atom.type", index, a1->type, a2->type);
 408         cmp_us(fp, "atom.ptype", index, a1->ptype, a2->ptype);
 409         cmp_int(fp, "atom.resind", index, a1->resind, a2->resind);
 410         cmp_int(fp, "atom.atomnumber", index, a1->atomnumber, a2->atomnumber);
 411         cmp_real(fp, "atom.m", index, a1->m, a2->m, ftol, abstol);
 412         cmp_real(fp, "atom.q", index, a1->q, a2->q, ftol, abstol);
 413         cmp_us(fp, "atom.typeB", index, a1->typeB, a2->typeB);
 414         cmp_real(fp, "atom.mB", index, a1->mB, a2->mB, ftol, abstol);
 415         cmp_real(fp, "atom.qB", index, a1->qB, a2->qB, ftol, abstol);
 416     }
 417     else
 418     {
 419         cmp_us(fp, "atom.type", index, a1->type, a1->typeB);
 420         cmp_real(fp, "atom.m", index, a1->m, a1->mB, ftol, abstol);
 421         cmp_real(fp, "atom.q", index, a1->q, a1->qB, ftol, abstol);
 422     }
 423 }
 424
 425 static void cmp_atoms(FILE *fp, t_atoms *a1, t_atoms *a2, real ftol, real abstol)
 426 {
 427     int i;
 428
 429     fprintf(fp, "comparing atoms\n");
 430
 431     if (a2)
 432     {
 433         cmp_int(fp, "atoms->nr", -1, a1->nr, a2->nr);
 434         for (i = 0; (i < a1->nr); i++)
 435         {
 436             cmp_atom(fp, i, &(a1->atom[i]), &(a2->atom[i]), ftol, abstol);
 437         }
 438     }
 439     else
 440     {
 441         for (i = 0; (i < a1->nr); i++)
 442         {
 443             cmp_atom(fp, i, &(a1->atom[i]), NULL, ftol, abstol);
 444         }
 445     }
 446 }
 447
 448 static void cmp_top(FILE *fp, t_topology *t1, t_topology *t2, real ftol, real abstol)
 449 {
 450     int i;
 451
 452     fprintf(fp, "comparing top\n");
 453     if (t2)
 454     {
 455         cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
 456         cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
 457         cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
 458         cmp_block(fp, &t1->mols, &t2->mols, "mols");
 459         cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
 460     }
 461     else
 462     {
 463         cmp_idef(fp, &(t1->idef), NULL, ftol, abstol);
 464         cmp_atoms(fp, &(t1->atoms), NULL, ftol, abstol);
 465     }
 466 }
 467
 468 static void cmp_groups(FILE *fp, gmx_groups_t *g0, gmx_groups_t *g1,
 469                        int natoms0, int natoms1)
 470 {
 471     int  i, j, ndiff;
 472     char buf[32];
 473
 474     fprintf(fp, "comparing groups\n");
 475
 476     for (i = 0; i < egcNR; i++)
 477     {
 478         sprintf(buf, "grps[%d].nr", i);
 479         cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
 480         if (g0->grps[i].nr == g1->grps[i].nr)
 481         {
 482             for (j = 0; j < g0->grps[i].nr; j++)
 483             {
 484                 sprintf(buf, "grps[%d].name[%d]", i, j);
 485                 cmp_str(fp, buf, -1,
 486                         *g0->grpname[g0->grps[i].nm_ind[j]],
 487                         *g1->grpname[g1->grps[i].nm_ind[j]]);
 488             }
 489         }
 490         cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
 491         if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
 492             (g0->grpnr[i] != NULL || g1->grpnr[i] != NULL))
 493         {
 494             for (j = 0; j < natoms0; j++)
 495             {
 496                 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
 497             }
 498         }
 499     }
 500     /* We have compared the names in the groups lists,
 501      * so we can skip the grpname list comparison.
 502      */
 503 }
 504
 505 static void cmp_rvecs(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
 506                       gmx_bool bRMSD, real ftol, real abstol)
 507 {
 508     int    i, m;
 509     double d, ssd;
 510
 511     if (bRMSD)
 512     {
 513         ssd = 0;
 514         for (i = 0; (i < n); i++)
 515         {
 516             for (m = 0; m < DIM; m++)
 517             {
 518                 d    = x1[i][m] - x2[i][m];
 519                 ssd += d*d;
 520             }
 521         }
 522         fprintf(fp, "%s RMSD %g\n", title, sqrt(ssd/n));
 523     }
 524     else
 525     {
 526         for (i = 0; (i < n); i++)
 527         {
 528             cmp_rvec(fp, title, i, x1[i], x2[i], ftol, abstol);
 529         }
 530     }
 531 }
 532
 533
 534 /* Similar to cmp_rvecs, but this routine scales the allowed absolute tolerance
 535  * by the RMS of the force components of x1.
 536  */
 537 static void cmp_rvecs_rmstol(FILE *fp, const char *title, int n, rvec x1[], rvec x2[],
 538                              real ftol, real abstol)
 539 {
 540     int    i, m;
 541     double d;
 542     double ave_x1, rms_x1;
 543
 544     /* It is tricky to compare real values, in particular forces that
 545      * are sums of lots of terms where the final value might be close to 0.0.
 546      * To get a reference magnitude we calculate the RMS value of each
 547      * component in x1, and then set the allowed absolute tolerance to the
 548      * relative tolerance times this RMS magnitude.
 549      */
 550     ave_x1 = 0.0;
 551     for (i = 0; i < n; i++)
 552     {
 553         for (m = 0; m < DIM; m++)
 554         {
 555             ave_x1 += x1[i][m];
 556         }
 557     }
 558     ave_x1 /= n*DIM;
 559
 560     rms_x1 = 0.0;
 561     for (i = 0; (i < n); i++)
 562     {
 563         for (m = 0; m < DIM; m++)
 564         {
 565             d       = x1[i][m] - ave_x1;
 566             rms_x1 += d*d;
 567         }
 568     }
 569     rms_x1 = sqrt(rms_x1/(DIM*n));
 570     /* And now do the actual comparision with a hopefully realistic abstol. */
 571     for (i = 0; (i < n); i++)
 572     {
 573         cmp_rvec(fp, title, i, x1[i], x2[i], ftol, abstol*rms_x1);
 574     }
 575 }
 576
 577 static void cmp_grpopts(FILE *fp, t_grpopts *opt1, t_grpopts *opt2, real ftol, real abstol)
 578 {
 579     int  i, j;
 580     char buf1[256], buf2[256];
 581
 582     cmp_int(fp, "inputrec->grpopts.ngtc", -1,  opt1->ngtc, opt2->ngtc);
 583     cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
 584     cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
 585     cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
 586     for (i = 0; (i < min(opt1->ngtc, opt2->ngtc)); i++)
 587     {
 588         cmp_real(fp, "inputrec->grpopts.nrdf", i, opt1->nrdf[i], opt2->nrdf[i], ftol, abstol);
 589         cmp_real(fp, "inputrec->grpopts.ref_t", i, opt1->ref_t[i], opt2->ref_t[i], ftol, abstol);
 590         cmp_real(fp, "inputrec->grpopts.tau_t", i, opt1->tau_t[i], opt2->tau_t[i], ftol, abstol);
 591         cmp_int(fp, "inputrec->grpopts.annealing", i, opt1->annealing[i], opt2->annealing[i]);
 592         cmp_int(fp, "inputrec->grpopts.anneal_npoints", i,
 593                 opt1->anneal_npoints[i], opt2->anneal_npoints[i]);
 594         if (opt1->anneal_npoints[i] == opt2->anneal_npoints[i])
 595         {
 596             sprintf(buf1, "inputrec->grpopts.anneal_time[%d]", i);
 597             sprintf(buf2, "inputrec->grpopts.anneal_temp[%d]", i);
 598             for (j = 0; j < opt1->anneal_npoints[i]; j++)
 599             {
 600                 cmp_real(fp, buf1, j, opt1->anneal_time[i][j], opt2->anneal_time[i][j], ftol, abstol);
 601                 cmp_real(fp, buf2, j, opt1->anneal_temp[i][j], opt2->anneal_temp[i][j], ftol, abstol);
 602             }
 603         }
 604     }
 605     if (opt1->ngener == opt2->ngener)
 606     {
 607         for (i = 0; i < opt1->ngener; i++)
 608         {
 609             for (j = i; j < opt1->ngener; j++)
 610             {
 611                 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
 612                 cmp_int(fp, buf1, j,
 613                         opt1->egp_flags[opt1->ngener*i+j],
 614                         opt2->egp_flags[opt1->ngener*i+j]);
 615             }
 616         }
 617     }
 618     for (i = 0; (i < min(opt1->ngacc, opt2->ngacc)); i++)
 619     {
 620         cmp_rvec(fp, "inputrec->grpopts.acc", i, opt1->acc[i], opt2->acc[i], ftol, abstol);
 621     }
 622     for (i = 0; (i < min(opt1->ngfrz, opt2->ngfrz)); i++)
 623     {
 624         cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
 625     }
 626 }
 627
 628 static void cmp_cosines(FILE *fp, const char *s, t_cosines c1[DIM], t_cosines c2[DIM], real ftol, real abstol)
 629 {
 630     int  i, m;
 631     char buf[256];
 632
 633     for (m = 0; (m < DIM); m++)
 634     {
 635         sprintf(buf, "inputrec->%s[%d]", s, m);
 636         cmp_int(fp, buf, 0, c1->n, c2->n);
 637         for (i = 0; (i < min(c1->n, c2->n)); i++)
 638         {
 639             cmp_real(fp, buf, i, c1->a[i], c2->a[i], ftol, abstol);
 640             cmp_real(fp, buf, i, c1->phi[i], c2->phi[i], ftol, abstol);
 641         }
 642     }
 643 }
 644 static void cmp_adress(FILE *fp, t_adress *ad1, t_adress *ad2,
 645                        real ftol, real abstol)
 646 {
 647     cmp_int(fp, "ir->adress->type", -1, ad1->type, ad2->type);
 648     cmp_real(fp, "ir->adress->const_wf", -1, ad1->const_wf, ad2->const_wf, ftol, abstol);
 649     cmp_real(fp, "ir->adress->ex_width", -1, ad1->ex_width, ad2->ex_width, ftol, abstol);
 650     cmp_real(fp, "ir->adress->hy_width", -1, ad1->hy_width, ad2->hy_width, ftol, abstol);
 651     cmp_int(fp, "ir->adress->icor", -1, ad1->icor, ad2->icor);
 652     cmp_int(fp, "ir->adress->site", -1, ad1->site, ad2->site);
 653     cmp_rvec(fp, "ir->adress->refs", -1, ad1->refs, ad2->refs, ftol, abstol);
 654     cmp_real(fp, "ir->adress->ex_forcecap", -1, ad1->ex_forcecap, ad2->ex_forcecap, ftol, abstol);
 655 }
 656
 657 static void cmp_pull(FILE *fp, t_pull *pull1, t_pull *pull2, real ftol, real abstol)
 658 {
 659     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");
 660 }
 661
 662 static void cmp_simtempvals(FILE *fp, t_simtemp *simtemp1, t_simtemp *simtemp2, int n_lambda, real ftol, real abstol)
 663 {
 664     int i;
 665     cmp_int(fp, "inputrec->simtempvals->eSimTempScale", -1, simtemp1->eSimTempScale, simtemp2->eSimTempScale);
 666     cmp_real(fp, "inputrec->simtempvals->simtemp_high", -1, simtemp1->simtemp_high, simtemp2->simtemp_high, ftol, abstol);
 667     cmp_real(fp, "inputrec->simtempvals->simtemp_low", -1, simtemp1->simtemp_low, simtemp2->simtemp_low, ftol, abstol);
 668     for (i = 0; i < n_lambda; i++)
 669     {
 670         cmp_real(fp, "inputrec->simtempvals->temperatures", -1, simtemp1->temperatures[i], simtemp2->temperatures[i], ftol, abstol);
 671     }
 672 }
 673
 674 static void cmp_expandedvals(FILE *fp, t_expanded *expand1, t_expanded *expand2, int n_lambda, real ftol, real abstol)
 675 {
 676     int i;
 677
 678     cmp_bool(fp, "inputrec->fepvals->bInit_weights", -1, expand1->bInit_weights, expand2->bInit_weights);
 679     cmp_bool(fp, "inputrec->fepvals->bWLoneovert", -1, expand1->bWLoneovert, expand2->bWLoneovert);
 680
 681     for (i = 0; i < n_lambda; i++)
 682     {
 683         cmp_real(fp, "inputrec->expandedvals->init_lambda_weights", -1,
 684                  expand1->init_lambda_weights[i], expand2->init_lambda_weights[i], ftol, abstol);
 685     }
 686
 687     cmp_int(fp, "inputrec->expandedvals->lambda-stats", -1, expand1->elamstats, expand2->elamstats);
 688     cmp_int(fp, "inputrec->expandedvals->lambda-mc-move", -1, expand1->elmcmove, expand2->elmcmove);
 689     cmp_int(fp, "inputrec->expandedvals->lmc-repeats", -1, expand1->lmc_repeats, expand2->lmc_repeats);
 690     cmp_int(fp, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1->gibbsdeltalam, expand2->gibbsdeltalam);
 691     cmp_int(fp, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1->lmc_forced_nstart, expand2->lmc_forced_nstart);
 692     cmp_int(fp, "inputrec->expandedvals->lambda-weights-equil", -1, expand1->elmceq, expand2->elmceq);
 693     cmp_int(fp, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1, expand1->equil_n_at_lam, expand2->equil_n_at_lam);
 694     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1->equil_samples, expand2->equil_samples);
 695     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1->equil_steps, expand2->equil_steps);
 696     cmp_real(fp, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1->equil_wl_delta, expand2->equil_wl_delta, ftol, abstol);
 697     cmp_real(fp, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1->equil_ratio, expand2->equil_ratio, ftol, abstol);
 698     cmp_bool(fp, "inputrec->expandedvals->symmetrized-transition-matrix", -1, expand1->bSymmetrizedTMatrix, expand2->bSymmetrizedTMatrix);
 699     cmp_int(fp, "inputrec->expandedvals->nstTij", -1, expand1->nstTij, expand2->nstTij);
 700     cmp_int(fp, "inputrec->expandedvals->mininum-var-min", -1, expand1->minvarmin, expand2->minvarmin); /*default is reasonable */
 701     cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range);      /* default is just C=0 */
 702     cmp_real(fp, "inputrec->expandedvals->wl-scale", -1, expand1->wl_scale, expand2->wl_scale, ftol, abstol);
 703     cmp_real(fp, "inputrec->expandedvals->init-wl-delta", -1, expand1->init_wl_delta, expand2->init_wl_delta, ftol, abstol);
 704     cmp_real(fp, "inputrec->expandedvals->wl-ratio", -1, expand1->wl_ratio, expand2->wl_ratio, ftol, abstol);
 705     cmp_int(fp, "inputrec->expandedvals->nstexpanded", -1, expand1->nstexpanded, expand2->nstexpanded);
 706     cmp_int(fp, "inputrec->expandedvals->lmc-seed", -1, expand1->lmc_seed, expand2->lmc_seed);
 707     cmp_real(fp, "inputrec->expandedvals->mc-temperature", -1, expand1->mc_temp, expand2->mc_temp, ftol, abstol);
 708 }
 709
 710 static void cmp_fepvals(FILE *fp, t_lambda *fep1, t_lambda *fep2, real ftol, real abstol)
 711 {
 712     int i, j;
 713     cmp_int(fp, "inputrec->nstdhdl", -1, fep1->nstdhdl, fep2->nstdhdl);
 714     cmp_double(fp, "inputrec->fepvals->init_fep_state", -1, fep1->init_fep_state, fep2->init_fep_state, ftol, abstol);
 715     cmp_double(fp, "inputrec->fepvals->delta_lambda", -1, fep1->delta_lambda, fep2->delta_lambda, ftol, abstol);
 716     cmp_int(fp, "inputrec->fepvals->n_lambda", -1, fep1->n_lambda, fep2->n_lambda);
 717     for (i = 0; i < efptNR; i++)
 718     {
 719         for (j = 0; j < min(fep1->n_lambda, fep2->n_lambda); j++)
 720         {
 721             cmp_double(fp, "inputrec->fepvals->all_lambda", -1, fep1->all_lambda[i][j], fep2->all_lambda[i][j], ftol, abstol);
 722         }
 723     }
 724     cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors,
 725             fep2->lambda_neighbors);
 726     cmp_real(fp, "inputrec->fepvals->sc_alpha", -1, fep1->sc_alpha, fep2->sc_alpha, ftol, abstol);
 727     cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
 728     cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
 729     cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
 730     cmp_bool(fp, "inputrec->fepvals->bPrintEnergy", -1, fep1->bPrintEnergy, fep1->bPrintEnergy);
 731     cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
 732     cmp_int(fp, "inputrec->separate_dhdl_file", -1, fep1->separate_dhdl_file, fep2->separate_dhdl_file);
 733     cmp_int(fp, "inputrec->dhdl_derivatives", -1, fep1->dhdl_derivatives, fep2->dhdl_derivatives);
 734     cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size);
 735     cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing, ftol, abstol);
 736 }
 737
 738 static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol, real abstol)
 739 {
 740     fprintf(fp, "comparing inputrec\n");
 741
 742     /* gcc 2.96 doesnt like these defines at all, but issues a huge list
 743      * of warnings. Maybe it will change in future versions, but for the
 744      * moment I've spelled them out instead. /EL 000820
 745      * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 746      * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 747      * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
 748      */
 749     cmp_int(fp, "inputrec->eI", -1, ir1->eI, ir2->eI);
 750     cmp_gmx_large_int(fp, "inputrec->nsteps", ir1->nsteps, ir2->nsteps);
 751     cmp_gmx_large_int(fp, "inputrec->init_step", ir1->init_step, ir2->init_step);
 752     cmp_int(fp, "inputrec->simulation_part", -1, ir1->simulation_part, ir2->simulation_part);
 753     cmp_int(fp, "inputrec->ePBC", -1, ir1->ePBC, ir2->ePBC);
 754     cmp_int(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
 755     cmp_int(fp, "inputrec->cutoff_scheme", -1, ir1->cutoff_scheme, ir2->cutoff_scheme);
 756     cmp_int(fp, "inputrec->ns_type", -1, ir1->ns_type, ir2->ns_type);
 757     cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
 758     cmp_int(fp, "inputrec->ndelta", -1, ir1->ndelta, ir2->ndelta);
 759     cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
 760     cmp_int(fp, "inputrec->comm_mode", -1, ir1->comm_mode, ir2->comm_mode);
 761     cmp_int(fp, "inputrec->nstcheckpoint", -1, ir1->nstcheckpoint, ir2->nstcheckpoint);
 762     cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
 763     cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
 764     cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
 765     cmp_int(fp, "inputrec->nstfout", -1, ir1->nstfout, ir2->nstfout);
 766     cmp_int(fp, "inputrec->nstcalcenergy", -1, ir1->nstcalcenergy, ir2->nstcalcenergy);
 767     cmp_int(fp, "inputrec->nstenergy", -1, ir1->nstenergy, ir2->nstenergy);
 768     cmp_int(fp, "inputrec->nstxtcout", -1, ir1->nstxtcout, ir2->nstxtcout);
 769     cmp_double(fp, "inputrec->init_t", -1, ir1->init_t, ir2->init_t, ftol, abstol);
 770     cmp_double(fp, "inputrec->delta_t", -1, ir1->delta_t, ir2->delta_t, ftol, abstol);
 771     cmp_real(fp, "inputrec->xtcprec", -1, ir1->xtcprec, ir2->xtcprec, ftol, abstol);
 772     cmp_real(fp, "inputrec->fourierspacing", -1, ir1->fourier_spacing, ir2->fourier_spacing, ftol, abstol);
 773     cmp_int(fp, "inputrec->nkx", -1, ir1->nkx, ir2->nkx);
 774     cmp_int(fp, "inputrec->nky", -1, ir1->nky, ir2->nky);
 775     cmp_int(fp, "inputrec->nkz", -1, ir1->nkz, ir2->nkz);
 776     cmp_int(fp, "inputrec->pme_order", -1, ir1->pme_order, ir2->pme_order);
 777     cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
 778     cmp_int(fp, "inputrec->ewald_geometry", -1, ir1->ewald_geometry, ir2->ewald_geometry);
 779     cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
 780     cmp_int(fp, "inputrec->bOptFFT", -1, ir1->bOptFFT, ir2->bOptFFT);
 781     cmp_int(fp, "inputrec->bContinuation", -1, ir1->bContinuation, ir2->bContinuation);
 782     cmp_int(fp, "inputrec->bShakeSOR", -1, ir1->bShakeSOR, ir2->bShakeSOR);
 783     cmp_int(fp, "inputrec->etc", -1, ir1->etc, ir2->etc);
 784     cmp_int(fp, "inputrec->bPrintNHChains", -1, ir1->bPrintNHChains, ir2->bPrintNHChains);
 785     cmp_int(fp, "inputrec->epc", -1, ir1->epc, ir2->epc);
 786     cmp_int(fp, "inputrec->epct", -1, ir1->epct, ir2->epct);
 787     cmp_real(fp, "inputrec->tau_p", -1, ir1->tau_p, ir2->tau_p, ftol, abstol);
 788     cmp_rvec(fp, "inputrec->ref_p(x)", -1, ir1->ref_p[XX], ir2->ref_p[XX], ftol, abstol);
 789     cmp_rvec(fp, "inputrec->ref_p(y)", -1, ir1->ref_p[YY], ir2->ref_p[YY], ftol, abstol);
 790     cmp_rvec(fp, "inputrec->ref_p(z)", -1, ir1->ref_p[ZZ], ir2->ref_p[ZZ], ftol, abstol);
 791     cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
 792     cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
 793     cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
 794     cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
 795     cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
 796     cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
 797     cmp_real(fp, "inputrec->verletbuf_drift", -1, ir1->verletbuf_drift, ir2->verletbuf_drift, ftol, abstol);
 798     cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
 799     cmp_real(fp, "inputrec->rlistlong", -1, ir1->rlistlong, ir2->rlistlong, ftol, abstol);
 800     cmp_int(fp, "inputrec->nstcalclr", -1, ir1->nstcalclr, ir2->nstcalclr);
 801     cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
 802     cmp_int(fp, "inputrec->coulombtype", -1, ir1->coulombtype, ir2->coulombtype);
 803     cmp_int(fp, "inputrec->coulomb_modifier", -1, ir1->coulomb_modifier, ir2->coulomb_modifier);
 804     cmp_real(fp, "inputrec->rcoulomb_switch", -1, ir1->rcoulomb_switch, ir2->rcoulomb_switch, ftol, abstol);
 805     cmp_real(fp, "inputrec->rcoulomb", -1, ir1->rcoulomb, ir2->rcoulomb, ftol, abstol);
 806     cmp_int(fp, "inputrec->vdwtype", -1, ir1->vdwtype, ir2->vdwtype);
 807     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);
 808     cmp_real(fp, "inputrec->rvdw", -1, ir1->rvdw, ir2->rvdw, ftol, abstol);
 809     cmp_real(fp, "inputrec->epsilon_r", -1, ir1->epsilon_r, ir2->epsilon_r, ftol, abstol);
 810     cmp_real(fp, "inputrec->epsilon_rf", -1, ir1->epsilon_rf, ir2->epsilon_rf, ftol, abstol);
 811     cmp_real(fp, "inputrec->tabext", -1, ir1->tabext, ir2->tabext, ftol, abstol);
 812     cmp_int(fp, "inputrec->implicit_solvent", -1, ir1->implicit_solvent, ir2->implicit_solvent);
 813     cmp_int(fp, "inputrec->gb_algorithm", -1, ir1->gb_algorithm, ir2->gb_algorithm);
 814     cmp_int(fp, "inputrec->nstgbradii", -1, ir1->nstgbradii, ir2->nstgbradii);
 815     cmp_real(fp, "inputrec->rgbradii", -1, ir1->rgbradii, ir2->rgbradii, ftol, abstol);
 816     cmp_real(fp, "inputrec->gb_saltconc", -1, ir1->gb_saltconc, ir2->gb_saltconc, ftol, abstol);
 817     cmp_real(fp, "inputrec->gb_epsilon_solvent", -1, ir1->gb_epsilon_solvent, ir2->gb_epsilon_solvent, ftol, abstol);
 818     cmp_real(fp, "inputrec->gb_obc_alpha", -1, ir1->gb_obc_alpha, ir2->gb_obc_alpha, ftol, abstol);
 819     cmp_real(fp, "inputrec->gb_obc_beta", -1, ir1->gb_obc_beta, ir2->gb_obc_beta, ftol, abstol);
 820     cmp_real(fp, "inputrec->gb_obc_gamma", -1, ir1->gb_obc_gamma, ir2->gb_obc_gamma, ftol, abstol);
 821     cmp_real(fp, "inputrec->gb_dielectric_offset", -1, ir1->gb_dielectric_offset, ir2->gb_dielectric_offset, ftol, abstol);
 822     cmp_int(fp, "inputrec->sa_algorithm", -1, ir1->sa_algorithm, ir2->sa_algorithm);
 823     cmp_real(fp, "inputrec->sa_surface_tension", -1, ir1->sa_surface_tension, ir2->sa_surface_tension, ftol, abstol);
 824
 825     cmp_int(fp, "inputrec->eDispCorr", -1, ir1->eDispCorr, ir2->eDispCorr);
 826     cmp_real(fp, "inputrec->shake_tol", -1, ir1->shake_tol, ir2->shake_tol, ftol, abstol);
 827     cmp_int(fp, "inputrec->efep", -1, ir1->efep, ir2->efep);
 828     cmp_fepvals(fp, ir1->fepvals, ir2->fepvals, ftol, abstol);
 829     cmp_int(fp, "inputrec->bSimTemp", -1, ir1->bSimTemp, ir2->bSimTemp);
 830     if ((ir1->bSimTemp == ir2->bSimTemp) && (ir1->bSimTemp))
 831     {
 832         cmp_simtempvals(fp, ir1->simtempvals, ir2->simtempvals, min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
 833     }
 834     cmp_int(fp, "inputrec->bExpanded", -1, ir1->bExpanded, ir2->bExpanded);
 835     if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
 836     {
 837         cmp_expandedvals(fp, ir1->expandedvals, ir2->expandedvals, min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
 838     }
 839     cmp_int(fp, "inputrec->nwall", -1, ir1->nwall, ir2->nwall);
 840     cmp_int(fp, "inputrec->wall_type", -1, ir1->wall_type, ir2->wall_type);
 841     cmp_int(fp, "inputrec->wall_atomtype[0]", -1, ir1->wall_atomtype[0], ir2->wall_atomtype[0]);
 842     cmp_int(fp, "inputrec->wall_atomtype[1]", -1, ir1->wall_atomtype[1], ir2->wall_atomtype[1]);
 843     cmp_real(fp, "inputrec->wall_density[0]", -1, ir1->wall_density[0], ir2->wall_density[0], ftol, abstol);
 844     cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
 845     cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
 846
 847     cmp_int(fp, "inputrec->ePull", -1, ir1->ePull, ir2->ePull);
 848     if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
 849     {
 850         cmp_pull(fp, ir1->pull, ir2->pull, ftol, abstol);
 851     }
 852
 853     cmp_int(fp, "inputrec->eDisre", -1, ir1->eDisre, ir2->eDisre);
 854     cmp_real(fp, "inputrec->dr_fc", -1, ir1->dr_fc, ir2->dr_fc, ftol, abstol);
 855     cmp_int(fp, "inputrec->eDisreWeighting", -1, ir1->eDisreWeighting, ir2->eDisreWeighting);
 856     cmp_int(fp, "inputrec->bDisreMixed", -1, ir1->bDisreMixed, ir2->bDisreMixed);
 857     cmp_int(fp, "inputrec->nstdisreout", -1, ir1->nstdisreout, ir2->nstdisreout);
 858     cmp_real(fp, "inputrec->dr_tau", -1, ir1->dr_tau, ir2->dr_tau, ftol, abstol);
 859     cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
 860     cmp_real(fp, "inputrec->orires_tau", -1, ir1->orires_tau, ir2->orires_tau, ftol, abstol);
 861     cmp_int(fp, "inputrec->nstorireout", -1, ir1->nstorireout, ir2->nstorireout);
 862     cmp_real(fp, "inputrec->dihre_fc", -1, ir1->dihre_fc, ir2->dihre_fc, ftol, abstol);
 863     cmp_real(fp, "inputrec->em_stepsize", -1, ir1->em_stepsize, ir2->em_stepsize, ftol, abstol);
 864     cmp_real(fp, "inputrec->em_tol", -1, ir1->em_tol, ir2->em_tol, ftol, abstol);
 865     cmp_int(fp, "inputrec->niter", -1, ir1->niter, ir2->niter);
 866     cmp_real(fp, "inputrec->fc_stepsize", -1, ir1->fc_stepsize, ir2->fc_stepsize, ftol, abstol);
 867     cmp_int(fp, "inputrec->nstcgsteep", -1, ir1->nstcgsteep, ir2->nstcgsteep);
 868     cmp_int(fp, "inputrec->nbfgscorr", 0, ir1->nbfgscorr, ir2->nbfgscorr);
 869     cmp_int(fp, "inputrec->eConstrAlg", -1, ir1->eConstrAlg, ir2->eConstrAlg);
 870     cmp_int(fp, "inputrec->nProjOrder", -1, ir1->nProjOrder, ir2->nProjOrder);
 871     cmp_real(fp, "inputrec->LincsWarnAngle", -1, ir1->LincsWarnAngle, ir2->LincsWarnAngle, ftol, abstol);
 872     cmp_int(fp, "inputrec->nLincsIter", -1, ir1->nLincsIter, ir2->nLincsIter);
 873     cmp_real(fp, "inputrec->bd_fric", -1, ir1->bd_fric, ir2->bd_fric, ftol, abstol);
 874     cmp_int(fp, "inputrec->ld_seed", -1, ir1->ld_seed, ir2->ld_seed);
 875     cmp_real(fp, "inputrec->cos_accel", -1, ir1->cos_accel, ir2->cos_accel, ftol, abstol);
 876     cmp_rvec(fp, "inputrec->deform(a)", -1, ir1->deform[XX], ir2->deform[XX], ftol, abstol);
 877     cmp_rvec(fp, "inputrec->deform(b)", -1, ir1->deform[YY], ir2->deform[YY], ftol, abstol);
 878     cmp_rvec(fp, "inputrec->deform(c)", -1, ir1->deform[ZZ], ir2->deform[ZZ], ftol, abstol);
 879
 880
 881     cmp_bool(fp, "ir->bAdress->type", -1, ir1->bAdress, ir2->bAdress);
 882     if (ir1->bAdress && ir2->bAdress)
 883     {
 884         cmp_adress(fp, ir1->adress, ir2->adress, ftol, abstol);
 885     }
 886
 887     cmp_int(fp, "inputrec->userint1", -1, ir1->userint1, ir2->userint1);
 888     cmp_int(fp, "inputrec->userint2", -1, ir1->userint2, ir2->userint2);
 889     cmp_int(fp, "inputrec->userint3", -1, ir1->userint3, ir2->userint3);
 890     cmp_int(fp, "inputrec->userint4", -1, ir1->userint4, ir2->userint4);
 891     cmp_real(fp, "inputrec->userreal1", -1, ir1->userreal1, ir2->userreal1, ftol, abstol);
 892     cmp_real(fp, "inputrec->userreal2", -1, ir1->userreal2, ir2->userreal2, ftol, abstol);
 893     cmp_real(fp, "inputrec->userreal3", -1, ir1->userreal3, ir2->userreal3, ftol, abstol);
 894     cmp_real(fp, "inputrec->userreal4", -1, ir1->userreal4, ir2->userreal4, ftol, abstol);
 895     cmp_grpopts(fp, &(ir1->opts), &(ir2->opts), ftol, abstol);
 896     cmp_cosines(fp, "ex", ir1->ex, ir2->ex, ftol, abstol);
 897     cmp_cosines(fp, "et", ir1->et, ir2->et, ftol, abstol);
 898 }
 899
 900 static void comp_pull_AB(FILE *fp, t_pull *pull, real ftol, real abstol)
 901 {
 902     int i;
 903
 904     for (i = 0; i < pull->ngrp+1; i++)
 905     {
 906         fprintf(fp, "comparing pull group %d\n", i);
 907         cmp_real(fp, "pullgrp->k", -1, pull->grp[i].k, pull->grp[i].kB, ftol, abstol);
 908     }
 909 }
 910
 911 static void comp_state(t_state *st1, t_state *st2,
 912                        gmx_bool bRMSD, real ftol, real abstol)
 913 {
 914     int i, j, nc;
 915
 916     fprintf(stdout, "comparing flags\n");
 917     cmp_int(stdout, "flags", -1, st1->flags, st2->flags);
 918     fprintf(stdout, "comparing box\n");
 919     cmp_rvecs(stdout, "box", DIM, st1->box, st2->box, FALSE, ftol, abstol);
 920     fprintf(stdout, "comparing box_rel\n");
 921     cmp_rvecs(stdout, "box_rel", DIM, st1->box_rel, st2->box_rel, FALSE, ftol, abstol);
 922     fprintf(stdout, "comparing boxv\n");
 923     cmp_rvecs(stdout, "boxv", DIM, st1->boxv, st2->boxv, FALSE, ftol, abstol);
 924     if (st1->flags & (1<<estSVIR_PREV))
 925     {
 926         fprintf(stdout, "comparing shake vir_prev\n");
 927         cmp_rvecs_rmstol(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, ftol, abstol);
 928     }
 929     if (st1->flags & (1<<estFVIR_PREV))
 930     {
 931         fprintf(stdout, "comparing force vir_prev\n");
 932         cmp_rvecs_rmstol(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, ftol, abstol);
 933     }
 934     if (st1->flags & (1<<estPRES_PREV))
 935     {
 936         fprintf(stdout, "comparing prev_pres\n");
 937         cmp_rvecs_rmstol(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, ftol, abstol);
 938     }
 939     cmp_int(stdout, "ngtc", -1, st1->ngtc, st2->ngtc);
 940     cmp_int(stdout, "nhchainlength", -1, st1->nhchainlength, st2->nhchainlength);
 941     if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength)
 942     {
 943         for (i = 0; i < st1->ngtc; i++)
 944         {
 945             nc = i*st1->nhchainlength;
 946             for (j = 0; j < nc; j++)
 947             {
 948                 cmp_real(stdout, "nosehoover_xi",
 949                          i, st1->nosehoover_xi[nc+j], st2->nosehoover_xi[nc+j], ftol, abstol);
 950             }
 951         }
 952     }
 953     cmp_int(stdout, "nnhpres", -1, st1->nnhpres, st2->nnhpres);
 954     if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength)
 955     {
 956         for (i = 0; i < st1->nnhpres; i++)
 957         {
 958             nc = i*st1->nhchainlength;
 959             for (j = 0; j < nc; j++)
 960             {
 961                 cmp_real(stdout, "nosehoover_xi",
 962                          i, st1->nhpres_xi[nc+j], st2->nhpres_xi[nc+j], ftol, abstol);
 963             }
 964         }
 965     }
 966
 967     cmp_int(stdout, "natoms", -1, st1->natoms, st2->natoms);
 968     if (st1->natoms == st2->natoms)
 969     {
 970         if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX)))
 971         {
 972             fprintf(stdout, "comparing x\n");
 973             cmp_rvecs(stdout, "x", st1->natoms, st1->x, st2->x, bRMSD, ftol, abstol);
 974         }
 975         if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV)))
 976         {
 977             fprintf(stdout, "comparing v\n");
 978             cmp_rvecs(stdout, "v", st1->natoms, st1->v, st2->v, bRMSD, ftol, abstol);
 979         }
 980     }
 981 }
 982
 983 void comp_tpx(const char *fn1, const char *fn2,
 984               gmx_bool bRMSD, real ftol, real abstol)
 985 {
 986     const char  *ff[2];
 987     t_tpxheader  sh[2];
 988     t_inputrec   ir[2];
 989     t_state      state[2];
 990     gmx_mtop_t   mtop[2];
 991     t_topology   top[2];
 992     int          i;
 993
 994     ff[0] = fn1;
 995     ff[1] = fn2;
 996     for (i = 0; i < (fn2 ? 2 : 1); i++)
 997     {
 998         read_tpx_state(ff[i], &(ir[i]), &state[i], NULL, &(mtop[i]));
 999     }
1000     if (fn2)
1001     {
1002         cmp_inputrec(stdout, &ir[0], &ir[1], ftol, abstol);
1003         /* Convert gmx_mtop_t to t_topology.
1004          * We should implement direct mtop comparison,
1005          * but it might be useful to keep t_topology comparison as an option.
1006          */
1007         top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
1008         top[1] = gmx_mtop_t_to_t_topology(&mtop[1]);
1009         cmp_top(stdout, &top[0], &top[1], ftol, abstol);
1010         cmp_groups(stdout, &mtop[0].groups, &mtop[1].groups,
1011                    mtop[0].natoms, mtop[1].natoms);
1012         comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
1013     }
1014     else
1015     {
1016         if (ir[0].efep == efepNO)
1017         {
1018             fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0].efep]);
1019         }
1020         else
1021         {
1022             if (ir[0].ePull != epullNO)
1023             {
1024                 comp_pull_AB(stdout, ir->pull, ftol, abstol);
1025             }
1026             /* Convert gmx_mtop_t to t_topology.
1027              * We should implement direct mtop comparison,
1028              * but it might be useful to keep t_topology comparison as an option.
1029              */
1030             top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
1031             cmp_top(stdout, &top[0], NULL, ftol, abstol);
1032         }
1033     }
1034 }
1035
1036 void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
1037                 gmx_bool bRMSD, real ftol, real abstol)
1038 {
1039     fprintf(fp, "\n");
1040     cmp_int(fp, "flags", -1, fr1->flags, fr2->flags);
1041     cmp_int(fp, "not_ok", -1, fr1->not_ok, fr2->not_ok);
1042     cmp_int(fp, "natoms", -1, fr1->natoms, fr2->natoms);
1043     cmp_real(fp, "t0", -1, fr1->t0, fr2->t0, ftol, abstol);
1044     if (cmp_bool(fp, "bTitle", -1, fr1->bTitle, fr2->bTitle))
1045     {
1046         cmp_str(fp, "title", -1, fr1->title, fr2->title);
1047     }
1048     if (cmp_bool(fp, "bStep", -1, fr1->bStep, fr2->bStep))
1049     {
1050         cmp_int(fp, "step", -1, fr1->step, fr2->step);
1051     }
1052     cmp_int(fp, "step", -1, fr1->step, fr2->step);
1053     if (cmp_bool(fp, "bTime", -1, fr1->bTime, fr2->bTime))
1054     {
1055         cmp_real(fp, "time", -1, fr1->time, fr2->time, ftol, abstol);
1056     }
1057     if (cmp_bool(fp, "bLambda", -1, fr1->bLambda, fr2->bLambda))
1058     {
1059         cmp_real(fp, "lambda", -1, fr1->lambda, fr2->lambda, ftol, abstol);
1060     }
1061     if (cmp_bool(fp, "bAtoms", -1, fr1->bAtoms, fr2->bAtoms))
1062     {
1063         cmp_atoms(fp, fr1->atoms, fr2->atoms, ftol, abstol);
1064     }
1065     if (cmp_bool(fp, "bPrec", -1, fr1->bPrec, fr2->bPrec))
1066     {
1067         cmp_real(fp, "prec", -1, fr1->prec, fr2->prec, ftol, abstol);
1068     }
1069     if (cmp_bool(fp, "bX", -1, fr1->bX, fr2->bX))
1070     {
1071         cmp_rvecs(fp, "x", min(fr1->natoms, fr2->natoms), fr1->x, fr2->x, bRMSD, ftol, abstol);
1072     }
1073     if (cmp_bool(fp, "bV", -1, fr1->bV, fr2->bV))
1074     {
1075         cmp_rvecs(fp, "v", min(fr1->natoms, fr2->natoms), fr1->v, fr2->v, bRMSD, ftol, abstol);
1076     }
1077     if (cmp_bool(fp, "bF", -1, fr1->bF, fr2->bF))
1078     {
1079         if (bRMSD)
1080         {
1081             cmp_rvecs(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, bRMSD, ftol, abstol);
1082         }
1083         else
1084         {
1085             cmp_rvecs_rmstol(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, ftol, abstol);
1086         }
1087     }
1088     if (cmp_bool(fp, "bBox", -1, fr1->bBox, fr2->bBox))
1089     {
1090         cmp_rvecs(fp, "box", 3, fr1->box, fr2->box, FALSE, ftol, abstol);
1091     }
1092 }
1093
1094 void comp_trx(const output_env_t oenv, const char *fn1, const char *fn2,
1095               gmx_bool bRMSD, real ftol, real abstol)
1096 {
1097     int          i;
1098     const char  *fn[2];
1099     t_trxframe   fr[2];
1100     t_trxstatus *status[2];
1101     gmx_bool     b[2];
1102
1103     fn[0] = fn1;
1104     fn[1] = fn2;
1105     fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
1106     for (i = 0; i < 2; i++)
1107     {
1108         b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X|TRX_READ_V|TRX_READ_F);
1109     }
1110
1111     if (b[0] && b[1])
1112     {
1113         do
1114         {
1115             comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
1116
1117             for (i = 0; i < 2; i++)
1118             {
1119                 b[i] = read_next_frame(oenv, status[i], &fr[i]);
1120             }
1121         }
1122         while (b[0] && b[1]);
1123
1124         for (i = 0; i < 2; i++)
1125         {
1126             if (b[i] && !b[1-i])
1127             {
1128                 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1-i], fn[i]);
1129             }
1130             close_trj(status[i]);
1131         }
1132     }
1133     if (!b[0] && !b[1])
1134     {
1135         fprintf(stdout, "\nBoth files read correctly\n");
1136     }
1137 }
1138
1139 static real ener_tensor_diag(int n, int *ind1, int *ind2,
1140                              gmx_enxnm_t *enm1,
1141                              int *tensi, int i,
1142                              t_energy e1[], t_energy e2[])
1143 {
1144     int  d1, d2;
1145     int  len;
1146     int  j;
1147     real prod1, prod2;
1148     int  nfound;
1149
1150     d1 = tensi[i]/DIM;
1151     d2 = tensi[i] - d1*DIM;
1152
1153     /* Find the diagonal elements d1 and d2 */
1154     len    = strlen(enm1[ind1[i]].name);
1155     prod1  = 1;
1156     prod2  = 1;
1157     nfound = 0;
1158     for (j = 0; j < n; j++)
1159     {
1160         if (tensi[j] >= 0 &&
1161             strlen(enm1[ind1[j]].name) == len &&
1162             strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len-2) == 0 &&
1163             (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2))
1164         {
1165             prod1 *= fabs(e1[ind1[j]].e);
1166             prod2 *= fabs(e2[ind2[j]].e);
1167             nfound++;
1168         }
1169     }
1170
1171     if (nfound == 2)
1172     {
1173         return 0.5*(sqrt(prod1) + sqrt(prod2));
1174     }
1175     else
1176     {
1177         return 0;
1178     }
1179 }
1180
1181 static gmx_bool enernm_equal(const char *nm1, const char *nm2)
1182 {
1183     int len1, len2;
1184
1185     len1 = strlen(nm1);
1186     len2 = strlen(nm2);
1187
1188     /* Remove " (bar)" at the end of a name */
1189     if (len1 > 6 && strcmp(nm1+len1-6, " (bar)") == 0)
1190     {
1191         len1 -= 6;
1192     }
1193     if (len2 > 6 && strcmp(nm2+len2-6, " (bar)") == 0)
1194     {
1195         len2 -= 6;
1196     }
1197
1198     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1199 }
1200
1201 static void cmp_energies(FILE *fp, int step1, int step2,
1202                          t_energy e1[], t_energy e2[],
1203                          gmx_enxnm_t *enm1, gmx_enxnm_t *enm2,
1204                          real ftol, real abstol,
1205                          int nre, int *ind1, int *ind2, int maxener)
1206 {
1207     int   i, ii;
1208     int  *tensi, len, d1, d2;
1209     real  ftol_i, abstol_i;
1210
1211     snew(tensi, maxener);
1212     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1213     for (i = 0; (i < maxener); i++)
1214     {
1215         ii       = ind1[i];
1216         tensi[i] = -1;
1217         len      = strlen(enm1[ii].name);
1218         if (len > 3 && enm1[ii].name[len-3] == '-')
1219         {
1220             d1 = enm1[ii].name[len-2] - 'X';
1221             d2 = enm1[ii].name[len-1] - 'X';
1222             if (d1 >= 0 && d1 < DIM &&
1223                 d2 >= 0 && d2 < DIM)
1224             {
1225                 tensi[i] = d1*DIM + d2;
1226             }
1227         }
1228     }
1229
1230     for (i = 0; (i < maxener); i++)
1231     {
1232         /* Check if this is an off-diagonal tensor element */
1233         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1234         {
1235             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1236             ftol_i = 5;
1237             /* Do the relative tolerance through an absolute tolerance times
1238              * the size of diagonal components of the tensor.
1239              */
1240             abstol_i = ftol*ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1241             if (debug)
1242             {
1243                 fprintf(debug, "tensor '%s' val %f diag %f\n",
1244                         enm1[i].name, e1[i].e, abstol_i/ftol);
1245             }
1246             if (abstol_i > 0)
1247             {
1248                 /* We found a diagonal, we need to check with the minimum tolerance */
1249                 abstol_i = min(abstol_i, abstol);
1250             }
1251             else
1252             {
1253                 /* We did not find a diagonal, ignore the relative tolerance check */
1254                 abstol_i = abstol;
1255             }
1256         }
1257         else
1258         {
1259             ftol_i   = ftol;
1260             abstol_i = abstol;
1261         }
1262         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1263         {
1264             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n",
1265                     enm1[ind1[i]].name,
1266                     step1, e1[ind1[i]].e,
1267                     step2, e2[ind2[i]].e);
1268         }
1269     }
1270
1271     sfree(tensi);
1272 }
1273
1274 #if 0
1275 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1276 {
1277     int  i;
1278     char bav[64], bt[64], bs[22];
1279
1280     cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1281     if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1282     {
1283         sprintf(bav, "step %s: disre rav", gmx_step_str(fr1->step, bs));
1284         sprintf(bt, "step %s: disre  rt", gmx_step_str(fr1->step, bs));
1285         for (i = 0; (i < fr1->ndisre); i++)
1286         {
1287             cmp_real(stdout, bav, i, fr1->disre_rm3tav[i], fr2->disre_rm3tav[i], ftol, abstol);
1288             cmp_real(stdout, bt, i, fr1->disre_rt[i], fr2->disre_rt[i], ftol, abstol);
1289         }
1290     }
1291 }
1292 #endif
1293
1294 static void cmp_eblocks(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1295 {
1296     int  i, j, k;
1297     char buf[64], bs[22];
1298
1299     cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1300     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1301     {
1302         for (j = 0; (j < fr1->nblock); j++)
1303         {
1304             t_enxblock *b1, *b2; /* convenience vars */
1305
1306             b1 = &(fr1->block[j]);
1307             b2 = &(fr2->block[j]);
1308
1309             sprintf(buf, "step %s: block[%d]", gmx_step_str(fr1->step, bs), j);
1310             cmp_int(stdout, buf, -1, b1->nsub, b2->nsub);
1311             cmp_int(stdout, buf, -1, b1->id, b2->id);
1312
1313             if ( (b1->nsub == b2->nsub) && (b1->id == b2->id) )
1314             {
1315                 for (i = 0; i < b1->nsub; i++)
1316                 {
1317                     t_enxsubblock *s1, *s2;
1318
1319                     s1 = &(b1->sub[i]);
1320                     s2 = &(b2->sub[i]);
1321
1322                     cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
1323                     cmp_gmx_large_int(stdout, buf, s1->nr, s2->nr);
1324
1325                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
1326                     {
1327                         switch (s1->type)
1328                         {
1329                             case xdr_datatype_float:
1330                                 for (k = 0; k < s1->nr; k++)
1331                                 {
1332                                     cmp_float(stdout, buf, i,
1333                                               s1->fval[k], s2->fval[k],
1334                                               ftol, abstol);
1335                                 }
1336                                 break;
1337                             case xdr_datatype_double:
1338                                 for (k = 0; k < s1->nr; k++)
1339                                 {
1340                                     cmp_double(stdout, buf, i,
1341                                                s1->dval[k], s2->dval[k],
1342                                                ftol, abstol);
1343                                 }
1344                                 break;
1345                             case xdr_datatype_int:
1346                                 for (k = 0; k < s1->nr; k++)
1347                                 {
1348                                     cmp_int(stdout, buf, i,
1349                                             s1->ival[k], s2->ival[k]);
1350                                 }
1351                                 break;
1352                             case xdr_datatype_large_int:
1353                                 for (k = 0; k < s1->nr; k++)
1354                                 {
1355                                     cmp_gmx_large_int(stdout, buf,
1356                                                       s1->lval[k], s2->lval[k]);
1357                                 }
1358                                 break;
1359                             case xdr_datatype_char:
1360                                 for (k = 0; k < s1->nr; k++)
1361                                 {
1362                                     cmp_uc(stdout, buf, i,
1363                                            s1->cval[k], s2->cval[k]);
1364                                 }
1365                                 break;
1366                             case xdr_datatype_string:
1367                                 for (k = 0; k < s1->nr; k++)
1368                                 {
1369                                     cmp_str(stdout, buf, i,
1370                                             s1->sval[k], s2->sval[k]);
1371                                 }
1372                                 break;
1373                             default:
1374                                 gmx_incons("Unknown data type!!");
1375                         }
1376                     }
1377                 }
1378             }
1379         }
1380     }
1381 }
1382
1383 void comp_enx(const char *fn1, const char *fn2, real ftol, real abstol, const char *lastener)
1384 {
1385     int            nre, nre1, nre2, block;
1386     ener_file_t    in1, in2;
1387     int            i, j, maxener, *ind1, *ind2, *have;
1388     char           buf[256];
1389     gmx_enxnm_t   *enm1 = NULL, *enm2 = NULL;
1390     t_enxframe    *fr1, *fr2;
1391     gmx_bool       b1, b2;
1392
1393     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1394
1395     in1 = open_enx(fn1, "r");
1396     in2 = open_enx(fn2, "r");
1397     do_enxnms(in1, &nre1, &enm1);
1398     do_enxnms(in2, &nre2, &enm2);
1399     if (nre1 != nre2)
1400     {
1401         fprintf(stdout, "There are %d and %d terms in the energy files\n\n",
1402                 nre1, nre2);
1403     }
1404     else
1405     {
1406         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1407     }
1408
1409     snew(ind1, nre1);
1410     snew(ind2, nre2);
1411     snew(have, nre2);
1412     nre = 0;
1413     for (i = 0; i < nre1; i++)
1414     {
1415         for (j = 0; j < nre2; j++)
1416         {
1417             if (enernm_equal(enm1[i].name, enm2[j].name))
1418             {
1419                 ind1[nre] = i;
1420                 ind2[nre] = j;
1421                 have[j]   = 1;
1422                 nre++;
1423                 break;
1424             }
1425         }
1426         if (nre == 0 || ind1[nre-1] != i)
1427         {
1428             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1429         }
1430     }
1431     for (i = 0; i < nre2; i++)
1432     {
1433         if (have[i] == 0)
1434         {
1435             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1436         }
1437     }
1438
1439     maxener = nre;
1440     for (i = 0; i < nre; i++)
1441     {
1442         if ((lastener != NULL) && (strstr(enm1[i].name, lastener) != NULL))
1443         {
1444             maxener = i+1;
1445             break;
1446         }
1447     }
1448
1449     fprintf(stdout, "There are %d terms to compare in the energy files\n\n",
1450             maxener);
1451
1452     for (i = 0; i < maxener; i++)
1453     {
1454         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1455     }
1456
1457     snew(fr1, 1);
1458     snew(fr2, 1);
1459     do
1460     {
1461         b1 = do_enx(in1, fr1);
1462         b2 = do_enx(in2, fr2);
1463         if (b1 && !b2)
1464         {
1465             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1466         }
1467         else if (!b1 && b2)
1468         {
1469             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1470         }
1471         else if (!b1 && !b2)
1472         {
1473             fprintf(stdout, "\nFiles read successfully\n");
1474         }
1475         else
1476         {
1477             cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1478             cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1479             /* We don't want to print the nre mismatch for every frame */
1480             /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1481             if ((fr1->nre >= nre) && (fr2->nre >= nre))
1482             {
1483                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener,
1484                              enm1, enm2, ftol, abstol, nre, ind1, ind2, maxener);
1485             }
1486             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1487             cmp_eblocks(fr1, fr2, ftol, abstol);
1488         }
1489     }
1490     while (b1 && b2);
1491
1492     close_enx(in1);
1493     close_enx(in2);
1494
1495     free_enxframe(fr2);
1496     sfree(fr2);
1497     free_enxframe(fr1);
1498     sfree(fr1);
1499 }