| | |
| | | #include "vertexmove.h" |
| | | #include <string.h> |
| | | #include "constvol.h" |
| | | #include "globals.h" |
| | | |
| | | ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){ |
| | | ts_uint i; |
| | | 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; |
| | |
| | | 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); |
| | | |
| | |
| | | } |
| | | /* 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+= |
| | |
| | | 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); |
| | |
| | | |
| | | //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){ |