| src/energy.c | ●●●●● patch | view | raw | blame | history | |
| src/energy.h | ●●●●● patch | view | raw | blame | history | |
| src/general.h | ●●●●● patch | view | raw | blame | history | |
| src/poly.c | ●●●●● patch | view | raw | blame | history | |
| src/tape | ●●●●● patch | view | raw | blame | history | |
| src/timestep.c | ●●●●● patch | view | raw | blame | history | |
| src/vertexmove.c | ●●●●● patch | view | raw | blame | history | |
| src/vertexmove.h | ●●●●● patch | view | raw | blame | history |
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! 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 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;
@@ -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; 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 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; } 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; } 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);
