Adding ending energy to filaments. In progress. Continue in single_filament...

mihaf

2014-03-18 58230a2591414fb38b9ec8d3a76439b290cb0a6f

Adding ending energy to filaments. In progress. Continue in single_filament_vertex_move

 src/bond.c

@@ -35,6 +35,17 @@
    return blist->bond[blist->n-1];
}


ts_bool bond_vector(ts_bond *bond){
	
    bond->x=bond->vtx1->x-bond->vtx2->x;
    bond->y=bond->vtx1->y-bond->vtx2->y;
    bond->z=bond->vtx1->z-bond->vtx2->z;

    return TS_SUCCESS;	
}


ts_bool bond_list_free(ts_bond_list *blist){
    ts_uint i;
    for(i=0;i<blist->n;i++){

 src/bond.h

@@ -20,6 +20,7 @@
 */
ts_bond *bond_add(ts_bond_list *blist, ts_vertex *vtx1, ts_vertex *vtx2);

ts_bool bond_vector(ts_bond *bond);
ts_bool bond_list_free(ts_bond_list *blist);



 src/frame.c

@@ -22,7 +22,7 @@
        vtx[i]->y-=vesicle->cm[1];
        vtx[i]->z-=vesicle->cm[2]; 
    } 

//move polymers for the same vector as we moved vesicle
    for(i=0;i<vesicle->poly_list->n;i++){
        for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            vesicle->poly_list->poly[i]->vlist->vtx[j]->x-=vesicle->cm[0];
@@ -30,22 +30,24 @@
            vesicle->poly_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
        }
    }
//move filaments for the same vector as we moved vesicle
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            vesicle->filament_list->poly[i]->vlist->vtx[j]->x-=vesicle->cm[0];
            vesicle->filament_list->poly[i]->vlist->vtx[j]->y-=vesicle->cm[1];
            vesicle->filament_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
        }
    }


    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
    vesicle->cm[0]=0.0;
    vesicle->cm[1]=0.0;
    vesicle->cm[2]=0.0;

    return TS_SUCCESS;
}

ts_bool cell_occupation(ts_vesicle *vesicle){
    ts_uint i,j,cellidx, n=vesicle->vlist->n;
    //ts_double shift;
    //ts_double dcell;
    //shift=(ts_double) vesicle->clist->ncmax[0]/2;
    //dcell=1.0/(1.0 + vesicle->stepsize);
    //`fprintf(stderr, "Bil sem tu\n"); 

    cell_list_cell_occupation_clear(vesicle->clist);
    for(i=0;i<n;i++){
@@ -56,17 +58,21 @@
    cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);
    }

//Add all polymers to cells
    for(i=0;i<vesicle->poly_list->n;i++){
    for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
        cellidx=vertex_self_avoidance(vesicle, vesicle->poly_list->poly[i]->vlist->vtx[j]);
        cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->poly_list->poly[i]->vlist->vtx[j]);
    }
    }
//Add all filaments to cells
     for(i=0;i<vesicle->filament_list->n;i++){
    for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
        cellidx=vertex_self_avoidance(vesicle, vesicle->filament_list->poly[i]->vlist->vtx[j]);
        cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->filament_list->poly[i]->vlist->vtx[j]);
    }
    }
   

    

    //fprintf(stderr, "Bil sem tu\n"); 
    //if(dcell);
    //if(shift);
    return TS_SUCCESS;
}

 src/general.h

@@ -166,10 +166,11 @@
        ts_uint idx;
    ts_vertex *vtx1;
    ts_vertex *vtx2;
    ts_double bond_length;
    ts_double bond_length_dual;
    ts_bool tainted;
        ts_double bond_length;
        ts_double bond_length_dual;
    ts_bool tainted; //TODO: remove
    ts_double energy;
    ts_double x,y,z;
};
typedef struct ts_bond ts_bond;


 src/initial_distribution.c

@@ -38,9 +38,24 @@
ts_vesicle *create_vesicle_from_tape(ts_tape *tape){
    ts_vesicle *vesicle;
    vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
    // Nucleus:
    vesicle->R_nucleus=tape->R_nucleus;
    //Initialize grafted polymers (brush):
    vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle);
    vesicle->spring_constant=tape->kspring;
    poly_assign_spring_const(vesicle);
    //Initialize filaments (polymers inside the vesicle):
    vesicle->filament_list=init_poly_list(tape->nfil,tape->nfono, NULL, vesicle);
