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