src/tools/gmx_disre.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, 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "copyrite.h"
  47 #include "mshift.h"
  48 #include "xvgr.h"
  49 #include "vec.h"
  50 #include "do_fit.h"
  51 #include "confio.h"
  52 #include "smalloc.h"
  53 #include "nrnb.h"
  54 #include "disre.h"
  55 #include "statutil.h"
  56 #include "force.h"
  57 #include "gstat.h"
  58 #include "main.h"
  59 #include "pdbio.h"
  60 #include "index.h"
  61 #include "mdatoms.h"
  62 #include "tpxio.h"
  63 #include "mdrun.h"
  64 #include "names.h"
  65 #include "matio.h"
  66 #include "mtop_util.h"
  67 #include "gmx_ana.h"
  68
  69
  70 typedef struct {
  71   int n;
  72   real v;
  73 } t_toppop;
  74
  75 t_toppop *top=NULL;
  76 int      ntop=0;
  77
  78 typedef struct {
  79   int  nv,nframes;
  80   real sumv,averv,maxv;
  81   real *aver1,*aver2,*aver_3,*aver_6;
  82 } t_dr_result;
  83
  84 static void init5(int n)
  85 {
  86   ntop=n;
  87   snew(top,ntop);
  88 }
  89
  90 static void reset5(void)
  91 {
  92   int i;
  93
  94   for(i=0; (i<ntop); i++) {
  95     top[i].n=-1;
  96     top[i].v= 0;
  97   }
  98 }
  99
 100 int tpcomp(const void *a,const void *b)
 101 {
 102   t_toppop *tpa;
 103   t_toppop *tpb;
 104
 105   tpa=(t_toppop *)a;
 106   tpb=(t_toppop *)b;
 107
 108   return  (1e7*(tpb->v-tpa->v));
 109 }
 110
 111 static void add5(int ndr,real viol)
 112 {
 113   int i,mini;
 114
 115   mini=0;
 116   for(i=1; (i<ntop); i++)
 117     if (top[i].v < top[mini].v)
 118       mini=i;
 119   if (viol > top[mini].v) {
 120     top[mini].v=viol;
 121     top[mini].n=ndr;
 122   }
 123 }
 124
 125 static void print5(FILE *fp)
 126 {
 127   int i;
 128
 129   qsort(top,ntop,sizeof(top[0]),tpcomp);
 130   fprintf(fp,"Index:");
 131   for(i=0; (i<ntop); i++)
 132     fprintf(fp," %6d",top[i].n);
 133   fprintf(fp,"\nViol: ");
 134   for(i=0; (i<ntop); i++)
 135     fprintf(fp," %6.3f",top[i].v);
 136   fprintf(fp,"\n");
 137 }
 138
 139 static void check_viol(FILE *log,t_commrec *cr,
 140                        t_ilist *disres,t_iparams forceparams[],
 141                        t_functype functype[],rvec x[],rvec f[],
 142                        t_forcerec *fr,t_pbc *pbc,t_graph *g,t_dr_result dr[],
 143                        int clust_id,int isize,atom_id index[],real vvindex[],
 144                        t_fcdata *fcd)
 145 {
 146   t_iatom *forceatoms;
 147   int     i,j,nat,n,type,nviol,ndr,label;
 148   real    ener,rt,mviol,tviol,viol,lam,dvdl,drt;
 149   static  gmx_bool bFirst=TRUE;
 150
 151   lam  =0;
 152   dvdl =0;
 153   tviol=0;
 154   nviol=0;
 155   mviol=0;
 156   ndr=0;
 157   if (ntop)
 158     reset5();
 159   forceatoms=disres->iatoms;
 160   for(j=0; (j<isize); j++) {
 161     vvindex[j]=0;
 162   }
 163   nat = interaction_function[F_DISRES].nratoms+1;
 164   for(i=0; (i<disres->nr); ) {
 165     type  = forceatoms[i];
 166     n     = 0;
 167     label = forceparams[type].disres.label;
 168     if (debug)
 169       fprintf(debug,"DISRE: ndr = %d, label = %d  i=%d, n =%d\n",
 170               ndr,label,i,n);
 171     if (ndr != label)
 172       gmx_fatal(FARGS,"tpr inconsistency. ndr = %d, label = %d\n",ndr,label);
 173     do {
 174       n += nat;
 175     } while (((i+n) < disres->nr) &&
 176              (forceparams[forceatoms[i+n]].disres.label == label));
 177
 178     calc_disres_R_6(cr->ms,n,&forceatoms[i],forceparams,
 179                     (const rvec*)x,pbc,fcd,NULL);
 180
 181     if (fcd->disres.Rt_6[0] <= 0)
 182       gmx_fatal(FARGS,"ndr = %d, rt_6 = %f",ndr,fcd->disres.Rt_6[0]);
 183
 184     rt = pow(fcd->disres.Rt_6[0],-1.0/6.0);
 185     dr[clust_id].aver1[ndr]  += rt;
 186     dr[clust_id].aver2[ndr]  += sqr(rt);
 187     drt = pow(rt,-3.0);
 188     dr[clust_id].aver_3[ndr] += drt;
 189     dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
 190
 191     ener=interaction_function[F_DISRES].ifunc(n,&forceatoms[i],forceparams,
 192                                               (const rvec*)x,f,fr->fshift,
 193                                               pbc,g,lam,&dvdl,NULL,fcd,NULL);
 194     viol = fcd->disres.sumviol;
 195
 196     if (viol > 0) {
 197       nviol++;
 198       if (ntop)
 199         add5(forceparams[type].disres.label,viol);
 200       if (viol > mviol)
 201         mviol = viol;
 202       tviol += viol;
 203       for(j=0; (j<isize); j++) {
 204         if (index[j] == forceparams[type].disres.label)
 205           vvindex[j]=pow(fcd->disres.Rt_6[0],-1.0/6.0);
 206         }
 207     }
 208     ndr ++;
 209     i   += n;
 210   }
 211   dr[clust_id].nv   = nviol;
 212   dr[clust_id].maxv = mviol;
 213   dr[clust_id].sumv = tviol;
 214   dr[clust_id].averv= tviol/ndr;
 215   dr[clust_id].nframes++;
 216
 217   if (bFirst) {
 218     fprintf(stderr,"\nThere are %d restraints and %d pairs\n",
 219             ndr,disres->nr/nat);
 220     bFirst = FALSE;
 221   }
 222   if (ntop)
 223     print5(log);
 224 }
 225
 226 typedef struct {
 227   int  index;
 228   gmx_bool bCore;
 229   real up1,r,rT3,rT6,viol,violT3,violT6;
 230 } t_dr_stats;
 231
 232 static int drs_comp(const void *a,const void *b)
 233 {
 234   t_dr_stats *da,*db;
 235
 236   da = (t_dr_stats *)a;
 237   db = (t_dr_stats *)b;
 238
 239   if (da->viol > db->viol)
 240     return -1;
 241   else if (da->viol < db->viol)
 242     return 1;
 243   else
 244     return 0;
 245 }
 246
 247 static void dump_dump(FILE *log,int ndr,t_dr_stats drs[])
 248 {
 249   static const char *core[] = { "All restraints", "Core restraints" };
 250   static const char *tp[]   = { "linear", "third power", "sixth power" };
 251   real viol_tot,viol_max,viol=0;
 252   gmx_bool bCore;
 253   int  nviol,nrestr;
 254   int  i,kkk;
 255
 256   for(bCore = FALSE; (bCore <= TRUE); bCore++) {
 257     for(kkk=0; (kkk<3); kkk++) {
 258       viol_tot  = 0;
 259       viol_max  = 0;
 260       nviol     = 0;
 261       nrestr    = 0;
 262       for(i=0; (i<ndr); i++) {
 263         if (!bCore || (bCore && drs[i].bCore)) {
 264           switch (kkk) {
 265           case 0:
 266             viol = drs[i].viol;
 267             break;
 268           case 1:
 269             viol = drs[i].violT3;
 270             break;
 271           case 2:
 272             viol = drs[i].violT6;
 273             break;
 274           default:
 275             gmx_incons("Dumping violations");
 276           }
 277           viol_max     = max(viol_max,viol);
 278           if (viol > 0)
 279             nviol++;
 280           viol_tot  += viol;
 281           nrestr++;
 282         }
 283       }
 284       if ((nrestr > 0) || (bCore && (nrestr < ndr))) {
 285         fprintf(log,"\n");
 286         fprintf(log,"+++++++ %s ++++++++\n",core[bCore]);
 287         fprintf(log,"+++++++ Using %s averaging: ++++++++\n",tp[kkk]);
 288         fprintf(log,"Sum of violations: %8.3f nm\n",viol_tot);
 289         if (nrestr > 0)
 290           fprintf(log,"Average violation: %8.3f nm\n",viol_tot/nrestr);
 291         fprintf(log,"Largest violation: %8.3f nm\n",viol_max);
 292         fprintf(log,"Number of violated restraints: %d/%d\n",nviol,nrestr);
 293       }
 294     }
 295   }
 296 }
 297
 298 static void dump_viol(FILE *log,int ndr,t_dr_stats *drs,gmx_bool bLinear)
 299 {
 300   int i;
 301
 302   fprintf(log,"Restr. Core     Up1     <r>   <rT3>   <rT6>  <viol><violT3><violT6>\n");
 303   for(i=0; (i<ndr); i++) {
 304     if (bLinear  && (drs[i].viol == 0))
 305       break;
 306     fprintf(log,"%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
 307             drs[i].index,yesno_names[drs[i].bCore],
 308             drs[i].up1,drs[i].r,drs[i].rT3,drs[i].rT6,
 309             drs[i].viol,drs[i].violT3,drs[i].violT6);
 310   }
 311 }
 312
 313 static gmx_bool is_core(int i,int isize,atom_id index[])
 314 {
 315   int kk;
 316   gmx_bool bIC = FALSE;
 317
 318   for(kk=0; !bIC && (kk<isize); kk++)
 319     bIC = (index[kk] == i);
 320
 321   return bIC;
 322 }
 323
 324 static void dump_stats(FILE *log,int nsteps,int ndr,t_ilist *disres,
 325                        t_iparams ip[],t_dr_result *dr,
 326                        int isize,atom_id index[],t_atoms *atoms)
 327 {
 328   int  i,j,nra;
 329   t_dr_stats *drs;
 330
 331   fprintf(log,"\n");
 332   fprintf(log,"++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
 333   snew(drs,ndr);
 334   j         = 0;
 335   nra       = interaction_function[F_DISRES].nratoms+1;
 336   for(i=0; (i<ndr); i++) {
 337     /* Search for the appropriate restraint */
 338     for( ; (j<disres->nr); j+=nra) {
 339       if (ip[disres->iatoms[j]].disres.label == i)
 340         break;
 341     }
 342     drs[i].index  = i;
 343     drs[i].bCore  = is_core(i,isize,index);
 344     drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 345     drs[i].r      = dr->aver1[i]/nsteps;
 346     drs[i].rT3    = pow(dr->aver_3[i]/nsteps,-1.0/3.0);
 347     drs[i].rT6    = pow(dr->aver_6[i]/nsteps,-1.0/6.0);
 348     drs[i].viol   = max(0,drs[i].r-drs[i].up1);
 349     drs[i].violT3 = max(0,drs[i].rT3-drs[i].up1);
 350     drs[i].violT6 = max(0,drs[i].rT6-drs[i].up1);
 351     if (atoms) {
 352       int j1 = disres->iatoms[j+1];
 353       int j2 = disres->iatoms[j+2];
 354       atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
 355       atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
 356     }
 357   }
 358   dump_viol(log,ndr,drs,FALSE);
 359
 360   fprintf(log,"+++ Sorted by linear averaged violations: +++\n");
 361   qsort(drs,ndr,sizeof(drs[0]),drs_comp);
 362   dump_viol(log,ndr,drs,TRUE);
 363
 364   dump_dump(log,ndr,drs);
 365
 366   sfree(drs);
 367 }
 368
 369 static void dump_clust_stats(FILE *fp,int ndr,t_ilist *disres,
 370                              t_iparams ip[],t_blocka *clust,t_dr_result dr[],
 371                              char *clust_name[],int isize,atom_id index[])
 372 {
 373   int    i,j,k,nra,mmm=0;
 374   double sumV,maxV,sumVT3,sumVT6,maxVT3,maxVT6;
 375   t_dr_stats *drs;
 376
 377   fprintf(fp,"\n");
 378   fprintf(fp,"++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
 379   fprintf(fp,"Cluster  NFrames    SumV      MaxV     SumVT     MaxVT     SumVS     MaxVS\n");
 380
 381   snew(drs,ndr);
 382
 383   for(k=0; (k<clust->nr); k++) {
 384     if (dr[k].nframes == 0)
 385       continue;
 386     if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
 387       gmx_fatal(FARGS,"Inconsistency in cluster %s.\n"
 388                 "Found %d frames in trajectory rather than the expected %d\n",
 389                 clust_name[k],dr[k].nframes,
 390                 clust->index[k+1]-clust->index[k]);
 391     if (!clust_name[k])
 392       gmx_fatal(FARGS,"Inconsistency with cluster %d. Invalid name",k);
 393     j         = 0;
 394     nra       = interaction_function[F_DISRES].nratoms+1;
 395     sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
 396     for(i=0; (i<ndr); i++) {
 397       /* Search for the appropriate restraint */
 398       for( ; (j<disres->nr); j+=nra) {
 399         if (ip[disres->iatoms[j]].disres.label == i)
 400           break;
 401       }
 402       drs[i].index  = i;
 403       drs[i].bCore  = is_core(i,isize,index);
 404       drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 405       drs[i].r      = dr[k].aver1[i]/dr[k].nframes;
 406       if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
 407         gmx_fatal(FARGS,"dr[%d].aver_3[%d] = %f",k,i,dr[k].aver_3[i]);
 408       drs[i].rT3    = pow(dr[k].aver_3[i]/dr[k].nframes,-1.0/3.0);
 409       drs[i].rT6    = pow(dr[k].aver_6[i]/dr[k].nframes,-1.0/6.0);
 410       drs[i].viol   = max(0,drs[i].r-drs[i].up1);
 411       drs[i].violT3 = max(0,drs[i].rT3-drs[i].up1);
 412       drs[i].violT6 = max(0,drs[i].rT6-drs[i].up1);
 413       sumV   += drs[i].viol;
 414       sumVT3 += drs[i].violT3;
 415       sumVT6 += drs[i].violT6;
 416       maxV    = max(maxV,drs[i].viol);
 417       maxVT3  = max(maxVT3,drs[i].violT3);
 418       maxVT6  = max(maxVT6,drs[i].violT6);
 419     }
 420     if (strcmp(clust_name[k],"1000") == 0) {
 421       mmm ++;
 422     }
 423     fprintf(fp,"%-10s%6d%8.3f  %8.3f  %8.3f  %8.3f  %8.3f  %8.3f\n",
 424             clust_name[k],
 425             dr[k].nframes,sumV,maxV,sumVT3,maxVT3,sumVT6,maxVT6);
 426
 427   }
 428   fflush(fp);
 429   sfree(drs);
 430 }
 431
 432 static void init_dr_res(t_dr_result *dr,int ndr)
 433 {
 434   snew(dr->aver1,ndr+1);
 435   snew(dr->aver2,ndr+1);
 436   snew(dr->aver_3,ndr+1);
 437   snew(dr->aver_6,ndr+1);
 438   dr->nv      = 0;
 439   dr->nframes = 0;
 440   dr->sumv    = 0;
 441   dr->maxv    = 0;
 442   dr->averv   = 0;
 443 }
 444
 445 static void dump_disre_matrix(const char *fn,t_dr_result *dr,int ndr,
 446                               int nsteps,t_idef *idef,gmx_mtop_t *mtop,
 447                               real max_dr,int nlevels,gmx_bool bThird)
 448 {
 449   FILE     *fp;
 450   int      *resnr;
 451   int      n_res,a_offset,mb,mol,a;
 452   t_atoms  *atoms;
 453   int      iii,i,j,nra,nratoms,tp,ri,rj,index,nlabel,label;
 454   atom_id  ai,aj,*ptr;
 455   real     **matrix,*t_res,hi,*w_dr,rav,rviol;
 456   t_rgb    rlo = { 1, 1, 1 };
 457   t_rgb    rhi = { 0, 0, 0 };
 458   if (fn == NULL)
 459     return;
 460
 461   snew(resnr,mtop->natoms);
 462   n_res = 0;
 463   a_offset = 0;
 464   for(mb=0; mb<mtop->nmolblock; mb++) {
 465     atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 466     for(mol=0; mol<mtop->molblock[mb].nmol; mol++) {
 467       for(a=0; a<atoms->nr; a++) {
 468         resnr[a_offset+a] = n_res + atoms->atom[a].resind;
 469       }
 470       n_res    += atoms->nres;
 471       a_offset += atoms->nr;
 472     }
 473   }
 474
 475   snew(t_res,n_res);
 476   for(i=0; (i<n_res); i++)
 477     t_res[i] = i+1;
 478   snew(matrix,n_res);
 479   for(i=0; (i<n_res); i++)
 480     snew(matrix[i],n_res);
 481   nratoms = interaction_function[F_DISRES].nratoms;
 482   nra = (idef->il[F_DISRES].nr/(nratoms+1));
 483   snew(ptr,nra+1);
 484   index   = 0;
 485   nlabel  = 0;
 486   ptr[0]  = 0;
 487   snew(w_dr,ndr);
 488   for(i=0; (i<idef->il[F_DISRES].nr); i+=nratoms+1) {
 489     tp       = idef->il[F_DISRES].iatoms[i];
 490     label    = idef->iparams[tp].disres.label;
 491
 492     if (label != index) {
 493       /* Set index pointer */
 494       ptr[index+1] = i;
 495       if (nlabel <= 0)
 496         gmx_fatal(FARGS,"nlabel is %d, label = %d",nlabel,label);
 497       if (index >= ndr)
 498         gmx_fatal(FARGS,"ndr = %d, index = %d",ndr,index);
 499       /* Update the weight */
 500       w_dr[index] = 1.0/nlabel;
 501       index  = label;
 502       nlabel = 1;
 503     }
 504     else {
 505       nlabel++;
 506     }
 507   }
 508   printf("nlabel = %d, index = %d, ndr = %d\n",nlabel,index,ndr);
 509   hi = 0;
 510   for(i=0; (i<ndr); i++) {
 511     for(j=ptr[i]; (j<ptr[i+1]); j+=nratoms+1) {
 512       tp  = idef->il[F_DISRES].iatoms[j];
 513       ai  = idef->il[F_DISRES].iatoms[j+1];
 514       aj  = idef->il[F_DISRES].iatoms[j+2];
 515
 516       ri = resnr[ai];
 517       rj = resnr[aj];
 518       if (bThird)
 519         rav = pow(dr->aver_3[i]/nsteps,-1.0/3.0);
 520       else
 521         rav = dr->aver1[i]/nsteps;
 522       if (debug)
 523         fprintf(debug,"DR %d, atoms %d, %d, distance %g\n",i,ai,aj,rav);
 524       rviol = max(0,rav-idef->iparams[tp].disres.up1);
 525       matrix[ri][rj] += w_dr[i]*rviol;
 526       matrix[rj][ri] += w_dr[i]*rviol;
 527       hi = max(hi,matrix[ri][rj]);
 528       hi = max(hi,matrix[rj][ri]);
 529     }
 530   }
 531
 532   sfree(resnr);
 533
 534   if (max_dr > 0) {
 535     if (hi > max_dr)
 536       printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n",max_dr,hi);
 537     hi = max_dr;
 538   }
 539   printf("Highest level in the matrix will be %g\n",hi);
 540   fp = ffopen(fn,"w");
 541   write_xpm(fp,0,"Distance Violations","<V> (nm)","Residue","Residue",
 542             n_res,n_res,t_res,t_res,matrix,0,hi,rlo,rhi,&nlevels);
 543   ffclose(fp);
 544 }
 545
 546 int gmx_disre(int argc,char *argv[])
 547 {
 548   const char *desc[] = {
 549     "[TT]g_disre[tt] computes violations of distance restraints.",
 550     "If necessary, all protons can be added to a protein molecule ",
 551     "using the [TT]g_protonate[tt] program.[PAR]",
 552     "The program always",
 553     "computes the instantaneous violations rather than time-averaged,",
 554     "because this analysis is done from a trajectory file afterwards",
 555     "it does not make sense to use time averaging. However,",
 556     "the time averaged values per restraint are given in the log file.[PAR]",
 557     "An index file may be used to select specific restraints for",
 558     "printing.[PAR]",
 559     "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
 560     "amount of average violations.[PAR]",
 561     "When the [TT]-c[tt] option is given, an index file will be read",
 562     "containing the frames in your trajectory corresponding to the clusters",
 563     "(defined in another manner) that you want to analyze. For these clusters",
 564     "the program will compute average violations using the third power",
 565     "averaging algorithm and print them in the log file."
 566   };
 567   static int  ntop      = 0;
 568   static int  nlevels   = 20;
 569   static real max_dr    = 0;
 570   static gmx_bool bThird    = TRUE;
 571   t_pargs pa[] = {
 572     { "-ntop", FALSE, etINT,  {&ntop},
 573       "Number of large violations that are stored in the log file every step" },
 574     { "-maxdr", FALSE, etREAL, {&max_dr},
 575       "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
 576     { "-nlevels", FALSE, etINT, {&nlevels},
 577       "Number of levels in the matrix output" },
 578     { "-third", FALSE, etBOOL, {&bThird},
 579       "Use inverse third power averaging or linear for matrix output" }
 580   };
 581
 582   FILE        *out=NULL,*aver=NULL,*numv=NULL,*maxxv=NULL,*xvg=NULL;
 583   t_tpxheader header;
 584   t_inputrec  ir;
 585   gmx_mtop_t  mtop;
 586   rvec        *xtop;
 587   gmx_localtop_t *top;
 588   t_atoms     *atoms=NULL;
 589   t_forcerec  *fr;
 590   t_fcdata    fcd;
 591   t_nrnb      nrnb;
 592   t_commrec   *cr;
 593   t_graph     *g;
 594   int         ntopatoms,natoms,i,j,kkk;
 595   t_trxstatus *status;
 596   real        t;
 597   rvec        *x,*f,*xav=NULL;
 598   matrix      box;
 599   gmx_bool        bPDB;
 600   int         isize;
 601   atom_id     *index=NULL,*ind_fit=NULL;
 602   char        *grpname;
 603   t_cluster_ndx *clust=NULL;
 604   t_dr_result dr,*dr_clust=NULL;
 605   char        **leg;
 606   real        *vvindex=NULL,*w_rls=NULL;
 607   t_mdatoms   *mdatoms;
 608   t_pbc       pbc,*pbc_null;
 609   int         my_clust;
 610   FILE        *fplog;
 611   output_env_t oenv;
 612   gmx_rmpbc_t  gpbc=NULL;
 613
 614   t_filenm fnm[] = {
 615     { efTPX, NULL, NULL, ffREAD },
 616     { efTRX, "-f", NULL, ffREAD },
 617     { efXVG, "-ds", "drsum",  ffWRITE },
 618     { efXVG, "-da", "draver", ffWRITE },
 619     { efXVG, "-dn", "drnum",  ffWRITE },
 620     { efXVG, "-dm", "drmax",  ffWRITE },
 621     { efXVG, "-dr", "restr",  ffWRITE },
 622     { efLOG, "-l",  "disres", ffWRITE },
 623     { efNDX, NULL,  "viol",   ffOPTRD },
 624     { efPDB, "-q",  "viol",   ffOPTWR },
 625     { efNDX, "-c",  "clust",  ffOPTRD },
 626     { efXPM, "-x",  "matrix", ffOPTWR }
 627   };
 628 #define NFILE asize(fnm)
 629
 630   cr  = init_par(&argc,&argv);
 631   CopyRight(stderr,argv[0]);
 632   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 633                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 634
 635   gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,FALSE,0,&fplog);
 636
 637   if (ntop)
 638     init5(ntop);
 639
 640   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header,FALSE,NULL,NULL);
 641   snew(xtop,header.natoms);
 642   read_tpx(ftp2fn(efTPX,NFILE,fnm),&ir,box,&ntopatoms,xtop,NULL,NULL,&mtop);
 643   bPDB = opt2bSet("-q",NFILE,fnm);
 644   if (bPDB) {
 645     snew(xav,ntopatoms);
 646     snew(ind_fit,ntopatoms);
 647     snew(w_rls,ntopatoms);
 648     for(kkk=0; (kkk<ntopatoms); kkk++) {
 649       w_rls[kkk] = 1;
 650       ind_fit[kkk] = kkk;
 651     }
 652
 653     snew(atoms,1);
 654     *atoms = gmx_mtop_global_atoms(&mtop);
 655
 656     if (atoms->pdbinfo == NULL) {
 657       snew(atoms->pdbinfo,atoms->nr);
 658     }
 659   }
 660
 661   top = gmx_mtop_generate_local_top(&mtop,&ir);
 662
 663   g = NULL;
 664   pbc_null = NULL;
 665   if (ir.ePBC != epbcNONE) {
 666     if (ir.bPeriodicMols)
 667       pbc_null = &pbc;
 668     else
 669       g = mk_graph(fplog,&top->idef,0,mtop.natoms,FALSE,FALSE);
 670   }
 671
 672   if (ftp2bSet(efNDX,NFILE,fnm)) {
 673     rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 674     xvg=xvgropen(opt2fn("-dr",NFILE,fnm),"Individual Restraints","Time (ps)",
 675                  "nm",oenv);
 676     snew(vvindex,isize);
 677     snew(leg,isize);
 678     for(i=0; (i<isize); i++) {
 679       index[i]++;
 680       snew(leg[i],12);
 681       sprintf(leg[i],"index %d",index[i]);
 682     }
 683     xvgr_legend(xvg,isize,(const char**)leg,oenv);
 684   }
 685   else
 686     isize=0;
 687
 688   ir.dr_tau=0.0;
 689   init_disres(fplog,&mtop,&ir,NULL,FALSE,&fcd,NULL,FALSE);
 690
 691   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 692   snew(f,5*natoms);
 693
 694   init_dr_res(&dr,fcd.disres.nres);
 695   if (opt2bSet("-c",NFILE,fnm)) {
 696     clust = cluster_index(fplog,opt2fn("-c",NFILE,fnm));
 697     snew(dr_clust,clust->clust->nr+1);
 698     for(i=0; (i<=clust->clust->nr); i++)
 699       init_dr_res(&dr_clust[i],fcd.disres.nres);
 700   }
 701   else {
 702     out =xvgropen(opt2fn("-ds",NFILE,fnm),
 703                   "Sum of Violations","Time (ps)","nm",oenv);
 704     aver=xvgropen(opt2fn("-da",NFILE,fnm),
 705                   "Average Violation","Time (ps)","nm",oenv);
 706     numv=xvgropen(opt2fn("-dn",NFILE,fnm),
 707                   "# Violations","Time (ps)","#",oenv);
 708     maxxv=xvgropen(opt2fn("-dm",NFILE,fnm),
 709                    "Largest Violation","Time (ps)","nm",oenv);
 710   }
 711
 712   mdatoms = init_mdatoms(fplog,&mtop,ir.efep!=efepNO);
 713   atoms2md(&mtop,&ir,0,NULL,0,mtop.natoms,mdatoms);
 714   update_mdatoms(mdatoms,ir.fepvals->init_lambda);
 715   fr      = mk_forcerec();
 716   fprintf(fplog,"Made forcerec\n");
 717   init_forcerec(fplog,oenv,fr,NULL,&ir,&mtop,cr,box,FALSE,
 718                 NULL,NULL,NULL,NULL,NULL,FALSE,-1);
 719   init_nrnb(&nrnb);
 720   if (ir.ePBC != epbcNONE)
 721     gpbc = gmx_rmpbc_init(&top->idef,ir.ePBC,natoms,box);
 722
 723   j=0;
 724   do {
 725     if (ir.ePBC != epbcNONE) {
 726       if (ir.bPeriodicMols)
 727         set_pbc(&pbc,ir.ePBC,box);
 728       else
 729         gmx_rmpbc(gpbc,natoms,box,x);
 730     }
 731
 732     if (clust) {
 733       if (j > clust->maxframe)
 734         gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file. t = %8f\n",t);
 735       my_clust = clust->inv_clust[j];
 736       range_check(my_clust,0,clust->clust->nr);
 737       check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
 738                  top->idef.iparams,top->idef.functype,
 739                  x,f,fr,pbc_null,g,dr_clust,my_clust,isize,index,vvindex,&fcd);
 740     }
 741     else
 742       check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
 743                  top->idef.iparams,top->idef.functype,
 744                  x,f,fr,pbc_null,g,&dr,0,isize,index,vvindex,&fcd);
 745     if (bPDB) {
 746       reset_x(atoms->nr,ind_fit,atoms->nr,NULL,x,w_rls);
 747       do_fit(atoms->nr,w_rls,x,x);
 748       if (j == 0) {
 749         /* Store the first frame of the trajectory as 'characteristic'
 750          * for colouring with violations.
 751          */
 752         for(kkk=0; (kkk<atoms->nr); kkk++)
 753           copy_rvec(x[kkk],xav[kkk]);
 754       }
 755     }
 756     if (!clust) {
 757       if (isize > 0) {
 758         fprintf(xvg,"%10g",t);
 759         for(i=0; (i<isize); i++)
 760           fprintf(xvg,"  %10g",vvindex[i]);
 761         fprintf(xvg,"\n");
 762       }
 763       fprintf(out,  "%10g  %10g\n",t,dr.sumv);
 764       fprintf(aver, "%10g  %10g\n",t,dr.averv);
 765       fprintf(maxxv,"%10g  %10g\n",t,dr.maxv);
 766       fprintf(numv, "%10g  %10d\n",t,dr.nv);
 767     }
 768     j++;
 769   } while (read_next_x(oenv,status,&t,natoms,x,box));
 770   close_trj(status);
 771   if (ir.ePBC != epbcNONE)
 772     gmx_rmpbc_done(gpbc);
 773
 774   if (clust) {
 775     dump_clust_stats(fplog,fcd.disres.nres,&(top->idef.il[F_DISRES]),
 776                      top->idef.iparams,clust->clust,dr_clust,
 777                      clust->grpname,isize,index);
 778   }
 779   else {
 780     dump_stats(fplog,j,fcd.disres.nres,&(top->idef.il[F_DISRES]),
 781                top->idef.iparams,&dr,isize,index,
 782                bPDB ? atoms : NULL);
 783     if (bPDB) {
 784       write_sto_conf(opt2fn("-q",NFILE,fnm),
 785                      "Coloured by average violation in Angstrom",
 786                      atoms,xav,NULL,ir.ePBC,box);
 787     }
 788     dump_disre_matrix(opt2fn_null("-x",NFILE,fnm),&dr,fcd.disres.nres,
 789                       j,&top->idef,&mtop,max_dr,nlevels,bThird);
 790     ffclose(out);
 791     ffclose(aver);
 792     ffclose(numv);
 793     ffclose(maxxv);
 794     if (isize > 0) {
 795       ffclose(xvg);
 796       do_view(oenv,opt2fn("-dr",NFILE,fnm),"-nxy");
 797     }
 798     do_view(oenv,opt2fn("-dn",NFILE,fnm),"-nxy");
 799     do_view(oenv,opt2fn("-da",NFILE,fnm),"-nxy");
 800     do_view(oenv,opt2fn("-ds",NFILE,fnm),"-nxy");
 801     do_view(oenv,opt2fn("-dm",NFILE,fnm),"-nxy");
 802   }
 803   thanx(stderr);
 804
 805   gmx_finalize_par();
 806
 807   gmx_log_close(fplog);
 808
 809   return 0;
 810 }