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