trisurf-ng.git - Gitblit

mihaf

2014-03-26 bc9583f9f05c86aa21da0d1b9081660a64d9ce31

commit \| author \| age
aec47d	1	#include<stdlib.h>
SP	2	#include<math.h>
	3	#include "general.h"
	4	#include "vertex.h"
	5	#include "bond.h"
	6	#include "triangle.h"
	7	#include "vesicle.h"
	8	#include "energy.h"
	9	#include "timestep.h"
	10	#include "cell.h"
	11	//#include "io.h"
	12	#include<stdio.h>
	13	#include "vertexmove.h"
1ad6d1	14	#include <string.h>
aec47d	15
fedf2b	16	ts_bool single_verticle_timestep(ts_vesicle vesicle,ts_vertex vtx,ts_double *rn){
aec47d	17	ts_uint i;
SP	18	ts_double dist;
	19	ts_bool retval;
	20	ts_uint cellidx;
414b8a	21	ts_double delta_energy,oenergy,dvol=0.0;
ed31fe	22	ts_double costheta,sintheta,phi,r;
1ad6d1	23	//This will hold all the information of vtx and its neighbours
dcd350	24	ts_vertex backupvtx[20];
SP	25	memcpy((void )&backupvtx[0],(void )vtx,sizeof(ts_vertex));
a63f17	26
SP	27	//Some stupid tests for debugging cell occupation!
	28	/* cellidx=vertex_self_avoidance(vesicle, vtx);
	29	if(vesicle->clist->cell[cellidx]==vtx->cell){
	30	fprintf(stderr,"Idx match!\n");
	31	} else {
	32	fprintf(stderr,"***** Idx don't match!\n");
	33	fatal("ENding.",1);
	34	}
	35	*/
	36
352fad	37	//temporarly moving the vertex
672ae4	38	// vtx->x=vtx->x+vesicle->stepsize(2.0rn[0]-1.0);
SP	39	// vtx->y=vtx->y+vesicle->stepsize(2.0rn[1]-1.0);
	40	// vtx->z=vtx->z+vesicle->stepsize(2.0rn[2]-1.0);
	41
ed31fe	42	//random move in a sphere with radius stepsize:
M	43	r=vesicle->stepsize*rn[0];
	44	phi=rn[1]2M_PI;
	45	costheta=2*rn[2]-1;
	46	sintheta=sqrt(1-pow(costheta,2));
672ae4	47	vtx->x=vtx->x+rsinthetacos(phi);
SP	48	vtx->y=vtx->y+rsinthetasin(phi);
	49	vtx->z=vtx->z+r*costheta;
	50
	51
a63f17	52	//distance with neighbours check
8f6a69	53	for(i=0;i<vtx->neigh_no;i++){
352fad	54	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
8f6a69	55	if(dist<1.0 \|\| dist>vesicle->dmax) {
dcd350	56	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
8f6a69	57	return TS_FAIL;
SP	58	}
aec47d	59	}
304510	60
M	61	// Distance with grafted poly-vertex check:
	62	if(vtx->grafted_poly!=NULL){
	63	dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
	64	if(dist<1.0 \|\| dist>vesicle->dmax) {
	65	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	66	return TS_FAIL;
	67	}
	68	}
	69
fe24d2	70	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	71	// Nucleus penetration check:
	72	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	73	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	74	return TS_FAIL;
	75	}
	76
	77	//self avoidance check with distant vertices
	78	cellidx=vertex_self_avoidance(vesicle, vtx);
	79	//check occupation number
	80	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	81
aec47d	82	if(retval==TS_FAIL){
dcd350	83	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
aec47d	84	return TS_FAIL;
SP	85	}
1ad6d1	86
SP	87
352fad	88	//if all the tests are successful, then energy for vtx and neighbours is calculated
1ad6d1	89	for(i=0;i<vtx->neigh_no;i++){
dcd350	90	memcpy((void )&backupvtx[i+1],(void )vtx->neigh[i],sizeof(ts_vertex));
1ad6d1	91	}
aec47d	92
414b8a	93	if(vesicle->pswitch == 1){
M	94	for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
	95	};
a63f17	96
aec47d	97	delta_energy=0;
SP	98	//update the normals of triangles that share bead i.
8f6a69	99	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17	100	oenergy=vtx->energy;
aec47d	101	energy_vertex(vtx);
a63f17	102	delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d	103	//the same is done for neighbouring vertices
8f6a69	104	for(i=0;i<vtx->neigh_no;i++){
SP	105	oenergy=vtx->neigh[i]->energy;
	106	energy_vertex(vtx->neigh[i]);
	107	delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d	108	}
