First working version, not debugged and energy assigned to vesicles on time...

Samo Penic

2018-05-15 7ec6fb63c690080d0624c7bf138d0fd84fe7664e

First working version, not debugged and energy assigned to vesicles on timesteps

 src/bondflip.c

@@ -176,7 +176,6 @@
/*    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);

@@ -189,6 +188,12 @@
  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){
@@ -334,6 +339,10 @@
//        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;
        }

 src/energy.c

@@ -238,3 +238,7 @@
    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);
}

 src/energy.h

@@ -8,4 +8,6 @@
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

 src/general.h

@@ -201,6 +201,7 @@

struct ts_triangle_list{
    ts_uint n;
    ts_double a0;
    ts_triangle **tria;
};
typedef struct ts_triangle_list ts_triangle_list;

 src/initial_distribution.c

@@ -112,7 +112,7 @@
        vesicle->sphHarmonics=NULL;
    }

    
    vesicle->tlist->a0=sqrt(3)/2*pow((vesicle->tape->dmax+1.0),2);  
    return TS_SUCCESS;

}

 src/io.c

@@ -826,6 +826,7 @@
    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;
@@ -973,7 +974,36 @@
            }
        }
    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");

 src/tape

@@ -26,7 +26,7 @@

#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
@@ -61,7 +61,7 @@
####### 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

 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;
@@ -113,7 +113,10 @@
        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);
@@ -191,6 +194,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/3.0;
            }
    }

    delta_energy+=dstretchenergy;	
		
/* No poly-bond energy for now!
    if(vtx->grafted_poly!=NULL){
        delta_energy+=
@@ -221,6 +235,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){