src/tools/gmx_lie.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,2013, 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
  42 #include <stdio.h>
  43 #include <stdlib.h>
  44 #include <math.h>
  45 #include <string.h>
  46
  47 #include "statutil.h"
  48 #include "sysstuff.h"
  49 #include "typedefs.h"
  50 #include "smalloc.h"
  51 #include "macros.h"
  52 #include "gmx_fatal.h"
  53 #include "vec.h"
  54 #include "copyrite.h"
  55 #include "futil.h"
  56 #include "statutil.h"
  57 #include "txtdump.h"
  58 #include "enxio.h"
  59 #include "gstat.h"
  60 #include "xvgr.h"
  61 #include "gmx_ana.h"
  62
  63
  64 typedef struct {
  65   int nlj,nqq;
  66   int *lj;
  67   int *qq;
  68 } t_liedata;
  69
  70 static t_liedata *analyze_names(int nre,gmx_enxnm_t *names,const char *ligand)
  71 {
  72   int       i;
  73   t_liedata *ld;
  74   char      self[256];
  75
  76   /* Skip until we come to pressure */
  77   for(i=0; (i<F_NRE); i++)
  78     if (strcmp(names[i].name,interaction_function[F_PRES].longname) == 0)
  79       break;
  80
  81   /* Now real analysis: find components of energies */
  82   sprintf(self,"%s-%s",ligand,ligand);
  83   snew(ld,1);
  84   for( ; (i<nre); i++) {
  85     if ((strstr(names[i].name,ligand) != NULL) &&
  86         (strstr(names[i].name,self) == NULL)) {
  87       if (strstr(names[i].name,"LJ") != NULL) {
  88         ld->nlj++;
  89         srenew(ld->lj,ld->nlj);
  90         ld->lj[ld->nlj-1] = i;
  91       }
  92       else if (strstr(names[i].name,"Coul") != NULL) {
  93         ld->nqq++;
  94         srenew(ld->qq,ld->nqq);
  95         ld->qq[ld->nqq-1] = i;
  96       }
  97     }
  98   }
  99   printf("Using the following energy terms:\n");
 100   printf("LJ:  ");
 101   for(i=0; (i<ld->nlj); i++)
 102     printf("  %12s",names[ld->lj[i]].name);
 103   printf("\nCoul:");
 104   for(i=0; (i<ld->nqq); i++)
 105     printf("  %12s",names[ld->qq[i]].name);
 106   printf("\n");
 107
 108   return ld;
 109 }
 110
 111 real calc_lie(t_liedata *ld,t_energy ee[],real lie_lj,real lie_qq,
 112               real fac_lj,real fac_qq)
 113 {
 114   int i;
 115   real lj_tot,qq_tot;
 116
 117   lj_tot = 0;
 118   for(i=0; (i<ld->nlj); i++)
 119     lj_tot += ee[ld->lj[i]].e;
 120   qq_tot = 0;
 121   for(i=0; (i<ld->nqq); i++)
 122     qq_tot += ee[ld->qq[i]].e;
 123
 124   /* And now the great LIE formula: */
 125   return fac_lj*(lj_tot-lie_lj)+fac_qq*(qq_tot-lie_qq);
 126 }
 127
 128 int gmx_lie(int argc,char *argv[])
 129 {
 130   const char *desc[] = {
 131     "[TT]g_lie[tt] computes a free energy estimate based on an energy analysis",
 132     "from. One needs an energy file with the following components:",
 133     "Coul (A-B) LJ-SR (A-B) etc."
 134   };
 135   static real lie_lj=0,lie_qq=0,fac_lj=0.181,fac_qq=0.5;
 136   static const char *ligand="none";
 137   t_pargs pa[] = {
 138     { "-Elj",  FALSE, etREAL, {&lie_lj},
 139       "Lennard-Jones interaction between ligand and solvent" },
 140     { "-Eqq",  FALSE, etREAL, {&lie_qq},
 141       "Coulomb interaction between ligand and solvent" },
 142     { "-Clj",  FALSE, etREAL, {&fac_lj},
 143       "Factor in the LIE equation for Lennard-Jones component of energy" },
 144     { "-Cqq",  FALSE, etREAL, {&fac_qq},
 145       "Factor in the LIE equation for Coulomb component of energy" },
 146     { "-ligand",  FALSE, etSTR, {&ligand},
 147       "Name of the ligand in the energy file" }
 148   };
 149 #define NPA asize(pa)
 150
 151   FILE      *out;
 152   int       nre,nframes=0,ct=0;
 153   ener_file_t fp;
 154   gmx_bool      bCont;
 155   t_liedata *ld;
 156   gmx_enxnm_t *enm=NULL;
 157   t_enxframe *fr;
 158   real      lie;
 159   double    lieaver=0,lieav2=0;
 160   output_env_t oenv;
 161
 162   t_filenm fnm[] = {
 163     { efEDR, "-f",    "ener",     ffREAD   },
 164     { efXVG, "-o",    "lie",      ffWRITE  }
 165   };
 166 #define NFILE asize(fnm)
 167
 168   CopyRight(stderr,argv[0]);
 169   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 170                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 171
 172   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
 173   do_enxnms(fp,&nre,&enm);
 174
 175   ld = analyze_names(nre,enm,ligand);
 176   snew(fr,1);
 177   out = xvgropen(ftp2fn(efXVG,NFILE,fnm),"LIE free energy estimate",
 178                  "Time (ps)","DGbind (kJ/mol)",oenv);
 179   while(do_enx(fp,fr))
 180   {
 181     ct = check_times(fr->t);
 182     if (ct == 0) {
 183       lie = calc_lie(ld,fr->ener,lie_lj,lie_qq,fac_lj,fac_qq);
 184       lieaver += lie;
 185       lieav2  += lie*lie;
 186       nframes ++;
 187       fprintf(out,"%10g  %10g\n",fr->t,lie);
 188     }
 189   }
 190   close_enx(fp);
 191   ffclose(out);
 192   fprintf(stderr,"\n");
 193
 194   if (nframes > 0)
 195     printf("DGbind = %.3f (%.3f)\n",lieaver/nframes,
 196            sqrt(lieav2/nframes-sqr(lieaver/nframes)));
 197
 198   do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
 199
 200   thanx(stderr);
 201
 202   return 0;
 203 }
 204