Isak's proposal for stretching energy calculation, take 1

Samo Penic

2018-06-21 620ba72765eda244ca3f3806b8ad709929e74b3c

Isak's proposal for stretching energy calculation,  take 1

 src/bondflip.c

@@ -40,6 +40,7 @@
    ts_double delta_energy_cv;
    ts_vertex *constvol_vtx_moved, *constvol_vtx_backup;
    ts_bool retval;
    ts_double temp_area=0.0;

    if(it->neigh_no< 3) return TS_FAIL;
    if(k->neigh_no< 3) return TS_FAIL;
@@ -176,6 +177,9 @@
/*    vesicle_volume(vesicle);
    fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
*/
    if(vesicle->tape->stretchswitch==1){
        temp_area=vesicle->area-lm->area-lp->area;
    }
/* fix data structure for flipped bond */
    ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1, vesicle->tape->w);

@@ -189,10 +193,13 @@
  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;
/*        oldenergy+=lm->energy+lp->energy;
        stretchenergy(vesicle,lm);
        stretchenergy(vesicle,lp);
        delta_energy+=lm->energy+lp->energy;
    */
        temp_area=temp_area+lm->area+lp->area;
        delta_energy+=stretchenergy2(temp_area,vesicle->tlist->a0*vesicle->tlist->n)-stretchenergy2(vesicle->area,vesicle->tlist->a0*vesicle->tlist->n);
    }

    delta_energy-=oldenergy;
@@ -339,11 +346,11 @@
//        fprintf(stderr,"Restoration complete!!!\n");
//    vesicle_volume(vesicle);
//    fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
    if(vesicle->tape->stretchswitch==1){
/*    if(vesicle->tape->stretchswitch==1){
        stretchenergy(vesicle,lm);
        stretchenergy(vesicle,lp);
    }

*/
        return TS_FAIL;
        }
    }
@@ -358,6 +365,10 @@
    if(vesicle->tape->constareaswitch==2){
        vesicle->area+=darea;
    }
    if(vesicle->tape->stretchswitch==1){
        vesicle->area=temp_area;
    }

    // delete all backups
    for(i=0;i<4;i++){
    free(bck_vtx[i]->neigh);

 src/energy.c

@@ -242,3 +242,8 @@
void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){
    triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2);
}

ts_double stretchenergy2(ts_double current_area,ts_double tensionless_area){
    return pow((current_area*tensionless_area),2)/tensionless_area;

}

 src/energy.h

@@ -10,4 +10,5 @@
ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old);

void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle);
ts_double stretchenergy2(ts_double current_area,ts_double tensionless_area);
#endif

 src/timestep.c

@@ -62,7 +62,7 @@
    centermass(vesicle);
    cell_occupation(vesicle);
    vesicle_volume(vesicle); //needed for constant volume at this moment
    vesicle_area(vesicle); //needed for constant area at this moment
    vesicle_area(vesicle); //needed for constant area and stretching energy at this moment
    if(V0<0.000001) 
        V0=vesicle->volume; 
    ts_fprintf(stdout,"Setting volume V0=%.17f\n",V0);

 src/vertexmove.c

@@ -22,6 +22,7 @@
    ts_bool retval; 
    ts_uint cellidx; 
    ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0, dstretchenergy=0.0;
    ts_double temp_area=0.0; //for stretching
    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;
@@ -115,7 +116,13 @@
    }
    //stretching energy 1 of 3
    if(vesicle->tape->stretchswitch==1){
        for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
        temp_area=vesicle->area;
        dstretchenergy-=stretchenergy2(vesicle->area,vesicle->tlist->a0*vesicle->tlist->n);
        for(i=0;i<vtx->tristar_no;i++){
            temp_area-=vtx->tristar[i]->area;
        }
        //0.5*vesicle->tape->xkA0*pow((vesicle->area-vesicle->tlist->a0*vesicle->tlist->n),2)/(vesicle->tlist->a0*vesicle->tlist->n);
        //for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
    }
    delta_energy=0;
    
@@ -197,10 +204,12 @@

    //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;
        for(i=0;i<vtx->tristar_no;i++){
            temp_area+=vtx->tristar[i]->area; 
            }
        dstretchenergy+=stretchenergy2(temp_area,vesicle->tlist->a0*vesicle->tlist->n);
        dstretchenergy=0.5*vesicle->tape->xkA0*dstretchenergy;
		
    }

    delta_energy+=dstretchenergy;    
@@ -235,13 +244,13 @@
    
    //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
/*    //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){
        constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
@@ -272,6 +281,10 @@
    if(vesicle->tape->constareaswitch==2){
        vesicle->area+=darea;
    }
    //stretching energy 3 of 3
    if(vesicle->tape->stretchswitch==1){
        vesicle->area=temp_area;
    }
//    if(oldcellidx);
    //END MONTE CARLOOOOOOO
//    vesicle_volume(vesicle);