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