src/kernel/tpbcmp.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <stdio.h>
  42 #include <string.h>
  43 #include "main.h"
  44 #include "macros.h"
  45 #include "smalloc.h"
  46 #include "futil.h"
  47 #include "statutil.h"
  48 #include "sysstuff.h"
  49 #include "txtdump.h"
  50 #include "gmx_fatal.h"
  51 #include "names.h"
  52 #include "tpxio.h"
  53 #include "enxio.h"
  54 #include "mtop_util.h"
  55 #include "string2.h"
  56
  57 static void cmp_int(FILE *fp,const char *s,int index,int i1,int i2)
  58 {
  59   if (i1 != i2) {
  60     if (index != -1)
  61       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  62     else
  63       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  64   }
  65 }
  66
  67 static void cmp_gmx_large_int(FILE *fp,const char *s,gmx_large_int_t i1,gmx_large_int_t i2)
  68 {
  69   if (i1 != i2) {
  70     fprintf(fp,"%s (",s);
  71     fprintf(fp,gmx_large_int_pfmt,i1);
  72     fprintf(fp," - ");
  73     fprintf(fp,gmx_large_int_pfmt,i2);
  74     fprintf(fp,")\n");
  75   }
  76 }
  77
  78 static void cmp_us(FILE *fp,const char *s,int index,unsigned short i1,unsigned short i2)
  79 {
  80   if (i1 != i2) {
  81     if (index != -1)
  82       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  83     else
  84       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  85   }
  86 }
  87
  88 static void cmp_uc(FILE *fp,const char *s,int index,unsigned char i1,unsigned char i2)
  89 {
  90   if (i1 != i2) {
  91     if (index != -1)
  92       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  93     else
  94       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  95   }
  96 }
  97
  98 static gmx_bool cmp_bool(FILE *fp, const char *s, int index, gmx_bool b1, gmx_bool b2)
  99 {
 100   if (b1) {
 101     b1 = 1;
 102   } else {
 103     b1 = 0;
 104   }
 105   if (b2) {
 106     b2 = 1;
 107   } else {
 108     b2 = 0;
 109   }
 110   if (b1 != b2) {
 111     if (index != -1)
 112       fprintf(fp,"%s[%d] (%s - %s)\n",s,index,
 113               bool_names[b1],bool_names[b2]);
 114     else
 115       fprintf(fp,"%s (%s - %s)\n",s,
 116               bool_names[b1],bool_names[b2]);
 117   }
 118   return b1 && b2;
 119 }
 120
 121 static void cmp_str(FILE *fp, const char *s, int index,
 122                     const char *s1, const char *s2)
 123 {
 124   if (strcmp(s1,s2) != 0) {
 125     if (index != -1)
 126       fprintf(fp,"%s[%d] (%s - %s)\n",s,index,s1,s2);
 127     else
 128       fprintf(fp,"%s (%s - %s)\n",s,s1,s2);
 129   }
 130 }
 131
 132 static gmx_bool equal_real(real i1,real i2,real ftol,real abstol)
 133 {
 134     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
 135 }
 136
 137 static gmx_bool equal_float(float i1,float i2,float ftol,float abstol)
 138 {
 139     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
 140 }
 141
 142 static gmx_bool equal_double(double i1,double i2,real ftol,real abstol)
 143 {
 144     return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
 145 }
 146
 147 static void
 148 cmp_real(FILE *fp,const char *s,int index,real i1,real i2,real ftol,real abstol)
 149 {
 150   if (!equal_real(i1,i2,ftol,abstol)) {
 151     if (index != -1)
 152       fprintf(fp,"%s[%2d] (%e - %e)\n",s,index,i1,i2);
 153     else
 154       fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
 155   }
 156 }
 157
 158 static void
 159 cmp_float(FILE *fp,const char *s,int index,float i1,float i2,float ftol,float abstol)
 160 {
 161   if (!equal_float(i1,i2,ftol,abstol)) {
 162     if (index != -1)
 163       fprintf(fp,"%s[%2d] (%e - %e)\n",s,index,i1,i2);
 164     else
 165       fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
 166   }
 167 }
 168
 169
 170
 171 static void
 172 cmp_double(FILE *fp,const char *s,int index,double i1,double i2,double ftol,double abstol)
 173 {
 174   if (!equal_double(i1,i2,ftol,abstol)) {
 175     if (index != -1)
 176       fprintf(fp,"%s[%2d] (%16.9e - %16.9e)\n",s,index,i1,i2);
 177     else
 178       fprintf(fp,"%s (%16.9e - %16.9e)\n",s,i1,i2);
 179   }
 180 }
 181
 182 static void cmp_rvec(FILE *fp,const char *s,int index,rvec i1,rvec i2,real ftol,real abstol)
 183 {
 184     if(!equal_real(i1[XX],i2[XX],ftol,abstol) ||
 185        !equal_real(i1[YY],i2[YY],ftol,abstol) ||
 186        !equal_real(i1[ZZ],i2[ZZ],ftol,abstol))
 187     {
 188         if (index != -1)
 189             fprintf(fp,"%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
 190                     s,index,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 191         else
 192             fprintf(fp,"%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
 193                     s,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 194     }
 195 }
 196
 197 static void cmp_ivec(FILE *fp,const char *s,int index,ivec i1,ivec i2)
 198 {
 199   if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ])) {
 200     if (index != -1)
 201       fprintf(fp,"%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,index,
 202               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 203     else
 204       fprintf(fp,"%s (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,
 205               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 206   }
 207 }
 208
 209 static void cmp_ilist(FILE *fp,int ftype,t_ilist *il1,t_ilist *il2)
 210 {
 211   int i;
 212   char buf[256];
 213
 214   fprintf(fp,"comparing ilist %s\n",interaction_function[ftype].name);
 215   sprintf(buf,"%s->nr",interaction_function[ftype].name);
 216   cmp_int(fp,buf,-1,il1->nr,il2->nr);
 217   sprintf(buf,"%s->iatoms",interaction_function[ftype].name);
 218   if (((il1->nr > 0) && (!il1->iatoms)) ||
 219       ((il2->nr > 0) && (!il2->iatoms)) ||
 220       ((il1->nr != il2->nr)))
 221     fprintf(fp,"Comparing radically different topologies - %s is different\n",
 222             buf);
 223   else
 224     for(i=0; (i<il1->nr); i++)
 225       cmp_int(fp,buf,i,il1->iatoms[i],il2->iatoms[i]);
 226 }
 227
 228 void cmp_iparm(FILE *fp,const char *s,t_functype ft,
 229                t_iparams ip1,t_iparams ip2,real ftol,real abstol)
 230 {
 231   int i;
 232   gmx_bool bDiff;
 233
 234   bDiff=FALSE;
 235   for(i=0; i<MAXFORCEPARAM && !bDiff; i++)
 236     bDiff = !equal_real(ip1.generic.buf[i],ip2.generic.buf[i],ftol,abstol);
 237   if (bDiff) {
 238     fprintf(fp,"%s1: ",s);
 239     pr_iparams(fp,ft,&ip1);
 240     fprintf(fp,"%s2: ",s);
 241     pr_iparams(fp,ft,&ip2);
 242   }
 243 }
 244
 245 void cmp_iparm_AB(FILE *fp,const char *s,t_functype ft,t_iparams ip1,real ftol,real abstol)
 246 {
 247   int nrfpA,nrfpB,p0,i;
 248   gmx_bool bDiff;
 249
 250   /* Normally the first parameter is perturbable */
 251   p0 = 0;
 252   nrfpA = interaction_function[ft].nrfpA;
 253   nrfpB = interaction_function[ft].nrfpB;
 254   if (ft == F_PDIHS) {
 255     nrfpB = 2;
 256   } else if (interaction_function[ft].flags & IF_TABULATED) {
 257     /* For tabulated interactions only the second parameter is perturbable */
 258     p0 = 1;
 259     nrfpB = 1;
 260   }
 261   bDiff=FALSE;
 262   for(i=0; i<nrfpB && !bDiff; i++) {
 263     bDiff = !equal_real(ip1.generic.buf[p0+i],ip1.generic.buf[nrfpA+i],ftol,abstol);
 264   }
 265   if (bDiff) {
 266     fprintf(fp,"%s: ",s);
 267     pr_iparams(fp,ft,&ip1);
 268   }
 269 }
 270
 271 static void cmp_idef(FILE *fp,t_idef *id1,t_idef *id2,real ftol,real abstol)
 272 {
 273   int i;
 274   char buf1[64],buf2[64];
 275
 276   fprintf(fp,"comparing idef\n");
 277   if (id2) {
 278     cmp_int(fp,"idef->ntypes",-1,id1->ntypes,id2->ntypes);
 279     cmp_int(fp,"idef->atnr",  -1,id1->atnr,id2->atnr);
 280     for(i=0; (i<id1->ntypes); i++) {
 281       sprintf(buf1,"idef->functype[%d]",i);
 282       sprintf(buf2,"idef->iparam[%d]",i);
 283       cmp_int(fp,buf1,i,(int)id1->functype[i],(int)id2->functype[i]);
 284       cmp_iparm(fp,buf2,id1->functype[i],
 285                 id1->iparams[i],id2->iparams[i],ftol,abstol);
 286     }
 287     cmp_real(fp,"fudgeQQ",-1,id1->fudgeQQ,id2->fudgeQQ,ftol,abstol);
 288     for(i=0; (i<F_NRE); i++)
 289       cmp_ilist(fp,i,&(id1->il[i]),&(id2->il[i]));
 290   } else {
 291     for(i=0; (i<id1->ntypes); i++)
 292       cmp_iparm_AB(fp,"idef->iparam",id1->functype[i],id1->iparams[i],ftol,abstol);
 293   }
 294 }
 295
 296 static void cmp_block(FILE *fp,t_block *b1,t_block *b2,const char *s)
 297 {
 298   int i,j,k;
 299   char buf[32];
 300
 301   fprintf(fp,"comparing block %s\n",s);
 302   sprintf(buf,"%s.nr",s);
 303   cmp_int(fp,buf,-1,b1->nr,b2->nr);
 304 }
 305
 306 static void cmp_blocka(FILE *fp,t_blocka *b1,t_blocka *b2,const char *s)
 307 {
 308   int i,j,k;
 309   char buf[32];
 310
 311   fprintf(fp,"comparing blocka %s\n",s);
 312   sprintf(buf,"%s.nr",s);
 313   cmp_int(fp,buf,-1,b1->nr,b2->nr);
 314   sprintf(buf,"%s.nra",s);
 315   cmp_int(fp,buf,-1,b1->nra,b2->nra);
 316 }
 317
 318 static void cmp_atom(FILE *fp,int index,t_atom *a1,t_atom *a2,real ftol,real abstol)
 319 {
 320   int  i;
 321   char buf[256];
 322
 323   if (a2) {
 324     cmp_us(fp,"atom.type",index,a1->type,a2->type);
 325     cmp_us(fp,"atom.ptype",index,a1->ptype,a2->ptype);
 326     cmp_int(fp,"atom.resind",index,a1->resind,a2->resind);
 327     cmp_int(fp,"atom.atomnumber",index,a1->atomnumber,a2->atomnumber);
 328     cmp_real(fp,"atom.m",index,a1->m,a2->m,ftol,abstol);
 329     cmp_real(fp,"atom.q",index,a1->q,a2->q,ftol,abstol);
 330     cmp_us(fp,"atom.typeB",index,a1->typeB,a2->typeB);
 331     cmp_real(fp,"atom.mB",index,a1->mB,a2->mB,ftol,abstol);
 332     cmp_real(fp,"atom.qB",index,a1->qB,a2->qB,ftol,abstol);
 333   } else {
 334     cmp_us(fp,"atom.type",index,a1->type,a1->typeB);
 335     cmp_real(fp,"atom.m",index,a1->m,a1->mB,ftol,abstol);
 336     cmp_real(fp,"atom.q",index,a1->q,a1->qB,ftol,abstol);
 337   }
 338 }
 339
 340 static void cmp_atoms(FILE *fp,t_atoms *a1,t_atoms *a2,real ftol, real abstol)
 341 {
 342   int i;
 343
 344   fprintf(fp,"comparing atoms\n");
 345
 346   if (a2) {
 347     cmp_int(fp,"atoms->nr",-1,a1->nr,a2->nr);
 348     for(i=0; (i<a1->nr); i++)
 349       cmp_atom(fp,i,&(a1->atom[i]),&(a2->atom[i]),ftol,abstol);
 350   } else {
 351     for(i=0; (i<a1->nr); i++)
 352       cmp_atom(fp,i,&(a1->atom[i]),NULL,ftol,abstol);
 353   }
 354 }
 355
 356 static void cmp_top(FILE *fp,t_topology *t1,t_topology *t2,real ftol, real abstol)
 357 {
 358   int i;
 359
 360   fprintf(fp,"comparing top\n");
 361   if (t2) {
 362     cmp_idef(fp,&(t1->idef),&(t2->idef),ftol,abstol);
 363     cmp_atoms(fp,&(t1->atoms),&(t2->atoms),ftol,abstol);
 364     cmp_block(fp,&t1->cgs,&t2->cgs,"cgs");
 365     cmp_block(fp,&t1->mols,&t2->mols,"mols");
 366     cmp_blocka(fp,&t1->excls,&t2->excls,"excls");
 367   } else {
 368     cmp_idef(fp,&(t1->idef),NULL,ftol,abstol);
 369     cmp_atoms(fp,&(t1->atoms),NULL,ftol,abstol);
 370   }
 371 }
 372
 373 static void cmp_groups(FILE *fp,gmx_groups_t *g0,gmx_groups_t *g1,
 374                        int natoms0,int natoms1)
 375 {
 376   int  i,j,ndiff;
 377   char buf[32];
 378
 379   fprintf(fp,"comparing groups\n");
 380
 381   for(i=0; i<egcNR; i++) {
 382     sprintf(buf,"grps[%d].nr",i);
 383     cmp_int(fp,buf,-1,g0->grps[i].nr,g1->grps[i].nr);
 384     if (g0->grps[i].nr == g1->grps[i].nr) {
 385       for(j=0; j<g0->grps[i].nr; j++) {
 386           sprintf(buf,"grps[%d].name[%d]",i,j);
 387           cmp_str(fp,buf,-1,
 388                   *g0->grpname[g0->grps[i].nm_ind[j]],
 389                   *g1->grpname[g1->grps[i].nm_ind[j]]);
 390       }
 391     }
 392     cmp_int(fp,"ngrpnr",i,g0->ngrpnr[i],g1->ngrpnr[i]);
 393     if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
 394         (g0->grpnr[i] != NULL || g1->grpnr[i] != NULL)) {
 395       for(j=0; j<natoms0; j++) {
 396         cmp_int(fp,gtypes[i],j,ggrpnr(g0,i,j),ggrpnr(g1,i,j));
 397       }
 398     }
 399   }
 400   /* We have compared the names in the groups lists,
 401    * so we can skip the grpname list comparison.
 402    */
 403 }
 404
 405 static void cmp_rvecs(FILE *fp,const char *title,int n,rvec x1[],rvec x2[],
 406                       gmx_bool bRMSD,real ftol,real abstol)
 407 {
 408   int i,m;
 409   double d,ssd;
 410
 411   if (bRMSD) {
 412     ssd = 0;
 413     for(i=0; (i<n); i++) {
 414       for(m=0; m<DIM; m++) {
 415         d = x1[i][m] - x2[i][m];
 416         ssd += d*d;
 417       }
 418     }
 419     fprintf(fp,"%s RMSD %g\n",title,sqrt(ssd/n));
 420   } else {
 421     for(i=0; (i<n); i++) {
 422       cmp_rvec(fp,title,i,x1[i],x2[i],ftol,abstol);
 423     }
 424   }
 425 }
 426
 427 static void cmp_grpopts(FILE *fp,t_grpopts *opt1,t_grpopts *opt2,real ftol, real abstol)
 428 {
 429   int i,j;
 430   char buf1[256],buf2[256];
 431
 432   cmp_int(fp,"inputrec->grpopts.ngtc",-1,  opt1->ngtc,opt2->ngtc);
 433   cmp_int(fp,"inputrec->grpopts.ngacc",-1, opt1->ngacc,opt2->ngacc);
 434   cmp_int(fp,"inputrec->grpopts.ngfrz",-1, opt1->ngfrz,opt2->ngfrz);
 435   cmp_int(fp,"inputrec->grpopts.ngener",-1,opt1->ngener,opt2->ngener);
 436   for(i=0; (i<min(opt1->ngtc,opt2->ngtc)); i++) {
 437     cmp_real(fp,"inputrec->grpopts.nrdf",i,opt1->nrdf[i],opt2->nrdf[i],ftol,abstol);
 438     cmp_real(fp,"inputrec->grpopts.ref_t",i,opt1->ref_t[i],opt2->ref_t[i],ftol,abstol);
 439     cmp_real(fp,"inputrec->grpopts.tau_t",i,opt1->tau_t[i],opt2->tau_t[i],ftol,abstol);
 440     cmp_int(fp,"inputrec->grpopts.annealing",i,opt1->annealing[i],opt2->annealing[i]);
 441     cmp_int(fp,"inputrec->grpopts.anneal_npoints",i,
 442             opt1->anneal_npoints[i],opt2->anneal_npoints[i]);
 443     if(opt1->anneal_npoints[i]==opt2->anneal_npoints[i]) {
 444       sprintf(buf1,"inputrec->grpopts.anneal_time[%d]",i);
 445       sprintf(buf2,"inputrec->grpopts.anneal_temp[%d]",i);
 446       for(j=0;j<opt1->anneal_npoints[i];j++) {
 447         cmp_real(fp,buf1,j,opt1->anneal_time[i][j],opt2->anneal_time[i][j],ftol,abstol);
 448         cmp_real(fp,buf2,j,opt1->anneal_temp[i][j],opt2->anneal_temp[i][j],ftol,abstol);
 449       }
 450     }
 451   }
 452   if (opt1->ngener == opt2->ngener) {
 453     for(i=0; i<opt1->ngener; i++)
 454       for(j=i; j<opt1->ngener; j++) {
 455         sprintf(buf1,"inputrec->grpopts.egp_flags[%d]",i);
 456         cmp_int(fp,buf1,j,
 457                 opt1->egp_flags[opt1->ngener*i+j],
 458                 opt2->egp_flags[opt1->ngener*i+j]);
 459       }
 460   }
 461   for(i=0; (i<min(opt1->ngacc,opt2->ngacc)); i++)
 462     cmp_rvec(fp,"inputrec->grpopts.acc",i,opt1->acc[i],opt2->acc[i],ftol,abstol);
 463   for(i=0; (i<min(opt1->ngfrz,opt2->ngfrz)); i++)
 464     cmp_ivec(fp,"inputrec->grpopts.nFreeze",i,opt1->nFreeze[i],opt2->nFreeze[i]);
 465 }
 466
 467 static void cmp_cosines(FILE *fp,const char *s,t_cosines c1[DIM],t_cosines c2[DIM],real ftol, real abstol)
 468 {
 469   int i,m;
 470   char buf[256];
 471
 472   for(m=0; (m<DIM); m++) {
 473     sprintf(buf,"inputrec->%s[%d]",s,m);
 474     cmp_int(fp,buf,0,c1->n,c2->n);
 475     for(i=0; (i<min(c1->n,c2->n)); i++) {
 476       cmp_real(fp,buf,i,c1->a[i],c2->a[i],ftol,abstol);
 477       cmp_real(fp,buf,i,c1->phi[i],c2->phi[i],ftol,abstol);
 478     }
 479   }
 480 }
 481
 482 static void cmp_pull(FILE *fp,t_pull *pull1,t_pull *pull2,real ftol, real abstol)
 483 {
 484   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");
 485 }
 486
 487 static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2,real ftol, real abstol)
 488 {
 489   int i;
 490
 491   fprintf(fp,"comparing inputrec\n");
 492
 493   /* gcc 2.96 doesnt like these defines at all, but issues a huge list
 494    * of warnings. Maybe it will change in future versions, but for the
 495    * moment I've spelled them out instead. /EL 000820
 496    * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 497    * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 498    * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
 499    */
 500   cmp_int(fp,"inputrec->eI",-1,ir1->eI,ir2->eI);
 501   cmp_gmx_large_int(fp,"inputrec->nsteps",ir1->nsteps,ir2->nsteps);
 502   cmp_gmx_large_int(fp,"inputrec->init_step",ir1->init_step,ir2->init_step);
 503   cmp_int(fp,"inputrec->simulation_part",-1,ir1->simulation_part,ir2->simulation_part);
 504   cmp_int(fp,"inputrec->ePBC",-1,ir1->ePBC,ir2->ePBC);
 505   cmp_int(fp,"inputrec->bPeriodicMols",-1,ir1->bPeriodicMols,ir2->bPeriodicMols);
 506   cmp_int(fp,"inputrec->ns_type",-1,ir1->ns_type,ir2->ns_type);
 507   cmp_int(fp,"inputrec->nstlist",-1,ir1->nstlist,ir2->nstlist);
 508   cmp_int(fp,"inputrec->ndelta",-1,ir1->ndelta,ir2->ndelta);
 509   cmp_int(fp,"inputrec->nstcomm",-1,ir1->nstcomm,ir2->nstcomm);
 510   cmp_int(fp,"inputrec->comm_mode",-1,ir1->comm_mode,ir2->comm_mode);
 511   cmp_int(fp,"inputrec->nstcheckpoint",-1,ir1->nstcheckpoint,ir2->nstcheckpoint);
 512   cmp_int(fp,"inputrec->nstlog",-1,ir1->nstlog,ir2->nstlog);
 513   cmp_int(fp,"inputrec->nstxout",-1,ir1->nstxout,ir2->nstxout);
 514   cmp_int(fp,"inputrec->nstvout",-1,ir1->nstvout,ir2->nstvout);
 515   cmp_int(fp,"inputrec->nstfout",-1,ir1->nstfout,ir2->nstfout);
 516   cmp_int(fp,"inputrec->nstcalcenergy",-1,ir1->nstcalcenergy,ir2->nstcalcenergy);
 517   cmp_int(fp,"inputrec->nstenergy",-1,ir1->nstenergy,ir2->nstenergy);
 518   cmp_int(fp,"inputrec->nstxtcout",-1,ir1->nstxtcout,ir2->nstxtcout);
 519   cmp_double(fp,"inputrec->init_t",-1,ir1->init_t,ir2->init_t,ftol,abstol);
 520   cmp_double(fp,"inputrec->delta_t",-1,ir1->delta_t,ir2->delta_t,ftol,abstol);
 521   cmp_real(fp,"inputrec->xtcprec",-1,ir1->xtcprec,ir2->xtcprec,ftol,abstol);
 522   cmp_int(fp,"inputrec->nkx",-1,ir1->nkx,ir2->nkx);
 523   cmp_int(fp,"inputrec->nky",-1,ir1->nky,ir2->nky);
 524   cmp_int(fp,"inputrec->nkz",-1,ir1->nkz,ir2->nkz);
 525   cmp_int(fp,"inputrec->pme_order",-1,ir1->pme_order,ir2->pme_order);
 526   cmp_real(fp,"inputrec->ewald_rtol",-1,ir1->ewald_rtol,ir2->ewald_rtol,ftol,abstol);
 527   cmp_int(fp,"inputrec->ewald_geometry",-1,ir1->ewald_geometry,ir2->ewald_geometry);
 528   cmp_real(fp,"inputrec->epsilon_surface",-1,ir1->epsilon_surface,ir2->epsilon_surface,ftol,abstol);
 529   cmp_int(fp,"inputrec->bOptFFT",-1,ir1->bOptFFT,ir2->bOptFFT);
 530   cmp_int(fp,"inputrec->bContinuation",-1,ir1->bContinuation,ir2->bContinuation);
 531   cmp_int(fp,"inputrec->bShakeSOR",-1,ir1->bShakeSOR,ir2->bShakeSOR);
 532   cmp_int(fp,"inputrec->etc",-1,ir1->etc,ir2->etc);
 533   cmp_int(fp,"inputrec->epc",-1,ir1->epc,ir2->epc);
 534   cmp_int(fp,"inputrec->epct",-1,ir1->epct,ir2->epct);
 535   cmp_real(fp,"inputrec->tau_p",-1,ir1->tau_p,ir2->tau_p,ftol,abstol);
 536   cmp_rvec(fp,"inputrec->ref_p(x)",-1,ir1->ref_p[XX],ir2->ref_p[XX],ftol,abstol);
 537   cmp_rvec(fp,"inputrec->ref_p(y)",-1,ir1->ref_p[YY],ir2->ref_p[YY],ftol,abstol);
 538   cmp_rvec(fp,"inputrec->ref_p(z)",-1,ir1->ref_p[ZZ],ir2->ref_p[ZZ],ftol,abstol);
 539   cmp_rvec(fp,"inputrec->compress(x)",-1,ir1->compress[XX],ir2->compress[XX],ftol,abstol);
 540   cmp_rvec(fp,"inputrec->compress(y)",-1,ir1->compress[YY],ir2->compress[YY],ftol,abstol);
 541   cmp_rvec(fp,"inputrec->compress(z)",-1,ir1->compress[ZZ],ir2->compress[ZZ],ftol,abstol);
 542   cmp_int(fp,"refcoord_scaling",-1,ir1->refcoord_scaling,ir2->refcoord_scaling);
 543    cmp_rvec(fp,"inputrec->posres_com",-1,ir1->posres_com,ir2->posres_com,ftol,abstol);
 544    cmp_rvec(fp,"inputrec->posres_comB",-1,ir1->posres_comB,ir2->posres_comB,ftol,abstol);
 545    cmp_int(fp,"inputrec->andersen_seed",-1,ir1->andersen_seed,ir2->andersen_seed);
 546   cmp_real(fp,"inputrec->rlist",-1,ir1->rlist,ir2->rlist,ftol,abstol);
 547   cmp_real(fp,"inputrec->rlistlong",-1,ir1->rlistlong,ir2->rlistlong,ftol,abstol);
 548   cmp_real(fp,"inputrec->rtpi",-1,ir1->rtpi,ir2->rtpi,ftol,abstol);
 549   cmp_int(fp,"inputrec->coulombtype",-1,ir1->coulombtype,ir2->coulombtype);
 550   cmp_real(fp,"inputrec->rcoulomb_switch",-1,ir1->rcoulomb_switch,ir2->rcoulomb_switch,ftol,abstol);
 551   cmp_real(fp,"inputrec->rcoulomb",-1,ir1->rcoulomb,ir2->rcoulomb,ftol,abstol);
 552   cmp_int(fp,"inputrec->vdwtype",-1,ir1->vdwtype,ir2->vdwtype);
 553   cmp_real(fp,"inputrec->rvdw_switch",-1,ir1->rvdw_switch,ir2->rvdw_switch,ftol,abstol);
 554   cmp_real(fp,"inputrec->rvdw",-1,ir1->rvdw,ir2->rvdw,ftol,abstol);
 555   cmp_real(fp,"inputrec->epsilon_r",-1,ir1->epsilon_r,ir2->epsilon_r,ftol,abstol);
 556   cmp_real(fp,"inputrec->epsilon_rf",-1,ir1->epsilon_rf,ir2->epsilon_rf,ftol,abstol);
 557   cmp_real(fp,"inputrec->tabext",-1,ir1->tabext,ir2->tabext,ftol,abstol);
 558   cmp_int(fp,"inputrec->implicit_solvent",-1,ir1->implicit_solvent,ir2->implicit_solvent);
 559   cmp_int(fp,"inputrec->gb_algorithm",-1,ir1->gb_algorithm,ir2->gb_algorithm);
 560   cmp_int(fp,"inputrec->nstgbradii",-1,ir1->nstgbradii,ir2->nstgbradii);
 561   cmp_real(fp,"inputrec->rgbradii",-1,ir1->rgbradii,ir2->rgbradii,ftol,abstol);
 562   cmp_real(fp,"inputrec->gb_saltconc",-1,ir1->gb_saltconc,ir2->gb_saltconc,ftol,abstol);
 563   cmp_real(fp,"inputrec->gb_epsilon_solvent",-1,ir1->gb_epsilon_solvent,ir2->gb_epsilon_solvent,ftol,abstol);
 564   cmp_real(fp,"inputrec->gb_obc_alpha",-1,ir1->gb_obc_alpha,ir2->gb_obc_alpha,ftol,abstol);
 565   cmp_real(fp,"inputrec->gb_obc_beta",-1,ir1->gb_obc_beta,ir2->gb_obc_beta,ftol,abstol);
 566   cmp_real(fp,"inputrec->gb_obc_gamma",-1,ir1->gb_obc_gamma,ir2->gb_obc_gamma,ftol,abstol);
 567   cmp_real(fp,"inputrec->gb_dielectric_offset",-1,ir1->gb_dielectric_offset,ir2->gb_dielectric_offset,ftol,abstol);
 568   cmp_int(fp,"inputrec->sa_algorithm",-1,ir1->sa_algorithm,ir2->sa_algorithm);
 569   cmp_real(fp,"inputrec->sa_surface_tension",-1,ir1->sa_surface_tension,ir2->sa_surface_tension,ftol,abstol);
 570
 571   cmp_int(fp,"inputrec->eDispCorr",-1,ir1->eDispCorr,ir2->eDispCorr);
 572   cmp_real(fp,"inputrec->shake_tol",-1,ir1->shake_tol,ir2->shake_tol,ftol,abstol);
 573   cmp_int(fp,"inputrec->efep",-1,ir1->efep,ir2->efep);
 574   cmp_double(fp,"inputrec->init_lambda",-1,ir1->init_lambda,ir2->init_lambda,ftol,abstol);
 575   cmp_double(fp,"inputrec->delta_lambda",-1,ir1->delta_lambda,ir2->delta_lambda,ftol,abstol);
 576   cmp_int(fp,"inputrec->n_foreign_lambda",-1,ir1->n_flambda,ir2->n_flambda);
 577   for(i=0; i<min(ir1->n_flambda,ir2->n_flambda); i++) {
 578     cmp_double(fp,"inputrec->foreign_lambda",-1,ir1->flambda[i],ir2->flambda[i],ftol,abstol);
 579   }
 580   cmp_real(fp,"inputrec->sc_alpha",-1,ir1->sc_alpha,ir2->sc_alpha,ftol,abstol);
 581   cmp_int(fp,"inputrec->sc_power",-1,ir1->sc_power,ir2->sc_power);
 582   cmp_real(fp,"inputrec->sc_sigma",-1,ir1->sc_sigma,ir2->sc_sigma,ftol,abstol);
 583   cmp_real(fp,"inputrec->sc_sigma_min",-1,ir1->sc_sigma_min,ir2->sc_sigma_min,ftol,abstol);
 584   cmp_int(fp,"inputrec->nstdhdl",-1,ir1->nstdhdl,ir2->nstdhdl);
 585   cmp_int(fp,"inputrec->separate_dhdl_file",-1,ir1->separate_dhdl_file,ir2->nstdhdl);
 586   cmp_int(fp,"inputrec->dhdl_derivatives",-1,ir1->dhdl_derivatives,ir2->dhdl_derivatives);
 587   cmp_int(fp,"inputrec->dh_hist_size",-1,ir1->dh_hist_size,ir2->dh_hist_size);
 588   cmp_double(fp,"inputrec->dh_hist_spacing",-1,ir1->dh_hist_spacing,ir2->dh_hist_spacing,ftol,abstol);
 589
 590   cmp_int(fp,"inputrec->nwall",-1,ir1->nwall,ir2->nwall);
 591   cmp_int(fp,"inputrec->wall_type",-1,ir1->wall_type,ir2->wall_type);
 592   cmp_int(fp,"inputrec->wall_atomtype[0]",-1,ir1->wall_atomtype[0],ir2->wall_atomtype[0]);
 593   cmp_int(fp,"inputrec->wall_atomtype[1]",-1,ir1->wall_atomtype[1],ir2->wall_atomtype[1]);
 594   cmp_real(fp,"inputrec->wall_density[0]",-1,ir1->wall_density[0],ir2->wall_density[0],ftol,abstol);
 595   cmp_real(fp,"inputrec->wall_density[1]",-1,ir1->wall_density[1],ir2->wall_density[1],ftol,abstol);
 596   cmp_real(fp,"inputrec->wall_ewald_zfac",-1,ir1->wall_ewald_zfac,ir2->wall_ewald_zfac,ftol,abstol);
 597
 598   cmp_int(fp,"inputrec->ePull",-1,ir1->ePull,ir2->ePull);
 599   if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
 600     cmp_pull(fp,ir1->pull,ir2->pull,ftol,abstol);
 601
 602   cmp_int(fp,"inputrec->eDisre",-1,ir1->eDisre,ir2->eDisre);
 603   cmp_real(fp,"inputrec->dr_fc",-1,ir1->dr_fc,ir2->dr_fc,ftol,abstol);
 604   cmp_int(fp,"inputrec->eDisreWeighting",-1,ir1->eDisreWeighting,ir2->eDisreWeighting);
 605   cmp_int(fp,"inputrec->bDisreMixed",-1,ir1->bDisreMixed,ir2->bDisreMixed);
 606   cmp_int(fp,"inputrec->nstdisreout",-1,ir1->nstdisreout,ir2->nstdisreout);
 607   cmp_real(fp,"inputrec->dr_tau",-1,ir1->dr_tau,ir2->dr_tau,ftol,abstol);
 608   cmp_real(fp,"inputrec->orires_fc",-1,ir1->orires_fc,ir2->orires_fc,ftol,abstol);
 609   cmp_real(fp,"inputrec->orires_tau",-1,ir1->orires_tau,ir2->orires_tau,ftol,abstol);
 610   cmp_int(fp,"inputrec->nstorireout",-1,ir1->nstorireout,ir2->nstorireout);
 611   cmp_real(fp,"inputrec->dihre_fc",-1,ir1->dihre_fc,ir2->dihre_fc,ftol,abstol);
 612   cmp_real(fp,"inputrec->em_stepsize",-1,ir1->em_stepsize,ir2->em_stepsize,ftol,abstol);
 613   cmp_real(fp,"inputrec->em_tol",-1,ir1->em_tol,ir2->em_tol,ftol,abstol);
 614   cmp_int(fp,"inputrec->niter",-1,ir1->niter,ir2->niter);
 615   cmp_real(fp,"inputrec->fc_stepsize",-1,ir1->fc_stepsize,ir2->fc_stepsize,ftol,abstol);
 616   cmp_int(fp,"inputrec->nstcgsteep",-1,ir1->nstcgsteep,ir2->nstcgsteep);
 617   cmp_int(fp,"inputrec->nbfgscorr",0,ir1->nbfgscorr,ir2->nbfgscorr);
 618   cmp_int(fp,"inputrec->eConstrAlg",-1,ir1->eConstrAlg,ir2->eConstrAlg);
 619   cmp_int(fp,"inputrec->nProjOrder",-1,ir1->nProjOrder,ir2->nProjOrder);
 620   cmp_real(fp,"inputrec->LincsWarnAngle",-1,ir1->LincsWarnAngle,ir2->LincsWarnAngle,ftol,abstol);
 621   cmp_int(fp,"inputrec->nLincsIter",-1,ir1->nLincsIter,ir2->nLincsIter);
 622   cmp_real(fp,"inputrec->bd_fric",-1,ir1->bd_fric,ir2->bd_fric,ftol,abstol);
 623   cmp_int(fp,"inputrec->ld_seed",-1,ir1->ld_seed,ir2->ld_seed);
 624   cmp_real(fp,"inputrec->cos_accel",-1,ir1->cos_accel,ir2->cos_accel,ftol,abstol);
 625   cmp_rvec(fp,"inputrec->deform(a)",-1,ir1->deform[XX],ir2->deform[XX],ftol,abstol);
 626   cmp_rvec(fp,"inputrec->deform(b)",-1,ir1->deform[YY],ir2->deform[YY],ftol,abstol);
 627   cmp_rvec(fp,"inputrec->deform(c)",-1,ir1->deform[ZZ],ir2->deform[ZZ],ftol,abstol);
 628   cmp_int(fp,"inputrec->userint1",-1,ir1->userint1,ir2->userint1);
 629   cmp_int(fp,"inputrec->userint2",-1,ir1->userint2,ir2->userint2);
 630   cmp_int(fp,"inputrec->userint3",-1,ir1->userint3,ir2->userint3);
 631   cmp_int(fp,"inputrec->userint4",-1,ir1->userint4,ir2->userint4);
 632   cmp_real(fp,"inputrec->userreal1",-1,ir1->userreal1,ir2->userreal1,ftol,abstol);
 633   cmp_real(fp,"inputrec->userreal2",-1,ir1->userreal2,ir2->userreal2,ftol,abstol);
 634   cmp_real(fp,"inputrec->userreal3",-1,ir1->userreal3,ir2->userreal3,ftol,abstol);
 635   cmp_real(fp,"inputrec->userreal4",-1,ir1->userreal4,ir2->userreal4,ftol,abstol);
 636   cmp_grpopts(fp,&(ir1->opts),&(ir2->opts),ftol,abstol);
 637   cmp_cosines(fp,"ex",ir1->ex,ir2->ex,ftol,abstol);
 638   cmp_cosines(fp,"et",ir1->et,ir2->et,ftol,abstol);
 639 }
 640
 641 static void comp_pull_AB(FILE *fp,t_pull *pull,real ftol,real abstol)
 642 {
 643   int i;
 644
 645   for(i=0; i<pull->ngrp+1; i++) {
 646     fprintf(fp,"comparing pull group %d\n",i);
 647     cmp_real(fp,"pullgrp->k",-1,pull->grp[i].k,pull->grp[i].kB,ftol,abstol);
 648   }
 649 }
 650
 651 static void comp_state(t_state *st1, t_state *st2,
 652                        gmx_bool bRMSD,real ftol,real abstol)
 653 {
 654   int i,j,nc;
 655
 656   fprintf(stdout,"comparing flags\n");
 657   cmp_int(stdout,"flags",-1,st1->flags,st2->flags);
 658   fprintf(stdout,"comparing box\n");
 659   cmp_rvecs(stdout,"box",DIM,st1->box,st2->box,FALSE,ftol,abstol);
 660   fprintf(stdout,"comparing box_rel\n");
 661   cmp_rvecs(stdout,"box_rel",DIM,st1->box_rel,st2->box_rel,FALSE,ftol,abstol);
 662   fprintf(stdout,"comparing boxv\n");
 663   cmp_rvecs(stdout,"boxv",DIM,st1->boxv,st2->boxv,FALSE,ftol,abstol);
 664   if (st1->flags & (1<<estSVIR_PREV)) {
 665     fprintf(stdout,"comparing shake vir_prev\n");
 666     cmp_rvecs(stdout,"svir_prev",DIM,st1->svir_prev,st2->svir_prev,FALSE,ftol,abstol);
 667   }
 668   if (st1->flags & (1<<estFVIR_PREV)) {
 669     fprintf(stdout,"comparing force vir_prev\n");
 670     cmp_rvecs(stdout,"fvir_prev",DIM,st1->fvir_prev,st2->fvir_prev,FALSE,ftol,abstol);
 671   }
 672   if (st1->flags & (1<<estPRES_PREV)) {
 673     fprintf(stdout,"comparing prev_pres\n");
 674     cmp_rvecs(stdout,"pres_prev",DIM,st1->pres_prev,st2->pres_prev,FALSE,ftol,abstol);
 675   }
 676   cmp_int(stdout,"ngtc",-1,st1->ngtc,st2->ngtc);
 677   cmp_int(stdout,"nhchainlength",-1,st1->nhchainlength,st2->nhchainlength);
 678   if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength){
 679     for(i=0; i<st1->ngtc; i++) {
 680       nc = i*st1->nhchainlength;
 681       for(j=0; j<nc; j++) {
 682         cmp_real(stdout,"nosehoover_xi",
 683                  i,st1->nosehoover_xi[nc+j],st2->nosehoover_xi[nc+j],ftol,abstol);
 684       }
 685     }
 686   }
 687   cmp_int(stdout,"nnhpres",-1,st1->nnhpres,st2->nnhpres);
 688   if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength) {
 689     for(i=0; i<st1->nnhpres; i++) {
 690       nc = i*st1->nhchainlength;
 691       for(j=0; j<nc; j++) {
 692         cmp_real(stdout,"nosehoover_xi",
 693                  i,st1->nhpres_xi[nc+j],st2->nhpres_xi[nc+j],ftol,abstol);
 694       }
 695     }
 696   }
 697
 698   cmp_int(stdout,"natoms",-1,st1->natoms,st2->natoms);
 699   if (st1->natoms == st2->natoms) {
 700     if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX))) {
 701       fprintf(stdout,"comparing x\n");
 702       cmp_rvecs(stdout,"x",st1->natoms,st1->x,st2->x,bRMSD,ftol,abstol);
 703     }
 704     if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV))) {
 705       fprintf(stdout,"comparing v\n");
 706       cmp_rvecs(stdout,"v",st1->natoms,st1->v,st2->v,bRMSD,ftol,abstol);
 707     }
 708   }
 709 }
 710
 711 void comp_tpx(const char *fn1,const char *fn2,
 712               gmx_bool bRMSD,real ftol,real abstol)
 713 {
 714   const char  *ff[2];
 715   t_tpxheader sh[2];
 716   t_inputrec  ir[2];
 717   t_state     state[2];
 718   gmx_mtop_t  mtop[2];
 719   t_topology  top[2];
 720   int         i;
 721
 722   ff[0]=fn1;
 723   ff[1]=fn2;
 724   for(i=0; i<(fn2 ? 2 : 1); i++) {
 725     read_tpx_state(ff[i],&(ir[i]),&state[i],NULL,&(mtop[i]));
 726   }
 727   if (fn2) {
 728     cmp_inputrec(stdout,&ir[0],&ir[1],ftol,abstol);
 729     /* Convert gmx_mtop_t to t_topology.
 730      * We should implement direct mtop comparison,
 731      * but it might be useful to keep t_topology comparison as an option.
 732      */
 733     top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
 734     top[1] = gmx_mtop_t_to_t_topology(&mtop[1]);
 735     cmp_top(stdout,&top[0],&top[1],ftol,abstol);
 736     cmp_groups(stdout,&mtop[0].groups,&mtop[1].groups,
 737                mtop[0].natoms,mtop[1].natoms);
 738     comp_state(&state[0],&state[1],bRMSD,ftol,abstol);
 739   } else {
 740     if (ir[0].efep == efepNO) {
 741       fprintf(stdout,"inputrec->efep = %s\n",efep_names[ir[0].efep]);
 742     } else {
 743       if (ir[0].ePull != epullNO) {
 744         comp_pull_AB(stdout,ir->pull,ftol,abstol);
 745       }
 746       /* Convert gmx_mtop_t to t_topology.
 747        * We should implement direct mtop comparison,
 748        * but it might be useful to keep t_topology comparison as an option.
 749        */
 750       top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
 751       cmp_top(stdout,&top[0],NULL,ftol,abstol);
 752     }
 753   }
 754 }
 755
 756 void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
 757                 gmx_bool bRMSD, real ftol,real abstol)
 758 {
 759   fprintf(fp,"\n");
 760   cmp_int(fp,"flags",-1,fr1->flags,fr2->flags);
 761   cmp_int(fp,"not_ok",-1,fr1->not_ok,fr2->not_ok);
 762   cmp_int(fp,"natoms",-1,fr1->natoms,fr2->natoms);
 763   cmp_real(fp,"t0",-1,fr1->t0,fr2->t0,ftol,abstol);
 764   if (cmp_bool(fp,"bTitle",-1,fr1->bTitle,fr2->bTitle))
 765     cmp_str(fp,"title", -1, fr1->title, fr2->title);
 766   if (cmp_bool(fp,"bStep",-1,fr1->bStep,fr2->bStep))
 767     cmp_int(fp,"step",-1,fr1->step,fr2->step);
 768   cmp_int(fp,"step",-1,fr1->step,fr2->step);
 769   if (cmp_bool(fp,"bTime",-1,fr1->bTime,fr2->bTime))
 770     cmp_real(fp,"time",-1,fr1->time,fr2->time,ftol,abstol);
 771   if (cmp_bool(fp,"bLambda",-1,fr1->bLambda,fr2->bLambda))
 772     cmp_real(fp,"lambda",-1,fr1->lambda,fr2->lambda,ftol,abstol);
 773   if (cmp_bool(fp,"bAtoms",-1,fr1->bAtoms,fr2->bAtoms))
 774     cmp_atoms(fp,fr1->atoms,fr2->atoms,ftol,abstol);
 775   if (cmp_bool(fp,"bPrec",-1,fr1->bPrec,fr2->bPrec))
 776     cmp_real(fp,"prec",-1,fr1->prec,fr2->prec,ftol,abstol);
 777   if (cmp_bool(fp,"bX",-1,fr1->bX,fr2->bX))
 778     cmp_rvecs(fp,"x",min(fr1->natoms,fr2->natoms),fr1->x,fr2->x,bRMSD,ftol,abstol);
 779   if (cmp_bool(fp,"bV",-1,fr1->bV,fr2->bV))
 780     cmp_rvecs(fp,"v",min(fr1->natoms,fr2->natoms),fr1->v,fr2->v,bRMSD,ftol,abstol);
 781   if (cmp_bool(fp,"bF",-1,fr1->bF,fr2->bF))
 782     cmp_rvecs(fp,"f",min(fr1->natoms,fr2->natoms),fr1->f,fr2->f,bRMSD,ftol,abstol);
 783   if (cmp_bool(fp,"bBox",-1,fr1->bBox,fr2->bBox))
 784     cmp_rvecs(fp,"box",3,fr1->box,fr2->box,FALSE,ftol,abstol);
 785 }
 786
 787 void comp_trx(const output_env_t oenv,const char *fn1, const char *fn2,
 788               gmx_bool bRMSD,real ftol,real abstol)
 789 {
 790   int i;
 791   const char *fn[2];
 792   t_trxframe fr[2];
 793   t_trxstatus *status[2];
 794   gmx_bool b[2];
 795
 796   fn[0]=fn1;
 797   fn[1]=fn2;
 798   fprintf(stderr,"Comparing trajectory files %s and %s\n",fn1,fn2);
 799   for (i=0; i<2; i++)
 800     b[i] = read_first_frame(oenv,&status[i],fn[i],&fr[i],TRX_READ_X|TRX_READ_V|TRX_READ_F);
 801
 802   if (b[0] && b[1]) {
 803     do {
 804       comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
 805
 806       for (i=0; i<2; i++)
 807         b[i] = read_next_frame(oenv,status[i],&fr[i]);
 808     } while (b[0] && b[1]);
 809
 810     for (i=0; i<2; i++) {
 811       if (b[i] && !b[1-i])
 812         fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn[i],fn[1-i]);
 813       close_trj(status[i]);
 814     }
 815   }
 816   if (!b[0] && !b[1])
 817     fprintf(stdout,"\nBoth files read correctly\n");
 818 }
 819
 820 static real ener_tensor_diag(int n,int *ind1,int *ind2,
 821                              gmx_enxnm_t *enm1,
 822                              int *tensi,int i,
 823                              t_energy e1[],t_energy e2[])
 824 {
 825   int  d1,d2;
 826   int  len;
 827   int  j;
 828   real prod1,prod2;
 829   int  nfound;
 830
 831   d1 = tensi[i]/DIM;
 832   d2 = tensi[i] - d1*DIM;
 833
 834   /* Find the diagonal elements d1 and d2 */
 835   len = strlen(enm1[ind1[i]].name);
 836   prod1 = 1;
 837   prod2 = 1;
 838   nfound = 0;
 839   for(j=0; j<n; j++) {
 840     if (tensi[j] >= 0 &&
 841         strlen(enm1[ind1[j]].name) == len &&
 842         strncmp(enm1[ind1[i]].name,enm1[ind1[j]].name,len-2) == 0 &&
 843         (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2)) {
 844       prod1 *= fabs(e1[ind1[j]].e);
 845       prod2 *= fabs(e2[ind2[j]].e);
 846       nfound++;
 847     }
 848   }
 849
 850   if (nfound == 2) {
 851     return 0.5*(sqrt(prod1) + sqrt(prod2));
 852   } else {
 853     return 0;
 854   }
 855 }
 856
 857 static gmx_bool enernm_equal(const char *nm1,const char *nm2)
 858 {
 859   int len1,len2;
 860
 861   len1 = strlen(nm1);
 862   len2 = strlen(nm2);
 863
 864   /* Remove " (bar)" at the end of a name */
 865   if (len1 > 6 && strcmp(nm1+len1-6," (bar)") == 0) {
 866     len1 -= 6;
 867   }
 868   if (len2 > 6 && strcmp(nm2+len2-6," (bar)") == 0) {
 869     len2 -= 6;
 870   }
 871
 872   return (len1 == len2 && gmx_strncasecmp(nm1,nm2,len1) == 0);
 873 }
 874
 875 static void cmp_energies(FILE *fp,int step1,int step2,
 876                          t_energy e1[],t_energy e2[],
 877                          gmx_enxnm_t *enm1,gmx_enxnm_t *enm2,
 878                          real ftol,real abstol,
 879                          int nre,int *ind1,int *ind2,int maxener)
 880 {
 881   int  i,ii;
 882   int  *tensi,len,d1,d2;
 883   real ftol_i,abstol_i;
 884
 885   snew(tensi,maxener);
 886   /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
 887   for(i=0; (i<maxener); i++) {
 888     ii = ind1[i];
 889     tensi[i] = -1;
 890     len = strlen(enm1[ii].name);
 891     if (len > 3 && enm1[ii].name[len-3] == '-') {
 892       d1 = enm1[ii].name[len-2] - 'X';
 893       d2 = enm1[ii].name[len-1] - 'X';
 894       if (d1 >= 0 && d1 < DIM &&
 895           d2 >= 0 && d2 < DIM) {
 896         tensi[i] = d1*DIM + d2;
 897       }
 898     }
 899   }
 900
 901   for(i=0; (i<maxener); i++) {
 902     /* Check if this is an off-diagonal tensor element */
 903     if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8) {
 904       /* Turn on the relative tolerance check (4 is maximum relative diff.) */
 905       ftol_i = 5;
 906       /* Do the relative tolerance through an absolute tolerance times
 907        * the size of diagonal components of the tensor.
 908        */
 909       abstol_i = ftol*ener_tensor_diag(nre,ind1,ind2,enm1,tensi,i,e1,e2);
 910       if (debug) {
 911         fprintf(debug,"tensor '%s' val %f diag %f\n",
 912                 enm1[i].name,e1[i].e,abstol_i/ftol);
 913       }
 914       if (abstol_i > 0) {
 915         /* We found a diagonal, we need to check with the minimum tolerance */
 916         abstol_i = min(abstol_i,abstol);
 917       } else {
 918         /* We did not find a diagonal, ignore the relative tolerance check */
 919         abstol_i = abstol;
 920       }
 921     } else {
 922       ftol_i   = ftol;
 923       abstol_i = abstol;
 924     }
 925     if (!equal_real(e1[ind1[i]].e,e2[ind2[i]].e,ftol_i,abstol_i)) {
 926       fprintf(fp,"%-15s  step %3d:  %12g,  step %3d: %12g\n",
 927               enm1[ind1[i]].name,
 928               step1,e1[ind1[i]].e,
 929               step2,e2[ind2[i]].e);
 930     }
 931   }
 932
 933   sfree(tensi);
 934 }
 935
 936 #if 0
 937 static void cmp_disres(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
 938 {
 939   int i;
 940   char bav[64],bt[64],bs[22];
 941
 942   cmp_int(stdout,"ndisre",-1,fr1->ndisre,fr2->ndisre);
 943   if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0)) {
 944     sprintf(bav,"step %s: disre rav",gmx_step_str(fr1->step,bs));
 945     sprintf(bt, "step %s: disre  rt",gmx_step_str(fr1->step,bs));
 946     for(i=0; (i<fr1->ndisre); i++) {
 947       cmp_real(stdout,bav,i,fr1->disre_rm3tav[i],fr2->disre_rm3tav[i],ftol,abstol);
 948       cmp_real(stdout,bt ,i,fr1->disre_rt[i]    ,fr2->disre_rt[i]    ,ftol,abstol);
 949     }
 950   }
 951 }
 952 #endif
 953
 954 static void cmp_eblocks(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
 955 {
 956     int i,j,k;
 957     char buf[64],bs[22];
 958
 959     cmp_int(stdout,"nblock",-1,fr1->nblock,fr2->nblock);
 960     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
 961     {
 962         for(j=0; (j<fr1->nblock); j++)
 963         {
 964             t_enxblock *b1, *b2; /* convenience vars */
 965
 966             b1=&(fr1->block[j]);
 967             b2=&(fr2->block[j]);
 968
 969             sprintf(buf,"step %s: block[%d]",gmx_step_str(fr1->step,bs),j);
 970             cmp_int(stdout,buf,-1,b1->nsub,b2->nsub);
 971             cmp_int(stdout,buf,-1,b1->id,b2->id);
 972
 973             if ( (b1->nsub==b2->nsub) && (b1->id==b2->id) )
 974             {
 975                 for(i=0;i<b1->nsub;i++)
 976                 {
 977                     t_enxsubblock *s1, *s2;
 978
 979                     s1=&(b1->sub[i]);
 980                     s2=&(b2->sub[i]);
 981
 982                     cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
 983                     cmp_gmx_large_int(stdout, buf, s1->nr, s2->nr);
 984
 985                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
 986                     {
 987                         switch(s1->type)
 988                         {
 989                             case xdr_datatype_float:
 990                                 for(k=0;k<s1->nr;k++)
 991                                 {
 992                                     cmp_float(stdout, buf, i,
 993                                              s1->fval[k], s2->fval[k],
 994                                              ftol, abstol);
 995                                 }
 996                                 break;
 997                             case xdr_datatype_double:
 998                                 for(k=0;k<s1->nr;k++)
 999                                 {
1000                                     cmp_double(stdout, buf, i,
1001                                              s1->dval[k], s2->dval[k],
1002                                              ftol, abstol);
1003                                 }
1004                                 break;
1005                             case xdr_datatype_int:
1006                                 for(k=0;k<s1->nr;k++)
1007                                 {
1008                                     cmp_int(stdout, buf, i,
1009                                             s1->ival[k], s2->ival[k]);
1010                                 }
1011                                 break;
1012                             case xdr_datatype_large_int:
1013                                 for(k=0;k<s1->nr;k++)
1014                                 {
1015                                     cmp_gmx_large_int(stdout, buf,
1016                                                       s1->lval[k], s2->lval[k]);
1017                                 }
1018                                 break;
1019                             case xdr_datatype_char:
1020                                 for(k=0;k<s1->nr;k++)
1021                                 {
1022                                     cmp_uc(stdout, buf, i,
1023                                            s1->cval[k], s2->cval[k]);
1024                                 }
1025                                 break;
1026                             case xdr_datatype_string:
1027                                 for(k=0;k<s1->nr;k++)
1028                                 {
1029                                     cmp_str(stdout, buf, i,
1030                                             s1->sval[k], s2->sval[k]);
1031                                 }
1032                                 break;
1033                             default:
1034                                 gmx_incons("Unknown data type!!");
1035                         }
1036                     }
1037                 }
1038             }
1039         }
1040     }
1041 }
1042
1043 void comp_enx(const char *fn1,const char *fn2,real ftol,real abstol,const char *lastener)
1044 {
1045   int        nre,nre1,nre2,block;
1046   ener_file_t in1, in2;
1047   int        i,j,maxener,*ind1,*ind2,*have;
1048   char       buf[256];
1049   gmx_enxnm_t *enm1=NULL,*enm2=NULL;
1050   t_enxframe *fr1,*fr2;
1051   gmx_bool       b1,b2;
1052
1053   fprintf(stdout,"comparing energy file %s and %s\n\n",fn1,fn2);
1054
1055   in1 = open_enx(fn1,"r");
1056   in2 = open_enx(fn2,"r");
1057   do_enxnms(in1,&nre1,&enm1);
1058   do_enxnms(in2,&nre2,&enm2);
1059   if (nre1 != nre2) {
1060     fprintf(stdout,"There are %d and %d terms in the energy files\n\n",
1061             nre1,nre2);
1062   } else {
1063     fprintf(stdout,"There are %d terms in the energy files\n\n",nre1);
1064   }
1065
1066   snew(ind1,nre1);
1067   snew(ind2,nre2);
1068   snew(have,nre2);
1069   nre = 0;
1070   for(i=0; i<nre1; i++) {
1071     for(j=0; j<nre2; j++) {
1072       if (enernm_equal(enm1[i].name,enm2[j].name)) {
1073         ind1[nre] = i;
1074         ind2[nre] = j;
1075         have[j] = 1;
1076         nre++;
1077         break;
1078       }
1079     }
1080     if (nre == 0 || ind1[nre-1] != i) {
1081       cmp_str(stdout,"enm",i,enm1[i].name,"-");
1082     }
1083   }
1084   for(i=0; i<nre2; i++) {
1085     if (have[i] == 0) {
1086       cmp_str(stdout,"enm",i,"-",enm2[i].name);
1087     }
1088   }
1089
1090   maxener = nre;
1091   for(i=0; i<nre; i++) {
1092     if ((lastener != NULL) && (strstr(enm1[i].name,lastener) != NULL)) {
1093       maxener=i+1;
1094       break;
1095     }
1096   }
1097
1098   fprintf(stdout,"There are %d terms to compare in the energy files\n\n",
1099           maxener);
1100
1101   for(i=0; i<maxener; i++) {
1102     cmp_str(stdout,"unit",i,enm1[ind1[i]].unit,enm2[ind2[i]].unit);
1103   }
1104
1105   snew(fr1,1);
1106   snew(fr2,1);
1107   do {
1108     b1 = do_enx(in1,fr1);
1109     b2 = do_enx(in2,fr2);
1110     if (b1 && !b2)
1111       fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn2,fn1);
1112     else if (!b1 && b2)
1113       fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn1,fn2);
1114     else if (!b1 && !b2)
1115       fprintf(stdout,"\nFiles read successfully\n");
1116     else {
1117       cmp_real(stdout,"t",-1,fr1->t,fr2->t,ftol,abstol);
1118       cmp_int(stdout,"step",-1,fr1->step,fr2->step);
1119       /* We don't want to print the nre mismatch for every frame */
1120       /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1121       if ((fr1->nre >= nre) && (fr2->nre >= nre))
1122         cmp_energies(stdout,fr1->step,fr1->step,fr1->ener,fr2->ener,
1123                      enm1,enm2,ftol,abstol,nre,ind1,ind2,maxener);
1124       /*cmp_disres(fr1,fr2,ftol,abstol);*/
1125       cmp_eblocks(fr1,fr2,ftol,abstol);
1126     }
1127   } while (b1 && b2);
1128
1129   close_enx(in1);
1130   close_enx(in2);
1131
1132   free_enxframe(fr2);
1133   sfree(fr2);
1134   free_enxframe(fr1);
1135   sfree(fr1);
1136 }