src/kernel/tpbcmp.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_tpbcmp_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdio.h>
  33 #include <string.h>
  34 #include "main.h"
  35 #include "macros.h"
  36 #include "smalloc.h"
  37 #include "futil.h"
  38 #include "statutil.h"
  39 #include "sysstuff.h"
  40 #include "txtdump.h"
  41 #include "fatal.h"
  42 #include "names.h"
  43 #include "tpxio.h"
  44 #include "enxio.h"
  45
  46 static void cmp_int(FILE *fp,char *s,int index,int i1,int i2)
  47 {
  48   if (i1 != i2) {
  49     if (index != -1)
  50       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  51     else
  52       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  53   }
  54 }
  55
  56 static void cmp_us(FILE *fp,char *s,int index,unsigned short i1,unsigned short i2)
  57 {
  58   if (i1 != i2) {
  59     if (index != -1)
  60       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  61     else
  62       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  63   }
  64 }
  65
  66 static void cmp_uc(FILE *fp,char *s,int index,unsigned char i1,unsigned char i2)
  67 {
  68   if (i1 != i2) {
  69     if (index != -1)
  70       fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
  71     else
  72       fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
  73   }
  74 }
  75
  76 static void cmp_real(FILE *fp,char *s,int index,real i1,real i2)
  77 {
  78   if (i1 != i2) {
  79     if (index != -1)
  80       fprintf(fp,"%s[%d] (%e - %e)\n",s,index,i1,i2);
  81     else
  82       fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
  83   }
  84 }
  85
  86 static void cmp_rvec(FILE *fp,char *s,int index,rvec i1,rvec i2)
  87 {
  88   if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ])) {
  89     if (index != -1)
  90       fprintf(fp,"%s[%d] (%e,%e,%e - %e,%e,%e)\n",s,index,
  91               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
  92     else
  93       fprintf(fp,"%s (%e,%e,%e - %e,%e,%e)\n",s,
  94               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
  95   }
  96 }
  97
  98 static void cmp_ivec(FILE *fp,char *s,int index,ivec i1,ivec i2)
  99 {
 100   if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ])) {
 101     if (index != -1)
 102       fprintf(fp,"%s[%d] (%d,%d,%d - %d,%d,%d)\n",s,index,
 103               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 104     else
 105       fprintf(fp,"%s (%d,%d,%d - %d,%d,%d)\n",s,
 106               i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
 107   }
 108 }
 109
 110 static void cmp_ilist(FILE *fp,int ftype,t_ilist *il1,t_ilist *il2)
 111 {
 112   int i;
 113   char buf[256];
 114
 115   fprintf(fp,"comparing ilist %s\n",interaction_function[ftype].name);
 116   sprintf(buf,"%s->nr",interaction_function[ftype].name);
 117   cmp_int(fp,buf,0,il1->nr,il2->nr);
 118   sprintf(buf,"%s->multinr",interaction_function[ftype].name);
 119   for(i=0; (i<MAXPROC); i++)
 120     cmp_int(fp,buf,i,il1->multinr[i],il2->multinr[i]);
 121   sprintf(buf,"%s->iatoms",interaction_function[ftype].name);
 122   for(i=0; (i<il1->nr); i++)
 123     cmp_int(fp,buf,i,il1->iatoms[i],il2->iatoms[i]);
 124 }
 125
 126 void cmp_iparm(FILE *fp,char *s,int index,t_functype ft,
 127                t_iparams ip1,t_iparams ip2)
 128 {
 129   if (memcmp(&ip1,&ip2,(size_t)sizeof(ip1)) != 0) {
 130     fprintf(fp,"%s1: ",s);
 131     pr_iparams(fp,ft,&ip1);
 132     fprintf(fp,"%s2: ",s);
 133     pr_iparams(fp,ft,&ip2);
 134   }
 135 }
 136
 137 static void cmp_idef(FILE *fp,t_idef *id1,t_idef *id2)
 138 {
 139   int i;
 140
 141   fprintf(fp,"comparing idef\n");
 142   cmp_int(fp,"idef->ntypes",-1,id1->ntypes,id2->ntypes);
 143   cmp_int(fp,"idef->pid",   -1,id1->pid,id2->pid);
 144   cmp_int(fp,"idef->atnr",  -1,id1->atnr,id2->atnr);
 145   for(i=0; (i<id1->ntypes); i++) {
 146     cmp_int(fp,"idef->functype",i,(int)id1->functype[i],(int)id2->functype[i]);
 147     cmp_iparm(fp,"idef->iparam",i,id1->functype[i],
 148               id1->iparams[i],id2->iparams[i]);
 149   }
 150   for(i=0; (i<F_NRE); i++)
 151     cmp_ilist(fp,i,&(id1->il[i]),&(id2->il[i]));
 152 }
 153
 154 static void cmp_block(FILE *fp,t_block *b1,t_block *b2,char *s)
 155 {
 156   int i,j,k;
 157   char buf[32];
 158
 159   fprintf(fp,"comparing block %s\n",s);
 160   sprintf(buf,"%s.nr",s);
 161   cmp_int(fp,buf,-1,b1->nr,b2->nr);
 162   sprintf(buf,"%s.nra",s);
 163   cmp_int(fp,buf,-1,b1->nra,b2->nra);
 164   sprintf(buf,"%s.multinr",s);
 165   for(i=0; (i<MAXPROC); i++)
 166     cmp_int(fp,buf,i,b1->multinr[i],b2->multinr[i]);
 167 }
 168
 169 static void cmp_atom(FILE *fp,int index,t_atom *a1,t_atom *a2)
 170 {
 171   int  i;
 172   char buf[256];
 173
 174   cmp_us(fp,"atom.type",index,a1->type,a2->type);
 175   cmp_us(fp,"atom.ptype",index,a1->ptype,a2->ptype);
 176   cmp_int(fp,"atom.resnr",index,a1->resnr,a2->resnr);
 177   cmp_real(fp,"atom.m",index,a1->m,a2->m);
 178   cmp_real(fp,"atom.q",index,a1->q,a2->q);
 179   for(i=0; (i<egcNR); i++) {
 180     sprintf(buf,"atom.grpnr(%d)",i);
 181     cmp_uc(fp,buf,index,a1->grpnr[i],a2->grpnr[i]);
 182   }
 183 }
 184
 185 static void cmp_atoms(FILE *fp,t_atoms *a1,t_atoms *a2)
 186 {
 187   int i;
 188
 189   fprintf(fp,"comparing atoms\n");
 190   cmp_int(fp,"atoms->nr",-1,a1->nr,a2->nr);
 191   for(i=0; (i<a1->nr); i++)
 192     cmp_atom(fp,i,&(a1->atom[i]),&(a2->atom[i]));
 193   cmp_block(fp,&a1->excl,&a2->excl,"excl");
 194 }
 195
 196 static void cmp_top(FILE *fp,t_topology *t1,t_topology *t2)
 197 {
 198   int i;
 199
 200   fprintf(fp,"comparing top\n");
 201   cmp_idef(fp,&(t1->idef),&(t2->idef));
 202   cmp_atoms(fp,&(t1->atoms),&(t2->atoms));
 203   for(i=0; (i<ebNR); i++)
 204     cmp_block(fp,&t1->blocks[i],&t2->blocks[i],EBLOCKS(i));
 205 }
 206
 207 static void cmp_rvecs(FILE *fp,char *title,int n,rvec x1[],rvec x2[])
 208 {
 209   int i;
 210
 211   fprintf(fp,"comparing %s\n",title);
 212   for(i=0; (i<n); i++)
 213     cmp_rvec(fp,title,i,x1[i],x2[i]);
 214 }
 215
 216 static void cmp_grpopts(FILE *fp,t_grpopts *opt1,t_grpopts *opt2)
 217 {
 218   int i;
 219
 220   cmp_int(fp,"inputrec->grpopts.ngtc",0,  opt1->ngtc,opt2->ngtc);
 221   cmp_int(fp,"inputrec->grpopts.ngacc",0, opt1->ngacc,opt2->ngacc);
 222   cmp_int(fp,"inputrec->grpopts.ngfrz",0, opt1->ngfrz,opt2->ngfrz);
 223   cmp_int(fp,"inputrec->grpopts.ngener",0,opt1->ngener,opt2->ngener);
 224   for(i=0; (i<min(opt1->ngtc,opt2->ngtc)); i++) {
 225     cmp_int(fp,"inputrec->grpopts.nrdf",i,opt1->nrdf[i],opt2->nrdf[i]);
 226     cmp_real(fp,"inputrec->grpopts.ref_t",i,opt1->ref_t[i],opt2->ref_t[i]);
 227     cmp_real(fp,"inputrec->grpopts.tau_t",i,opt1->tau_t[i],opt2->tau_t[i]);
 228   }
 229   for(i=0; (i<min(opt1->ngacc,opt2->ngacc)); i++)
 230     cmp_rvec(fp,"inputrec->grpopts.acc",i,opt1->acc[i],opt2->acc[i]);
 231   for(i=0; (i<min(opt1->ngfrz,opt2->ngfrz)); i++)
 232     cmp_ivec(fp,"inputrec->grpopts.nFreeze",i,opt1->nFreeze[i],opt2->nFreeze[i]);
 233 }
 234
 235 static void cmp_cosines(FILE *fp,char *s,t_cosines c1[DIM],t_cosines c2[DIM])
 236 {
 237   int i,m;
 238   char buf[256];
 239
 240   for(m=0; (m<DIM); m++) {
 241     sprintf(buf,"inputrec->%s[%d]",s,m);
 242     cmp_int(fp,buf,0,c1->n,c2->n);
 243     for(i=0; (i<min(c1->n,c2->n)); i++) {
 244       cmp_real(fp,buf,i,c1->a[i],c2->a[i]);
 245       cmp_real(fp,buf,i,c1->phi[i],c2->phi[i]);
 246     }
 247   }
 248 }
 249
 250 static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2)
 251 {
 252   fprintf(fp,"comparing inputrec\n");
 253 #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 254 #define CII(s) cmp_int (fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 255 #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
 256   CII(eI);
 257   CII(nsteps);
 258   CII(ns_type);
 259   CII(nstlist);
 260   CII(ndelta);
 261   CIB(bDomDecomp);
 262   CII(decomp_dir);
 263   CII(nstcomm);
 264   CII(nstlog);
 265   CII(nstxout);
 266   CII(nstvout);
 267   CII(nstfout);
 268   CII(nstenergy);
 269   CII(nstxtcout);
 270   CIR(init_t);
 271   CIR(delta_t);
 272   CIR(xtcprec);
 273   CII(solvent_opt);
 274   CII(nsatoms);
 275   CII(nkx);
 276   CII(nky);
 277   CII(nkz);
 278   CII(pme_order);
 279   CIR(ewald_rtol);
 280   CIB(bOptFFT);
 281   CII(eBox);
 282   CIB(bUncStart);
 283   CIB(btc);
 284   CII(ntcmemory);
 285   CII(epc);
 286   CII(npcmemory);
 287   CIR(tau_p);
 288   cmp_rvec(fp,"inputrec->ref_p",0,ir1->ref_p,ir2->ref_p);
 289   cmp_rvec(fp,"inputrec->compress",0,ir1->compress,ir2->compress);
 290   CIB(bSimAnn);
 291   CIR(zero_temp_time);
 292   CIR(rlist);
 293   CII(coulombtype);
 294   CIR(rcoulomb_switch);
 295   CIR(rcoulomb);
 296   CII(vdwtype);
 297   CIR(rvdw_switch);
 298   CIR(rvdw);
 299   CIR(epsilon_r);
 300   CII(bDispCorr);
 301   CIR(shake_tol);
 302   CIR(fudgeQQ);
 303   CIB(bPert);
 304   CIR(init_lambda);
 305   CIR(delta_lambda);
 306   CIR(dr_fc);
 307   CII(eDisreWeighting);
 308   CIB(bDisreMixed);
 309   CII(nstdisreout);
 310   CIR(dr_tau);
 311   CIR(em_stepsize);
 312   CIR(em_tol);
 313   CII(nstcgsteep);
 314   CII(eConstrAlg);
 315   CII(nProjOrder);
 316   CIR(LincsWarnAngle);
 317   CII(nstLincsout);
 318   CIR(ld_temp);
 319   CIR(ld_fric);
 320   CII(ld_seed);
 321   CII(userint1);
 322   CII(userint2);
 323   CII(userint3);
 324   CII(userint4);
 325   CIR(userreal1);
 326   CIR(userreal2);
 327   CIR(userreal3);
 328   CIR(userreal4);
 329   cmp_grpopts(fp,&(ir1->opts),&(ir2->opts));
 330   cmp_cosines(fp,"ex",ir1->ex,ir2->ex);
 331   cmp_cosines(fp,"et",ir1->et,ir2->et);
 332 }
 333
 334
 335 void comp_tpx(char *fn1,char *fn2)
 336 {
 337   char        *ff[2];
 338   t_tpxheader sh[2];
 339   t_inputrec  ir[2];
 340   rvec        *xx[2],*vv[2];
 341   t_topology  top[2];
 342   matrix      box[2];
 343   int         i,step,natoms;
 344   real        t,lambda;
 345
 346   ff[0]=fn1;
 347   ff[1]=fn2;
 348   for(i=0; (i<2); i++) {
 349     read_tpxheader(ff[i],&(sh[i]));
 350     snew(xx[i],sh[i].natoms);
 351     snew(vv[i],sh[i].natoms);
 352     read_tpx(ff[i],&step,&t,&lambda,&(ir[i]),box[i],&natoms,
 353              xx[i],vv[i],NULL,&(top[i]));
 354   }
 355   cmp_inputrec(stdout,&ir[0],&ir[1]);
 356   cmp_top(stdout,&top[0],&top[1]);
 357   cmp_rvecs(stdout,"box",DIM,box[0],box[1]);
 358   cmp_rvecs(stdout,"x",natoms,xx[0],xx[1]);
 359   cmp_rvecs(stdout,"v",natoms,vv[0],vv[1]);
 360 }
 361
 362 static void cmp_energies(FILE *fp,int frame1,int frame2,int nre,
 363                          t_energy e1[],t_energy e2[],char *enm1[],char *enm2[],
 364                          real ftol)
 365 {
 366   int  i;
 367
 368   for(i=0; (i<nre); i++)
 369     if (fabs(e1[i].e - e2[i].e) > ftol)
 370       fprintf(fp,"%-15s  frame %3d:  %12g,  frame %3d: %12g\n",
 371               enm1[i],frame1,e1[i].e,frame2,e2[i].e);
 372 }
 373
 374 void comp_enx(char *fn1,char *fn2,real ftol)
 375 {
 376   int       in1,in2,nre,nre1,nre2,frame1,frame2,ndr1,ndr2;
 377   int       i;
 378   char      **enm1=NULL,**enm2=NULL;
 379   t_energy  *ee1,*ee2;
 380   t_drblock dr1,dr2;
 381   bool      b1,b2;
 382   real      t1,t2;
 383
 384   fprintf(stdout,"comparing energy file %s and %s\n\n",fn1,fn2);
 385
 386   in1 = open_enx(fn1,"r");
 387   in2 = open_enx(fn2,"r");
 388   do_enxnms(in1,&nre1,&enm1);
 389   do_enxnms(in2,&nre2,&enm2);
 390   if (nre1 != nre2) {
 391     fprintf(stdout,"%s: nre=%d, %s: nre=%d\n",fn1,nre1,fn2,nre2);
 392     return;
 393   }
 394   nre = nre1;
 395   fprintf(stdout,"There are %d terms in the energy files\n\n",nre);
 396
 397   for(i=0; (i<nre); i++) {
 398     if (strcmp(enm1[i],enm2[i]) != 0)
 399       fprintf(stdout,"%s: %16s - %s: %16s\n",fn1,enm1[i],fn2,enm2[i]);
 400   }
 401
 402   snew(ee1,nre);
 403   snew(ee2,nre);
 404   do {
 405     b1 = do_enx(in1,&t1,&frame1,&nre1,ee1,&ndr1,&dr1);
 406     b2 = do_enx(in2,&t2,&frame2,&nre2,ee2,&ndr2,&dr2);
 407     if (b1 && !b2)
 408       fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn2,fn1);
 409     else if (!b1 && b2)
 410       fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn1,fn2);
 411     else if (!b1 && !b2)
 412       fprintf(stdout,"\nFiles read succesfully\n");
 413     else {
 414       cmp_real(stdout,"t",-1,t1,t2);
 415       cmp_int(stdout,"frame",-1,frame1,frame2);
 416       cmp_int(stdout,"nre",-1,nre1,nre2);
 417       cmp_int(stdout,"ndr",-1,ndr1,ndr2);
 418
 419       if (nre1 != nre) {
 420         fprintf(stdout,
 421                 "file %s internally inconsistent: nre changed from %d to %d\n",
 422                 fn1,nre,nre1);
 423         break;
 424       }
 425       if (nre2 != nre) {
 426         fprintf(stdout,
 427                 "file %s internally inconsistent: nre changed from %d to %d\n",
 428                 fn2,nre,nre2);
 429         break;
 430       }
 431       cmp_energies(stdout,frame1,frame2,nre,ee1,ee2,enm1,enm2,ftol);
 432       /*if ((ndr1 == ndr2) && (ndr1 > 0))
 433         cmp_disres(stdout,frame1,frame2,ndr1,dr1,dr2);*/
 434     }
 435   } while (b1 && b2);
 436
 437   close_enx(in1);
 438   close_enx(in2);
 439
 440 }