Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2018-06-21 620ba72765eda244ca3f3806b8ad709929e74b3c
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 */
168             //fprintf(stderr,"DVOL=%1.16e\n",dvol);
169         if(fabs(vesicle->volume+dvol-V0)>epsvol){
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]); 
176                     //fprintf(stderr,"fajlam!\n");
177                     return TS_FAIL;
178         }
179
180     } else
fe5069 181 //    vesicle_volume(vesicle);
SP 182 //    fprintf(stderr,"Volume before=%1.16e\n", vesicle->volume);
fbcbdf 183    if(vesicle->tape->constvolswitch == 1){
fe5069 184         retval=constvolume(vesicle, vtx, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
fbcbdf 185         if(retval==TS_FAIL){ // if we couldn't move the vertex to assure constant volume
SP 186             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
187             for(i=0;i<vtx->neigh_no;i++){
188                 vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
189             }
190             for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
fe5069 191  //           fprintf(stderr,"fajlam!\n");
fbcbdf 192             return TS_FAIL;
SP 193         }
fe5069 194 //    vesicle_volume(vesicle);
SP 195 //    fprintf(stderr,"Volume after=%1.16e\n", vesicle->volume);
196 //    fprintf(stderr,"Volume after-dvol=%1.16e\n", vesicle->volume-dvol);
197 //    fprintf(stderr,"Denergy before=%e\n",delta_energy);
198     
fbcbdf 199     delta_energy+=delta_energy_cv;
fe5069 200 //    fprintf(stderr,"Denergy after=%e\n",delta_energy);
fbcbdf 201     }
250de4 202 /* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */
SP 203     delta_energy+=direct_force_energy(vesicle,vtx,backupvtx);
7ec6fb 204
SP 205     //stretching energy 2 of 3
206     if(vesicle->tape->stretchswitch==1){
620ba7 207         for(i=0;i<vtx->tristar_no;i++){
SP 208             temp_area+=vtx->tristar[i]->area; 
7ec6fb 209             }
620ba7 210         dstretchenergy+=stretchenergy2(temp_area,vesicle->tlist->a0*vesicle->tlist->n);
SP 211         dstretchenergy=0.5*vesicle->tape->xkA0*dstretchenergy;
212         
7ec6fb 213     }
SP 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 */
314f2d 224 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 225     //MONTE CARLOOOOOOOO
e5858f 226 //    if(vtx->c!=0.0) printf("DE=%f\n",delta_energy);
aec47d 227     if(delta_energy>=0){
SP 228 #ifdef TS_DOUBLE_DOUBLE
3de289 229         if(exp(-delta_energy)< drand48())
aec47d 230 #endif
SP 231 #ifdef TS_DOUBLE_FLOAT
232         if(expf(-delta_energy)< (ts_float)drand48())
233 #endif
234 #ifdef TS_DOUBLE_LONGDOUBLE
235         if(expl(-delta_energy)< (ts_ldouble)drand48())
236 #endif
237     {
238     //not accepted, reverting changes
fbcbdf 239   //  fprintf(stderr,"MC failed\n");
dcd350 240     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 241     for(i=0;i<vtx->neigh_no;i++){
a63f17 242         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 243     }
SP 244     
aec47d 245     //update the normals of triangles that share bead i.
dcd350 246    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
620ba7 247 /*    //stretching energy 3 of 3
7ec6fb 248     if(vesicle->tape->stretchswitch==1){
SP 249         for(i=0;i<vtx->tristar_no;i++){ 
250             stretchenergy(vesicle,vtx->tristar[i]);
251             }
252     }
620ba7 253 */
fe5069 254 //    fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
43c042 255     if(vesicle->tape->constvolswitch == 1){
958e0e 256         constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
43c042 257     }
fe5069 258 //    fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
dd5aca 259 //    vesicle_volume(vesicle);
SP 260 //    fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
aec47d 261     return TS_FAIL; 
SP 262     }
263 }
2b14da 264     //accepted    
fbcbdf 265  //   fprintf(stderr,"MC accepted\n");
a63f17 266 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
SP 267     if(vtx->cell!=vesicle->clist->cell[cellidx]){
268         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
269 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
270         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
271         
272     }
2b14da 273
1121fa 274     if(vesicle->tape->constvolswitch == 2){
SP 275     vesicle->volume+=dvol;
276     } else
43c042 277     if(vesicle->tape->constvolswitch == 1){
fbcbdf 278         constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
43c042 279     }
c0ae90 280
SP 281     if(vesicle->tape->constareaswitch==2){
282         vesicle->area+=darea;
283     }
620ba7 284     //stretching energy 3 of 3
SP 285     if(vesicle->tape->stretchswitch==1){
286         vesicle->area=temp_area;
287     }
a63f17 288 //    if(oldcellidx);
aec47d 289     //END MONTE CARLOOOOOOO
dd5aca 290 //    vesicle_volume(vesicle);
SP 291 //    fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
aec47d 292     return TS_SUCCESS;
SP 293 }
294
fedf2b 295
M 296 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
297     ts_uint i;
298     ts_bool retval; 
299     ts_uint cellidx; 
304510 300 //    ts_double delta_energy;
fedf2b 301     ts_double costheta,sintheta,phi,r;
304510 302     ts_double dist;
fedf2b 303     //This will hold all the information of vtx and its neighbours
M 304     ts_vertex backupvtx;
304510 305 //    ts_bond backupbond[2];
fedf2b 306     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 307
308     //random move in a sphere with radius stepsize:
309     r=vesicle->stepsize*rn[0];
310     phi=rn[1]*2*M_PI;
311     costheta=2*rn[2]-1;
312     sintheta=sqrt(1-pow(costheta,2));
313     vtx->x=vtx->x+r*sintheta*cos(phi);
314     vtx->y=vtx->y+r*sintheta*sin(phi);
315     vtx->z=vtx->z+r*costheta;
316
317
318     //distance with neighbours check
304510 319     for(i=0;i<vtx->neigh_no;i++){
M 320         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
321         if(dist<1.0 || dist>vesicle->dmax) {
322             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
323             return TS_FAIL;
324         }
325     }
326
327 // Distance with grafted vesicle-vertex check:    
328     if(vtx==poly->vlist->vtx[0]){
329         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
330         if(dist<1.0 || dist>vesicle->dmax) {
331         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
332         return TS_FAIL;
333         }
334     }
335
fedf2b 336
M 337     //self avoidance check with distant vertices
338     cellidx=vertex_self_avoidance(vesicle, vtx);
339     //check occupation number
340     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
341     
342     if(retval==TS_FAIL){
343         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
344         return TS_FAIL;
345     } 
346
347
348     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 349 /* Energy ignored for now!
fedf2b 350     delta_energy=0;
M 351     for(i=0;i<vtx->bond_no;i++){
352         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
353
354         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
355         bond_energy(vtx->bond[i],poly);
356         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
357     }
358
359     if(vtx==poly->vlist->vtx[0]){
360         delta_energy+=
361             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
362             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
363         
364     }
365
366
367     if(delta_energy>=0){
368 #ifdef TS_DOUBLE_DOUBLE
369         if(exp(-delta_energy)< drand48() )
370 #endif
371 #ifdef TS_DOUBLE_FLOAT
372         if(expf(-delta_energy)< (ts_float)drand48())
373 #endif
374 #ifdef TS_DOUBLE_LONGDOUBLE
375         if(expl(-delta_energy)< (ts_ldouble)drand48())
376 #endif
377         {
378     //not accepted, reverting changes
379     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
380     for(i=0;i<vtx->bond_no;i++){
381     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
382     }
383
384     return TS_FAIL; 
385     }
386     }
304510 387 */
fedf2b 388         
M 389 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
390     if(vtx->cell!=vesicle->clist->cell[cellidx]){
391         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
392 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
393         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
394     }
395 //    if(oldcellidx);
396     //END MONTE CARLOOOOOOO
397     return TS_SUCCESS;
398 }
58230a 399
M 400
401
402
403 ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
404     ts_uint i;
405     ts_bool retval; 
406     ts_uint cellidx; 
b30f45 407     ts_double delta_energy;
58230a 408     ts_double costheta,sintheta,phi,r;
M 409     ts_double dist[2];
410     //This will hold all the information of vtx and its neighbours
b30f45 411     ts_vertex backupvtx,backupneigh[2];
58230a 412     ts_bond backupbond[2];
b30f45 413
M 414     //backup vertex:        
58230a 415     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 416
417     //random move in a sphere with radius stepsize:
418     r=vesicle->stepsize*rn[0];
419     phi=rn[1]*2*M_PI;
420     costheta=2*rn[2]-1;
421     sintheta=sqrt(1-pow(costheta,2));
422     vtx->x=vtx->x+r*sintheta*cos(phi);
423     vtx->y=vtx->y+r*sintheta*sin(phi);
424     vtx->z=vtx->z+r*costheta;
425
426
427     //distance with neighbours check
428     for(i=0;i<vtx->bond_no;i++){
429         dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
430         if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
431             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
432             return TS_FAIL;
433         }
434     }
435
fe24d2 436 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 437 // Nucleus penetration check:
438     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
439         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
440         return TS_FAIL;
441     }
442
58230a 443
M 444     //self avoidance check with distant vertices
445     cellidx=vertex_self_avoidance(vesicle, vtx);
446     //check occupation number
447     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
448     if(retval==TS_FAIL){
449         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
450         return TS_FAIL;
451     } 
452
453     //backup bonds
454     for(i=0;i<vtx->bond_no;i++){
455         memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
456         vtx->bond[i]->bond_length=sqrt(dist[i]);
457         bond_vector(vtx->bond[i]);
b30f45 458     }
M 459
460     //backup neighboring vertices:
461     for(i=0;i<vtx->neigh_no;i++){
462         memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a 463     }
M 464     
465     //if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45 466     delta_energy=0;
M 467     
468     if(vtx->bond_no == 2){
469         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;
470         delta_energy += vtx->energy - backupvtx.energy;
58230a 471     }
M 472
b30f45 473     for(i=0;i<vtx->neigh_no;i++){
M 474         if(vtx->neigh[i]->bond_no == 2){
475             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;
476             delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
477         }
58230a 478     }
M 479
b30f45 480     // poly->k is filament persistence length (in units l_min)
M 481     delta_energy *= poly->k;
58230a 482
M 483     if(delta_energy>=0){
484 #ifdef TS_DOUBLE_DOUBLE
485         if(exp(-delta_energy)< drand48() )
486 #endif
487 #ifdef TS_DOUBLE_FLOAT
488         if(expf(-delta_energy)< (ts_float)drand48())
489 #endif
490 #ifdef TS_DOUBLE_LONGDOUBLE
491         if(expl(-delta_energy)< (ts_ldouble)drand48())
492 #endif
493         {
494     //not accepted, reverting changes
495     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
b30f45 496     for(i=0;i<vtx->neigh_no;i++){
M 497         memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
498     }
58230a 499     for(i=0;i<vtx->bond_no;i++){
b30f45 500         vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
58230a 501     }
M 502
503     return TS_FAIL; 
504     }
505     }
506     
b30f45 507     
58230a 508 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M 509     if(vtx->cell!=vesicle->clist->cell[cellidx]){
510         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
511 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
512         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
513     }
514 //    if(oldcellidx);
515     //END MONTE CARLOOOOOOO
516     return TS_SUCCESS;
517 }