ts_uint i,j;
    for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){

    fprintf(stderr,"%e %e %e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
    }
    }
//    vesicle->spring_constant=tape->kspring;
//    poly_assign_spring_const(vesicle);

    
    vesicle->nshell=tape->nshell;
    vesicle->dmax=tape->dmax*tape->dmax; /* dmax^2 in the vesicle dmax variable */

 src/io.c

@@ -664,8 +664,8 @@
    /* Here comes header of the file */

    //find number of extra vtxs and bonds of polymeres
    ts_uint monono=0, polyno=0, poly_idx=0;
    ts_bool poly=0;
    ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
    ts_bool poly=0, fil=0;
    if(vesicle->poly_list!=NULL){
        if(vesicle->poly_list->poly[0]!=NULL){
        polyno=vesicle->poly_list->n;
@@ -673,8 +673,17 @@
        poly=1;
        }
    }

    if(vesicle->filament_list!=NULL){
        if(vesicle->filament_list->poly[0]!=NULL){
        filno=vesicle->filament_list->n;
        fonono=vesicle->filament_list->poly[0]->vlist->n;
        fil=1;
        }
    }

    fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n");
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno, blist->n+monono*polyno);
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1));
    fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
       for(i=0;i<vlist->n;i++){
        fprintf(fh,"%u ",vtx[i]->idx);
@@ -684,6 +693,16 @@
        poly_idx=vlist->n;
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
                fprintf(fh,"%u ", poly_idx);
            }
        }
    }
    //filaments
    if(fil){
        poly_idx=vlist->n+monono*polyno;
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
    //    fprintf(stderr,"was here\n");
                fprintf(fh,"%u ", poly_idx);
            }
        }
@@ -698,6 +717,14 @@
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
                fprintf(fh,"%e %e %e\n", vesicle->poly_list->poly[i]->vlist->vtx[j]->x,vesicle->poly_list->poly[i]->vlist->vtx[j]->y, vesicle->poly_list->poly[i]->vlist->vtx[j]->z );
            }
        }
    }
    //filaments
    if(fil){
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
                fprintf(fh,"%e %e %e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
            }
        }
    }
@@ -720,14 +747,26 @@

    }
    
    //filaments
    if(fil){
        poly_idx=vlist->n+monono*polyno;
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
                fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono);
//        fprintf(stderr,"was here\n");
			
            }
        }

    }

    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
    for (i=2;i<(blist->n+monono*polyno)*2+1;i+=2){
    for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){
    fprintf(fh,"%u ",i);
    }
    fprintf(fh,"\n");
    fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
     for (i=0;i<blist->n+monono*polyno;i++){
     for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){
        fprintf(fh,"3 ");
    }

@@ -797,7 +836,10 @@
        CFG_SIMPLE_INT("nshell", &tape->nshell),
        CFG_SIMPLE_INT("npoly", &tape->npoly),
        CFG_SIMPLE_INT("nmono", &tape->nmono),
        CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
    CFG_SIMPLE_INT("nfil",&tape->nfil),
    CFG_SIMPLE_INT("nfono",&tape->nfono),
    CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
    CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
        CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
    CFG_SIMPLE_INT("pswitch",&tape->pswitch),
    CFG_SIMPLE_FLOAT("pressure",&tape->pressure),

 src/io.h

@@ -16,6 +16,9 @@
    long int nczmax;
    long int npoly;
    long int nmono;
    long int nfil;
    long int nfono;
    long int R_nucleus;
    long int pswitch;
        char *multiprocessing;
       long int brezveze0;

 src/poly.c

@@ -97,14 +97,14 @@
    else
    {
    /* Make filaments inside the vesicle. Helix with radius... Dist. between poly vertices put to 1*/
        dphi = 2*asin(1/2/vesicle->R_nucleus)*1.001;
        dh = dphi/2/M_PI*1.001;
        dphi = 2.0*asin(1.0/2.0/vesicle->R_nucleus)*1.001;
        dh = dphi/2.0/M_PI*1.001;
        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;
                poly_list->poly[i]->vlist->vtx[j]->z = ji*dh - (dh*poly_list->n*poly_list->poly[i]->vlist->n/2.0);
            }
        }
    }

 src/tape

@@ -23,12 +23,15 @@
k_spring=800

####### Filament (inside the vesicle) definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
nfil=0
# nmono is a number of monomers in each filament
nfono=10
# nfil is a number of filaments inside the vesicle
nfil=1
# nfono is a number of monomers in each filament
nfono=300
# Spring constant between monomers of the filament
k_sfring=800

