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