src/tools/mcprop.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.2.0
  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-2004, 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  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "typedefs.h"
  41 #include "random.h"
  42 #include "mcprop.h"
  43 #include "smalloc.h"
  44 #include "vec.h"
  45 #include "futil.h"
  46
  47 void normalise_vec(int nx,real x[])
  48 {
  49   real   fac,nnorm;
  50   int    j;
  51
  52   /* Normalise vector */
  53   nnorm=0.0;
  54   for(j=0; (j<nx); j++)
  55     nnorm+=x[j]*x[j];
  56   fac=1.0/sqrt(nnorm);
  57   for(j=0; (j<nx); j++)
  58     x[j]*=fac;
  59 }
  60
  61 static real do_step(int nx,real x[],int i,int *ig,real step,bool bPlus)
  62 {
  63   static real   r=0;
  64
  65   /* Modify a coordinate */
  66   if (bPlus)
  67     r = rando(ig)*step;
  68   else
  69     r = -r;
  70   x[i] += r;
  71
  72   normalise_vec(nx,x);
  73
  74   return r;
  75 }
  76
  77 real find_min(real x0,real x1,real x2,real y0,real y1,real y2)
  78 {
  79   matrix X,X_1;
  80   rvec   y;
  81   rvec   abc;
  82   real   a2_1,wortel,xx0,yy0,XXX;
  83
  84   X[0][0]=x0*x0, X[0][1]=x0, X[0][2]=1;
  85   X[1][0]=x1*x1, X[1][1]=x1, X[1][2]=1;
  86   X[2][2]=x2*x2, X[2][1]=x2, X[2][2]=1;
  87   y[0]=y0, y[1]=y1, y[2]=y2;
  88
  89   m_inv(X,X_1);
  90   mvmul(X_1,y,abc);
  91
  92   if (abc[0] > 0) {
  93     /* There IS a minimum */
  94     xx0 = -abc[1]/(2.0*abc[0]);
  95     if ((x0 < xx0) && (xx0 < x1))
  96       XXX = xx0;
  97     else {
  98       /* The minimum is not on our interval */
  99       if (xx0 < x0)
 100         XXX = x0;
 101       else
 102         XXX = x2;
 103     }
 104   }
 105   else if (abc[0] < 0) {
 106     if (y0 < y2)
 107       XXX = x0;
 108     else
 109       XXX = x2;
 110   }
 111   else
 112     XXX = x1;
 113
 114   return XXX;
 115 }
 116
 117 void do_mc(FILE *fp,int nx,real x[],real step,real v0,real tol,
 118            int maxsteps,t_propfunc *func)
 119 {
 120   FILE   *ffp[2];
 121   FILE   *ftrj;
 122   int    i,j,k,m,f,n,ig,cur=0;
 123   bool   bConv,bUp;
 124   real   vtol,r,bmf,*rx[2],valmin,vplusmin[2],stepsize;
 125   double dv,val[2];
 126 #define next (1-cur)
 127
 128   snew(rx[cur], nx);
 129   snew(rx[next],nx);
 130   ffp[0]=fp;
 131   ffp[1]=stderr;
 132
 133   ftrj=ffopen("ftrj.out","w");
 134
 135   for(j=0; (j<nx); j++)
 136     rx[cur][j]=x[j];
 137
 138   /* Random seed */
 139   ig       = 1993;
 140
 141   /* Initial value */
 142   val[cur] = func(nx,x);
 143   vtol     = tol*v0;
 144
 145   for(f=0; (f<2); f++) {
 146     fprintf(ffp[f],"Starting MC in property space, YES!\n\n");
 147     fprintf(ffp[f],"Initial value: %10.3e\n",val[cur]);
 148     fprintf(ffp[f],"Going to do %d steps in %dD space\n",maxsteps,nx);
 149   }
 150   bConv=FALSE;
 151   valmin=val[cur];
 152   for(n=0; (n<maxsteps) && !bConv; ) {
 153     for(i=0; (i<nx)  && !bConv; i++,n++) {
 154
 155       for(m=0; (m<2); m++) {
 156         for(j=0; (j<nx); j++)
 157           rx[next][j]=rx[cur][j];
 158         stepsize=do_step(nx,rx[next],i,&ig,step,1-m);
 159         vplusmin[m]=func(nx,rx[next]);
 160       }
 161       for(j=0; (j<nx); j++)
 162         rx[next][j]=rx[cur][j];
 163       rx[next][i]=find_min(rx[cur][i]+stepsize,rx[cur][i],rx[cur][i]-stepsize,
 164                            vplusmin[0],val[cur],vplusmin[1]);
 165       normalise_vec(nx,rx[next]);
 166       val[next]=func(nx,rx[next]);
 167
 168       bmf=0;
 169       bUp=FALSE;
 170       dv=val[next]-val[cur];
 171       if (dv < 0) {
 172         cur=next;
 173         if (val[cur] < valmin)
 174           valmin=val[cur];
 175         for(k=0; (k<nx); k++) {
 176           x[k]=rx[cur][k];
 177           fprintf(ftrj,"%6.3f  ",x[k]);
 178         }
 179         fprintf(ftrj,"\n");
 180       }
 181       if ((fabs(dv) < vtol) && (val[cur]<=valmin))
 182         bConv=TRUE;
 183       else if ((dv >= 0) && (v0 > 0)) {
 184         r=rando(&ig);
 185         bmf=exp(-dv/v0);
 186         if (bmf < r) {
 187           cur=next;
 188           bUp=TRUE;
 189         }
 190       }
 191       for(f=0; (f<2); f++)
 192         fprintf(ffp[f],"Step %5d, Min: %10.3e,  Cur: %10.3e,  BMF %6.3f %s\n",
 193                 n,valmin,val[cur],bmf,bUp ? "+" : "");
 194     }
 195   }
 196   if (bConv) {
 197     fprintf(fp,"Converged !\n");
 198     fprintf(stderr,"Converged !\n");
 199   }
 200   ffclose(ftrj);
 201 }
 202