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