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