trisurf-ng.git - Gitblit

Samo Penic

2014-03-05 719c9febac2eaff9483fda487b57684afbb59bb2

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;
	21	ts_double delta_energy,oenergy;
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
6fa0c9	53	for(i=0;i<vtx->neigh->n;i++){
SP	54	dist=vtx_distance_sq(vtx,vtx->neigh->vtx[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
aec47d	70	//self avoidance check with distant vertices
352fad	71	cellidx=vertex_self_avoidance(vesicle, vtx);
aec47d	72	//check occupation number
a63f17	73	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
SP	74
aec47d	75	if(retval==TS_FAIL){
dcd350	76	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
aec47d	77	return TS_FAIL;
SP	78	}
1ad6d1	79
SP	80
352fad	81	//if all the tests are successful, then energy for vtx and neighbours is calculated
6fa0c9	82	for(i=0;i<vtx->neigh->n;i++){
SP	83	memcpy((void )&backupvtx[i+1],(void )vtx->neigh->vtx[i],sizeof(ts_vertex));
1ad6d1	84	}
aec47d	85
a63f17	86
SP	87
aec47d	88	delta_energy=0;
SP	89	//update the normals of triangles that share bead i.
8f6a69	90	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17	91	oenergy=vtx->energy;
aec47d	92	energy_vertex(vtx);
a63f17	93	delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d	94	//the same is done for neighbouring vertices
6fa0c9	95	for(i=0;i<vtx->neigh->n;i++){
SP	96	oenergy=vtx->neigh->vtx[i]->energy;
	97	energy_vertex(vtx->neigh->vtx[i]);
	98	delta_energy+=vtx->neigh->vtx[i]->xk*(vtx->neigh->vtx[i]->energy-oenergy);
aec47d	99	}
304510	100	/* No poly-bond energy for now!
fedf2b	101	if(vtx->grafted_poly!=NULL){
M	102	delta_energy+=
	103	(pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
	104	pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
	105	}
304510	106	*/
314f2d	107	// fprintf(stderr, "DE=%f\n",delta_energy);
aec47d	108	//MONTE CARLOOOOOOOO
SP	109	if(delta_energy>=0){
	110	#ifdef TS_DOUBLE_DOUBLE
	111	if(exp(-delta_energy)< drand48() )
	112	#endif
	113	#ifdef TS_DOUBLE_FLOAT
	114	if(expf(-delta_energy)< (ts_float)drand48())
	115	#endif
	116	#ifdef TS_DOUBLE_LONGDOUBLE
	117	if(expl(-delta_energy)< (ts_ldouble)drand48())
	118	#endif
	119	{
	120	//not accepted, reverting changes
dcd350	121	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
6fa0c9	122	for(i=0;i<vtx->neigh->n;i++){
SP	123	vtx->neigh->vtx[i]=memcpy((void )vtx->neigh->vtx[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1	124	}
SP	125
aec47d	126	//update the normals of triangles that share bead i.
dcd350	127	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1	128
aec47d	129	return TS_FAIL;
SP	130	}
	131	}
a63f17	132
SP	133	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	134	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	135	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	136	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	137	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
	138
	139	}
	140	// if(oldcellidx);
aec47d	141	//END MONTE CARLOOOOOOO
SP	142	return TS_SUCCESS;
	143	}
	144
fedf2b	145
M	146	ts_bool single_poly_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	147	ts_uint i;
	148	ts_bool retval;
	149	ts_uint cellidx;
304510	150	// ts_double delta_energy;
fedf2b	151	ts_double costheta,sintheta,phi,r;
304510	152	ts_double dist;
fedf2b	153	//This will hold all the information of vtx and its neighbours
M	154	ts_vertex backupvtx;
304510	155	// ts_bond backupbond[2];
fedf2b	156	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	157
	158	//random move in a sphere with radius stepsize:
	159	r=vesicle->stepsize*rn[0];
	160	phi=rn[1]2M_PI;
	161	costheta=2*rn[2]-1;
	162	sintheta=sqrt(1-pow(costheta,2));
	163	vtx->x=vtx->x+rsinthetacos(phi);
	164	vtx->y=vtx->y+rsinthetasin(phi);
	165	vtx->z=vtx->z+r*costheta;
	166
	167
	168	//distance with neighbours check
304510	169	for(i=0;i<vtx->neigh_no;i++){
M	170	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	171	if(dist<1.0 \|\| dist>vesicle->dmax) {
	172	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	173	return TS_FAIL;
	174	}
	175	}
	176
	177	// Distance with grafted vesicle-vertex check:
	178	if(vtx==poly->vlist->vtx[0]){
	179	dist=vtx_distance_sq(vtx,poly->grafted_vtx);
	180	if(dist<1.0 \|\| dist>vesicle->dmax) {
	181	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	182	return TS_FAIL;
	183	}
	184	}
	185
fedf2b	186
M	187	//self avoidance check with distant vertices
	188	cellidx=vertex_self_avoidance(vesicle, vtx);
	189	//check occupation number
	190	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	191
	192	if(retval==TS_FAIL){
	193	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	194	return TS_FAIL;
	195	}
	196
	197
	198	//if all the tests are successful, then energy for vtx and neighbours is calculated
304510	199	/* Energy ignored for now!
fedf2b	200	delta_energy=0;
M	201	for(i=0;i<vtx->bond_no;i++){
	202	memcpy((void )&backupbond[i],(void )vtx->bond[i],sizeof(ts_bond));
	203
	204	vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
	205	bond_energy(vtx->bond[i],poly);
	206	delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
	207	}
	208
	209	if(vtx==poly->vlist->vtx[0]){
	210	delta_energy+=
	211	(pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
	212	pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
	213
	214	}
	215
	216
	217	if(delta_energy>=0){
	218	#ifdef TS_DOUBLE_DOUBLE
	219	if(exp(-delta_energy)< drand48() )
	220	#endif
	221	#ifdef TS_DOUBLE_FLOAT
	222	if(expf(-delta_energy)< (ts_float)drand48())
	223	#endif
	224	#ifdef TS_DOUBLE_LONGDOUBLE
	225	if(expl(-delta_energy)< (ts_ldouble)drand48())
	226	#endif
	227	{
	228	//not accepted, reverting changes
	229	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	230	for(i=0;i<vtx->bond_no;i++){
	231	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
	232	}
	233
	234	return TS_FAIL;
	235	}
	236	}
304510	237	*/
fedf2b	238
M	239	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	240	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	241	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	242	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	243	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	244	}
	245	// if(oldcellidx);
	246	//END MONTE CARLOOOOOOO
	247	return TS_SUCCESS;
	248	}