414b8a	109
M	110	if(vesicle->pswitch == 1){
	111	for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
	112	delta_energy-=vesicle->pressure*dvol;
	113	};
	114
304510	115	/* No poly-bond energy for now!
fedf2b	116	if(vtx->grafted_poly!=NULL){
M	117	delta_energy+=
	118	(pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
	119	pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
	120	}
304510	121	*/
314f2d	122	// fprintf(stderr, "DE=%f\n",delta_energy);
aec47d	123	//MONTE CARLOOOOOOOO
SP	124	if(delta_energy>=0){
	125	#ifdef TS_DOUBLE_DOUBLE
	126	if(exp(-delta_energy)< drand48() )
	127	#endif
	128	#ifdef TS_DOUBLE_FLOAT
	129	if(expf(-delta_energy)< (ts_float)drand48())
	130	#endif
	131	#ifdef TS_DOUBLE_LONGDOUBLE
	132	if(expl(-delta_energy)< (ts_ldouble)drand48())
	133	#endif
	134	{
	135	//not accepted, reverting changes
dcd350	136	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
1ad6d1	137	for(i=0;i<vtx->neigh_no;i++){
a63f17	138	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1	139	}
SP	140
aec47d	141	//update the normals of triangles that share bead i.
dcd350	142	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1	143
aec47d	144	return TS_FAIL;
SP	145	}
	146	}
a63f17	147
SP	148	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	149	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	150	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	151	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	152	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
	153
	154	}
	155	// if(oldcellidx);
aec47d	156	//END MONTE CARLOOOOOOO
SP	157	return TS_SUCCESS;
	158	}
	159
fedf2b	160
M	161	ts_bool single_poly_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	162	ts_uint i;
	163	ts_bool retval;
	164	ts_uint cellidx;
304510	165	// ts_double delta_energy;
fedf2b	166	ts_double costheta,sintheta,phi,r;
304510	167	ts_double dist;
fedf2b	168	//This will hold all the information of vtx and its neighbours
M	169	ts_vertex backupvtx;
304510	170	// ts_bond backupbond[2];
fedf2b	171	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	172
	173	//random move in a sphere with radius stepsize:
	174	r=vesicle->stepsize*rn[0];
	175	phi=rn[1]2M_PI;
	176	costheta=2*rn[2]-1;
	177	sintheta=sqrt(1-pow(costheta,2));
	178	vtx->x=vtx->x+rsinthetacos(phi);
	179	vtx->y=vtx->y+rsinthetasin(phi);
	180	vtx->z=vtx->z+r*costheta;
	181
	182
	183	//distance with neighbours check
304510	184	for(i=0;i<vtx->neigh_no;i++){
M	185	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	186	if(dist<1.0 \|\| dist>vesicle->dmax) {
	187	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	188	return TS_FAIL;
	189	}
	190	}
	191
	192	// Distance with grafted vesicle-vertex check:
	193	if(vtx==poly->vlist->vtx[0]){
	194	dist=vtx_distance_sq(vtx,poly->grafted_vtx);
	195	if(dist<1.0 \|\| dist>vesicle->dmax) {
	196	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	197	return TS_FAIL;
	198	}
	199	}
	200
fedf2b	201
M	202	//self avoidance check with distant vertices
	203	cellidx=vertex_self_avoidance(vesicle, vtx);
	204	//check occupation number
	205	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	206
	207	if(retval==TS_FAIL){
	208	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	209	return TS_FAIL;
	210	}
	211
	212
	213	//if all the tests are successful, then energy for vtx and neighbours is calculated
304510	214	/* Energy ignored for now!
fedf2b	215	delta_energy=0;
M	216	for(i=0;i<vtx->bond_no;i++){
	217	memcpy((void )&backupbond[i],(void )vtx->bond[i],sizeof(ts_bond));
	218
	219	vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
	220	bond_energy(vtx->bond[i],poly);
	221	delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
	222	}
	223
	224	if(vtx==poly->vlist->vtx[0]){
	225	delta_energy+=
	226	(pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
	227	pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
	228
	229	}
	230
	231
	232	if(delta_energy>=0){
	233	#ifdef TS_DOUBLE_DOUBLE
	234	if(exp(-delta_energy)< drand48() )
	235	#endif
	236	#ifdef TS_DOUBLE_FLOAT
	237	if(expf(-delta_energy)< (ts_float)drand48())
	238	#endif
	239	#ifdef TS_DOUBLE_LONGDOUBLE
	240	if(expl(-delta_energy)< (ts_ldouble)drand48())
	241	#endif
	242	{
	243	//not accepted, reverting changes
	244	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	245	for(i=0;i<vtx->bond_no;i++){
	246	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
	247	}
	248
	249	return TS_FAIL;
	250	}
	251	}
304510	252	*/
fedf2b	253
M	254	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	255	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	256	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	257	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	258	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	259	}
	260	// if(oldcellidx);
	261	//END MONTE CARLOOOOOOO
	262	return TS_SUCCESS;
	263	}
