Confinement planes, not debugged

Samo Penic

2018-09-23 267cf14b1ced6fce480e4292131a987a1802e28a

Confinement planes, not debugged

 src/general.h

@@ -289,8 +289,15 @@
    ts_double c0;
    ts_double w;
    ts_double F;
    long int plane_confinement_switch;
    ts_double plane_d;
    ts_double plane_F;
} ts_tape;

typedef struct{
    ts_float z_max;
    ts_float z_min;
} ts_confinement_plane;



@@ -321,6 +328,7 @@
    ts_double R_nucleusZ;
    ts_double nucleus_center[3];
    ts_double area;
    ts_confinement_plane confinement_plane;
} ts_vesicle;


@@ -339,6 +347,8 @@
} ts_cluster_list;




/* GLOBAL VARIABLES */

int quiet;

 src/io.c

@@ -1185,6 +1185,9 @@
    CFG_SIMPLE_FLOAT("c0",&tape->c0),
    CFG_SIMPLE_FLOAT("w",&tape->w),
    CFG_SIMPLE_FLOAT("F",&tape->F),
    CFG_SIMPLE_INT("plane_confinement_switch", &tape->plane_confinement_switch),
    CFG_SIMPLE_FLOAT("plane_d", &tape->plane_d),
    CFG_SIMPLE_FLOAT("plane_F", &tape->plane_F),
        CFG_END()
    };
    cfg_t *cfg;    

 src/tape

@@ -94,3 +94,11 @@
w=10.0
#direct force on vesicles with spontaneous curvature (float)
F=2.0


#plane confinement
plane_confinement_switch=1
#final plane distance (float in lmin)
plane_d=10.0
#plane to vesicle repulsion force while closing
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;
@@ -71,11 +72,33 @@
    ts_fprintf(stdout,"Setting area A0=%.17f\n",A0);
    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;


  //  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){
            if(vesicle->confinement_plane.z_max-vesicle->confinement_plane.z_min<=vesicle->tape->plane_d){
                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;
                }
                if(max_z-min_z<=vesicle->tape->plane_d) {
                    vesicle->confinement_plane.z_max=max_z;
                    vesicle->confinement_plane.z_min=max_z-vesicle->tape->plane_d;
                }  else {
                    vesicle->confinement_plane.z_min=min_z-1e-5;
                    vesicle->confinement_plane.z_max=max_z+1e-5;
                }
            }	
        }

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

 src/vertexmove.c

@@ -99,7 +99,13 @@
        return TS_FAIL;
    } 
   
 
// plane confinement test
    if(vesicle->tape->plane_confinement_switch){
        if(vtx->z >= vesicle->confinement_plane.z_max || vtx->z <= vesicle->confinement_plane.z_min){
                vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
                return TS_FAIL;
        }
    }
//if all the tests are successful, then energy for vtx and neighbours is calculated
    for(i=0;i<vtx->neigh_no;i++){
    memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
@@ -201,6 +207,16 @@
//   fprintf(stderr, "DE=%f\n",delta_energy);
    //MONTE CARLOOOOOOOO
//    if(vtx->c!=0.0) printf("DE=%f\n",delta_energy);
// plane confinement
    if(vesicle->tape->plane_confinement_switch){
        //if planes are not close enough, then repusion force is on
        if(vesicle->confinement_plane.z_max-vesicle->confinement_plane.z_min > vesicle->tape->plane_d){	
            delta_energy+=vesicle->tape->plane_F * 1.0/( (backupvtx->z-vesicle->confinement_plane.z_min) + (backupvtx->z-vesicle->confinement_plane.z_max) );
            delta_energy+=-(vesicle->tape->plane_F * 1.0/( (vtx->z-vesicle->confinement_plane.z_min) + (vtx->z-vesicle->confinement_plane.z_max) ) );

        }
    }
// end plane confinement
    if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
        if(exp(-delta_energy)< drand48())