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