Trisurf Monte Carlo simulator
Samo Penic
2016-05-15 5f4aea824bd6785dfa67f922bbc67a57178ae0d7
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
aec47d 2 #include<stdlib.h>
SP 3 #include<math.h>
4 #include "general.h"
5 #include "vertex.h"
6 #include "bond.h"
7 #include "triangle.h"
8 #include "vesicle.h"
9 #include "energy.h"
10 #include "timestep.h"
11 #include "cell.h"
12 //#include "io.h"
9166cb 13 #include "io.h"
aec47d 14 #include<stdio.h>
SP 15 #include "vertexmove.h"
1ad6d1 16 #include <string.h>
43c042 17 #include "constvol.h"
aec47d 18
fedf2b 19 ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
aec47d 20     ts_uint i;
SP 21     ts_double dist;
22     ts_bool retval; 
23     ts_uint cellidx; 
c0ae90 24     ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0;
ed31fe 25     ts_double costheta,sintheta,phi,r;
1ad6d1 26     //This will hold all the information of vtx and its neighbours
958e0e 27     ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL;
dcd350 28     memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
a63f17 29
SP 30     //Some stupid tests for debugging cell occupation!
31 /*         cellidx=vertex_self_avoidance(vesicle, vtx);
32     if(vesicle->clist->cell[cellidx]==vtx->cell){
33         fprintf(stderr,"Idx match!\n");
34     } else {
35         fprintf(stderr,"***** Idx don't match!\n");
36         fatal("ENding.",1);
37     }
38 */
39
352fad 40         //temporarly moving the vertex
672ae4 41 //    vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
SP 42 //        vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
43 //        vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
44
c0ae90 45 //random move in a sphere with radius stepsize:
ed31fe 46     r=vesicle->stepsize*rn[0];
M 47     phi=rn[1]*2*M_PI;
48     costheta=2*rn[2]-1;
49     sintheta=sqrt(1-pow(costheta,2));
672ae4 50     vtx->x=vtx->x+r*sintheta*cos(phi);
SP 51     vtx->y=vtx->y+r*sintheta*sin(phi);
52     vtx->z=vtx->z+r*costheta;
53
54
c0ae90 55 //distance with neighbours check
8f6a69 56     for(i=0;i<vtx->neigh_no;i++){
352fad 57         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
8f6a69 58         if(dist<1.0 || dist>vesicle->dmax) {
c0ae90 59             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
SP 60             return TS_FAIL;
8f6a69 61         }
aec47d 62     }
304510 63
M 64 // Distance with grafted poly-vertex check:    
65     if(vtx->grafted_poly!=NULL){
66         dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
67         if(dist<1.0 || dist>vesicle->dmax) {
68         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
69         return TS_FAIL;
70         }
71     }
72
fe24d2 73 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 74 // Nucleus penetration check:
75     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
76         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
77         return TS_FAIL;
78     }
79
80 //self avoidance check with distant vertices
81     cellidx=vertex_self_avoidance(vesicle, vtx);
82     //check occupation number
83     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
84
aec47d 85     if(retval==TS_FAIL){
dcd350 86         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
aec47d 87         return TS_FAIL;
SP 88     } 
1ad6d1 89    
SP 90  
c0ae90 91 //if all the tests are successful, then energy for vtx and neighbours is calculated
1ad6d1 92     for(i=0;i<vtx->neigh_no;i++){
dcd350 93     memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
1ad6d1 94     }
aec47d 95
1121fa 96     if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
414b8a 97         for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
c0ae90 98     }
SP 99
100     if(vesicle->tape->constareaswitch==2){
101         for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
102     
103     }
a63f17 104
aec47d 105     delta_energy=0;
fe5069 106     
SP 107 //    vesicle_volume(vesicle);
108 //    fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
43c042 109
aec47d 110     //update the normals of triangles that share bead i.
8f6a69 111     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17 112     oenergy=vtx->energy;
aec47d 113     energy_vertex(vtx);
a63f17 114     delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d 115     //the same is done for neighbouring vertices
8f6a69 116     for(i=0;i<vtx->neigh_no;i++){
SP 117         oenergy=vtx->neigh[i]->energy;
118         energy_vertex(vtx->neigh[i]);
119         delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d 120     }
414b8a 121
1121fa 122     if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch >0){
414b8a 123         for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
fbcbdf 124         if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
414b8a 125     };
43c042 126
c0ae90 127     if(vesicle->tape->constareaswitch==2){
SP 128         /* check whether the darea is gt epsarea */
129         for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
130         if(fabs(vesicle->area+darea-A0)>epsarea){
131             //restore old state.
132              vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
133                 for(i=0;i<vtx->neigh_no;i++){
134                     vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
135                 }
136                     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
137                     //fprintf(stderr,"fajlam!\n");
138                     return TS_FAIL;
139         }
140
141
142     }
1121fa 143
SP 144     if(vesicle->tape->constvolswitch==2){
145         /*check whether the dvol is gt than epsvol */
146             //fprintf(stderr,"DVOL=%1.16e\n",dvol);
147         if(fabs(vesicle->volume+dvol-V0)>epsvol){
148             //restore old state.
149              vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
150                 for(i=0;i<vtx->neigh_no;i++){
151                     vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
152                 }
153                     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
154                     //fprintf(stderr,"fajlam!\n");
155                     return TS_FAIL;
156         }
157
158     } else
fe5069 159 //    vesicle_volume(vesicle);
SP 160 //    fprintf(stderr,"Volume before=%1.16e\n", vesicle->volume);
fbcbdf 161    if(vesicle->tape->constvolswitch == 1){
fe5069 162         retval=constvolume(vesicle, vtx, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
fbcbdf 163         if(retval==TS_FAIL){ // if we couldn't move the vertex to assure constant volume
SP 164             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
165             for(i=0;i<vtx->neigh_no;i++){
166                 vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
167             }
168             for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
fe5069 169  //           fprintf(stderr,"fajlam!\n");
fbcbdf 170             return TS_FAIL;
SP 171         }
fe5069 172 //    vesicle_volume(vesicle);
SP 173 //    fprintf(stderr,"Volume after=%1.16e\n", vesicle->volume);
174 //    fprintf(stderr,"Volume after-dvol=%1.16e\n", vesicle->volume-dvol);
175 //    fprintf(stderr,"Denergy before=%e\n",delta_energy);
176     
fbcbdf 177     delta_energy+=delta_energy_cv;
fe5069 178 //    fprintf(stderr,"Denergy after=%e\n",delta_energy);
fbcbdf 179     }
304510 180 /* No poly-bond energy for now!
fedf2b 181     if(vtx->grafted_poly!=NULL){
M 182         delta_energy+=
183             (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
184             pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
185     }
304510 186 */
314f2d 187 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 188     //MONTE CARLOOOOOOOO
SP 189     if(delta_energy>=0){
190 #ifdef TS_DOUBLE_DOUBLE
3de289 191         if(exp(-delta_energy)< drand48())
aec47d 192 #endif
SP 193 #ifdef TS_DOUBLE_FLOAT
194         if(expf(-delta_energy)< (ts_float)drand48())
195 #endif
196 #ifdef TS_DOUBLE_LONGDOUBLE
197         if(expl(-delta_energy)< (ts_ldouble)drand48())
198 #endif
199     {
200     //not accepted, reverting changes
fbcbdf 201   //  fprintf(stderr,"MC failed\n");
dcd350 202     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 203     for(i=0;i<vtx->neigh_no;i++){
a63f17 204         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 205     }
SP 206     
aec47d 207     //update the normals of triangles that share bead i.
dcd350 208    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1 209
fe5069 210 //    fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
43c042 211     if(vesicle->tape->constvolswitch == 1){
958e0e 212         constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
43c042 213     }
fe5069 214 //    fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
dd5aca 215 //    vesicle_volume(vesicle);
SP 216 //    fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
aec47d 217     return TS_FAIL; 
SP 218     }
219 }
2b14da 220     //accepted    
fbcbdf 221  //   fprintf(stderr,"MC accepted\n");
a63f17 222 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
SP 223     if(vtx->cell!=vesicle->clist->cell[cellidx]){
224         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
225 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
226         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
227         
228     }
2b14da 229
1121fa 230     if(vesicle->tape->constvolswitch == 2){
SP 231     vesicle->volume+=dvol;
232     } else
43c042 233     if(vesicle->tape->constvolswitch == 1){
fbcbdf 234         constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
43c042 235     }
c0ae90 236
SP 237     if(vesicle->tape->constareaswitch==2){
238         vesicle->area+=darea;
239     }
a63f17 240 //    if(oldcellidx);
aec47d 241     //END MONTE CARLOOOOOOO
dd5aca 242 //    vesicle_volume(vesicle);
SP 243 //    fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
aec47d 244     return TS_SUCCESS;
SP 245 }
246
fedf2b 247
M 248 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
249     ts_uint i;
250     ts_bool retval; 
251     ts_uint cellidx; 
304510 252 //    ts_double delta_energy;
fedf2b 253     ts_double costheta,sintheta,phi,r;
304510 254     ts_double dist;
fedf2b 255     //This will hold all the information of vtx and its neighbours
M 256     ts_vertex backupvtx;
304510 257 //    ts_bond backupbond[2];
fedf2b 258     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 259
260     //random move in a sphere with radius stepsize:
261     r=vesicle->stepsize*rn[0];
262     phi=rn[1]*2*M_PI;
263     costheta=2*rn[2]-1;
264     sintheta=sqrt(1-pow(costheta,2));
265     vtx->x=vtx->x+r*sintheta*cos(phi);
266     vtx->y=vtx->y+r*sintheta*sin(phi);
267     vtx->z=vtx->z+r*costheta;
268
269
270     //distance with neighbours check
304510 271     for(i=0;i<vtx->neigh_no;i++){
M 272         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
273         if(dist<1.0 || dist>vesicle->dmax) {
274             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
275             return TS_FAIL;
276         }
277     }
278
279 // Distance with grafted vesicle-vertex check:    
280     if(vtx==poly->vlist->vtx[0]){
281         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
282         if(dist<1.0 || dist>vesicle->dmax) {
283         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
284         return TS_FAIL;
285         }
286     }
287
fedf2b 288
M 289     //self avoidance check with distant vertices
290     cellidx=vertex_self_avoidance(vesicle, vtx);
291     //check occupation number
292     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
293     
294     if(retval==TS_FAIL){
295         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
296         return TS_FAIL;
297     } 
298
299
300     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 301 /* Energy ignored for now!
fedf2b 302     delta_energy=0;
M 303     for(i=0;i<vtx->bond_no;i++){
304         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
305
306         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
307         bond_energy(vtx->bond[i],poly);
308         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
309     }
310
311     if(vtx==poly->vlist->vtx[0]){
312         delta_energy+=
313             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
314             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
315         
316     }
317
318
319     if(delta_energy>=0){
320 #ifdef TS_DOUBLE_DOUBLE
321         if(exp(-delta_energy)< drand48() )
322 #endif
323 #ifdef TS_DOUBLE_FLOAT
324         if(expf(-delta_energy)< (ts_float)drand48())
325 #endif
326 #ifdef TS_DOUBLE_LONGDOUBLE
327         if(expl(-delta_energy)< (ts_ldouble)drand48())
328 #endif
329         {
330     //not accepted, reverting changes
331     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
332     for(i=0;i<vtx->bond_no;i++){
333     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
334     }
335
336     return TS_FAIL; 
337     }
338     }
304510 339 */
fedf2b 340         
M 341 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
342     if(vtx->cell!=vesicle->clist->cell[cellidx]){
343         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
344 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
345         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
346     }
347 //    if(oldcellidx);
348     //END MONTE CARLOOOOOOO
349     return TS_SUCCESS;
350 }
58230a 351
M 352
353
354
355 ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
356     ts_uint i;
357     ts_bool retval; 
358     ts_uint cellidx; 
b30f45 359     ts_double delta_energy;
58230a 360     ts_double costheta,sintheta,phi,r;
M 361     ts_double dist[2];
362     //This will hold all the information of vtx and its neighbours
b30f45 363     ts_vertex backupvtx,backupneigh[2];
58230a 364     ts_bond backupbond[2];
b30f45 365
M 366     //backup vertex:        
58230a 367     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 368
369     //random move in a sphere with radius stepsize:
370     r=vesicle->stepsize*rn[0];
371     phi=rn[1]*2*M_PI;
372     costheta=2*rn[2]-1;
373     sintheta=sqrt(1-pow(costheta,2));
374     vtx->x=vtx->x+r*sintheta*cos(phi);
375     vtx->y=vtx->y+r*sintheta*sin(phi);
376     vtx->z=vtx->z+r*costheta;
377
378
379     //distance with neighbours check
380     for(i=0;i<vtx->bond_no;i++){
381         dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
382         if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
383             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
384             return TS_FAIL;
385         }
386     }
387
fe24d2 388 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 389 // Nucleus penetration check:
390     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
391         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
392         return TS_FAIL;
393     }
394
58230a 395
M 396     //self avoidance check with distant vertices
397     cellidx=vertex_self_avoidance(vesicle, vtx);
398     //check occupation number
399     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
400     if(retval==TS_FAIL){
401         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
402         return TS_FAIL;
403     } 
404
405     //backup bonds
406     for(i=0;i<vtx->bond_no;i++){
407         memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
408         vtx->bond[i]->bond_length=sqrt(dist[i]);
409         bond_vector(vtx->bond[i]);
b30f45 410     }
M 411
412     //backup neighboring vertices:
413     for(i=0;i<vtx->neigh_no;i++){
414         memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a 415     }
M 416     
417     //if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45 418     delta_energy=0;
M 419     
420     if(vtx->bond_no == 2){
421         vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
422         delta_energy += vtx->energy - backupvtx.energy;
58230a 423     }
M 424
b30f45 425     for(i=0;i<vtx->neigh_no;i++){
M 426         if(vtx->neigh[i]->bond_no == 2){
427             vtx->neigh[i]->energy = -(vtx->neigh[i]->bond[0]->x*vtx->neigh[i]->bond[1]->x + vtx->neigh[i]->bond[0]->y*vtx->neigh[i]->bond[1]->y + vtx->neigh[i]->bond[0]->z*vtx->neigh[i]->bond[1]->z)/vtx->neigh[i]->bond[0]->bond_length/vtx->neigh[i]->bond[1]->bond_length;
428             delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
429         }
58230a 430     }
M 431
b30f45 432     // poly->k is filament persistence length (in units l_min)
M 433     delta_energy *= poly->k;
58230a 434
M 435     if(delta_energy>=0){
436 #ifdef TS_DOUBLE_DOUBLE
437         if(exp(-delta_energy)< drand48() )
438 #endif
439 #ifdef TS_DOUBLE_FLOAT
440         if(expf(-delta_energy)< (ts_float)drand48())
441 #endif
442 #ifdef TS_DOUBLE_LONGDOUBLE
443         if(expl(-delta_energy)< (ts_ldouble)drand48())
444 #endif
445         {
446     //not accepted, reverting changes
447     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
b30f45 448     for(i=0;i<vtx->neigh_no;i++){
M 449         memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
450     }
58230a 451     for(i=0;i<vtx->bond_no;i++){
b30f45 452         vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
58230a 453     }
M 454
455     return TS_FAIL; 
456     }
457     }
458     
b30f45 459     
58230a 460 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M 461     if(vtx->cell!=vesicle->clist->cell[cellidx]){
462         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
463 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
464         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
465     }
466 //    if(oldcellidx);
467     //END MONTE CARLOOOOOOO
468     return TS_SUCCESS;
469 }