Trisurf Monte Carlo simulator
mihaf
2014-03-25 fe24d29b3c8684f08dccd01f5785aa48b7137322
Added nucleus. Better initial configuration of filaments.
5 files modified
61 ■■■■ changed files
src/initial_distribution.c 5 ●●●●● patch | view | raw | blame | history
src/poly.c 26 ●●●● patch | view | raw | blame | history
src/tape 10 ●●●● patch | view | raw | blame | history
src/timestep.c 6 ●●●●● patch | view | raw | blame | history
src/vertexmove.c 14 ●●●●● patch | view | raw | blame | history
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);
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);
            }
        }
    }
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!!!
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;
src/vertexmove.c
@@ -67,6 +67,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[0],sizeof(ts_vertex));
        return TS_FAIL;
    }
    //self avoidance check with distant vertices
     cellidx=vertex_self_avoidance(vesicle, vtx);
    //check occupation number
@@ -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);