Debugged

Samo Penic

2018-09-25 5eefcc281fda031d322599fd7094da4fc8450948

Debugged

 src/tape

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

 src/timestep.c

@@ -24,7 +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;
    ts_double max_z=0, min_z=0;
    FILE *fd1,*fd2=NULL,*fd3=NULL;
     char filename[10000];
    //struct stat st;
@@ -72,7 +72,21 @@
    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;

    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;
                }
                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;
                }
    }

  //  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);
@@ -83,6 +97,10 @@
        // plane confinement
        if(vesicle->tape->plane_confinement_switch){
            if(vesicle->confinement_plane.z_max-vesicle->confinement_plane.z_min<=vesicle->tape->plane_d){
                ts_fprintf(stdout,"Force is off\n");
            } else {
                ts_fprintf(stdout,"Force is on\n");
			
                min_z=1e10;
                max_z=-1e10;
                for(k=0;k<vesicle->vlist->n;k++){
@@ -93,10 +111,11 @@
                    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->confinement_plane.z_min=min_z;//-1e-5;
                    vesicle->confinement_plane.z_max=max_z;//+1e-5;
                }
            }	
            }
            ts_fprintf(stdout,"Planes (zmin, zmax)=(%e,%e)\n",vesicle->confinement_plane.z_min,vesicle->confinement_plane.z_max);	
        }

/*    vesicle_volume(vesicle);

 src/vertexmove.c

@@ -101,7 +101,7 @@
   
// plane confinement test
    if(vesicle->tape->plane_confinement_switch){
        if(vtx->z >= vesicle->confinement_plane.z_max || vtx->z <= vesicle->confinement_plane.z_min){
        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;
        }
@@ -211,8 +211,9 @@
    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) ) );
            //ts_fprintf(stderr,"force!!!\n");
            delta_energy-=vesicle->tape->plane_F * ( 1.0/pow(backupvtx->z-vesicle->confinement_plane.z_min,2) + 1.0/pow(-backupvtx->z+vesicle->confinement_plane.z_max,2) );
            delta_energy+=vesicle->tape->plane_F * ( 1.0/pow(vtx->z-vesicle->confinement_plane.z_min,2) + +1.0/pow(-vtx->z+vesicle->confinement_plane.z_max,2) );

        }
    }