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