trisurf-ng.git - Gitblit

Samo Penic

2018-12-09 88bdd70987e76a58ea0fd917f63aa0c682848116

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