src/tools/g_disre.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_g_disre_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "typedefs.h"
  34 #include "macros.h"
  35 #include "copyrite.h"
  36 #include "mshift.h"
  37 #include "xvgr.h"
  38 #include "smalloc.h"
  39 #include "nrnb.h"
  40 #include "disre.h"
  41 #include "statutil.h"
  42 #include "force.h"
  43 #include "gstat.h"
  44 #include "main.h"
  45 #include "genhydro.h"
  46 #include "pdbio.h"
  47 #include "rdgroup.h"
  48 #include "mdatoms.h"
  49 #include "nsb.h"
  50 #include "tpxio.h"
  51 #include "init.h"
  52
  53 typedef struct {
  54   int n;
  55   real v;
  56 } t_toppop;
  57
  58 t_toppop *top=NULL;
  59 int      ntop=0;
  60
  61 static void init5(int n)
  62 {
  63   ntop=n;
  64   snew(top,ntop);
  65 }
  66
  67 static void reset5(void)
  68 {
  69   int i;
  70
  71   for(i=0; (i<ntop); i++) {
  72     top[i].n=-1;
  73     top[i].v= 0;
  74   }
  75 }
  76
  77 int tpcomp(const void *a,const void *b)
  78 {
  79   t_toppop *tpa;
  80   t_toppop *tpb;
  81
  82   tpa=(t_toppop *)a;
  83   tpb=(t_toppop *)b;
  84
  85   return  (1e7*(tpb->v-tpa->v));
  86 }
  87
  88 static void add5(int ndr,real viol)
  89 {
  90   int i,mini;
  91
  92   mini=0;
  93   for(i=1; (i<ntop); i++)
  94     if (top[i].v < top[mini].v)
  95       mini=i;
  96   if (viol > top[mini].v) {
  97     top[mini].v=viol;
  98     top[mini].n=ndr;
  99   }
 100 }
 101
 102 static void print5(FILE *fp)
 103 {
 104   int i;
 105
 106   qsort(top,ntop,sizeof(top[0]),tpcomp);
 107   fprintf(fp,"Index:");
 108   for(i=0; (i<ntop); i++)
 109     fprintf(fp," %6d",top[i].n);
 110   fprintf(fp,"\nViol: ");
 111   for(i=0; (i<ntop); i++)
 112     fprintf(fp," %6.3f",top[i].v);
 113   fprintf(fp,"\n");
 114 }
 115
 116 void check_viol(FILE *log,
 117                 t_ilist *bonds,t_iparams forceparams[],
 118                 t_functype functype[],
 119                 rvec x[],rvec f[],
 120                 t_forcerec *fr,matrix box,t_graph *g,
 121                 real *sumv,real *averv,
 122                 real *maxv,int *nv,
 123                 int isize,atom_id index[],real vvindex[])
 124 {
 125   t_iatom *forceatoms;
 126   int     i,j,nat,n,type,ftype,nviol,ndr;
 127   real    mviol,tviol,viol,lam,dvdl;
 128   static  bool bFirst=TRUE;
 129
 130   lam  =0;
 131   dvdl =0;
 132   tviol=0;
 133   nviol=0;
 134   mviol=0;
 135   ndr=0;
 136   reset5();
 137   forceatoms=bonds->iatoms;
 138   for(j=0; (j<isize); j++) {
 139     vvindex[j]=0;
 140   }
 141   for(i=0; (i<bonds->nr); ) {
 142     type=forceatoms[i];
 143     ftype=functype[type];
 144     nat=interaction_function[ftype].nratoms;
 145     n=0;
 146     do
 147       n += 1+nat;
 148     while ((i+n < bonds->nr) &&
 149            (forceparams[forceatoms[i+n]].disres.index ==
 150             forceparams[forceatoms[i]].disres.index));
 151
 152     viol=interaction_function[ftype].ifunc(n,&forceatoms[i],
 153                                            forceparams,
 154                                            x,f,fr,g,box,lam,&dvdl,
 155                                            NULL,0,NULL,NULL);
 156     if (viol > 0) {
 157       nviol++;
 158       add5(forceparams[type].disres.index,viol);
 159       if (viol > mviol)
 160         mviol=viol;
 161       tviol+=viol;
 162       for(j=0; (j<isize); j++) {
 163         if (index[j] == forceparams[type].disres.index)
 164           vvindex[j]=viol;
 165         }
 166     }
 167     ndr ++;
 168     i   += n;
 169   }
 170   *nv   = nviol;
 171   *maxv = mviol;
 172   *sumv = tviol;
 173   *averv= tviol/ndr;
 174
 175   if (bFirst) {
 176     fprintf(stderr,"\nThere are %d restraints and %d pairs\n",ndr,bonds->nr/3);
 177     bFirst = FALSE;
 178   }
 179
 180   print5(log);
 181 }
 182
 183 void patch_viol(t_ilist *bonds,t_iparams forceparams[],
 184                 t_functype functype[])
 185 {
 186   t_iatom *forceatoms;
 187   int     i,j,type,ftype,nat;
 188
 189   forceatoms=bonds->iatoms;
 190   for(i=j=0; (i<bonds->nr); ) {
 191     type=forceatoms[i++];
 192     ftype=functype[type];
 193     if (ftype == F_DISRES)
 194       forceparams[ftype].disres.index=j++;
 195     nat=interaction_function[ftype].nratoms;
 196     i+=nat;
 197   }
 198 }
 199
 200 int main (int argc,char *argv[])
 201 {
 202   static char *desc[] = {
 203     "g_disre computes violations of distance restraints. If necessary",
 204     "all protons can be added to a protein molecule. The program allways",
 205     "computes the instantaneous violations rather than time-averaged,",
 206     "because this analysis is done from a trajectory file afterwards",
 207     "it does not make sense to use time averaging.[PAR]",
 208     "An index file may be used to select specific restraints for",
 209     "printing."
 210   };
 211   static bool bProt=FALSE;
 212   static int  ntop = 6;
 213   t_pargs pa[] = {
 214     { "-prot", FALSE, etBOOL, {&bProt},
 215       "Protonate protein every step. This currently does not add terminal hydrogens, and therefore works only when the termini are capped." },
 216     { "-ntop", FALSE, etINT,  {&ntop},
 217       "Number of large violations that are stored in the log file every step" }
 218   };
 219
 220   FILE        *out,*aver,*numv,*maxxv,*xvg=NULL;
 221   t_tpxheader header;
 222   t_inputrec  ir;
 223   t_topology  top;
 224   rvec        *xtop;
 225   t_atoms     *atoms=NULL;
 226   t_forcerec  *fr;
 227   t_nrnb      nrnb;
 228   t_nsborder  *nsb;
 229   t_commrec   *cr;
 230   t_graph     *g;
 231   int         status,ntopatoms,natoms,i,j,nv,step;
 232   real        t,sumv,averv,maxv,lambda;
 233   rvec        *x,*f;
 234   matrix      box;
 235   int         isize;
 236   atom_id     *index=NULL;
 237   char        *grpname;
 238   char        **leg;
 239   real        *vvindex=NULL;
 240   t_mdatoms   *mdatoms;
 241
 242   t_filenm fnm[] = {
 243     { efTPX, NULL, NULL, ffREAD },
 244     { efTRX, "-f", NULL, ffREAD },
 245     { efXVG, "-ds", "drsum",  ffWRITE },
 246     { efXVG, "-da", "draver", ffWRITE },
 247     { efXVG, "-dn", "drnum",  ffWRITE },
 248     { efXVG, "-dm", "drmax",  ffWRITE },
 249     { efXVG, "-dr", "restr",  ffWRITE },
 250     { efLOG, "-l", "disres", ffWRITE },
 251     { efNDX, NULL, "viol",   ffOPTRD }
 252   };
 253 #define NFILE asize(fnm)
 254
 255   CopyRight(stderr,argv[0]);
 256   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,TRUE,
 257                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 258   init5(ntop);
 259
 260   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header);
 261   snew(xtop,header.natoms);
 262   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,
 263            box,&ntopatoms,xtop,NULL,NULL,&top);
 264
 265   check_nprocs_top(ftp2fn(efTPX,NFILE,fnm),&top,1);
 266
 267   g   = mk_graph(&top.idef,top.atoms.nr,0);
 268   cr  = init_par(&argc,&argv);
 269   open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 270
 271   if (ftp2bSet(efNDX,NFILE,fnm)) {
 272     rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 273     xvg=xvgropen(opt2fn("-dr",NFILE,fnm),"Inidividual Restraints","Time (ps)",
 274                  "nm");
 275     snew(vvindex,isize);
 276     snew(leg,isize);
 277     for(i=0; (i<isize); i++) {
 278       index[i]++;
 279       snew(leg[i],12);
 280       sprintf(leg[i],"index %u",index[i]);
 281     }
 282     xvgr_legend(xvg,isize,leg);
 283   }
 284   else
 285     isize=0;
 286
 287   ir.dr_tau=0.0;
 288   init_disres(stdlog,top.idef.il[F_DISRES].nr,&ir);
 289
 290   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 291   snew(f,5*natoms);
 292
 293   out =xvgropen(opt2fn("-ds",NFILE,fnm),"Sum of Violations","Time (ps)","nm");
 294   aver=xvgropen(opt2fn("-da",NFILE,fnm),"Average Violation","Time (ps)","nm");
 295   numv=xvgropen(opt2fn("-dn",NFILE,fnm),"# Violations","Time (ps)","#");
 296   maxxv=xvgropen(opt2fn("-dm",NFILE,fnm),"Largest Violation","Time (ps)","nm");
 297
 298   snew(nsb,1);
 299   snew(atoms,1);
 300   atoms->nr=top.atoms.nr;
 301   atoms->nres=top.atoms.nres;
 302   snew(atoms->atomname,atoms->nr);
 303   snew(atoms->resname,atoms->nres);
 304   snew(atoms->atom,atoms->nr);
 305   memcpy(atoms->atom,top.atoms.atom,atoms->nr*sizeof(atoms->atom[0]));
 306   memcpy(atoms->atomname,top.atoms.atomname,
 307          atoms->nr*sizeof(atoms->atomname[0]));
 308   memcpy(atoms->resname,top.atoms.resname,
 309          atoms->nres*sizeof(atoms->resname[0]));
 310   mdatoms = atoms2md(&top.atoms,ir.opts.nFreeze,FALSE,FALSE,FALSE);
 311   fr      = mk_forcerec();
 312   fprintf(stdlog,"Made forcerec...\n");
 313   calc_nsb(&(top.blocks[ebCGS]),1,nsb,0);
 314   init_forcerec(stdlog,fr,&ir,&(top.blocks[ebMOLS]),cr,
 315                 &(top.blocks[ebCGS]),&(top.idef),mdatoms,nsb,box,FALSE);
 316   init_nrnb(&nrnb);
 317   j=0;
 318   do {
 319     rm_pbc(&top.idef,natoms,box,x,x);
 320
 321     if (bProt) {
 322       protonate(&atoms,&x);
 323     }
 324
 325     check_viol(stdlog,
 326                &(top.idef.il[F_DISRES]),
 327                top.idef.iparams,top.idef.functype,
 328                x,f,fr,box,g,&sumv,&averv,&maxv,&nv,
 329                isize,index,vvindex);
 330     if (isize > 0) {
 331       fprintf(xvg,"%10g",t);
 332       for(i=0; (i<isize); i++)
 333         fprintf(xvg,"  %10g",vvindex[i]);
 334       fprintf(xvg,"\n");
 335     }
 336     fprintf(out,  "%10g  %10g\n",t,sumv);
 337     fprintf(aver, "%10g  %10g\n",t,averv);
 338     fprintf(maxxv,"%10g  %10g\n",t,maxv);
 339     fprintf(numv, "%10g  %10d\n",t,nv);
 340
 341     j++;
 342   } while (read_next_x(status,&t,natoms,x,box));
 343
 344   close_trj(status);
 345   fclose(out);
 346   fclose(aver);
 347   fclose(numv);
 348   fclose(maxxv);
 349   if (isize > 0) {
 350     fclose(xvg);
 351     xvgr_file(opt2fn("-dr",NFILE,fnm),"-nxy");
 352   }
 353   xvgr_file(opt2fn("-dn",NFILE,fnm),NULL);
 354   xvgr_file(opt2fn("-da",NFILE,fnm),NULL);
 355   xvgr_file(opt2fn("-ds",NFILE,fnm),NULL);
 356   xvgr_file(opt2fn("-dm",NFILE,fnm),NULL);
 357
 358   thanx(stdout);
 359
 360 #ifdef USE_MPI
 361   MPI_Finalize();
 362 #endif
 363
 364   return 0;
 365 }
 366
 367