src/tools/g_enemat.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_enemat_c = "$Id$";
  30
  31 #include <string.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "fatal.h"
  35 #include "vec.h"
  36 #include "smalloc.h"
  37 #include "enxio.h"
  38 #include "statutil.h"
  39 #include "assert.h"
  40 #include "names.h"
  41 #include "copyrite.h"
  42 #include "macros.h"
  43 #include "xvgr.h"
  44 #include "gstat.h"
  45 #include "physics.h"
  46 #include "matio.h"
  47 #include "strdb.h"
  48
  49 static int search_str2(int nstr,char **str,char *key)
  50 {
  51   int  i,n;
  52
  53   /* Linear search */
  54   n=0;
  55   while( (n<strlen(key)) && ((key[n]<'0') || (key[n]>'9')) )
  56     n++;
  57   for(i=0; (i<nstr); i++)
  58     if (strncasecmp(str[i],key,n)==0)
  59       return i;
  60
  61   return -1;
  62 }
  63
  64 int main(int argc,char *argv[])
  65 {
  66   static char *desc[] = {
  67     "g_enemat extracts an energy matrix from an energy file.",
  68     "With [BB]-groups[bb] a file must be supplied with on each",
  69     "line a group to be used. For these groups a matrices of",
  70     "interaction energies will be calculated. Also the total",
  71     "interaction energy energy per group is calculated.[PAR]",
  72     "An approximation of the free energy is calculated using:",
  73     "E(free) = E0 + kT log( <exp((E-E0)/kT)> ), where '<>'",
  74     "stands for time-average. A file with reference free energies",
  75     "can be supplied to calculate the free energy difference",
  76     "with some reference state. Group names (e.g. residue names",
  77     "in the reference file should correspond to the group names",
  78     "as used in the [BB]-groups[bb] file, but a appended number",
  79     "(e.g. residue number)in the [BB]-groups[bb] will be ignored",
  80     "in the comparison."
  81   };
  82   static bool bSum=FALSE;
  83   static bool bMeanEmtx=TRUE;
  84   static int  skip=0,nlevels=20;
  85   static real cutmax=1e20,cutmin=-1e20,reftemp=300.0;
  86   static bool bCoul=TRUE,bCoulLR=FALSE,bCoul14=FALSE;
  87   static bool bLJ=TRUE,bLJ14=FALSE,bBham=FALSE,bFree=TRUE;
  88   t_pargs pa[] = {
  89     { "-sum",  FALSE, etBOOL, {&bSum},
  90       "Sum the energy terms selected rather than display them all" },
  91     { "-skip", FALSE, etINT,  {&skip},
  92       "Skip number of frames between data points" },
  93     { "-mean", FALSE, etBOOL, {&bMeanEmtx},
  94       "with -groups calculates matrix of mean energies in stead of "
  95       "matrix for each timestep" },
  96     { "-nlevels", FALSE, etINT, {&nlevels},"number of levels for matrix colors"},
  97     { "-max",FALSE, etREAL, {&cutmax},"max value for energies"},
  98     { "-min",FALSE, etREAL, {&cutmin},"min value for energies"},
  99     { "-coul", FALSE, etBOOL, {&bCoul},"calculate Coulomb SR energies"},
 100     { "-coulr", FALSE, etBOOL, {&bCoulLR},"calculate Coulomb LR energies"},
 101     { "-coul14",FALSE, etBOOL, {&bCoul14},"calculate Coulomb 1-4 energies"},
 102     { "-lj", FALSE, etBOOL, {&bLJ},"calculate Lennard-Jones SR energies"},
 103     { "-lj14",FALSE, etBOOL, {&bLJ14},"calculate Lennard-Jones 1-4 energies"},
 104     { "-bham",FALSE, etBOOL, {&bBham},"calculate Buckingham energies"},
 105     { "-free",FALSE, etBOOL, {&bFree},"calculate free energy"},
 106     { "-temp",FALSE, etREAL, {&reftemp},
 107       "reference temperature for free energy calculation"}
 108   };
 109   /* We will define egSP more energy-groups:
 110      egTotal (total energy) */
 111 #define egTotal egNR
 112 #define egSP 1
 113   bool       egrp_use[egNR+egSP];
 114   int        in;
 115   FILE       *out;
 116   int        timecheck=0;
 117   t_energy   *ee;
 118   t_drblock  dr;
 119   int        teller=0,nre,step;
 120   real       t,sum;
 121   bool       bCont,bRef;
 122   bool       bCutmax,bCutmin;
 123   real       **eneset,*time=NULL;
 124   int        *set,i,j,k,prevk,m=0,n,nset,nenergy,ndr;
 125   char       **enm,**groups;
 126   char       groupname[255],fn[255];
 127   int        ngroups;
 128   t_rgb      rlo,rhi,rmid;
 129   real       emax,emid,emin;
 130   real       ***emat,**etot,*groupnr;
 131   double     beta,expE,**e,*eaver,*efree=NULL,edum;
 132   char       label[234];
 133   char       **ereflines,**erefres=NULL;
 134   real       *eref=NULL,*edif=NULL;
 135   int        neref=0;
 136
 137   t_filenm   fnm[] = {
 138     { efENX, "-f", NULL, ffOPTRD },
 139     { efDAT, "-groups", "groups.dat", ffREAD },
 140     { efDAT, "-eref",   "eref.dat", ffOPTRD },
 141     { efXPM, "-emat",   "emat", ffWRITE },
 142     { efXVG, "-etot",   "energy", ffWRITE }
 143   };
 144 #define NFILE asize(fnm)
 145
 146   CopyRight(stderr,argv[0]);
 147   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 148                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 149
 150   egrp_use[egCOUL]=bCoul;
 151   egrp_use[egLJ]=bLJ;
 152   egrp_use[egBHAM]=bBham;
 153   egrp_use[egLR]=bCoulLR;
 154   egrp_use[egCOUL14]=bCoul14;
 155   egrp_use[egLJ14]=bLJ14;
 156   egrp_use[egTotal]=TRUE;
 157
 158   bRef=opt2bSet("-eref",NFILE,fnm);
 159   in=open_enx(ftp2fn(efENE,NFILE,fnm),"r");
 160   do_enxnms(in,&nre,&enm);
 161
 162   if (nre == 0)
 163     fatal_error(0,"No energies!\n");
 164
 165   bCutmax=opt2parg_bSet("-max",asize(pa),pa);
 166   bCutmin=opt2parg_bSet("-min",asize(pa),pa);
 167
 168   snew(ee,nre);
 169   nenergy = 0;
 170
 171   /* Read groupnames from input file and construct selection of
 172      energy groups from it*/
 173
 174   fprintf(stderr,"Will read groupnames from inputfile\n");
 175   ngroups = get_lines(opt2fn("-groups",NFILE,fnm), &groups);
 176   fprintf(stderr,"Read %d groups\n",ngroups);
 177   snew(set,sqr(ngroups)*egNR/2);
 178   n=0;
 179   prevk=0;
 180   for (i=0; (i<ngroups); i++) {
 181     fprintf(stderr,"\rgroup %d",i);
 182     for (j=i; (j<ngroups); j++)
 183       for(m=0; (m<egNR); m++)
 184         if (egrp_use[m]) {
 185           sprintf(groupname,"%s:%s-%s",egrp_nm[m],groups[i],groups[j]);
 186 #ifdef DEBUG
 187           fprintf(stderr,"\r%-15s %5d",groupname,n);
 188 #endif
 189           for(k=prevk; (k<prevk+nre); k++)
 190             if (strcmp(enm[k%nre],groupname) == 0) {
 191               set[n++] = k;
 192               break;
 193             }
 194           if (k==prevk+nre)
 195             fprintf(stderr,"WARNING! could not find group %s (%d,%d)"
 196                     "in energy file\n",groupname,i,j);
 197           else
 198             prevk = k;
 199         }
 200   }
 201   fprintf(stderr,"\n");
 202   nset=n;
 203   snew(eneset,nset+1);
 204   fprintf(stderr,"Will select half-matrix of energies with %d elements\n",n);
 205
 206   /* Start reading energy frames */
 207   step     = 0;
 208   t        = 0;
 209   do {
 210     do {
 211       bCont=do_enx(in,&t,&step,&nre,ee,&ndr,&dr);
 212       if (bCont) {
 213         timecheck=check_times(t);
 214
 215         /* It is necessary for statistics to start counting from 1 */
 216         step += 1;
 217       }
 218
 219     } while (bCont && (timecheck < 0));
 220
 221     if (timecheck == 0) {
 222 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
 223
 224       if (bCont) {
 225         fprintf(stderr,"\rRead frame: %d, Time: %.3f",teller,t);
 226
 227         if ((nenergy % 1000) == 0) {
 228           srenew(time,nenergy+1000);
 229           for(i=0; (i<=nset); i++)
 230             srenew(eneset[i],nenergy+1000);
 231         }
 232         time[nenergy] = t;
 233         sum=0;
 234         for(i=0; (i<nset); i++) {
 235           eneset[i][nenergy] = ee[set[i]].e;
 236           sum+=ee[set[i]].e;
 237         }
 238         if (bSum)
 239           eneset[nset][nenergy] = sum;
 240         nenergy++;
 241       }
 242       teller++;
 243     }
 244   } while (bCont && (timecheck == 0));
 245
 246   fprintf(stderr,"\n");
 247
 248   fprintf(stderr,"Will build energy half-matrix of %d groups, %d elements, "
 249           "over %d frames\n",ngroups,nset,nenergy);
 250
 251   snew(emat,egNR+egSP);
 252   for(j=0; (j<egNR+egSP); j++)
 253     if (egrp_use[m]) {
 254       snew(emat[j],ngroups);
 255       for (i=0; (i<ngroups); i++)
 256         snew(emat[j][i],ngroups);
 257     }
 258   snew(groupnr,ngroups);
 259   for (i=0; (i<ngroups); i++)
 260     groupnr[i] = i+1;
 261   rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 262   rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 263   rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 264   if (bMeanEmtx) {
 265     snew(e,ngroups);
 266     for (i=0; (i<ngroups); i++)
 267       snew(e[i],nenergy);
 268     n = 0;
 269     for (i=0; (i<ngroups); i++) {
 270       for (j=i; (j<ngroups); j++) {
 271         for (m=0; (m<egNR); m++) {
 272           if (egrp_use[m]) {
 273             for (k=0; (k<nenergy); k++) {
 274               emat[m][i][j] += eneset[n][k];
 275               e[i][k] += eneset[n][k];/* *0.5; */
 276               e[j][k] += eneset[n][k];/* *0.5; */
 277             }
 278             n++;
 279             emat[egTotal][i][j]+=emat[m][i][j];
 280             emat[m][i][j]/=nenergy;
 281             emat[m][j][i]=emat[m][i][j];
 282           }
 283         }
 284         emat[egTotal][i][j]/=nenergy;
 285         emat[egTotal][j][i]=emat[egTotal][i][j];
 286       }
 287     }
 288     if (bFree) {
 289       if (bRef) {
 290         fprintf(stderr,"Will read reference energies from inputfile\n");
 291         neref = get_lines(opt2fn("-eref",NFILE,fnm), &ereflines);
 292         fprintf(stderr,"Read %d reference energies\n",neref);
 293         snew(eref, neref);
 294         snew(erefres, neref);
 295         for(i=0; (i<neref); i++) {
 296           snew(erefres[i], 5);
 297           sscanf(ereflines[i],"%s %lf",erefres[i],&edum);
 298           eref[i] = edum;
 299         }
 300       }
 301       snew(eaver,ngroups);
 302       for (i=0; (i<ngroups); i++) {
 303         for (k=0; (k<nenergy); k++)
 304           eaver[i] += e[i][k];
 305         eaver[i] /= nenergy;
 306       }
 307       beta = 1.0/(BOLTZ*reftemp);
 308       snew(efree,ngroups);
 309       snew(edif,ngroups);
 310       for (i=0; (i<ngroups); i++) {
 311         expE=0;
 312         for (k=0; (k<nenergy); k++) {
 313           expE += exp(beta*(e[i][k]-eaver[i]));
 314         }
 315         efree[i] = log(expE/nenergy)/beta + eaver[i];
 316         if (bRef) {
 317           n = search_str2(neref,erefres,groups[i]);
 318           if (n != -1) {
 319             edif[i] = efree[i]-eref[n];
 320           } else {
 321             edif[i] = efree[i];
 322             fprintf(stderr,"WARNING: group %s not found "
 323                     "in reference energies.\n",groups[i]);
 324           }
 325         } else
 326           edif[i]=0;
 327       }
 328     }
 329
 330     emid = 0.0;/*(emin+emax)*0.5;*/
 331     for(m=0; (m<egNR); m++)
 332       egrp_nm[m]=egrp_nm[m];
 333     egrp_nm[egTotal]="total";
 334     for (m=0; (m<egNR+egSP); m++)
 335       if (egrp_use[m]) {
 336         emin=1e10;
 337         emax=-1e10;
 338         for (i=0; (i<ngroups); i++) {
 339           for (j=i; (j<ngroups); j++) {
 340             if (emat[m][i][j] > emax)
 341               emax=emat[m][i][j];
 342             else if (emat[m][i][j] < emin)
 343               emin=emat[m][i][j];
 344           }
 345         }
 346         if (emax==emin)
 347           fprintf(stderr,"Matrix of %s energy is uniform at %f "
 348                   "(will not produce output).\n",egrp_nm[m],emax);
 349         else {
 350           fprintf(stderr,"Matrix of %s energy ranges from %f to %f\n",
 351                   egrp_nm[m],emin,emax);
 352           if ((bCutmax) || (emax>cutmax)) emax=cutmax;
 353           if ((bCutmin) || (emin<cutmin)) emin=cutmin;
 354           if ((emax==cutmax) || (emin==cutmin))
 355             fprintf(stderr,"Energy range adjusted: %f to %f\n",emin,emax);
 356
 357           sprintf(fn,"%s%s",egrp_nm[m],ftp2fn(efXPM,NFILE,fnm));
 358           sprintf(label,"%s Interaction Energies",egrp_nm[m]);
 359           out=ffopen(fn,"w");
 360           if (emin>=emid)
 361             write_xpm(out,label,"Energy (kJ/mol)",
 362                       "Residue Index","Residue Index",
 363                       ngroups,ngroups,groupnr,groupnr,emat[m],
 364                       emid,emax,rmid,rhi,&nlevels);
 365           else if (emax<=emid)
 366             write_xpm(out,label,"Energy (kJ/mol)",
 367                       "Residue Index","Residue Index",
 368                       ngroups,ngroups,groupnr,groupnr,emat[m],
 369                       emin,emid,rlo,rmid,&nlevels);
 370           else
 371             write_xpm3(out,label,"Energy (kJ/mol)",
 372                        "Residue Index","Residue Index",
 373                        ngroups,ngroups,groupnr,groupnr,emat[m],
 374                        emin,emid,emax,rlo,rmid,rhi,&nlevels);
 375           ffclose(out);
 376         }
 377       }
 378     snew(etot,egNR+egSP);
 379     for (m=0; (m<egNR+egSP); m++) {
 380       snew(etot[m],ngroups);
 381       for (i=0; (i<ngroups); i++) {
 382         for (j=0; (j<ngroups); j++)
 383           etot[m][i]+=emat[m][i][j];
 384       }
 385     }
 386
 387     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Mean Energy","Residue","kJ/mol");
 388     xvgr_legend(out,0,NULL);
 389     j=0;
 390     for (m=0; (m<egNR+egSP); m++)
 391       if (egrp_use[m])
 392         fprintf(out,"@ legend string %d \"%s\"\n",j++,egrp_nm[m]);
 393     if (bFree)
 394       fprintf(out,"@ legend string %d \"%s\"\n",j++,"Free");
 395     if (bFree)
 396       fprintf(out,"@ legend string %d \"%s\"\n",j++,"Diff");
 397     fprintf(out,"@TYPE nxy\n");
 398     fprintf(out,"#%3s","grp");
 399     for (m=0; (m<egNR+egSP); m++)
 400       if (egrp_use[m])
 401         fprintf(out," %9s",egrp_nm[m]);
 402     if (bFree)
 403       fprintf(out," %9s","Free");
 404     if (bFree)
 405       fprintf(out," %9s","Diff");
 406     fprintf(out,"\n");
 407     for (i=0; (i<ngroups); i++) {
 408       fprintf(out,"%3.0f",groupnr[i]);
 409       for (m=0; (m<egNR+egSP); m++)
 410         if (egrp_use[m])
 411           fprintf(out," %9.5g",etot[m][i]);
 412       if (bFree)
 413         fprintf(out," %9.5g",efree[i]);
 414       if (bRef)
 415         fprintf(out," %9.5g",edif[i]);
 416       fprintf(out,"\n");
 417     }
 418     fclose(out);
 419   } else {
 420     fprintf(stderr,"While typing at your keyboard, suddenly...\n"
 421             "...nothing happens.\nWARNING: Not Implemented Yet\n");
 422 /*
 423     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 424     n=0;
 425     emin=emax=0.0;
 426     for (k=0; (k<nenergy); k++) {
 427       for (i=0; (i<ngroups); i++)
 428         for (j=i+1; (j<ngroups); j++)
 429           emat[i][j]=eneset[n][k];
 430       sprintf(label,"t=%.0f ps",time[k]);
 431       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 432       n++;
 433     }
 434     fclose(out);
 435 */
 436   }
 437   close_enx(in);
 438
 439   thanx(stdout);
 440
 441   return 0;
 442 }