Added nucleus. Better initial configuration of filaments.

mihaf

2014-03-25 fe24d29b3c8684f08dccd01f5785aa48b7137322

Added nucleus. Better initial configuration of filaments.

 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,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);