58230a	264
M	265
	266
	267
	268	ts_bool single_filament_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	269	ts_uint i;
	270	ts_bool retval;
	271	ts_uint cellidx;
b30f45	272	ts_double delta_energy;
58230a	273	ts_double costheta,sintheta,phi,r;
M	274	ts_double dist[2];
	275	//This will hold all the information of vtx and its neighbours
b30f45	276	ts_vertex backupvtx,backupneigh[2];
58230a	277	ts_bond backupbond[2];
b30f45	278
M	279	//backup vertex:
58230a	280	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	281
	282	//random move in a sphere with radius stepsize:
	283	r=vesicle->stepsize*rn[0];
	284	phi=rn[1]2M_PI;
	285	costheta=2*rn[2]-1;
	286	sintheta=sqrt(1-pow(costheta,2));
	287	vtx->x=vtx->x+rsinthetacos(phi);
	288	vtx->y=vtx->y+rsinthetasin(phi);
	289	vtx->z=vtx->z+r*costheta;
	290
	291
	292	//distance with neighbours check
	293	for(i=0;i<vtx->bond_no;i++){
	294	dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
	295	if(dist[i]<1.0 \|\| dist[i]>vesicle->dmax) {
	296	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	297	return TS_FAIL;
	298	}
	299	}
	300
fe24d2	301	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	302	// Nucleus penetration check:
	303	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	304	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	305	return TS_FAIL;
	306	}
	307
58230a	308
M	309	//self avoidance check with distant vertices
	310	cellidx=vertex_self_avoidance(vesicle, vtx);
	311	//check occupation number
	312	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	313
	314	if(retval==TS_FAIL){
	315	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	316	return TS_FAIL;
	317	}
	318
	319	//backup bonds
	320	for(i=0;i<vtx->bond_no;i++){
	321	memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
	322	vtx->bond[i]->bond_length=sqrt(dist[i]);
	323	bond_vector(vtx->bond[i]);
b30f45	324	}
M	325
	326	//backup neighboring vertices:
	327	for(i=0;i<vtx->neigh_no;i++){
	328	memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a	329	}
M	330
	331	//if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45	332	delta_energy=0;
M	333
	334	if(vtx->bond_no == 2){
	335	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;
	336	delta_energy += vtx->energy - backupvtx.energy;
58230a	337	}
M	338
b30f45	339	for(i=0;i<vtx->neigh_no;i++){
M	340	if(vtx->neigh[i]->bond_no == 2){
	341	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;
	342	delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
	343	}
58230a	344	}
M	345
b30f45	346	// poly->k is filament persistence length (in units l_min)
M	347	delta_energy *= poly->k;
58230a	348
M	349	if(delta_energy>=0){
	350	#ifdef TS_DOUBLE_DOUBLE
	351	if(exp(-delta_energy)< drand48() )
	352	#endif
	353	#ifdef TS_DOUBLE_FLOAT
	354	if(expf(-delta_energy)< (ts_float)drand48())
	355	#endif
	356	#ifdef TS_DOUBLE_LONGDOUBLE
	357	if(expl(-delta_energy)< (ts_ldouble)drand48())
	358	#endif
	359	{
	360	//not accepted, reverting changes
	361	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
b30f45	362	for(i=0;i<vtx->neigh_no;i++){
M	363	memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
	364	}
58230a	365	for(i=0;i<vtx->bond_no;i++){
b30f45	366	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
58230a	367	}
M	368
	369	return TS_FAIL;
	370	}
	371	}
	372
b30f45	373
58230a	374	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M	375	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	376	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	377	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	378	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	379	}
	380	// if(oldcellidx);
	381	//END MONTE CARLOOOOOOO
	382	return TS_SUCCESS;
	383	}