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