Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2010-12-29 aec47dec1b43faa9cd7ad12821d85882a79aa5a6
commit | author | age
aec47d 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_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn){
15 /*c  Vertex and triangle (lm and lp) indexing for bond flip:
16 c      +----- k-------+              +----- k ------+
17 c      |lm1 / | \ lp1 |              |lm1 /   \ lp1 |
18 c      |  /   |   \   |              |  /       \   |
19 c      |/     |     \ |     FLIP     |/    lm     \ |
20 c     km  lm  | lp   kp    --->      km ---------- kp  
21 c      |\     |     / |              |\    lp     / |  
22 c      |  \   |   /   |              |  \       /   |
23 c      |lm2 \ | / lp2 |              |lm2 \   / lp2 |
24 c      +------it------+              +----- it -----+
25 c
26 */
27     ts_vertex *it=bond->vtx1;
28     ts_vertex *k=bond->vtx2;
29     ts_uint nei,neip,neim;
30     ts_uint i,j;
31     ts_double oldenergy, delta_energy;
32  //   ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lp2=NULL, *lm1=NULL, *lm2=NULL;
33
34     ts_vertex *kp,*km;
35
36     if(it->neigh_no< 3) return TS_FAIL;
37     if(k->neigh_no< 3) return TS_FAIL;
38     if(k==NULL || it==NULL){
39         fatal("In bondflip, number of neighbours of k or it is less than 3!",999);
40     }
41
42
43     for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k 
44         if(it->neigh[i]==k){
45             nei=i;
46             break;
47         }
48     }
49     neip=nei+1;  // I don't like it.. Smells like I must have it in correct order
50     neim=nei-1;
51     if(neip>=it->neigh_no) neip=0;
52     if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not
53 there the neim is never <0 !!! */
54   //  fprintf(stderr,"The numbers are: %u %u\n",neip, neim);
55     km=it->neigh[neim];  // We located km and kp
56     kp=it->neigh[neip];
57
58     if(km==NULL || kp==NULL){
59         fatal("In bondflip, cannot determine km and kp!",999);
60     }
61
62   //  fprintf(stderr,"I WAS HERE! after the 4 vertices are known!\n");
63
64 /* test if the membrane is wrapped too much, so that kp is nearest neighbour of
65  * km. If it is true, then don't flip! */
66     for(i=0;i<km->neigh_no;i++){
67         if(km->neigh[i] == kp) return TS_FAIL;
68     }
69  //   fprintf(stderr,"Membrane didn't wrap too much.. Continue.\n");
70 /* if bond would be too long, return... */
71     if(vertex_distance_sq(km,kp) > vesicle->dmax ) return TS_FAIL;
72  //   fprintf(stderr,"Bond will not be too long.. Continue.\n");
73
74 /* we make a bond flip. this is different than in original fortran */
75 // 0. step. Get memory prior the flip
76   oldenergy=0;
77   oldenergy+=k->xk* k->energy;
78   oldenergy+=kp->xk* kp->energy;
79   oldenergy+=km->xk* km->energy;
80   oldenergy+=it->xk* it->energy;
81 //  for(i=0;i<k->neigh_no;i++) oldenergy+=k->neigh[i]->xk*k->neigh[i]->energy;
82 //  for(i=0;i<kp->neigh_no;i++) oldenergy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
83 //  for(i=0;i<km->neigh_no;i++) oldenergy+=km->neigh[i]->xk*km->neigh[i]->energy;
84 //  for(i=0;i<it->neigh_no;i++) oldenergy+=it->neigh[i]->xk*it->neigh[i]->energy;
85 /*
86 fprintf(stderr,"*** Naslov k=%d\n",k);
87 fprintf(stderr,"*** Naslov it=%d\n",it);
88 fprintf(stderr,"*** Naslov km=%d\n",km);
89 fprintf(stderr,"*** Naslov kp=%d\n",kp);
90
91 for(i=0;i<k->neigh_no;i++)
92     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
93 for(i=0;i<it->neigh_no;i++)
94     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
95
96 for(i=0;i<km->neigh_no;i++)
97     fprintf(stderr,"km sosed=%d\n",km->neigh[i]);
98 for(i=0;i<kp->neigh_no;i++)
99     fprintf(stderr,"kp sosed=%d\n",kp->neigh[i]);
100
101
102 */
103   //  fprintf(stderr,"I WAS HERE! Before bondflip!\n");
104     ts_flip_bond(k,it,km,kp, bond);
105    // fprintf(stderr,"I WAS HERE! Bondflip successful!\n");
106
107 /* Calculating the new energy */
108   delta_energy=0;
109   for(i=0;i<k->neigh_no;i++) energy_vertex(k->neigh[i]);
110   for(i=0;i<kp->neigh_no;i++) energy_vertex(kp->neigh[i]);
111   for(i=0;i<km->neigh_no;i++) energy_vertex(km->neigh[i]);
112   for(i=0;i<it->neigh_no;i++) energy_vertex(it->neigh[i]);
113   delta_energy+=k->xk* k->energy;
114   delta_energy+=kp->xk* kp->energy;
115   delta_energy+=km->xk* km->energy;
116   delta_energy+=it->xk* it->energy;
117 //  for(i=0;i<k->neigh_no;i++) delta_energy+=k->neigh[i]->xk*k->neigh[i]->energy;
118 //  for(i=0;i<kp->neigh_no;i++) delta_energy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
119 //  for(i=0;i<km->neigh_no;i++) delta_energy+=km->neigh[i]->xk*km->neigh[i]->energy;
120 //  for(i=0;i<it->neigh_no;i++) delta_energy+=it->neigh[i]->xk*it->neigh[i]->energy;
121   delta_energy-=oldenergy;
122  // fprintf(stderr,"I WAS HERE! Got energy!\n");
123 /* MONTE CARLO */
124     if(delta_energy>=0){
125 #ifdef TS_DOUBLE_DOUBLE
126         if(exp(-delta_energy)< drand48() )
127 #endif
128 #ifdef TS_DOUBLE_FLOAT
129         if(expf(-delta_energy)< (ts_float)drand48())
130 #endif
131 #ifdef TS_DOUBLE_LONGDOUBLE
132         if(expl(-delta_energy)< (ts_ldouble)drand48())
133 #endif
134         {
135             //not accepted, reverting changes
136        //     fprintf(stderr,"Failed to move, due to MC\n");
137
138 //            ts_flip_bond(km,kp,it,k, bond);
139             ts_flip_bond(kp,km,k,it, bond);
140                 
141
142 /*
143 fprintf(stderr,"*** Naslov k=%d\n",k);
144 fprintf(stderr,"*** Naslov it=%d\n",it);
145 fprintf(stderr,"*** Naslov km=%d\n",km);
146 fprintf(stderr,"*** Naslov kp=%d\n",kp);
147 for(i=0;i<k->neigh_no;i++)
148     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
149 for(i=0;i<it->neigh_no;i++)
150     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
151
152
153 for(i=0;i<km->neigh_no;i++)
154     fprintf(stderr,"km sosed=%d\n",km->neigh[i]);
155 for(i=0;i<kp->neigh_no;i++)
156     fprintf(stderr,"kp sosed=%d\n",kp->neigh[i]);
157 */
158
159
160
161         //    fprintf(stderr,"Reverted condition!\n");
162             return TS_FAIL;
163         }
164     }
165         //    fprintf(stderr,"Success\n");
166
167
168 /* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
169     return TS_SUCCESS;
170 }
171
172
173 ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp,
174 ts_bond *bond){
175
176     ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
177     ts_uint i,j, lmidx, lpidx;
178 if(k==NULL || it==NULL || km==NULL || kp==NULL){
179     fatal("ts_flip_bond: You called me with invalid pointers to vertices",999);
180 }
181 // 1. step. We find lm and lp from k->tristar !
182     for(i=0;i<it->tristar_no;i++){
183         for(j=0;j<k->tristar_no;j++){
184             if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
185                 if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
186 == km || it->tristar[i]->vertex[2]== km )){
187                 lm=it->tristar[i];
188          //       lmidx=i;
189                 }
190                 else
191                 {
192                 lp=it->tristar[i];
193          //       lpidx=i;
194                 }
195
196             }
197         }
198     }
199 if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999);
200
201 //we look for important triangles lp1 and lm2.
202
203  for(i=0;i<k->tristar_no;i++){
204         for(j=0;j<kp->tristar_no;j++){
205                 if((k->tristar[i] == kp->tristar[j]) && k->tristar[i]!=lp){ //ce gre za skupen trikotnik
206                     lp1=k->tristar[i];
207             }
208         }
209 }
210
211  for(i=0;i<it->tristar_no;i++){
212         for(j=0;j<km->tristar_no;j++){
213             if((it->tristar[i] == km->tristar[j]) && it->tristar[i]!=lm){ //ce gre za skupen trikotnik
214                     lm2=it->tristar[i];
215             } 
216         }
217     }
218 /*
219 // DEBUG TESTING!
220 fprintf(stderr,"*** Naslov k=%d\n",k);
221 fprintf(stderr,"*** Naslov it=%d\n",it);
222 fprintf(stderr,"*** Naslov km=%d\n",km);
223 fprintf(stderr,"*** Naslov kp=%d\n",kp);
224
225 for(i=0;i<k->neigh_no;i++)
226     fprintf(stderr,"k sosed=%d\n",k->neigh[i]);
227 for(i=0;i<it->neigh_no;i++)
228     fprintf(stderr,"it sosed=%d\n",it->neigh[i]);
229
230
231 // END DEBUG TESTING!
232 */
233 if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999);
234
235
236 //fprintf(stderr,"1. step: lm, lm2, lp1 and lp found!\n");
237
238 /*
239 //DEBUG TESTING
240 fprintf(stderr,"--- Naslov lm=%d",lm);
241
242
243 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm->vertex[0],lm->vertex[1], lm->vertex[2]);
244 fprintf(stderr,"--- Naslov lp=%d",lp);
245 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp->vertex[0],lp->vertex[1], lp->vertex[2]);
246 fprintf(stderr,"--- Naslov lm2=%d",lm2);
247 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm2->vertex[0],lm2->vertex[1], lm2->vertex[2]);
248 fprintf(stderr,"--- Naslov lp1=%d",lp1);
249 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp1->vertex[0],lp1->vertex[1], lp1->vertex[2]);
250
251 for(i=0;i<lm->neigh_no;i++)
252     fprintf(stderr,"lm sosed=%d\n",lm->neigh[i]);
253 for(i=0;i<lp->neigh_no;i++)
254     fprintf(stderr,"lp sosed=%d\n",lp->neigh[i]);
255 // END DEBUG TESTING
256 */
257 /*
258 // DEBUG TESTING!
259
260 for(i=0;i<3;i++){
261
262     if(lp1->neigh[i]==lp) fprintf(stderr,"Nasel sem par lp1->lp\n");
263     if(lp->neigh[i]==lp1) fprintf(stderr,"Nasel sem par lp->lp1\n");
264     if(lm2->neigh[i]==lm) fprintf(stderr,"Nasel sem par lm2->lm\n");
265     if(lm->neigh[i]==lm2) fprintf(stderr,"Nasel sem par lm->lm2\n");
266 }
267 // END DEBUG TESTING!
268 */
269
270
271 // 2. step. We change the triangle vertices... (actual bond flip)
272     for(i=0;i<3;i++) if(lm->vertex[i]== it) lm->vertex[i]= kp;
273     for(i=0;i<3;i++) if(lp->vertex[i]== k) lp->vertex[i]= km;
274 //fprintf(stderr,"2. step: actual bondflip made\n");
275 // 2a. step. If any changes in triangle calculations must be done, do it here!
276 //   * normals are recalculated here
277     triangle_normal_vector(lp);
278     triangle_normal_vector(lm);
279 // 3. step. Correct neighbours in vertex_list
280
281
282             vertex_remove_neighbour(k,it);
283             vertex_remove_neighbour(it,k);
284             //Tukaj pa nastopi tezava... Kam dodati soseda?
285             vertex_insert_neighbour(km,kp,k);
286             vertex_insert_neighbour(kp,km,it);
287 //            vertex_add_neighbour(km,kp); //pazi na vrstni red.
288 //            vertex_add_neighbour(kp,km);
289 //fprintf(stderr,"3. step: vertex neighbours corrected\n");
290
291 // 3a. step. If any changes to ts_vertex, do it here!
292 //   bond_length calculatons not required for it is done in energy.c
293
294 // 4. step. Correct bond_list (don't know why I still have it!)
295             bond->vtx1=km;
296             bond->vtx2=kp;
297 //fprintf(stderr,"4. step: bondlist corrected\n");
298
299
300 // 5. step. Correct neighbouring triangles 
301    
302     triangle_remove_neighbour(lp,lp1);
303    // fprintf(stderr,".\n");
304     triangle_remove_neighbour(lp1,lp);
305   //  fprintf(stderr,".\n");
306     triangle_remove_neighbour(lm,lm2);
307   //  fprintf(stderr,".\n");
308     triangle_remove_neighbour(lm2,lm);
309    
310     triangle_add_neighbour(lm,lp1);    
311     triangle_add_neighbour(lp1,lm);
312     triangle_add_neighbour(lp,lm2);  //Vrstni red?!
313     triangle_add_neighbour(lm2,lp);
314
315 //fprintf(stderr,"5. step: triangle neigbours corrected\n");
316
317
318 // 6. step. Correct tristar for vertices km, kp, k and it
319             vertex_add_tristar(km,lp);  // Preveri vrstni red!
320             vertex_add_tristar(kp,lm);
321             vertex_remove_tristar(it,lm);
322             vertex_remove_tristar(k,lp);
323 //fprintf(stderr,"6. step: tristar corrected\n");
324
325 /*
326 //DEBUG TESTING
327 fprintf(stderr,"--- Naslov lm=%d",lm);
328
329
330 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm->vertex[0],lm->vertex[1], lm->vertex[2]);
331 fprintf(stderr,"--- Naslov lp=%d",lp);
332 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp->vertex[0],lp->vertex[1], lp->vertex[2]);
333 fprintf(stderr,"--- Naslov lm2=%d",lm2);
334 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lm2->vertex[0],lm2->vertex[1], lm2->vertex[2]);
335 fprintf(stderr,"--- Naslov lp1=%d",lp1);
336 fprintf(stderr,"   vtxs(%d, %d, %d)\n",lp1->vertex[0],lp1->vertex[1], lp1->vertex[2]);
337
338 for(i=0;i<lm->neigh_no;i++)
339     fprintf(stderr,"lm sosed=%d\n",lm->neigh[i]);
340 for(i=0;i<lp->neigh_no;i++)
341     fprintf(stderr,"lp sosed=%d\n",lp->neigh[i]);
342 // END DEBUG TESTING
343 */
344   energy_vertex(k);
345   energy_vertex(kp);
346   energy_vertex(km);
347   energy_vertex(it);
348
349
350 // END modifications to data structure!
351
352
353     return TS_SUCCESS;
354 }