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