From 58230a2591414fb38b9ec8d3a76439b290cb0a6f Mon Sep 17 00:00:00 2001
From: mihaf <miha.fosnaric@gmail.com>
Date: Tue, 18 Mar 2014 13:09:23 +0000
Subject: [PATCH] Adding ending energy to filaments. In progress. Continue in single_filament_vertex_move
---
src/timestep.c | 22 +++-
src/poly.c | 6
src/vertexmove.h | 2
src/bond.h | 1
src/io.c | 54 +++++++++-
src/tape | 19 ++-
src/bond.c | 11 ++
src/initial_distribution.c | 15 +++
src/general.h | 7
src/frame.c | 36 ++++---
src/vertexmove.c | 104 ++++++++++++++++++++
src/io.h | 3
12 files changed, 239 insertions(+), 41 deletions(-)
diff --git a/src/bond.c b/src/bond.c
index a348d61..4d11c97 100644
--- a/src/bond.c
+++ b/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++){
diff --git a/src/bond.h b/src/bond.h
index e2be586..34cdba5 100644
--- a/src/bond.h
+++ b/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);
diff --git a/src/frame.c b/src/frame.c
index d2d7e10..aeaef21 100644
--- a/src/frame.c
+++ b/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;
}
diff --git a/src/general.h b/src/general.h
index 26e8b55..79da2ba 100644
--- a/src/general.h
+++ b/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;
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index f450dbf..c7eed90 100644
--- a/src/initial_distribution.c
+++ b/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 */
diff --git a/src/io.c b/src/io.c
index db4e64f..9e23560 100644
--- a/src/io.c
+++ b/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),
diff --git a/src/io.h b/src/io.h
index eb0b737..71390bc 100644
--- a/src/io.h
+++ b/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;
diff --git a/src/poly.c b/src/poly.c
index aa41184..d99d2ad 100644
--- a/src/poly.c
+++ b/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);
}
}
}
diff --git a/src/tape b/src/tape
index 38871bd..ac316ef 100644
--- a/src/tape
+++ b/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!!!
diff --git a/src/timestep.c b/src/timestep.c
index 92a4a25..d631d42 100644
--- a/src/timestep.c
+++ b/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;
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 7f6391f..04325ef 100644
--- a/src/vertexmove.c
+++ b/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;
+}
diff --git a/src/vertexmove.h b/src/vertexmove.h
index 6387ab8..c58ee64 100644
--- a/src/vertexmove.h
+++ b/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);
--
Gitblit v1.9.3