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