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/vertexmove.c | 104 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 files changed, 104 insertions(+), 0 deletions(-) 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; +} -- Gitblit v1.9.3