src/contrib/optwat.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2006, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <string.h>
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "vec.h"
  46 #include "macros.h"
  47 #include "txtdump.h"
  48 #include "tpxio.h"
  49 #include "enxio.h"
  50 #include "names.h"
  51 #include "statutil.h"
  52 #include "copyrite.h"
  53 #include "random.h"
  54
  55 real ener(matrix P,real e,real e0,int nmol,real kp,real ke,gmx_bool bPScal)
  56 {
  57   if (bPScal)
  58     return (kp*(sqr(P[XX][XX]+P[YY][YY]+P[ZZ][ZZ]-3))+
  59             ke*sqr(e/nmol-e0));
  60   else
  61     return (kp*(sqr(P[XX][XX]-1)+sqr(P[YY][YY]-1)+sqr(P[ZZ][ZZ]-1)+
  62                 sqr(P[XX][YY])+sqr(P[XX][ZZ])+sqr(P[YY][ZZ])) +
  63             ke*sqr(e/nmol-e0));
  64 }
  65
  66 void do_sim(char *enx,
  67             t_topology *top,rvec *x,rvec *v,t_inputrec *ir,matrix box)
  68 {
  69   char *tpx = "optwat.tpr";
  70   char buf[128];
  71
  72   write_tpx(tpx,0,0.0,0.0,ir,box,top->atoms.nr,x,v,NULL,top);
  73   sprintf(buf,"xmdrun -s %s -e %s >& /dev/null",tpx,enx);
  74   system(buf);
  75 }
  76
  77 void get_results(char *enx,real P[],real *epot,int pindex,int eindex)
  78 {
  79   ener_file_t fp_ene;
  80   int      nre,step,ndr,i;
  81   real     t;
  82   gmx_enxnm_t *nms=NULL;
  83   t_enxframe fr;
  84
  85   fp_ene = open_enx(enx,"r");
  86
  87   do_enxnms(fp_ene,&nre,&nms);
  88
  89   /* Read until the last frame */
  90   while (do_enx(fp_ene,&fr));
  91
  92   close_enx(fp_ene);
  93
  94   *epot = fr.ener[eindex].e;
  95   for(i=pindex; (i<pindex+9); i++)
  96     P[i-pindex] = fr.ener[i].e;
  97
  98   sfree(ener);
  99 }
 100
 101 void copy_iparams(int nr,t_iparams dest[],t_iparams src[])
 102 {
 103   memcpy(dest,src,nr*sizeof(dest[0]));
 104 }
 105
 106 void rand_step(FILE *fp,int nr,t_iparams ip[],int *seed,real frac)
 107 {
 108   int  i;
 109   real ff;
 110
 111   do {
 112     i = (int) (rando(seed)*nr);
 113   } while (ip[i].lj.c12 == 0.0);
 114
 115   do {
 116     ff = frac*rando(seed);
 117   } while (ff == 0.0);
 118
 119   if (rando(seed) > 0.5) {
 120     ip[i].lj.c12 *= 1.0+ff;
 121     fprintf(fp,"Increasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
 122   }
 123   else {
 124     ip[i].lj.c12 *= 1.0-ff;
 125     fprintf(fp,"Decreasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
 126   }
 127 }
 128
 129 void pr_progress(FILE *fp,int nit,tensor P,real epot,real eFF,
 130                  double mc_crit,gmx_bool bConv,gmx_bool bAccept)
 131 {
 132   fprintf(fp,"Iter %3d, eFF = %g, Converged = %s,  Accepted = %s\n",
 133           nit,eFF,yesno_names[bConv],yesno_names[bAccept]);
 134   fprintf(fp,"Epot = %g  Pscal = %g,  mc_crit = %g\n",epot,
 135           trace(P)/3,mc_crit);
 136   pr_rvecs(fp,0,"Pres",P,DIM);
 137   fprintf(fp,"-----------------------------------------------------\n");
 138   fflush(fp);
 139 }
 140
 141 int main(int argc,char *argv[])
 142 {
 143   static char *desc[] = {
 144     "[TT]optwat[tt] optimizes the force field parameter set of a molecular crystal",
 145     "to reproduce the pressure tensor and experimental energy.[PAR]",
 146     "Note that for good results the [TT].tpx[tt] file must contain input for a",
 147     "simulated annealing run, or a single point energy calculation at 0 K"
 148   };
 149   t_filenm fnm[] = {
 150     { efTPX, NULL,      NULL,       ffREAD },
 151     { efEDR, "-e",      NULL,       ffRW },
 152     { efLOG, "-g",      NULL,       ffWRITE }
 153   };
 154 #define NFILE asize(fnm)
 155
 156   static real epot0    = -57, tol    = 1,   kT     = 0.0;
 157   static real kp       = 1,   ke     = 100, frac   = 0.1;
 158   static int  maxnit   = 100, eindex = 5,   pindex = 19;
 159   static int  seed     = 1993;
 160   static gmx_bool bPScal   = FALSE;
 161   static t_pargs pa[] = {
 162     { "-epot0",  FALSE, etREAL, {&epot0},
 163       "Potential energy in kJ/mol" },
 164     { "-tol",    FALSE, etREAL, {&tol},
 165       "Tolerance for converging" },
 166     { "-nit",    FALSE, etINT,  {&maxnit},
 167       "Max number of iterations" },
 168     { "-seed",   FALSE, etINT,  {&seed},
 169       "Random seed for MC steps" },
 170     { "-frac",   FALSE, etREAL, {&frac},
 171       "Maximum fraction by which to change parameters. Actual fraction is random between 0 and this parameter" },
 172     { "-pindex", FALSE, etINT,  {&pindex},
 173       "Index of P[X][X] in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
 174     { "-eindex", FALSE, etINT,  {&pindex},
 175       "Index of Epot in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
 176     { "-kp",     FALSE, etREAL, {&kp},
 177       "Force constant for pressure components"},
 178     { "-ke",     FALSE, etREAL, {&ke},
 179       "Force constant for energy component"},
 180     { "-kT",     FALSE, etREAL, {&kT},
 181       "Boltzmann Energy for Monte Carlo" },
 182     { "-pscal",  FALSE, etBOOL, {&bPScal},
 183       "Optimize params for scalar pressure, instead of tensor" }
 184   };
 185
 186   FILE        *fp;
 187   t_topology  top;
 188   t_tpxheader sh;
 189   t_inputrec  ir;
 190   t_iparams   *ip[2];
 191   int         cur=0;
 192 #define next  (1-cur)
 193   rvec        *xx,*vv;
 194   matrix      box;
 195   int         i,step,natoms,nmol,nit,atnr2;
 196   real        t,lambda,epot,eFF[2];
 197   double      mc_crit=0;
 198   gmx_bool        bConverged,bAccept;
 199   tensor      P;
 200
 201   CopyRight(stdout,argv[0]);
 202   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 203                     asize(desc),desc,0,NULL);
 204
 205   /* Read initial topology and coordaintes etc. */
 206   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&sh,TRUE,NULL,NULL);
 207   snew(xx,sh.natoms);
 208   snew(vv,sh.natoms);
 209   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,box,&natoms,
 210            xx,vv,NULL,&top);
 211
 212   /* Open log file and print options */
 213   fp    = ftp2FILE(efLOG,NFILE,fnm,"w");
 214   fprintf(fp,"%s started with the following parameters\n",argv[0]);
 215   fprintf(fp,"epot   = %8g  ke     = %8g  kp     = %8g\n",epot0,ke,kp);
 216   fprintf(fp,"maxnit = %8d  tol    = %8g  seed   = %8d\n",maxnit,tol,seed);
 217   fprintf(fp,"frac   = %8g  pindex = %8d  eindex = %8d\n",frac,pindex,eindex);
 218   fprintf(fp,"kT     = %8g  pscal  = %8s\n",kT,gmx_bool_names[bPScal]);
 219
 220   /* Unpack some topology numbers */
 221   nmol  = top.blocks[ebMOLS].nr;
 222   atnr2 = top.idef.atnr*top.idef.atnr;
 223
 224   /* Copy input params */
 225   snew(ip[cur],atnr2);
 226   snew(ip[next],atnr2);
 227   copy_iparams(atnr2,ip[cur],top.idef.iparams);
 228   copy_iparams(atnr2,ip[next],top.idef.iparams);
 229
 230   /* Loop over iterations */
 231   nit = 0;
 232   do {
 233     if (nit > 0) {
 234       /* Do random step */
 235       rand_step(fp,atnr2,ip[next],&seed,frac);
 236       copy_iparams(atnr2,top.idef.iparams,ip[next]);
 237     }
 238     do_sim(ftp2fn(efEDR,NFILE,fnm),&top,xx,vv,&ir,box);
 239
 240     get_results(ftp2fn(efEDR,NFILE,fnm),P[0],&epot,pindex,eindex);
 241
 242     /* Calculate penalty */
 243     eFF[(nit > 0) ? next : cur]  = ener(P,epot,epot0,nmol,kp,ke,bPScal);
 244
 245     bConverged = (eFF[(nit > 0) ? next : cur] < tol);
 246
 247     if (nit > 0) {
 248       /* Do Metropolis criterium */
 249       if (kT > 0)
 250         mc_crit = exp(-(eFF[next]-eFF[cur])/kT);
 251       bAccept = ((eFF[next] < eFF[cur]) ||
 252                  ((kT > 0) && (mc_crit > rando(&seed))));
 253       pr_progress(fp,nit,P,epot/nmol,eFF[next],mc_crit,
 254                   bConverged,bAccept);
 255       if (bAccept) {
 256         /* Better params! */
 257         cur = next;
 258       }
 259       else {
 260         /* Restore old parameters */
 261         copy_iparams(atnr2,ip[next],ip[cur]);
 262       }
 263     }
 264     else
 265       pr_progress(fp,nit,P,epot/nmol,eFF[cur],mc_crit,bConverged,FALSE);
 266
 267     nit++;
 268   } while (!bConverged && (nit < maxnit));
 269
 270   for(i=0; (i<atnr2); i++)
 271     pr_iparams(fp,F_LJ,&ip[cur][i]);
 272
 273   ffclose(fp);
 274
 275   thanx(stderr);
 276
 277   return 0;
 278 }