From bd826de2f539f2e48c8c01d2d7f9f34c7e97104a Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Fri, 13 May 2016 07:43:27 +0000
Subject: [PATCH] Fix in trisurf output, inhibiting print of successful reconstruction. Multiple fixes and improvements in python module. Added symlinking of tapes into the running directories and dumping tapes from snapshots into tape files.
---
src/bondflip.c | 231 ++++++++++++++++++++++++++++++++++++++++++++-------------
1 files changed, 177 insertions(+), 54 deletions(-)
diff --git a/src/bondflip.c b/src/bondflip.c
index fae59d9..858dc74 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -1,3 +1,4 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdlib.h>
#include<math.h>
#include "general.h"
@@ -12,6 +13,7 @@
//#include "io.h"
#include<stdio.h>
#include<string.h>
+#include "constvol.h"
ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn){
/*c Vertex and triangle (lm and lp) indexing for bond flip:
@@ -30,10 +32,14 @@
ts_vertex *k=bond->vtx2;
ts_uint nei,neip,neim;
ts_uint i,j;
- ts_double oldenergy, delta_energy, dvol=0.0;
+ ts_double oldenergy, delta_energy, dvol=0.0, darea=0.0;
ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
ts_vertex *kp,*km;
+
+ ts_double delta_energy_cv;
+ ts_vertex *constvol_vtx_moved, *constvol_vtx_backup;
+ ts_bool retval;
if(it->neigh_no< 3) return TS_FAIL;
if(k->neigh_no< 3) return TS_FAIL;
@@ -126,32 +132,36 @@
ts_triangle *bck_tria[4];
ts_bond *bck_bond;
ts_vertex *orig_vtx[]={k,it,kp,km};
-ts_triangle *orig_tria[]={lm,lp, lm2,lp1};
+ts_triangle *orig_tria[]={lm,lp,lm2,lp1};
-fprintf(stderr,"Backuping!!!\n");
+//fprintf(stderr,"Backuping!!!\n");
+ bck_bond=(ts_bond *)malloc(sizeof(ts_bond));
for(i=0;i<4;i++){
+/* fprintf(stderr,"vtx neigh[%d]=",i);
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n");
+*/
bck_vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
bck_tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
- bck_bond=(ts_bond *)malloc(sizeof(ts_bond));
+ memcpy((void *)bck_vtx[i],(void *)orig_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)bck_tria[i],(void *)orig_tria[i],sizeof(ts_triangle));
/* level 2 pointers */
+
bck_vtx[i]->neigh=(ts_vertex **)malloc(orig_vtx[i]->neigh_no*sizeof(ts_vertex *));
bck_vtx[i]->tristar=(ts_triangle **)malloc(orig_vtx[i]->tristar_no*sizeof(ts_triangle *));
bck_vtx[i]->bond=(ts_bond **)malloc(orig_vtx[i]->bond_no*sizeof(ts_bond *));
bck_tria[i]->neigh=(ts_triangle **)malloc(orig_tria[i]->neigh_no*sizeof(ts_triangle *));
- bck_vtx[i]=memcpy(bck_vtx[i],orig_vtx[i],sizeof(ts_vertex));
- bck_vtx[i]->neigh=memcpy(bck_vtx[i]->neigh,orig_vtx[i]->neigh,orig_vtx[i]->neigh_no*sizeof(ts_vertex *));
- bck_vtx[i]->tristar=memcpy(bck_vtx[i]->tristar,orig_vtx[i]->tristar,orig_vtx[i]->tristar_no*sizeof(ts_triangle *));
- bck_vtx[i]->bond=memcpy(bck_vtx[i]->bond,orig_vtx[i]->bond,orig_vtx[i]->bond_no*sizeof(ts_bond *));
-
-
- bck_tria[i]=memcpy(bck_tria[i],orig_tria[i],sizeof(ts_triangle));
- bck_tria[i]->neigh=memcpy(bck_tria[i]->neigh,orig_tria[i]->neigh,orig_tria[i]->neigh_no*sizeof(ts_triangle *));
-}
- bck_bond=memcpy(bck_bond,bond,sizeof(ts_bond));
-fprintf(stderr,"Backup complete!!!\n");
-/* end backup vertex */
+ memcpy((void *)bck_vtx[i]->neigh,(void *)orig_vtx[i]->neigh,orig_vtx[i]->neigh_no*sizeof(ts_vertex *));
+ memcpy((void *)bck_vtx[i]->tristar,(void *)orig_vtx[i]->tristar,orig_vtx[i]->tristar_no*sizeof(ts_triangle *));
+ memcpy((void *)bck_vtx[i]->bond,(void *)orig_vtx[i]->bond,orig_vtx[i]->bond_no*sizeof(ts_bond *));
+ memcpy((void *)bck_tria[i]->neigh,(void *)orig_tria[i]->neigh,orig_tria[i]->neigh_no*sizeof(ts_triangle *));
+}
+ memcpy(bck_bond,bond,sizeof(ts_bond));
+//fprintf(stderr,"Backup complete!!!\n");
+/* end backup vertex */
+
/* Save old energy */
oldenergy=0;
oldenergy+=k->xk* k->energy;
@@ -160,7 +170,11 @@
oldenergy+=it->xk* it->energy;
//Neigbours of k, it, km, kp don't change its energy.
- if(vesicle->pswitch == 1){dvol = -lm->volume - lp->volume;}
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){dvol = -lm->volume - lp->volume;}
+ if(vesicle->tape->constareaswitch==2){darea=-lm->area-lp->area;}
+/* vesicle_volume(vesicle);
+ fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
+*/
/* fix data structure for flipped bond */
ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1);
@@ -174,16 +188,109 @@
delta_energy+=it->xk* it->energy;
//Neigbours of k, it, km, kp don't change its energy.
- delta_energy-=oldenergy;
- if(vesicle->pswitch == 1){
+ delta_energy-=oldenergy;
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
dvol = dvol + lm->volume + lp->volume;
- delta_energy-= vesicle->pressure*dvol;
+ if(vesicle->pswitch==1) delta_energy-= vesicle->pressure*dvol;
}
+ if(vesicle->tape->constareaswitch==2){
+ darea=darea+lm->area+lp->area;
+/*check whether the dvol is gt than epsvol */
+ if(fabs(vesicle->area+darea-A0)>epsarea){
+ //restore old state.
+ /* restoration procedure copied from few lines below */
+ for(i=0;i<4;i++){
+ // fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ free(orig_vtx[i]->neigh);
+ free(orig_vtx[i]->tristar);
+ free(orig_vtx[i]->bond);
+ free(orig_tria[i]->neigh);
+ memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle));
+ // fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ /* level 2 pointers are redirected*/
+ }
+ memcpy(bond,bck_bond,sizeof(ts_bond));
+ for(i=0;i<4;i++){
+ free(bck_vtx[i]);
+ free(bck_tria[i]);
+ /* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no );
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n"); */
+ }
+ free(bck_bond);
+ return TS_FAIL;
+
+ }
+ }
+
+
+
+
+ if(vesicle->tape->constvolswitch == 2){
+ /*check whether the dvol is gt than epsvol */
+ if(fabs(vesicle->volume+dvol-V0)>epsvol){
+ //restore old state.
+ /* restoration procedure copied from few lines below */
+ for(i=0;i<4;i++){
+ // fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ free(orig_vtx[i]->neigh);
+ free(orig_vtx[i]->tristar);
+ free(orig_vtx[i]->bond);
+ free(orig_tria[i]->neigh);
+ memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle));
+ // fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ /* level 2 pointers are redirected*/
+ }
+ memcpy(bond,bck_bond,sizeof(ts_bond));
+ for(i=0;i<4;i++){
+ free(bck_vtx[i]);
+ free(bck_tria[i]);
+ /* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no );
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n"); */
+ }
+ free(bck_bond);
+ return TS_FAIL;
+
+ }
+
+ } else
+ if(vesicle->tape->constvolswitch == 1){
+ retval=constvolume(vesicle, it, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
+ if(retval==TS_FAIL){
+/* restoration procedure copied from few lines below */
+ for(i=0;i<4;i++){
+ // fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ free(orig_vtx[i]->neigh);
+ free(orig_vtx[i]->tristar);
+ free(orig_vtx[i]->bond);
+ free(orig_tria[i]->neigh);
+ memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle));
+ // fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ /* level 2 pointers are redirected*/
+ }
+ memcpy(bond,bck_bond,sizeof(ts_bond));
+ for(i=0;i<4;i++){
+ free(bck_vtx[i]);
+ free(bck_tria[i]);
+ /* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no );
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n"); */
+ }
+ free(bck_bond);
+ return TS_FAIL;
+ }
+ delta_energy+=delta_energy_cv;
+ }
+
/* MONTE CARLO */
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())
@@ -193,54 +300,70 @@
#endif
{
//not accepted, reverting changes
-// fprintf(stderr,"Failed to move, due to MC\n");
+ //restore all backups
+// fprintf(stderr,"Restoring!!!\n");
+ if(vesicle->tape->constvolswitch == 1){
+ constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
+ }
- // ts_flip_bond(kp,km,k,it, bond, lm,lp,lm2,lp1);
-// fprintf(stderr,"%e, %e, %e\n", lp->xnorm, lp->ynorm, lp->znorm);
+ for(i=0;i<4;i++){
+// fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ free(orig_vtx[i]->neigh);
+ free(orig_vtx[i]->tristar);
+ free(orig_vtx[i]->bond);
+ free(orig_tria[i]->neigh);
+ memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex));
+ memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle));
+// fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no );
+ /* level 2 pointers are redirected*/
+ }
+ memcpy(bond,bck_bond,sizeof(ts_bond));
- //restore all backups
- fprintf(stderr,"Restoring!!!\n");
+ for(i=0;i<4;i++){
+ free(bck_vtx[i]);
+ free(bck_tria[i]);
+/* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no );
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n"); */
+ }
-for(i=0;i<4;i++){
- /* level 2 pointers */
- orig_vtx[i]=memcpy(orig_vtx[i],bck_vtx[i],sizeof(ts_vertex));
- orig_vtx[i]->neigh=memcpy(orig_vtx[i]->neigh,bck_vtx[i]->neigh,bck_vtx[i]->neigh_no*sizeof(ts_vertex *));
- orig_vtx[i]->tristar=memcpy(orig_vtx[i]->tristar,bck_vtx[i]->tristar,bck_vtx[i]->tristar_no*sizeof(ts_triangle *));
- orig_vtx[i]->bond=memcpy(orig_vtx[i]->bond,bck_vtx[i]->bond,bck_vtx[i]->bond_no*sizeof(ts_bond *));
-
-
- orig_tria[i]=memcpy(orig_tria[i],bck_tria[i],sizeof(ts_triangle));
- orig_tria[i]->neigh=memcpy(orig_tria[i]->neigh,bck_tria[i]->neigh,bck_tria[i]->neigh_no*sizeof(ts_triangle *));
-}
- bond=memcpy(bond,bck_bond,sizeof(ts_bond));
-
-
- for(i=0;i<4;i++){
- vtx_free(bck_vtx[i]);
-
- free(bck_tria[i]->neigh);
- free(bck_tria[i]);
-
- }
free(bck_bond);
- fprintf(stderr,"Restoration complete!!!\n");
+// fprintf(stderr,"Restoration complete!!!\n");
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
- return TS_FAIL;
+ return TS_FAIL;
}
}
/* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
- fprintf(stderr,"SUCCESS!!!\n");
+// fprintf(stderr,"SUCCESS!!!\n");
+ if(vesicle->tape->constvolswitch == 2){
+ vesicle->volume+=dvol;
+ } else if(vesicle->tape->constvolswitch == 1){
+ constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
+ }
+ if(vesicle->tape->constareaswitch==2){
+ vesicle->area+=darea;
+ }
// delete all backups
for(i=0;i<4;i++){
- vtx_free(bck_vtx[i]);
-
+ free(bck_vtx[i]->neigh);
+ free(bck_vtx[i]->bond);
+ free(bck_vtx[i]->tristar);
+ free(bck_vtx[i]);
free(bck_tria[i]->neigh);
free(bck_tria[i]);
-
+/* fprintf(stderr,"Afret backup deletion vtx neigh[%d]=",i);
+ for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx);
+ fprintf(stderr,"\n");
+*/
}
- free(bck_bond);
+ free(bck_bond);
+
+// vesicle_volume(vesicle);
+// fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
return TS_SUCCESS;
}
--
Gitblit v1.9.3