source: sans/XOP_Dev/MonteCarlo/Metropolis.c

Last change on this file was 997, checked in by srkline, 7 years ago

Updates for Xcode 7, Toolkit 7, and 64-bit version

File size: 7.9 KB
Line 
1/*
2 *  Metropolis.c
3 *  SANSAnalysis
4 *
5 *  Created by Steve Kline on 10/16/08.
6 *  Copyright 2008 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10
11#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
12#include "Metropolis.h"
13#include "MonteCarlo.h"
14
15
16int
17MetropolisX(Metropolis_ParamsPtr p) {
18        double *UofR;                           /* pointer to double precision wave data */
19        double *xx;                             /* pointer to double precision wave data */
20        double *yy;                             /* pointer to double precision wave data */
21        double *zz;                             /* pointer to double precision wave data */
22        double *energy;                         /* pointer to double precision wave data */
23//      double radius;                         
24        double tVox;                           
25        CountInt nVox;                  // these are declared double in the struct definition!
26        CountInt numIter;
27        CountInt gCount;
28//      double retVal;                          //return value
29       
30        // for accessing the wave data, direct method (see the WaveAccess example XOP)
31//      waveHndl wavH;
32       
33        CountInt jj,kk,num,good;
34        double Uold,Unew;
35        double step,delta,bolt;
36        double ratio,finalEnergy;
37    CountInt xnew,ynew,znew;
38        float ran;
39        SInt32 seed;
40       
41        char buf[256];
42               
43
44        /* check that wave handles are all valid */
45        if (p->UofRH == NIL) {
46                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
47                return(NON_EXISTENT_WAVE);
48        }
49        if (p->xxH == NIL) {
50                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
51                return(NON_EXISTENT_WAVE);
52        }       
53        if (p->yyH == NIL) {
54                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
55                return(NON_EXISTENT_WAVE);
56        }
57        if (p->zzH == NIL) {
58                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
59                return(NON_EXISTENT_WAVE);
60        }
61        if (p->energyH == NIL) {
62                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
63                return(NON_EXISTENT_WAVE);
64        }
65       
66        p->retVal = 0.0;
67       
68// trusting that all inputs are double PRECISION WAVES!!!
69        UofR = (double*)WaveData(p->UofRH);
70        xx = (double*)WaveData(p->xxH);
71        yy = (double*)WaveData(p->yyH);
72        zz = (double*)WaveData(p->zzH);
73        energy = (double*)WaveData(p->energyH);
74       
75        //input variables
76        tVox = (double)p->tVox;
77        nVox = (CountInt)(p->nVox);
78        numIter = (CountInt)(p->nIter);
79        gCount = (CountInt)(p->gCount);
80       
81        num = WavePoints(p->xxH);               //number of points in x,y,z
82       
83        seed = (SInt32)energy[0];                       // pick up a seed from somewhere
84       
85//              sprintf(buf, "XOP number of points = %lld\r", (SInt64)num);
86//              XOPNotice(buf);
87       
88        if(seed >= 0) {
89                seed = -1234509876;
90        }
91       
92        ran = ran3(&seed);              //initialize the random sequence by passing in a negative value
93        seed = 12348765;                //non-negative after that does nothing
94       
95        good = 0;               // # of accepted moves
96       
97        step = tVox;            //default step?? should this be in terms of radius?
98       
99        // do some number of iterations (kk)
100        for(kk=0;kk<numIter;kk+=1) {
101        //at each iteration, loop over all of the particles
102                if( (kk*10 % numIter) == 0) {
103                        sprintf(buf, "Iteration %lld of %lld\r", (SInt64)kk,(SInt64)numIter);
104                        XOPNotice(buf);
105                }
106                       
107                good = 0;               //re-initialize the number of good moves
108                       
109                for(jj=0;jj<num;jj+=1) {
110                               
111                        //if(mod(jj, num/100 ) == 0)
112                        //      Print "sphere = ",jj," of ",num
113                        //endif
114                               
115                        // 2- calculate the initial energy
116                        Uold = M_energy(UofR,xx,yy,zz,tVox,nVox,xx[jj],yy[jj],zz[jj],jj,num);
117                               
118                        // 3- move a sphere randomly (allow fractional positions at this point)
119                        //displacement d=(ran - 0.5)*dmax in each xyz. this gives forward or backward steps.
120                        // ran is from 0,1
121                        ran = ran3(&seed);
122                        xnew = (CountInt)MC_round(xx[jj] + ((double)ran-0.5)*step); // double new x coordinate for j-particle
123                        ran = ran3(&seed);     
124                        ynew = (CountInt)MC_round(yy[jj] + ((double)ran-0.5)*step); // new y coordinate for j-particle
125                        ran = ran3(&seed);
126                        znew = (CountInt)MC_round(zz[jj] + ((double)ran-0.5)*step); // new z coordinate for j-particle
127                               
128                        // nix the fractional positions -- gives misleading energy calculations, esp. with hard spheres                 
129                        // test point converted to round above
130                               
131                        // enforce in-bounds (periodic)
132                        if( (xnew>=nVox) || (xnew<0) ) {
133                                xnew = labs(nVox - labs(xnew));
134                                }
135                        if( (ynew>=nVox) || (ynew<0) ) {
136                                ynew = labs(nVox - labs(ynew));
137                        }
138                        if( (znew>=nVox) || (znew<0) ) {
139                                znew = labs(nVox - labs(znew));
140                        }
141                                                       
142                        // 4- calculate the new energy
143                        Unew = M_energy(UofR,xx,yy,zz,tVox,nVox,xnew,ynew,znew,jj,num);
144                       
145                        // 5- decide whether to accept the move:
146                        //              if dU <(=?) 0, accept the move
147                        //              if dU > 0, calculate W=exp(-dU)
148                        //              then if random num > W, accept new configuration       
149                        // 6- update the current energy, go back to step 3 and repeat
150                       
151                        delta = Unew-Uold;
152                        /*
153                        if( (kk*10 % numIter) == 0 && jj==0) {
154                                sprintf(buf, "XOP delta = %lf\r", delta);
155                                XOPNotice(buf);
156                        }
157                         */
158                        if(delta <= 0) {
159                                //move is good
160                                Uold = Unew;
161                                xx[jj] = (double)xnew;
162                                yy[jj] = (double)ynew;
163                                zz[jj] = (double)znew;
164                                good += 1;
165                        } else {
166                                bolt = exp(-delta);
167                                ran = ran3(&seed);
168                                if ((double)ran < bolt)  {      // accept anyways
169                                        xx[jj] = (double)xnew; 
170                                        yy[jj] = (double)ynew; 
171                                        zz[jj] = (double)znew;
172                                        good += 1;
173                                        Uold = Unew;
174                                }
175                        }
176                                                                               
177                }               //jj loop over particles
178                                                                               
179                // calculate the final total energy after each pass, and save it
180                if( (kk*10 % numIter) == 0) {
181                        finalEnergy=M_TotalEnergy(UofR,xx,yy,zz,tVox,nVox,num);
182                        energy[gCount+kk+1] = finalEnergy;
183                }
184                                                                                       
185                // calculate the fraction of accepted moves
186                ratio = (double)good/(double)num*100;
187                                                                                       
188                if (ratio > 50.0) { //! adjust the step for randomization
189                        step = step*1.05;
190                        if (step > 2*tVox) {
191                                step = 2*tVox;
192                        }
193                } else {
194                        step = step*0.95;
195                        if (step < 0.1*tVox) {
196                                step = 0.1*tVox;
197                        }
198                } //! end adjust the step for randomization
199
200                /*
201                if( (kk*10 % numIter) == 0) {
202                        sprintf(buf, "XOP percent of accepted moves = %lf\r", ratio);
203                        XOPNotice(buf);
204                        sprintf(buf, "XOP step = %lf\r", step);
205                        XOPNotice(buf);
206                }
207                 */
208
209                                                       
210        }               //kk total iterations
211       
212        finalEnergy=M_TotalEnergy(UofR,xx,yy,zz,tVox,nVox,num);
213        sprintf(buf, "XOP final energy = %lf\r", finalEnergy);
214        XOPNotice(buf);
215        return(0);
216}
217
218
219
220// total energy of the system - without double counting
221double
222M_TotalEnergy(double ur[], double xx[], double yy[], double zz[], double tVox, CountInt nVox, CountInt num) {
223
224        IndexInt ii,jj;
225        double energy,dx,dy,dz,dist;
226
227        energy=0;
228
229        for(ii=0;ii<num;ii+=1) {
230                for(jj=ii+1;jj<num;jj+=1) {
231                        dx = xx[jj] - xx[ii]; //! calculates distance
232                        if(fabs(dx) > nVox/2) {
233                                dx = (nVox-1) - fabs(dx);               //enforce periodic conditions
234                        }
235
236                        dy = yy[jj] - yy[ii]; //! calculates distance
237                        if(fabs(dx) > nVox/2) {
238                                dy = (nVox-1) - fabs(dy);               //enforce periodic conditions
239                        }
240
241                        dz = zz[jj] - zz[ii]; //! calculates distance
242                        if(fabs(dz) > nVox/2) {
243                                dz = (nVox-1) - fabs(dz);               //enforce periodic conditions
244                        }
245
246                        dist = sqrt(dx*dx+dy*dy+dz*dz);// ! distance (in box units)
247                        dist *= tVox;           //real units
248
249        //                      energy += UofR[round(dist/tVox)]                 //! calc. the energy           
250                        energy += ur[MC_round(dist)];            //! calc. the energy           
251
252                }
253        }
254
255        return(energy);
256}
257
258
259
260
261double
262M_energy(double ur[], double xx[], double yy[], double zz[], double tVox, CountInt nVox, double rx, double ry, double rz, CountInt jj,CountInt num)
263{
264        IndexInt ii;
265        double dx,dy,dz,dist,energy;
266
267        energy = 0.0; //! initial energy of 0.0
268
269        for(ii=0;ii<num;ii+=1) {
270                if (ii != jj)  {                //! if i and j are the same point then skip
271                       
272                        dx = rx - xx[ii]; //! calculates distance
273                        if(fabs(dx) > nVox/2) {
274                                dx = (nVox-1) - fabs(dx);               //enforce periodic conditions
275                        }
276
277                        dy = ry - yy[ii]; //! calculates distance
278                        if(fabs(dy) > nVox/2) {
279                                dy = (nVox-1) - fabs(dy);       //enforce periodic conditions
280                        }
281
282                        dz = rz - zz[ii]; //! calculates distance
283                        if(fabs(dz) > nVox/2) {
284                                dz = (nVox-1) - fabs(dz);               //enforce periodic conditions
285                        }
286
287                        dist = sqrt(dx*dx+dy*dy+dz*dz);// ! distance (in box units)
288                        dist *= tVox;           //real units
289
290                        //                      energy += UofR[round(dist/tVox)]                 //! calc. the energy
291                        energy += ur[MC_round(dist)];            //! calc. the energy
292                }
293        }
294
295        return(energy);
296} 
Note: See TracBrowser for help on using the repository browser.