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