finished constant volume. Problems still exists.

mihaf

2014-03-07 414b8ab65b13c32594178ba3d09b16d34f14bb19

finished constant volume. Problems still exists.

 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

 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!!!

 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+=