trisurf-ng.git - Gitblit

Merge branch 'stretching-devel' into plane-new

Samo Penic

2019-02-08 4019ca8157136a68b6e43b31f166dc93acd81b4d

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