3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
47 void normalise_vec(int nx,real x[])
52 /* Normalise vector */
61 static real do_step(int nx,real x[],int i,int *ig,real step,bool bPlus)
65 /* Modify a coordinate */
77 real find_min(real x0,real x1,real x2,real y0,real y1,real y2)
82 real a2_1,wortel,xx0,yy0,XXX;
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;
93 /* There IS a minimum */
94 xx0 = -abc[1]/(2.0*abc[0]);
95 if ((x0 < xx0) && (xx0 < x1))
98 /* The minimum is not on our interval */
105 else if (abc[0] < 0) {
117 void do_mc(FILE *fp,int nx,real x[],real step,real v0,real tol,
118 int maxsteps,t_propfunc *func)
122 int i,j,k,m,f,n,ig,cur=0;
124 real vtol,r,bmf,*rx[2],valmin,vplusmin[2],stepsize;
133 ftrj=ffopen("ftrj.out","w");
135 for(j=0; (j<nx); j++)
142 val[cur] = func(nx,x);
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);
152 for(n=0; (n<maxsteps) && !bConv; ) {
153 for(i=0; (i<nx) && !bConv; i++,n++) {
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]);
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]);
170 dv=val[next]-val[cur];
173 if (val[cur] < valmin)
175 for(k=0; (k<nx); k++) {
177 fprintf(ftrj,"%6.3f ",x[k]);
181 if ((fabs(dv) < vtol) && (val[cur]<=valmin))
183 else if ((dv >= 0) && (v0 > 0)) {
192 fprintf(ffp[f],"Step %5d, Min: %10.3e, Cur: %10.3e, BMF %6.3f %s\n",
193 n,valmin,val[cur],bmf,bUp ? "+" : "");
197 fprintf(fp,"Converged !\n");
198 fprintf(stderr,"Converged !\n");