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     "E(free) = E0 + kT log( <exp((E-E0)/kT)> ), where '<>'",
  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     { "-coul", FALSE, etBOOL, {&bCoulSR},"extract Coulomb SR energies"},
 126     { "-coulr", FALSE, etBOOL, {&bCoulLR},"extract Coulomb LR energies"},
 127     { "-coul14",FALSE, etBOOL, {&bCoul14},"extract Coulomb 1-4 energies"},
 128     { "-lj", FALSE, etBOOL, {&bLJSR},"extract Lennard-Jones SR energies"},
 129     { "-lj", 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     for(m=0; (m<egNR); m++)
 356       egrp_nm[m]=egrp_nm[m];
 357     egrp_nm[egTotal]="total";
 358     for (m=0; (m<egNR+egSP); m++)
 359       if (egrp_use[m]) {
 360         emin=1e10;
 361         emax=-1e10;
 362         for (i=0; (i<ngroups); i++) {
 363           for (j=i; (j<ngroups); j++) {
 364             if (emat[m][i][j] > emax)
 365               emax=emat[m][i][j];
 366             else if (emat[m][i][j] < emin)
 367               emin=emat[m][i][j];
 368           }
 369         }
 370         if (emax==emin)
 371           fprintf(stderr,"Matrix of %s energy is uniform at %f "
 372                   "(will not produce output).\n",egrp_nm[m],emax);
 373         else {
 374           fprintf(stderr,"Matrix of %s energy ranges from %f to %f\n",
 375                   egrp_nm[m],emin,emax);
 376           if ((bCutmax) || (emax>cutmax)) emax=cutmax;
 377           if ((bCutmin) || (emin<cutmin)) emin=cutmin;
 378           if ((emax==cutmax) || (emin==cutmin))
 379             fprintf(stderr,"Energy range adjusted: %f to %f\n",emin,emax);
 380
 381           sprintf(fn,"%s%s",egrp_nm[m],ftp2fn(efXPM,NFILE,fnm));
 382           sprintf(label,"%s Interaction Energies",egrp_nm[m]);
 383           out=ffopen(fn,"w");
 384           if (emin>=emid)
 385             write_xpm(out,0,label,"Energy (kJ/mol)",
 386                       "Residue Index","Residue Index",
 387                       ngroups,ngroups,groupnr,groupnr,emat[m],
 388                       emid,emax,rmid,rhi,&nlevels);
 389           else if (emax<=emid)
 390             write_xpm(out,0,label,"Energy (kJ/mol)",
 391                       "Residue Index","Residue Index",
 392                       ngroups,ngroups,groupnr,groupnr,emat[m],
 393                       emin,emid,rlo,rmid,&nlevels);
 394           else
 395             write_xpm3(out,0,label,"Energy (kJ/mol)",
 396                        "Residue Index","Residue Index",
 397                        ngroups,ngroups,groupnr,groupnr,emat[m],
 398                        emin,emid,emax,rlo,rmid,rhi,&nlevels);
 399           ffclose(out);
 400         }
 401       }
 402     snew(etot,egNR+egSP);
 403     for (m=0; (m<egNR+egSP); m++) {
 404       snew(etot[m],ngroups);
 405       for (i=0; (i<ngroups); i++) {
 406         for (j=0; (j<ngroups); j++)
 407           etot[m][i]+=emat[m][i][j];
 408       }
 409     }
 410
 411     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Mean Energy","Residue","kJ/mol",
 412                  oenv);
 413     xvgr_legend(out,0,NULL, oenv);
 414     j=0;
 415     for (m=0; (m<egNR+egSP); m++)
 416       if (egrp_use[m])
 417         fprintf(out,"@ legend string %d \"%s\"\n",j++,egrp_nm[m]);
 418     if (bFree)
 419       fprintf(out,"@ legend string %d \"%s\"\n",j++,"Free");
 420     if (bFree)
 421       fprintf(out,"@ legend string %d \"%s\"\n",j++,"Diff");
 422     fprintf(out,"@TYPE xy\n");
 423     fprintf(out,"#%3s","grp");
 424     for (m=0; (m<egNR+egSP); m++)
 425       if (egrp_use[m])
 426         fprintf(out," %9s",egrp_nm[m]);
 427     if (bFree)
 428       fprintf(out," %9s","Free");
 429     if (bFree)
 430       fprintf(out," %9s","Diff");
 431     fprintf(out,"\n");
 432     for (i=0; (i<ngroups); i++) {
 433       fprintf(out,"%3.0f",groupnr[i]);
 434       for (m=0; (m<egNR+egSP); m++)
 435         if (egrp_use[m])
 436           fprintf(out," %9.5g",etot[m][i]);
 437       if (bFree)
 438         fprintf(out," %9.5g",efree[i]);
 439       if (bRef)
 440         fprintf(out," %9.5g",edif[i]);
 441       fprintf(out,"\n");
 442     }
 443     ffclose(out);
 444   } else {
 445     fprintf(stderr,"While typing at your keyboard, suddenly...\n"
 446             "...nothing happens.\nWARNING: Not Implemented Yet\n");
 447 /*
 448     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 449     n=0;
 450     emin=emax=0.0;
 451     for (k=0; (k<nenergy); k++) {
 452       for (i=0; (i<ngroups); i++)
 453         for (j=i+1; (j<ngroups); j++)
 454           emat[i][j]=eneset[n][k];
 455       sprintf(label,"t=%.0f ps",time[k]);
 456       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 457       n++;
 458     }
 459     ffclose(out);
 460 */
 461   }
 462   close_enx(in);
 463
 464   thanx(stderr);
 465
 466   return 0;
 467 }