src/tools/g_lie.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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include <stdio.h>
  34 #include <stdlib.h>
  35 #include <math.h>
  36 #include <string.h>
  37 #include "statutil.h"
  38 #include "sysstuff.h"
  39 #include "typedefs.h"
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "fatal.h"
  43 #include "vec.h"
  44 #include "copyrite.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "txtdump.h"
  48 #include "enxio.h"
  49 #include "gstat.h"
  50 #include "xvgr.h"
  51
  52 typedef struct {
  53   int nlj,nqq;
  54   int *lj;
  55   int *qq;
  56 } t_liedata;
  57
  58 static t_liedata *analyze_names(int nre,char *names[],char *ligand)
  59 {
  60   int       i;
  61   t_liedata *ld;
  62   char      self[256];
  63
  64   /* Skip until we come to pressure */
  65   for(i=0; (i<F_NRE); i++)
  66     if (strcmp(names[i],interaction_function[F_PRES].longname) == 0)
  67       break;
  68
  69   /* Now real analysis: find components of energies */
  70   sprintf(self,"%s-%s",ligand,ligand);
  71   snew(ld,1);
  72   for( ; (i<nre); i++) {
  73     if ((strstr(names[i],ligand) != NULL) &&
  74         (strstr(names[i],self) == NULL)) {
  75       if (strstr(names[i],"LJ") != NULL) {
  76         ld->nlj++;
  77         srenew(ld->lj,ld->nlj);
  78         ld->lj[ld->nlj-1] = i;
  79       }
  80       else if (strstr(names[i],"Coul") != NULL) {
  81         ld->nqq++;
  82         srenew(ld->qq,ld->nqq);
  83         ld->qq[ld->nqq-1] = i;
  84       }
  85     }
  86   }
  87   printf("Using the following energy terms:\n");
  88   printf("LJ:  ");
  89   for(i=0; (i<ld->nlj); i++)
  90     printf("  %12s",names[ld->lj[i]]);
  91   printf("\nCoul:");
  92   for(i=0; (i<ld->nqq); i++)
  93     printf("  %12s",names[ld->qq[i]]);
  94   printf("\n");
  95
  96   return ld;
  97 }
  98
  99 real calc_lie(t_liedata *ld,t_energy ee[],real lie_lj,real lie_qq,
 100               real fac_lj,real fac_qq)
 101 {
 102   int i;
 103   real lj_tot,qq_tot;
 104
 105   lj_tot = 0;
 106   for(i=0; (i<ld->nlj); i++)
 107     lj_tot += ee[ld->lj[i]].e;
 108   qq_tot = 0;
 109   for(i=0; (i<ld->nqq); i++)
 110     qq_tot += ee[ld->qq[i]].e;
 111
 112   /* And now the great LIE formula: */
 113   return fac_lj*(lj_tot-lie_lj)+fac_qq*(qq_tot-lie_qq);
 114 }
 115
 116 int main(int argc,char *argv[])
 117 {
 118   static char *desc[] = {
 119     "g_lie computes a free energy estimate based on an energy analysis",
 120     "from. One needs an energy file with the following components:",
 121     "Coul (A-B) LJ-SR (A-B) etc."
 122   };
 123   static real lie_lj=0,lie_qq=0,fac_lj=0.181,fac_qq=0.5;
 124   static char *ligand="none";
 125   t_pargs pa[] = {
 126     { "-Elj",  FALSE, etREAL, {&lie_lj},
 127       "Lennard-Jones interaction between ligand and solvent" },
 128     { "-Eqq",  FALSE, etREAL, {&lie_qq},
 129       "Coulomb interaction between ligand and solvent" },
 130     { "-Clj",  FALSE, etREAL, {&fac_lj},
 131       "Factor in the LIE equation for Lennard-Jones component of energy" },
 132     { "-Cqq",  FALSE, etREAL, {&lie_qq},
 133       "Factor in the LIE equation for Coulomb component of energy" },
 134     { "-ligand",  FALSE, etSTR, {&ligand},
 135       "Name of the ligand in the energy file" }
 136   };
 137 #define NPA asize(pa)
 138
 139   FILE      *out;
 140   int       fp,nre,nframes=0,ct=0;
 141   char      **enm=NULL;
 142   bool      bCont;
 143   t_liedata *ld;
 144   t_enxframe *fr;
 145   real      lie;
 146   double    lieaver=0,lieav2=0;
 147
 148   t_filenm fnm[] = {
 149     { efENX, "-f",    "ener",     ffREAD   },
 150     { efXVG, "-o",    "lie",      ffWRITE  }
 151   };
 152 #define NFILE asize(fnm)
 153
 154   CopyRight(stderr,argv[0]);
 155   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 156                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 157
 158   fp = open_enx(ftp2fn(efENX,NFILE,fnm),"r");
 159   do_enxnms(fp,&nre,&enm);
 160
 161   ld = analyze_names(nre,enm,ligand);
 162   snew(fr,1);
 163   out = xvgropen(ftp2fn(efXVG,NFILE,fnm),"LIE free energy estimate",
 164                  "Time (ps)","DGbind (kJ/mol)");
 165   do {
 166     bCont = do_enx(fp,fr);
 167     ct    = check_times(fr->t);
 168     if (ct == 0) {
 169       lie = calc_lie(ld,fr->ener,lie_lj,lie_qq,fac_lj,fac_qq);
 170       lieaver += lie;
 171       lieav2  += lie*lie;
 172       nframes ++;
 173       fprintf(out,"%10g  %10g\n",fr->t,lie);
 174     }
 175   } while (bCont);
 176   close_enx(fp);
 177   fclose(out);
 178   fprintf(stderr,"\n");
 179
 180   if (nframes > 0)
 181     printf("DGbind = %.3f (%.3f)\n",lieaver/nframes,
 182            sqrt(lieav2/nframes-sqr(lieaver/nframes)));
 183
 184   do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
 185
 186   thanx(stderr);
 187
 188   return 0;
 189 }
 190