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