From 414b8ab65b13c32594178ba3d09b16d34f14bb19 Mon Sep 17 00:00:00 2001
From: mihaf <miha.fosnaric@gmail.com>
Date: Fri, 07 Mar 2014 15:15:08 +0000
Subject: [PATCH] finished constant volume. Problems still exists.
---
src/tape | 10 +++++-----
src/bondflip.c | 23 ++++++++++++++---------
src/vertexmove.c | 12 ++++++++++--
3 files changed, 29 insertions(+), 16 deletions(-)
diff --git a/src/bondflip.c b/src/bondflip.c
index da2c664..4935901 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -29,7 +29,7 @@
ts_vertex *k=bond->vtx2;
ts_uint nei,neip,neim;
ts_uint i,j;
- ts_double oldenergy, delta_energy;
+ ts_double oldenergy, delta_energy, dvol=0.0;
ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
ts_vertex *kp,*km;
@@ -113,19 +113,19 @@
}
if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999);
- // fprintf(stderr,"I WAS HERE! Before bondflip!\n");
- // fprintf(stderr,"%e, %e, %e\n", lm->xnorm, lm->ynorm, lm->znorm);
- ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1);
-// fprintf(stderr,"I WAS HERE! Bondflip successful!\n");
-
+/* Save old energy */
oldenergy=0;
oldenergy+=k->xk* k->energy;
oldenergy+=kp->xk* kp->energy;
oldenergy+=km->xk* km->energy;
oldenergy+=it->xk* it->energy;
- //Neigbours don't change its energy.
+ //Neigbours of k, it, km, kp don't change its energy.
+ if(vesicle->pswitch == 1){dvol = -lm->volume - lp->volume;}
+
+/* fix data structure for flipped bond */
+ ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1);
/* Calculating the new energy */
@@ -134,9 +134,14 @@
delta_energy+=kp->xk* kp->energy;
delta_energy+=km->xk* km->energy;
delta_energy+=it->xk* it->energy;
- //Neigbours don't change its energy.
+ //Neigbours of k, it, km, kp don't change its energy.
+
delta_energy-=oldenergy;
- // fprintf(stderr,"I WAS HERE! Got energy!\n");
+ if(vesicle->pswitch == 1){
+ dvol = dvol + lm->volume + lp->volume;
+ delta_energy-= vesicle->pressure*dvol;
+ }
+
/* MONTE CARLO */
if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
diff --git a/src/tape b/src/tape
index 1c0d279..2d9f978 100644
--- a/src/tape
+++ b/src/tape
@@ -4,19 +4,19 @@
# dmax is the max. bond length (in units l_min)
dmax=1.7
# bending rigidity of the membrane (in units kT)
-xk0=25.0
+xk0=1.0
# max step size (in units l_min)
stepsize=0.15
# Pressure calculations
# (pswitch=1: calc. p*dV energy contribution)
-pswitch = 1
+pswitch = 0
# pressure difference: p_inside - p_outside (in units l_min^3/kT):
pressure=0.0
####### Polymer definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
-npoly=30
+npoly=0
# nmono is a number of monomers in each polymer
nmono=10
# Spring constant between monomers of the polymer
@@ -30,11 +30,11 @@
####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
-mcsweeps=5000
+mcsweeps=50
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=1
#how many records do you want on the disk iteration are there in a run?
-iterations=2000
+iterations=10
#shut up if we are using cluster!!!
diff --git a/src/vertexmove.c b/src/vertexmove.c
index 86b8533..7f6391f 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -18,7 +18,7 @@
ts_double dist;
ts_bool retval;
ts_uint cellidx;
- ts_double delta_energy,oenergy;
+ ts_double delta_energy,oenergy,dvol=0.0;
ts_double costheta,sintheta,phi,r;
//This will hold all the information of vtx and its neighbours
ts_vertex backupvtx[20];
@@ -83,7 +83,9 @@
memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
}
-
+ if(vesicle->pswitch == 1){
+ for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
+ };
delta_energy=0;
//update the normals of triangles that share bead i.
@@ -97,6 +99,12 @@
energy_vertex(vtx->neigh[i]);
delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
}
+
+ if(vesicle->pswitch == 1){
+ for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
+ delta_energy-=vesicle->pressure*dvol;
+ };
+
/* No poly-bond energy for now!
if(vtx->grafted_poly!=NULL){
delta_energy+=
--
Gitblit v1.9.3