src/local/optwat.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_optwat_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "smalloc.h"
  33 #include "vec.h"
  34 #include "macros.h"
  35 #include "txtdump.h"
  36 #include "tpxio.h"
  37 #include "enxio.h"
  38 #include "names.h"
  39 #include "statutil.h"
  40 #include "copyrite.h"
  41 #include "random.h"
  42
  43 real ener(matrix P,real e,real e0,int nmol,real kp,real ke,bool bPScal)
  44 {
  45   if (bPScal)
  46     return (kp*(sqr(P[XX][XX]+P[YY][YY]+P[ZZ][ZZ]-3))+
  47             ke*sqr(e/nmol-e0));
  48   else
  49     return (kp*(sqr(P[XX][XX]-1)+sqr(P[YY][YY]-1)+sqr(P[ZZ][ZZ]-1)+
  50                 sqr(P[XX][YY])+sqr(P[XX][ZZ])+sqr(P[YY][ZZ])) +
  51             ke*sqr(e/nmol-e0));
  52 }
  53
  54 void do_sim(char *enx,
  55             t_topology *top,rvec *x,rvec *v,t_inputrec *ir,matrix box)
  56 {
  57   char *tpx = "optwat.tpr";
  58   char buf[128];
  59
  60   write_tpx(tpx,0,0.0,0.0,ir,box,top->atoms.nr,x,v,NULL,top);
  61   sprintf(buf,"xmdrun -s %s -e %s >& /dev/null",tpx,enx);
  62   system(buf);
  63 }
  64
  65 void get_results(char *enx,real P[],real *epot,int pindex,int eindex)
  66 {
  67   int      fp_ene;
  68   char     **nms=NULL;
  69   int      nre,step,ndr,i;
  70   real     t;
  71   t_energy *ener;
  72
  73   fp_ene = open_enx(enx,"r");
  74
  75   do_enxnms(fp_ene,&nre,&nms);
  76   snew(ener,nre);
  77
  78   /* Read until the last frame */
  79   while (do_enx(fp_ene,&t,&step,&nre,ener,&ndr,NULL));
  80
  81   close_enx(fp_ene);
  82
  83   *epot = ener[eindex].e;
  84   for(i=pindex; (i<pindex+9); i++)
  85     P[i-pindex] = ener[i].e;
  86
  87   sfree(ener);
  88 }
  89
  90 void copy_iparams(int nr,t_iparams dest[],t_iparams src[])
  91 {
  92   memcpy(dest,src,nr*sizeof(dest[0]));
  93 }
  94
  95 void rand_step(FILE *fp,int nr,t_iparams ip[],int *seed,real frac)
  96 {
  97   int  i;
  98   real ff;
  99
 100   do {
 101     i = (int) (rando(seed)*nr);
 102   } while (ip[i].lj.c12 == 0.0);
 103
 104   do {
 105     ff = frac*rando(seed);
 106   } while (ff == 0.0);
 107
 108   if (rando(seed) > 0.5) {
 109     ip[i].lj.c12 *= 1.0+ff;
 110     fprintf(fp,"Increasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
 111   }
 112   else {
 113     ip[i].lj.c12 *= 1.0-ff;
 114     fprintf(fp,"Decreasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
 115   }
 116 }
 117
 118 void pr_progress(FILE *fp,int nit,tensor P,real epot,real eFF,
 119                  double mc_crit,bool bConv,bool bAccept)
 120 {
 121   fprintf(fp,"Iter %3d, eFF = %g, Converged = %s,  Accepted = %s\n",
 122           nit,eFF,yesno_names[bConv],yesno_names[bAccept]);
 123   fprintf(fp,"Epot = %g  Pscal = %g,  mc_crit = %g\n",epot,
 124           trace(P)/3,mc_crit);
 125   pr_rvecs(fp,0,"Pres",P,DIM);
 126   fprintf(fp,"-----------------------------------------------------\n");
 127   fflush(fp);
 128 }
 129
 130 int main(int argc,char *argv[])
 131 {
 132   static char *desc[] = {
 133     "optwat optimizes the force field parameter set of a molecular crystal",
 134     "to reproduce the pressure tensor and experimental energy.[PAR]",
 135     "Note that for good results the tpx file must contain input for a",
 136     "simulated annealing run, or a single point energy calculation at 0 K"
 137   };
 138   t_filenm fnm[] = {
 139     { efTPX, NULL,      NULL,       ffREAD },
 140     { efENX, "-e",      NULL,       ffRW },
 141     { efLOG, "-g",      NULL,       ffWRITE }
 142   };
 143 #define NFILE asize(fnm)
 144
 145   static real epot0    = -57, tol    = 1,   kT     = 0.0;
 146   static real kp       = 1,   ke     = 100, frac   = 0.1;
 147   static int  maxnit   = 100, eindex = 5,   pindex = 19;
 148   static int  seed     = 1993;
 149   static bool bPScal   = FALSE;
 150   static t_pargs pa[] = {
 151     { "-epot0",  FALSE, etREAL, {&epot0},
 152       "Potential energy in kJ/mol" },
 153     { "-tol",    FALSE, etREAL, {&tol},
 154       "Tolerance for converging" },
 155     { "-nit",    FALSE, etINT,  {&maxnit},
 156       "Max number of iterations" },
 157     { "-seed",   FALSE, etINT,  {&seed},
 158       "Random seed for MC steps" },
 159     { "-frac",   FALSE, etREAL, {&frac},
 160       "Maximum fraction by which to change parameters. Actual fraction is random between 0 and this parameter" },
 161     { "-pindex", FALSE, etINT,  {&pindex},
 162       "Index of P[X][X] in the energy file (check with g_energy and subtract 1)" },
 163     { "-eindex", FALSE, etINT,  {&pindex},
 164       "Index of Epot in the energy file (check with g_energy and subtract 1)" },
 165     { "-kp",     FALSE, etREAL, {&kp},
 166       "Force constant for pressure components"},
 167     { "-ke",     FALSE, etREAL, {&ke},
 168       "Force constant for energy component"},
 169     { "-kT",     FALSE, etREAL, {&kT},
 170       "Boltzmann Energy for Monte Carlo" },
 171     { "-pscal",  FALSE, etBOOL, {&bPScal},
 172       "Optimize params for scalar pressure, instead of tensor" }
 173   };
 174
 175   FILE        *fp;
 176   t_topology  top;
 177   t_tpxheader sh;
 178   t_inputrec  ir;
 179   t_iparams   *ip[2];
 180   int         cur=0;
 181 #define next  (1-cur)
 182   rvec        *xx,*vv;
 183   matrix      box;
 184   int         i,step,natoms,nmol,nit,atnr2;
 185   real        t,lambda,epot,eFF[2];
 186   double      mc_crit=0;
 187   bool        bConverged,bAccept;
 188   tensor      P;
 189
 190   CopyRight(stdout,argv[0]);
 191   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
 192                     asize(desc),desc,0,NULL);
 193
 194   /* Read initial topology and coordaintes etc. */
 195   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&sh);
 196   snew(xx,sh.natoms);
 197   snew(vv,sh.natoms);
 198   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,box,&natoms,
 199            xx,vv,NULL,&top);
 200
 201   /* Open log file and print options */
 202   fp    = ftp2FILE(efLOG,NFILE,fnm,"w");
 203   fprintf(fp,"%s started with the following parameters\n",argv[0]);
 204   fprintf(fp,"epot   = %8g  ke     = %8g  kp     = %8g\n",epot0,ke,kp);
 205   fprintf(fp,"maxnit = %8d  tol    = %8g  seed   = %8d\n",maxnit,tol,seed);
 206   fprintf(fp,"frac   = %8g  pindex = %8d  eindex = %8d\n",frac,pindex,eindex);
 207   fprintf(fp,"kT     = %8g  pscal  = %8s\n",kT,bool_names[bPScal]);
 208
 209   /* Unpack some topology numbers */
 210   nmol  = top.blocks[ebMOLS].nr;
 211   atnr2 = top.idef.atnr*top.idef.atnr;
 212
 213   /* Copy input params */
 214   snew(ip[cur],atnr2);
 215   snew(ip[next],atnr2);
 216   copy_iparams(atnr2,ip[cur],top.idef.iparams);
 217   copy_iparams(atnr2,ip[next],top.idef.iparams);
 218
 219   /* Loop over iterations */
 220   nit = 0;
 221   do {
 222     if (nit > 0) {
 223       /* Do random step */
 224       rand_step(fp,atnr2,ip[next],&seed,frac);
 225       copy_iparams(atnr2,top.idef.iparams,ip[next]);
 226     }
 227     do_sim(ftp2fn(efENX,NFILE,fnm),&top,xx,vv,&ir,box);
 228
 229     get_results(ftp2fn(efENX,NFILE,fnm),P[0],&epot,pindex,eindex);
 230
 231     /* Calculate penalty */
 232     eFF[(nit > 0) ? next : cur]  = ener(P,epot,epot0,nmol,kp,ke,bPScal);
 233
 234     bConverged = (eFF[(nit > 0) ? next : cur] < tol);
 235
 236     if (nit > 0) {
 237       /* Do Metropolis criterium */
 238       if (kT > 0)
 239         mc_crit = exp(-(eFF[next]-eFF[cur])/kT);
 240       bAccept = ((eFF[next] < eFF[cur]) ||
 241                  ((kT > 0) && (mc_crit > rando(&seed))));
 242       pr_progress(fp,nit,P,epot/nmol,eFF[next],mc_crit,
 243                   bConverged,bAccept);
 244       if (bAccept) {
 245         /* Better params! */
 246         cur = next;
 247       }
 248       else {
 249         /* Restore old parameters */
 250         copy_iparams(atnr2,ip[next],ip[cur]);
 251       }
 252     }
 253     else
 254       pr_progress(fp,nit,P,epot/nmol,eFF[cur],mc_crit,bConverged,FALSE);
 255
 256     nit++;
 257   } while (!bConverged && (nit < maxnit));
 258
 259   for(i=0; (i<atnr2); i++)
 260     pr_iparams(fp,F_LJ,&ip[cur][i]);
 261
 262   fclose(fp);
 263
 264   thanx(stdout);
 265
 266   return 0;
 267 }