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 | |
---|
16 | int |
---|
17 | MetropolisX(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 |
---|
221 | double |
---|
222 | M_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 | |
---|
261 | double |
---|
262 | M_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 | } |
---|