trisurf-ng.git - Gitblit

Before adding constant volume subrutines

Samo Penic

2014-04-28 2b14da2e8f012d7b1bf669af085ab75d87f378f1

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