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