First test successful

Samo Penic

2018-12-11 0dd5baa7166ab9abd7ef2d6b374e72beab03ef2a

First test successful

 src/frame.c

@@ -9,6 +9,7 @@
ts_bool centermass(ts_vesicle *vesicle){
    ts_uint i,j, n=vesicle->vlist->n;
    ts_vertex **vtx=vesicle->vlist->vtx;
    ts_double temp_z_cm=0;
    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
@@ -20,6 +21,12 @@
    vesicle->cm[0]/=(ts_float)n;
    vesicle->cm[1]/=(ts_float)n;
    vesicle->cm[2]/=(ts_float)n;

    //center mass for z component does not  change if we confine the vesicle with plates
    if(vesicle->tape->plane_confinement_switch){
        temp_z_cm=vesicle->cm[2];
        vesicle->cm[2]=0;
    }

    for(i=0;i<n;i++){
        vtx[i]->x-=vesicle->cm[0];
@@ -49,7 +56,9 @@

    vesicle->cm[0]=0.0;
    vesicle->cm[1]=0.0;
    vesicle->cm[2]=0.0;
    if(vesicle->tape->plane_confinement_switch){
        vesicle->cm[2]=temp_z_cm;
    }

    for(i=0;i<vesicle->tlist->n;i++){
        triangle_normal_vector(vesicle->tlist->tria[i]);

 src/tape

@@ -100,6 +100,6 @@
#plane confinement
plane_confinement_switch=1
#final plane distance (float in lmin)
plane_d=15
plane_d=10
#plane to vesicle repulsion force while closing
plane_F=1000
plane_F=10

 src/timestep.c

@@ -24,6 +24,7 @@
    ts_double r0,kc1=0,kc2=0,kc3=0,kc4=0;
    ts_double l1,l2,l3,vmsr,bfsr, vmsrt, bfsrt;
    ts_ulong epochtime;
    ts_double max_z,min_z;
    FILE *fd1,*fd2=NULL,*fd3=NULL;
     char filename[10000];
    //struct stat st;
@@ -72,19 +73,40 @@
    epsvol=4.0*sqrt(2.0*M_PI)/pow(3.0,3.0/4.0)*V0/pow(vesicle->tlist->n,3.0/2.0);
    epsarea=A0/(ts_double)vesicle->tlist->n;

    //plane confinement part 1

    if(vesicle->tape->plane_confinement_switch){
        vesicle->confinement_plane.z_min=-vesicle->tape->plane_d/2.0;
        vesicle->confinement_plane.z_max=vesicle->tape->plane_d/2.0;
        ts_fprintf(stderr,"Vesicle confinement by plane set to (zmin, zmax)=(%e,%e).\n",vesicle->confinement_plane.z_min,vesicle->confinement_plane.z_max);
    }
	
  //  fprintf(stderr, "DVol=%1.16f (%1.16f), V0=%1.16f\n", epsvol,0.003e-2*V0,V0);
    if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
    for(i=start_iteration;i<inititer+iterations;i++){
        vmsr=0.0;
        bfsr=0.0;

    //plane confinement
    if(vesicle->tape->plane_confinement_switch){
        min_z=1e10;
        max_z=-1e10;
        for(k=0;k<vesicle->vlist->n;k++){
            if(vesicle->vlist->vtx[k]->z > max_z) max_z=vesicle->vlist->vtx[k]->z;
            if(vesicle->vlist->vtx[k]->z < min_z) min_z=vesicle->vlist->vtx[k]->z;
        }
        vesicle->confinement_plane.force_switch=0;
        if(max_z>=vesicle->tape->plane_d/2.0){
            ts_fprintf(stdout, "Max vertex out of bounds (z>=%e). Plane set to max_z = %e.\n",vesicle->tape->plane_d/2.0,max_z);
            vesicle->confinement_plane.z_max = max_z;
            vesicle->confinement_plane.force_switch=1;
        } else {
            vesicle->confinement_plane.z_max=vesicle->tape->plane_d/2.0;
        }
        if(min_z<=-vesicle->tape->plane_d/2.0){
            ts_fprintf(stdout, "Min vertex out of bounds (z<=%e). Plane set to min_z = %e.\n",-vesicle->tape->plane_d/2.0,min_z);
            vesicle->confinement_plane.z_min = min_z;
            vesicle->confinement_plane.force_switch=1;
        } else {
            vesicle->confinement_plane.z_min=-vesicle->tape->plane_d/2.0;
        }
        ts_fprintf(stdout,"Vesicle confinement by plane set to (zmin, zmax)=(%e,%e).\n",vesicle->confinement_plane.z_min,vesicle->confinement_plane.z_max);
        if(vesicle->confinement_plane.force_switch) ts_fprintf(stdout,"Squeezing with force %e.\n",vesicle->tape->plane_F);
    }

    //end plane confinement

/*    vesicle_volume(vesicle);
    fprintf(stderr,"Volume before TS=%1.16e\n", vesicle->volume); */
        for(j=0;j<mcsweeps;j++){

 src/vertexmove.c

@@ -124,6 +124,7 @@

    delta_energy=0;
    

//    vesicle_volume(vesicle);
//    fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);

@@ -206,6 +207,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);