trisurf-ng.git - Gitblit

Fixes in a0 calculation and calculation of stretching energy in vertexmove

Samo Penic

2018-05-21 699ac44864e75549fb3d238d36f0bc1dd645ddac

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.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	}
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){
699ac4	118	for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
7ec6fb	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]);
699ac4	202	dstretchenergy+=vtx->tristar[i]->energy;
7ec6fb	203	}
SP	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]2M_PI;
	298	costheta=2*rn[2]-1;
	299	sintheta=sqrt(1-pow(costheta,2));
	300	vtx->x=vtx->x+rsinthetacos(phi);
	301	vtx->y=vtx->y+rsinthetasin(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]2M_PI;
	407	costheta=2*rn[2]-1;
	408	sintheta=sqrt(1-pow(costheta,2));
	409	vtx->x=vtx->x+rsinthetacos(phi);
	410	vtx->y=vtx->y+rsinthetasin(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->xvtx->x + vtx->yvtx->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]->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;
	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]->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;
	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	}