####### Nucleus (inside the vesicle) ###########
# Radius of an impenetrable hard sphere inside the vesicle
R_nucleus=5

#######  Cell definitions ############
nxmax=60
@@ -38,11 +41,11 @@

####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
mcsweeps=5000
mcsweeps=300
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=1
inititer=0
#how many records do you want on the disk iteration are there in a run?
iterations=1000
iterations=300


#shut up if we are using cluster!!!

 src/timestep.c

@@ -23,7 +23,7 @@
        cell_occupation(vesicle);
        ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps);
            dump_state(vesicle,i);
        if(i>inititer){
        if(i>=inititer){
            write_vertex_xml_file(vesicle,i-inititer);
        }
    }
@@ -56,15 +56,25 @@
    }

    for(i=0;i<vesicle->poly_list->n;i++){
    for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
        rnvec[0]=drand48();
        rnvec[1]=drand48();
        rnvec[2]=drand48();
        retval=single_poly_vertex_move(vesicle,vesicle->poly_list->poly[i],vesicle->poly_list->poly[i]->vlist->vtx[j],rnvec);	
        for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            rnvec[0]=drand48();
            rnvec[1]=drand48();
            rnvec[2]=drand48();
            retval=single_poly_vertex_move(vesicle,vesicle->poly_list->poly[i],vesicle->poly_list->poly[i]->vlist->vtx[j],rnvec);	
        }
    }


    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            rnvec[0]=drand48();
            rnvec[1]=drand48();
            rnvec[2]=drand48();
            retval=single_filament_vertex_move(vesicle,vesicle->filament_list->poly[i],vesicle->filament_list->poly[i]->vlist->vtx[j],rnvec);	
        }
    }
 

//    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,3*vesicle->blist->n,(double)cnt/(double)vesicle->blist->n/3.0);
    if(retval);
    return TS_SUCCESS;

 src/vertexmove.c

@@ -254,3 +254,107 @@
    //END MONTE CARLOOOOOOO
    return TS_SUCCESS;
}




ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
    ts_uint i;
    ts_bool retval; 
    ts_uint cellidx; 
//    ts_double delta_energy;
    ts_double costheta,sintheta,phi,r;
    ts_double dist[2];
    //This will hold all the information of vtx and its neighbours
    ts_vertex backupvtx;
    ts_bond backupbond[2];
    memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));

    //random move in a sphere with radius stepsize:
    r=vesicle->stepsize*rn[0];
    phi=rn[1]*2*M_PI;
    costheta=2*rn[2]-1;
    sintheta=sqrt(1-pow(costheta,2));
    vtx->x=vtx->x+r*sintheta*cos(phi);
    vtx->y=vtx->y+r*sintheta*sin(phi);
    vtx->z=vtx->z+r*costheta;


    //distance with neighbours check
    for(i=0;i<vtx->bond_no;i++){
        dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
        if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
            vtx=memcpy((void *)vtx,(void *)&backupvtx,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,sizeof(ts_vertex));
        return TS_FAIL;
    } 

    //backup bonds
    for(i=0;i<vtx->bond_no;i++){
        memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
        vtx->bond[i]->bond_length=sqrt(dist[i]);
        bond_vector(vtx->bond[i]);
		
    }
	
    //if all the tests are successful, then energy for vtx and neighbours is calculated
//    delta_energy=0;
/*    for(i=0;i<vtx->neigh_no;i++){
//        memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
        xp = vtx->neigh[i]
        vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
        bond_energy(vtx->bond[i],poly);
        delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
    }

    if(vtx==poly->vlist->vtx[0]){
        delta_energy+=
            (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
            pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
		
    }


    if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
        if(exp(-delta_energy)< drand48() )
#endif
#ifdef TS_DOUBLE_FLOAT
        if(expf(-delta_energy)< (ts_float)drand48())
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
        if(expl(-delta_energy)< (ts_ldouble)drand48())
#endif
        {
    //not accepted, reverting changes
    vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
    for(i=0;i<vtx->bond_no;i++){
    vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
    }

    return TS_FAIL; 
    }
    }
	
*/	
//    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
    if(vtx->cell!=vesicle->clist->cell[cellidx]){
        retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
//        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
        if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);	
    }
//    if(oldcellidx);
    //END MONTE CARLOOOOOOO
    return TS_SUCCESS;
}

 src/vertexmove.h

@@ -2,3 +2,5 @@
ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn);

ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn);

ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn);