trisurf-ng.git - Gitblit

Samo Penic

2014-04-28 43c0420ae7c32fbcd3d784803bdbcf56b25b69e4

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;
414b8a	23	ts_double delta_energy,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
43c042	26	ts_vertex backupvtx[20], constvol_vtx_moved, constvol_vtx_backup;
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
43c042	95	if(vesicle->pswitch == 1){
414b8a	96	for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
M	97	};
a63f17	98
aec47d	99	delta_energy=0;
43c042	100
SP	101	if(vesicle->tape->constvolswitch == 1){
	102	retval=constvolume(vesicle, vtx, dvol, &delta_energy, constvol_vtx_moved,constvol_vtx_backup);
	103	if(retval==TS_FAIL){ // if we couldn't move the vertex to assure constant volume
	104	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	105	for(i=0;i<vtx->neigh_no;i++){
	106	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	107	}
	108	return TS_FAIL;
	109	}
	110	}
	111
aec47d	112	//update the normals of triangles that share bead i.
8f6a69	113	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17	114	oenergy=vtx->energy;
aec47d	115	energy_vertex(vtx);
a63f17	116	delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d	117	//the same is done for neighbouring vertices
8f6a69	118	for(i=0;i<vtx->neigh_no;i++){
SP	119	oenergy=vtx->neigh[i]->energy;
	120	energy_vertex(vtx->neigh[i]);
	121	delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d	122	}
414b8a	123
43c042	124	if(vesicle->pswitch == 1){
414b8a	125	for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
43c042	126	delta_energy-=vesicle->pressure*dvol;
414b8a	127	};
43c042	128
414b8a	129
304510	130	/* No poly-bond energy for now!
fedf2b	131	if(vtx->grafted_poly!=NULL){
M	132	delta_energy+=
	133	(pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
	134	pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
	135	}
304510	136	*/
314f2d	137	// fprintf(stderr, "DE=%f\n",delta_energy);
aec47d	138	//MONTE CARLOOOOOOOO
SP	139	if(delta_energy>=0){
	140	#ifdef TS_DOUBLE_DOUBLE
	141	if(exp(-delta_energy)< drand48() )
	142	#endif
	143	#ifdef TS_DOUBLE_FLOAT
	144	if(expf(-delta_energy)< (ts_float)drand48())
	145	#endif
	146	#ifdef TS_DOUBLE_LONGDOUBLE
	147	if(expl(-delta_energy)< (ts_ldouble)drand48())
	148	#endif
	149	{
	150	//not accepted, reverting changes
dcd350	151	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
1ad6d1	152	for(i=0;i<vtx->neigh_no;i++){
a63f17	153	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1	154	}
SP	155
aec47d	156	//update the normals of triangles that share bead i.
dcd350	157	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1	158
43c042	159	if(vesicle->tape->constvolswitch == 1){
SP	160	ts_bool constvolumerestore(constvol_vtx_backup);
	161	}
	162
aec47d	163	return TS_FAIL;
SP	164	}
	165	}
2b14da	166	//accepted
a63f17	167	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
SP	168	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	169	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	170	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	171	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
	172
	173	}
2b14da	174
43c042	175	if(vesicle->tape->constvolswitch == 1){
SP	176	ts_bool constvolumeaccept(constvol_vtx_backup);
	177	}
a63f17	178	// if(oldcellidx);
aec47d	179	//END MONTE CARLOOOOOOO
SP	180	return TS_SUCCESS;
	181	}
	182
fedf2b	183
M	184	ts_bool single_poly_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	185	ts_uint i;
	186	ts_bool retval;
	187	ts_uint cellidx;
304510	188	// ts_double delta_energy;
fedf2b	189	ts_double costheta,sintheta,phi,r;
304510	190	ts_double dist;
fedf2b	191	//This will hold all the information of vtx and its neighbours
M	192	ts_vertex backupvtx;
304510	193	// ts_bond backupbond[2];
fedf2b	194	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	195
	196	//random move in a sphere with radius stepsize:
	197	r=vesicle->stepsize*rn[0];
	198	phi=rn[1]2M_PI;
	199	costheta=2*rn[2]-1;
	200	sintheta=sqrt(1-pow(costheta,2));
	201	vtx->x=vtx->x+rsinthetacos(phi);
	202	vtx->y=vtx->y+rsinthetasin(phi);
	203	vtx->z=vtx->z+r*costheta;
	204
	205
	206	//distance with neighbours check
304510	207	for(i=0;i<vtx->neigh_no;i++){
M	208	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	209	if(dist<1.0 \|\| dist>vesicle->dmax) {
	210	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	211	return TS_FAIL;
	212	}
	213	}
	214
	215	// Distance with grafted vesicle-vertex check:
	216	if(vtx==poly->vlist->vtx[0]){
	217	dist=vtx_distance_sq(vtx,poly->grafted_vtx);
	218	if(dist<1.0 \|\| dist>vesicle->dmax) {
	219	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	220	return TS_FAIL;
	221	}
	222	}
	223
fedf2b	224
M	225	//self avoidance check with distant vertices
	226	cellidx=vertex_self_avoidance(vesicle, vtx);
	227	//check occupation number
	228	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	229
	230	if(retval==TS_FAIL){
	231	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	232	return TS_FAIL;
	233	}
	234
	235
	236	//if all the tests are successful, then energy for vtx and neighbours is calculated
304510	237	/* Energy ignored for now!
fedf2b	238	delta_energy=0;
M	239	for(i=0;i<vtx->bond_no;i++){
	240	memcpy((void )&backupbond[i],(void )vtx->bond[i],sizeof(ts_bond));
	241
	242	vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
	243	bond_energy(vtx->bond[i],poly);
	244	delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
	245	}
	246
	247	if(vtx==poly->vlist->vtx[0]){
	248	delta_energy+=
	249	(pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
	250	pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
	251
	252	}
	253
	254
	255	if(delta_energy>=0){
	256	#ifdef TS_DOUBLE_DOUBLE
	257	if(exp(-delta_energy)< drand48() )
	258	#endif
	259	#ifdef TS_DOUBLE_FLOAT
	260	if(expf(-delta_energy)< (ts_float)drand48())
	261	#endif
	262	#ifdef TS_DOUBLE_LONGDOUBLE
	263	if(expl(-delta_energy)< (ts_ldouble)drand48())
	264	#endif
	265	{
	266	//not accepted, reverting changes
	267	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	268	for(i=0;i<vtx->bond_no;i++){
	269	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
	270	}
	271
	272	return TS_FAIL;
	273	}
	274	}
304510	275	*/
fedf2b	276
M	277	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	278	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	279	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	280	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	281	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	282	}
	283	// if(oldcellidx);
	284	//END MONTE CARLOOOOOOO
	285	return TS_SUCCESS;
	286	}
58230a	287
M	288
	289
	290
	291	ts_bool single_filament_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	292	ts_uint i;
	293	ts_bool retval;
	294	ts_uint cellidx;
b30f45	295	ts_double delta_energy;
58230a	296	ts_double costheta,sintheta,phi,r;
M	297	ts_double dist[2];
	298	//This will hold all the information of vtx and its neighbours
b30f45	299	ts_vertex backupvtx,backupneigh[2];
58230a	300	ts_bond backupbond[2];
b30f45	301
M	302	//backup vertex:
58230a	303	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	304
	305	//random move in a sphere with radius stepsize:
	306	r=vesicle->stepsize*rn[0];
	307	phi=rn[1]2M_PI;
	308	costheta=2*rn[2]-1;
	309	sintheta=sqrt(1-pow(costheta,2));
	310	vtx->x=vtx->x+rsinthetacos(phi);
	311	vtx->y=vtx->y+rsinthetasin(phi);
	312	vtx->z=vtx->z+r*costheta;
	313
	314
	315	//distance with neighbours check
	316	for(i=0;i<vtx->bond_no;i++){
	317	dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
	318	if(dist[i]<1.0 \|\| dist[i]>vesicle->dmax) {
	319	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	320	return TS_FAIL;
	321	}
	322	}
	323
fe24d2	324	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	325	// Nucleus penetration check:
	326	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	327	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	328	return TS_FAIL;
	329	}
	330
58230a	331
M	332	//self avoidance check with distant vertices
	333	cellidx=vertex_self_avoidance(vesicle, vtx);
	334	//check occupation number
	335	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	336	if(retval==TS_FAIL){
	337	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	338	return TS_FAIL;
	339	}
	340
	341	//backup bonds
	342	for(i=0;i<vtx->bond_no;i++){
	343	memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
	344	vtx->bond[i]->bond_length=sqrt(dist[i]);
	345	bond_vector(vtx->bond[i]);
b30f45	346	}
M	347
	348	//backup neighboring vertices:
	349	for(i=0;i<vtx->neigh_no;i++){
	350	memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a	351	}
M	352
	353	//if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45	354	delta_energy=0;
M	355
	356	if(vtx->bond_no == 2){
	357	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;
	358	delta_energy += vtx->energy - backupvtx.energy;
58230a	359	}
M	360
b30f45	361	for(i=0;i<vtx->neigh_no;i++){
M	362	if(vtx->neigh[i]->bond_no == 2){
	363	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;
	364	delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
	365	}
58230a	366	}
M	367
b30f45	368	// poly->k is filament persistence length (in units l_min)
M	369	delta_energy *= poly->k;
58230a	370
M	371	if(delta_energy>=0){
	372	#ifdef TS_DOUBLE_DOUBLE
	373	if(exp(-delta_energy)< drand48() )
	374	#endif
	375	#ifdef TS_DOUBLE_FLOAT
	376	if(expf(-delta_energy)< (ts_float)drand48())
	377	#endif
	378	#ifdef TS_DOUBLE_LONGDOUBLE
	379	if(expl(-delta_energy)< (ts_ldouble)drand48())
	380	#endif
	381	{
	382	//not accepted, reverting changes
	383	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
b30f45	384	for(i=0;i<vtx->neigh_no;i++){
M	385	memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
	386	}
58230a	387	for(i=0;i<vtx->bond_no;i++){
b30f45	388	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
58230a	389	}
M	390
	391	return TS_FAIL;
	392	}
	393	}
	394
b30f45	395
58230a	396	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M	397	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	398	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	399	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	400	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	401	}
	402	// if(oldcellidx);
	403	//END MONTE CARLOOOOOOO
	404	return TS_SUCCESS;
	405	}