Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2010-11-27 73771431acbfbd1ffa2f3aeabb21fed9016fca42
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include<math.h>
3 #include "general.h"
4 #include "vertex.h"
5 #include "bond.h"
6 #include "triangle.h"
7 #include "vesicle.h"
8 #include "energy.h"
9 #include "timestep.h"
10 #include "cell.h"
11 //#include "io.h"
12 #include<stdio.h>
13
14 ts_bool single_timestep(ts_vesicle *vesicle){
15     ts_bool retval;
16     ts_double rnvec[3];
17     ts_uint i;
18     for(i=0;i<vesicle->vlist.n;i++){
19         rnvec[0]=drand48();
20         rnvec[1]=drand48();
21         rnvec[2]=drand48();
22         retval=single_verticle_timestep(vesicle,&vesicle->vlist.vertex[i],rnvec);
23     }
24
25     for(i=0;i<vesicle->blist.n;i++){
26         rnvec[0]=drand48();
27         rnvec[1]=drand48();
28         rnvec[2]=drand48();
29         //find a bond and return a pointer to a bond...
30         //call single_bondflip_timestep...
31         retval=single_bondflip_timestep(vesicle,&vesicle->blist.bond[i],rnvec);
32         
33     } 
34
35     return TS_SUCCESS;
36 }
37
38
39
40
41
42 ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double
43 *rn){
44     ts_uint i;
45     ts_double dist;
46     ts_vertex tvtx;
47     ts_bool retval; 
48     ts_uint cellidx; 
49     ts_double xold,yold,zold;
50     ts_double delta_energy,oenergy;
51     ts_vertex *ovtx;
52
53     //randomly we move the temporary vertex
54     tvtx.x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
55     tvtx.y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
56     tvtx.z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
57     //check we if some length to neighbours are too much
58
59     for(i=0;i<vtx->neigh_no;i++){
60         dist=vertex_distance_sq(&tvtx,vtx->neigh[i]);
61         if(dist<1.0 || dist>vesicle->dmax) return TS_FAIL;
62     }
63     //self avoidance check with distant vertices
64      cellidx=vertex_self_avoidance(vesicle, &tvtx);
65     //check occupation number
66      retval=cell_occupation_number_and_internal_proximity(&vesicle->clist,cellidx,vtx,&tvtx);
67     if(retval==TS_FAIL){
68         return TS_FAIL;
69     }
70     
71     //if all the tests are successful, then we update the vertex position
72     xold=vtx->x;
73     yold=vtx->y;
74     zold=vtx->z;
75     ovtx=malloc(sizeof(ts_vertex));
76     vertex_full_copy(ovtx,vtx);
77     vtx->x=tvtx.x;
78     vtx->y=tvtx.y;
79     vtx->z=tvtx.z;
80
81     delta_energy=0;
82     //update the normals of triangles that share bead i.
83     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
84     //energy and curvature
85     energy_vertex(vtx);
86     delta_energy=vtx->xk*(vtx->energy - ovtx->energy);
87     //the same is done for neighbouring vertices
88     for(i=0;i<vtx->neigh_no;i++){
89         oenergy=vtx->neigh[i]->energy;
90         energy_vertex(vtx->neigh[i]);
91         delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
92     }
93   //  fprintf(stderr, "DE=%f\n",delta_energy);
94     //MONTE CARLOOOOOOOO
95     if(delta_energy>=0){
96 #ifdef TS_DOUBLE_DOUBLE
97         if(exp(-delta_energy)< drand48() )
98 #endif
99 #ifdef TS_DOUBLE_FLOAT
100         if(expf(-delta_energy)< (ts_float)drand48())
101 #endif
102 #ifdef TS_DOUBLE_LONGDOUBLE
103         if(expl(-delta_energy)< (ts_ldouble)drand48())
104 #endif
105     {
106     //not accepted, reverting changes
107     vtx->x=xold;
108     vtx->y=yold;
109     vtx->z=zold;
110     //update the normals of triangles that share bead i.
111     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
112     //energy and curvature
113     energy_vertex(vtx);
114     //the same is done for neighbouring vertices
115     for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]);
116   free(ovtx->bond_length);
117     free(ovtx->bond_length_dual);
118     free(ovtx);
119     return TS_FAIL; 
120     }
121 }
122     //END MONTE CARLOOOOOOO
123
124     //TODO: change cell occupation if necessary!
125
126     free(ovtx->bond_length);
127     free(ovtx->bond_length_dual);
128     free(ovtx);
129     return TS_SUCCESS;
130 }
131
132 ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn){
133 /*c  Vertex and triangle (lm and lp) indexing for bond flip:
134 c      +----- k-------+              +----- k ------+
135 c      |lm1 / | \ lp1 |              |lm1 /   \ lp1 |
136 c      |  /   |   \   |              |  /       \   |
137 c      |/     |     \ |     FLIP     |/    lm     \ |
138 c     km  lm  | lp   kp    --->      km ---------- kp  
139 c      |\     |     / |              |\    lp     / |  
140 c      |  \   |   /   |              |  \       /   |
141 c      |lm2 \ | / lp2 |              |lm2 \   / lp2 |
142 c      +------it------+              +----- it -----+
143 c
144 */
145     ts_vertex *it=bond->vtx1;
146     ts_vertex *k=bond->vtx2;
147     ts_uint nei,neip,neim;
148     ts_uint i,j;
149     ts_double oldenergy, delta_energy;
150  //   ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lp2=NULL, *lm1=NULL, *lm2=NULL;
151
152     ts_vertex *kp,*km;
153
154     if(it->neigh_no< 3) return TS_FAIL;
155     if(k->neigh_no< 3) return TS_FAIL;
156     if(k==NULL || it==NULL){
157         fatal("In bondflip, number of neighbours of k or it is less than 3!",999);
158     }
159
160
161     for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k 
162         if(it->neigh[i]==k){
163             nei=i;
164             break;
165         }
166     }
167     neip=nei+1;  // I don't like it.. Smells like I must have it in correct order
168     neim=nei-1;
169     if(neip>=it->neigh_no) neip=0;
170     if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not
171 there the neim is never <0 !!! */
172   //  fprintf(stderr,"The numbers are: %u %u\n",neip, neim);
173     km=it->neigh[neim];  // We located km and kp
174     kp=it->neigh[neip];
175
176     if(km==NULL || kp==NULL){
177         fatal("In bondflip, cannot determine km and kp!",999);
178     }
179
180   //  fprintf(stderr,"I WAS HERE! after the 4 vertices are known!\n");
181
182 /* test if the membrane is wrapped too much, so that kp is nearest neighbour of
183  * km. If it is true, then don't flip! */
184     for(i=0;i<km->neigh_no;i++){
185         if(km->neigh[i] == kp) return TS_FAIL;
186     }
187  //   fprintf(stderr,"Membrane didn't wrap too much.. Continue.\n");
188 /* if bond would be too long, return... */
189     if(vertex_distance_sq(km,kp) > vesicle->dmax ) return TS_FAIL;
190  //   fprintf(stderr,"Bond will not be too long.. Continue.\n");
191
192 /* we make a bond flip. this is different than in original fortran */
193 // 0. step. Get memory prior the flip
194   oldenergy=0;
195   oldenergy+=k->xk* k->energy;
196   oldenergy+=kp->xk* kp->energy;
197   oldenergy+=km->xk* km->energy;
198   oldenergy+=it->xk* it->energy;
199 //  for(i=0;i<k->neigh_no;i++) oldenergy+=k->neigh[i]->xk*k->neigh[i]->energy;
200 //  for(i=0;i<kp->neigh_no;i++) oldenergy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
201 //  for(i=0;i<km->neigh_no;i++) oldenergy+=km->neigh[i]->xk*km->neigh[i]->energy;
202 //  for(i=0;i<it->neigh_no;i++) oldenergy+=it->neigh[i]->xk*it->neigh[i]->energy;
203 /*
204 fprintf(stderr,"*** Naslov k=%d\n",k);
205 fprintf(stderr,"*** Naslov it=%d\n",it);
206 fprintf(stderr,"*** Naslov km=%d\n",km);
207 fprintf(stderr,"*** Naslov kp=%d\n",kp);
208
209 for(i=0;i<k->neigh_no;i++)
210     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
211 for(i=0;i<it->neigh_no;i++)
212     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
213
214 for(i=0;i<km->neigh_no;i++)
215     fprintf(stderr,"km sosed=%d\n",km->neigh[i]);
216 for(i=0;i<kp->neigh_no;i++)
217     fprintf(stderr,"kp sosed=%d\n",kp->neigh[i]);
218
219
220 */
221   //  fprintf(stderr,"I WAS HERE! Before bondflip!\n");
222     ts_flip_bond(k,it,km,kp, bond);
223    // fprintf(stderr,"I WAS HERE! Bondflip successful!\n");
224
225 /* Calculating the new energy */
226   delta_energy=0;
227   for(i=0;i<k->neigh_no;i++) energy_vertex(k->neigh[i]);
228   for(i=0;i<kp->neigh_no;i++) energy_vertex(kp->neigh[i]);
229   for(i=0;i<km->neigh_no;i++) energy_vertex(km->neigh[i]);
230   for(i=0;i<it->neigh_no;i++) energy_vertex(it->neigh[i]);
231   delta_energy+=k->xk* k->energy;
232   delta_energy+=kp->xk* kp->energy;
233   delta_energy+=km->xk* km->energy;
234   delta_energy+=it->xk* it->energy;
235 //  for(i=0;i<k->neigh_no;i++) delta_energy+=k->neigh[i]->xk*k->neigh[i]->energy;
236 //  for(i=0;i<kp->neigh_no;i++) delta_energy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
237 //  for(i=0;i<km->neigh_no;i++) delta_energy+=km->neigh[i]->xk*km->neigh[i]->energy;
238 //  for(i=0;i<it->neigh_no;i++) delta_energy+=it->neigh[i]->xk*it->neigh[i]->energy;
239   delta_energy-=oldenergy;
240  // fprintf(stderr,"I WAS HERE! Got energy!\n");
241 /* MONTE CARLO */
242     if(delta_energy>=0){
243 #ifdef TS_DOUBLE_DOUBLE
244         if(exp(-delta_energy)< drand48() )
245 #endif
246 #ifdef TS_DOUBLE_FLOAT
247         if(expf(-delta_energy)< (ts_float)drand48())
248 #endif
249 #ifdef TS_DOUBLE_LONGDOUBLE
250         if(expl(-delta_energy)< (ts_ldouble)drand48())
251 #endif
252         {
253             //not accepted, reverting changes
254        //     fprintf(stderr,"Failed to move, due to MC\n");
255
256 //            ts_flip_bond(km,kp,it,k, bond);
257             ts_flip_bond(kp,km,k,it, bond);
258                 
259
260 /*
261 fprintf(stderr,"*** Naslov k=%d\n",k);
262 fprintf(stderr,"*** Naslov it=%d\n",it);
263 fprintf(stderr,"*** Naslov km=%d\n",km);
264 fprintf(stderr,"*** Naslov kp=%d\n",kp);
265 for(i=0;i<k->neigh_no;i++)
266     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
267 for(i=0;i<it->neigh_no;i++)
268     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
269
270
271 for(i=0;i<km->neigh_no;i++)
272     fprintf(stderr,"km sosed=%d\n",km->neigh[i]);
273 for(i=0;i<kp->neigh_no;i++)
274     fprintf(stderr,"kp sosed=%d\n",kp->neigh[i]);
275 */
276
277
278
279         //    fprintf(stderr,"Reverted condition!\n");
280             return TS_FAIL;
281         }
282     }
283         //    fprintf(stderr,"Success\n");
284
285
286 /* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
287     return TS_SUCCESS;
288 }
289
290
291 ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp,
292 ts_bond *bond){
293
294     ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
295     ts_uint i,j, lmidx, lpidx;
296 if(k==NULL || it==NULL || km==NULL || kp==NULL){
297     fatal("ts_flip_bond: You called me with invalid pointers to vertices",999);
298 }
299 // 1. step. We find lm and lp from k->tristar !
300     for(i=0;i<it->tristar_no;i++){
301         for(j=0;j<k->tristar_no;j++){
302             if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
303                 if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
304 == km || it->tristar[i]->vertex[2]== km )){
305                 lm=it->tristar[i];
306          //       lmidx=i;
307                 }
308                 else
309                 {
310                 lp=it->tristar[i];
311          //       lpidx=i;
312                 }
313
314             }
315         }
316     }
317 if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999);
318
319 //we look for important triangles lp1 and lm2.
320
321  for(i=0;i<k->tristar_no;i++){
322         for(j=0;j<kp->tristar_no;j++){
323                 if((k->tristar[i] == kp->tristar[j]) && k->tristar[i]!=lp){ //ce gre za skupen trikotnik
324                     lp1=k->tristar[i];
325             }
326         }
327 }
328
329  for(i=0;i<it->tristar_no;i++){
330         for(j=0;j<km->tristar_no;j++){
331             if((it->tristar[i] == km->tristar[j]) && it->tristar[i]!=lm){ //ce gre za skupen trikotnik
332                     lm2=it->tristar[i];
333             } 
334         }
335     }
336 /*
337 // DEBUG TESTING!
338 fprintf(stderr,"*** Naslov k=%d\n",k);
339 fprintf(stderr,"*** Naslov it=%d\n",it);
340 fprintf(stderr,"*** Naslov km=%d\n",km);
341 fprintf(stderr,"*** Naslov kp=%d\n",kp);
342
343 for(i=0;i<k->neigh_no;i++)
344     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
345 for(i=0;i<it->neigh_no;i++)
346     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
347
348
349 // END DEBUG TESTING!
350 */
351 if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999);
352
353
354 //fprintf(stderr,"1. step: lm, lm2, lp1 and lp found!\n");
355
356 /*
357 //DEBUG TESTING
358 fprintf(stderr,"--- Naslov lm=%d",lm);
359
360
361 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm->vertex[0],lm->vertex[1], lm->vertex[2]);
362 fprintf(stderr,"--- Naslov lp=%d",lp);
363 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp->vertex[0],lp->vertex[1], lp->vertex[2]);
364 fprintf(stderr,"--- Naslov lm2=%d",lm2);
365 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm2->vertex[0],lm2->vertex[1], lm2->vertex[2]);
366 fprintf(stderr,"--- Naslov lp1=%d",lp1);
367 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp1->vertex[0],lp1->vertex[1], lp1->vertex[2]);
368
369 for(i=0;i<lm->neigh_no;i++)
370     fprintf(stderr,"lm sosed=%d\n",lm->neigh[i]);
371 for(i=0;i<lp->neigh_no;i++)
372     fprintf(stderr,"lp sosed=%d\n",lp->neigh[i]);
373 // END DEBUG TESTING
374 */
375 /*
376 // DEBUG TESTING!
377
378 for(i=0;i<3;i++){
379
380     if(lp1->neigh[i]==lp) fprintf(stderr,"Nasel sem par lp1->lp\n");
381     if(lp->neigh[i]==lp1) fprintf(stderr,"Nasel sem par lp->lp1\n");
382     if(lm2->neigh[i]==lm) fprintf(stderr,"Nasel sem par lm2->lm\n");
383     if(lm->neigh[i]==lm2) fprintf(stderr,"Nasel sem par lm->lm2\n");
384 }
385 // END DEBUG TESTING!
386 */
387
388
389 // 2. step. We change the triangle vertices... (actual bond flip)
390     for(i=0;i<3;i++) if(lm->vertex[i]== it) lm->vertex[i]= kp;
391     for(i=0;i<3;i++) if(lp->vertex[i]== k) lp->vertex[i]= km;
392 //fprintf(stderr,"2. step: actual bondflip made\n");
393 // 2a. step. If any changes in triangle calculations must be done, do it here!
394 //   * normals are recalculated here
395     triangle_normal_vector(lp);
396     triangle_normal_vector(lm);
397 // 3. step. Correct neighbours in vertex_list
398
399
400             vertex_remove_neighbour(k,it);
401             vertex_remove_neighbour(it,k);
402             //Tukaj pa nastopi tezava... Kam dodati soseda?
403             vertex_insert_neighbour(km,kp,k);
404             vertex_insert_neighbour(kp,km,it);
405 //            vertex_add_neighbour(km,kp); //pazi na vrstni red.
406 //            vertex_add_neighbour(kp,km);
407 //fprintf(stderr,"3. step: vertex neighbours corrected\n");
408
409 // 3a. step. If any changes to ts_vertex, do it here!
410 //   bond_length calculatons not required for it is done in energy.c
411
412 // 4. step. Correct bond_list (don't know why I still have it!)
413             bond->vtx1=km;
414             bond->vtx2=kp;
415 //fprintf(stderr,"4. step: bondlist corrected\n");
416
417
418 // 5. step. Correct neighbouring triangles 
419    
420     triangle_remove_neighbour(lp,lp1);
421    // fprintf(stderr,".\n");
422     triangle_remove_neighbour(lp1,lp);
423   //  fprintf(stderr,".\n");
424     triangle_remove_neighbour(lm,lm2);
425   //  fprintf(stderr,".\n");
426     triangle_remove_neighbour(lm2,lm);
427    
428     triangle_add_neighbour(lm,lp1);    
429     triangle_add_neighbour(lp1,lm);
430     triangle_add_neighbour(lp,lm2);  //Vrstni red?!
431     triangle_add_neighbour(lm2,lp);
432
433 //fprintf(stderr,"5. step: triangle neigbours corrected\n");
434
435
436 // 6. step. Correct tristar for vertices km, kp, k and it
437             vertex_add_tristar(km,lp);  // Preveri vrstni red!
438             vertex_add_tristar(kp,lm);
439             vertex_remove_tristar(it,lm);
440             vertex_remove_tristar(k,lp);
441 //fprintf(stderr,"6. step: tristar corrected\n");
442
443 /*
444 //DEBUG TESTING
445 fprintf(stderr,"--- Naslov lm=%d",lm);
446
447
448 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm->vertex[0],lm->vertex[1], lm->vertex[2]);
449 fprintf(stderr,"--- Naslov lp=%d",lp);
450 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp->vertex[0],lp->vertex[1], lp->vertex[2]);
451 fprintf(stderr,"--- Naslov lm2=%d",lm2);
452 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm2->vertex[0],lm2->vertex[1], lm2->vertex[2]);
453 fprintf(stderr,"--- Naslov lp1=%d",lp1);
454 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp1->vertex[0],lp1->vertex[1], lp1->vertex[2]);
455
456 for(i=0;i<lm->neigh_no;i++)
457     fprintf(stderr,"lm sosed=%d\n",lm->neigh[i]);
458 for(i=0;i<lp->neigh_no;i++)
459     fprintf(stderr,"lp sosed=%d\n",lp->neigh[i]);
460 // END DEBUG TESTING
461 */
462   energy_vertex(k);
463   energy_vertex(kp);
464   energy_vertex(km);
465   energy_vertex(it);
466
467
468 // END modifications to data structure!
469
470
471     return TS_SUCCESS;
472 }