From 267cf14b1ced6fce480e4292131a987a1802e28a Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Sun, 23 Sep 2018 08:32:24 +0000
Subject: [PATCH] Confinement planes, not debugged
---
src/timestep.c | 23 +++++++++++
src/io.c | 3 +
src/tape | 8 ++++
src/general.h | 12 +++++
src/vertexmove.c | 18 ++++++++
5 files changed, 62 insertions(+), 2 deletions(-)
diff --git a/src/general.h b/src/general.h
index a0e31a5..503486f 100644
--- a/src/general.h
+++ b/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;
@@ -320,7 +327,8 @@
ts_double R_nucleusY;
ts_double R_nucleusZ;
ts_double nucleus_center[3];
- ts_double area;
+ ts_double area;
+ ts_confinement_plane confinement_plane;
} ts_vesicle;
@@ -339,6 +347,8 @@
} ts_cluster_list;
+
+
/* GLOBAL VARIABLES */
int quiet;
diff --git a/src/io.c b/src/io.c
index 5a5e9a2..1710d17 100644
--- a/src/io.c
+++ b/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;
diff --git a/src/tape b/src/tape
index 8f38645..8e4cc9a 100644
--- a/src/tape
+++ b/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
diff --git a/src/timestep.c b/src/timestep.c
index 27400e7..204e9ba 100644
--- a/src/timestep.c
+++ b/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++){
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 2517f69..a6fb0e2 100644
--- a/src/vertexmove.c
+++ b/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())
--
Gitblit v1.9.3