From e86357faa6caf28f5032c624fa0cba24e5dc4686 Mon Sep 17 00:00:00 2001
From: mihaf <miha.fosnaric@gmail.com>
Date: Fri, 07 Mar 2014 12:13:21 +0000
Subject: [PATCH] Many fixes, before changing bondflip. triangles are still calculated in ts_bond_flip
---
src/io.c | 9 +++-
src/tape | 14 +++++--
src/bondflip.c | 37 ++++++++++++------
src/initial_distribution.c | 3 -
src/general.h | 1
5 files changed, 43 insertions(+), 21 deletions(-)
diff --git a/src/bondflip.c b/src/bondflip.c
index 478027d..2d41849 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -79,10 +79,7 @@
oldenergy+=kp->xk* kp->energy;
oldenergy+=km->xk* km->energy;
oldenergy+=it->xk* it->energy;
-// for(i=0;i<k->neigh_no;i++) oldenergy+=k->neigh[i]->xk*k->neigh[i]->energy;
-// for(i=0;i<kp->neigh_no;i++) oldenergy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
-// for(i=0;i<km->neigh_no;i++) oldenergy+=km->neigh[i]->xk*km->neigh[i]->energy;
-// for(i=0;i<it->neigh_no;i++) oldenergy+=it->neigh[i]->xk*it->neigh[i]->energy;
+ //Neigbours don't change its energy.
/*
fprintf(stderr,"*** Naslov k=%ld\n",(long)k);
fprintf(stderr,"*** Naslov it=%ld\n",(long)it);
@@ -107,18 +104,34 @@
/* Calculating the new energy */
delta_energy=0;
- for(i=0;i<k->neigh_no;i++) energy_vertex(k->neigh[i]);
- for(i=0;i<kp->neigh_no;i++) energy_vertex(kp->neigh[i]);
- for(i=0;i<km->neigh_no;i++) energy_vertex(km->neigh[i]);
- for(i=0;i<it->neigh_no;i++) energy_vertex(it->neigh[i]);
+/* ts_double dbg_energy=0.0;
+ for(i=0;i<k->neigh_no;i++) {
+ dbg_energy=k->neigh[i]->energy;
+ energy_vertex(k->neigh[i]);
+ if(abs(dbg_energy-k->neigh[i]->energy)>1e-100) fatal("Energy changes!",1);
+ }
+ for(i=0;i<kp->neigh_no;i++) {
+ dbg_energy=kp->neigh[i]->energy;
+ energy_vertex(kp->neigh[i]);
+ if(abs(dbg_energy-kp->neigh[i]->energy)>1e-100) fatal("Energy changes!",1);
+ }
+ for(i=0;i<km->neigh_no;i++){
+ dbg_energy=km->neigh[i]->energy;
+ energy_vertex(km->neigh[i]);
+ if(abs(dbg_energy-km->neigh[i]->energy)>1e-100) fatal("Energy changes!",1);
+ }
+
+ for(i=0;i<it->neigh_no;i++){
+ dbg_energy=it->neigh[i]->energy;
+ energy_vertex(it->neigh[i]);
+ if(abs(dbg_energy-it->neigh[i]->energy)>1e-100) fatal("Energy changes!",1);
+ }
+*/
delta_energy+=k->xk* k->energy;
delta_energy+=kp->xk* kp->energy;
delta_energy+=km->xk* km->energy;
delta_energy+=it->xk* it->energy;
-// for(i=0;i<k->neigh_no;i++) delta_energy+=k->neigh[i]->xk*k->neigh[i]->energy;
-// for(i=0;i<kp->neigh_no;i++) delta_energy+=kp->neigh[i]->xk*kp->neigh[i]->energy;
-// for(i=0;i<km->neigh_no;i++) delta_energy+=km->neigh[i]->xk*km->neigh[i]->energy;
-// for(i=0;i<it->neigh_no;i++) delta_energy+=it->neigh[i]->xk*it->neigh[i]->energy;
+ //Neigbours don't change its energy.
delta_energy-=oldenergy;
// fprintf(stderr,"I WAS HERE! Got energy!\n");
/* MONTE CARLO */
diff --git a/src/general.h b/src/general.h
index e7e9545..f4e5833 100644
--- a/src/general.h
+++ b/src/general.h
@@ -260,6 +260,7 @@
ts_poly_list *poly_list;
ts_double spring_constant;
+ ts_double pressure;
} ts_vesicle;
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 382750b..b8f5965 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -19,9 +19,6 @@
ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize);
-//TODO: debugging only. Please remove ASAP!
- vesicle->bending_rigidity=25.0;
-
vesicle->nshell=nshell;
retval = vtx_set_global_values(vesicle);
retval = pentagonal_dipyramid_vertex_distribution(vesicle->vlist);
diff --git a/src/io.c b/src/io.c
index 8b51793..4db7f56 100644
--- a/src/io.c
+++ b/src/io.c
@@ -341,7 +341,7 @@
long int brezveze0=1;
long int brezveze1=1;
long int brezveze2=1;
- ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0;
+ ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0,pressure=0.0;
long int iter=1000, init=1000, mcsw=1000;
@@ -351,7 +351,8 @@
CFG_SIMPLE_INT("nmono", &nmono),
CFG_SIMPLE_FLOAT("dmax", &dmax),
CFG_SIMPLE_FLOAT("xk0",&xk0),
- CFG_SIMPLE_FLOAT("k_spring",&kspring),
+ CFG_SIMPLE_FLOAT("pressure",&pressure),
+ CFG_SIMPLE_FLOAT("k_spring",&kspring),
CFG_SIMPLE_FLOAT("stepsize",&stepsize),
CFG_SIMPLE_INT("nxmax", &ncxmax),
CFG_SIMPLE_INT("nymax", &ncymax),
@@ -389,11 +390,15 @@
vesicle->nshell=nshell;
vesicle->dmax=dmax*dmax;
vesicle->bending_rigidity=xk0;
+ vtx_set_global_values(vesicle); //copies xk0 to every vertex
+
+
vesicle->stepsize=stepsize;
vesicle->clist->ncmax[0]=ncxmax;
vesicle->clist->ncmax[1]=ncymax;
vesicle->clist->ncmax[2]=nczmax;
vesicle->clist->max_occupancy=8;
+ vesicle->pressure=pressure/vesicle->bending_rigidity; //all energy contributions need to be divided by bending_rigidity!
cfg_free(cfg);
free(buf);
diff --git a/src/tape b/src/tape
index 9a44ef3..1c0d279 100644
--- a/src/tape
+++ b/src/tape
@@ -1,12 +1,18 @@
####### Vesicle definitions ###########
# nshell is a number of divisions of dipyramid
nshell=17
-# dmax is the max. bond length
+# dmax is the max. bond length (in units l_min)
dmax=1.7
-# bending rigidity of the membrane
+# bending rigidity of the membrane (in units kT)
xk0=25.0
-# max step size
+# max step size (in units l_min)
stepsize=0.15
+
+# Pressure calculations
+# (pswitch=1: calc. p*dV energy contribution)
+pswitch = 1
+# 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
@@ -28,7 +34,7 @@
#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=10000
+iterations=2000
#shut up if we are using cluster!!!
--
Gitblit v1.9.3