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/vertexmove.c | 112 ++++++++++++++++++++++++++++++++++++++++++++------------
1 files changed, 88 insertions(+), 24 deletions(-)
diff --git a/src/vertexmove.c b/src/vertexmove.c
index a27a776..b065cbb 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -1,3 +1,4 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<stdlib.h>
#include<math.h>
#include "general.h"
@@ -20,7 +21,7 @@
ts_double dist;
ts_bool retval;
ts_uint cellidx;
- ts_double delta_energy,oenergy,dvol=0.0;
+ ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0;
ts_double costheta,sintheta,phi,r;
//This will hold all the information of vtx and its neighbours
ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL;
@@ -41,7 +42,7 @@
// vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
// vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
- //random move in a sphere with radius stepsize:
+//random move in a sphere with radius stepsize:
r=vesicle->stepsize*rn[0];
phi=rn[1]*2*M_PI;
costheta=2*rn[2]-1;
@@ -51,12 +52,12 @@
vtx->z=vtx->z+r*costheta;
- //distance with neighbours check
+//distance with neighbours check
for(i=0;i<vtx->neigh_no;i++){
dist=vtx_distance_sq(vtx,vtx->neigh[i]);
if(dist<1.0 || dist>vesicle->dmax) {
- vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
- return TS_FAIL;
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ return TS_FAIL;
}
}
@@ -87,27 +88,24 @@
}
- //if all the tests are successful, then energy for vtx and neighbours is calculated
+//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));
}
- if(vesicle->pswitch == 1){
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
- };
+ }
+
+ if(vesicle->tape->constareaswitch==2){
+ for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
+
+ }
delta_energy=0;
-
- if(vesicle->tape->constvolswitch == 1){
- retval=constvolume(vesicle, vtx, dvol, &delta_energy, 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));
- }
- return TS_FAIL;
- }
- }
+
+// 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]);
@@ -121,12 +119,64 @@
delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
}
- if(vesicle->pswitch == 1){
+ if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch >0){
for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
- delta_energy-=vesicle->pressure*dvol;
+ if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
};
+ if(vesicle->tape->constareaswitch==2){
+ /* check whether the darea is gt epsarea */
+ for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
+ if(fabs(vesicle->area+darea-A0)>epsarea){
+ //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;
+ }
+
+ }
+
+ 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);
+ 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){
delta_energy+=
@@ -138,7 +188,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())
@@ -148,6 +198,7 @@
#endif
{
//not accepted, reverting changes
+ // fprintf(stderr,"MC failed\n");
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));
@@ -156,14 +207,18 @@
//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;
}
}
//accepted
+ // fprintf(stderr,"MC accepted\n");
// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
if(vtx->cell!=vesicle->clist->cell[cellidx]){
retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
@@ -172,11 +227,20 @@
}
+ if(vesicle->tape->constvolswitch == 2){
+ vesicle->volume+=dvol;
+ } else
if(vesicle->tape->constvolswitch == 1){
- constvolumeaccept(constvol_vtx_moved,constvol_vtx_backup);
+ constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
+ }
+
+ if(vesicle->tape->constareaswitch==2){
+ vesicle->area+=darea;
}
// 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