From fe24d29b3c8684f08dccd01f5785aa48b7137322 Mon Sep 17 00:00:00 2001
From: mihaf <miha.fosnaric@gmail.com>
Date: Tue, 25 Mar 2014 12:49:06 +0000
Subject: [PATCH] Added nucleus. Better initial configuration of filaments.
---
src/timestep.c | 6 ++-
src/poly.c | 26 ++++++++++---
src/tape | 10 ++--
src/initial_distribution.c | 5 +-
src/vertexmove.c | 24 +++++++++--
5 files changed, 51 insertions(+), 20 deletions(-)
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index a93e2ee..0587d18 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -41,7 +41,9 @@
vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
// Nucleus:
- vesicle->R_nucleus=tape->R_nucleus;
+ vesicle->R_nucleus=tape->R_nucleus*tape->R_nucleus;
+
+ vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
//Initialize grafted polymers (brush):
vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle);
@@ -79,7 +81,6 @@
vesicle->nshell=tape->nshell;
vesicle->dmax=tape->dmax*tape->dmax; /* dmax^2 in the vesicle dmax variable */
- vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
vesicle->bending_rigidity=tape->xk0;
vtx_set_global_values(vesicle); /* make xk0 default value for every vertex */
ts_fprintf(stdout, "Tape setting: xk0=%e\n",tape->xk0);
diff --git a/src/poly.c b/src/poly.c
index 5cfc6e0..b2b13f9 100644
--- a/src/poly.c
+++ b/src/poly.c
@@ -58,7 +58,6 @@
ts_uint gvtxi;
ts_double xnorm,ynorm,znorm,normlength;
ts_double dphi,dh;
- ts_uint ji;
// Grafting polymers:
if (vlist!=NULL){
@@ -108,14 +107,29 @@
else
{
/* Make filaments inside the vesicle. Helix with radius... Dist. between poly vertices put to 1*/
- dphi = 2.0*asin(1.0/2.0/vesicle->R_nucleus)*1.001;
+ ts_double a,R,H,tantheta,h,r,phi,A0=1.2;
+
+ a = A0*(ts_double)vesicle->nshell;
+ R = A0*((ts_double)vesicle->nshell)/(2.0*sin(M_PI/5.0));
+ H = sqrt(a*a - R*R);
+ tantheta = sqrt(R*R - a*a/4.0)/H;
+
+ h = -H + sqrt(vesicle->clist->dmin_interspecies)*1.5;
+ r = (H-fabs(h))*tantheta - sqrt(vesicle->clist->dmin_interspecies)*1.5;
+ dphi = 2.0*asin(1.0/2.0/r)*1.001;
dh = dphi/2.0/M_PI*1.001;
+ phi=0.0;
for(i=0;i<poly_list->n;i++){
for (j=0;j<poly_list->poly[i]->vlist->n;j++){
- ji = j + i*poly_list->poly[i]->vlist->n;
- poly_list->poly[i]->vlist->vtx[j]->x = vesicle->R_nucleus*cos(ji*dphi);
- poly_list->poly[i]->vlist->vtx[j]->y = vesicle->R_nucleus*sin(ji*dphi);
- poly_list->poly[i]->vlist->vtx[j]->z = ji*dh - (dh*poly_list->n*poly_list->poly[i]->vlist->n/2.0);
+ h = h + dh;
+ r = (H-fabs(h))*tantheta - sqrt(vesicle->clist->dmin_interspecies)*1.5;
+ dphi = 2.0*asin(1.0/2.0/r)*1.001;
+ dh = dphi/2.0/M_PI*1.001;
+ phi+=dphi;
+ //ji = j + i*poly_list->poly[i]->vlist->n;
+ poly_list->poly[i]->vlist->vtx[j]->x = r*cos(phi);
+ poly_list->poly[i]->vlist->vtx[j]->y = r*sin(phi);
+ poly_list->poly[i]->vlist->vtx[j]->z = h;// ji*dh - (dh*poly_list->n*poly_list->poly[i]->vlist->n/2.0);
}
}
}
diff --git a/src/tape b/src/tape
index 8d589d0..e7e42e7 100644
--- a/src/tape
+++ b/src/tape
@@ -4,7 +4,7 @@
# dmax is the max. bond length (in units l_min)
dmax=1.7
# dmin_interspecies in the min. dist. between different vertex species (in units l_min)
-dmin_interspecies=1
+dmin_interspecies=1.2
# bending rigidity of the membrane (in units kT)
xk0=10.0
# max step size (in units l_min)
@@ -12,7 +12,7 @@
# Pressure calculations
# (pswitch=1: calc. p*dV energy contribution)
-pswitch = 1
+pswitch = 0
# pressure difference: p_inside - p_outside (in units kT/l_min^3):
pressure=-2.0
@@ -26,7 +26,7 @@
####### Filament (inside the vesicle) definitions ###########
# nfil is a number of filaments inside the vesicle
-nfil=1
+nfil=2
# nfono is a number of monomers in each filament
nfono=300
# Persistence lenght of the filaments (in units l_min)
@@ -44,11 +44,11 @@
####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
-mcsweeps=100
+mcsweeps=1000
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=0
#how many records do you want on the disk iteration are there in a run?
-iterations=20
+iterations=100
#shut up if we are using cluster!!!
diff --git a/src/timestep.c b/src/timestep.c
index 3369b9f..af152e4 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -11,7 +11,7 @@
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
ts_uint i, j;
- char filename[255];
+// char filename[255];
centermass(vesicle);
cell_occupation(vesicle);
@@ -26,8 +26,10 @@
dump_state(vesicle,i);
if(i>=inititer){
write_vertex_xml_file(vesicle,i-inititer);
+ write_master_xml_file("test.pvd");
+
// sprintf(filename,"timestep-%05d.pov",i-inititer);
- write_pov_file(vesicle,filename);
+ // write_pov_file(vesicle,filename);
}
}
return TS_SUCCESS;
diff --git a/src/vertexmove.c b/src/vertexmove.c
index da8e70f..9139ae2 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -67,11 +67,18 @@
}
}
- //self avoidance check with distant vertices
- cellidx=vertex_self_avoidance(vesicle, vtx);
- //check occupation number
- retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
-
+// TODO: Maybe faster if checks only nucleus-neighboring cells
+// Nucleus penetration check:
+ if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ return TS_FAIL;
+ }
+
+//self avoidance check with distant vertices
+ cellidx=vertex_self_avoidance(vesicle, vtx);
+ //check occupation number
+ retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
+
if(retval==TS_FAIL){
vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
return TS_FAIL;
@@ -291,6 +298,13 @@
}
}
+// TODO: Maybe faster if checks only nucleus-neighboring cells
+// Nucleus penetration check:
+ if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
+ vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
+ return TS_FAIL;
+ }
+
//self avoidance check with distant vertices
cellidx=vertex_self_avoidance(vesicle, vtx);
--
Gitblit v1.9.3