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