trisurf-ng.git - Gitblit

Samo Penic

2014-09-11 978e52dad1ade8ab82871846e52449a064c201ac

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