First working version, not debugged and energy assigned to vesicles on timesteps
| | |
| | | /* vesicle_volume(vesicle); |
| | | fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume); |
| | | */ |
| | | |
| | | /* fix data structure for flipped bond */ |
| | | ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1, vesicle->tape->w); |
| | | |
| | |
| | | delta_energy+=it->xk* it->energy; |
| | | delta_energy+=bond->energy; /* attraction with neighboring vertices, that have spontaneous curvature */ |
| | | //Neigbours of k, it, km, kp don't change its energy. |
| | | if(vesicle->tape->stretchswitch==1){ |
| | | oldenergy+=lm->energy+lp->energy; |
| | | stretchenergy(vesicle,lm); |
| | | stretchenergy(vesicle,lp); |
| | | delta_energy+=lm->energy+lp->energy; |
| | | } |
| | | |
| | | delta_energy-=oldenergy; |
| | | if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){ |
| | |
| | | // fprintf(stderr,"Restoration complete!!!\n"); |
| | | // vesicle_volume(vesicle); |
| | | // fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume); |
| | | if(vesicle->tape->stretchswitch==1){ |
| | | stretchenergy(vesicle,lm); |
| | | stretchenergy(vesicle,lp); |
| | | } |
| | | |
| | | return TS_FAIL; |
| | | } |
| | |
| | | return vesicle->tape->F*ddp; |
| | | |
| | | } |
| | | |
| | | void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){ |
| | | triangle->energy=vesicle->tape->xkA0/vesicle->tlist->n*pow((triangle->area/vesicle->tlist->a0-1),2); |
| | | } |
| | |
| | | ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle); |
| | | inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w); |
| | | ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old); |
| | | |
| | | void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle); |
| | | #endif |
| | |
| | | |
| | | struct ts_triangle_list{ |
| | | ts_uint n; |
| | | ts_double a0; |
| | | ts_triangle **tria; |
| | | }; |
| | | typedef struct ts_triangle_list ts_triangle_list; |
| | |
| | | vesicle->sphHarmonics=NULL; |
| | | } |
| | | |
| | | |
| | | vesicle->tlist->a0=sqrt(3)/2*pow((vesicle->tape->dmax+1.0),2); |
| | | return TS_SUCCESS; |
| | | |
| | | } |
| | |
| | | ts_bond_list *blist=vesicle->blist; |
| | | ts_vertex **vtx=vlist->vtx; |
| | | ts_uint i,j; |
| | | ts_double senergy=0.0; |
| | | char filename[10000]; |
| | | char just_name[255]; |
| | | FILE *fh; |
| | |
| | | } |
| | | } |
| | | fprintf(fh,"</DataArray>\n"); |
| | | |
| | | if(vesicle->tape->stretchswitch==1){ |
| | | fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"ascii\">"); |
| | | for(i=0;i<vlist->n;i++){ |
| | | senergy=0.0; |
| | | for(j=0;j<vtx[i]->tristar_no;j++){ |
| | | senergy+=vtx[i]->tristar[j]->energy; |
| | | } |
| | | fprintf(fh,"%.17e ",senergy); |
| | | } |
| | | //polymeres |
| | | if(poly){ |
| | | poly_idx=vlist->n; |
| | | for(i=0;i<vesicle->poly_list->n;i++){ |
| | | for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){ |
| | | fprintf(fh,"%.17e ", 0.0); |
| | | } |
| | | } |
| | | } |
| | | //filaments |
| | | if(fil){ |
| | | poly_idx=vlist->n+monono*polyno; |
| | | for(i=0;i<vesicle->filament_list->n;i++){ |
| | | for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){ |
| | | // fprintf(stderr,"was here\n"); |
| | | fprintf(fh,"%.17e ", 0.0); |
| | | } |
| | | } |
| | | } |
| | | fprintf(fh,"</DataArray>\n"); |
| | | } |
| | | |
| | | |
| | | fprintf(fh,"</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n"); |
| | |
| | | |
| | | #Stretching |
| | | stretchswitch=1 |
| | | xkA0=10.0 |
| | | xkA0=1.0 |
| | | |
| | | ####### Polymer (brush) definitions ########### |
| | | # npoly is a number of polymers attached to npoly distinct vertices on vesicle |
| | |
| | | ####### Program Control ############ |
| | | #how many MC sweeps between subsequent records of states to disk |
| | | #200000 |
| | | mcsweeps=2000 |
| | | mcsweeps=200 |
| | | #how many initial mcsweeps*inititer MC sweeps before recording to disk? |
| | | #2 |
| | | inititer=1 |
| | |
| | | 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/3.0; |
| | | } |
| | | delta_energy=0; |
| | | |
| | | // vesicle_volume(vesicle); |
| | |
| | | } |
| | | /* 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/3.0; |
| | | } |
| | | } |
| | | |
| | | delta_energy+=dstretchenergy; |
| | | |
| | | /* No poly-bond energy for now! |
| | | if(vtx->grafted_poly!=NULL){ |
| | | 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){ |