From 1121faa13a2038facad22073f0fc610903d98691 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo@CAE-linux.(none)>
Date: Fri, 05 Sep 2014 20:18:05 +0000
Subject: [PATCH] First variant of constant volume constrant (a new one proposed by Miha after reading his article). It seems to work, however there are still some things to be done, such as Miha's derivation of the epsvol (0.003% is used at the moment) and solving the problem of additional global variables.
---
src/vertexmove.c | 48 +++++++++++++++++++++++++++++++++++++++++-------
1 files changed, 41 insertions(+), 7 deletions(-)
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 9f263c8..28a8880 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -92,14 +92,14 @@
memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
}
- if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch==1){
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
};
delta_energy=0;
-
-
-// fprintf(stderr,"Success for now.\n");
+
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
//update the normals of triangles that share bead i.
for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
@@ -113,22 +113,47 @@
delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
}
- if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch == 1){
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch >0){
for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
};
+
+ if(vesicle->tape->constvolswitch==2){
+ /*check whether the dvol is gt than epsvol */
+ //fprintf(stderr,"DVOL=%1.16e\n",dvol);
+ if(fabs(vesicle->volume+dvol-V0)>epsvol){
+ //restore old state.
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ for(i=0;i<vtx->neigh_no;i++){
+ vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
+ }
+ for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
+ //fprintf(stderr,"fajlam!\n");
+ return TS_FAIL;
+ }
+
+ } else
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume before=%1.16e\n", vesicle->volume);
if(vesicle->tape->constvolswitch == 1){
- retval=constvolume(vesicle, vtx, dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
+ retval=constvolume(vesicle, vtx, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
if(retval==TS_FAIL){ // if we couldn't move the vertex to assure constant volume
vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
for(i=0;i<vtx->neigh_no;i++){
vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
}
for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
+ // fprintf(stderr,"fajlam!\n");
return TS_FAIL;
}
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume after=%1.16e\n", vesicle->volume);
+// fprintf(stderr,"Volume after-dvol=%1.16e\n", vesicle->volume-dvol);
+// fprintf(stderr,"Denergy before=%e\n",delta_energy);
+
delta_energy+=delta_energy_cv;
+// fprintf(stderr,"Denergy after=%e\n",delta_energy);
}
/* No poly-bond energy for now!
if(vtx->grafted_poly!=NULL){
@@ -141,7 +166,7 @@
//MONTE CARLOOOOOOOO
if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
- if(exp(-delta_energy)< drand48() )
+ if(exp(-delta_energy)< drand48())
#endif
#ifdef TS_DOUBLE_FLOAT
if(expf(-delta_energy)< (ts_float)drand48())
@@ -160,9 +185,13 @@
//update the normals of triangles that share bead i.
for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(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){
constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
}
+// fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
return TS_FAIL;
}
}
@@ -176,11 +205,16 @@
}
+ if(vesicle->tape->constvolswitch == 2){
+ vesicle->volume+=dvol;
+ } else
if(vesicle->tape->constvolswitch == 1){
constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
}
// if(oldcellidx);
//END MONTE CARLOOOOOOO
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
return TS_SUCCESS;
}
--
Gitblit v1.9.3