From fedf2b217113800a15f367e111e2dbcad4a3c92c Mon Sep 17 00:00:00 2001 From: Miha <mihaf@hapi.(none)> Date: Thu, 06 Feb 2014 17:39:37 +0000 Subject: [PATCH] Polymers thermicaly move now! HArmonic potential added. --- src/timestep.c | 20 +++++- src/poly.c | 2 src/vertexmove.h | 1 src/tape | 6 +- src/energy.c | 9 ++- src/general.h | 3 src/energy.h | 1 src/vertexmove.c | 103 +++++++++++++++++++++++++++++++++ 8 files changed, 132 insertions(+), 13 deletions(-) diff --git a/src/energy.c b/src/energy.c index 3b1779c..1e37e08 100644 --- a/src/energy.c +++ b/src/energy.c @@ -19,9 +19,12 @@ return TS_SUCCESS; } -inline ts_bool energy_poly_vertex(ts_vertex *vtx,ts_poly *poly){ - -} + +inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){ + bond->energy=poly->k*pow(bond->bond_length-1,2); + return TS_SUCCESS; +}; + inline ts_bool energy_vertex(ts_vertex *vtx){ // ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value! diff --git a/src/energy.h b/src/energy.h index 90463d8..8a37b85 100644 --- a/src/energy.h +++ b/src/energy.h @@ -2,4 +2,5 @@ #define _ENERGY_H ts_bool mean_curvature_and_energy(ts_vesicle *vesicle); inline ts_bool energy_vertex(ts_vertex *vtx); +inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly); #endif diff --git a/src/general.h b/src/general.h index 0372e7a..e7e9545 100644 --- a/src/general.h +++ b/src/general.h @@ -164,12 +164,13 @@ } ts_vertex_list; struct ts_bond { - ts_uint idx; + ts_uint idx; ts_vertex *vtx1; ts_vertex *vtx2; ts_double bond_length; ts_double bond_length_dual; ts_bool tainted; + ts_double energy; }; typedef struct ts_bond ts_bond; diff --git a/src/poly.c b/src/poly.c index edab18f..1ddcec7 100644 --- a/src/poly.c +++ b/src/poly.c @@ -4,6 +4,7 @@ #include"vertex.h" #include"bond.h" #include<math.h> +#include"energy.h" ts_bool poly_assign_spring_const(ts_vesicle *vesicle){ ts_uint i; @@ -30,6 +31,7 @@ for(i=0;i<poly->blist->n;i++){ poly->blist->bond[i]->bond_length=sqrt(vtx_distance_sq(poly->blist->bond[i]->vtx1,poly->blist->bond[i]->vtx2)); + bond_energy(poly->blist->bond[i],poly); } return poly; diff --git a/src/tape b/src/tape index c8c8cdb..02c85b8 100644 --- a/src/tape +++ b/src/tape @@ -1,8 +1,8 @@ ####### Vesicle definitions ########### # nshell is a number of divisions of dipyramid nshell=17 -# dmax is the square of the bond length -dmax=1.67 +# dmax is the max. bond length +dmax=1.7 # bending rigidity of the membrane xk0=25.0 # max step size @@ -10,7 +10,7 @@ ####### Polymer definitions ########### # npoly is a number of polymers attached to npoly distinct vertices on vesicle -npoly=10 +npoly=20 # nmono is a number of monomers in each polymer nmono=15 # Spring constant between monomers of the polymer diff --git a/src/timestep.c b/src/timestep.c index 9e5e61c..0ce9372 100644 --- a/src/timestep.c +++ b/src/timestep.c @@ -8,6 +8,7 @@ #include "bondflip.h" #include "frame.h" #include "io.h" + ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){ ts_uint i, j; @@ -30,7 +31,7 @@ ts_bool single_timestep(ts_vesicle *vesicle){ ts_bool retval; ts_double rnvec[3]; - ts_uint i, b; + ts_uint i,j,b; for(i=0;i<vesicle->vlist->n;i++){ rnvec[0]=drand48(); rnvec[1]=drand48(); @@ -39,7 +40,7 @@ } // ts_int cnt=0; - for(i=0;i<vesicle->vlist->n;i++){ + for(i=0;i<3*vesicle->vlist->n;i++){ //why is rnvec needed in bondflip? /* rnvec[0]=drand48(); rnvec[1]=drand48(); @@ -50,8 +51,19 @@ //call single_bondflip_timestep... retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec); // if(retval==TS_SUCCESS) cnt++; - } -// printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n); + } + + 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); + } + + } + +// 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 a79b2f0..f78fd61 100644 --- a/src/vertexmove.c +++ b/src/vertexmove.c @@ -13,8 +13,7 @@ #include "vertexmove.h" #include <string.h> -ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double -*rn){ +ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){ ts_uint i; ts_double dist; ts_bool retval; @@ -88,6 +87,15 @@ energy_vertex(vtx->neigh[i]); delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy); } + + if(vtx->grafted_poly!=NULL){ + delta_energy+= + (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)- + pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k; + + + } + // fprintf(stderr, "DE=%f\n",delta_energy); //MONTE CARLOOOOOOOO if(delta_energy>=0){ @@ -126,3 +134,94 @@ return TS_SUCCESS; } + +ts_bool single_poly_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; + //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->neigh_no;i++){ +// dist=vtx_distance_sq(vtx,vtx->neigh[i]); +// if(dist<1.0 || dist>vesicle->dmax) { +// 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,sizeof(ts_vertex)); + return TS_FAIL; + } + + + //if all the tests are successful, then energy for vtx and neighbours is calculated + delta_energy=0; + for(i=0;i<vtx->bond_no;i++){ + memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond)); + + 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 0c24d54..6387ab8 100644 --- a/src/vertexmove.h +++ b/src/vertexmove.h @@ -1,3 +1,4 @@ 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); -- Gitblit v1.9.3