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