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