trisurf-ng.git - Gitblit

Moved all global variables to separate file and defined extern keyword wher...

Samo Penic

2022-05-04 06afc729f9061c1cfa14c78728d61d518324c2f0

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