From 4f5ffcfd5ba38b6b7dd6d3b15de8bb8677537b9f Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Thu, 09 Apr 2020 16:15:56 +0000 Subject: [PATCH] Commiting old changes --- src/vertexmove.c | 51 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 files changed, 48 insertions(+), 3 deletions(-) diff --git a/src/vertexmove.c b/src/vertexmove.c index 2517f69..e74fa74 100644 --- a/src/vertexmove.c +++ b/src/vertexmove.c @@ -21,7 +21,7 @@ ts_double dist; ts_bool retval; ts_uint cellidx; - ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0; + ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0, dstretchenergy=0.0; ts_double costheta,sintheta,phi,r; //This will hold all the information of vtx and its neighbours ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL; @@ -88,6 +88,14 @@ } } + + // plane confinement check whether the new position of vertex will be out of bounds + if(vesicle->tape->plane_confinement_switch){ + if(vtx->z>vesicle->confinement_plane.z_max || vtx->z<vesicle->confinement_plane.z_min){ + vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); + return TS_FAIL; + } + } //#undef SQ //self avoidance check with distant vertices cellidx=vertex_self_avoidance(vesicle, vtx); @@ -113,9 +121,13 @@ for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area; } - + //stretching energy 1 of 3 + if(vesicle->tape->stretchswitch==1){ + for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy; + } delta_energy=0; - + + // vesicle_volume(vesicle); // fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume); @@ -191,6 +203,17 @@ } /* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */ delta_energy+=direct_force_energy(vesicle,vtx,backupvtx); + + //stretching energy 2 of 3 + if(vesicle->tape->stretchswitch==1){ + for(i=0;i<vtx->tristar_no;i++){ + stretchenergy(vesicle, vtx->tristar[i]); + dstretchenergy+=vtx->tristar[i]->energy; + } + } + + delta_energy+=dstretchenergy; + /* No poly-bond energy for now! if(vtx->grafted_poly!=NULL){ delta_energy+= @@ -198,6 +221,22 @@ pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k; } */ + +// plane confinement energy due to compressing force + if(vesicle->tape->plane_confinement_switch){ + if(vesicle->confinement_plane.force_switch){ + //substract old energy + if(abs(vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_max)>1e-10) { + delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-backupvtx[0].z,2); + delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-vtx->z,2); + } + if(abs(-vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_min)>1e-10) { + delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-backupvtx[0].z,2); + delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-vtx->z,2); + } + } + } + // fprintf(stderr, "DE=%f\n",delta_energy); //MONTE CARLOOOOOOOO // if(vtx->c!=0.0) printf("DE=%f\n",delta_energy); @@ -221,6 +260,12 @@ //update the normals of triangles that share bead i. for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); + //stretching energy 3 of 3 + if(vesicle->tape->stretchswitch==1){ + for(i=0;i<vtx->tristar_no;i++){ + stretchenergy(vesicle,vtx->tristar[i]); + } + } // fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z); if(vesicle->tape->constvolswitch == 1){ -- Gitblit v1.9.3