dont know

Samo Penic

2014-03-05 719c9febac2eaff9483fda487b57684afbb59bb2

dont know

 aclocal.m4

@@ -1,4 +1,4 @@
# generated automatically by aclocal 1.11.6 -*- Autoconf -*-
# generated automatically by aclocal 1.11.3 -*- Autoconf -*-

# Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
# 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation,
@@ -14,8 +14,8 @@

m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.69],,
[m4_warning([this file was generated for autoconf 2.69.
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],,
[m4_warning([this file was generated for autoconf 2.68.
You have another version of autoconf.  It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])
@@ -38,7 +38,7 @@
[am__api_version='1.11'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version.  Point them to the right macro.
m4_if([$1], [1.11.6], [],
m4_if([$1], [1.11.3], [],
      [AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])

@@ -54,7 +54,7 @@
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
[AM_AUTOMAKE_VERSION([1.11.6])dnl
[AM_AUTOMAKE_VERSION([1.11.3])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])

 src/Makefile.am

@@ -1,21 +1,25 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c
#trisurf_LDFLAGS = -lm -lconfuse
shdiscoverdir=../
shdiscover_PROGRAMS= shdiscover
shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c
shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c poly.c
AM_CFLAGS = -Wall -Werror
co_testdir=../
co_test_PROGRAMS=co_test
co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c
co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c poly.c
spherical_trisurfdir=../
spherical_trisurf_PROGRAMS = spherical_trisurf
spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c
spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c poly.c
spherical_trisurf_ffdir=../
spherical_trisurf_ff_PROGRAMS = spherical_trisurf_ff
<<<<<<< HEAD
spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c
neigh_testdir=../
neigh_test_PROGRAMS = neigh_test
neigh_test_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c neigh_test.c

=======
spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c poly.c
>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768

 src/bondflip.c

@@ -15,10 +15,10 @@
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:
c      +----- k-------+              +----- k ------+
c      |lm1 / | \ lp1 |              |lm1 /   \ lp1 |
c      |lm1 / ^ \ lp1 |              |lm1 /   \ lp1 |
c      |  /   |   \   |              |  /       \   |
c      |/     |     \ |     FLIP     |/    lm     \ |
c     km  lm  | lp   kp    --->      km ---------- kp  
c     km  lm  | lp   kp    --->      km ---------> kp  
c      |\     |     / |              |\    lp     / |  
c      |  \   |   /   |              |  \       /   |
c      |lm2 \ | / lp2 |              |lm2 \   / lp2 |
@@ -113,7 +113,7 @@
#endif
        {
            //not accepted, reverting changes
       //     fprintf(stderr,"Failed to move, due to MC\n");
            fprintf(stderr,"Failed to move, due to MC\n");

//            ts_flip_bond(km,kp,it,k, bond);
            ts_flip_bond(bond);
@@ -142,7 +142,7 @@
            return TS_FAIL;
        }
    }
        //    fprintf(stderr,"Success\n");
            fprintf(stderr,"Success\n");


/* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
@@ -153,49 +153,8 @@
ts_bool ts_flip_bond(ts_bond *bond){
fprintf(stderr,"Called!\n");
    ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
    ts_uint i; //lmidx, lpidx;
/*This is no longer necessary! if(k==NULL || it==NULL || km==NULL || kp==NULL){
    fatal("ts_flip_bond: You called me with invalid pointers to vertices",999);
}*/
// 1. step. We find lm and lp from k->tristar !
/*    for(i=0;i<it->tristar_no;i++){
        for(j=0;j<k->tristar_no;j++){
            if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
                if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
== km || it->tristar[i]->vertex[2]== km )){
                lm=it->tristar[i];
         //       lmidx=i;
                }
                else
                {
                lp=it->tristar[i];
         //       lpidx=i;
                }
    ts_uint i,j; //lmidx, lpidx;

            }
        }
    }
if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999);

//we look for important triangles lp1 and lm2.

 for(i=0;i<k->tristar_no;i++){
        for(j=0;j<kp->tristar_no;j++){
                if((k->tristar[i] == kp->tristar[j]) && k->tristar[i]!=lp){ //ce gre za skupen trikotnik
                    lp1=k->tristar[i];
            }
        }
}

 for(i=0;i<it->tristar_no;i++){
        for(j=0;j<km->tristar_no;j++){
            if((it->tristar[i] == km->tristar[j]) && it->tristar[i]!=lm){ //ce gre za skupen trikotnik
                    lm2=it->tristar[i];
            } 
        }
    }
*/
//Upper block no longer neccessary.
ts_vertex *it=bond->vtx1;
ts_vertex *k=bond->vtx2;
ts_vertex *km=bond->adjvtx[0];
@@ -219,6 +178,22 @@
*/
//find bonds between k and kp (adj[1] vtx)
//find bonds between it and km (adj[0] vtx)

for(i=0;i<3;i++){
if(lm->neigh[i]==lp) break;
    for(j=0;j<3;j++){
        if(lm->neigh[i]->vertex[j]==it) lm2=lm->neigh[i];
    }
}

for(i=0;i<3;i++){
if(lp->neigh[i]==lm) break;
    for(j=0;j<3;j++){
        if(lp->neigh[i]->vertex[j]==k) lp1=lp->neigh[i];
    }
}

/*
for(i=0;i<it->tristar_no;i++){
    if((it->tristar[i]->vertex[0]==km || it->tristar[i]->vertex[1]==km || it->tristar[i]->vertex[2]==km) && (it->tristar[i]->vertex[0]==it || it->tristar[i]->vertex[1]==it || it->tristar[i]->vertex[2]==it)){
        lm2=it->tristar[i];
@@ -227,7 +202,11 @@
        lp1=it->tristar[i];
    }
}
*/
if(lm2==NULL ) fatal("ts_flip_bond: Cannot find triangles lm2!",999);
if(lp1==NULL ) fatal("ts_flip_bond: Cannot find triangles lp1!",999);
if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999);

/*
//DEBUG TESTING
fprintf(stderr,"1. step: lm, lm2, lp1 and lp found!\n");
@@ -347,9 +326,9 @@
  energy_vertex(it);

//Extra steps for new bond data structure

bond->adjvtx[0]=k;
bond->adjvtx[1]=it;

// END modifications to data structure!



 src/cell.c

@@ -147,7 +147,7 @@
                }
// Now we check whether we didn't come close to some other vertices in the same
// cell!
                if(cell_occupation>1){
                if(cell_occupation>0){
                    for(l=0;l<cell_occupation;l++){

                //carefull with this checks!

 src/energy.c

@@ -20,6 +20,14 @@
}


inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){
//TODO: This value to be changed and implemented in data structure:
    ts_double d_relaxed=1.0;
    bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2);
    return TS_SUCCESS;
};


inline ts_bool energy_vertex(ts_vertex *vtx){
//    ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
//    ts_triangle *tristar=vtx->tristar-1;

 src/energy.h

@@ -2,4 +2,5 @@
#define _ENERGY_H
ts_bool mean_curvature_and_energy(ts_vesicle *vesicle);
inline ts_bool energy_vertex(ts_vertex *vtx);
inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly);
#endif

 src/frame.c

@@ -3,7 +3,7 @@
#include "cell.h"
#include "frame.h"
ts_bool centermass(ts_vesicle *vesicle){
    ts_uint i, n=vesicle->vlist->n;
    ts_uint i,j, n=vesicle->vlist->n;
    ts_vertex **vtx=vesicle->vlist->vtx;
    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
@@ -23,6 +23,15 @@
        vtx[i]->z-=vesicle->cm[2]; 
    } 

    for(i=0;i<vesicle->poly_list->n;i++){
        for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            vesicle->poly_list->poly[i]->vlist->vtx[j]->x-=vesicle->cm[0];
            vesicle->poly_list->poly[i]->vlist->vtx[j]->y-=vesicle->cm[1];
            vesicle->poly_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
        }
    }


    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
@@ -31,11 +40,11 @@
}

ts_bool cell_occupation(ts_vesicle *vesicle){
    ts_uint i,cellidx, n=vesicle->vlist->n;
    ts_double shift;
    ts_double dcell;
    shift=(ts_double) vesicle->clist->ncmax[0]/2;
    dcell=1.0/(1.0 + vesicle->stepsize);
    ts_uint i,j,cellidx, n=vesicle->vlist->n;
    //ts_double shift;
    //ts_double dcell;
    //shift=(ts_double) vesicle->clist->ncmax[0]/2;
    //dcell=1.0/(1.0 + vesicle->stepsize);
    //`fprintf(stderr, "Bil sem tu\n"); 

    cell_list_cell_occupation_clear(vesicle->clist);
@@ -45,10 +54,19 @@
//    vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];

    cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);

    }

    for(i=0;i<vesicle->poly_list->n;i++){
    for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
        cellidx=vertex_self_avoidance(vesicle, vesicle->poly_list->poly[i]->vlist->vtx[j]);
        cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->poly_list->poly[i]->vlist->vtx[j]);
    }
    }

    

    //fprintf(stderr, "Bil sem tu\n"); 
    if(dcell);
    if(shift);
    //if(dcell);
    //if(shift);
    return TS_SUCCESS;
}

 src/general.h

@@ -12,6 +12,7 @@
  * Header file for general inclusion in all the code, defining data structures
  * and general constans. All datatypes used in the code is also defined here.
  *
  * Miha: branch trisurf-polyel
  */

/* Defines */
@@ -130,9 +131,12 @@
    ts_uint coord_type;
} ts_coord;

<<<<<<< HEAD



=======
>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768
/** @brief Data structure of all data connected to a vertex
 *
 *  ts_vertex holds the data for one single point (bead, vertex). To understand how to use it
@@ -161,6 +165,7 @@
        ts_double projArea;
        ts_double relR;
        ts_double solAngle;
    struct ts_poly *grafted_poly;
};
typedef struct ts_vertex ts_vertex;

@@ -181,6 +186,7 @@
    ts_double bond_length;
    ts_double bond_length_dual;
    ts_bool tainted;
    ts_double energy;
};
typedef struct ts_bond ts_bond;

@@ -237,6 +243,24 @@



struct ts_poly {
    ts_vertex_list *vlist;
    ts_bond_list *blist;
    ts_vertex *grafted_vtx;
    ts_double k;
};
typedef struct ts_poly ts_poly;


struct ts_poly_list {
    ts_uint    n;
    ts_poly **poly;
};
typedef struct ts_poly_list ts_poly_list;




typedef struct {
    ts_vertex_list *vlist;
    ts_bond_list *blist;
@@ -249,6 +273,9 @@
    ts_double cm[3];
    ts_double volume;
    ts_spharm *sphHarmonics;

    ts_poly_list *poly_list;
    ts_double spring_constant;
} ts_vesicle;



 src/io.c

@@ -9,6 +9,329 @@
#include <sys/types.h>
#include <dirent.h>
#include "initial_distribution.h"
#include "poly.h"



/** DUMP STATE TO DISK DRIVE **/

ts_bool dump_state(ts_vesicle *vesicle){

    /* save current state with wrong pointers. Will fix that later */
    ts_uint i,j,k;
    ts_ulong off;
    FILE *fh=fopen("dump.bin","wb");

    /* dump vesicle */
    fwrite(vesicle, sizeof(vesicle),1,fh);
    /* dump vertex list */
    fwrite(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
    /* dump bond list */
    fwrite(vesicle->blist, sizeof(ts_bond_list),1,fh);
    /* dump triangle list */
    fwrite(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
    /* dump cell list */
    fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
    /* dump poly list */
    fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
    /* level 1 complete */

    /*dump vertices*/
    for(i=0;i<vesicle->vlist->n;i++){
        fwrite(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
        /* dump pointer offsets for:
                    neigh
                    bond
                    tria
                    cell is ignored
        */
        for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
            off=(ts_ulong)(vesicle->vlist->vtx[i]->neigh[j]-vesicle->vlist->vtx[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
        }
        for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
            off=(ts_ulong)(vesicle->vlist->vtx[i]->bond[j]-vesicle->blist->bond[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
        }
        for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
            off=(ts_ulong)(vesicle->vlist->vtx[i]->tristar[j]-vesicle->tlist->tria[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
        }
//        off=(ts_ulong)(vesicle->vlist->vtx[i]->cell-vesicle->clist->cell[0]);
//        fwrite(&off,sizeof(ts_ulong),1,fh); 

    }

    /*dump bonds*/
    for(i=0;i<vesicle->blist->n;i++){
        fwrite(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
        /* dump pointer offsets for vtx1 and vtx2 */
        off=(ts_ulong)(vesicle->blist->bond[i]->vtx1-vesicle->vlist->vtx[0]);
        fwrite(&off,sizeof(ts_ulong),1,fh); 
        off=(ts_ulong)(vesicle->blist->bond[i]->vtx2-vesicle->vlist->vtx[0]);
        fwrite(&off,sizeof(ts_ulong),1,fh); 
    }

    /*dump triangles*/
    for(i=0;i<vesicle->tlist->n;i++){
        fwrite(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
        /* dump pointer offsets for vertex */
        off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[0]-vesicle->vlist->vtx[0]);
        fwrite(&off,sizeof(ts_ulong),1,fh); 
        off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[1]-vesicle->vlist->vtx[0]);
        fwrite(&off,sizeof(ts_ulong),1,fh); 
        off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[2]-vesicle->vlist->vtx[0]);
        fwrite(&off,sizeof(ts_ulong),1,fh); 
        /* dump pointer offsets for neigh */
        for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
            off=(ts_ulong)(vesicle->tlist->tria[i]->neigh[j]-vesicle->tlist->tria[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
        }
    }
   
    /*dump cells and its contents is not dumped but recalculated on restore*/

/*    for(i=0;i<vesicle->clist->cellno;i++){
        fwrite(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
        for(j=0;j<vesicle->clist->cell[i]->nvertex;j++){
            off=(ts_ulong)(vesicle->clist->cell[i]->vertex[j]-vesicle->vlist->vtx[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh);
        }
    }
*/


    /*dump polymeres */
    for(i=0;i<vesicle->poly_list->n;i++){
        fwrite(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
    } 
     
    /* dump poly vertex(monomer) list*/
    for(i=0;i<vesicle->poly_list->n;i++){
        for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            fwrite(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
            /* dump offset for neigh and bond */
            for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
                off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]-vesicle->poly_list->poly[i]->vlist->vtx[0]);
                fwrite(&off,sizeof(ts_ulong),1,fh); 
            }
            for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
                off=(ts_ulong)(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]-vesicle->poly_list->poly[i]->blist->bond[0]);
                fwrite(&off,sizeof(ts_ulong),1,fh); 
            }
        }
    }
    /* dump poly bonds between monomers list*/
    for(i=0;i<vesicle->poly_list->n;i++){
        for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
            fwrite(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
            /* dump vtx1 and vtx2 offsets */
            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx1-vesicle->poly_list->poly[i]->vlist->vtx[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
            fwrite(&off,sizeof(ts_ulong),1,fh); 
        }
    }

/* pointer offsets for fixing the restored pointers */
/* need pointers for 
    vlist->vtx
    blist->bond
    tlist->tria
    clist->cell
    poly_list->poly
    and for each poly:
        poly_list->poly->vtx
        poly_list->poly->bond
*/

    fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
    fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
    fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
    fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
    fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
    for(i=0;i<vesicle->poly_list->n;i++){
        fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
        fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
    }
    fclose(fh);
    return TS_SUCCESS;
}


/** RESTORE DUMP FROM DISK **/
ts_vesicle *restore_state(){
    ts_uint i,j,k;
    FILE *fh=fopen("dump.bin","rb");
    ts_uint retval;
    ts_ulong off;

/* we restore all the data from the dump */
    /* restore vesicle */
    ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
    retval=fread(vesicle, sizeof(vesicle),1,fh);
    /* restore vertex list */
    vesicle->vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
    retval=fread(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
    /* restore bond list */
    vesicle->blist=(ts_bond_list *)malloc(sizeof(ts_bond_list));
    retval=fread(vesicle->blist, sizeof(ts_bond_list),1,fh);
    /* restore triangle list */
    vesicle->tlist=(ts_triangle_list *)malloc(sizeof(ts_triangle_list));
    retval=fread(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
    /* restore cell list */
    vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
    retval=fread(vesicle->clist, sizeof(ts_cell_list),1,fh);
    /* restore poly list */
    vesicle->poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
    retval=fread(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
    /* level 1 complete */

  /*restore vertices*/
    vesicle->vlist->vtx=(ts_vertex **)calloc(vesicle->vlist->n,sizeof(ts_vertex *));
    for(i=0;i<vesicle->vlist->n;i++){
        vesicle->vlist->vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
        retval=fread(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
        /*restore neigh, bond, tristar. Ignoring cell */
        vesicle->vlist->vtx[i]->neigh=(ts_vertex **)calloc(vesicle->vlist->vtx[i]->neigh_no, sizeof(ts_vertex *));
        for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
            retval=fread(&off,sizeof(ts_ulong),1,fh);
            vesicle->vlist->vtx[i]->neigh[j]=vesicle->vlist->vtx[0]+off;
        }
        vesicle->vlist->vtx[i]->bond=(ts_bond **)calloc(vesicle->vlist->vtx[i]->bond_no, sizeof(ts_bond *));
        for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
            retval=fread(&off,sizeof(ts_ulong),1,fh);
            vesicle->vlist->vtx[i]->bond[j]=vesicle->blist->bond[0]+off;    
        }

        vesicle->vlist->vtx[i]->tristar=(ts_triangle **)calloc(vesicle->vlist->vtx[i]->tristar_no, sizeof(ts_triangle *));
        for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
            retval=fread(&off,sizeof(ts_ulong),1,fh);
            vesicle->vlist->vtx[i]->tristar[j]=vesicle->tlist->tria[0]+off;
        }

    }

    /*restore bonds*/
    vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *));
    for(i=0;i<vesicle->blist->n;i++){
        vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond));
        retval=fread(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
        /* restore vtx1 and vtx2 */
        retval=fread(&off,sizeof(ts_ulong),1,fh);
        vesicle->blist->bond[i]->vtx1=vesicle->vlist->vtx[0]+off;
        retval=fread(&off,sizeof(ts_ulong),1,fh);
        vesicle->blist->bond[i]->vtx2=vesicle->vlist->vtx[0]+off;
    }

    /*restore triangles*/
    vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *));
    for(i=0;i<vesicle->tlist->n;i++){
        vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
        retval=fread(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
        /* restore pointers for vertices */
        retval=fread(&off,sizeof(ts_ulong),1,fh);
        vesicle->tlist->tria[i]->vertex[0]=vesicle->vlist->vtx[0]+off;
        retval=fread(&off,sizeof(ts_ulong),1,fh);
        vesicle->tlist->tria[i]->vertex[1]=vesicle->vlist->vtx[0]+off;
        retval=fread(&off,sizeof(ts_ulong),1,fh);
        vesicle->tlist->tria[i]->vertex[2]=vesicle->vlist->vtx[0]+off;
        /* restore pointers for neigh */
        for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
            retval=fread(&off,sizeof(ts_ulong),1,fh);
            vesicle->tlist->tria[i]->neigh[j]=vesicle->tlist->tria[0]+off;
        }

    }
   
    /*restore cells */
/*TODO: do we need to recalculate cells here? */
/*    vesicle->clist->cell=(ts_cell **)malloc(vesicle->clist->cellno*sizeof(ts_cell *));
    for(i=0;i<vesicle->clist->cellno;i++){
        vesicle->clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
        retval=fread(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
    }
*/
    /*restore polymeres */
    vesicle->poly_list->poly = (ts_poly **)calloc(vesicle->poly_list->n,sizeof(ts_poly *));
    for(i=0;i<vesicle->poly_list->n;i++){
        vesicle->poly_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
        retval=fread(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
    } 
     
    /* restore poly vertex(monomer) list*/
    for(i=0;i<vesicle->poly_list->n;i++){
        vesicle->poly_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->n,sizeof(ts_vertex *));
        for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
            vesicle->poly_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
            retval=fread(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
            
            /* restore neigh and bonds */
            vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *));
            for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
                retval=fread(&off,sizeof(ts_ulong),1,fh);
                vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
            }
            vesicle->poly_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *));
            for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
                retval=fread(&off,sizeof(ts_ulong),1,fh);
                vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->poly_list->poly[i]->blist->bond[0]+off;
            }

        }
    }
    /* restore poly bonds between monomers list*/
    for(i=0;i<vesicle->poly_list->n;i++){
        vesicle->poly_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->blist->n,sizeof(ts_bond *));
        for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
            vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
            retval=fread(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
            /* restore vtx1 and vtx2 */
                retval=fread(&off,sizeof(ts_ulong),1,fh);
                vesicle->poly_list->poly[i]->blist->bond[j]->vtx1=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
                retval=fread(&off,sizeof(ts_ulong),1,fh);
                vesicle->poly_list->poly[i]->blist->bond[j]->vtx2=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
        }
    }


    /* recalculate pointers using saved information on pointer addresses */
    /* following pointers need to be recalculated:
        all vtx->neigh
        all vtx->tristar
        all vtx->bond
        vtx->cell
        vtx->grafted_poly
        bond->vtx1
        bond->vtx2
        tria->vertex[*]
        tria->neigh
        cell->vertex
        poly->grafted_vtx
        poly->vlist->vtx->neigh
        poly->vlist->vtx->bond
        poly->vlist->vtx->cell
    */

    /* read pointers from disk */
/*    fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
    fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
    fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
    fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
    fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
    for(i=0;i<vesicle->poly_list->n;i++){
        fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
        fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
    }
*/


    if(retval); 
    fclose(fh);
    return vesicle;
}



ts_bool print_vertex_list(ts_vertex_list *vlist){
    ts_uint i;
@@ -206,7 +529,7 @@
    ts_vertex_list *vlist=vesicle->vlist;
    ts_bond_list *blist=vesicle->blist;
    ts_vertex **vtx=vlist->vtx;
    ts_uint i;
    ts_uint i,j;
        char filename[255];
    FILE *fh;

@@ -217,29 +540,69 @@
        return TS_FAIL;
    }
    /* Here comes header of the file */

    //find number of extra vtxs and bonds of polymeres
    ts_uint monono=0, polyno=0;
    ts_bool poly=0;
    if(vesicle->poly_list!=NULL){
        if(vesicle->poly_list->poly[0]!=NULL){
        polyno=vesicle->poly_list->n;
        monono=vesicle->poly_list->poly[0]->vlist->n;
        poly=1;
        }
    }
    fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n <UnstructuredGrid>\n");
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n, blist->n);
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno, blist->n+monono*polyno);
    fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
       for(i=0;i<vlist->n;i++){
        fprintf(fh,"%u ",vtx[i]->idx);
    }
    //polymeres
    if(poly){
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
                fprintf(fh,"%u ", vesicle->poly_list->poly[i]->vlist->vtx[j]->idx);
            }
        }
    }

    fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
    for(i=0;i<vlist->n;i++){
        fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
    }
    //polymeres
    if(poly){
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
                fprintf(fh,"%e %e %e\n", vesicle->poly_list->poly[i]->vlist->vtx[j]->x,vesicle->poly_list->poly[i]->vlist->vtx[j]->y, vesicle->poly_list->poly[i]->vlist->vtx[j]->z );
            }
        }
    }

    fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
    for(i=0;i<blist->n;i++){
            fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
    }
    //polymeres
    if(poly){
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
                fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx);
            }
    //grafted bonds
        fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx);
        }

    }
	

    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
    for (i=2;i<blist->n*2+1;i+=2){
    for (i=2;i<(blist->n+monono*polyno)*2+1;i+=2){
    fprintf(fh,"%u ",i);
    }
    fprintf(fh,"\n");
    fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
     for (i=0;i<blist->n;i++){
     for (i=0;i<blist->n+monono*polyno;i++){
        fprintf(fh,"3 ");
    }

@@ -295,17 +658,22 @@


ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){
    long int nshell=17,ncxmax=60, ncymax=60, nczmax=60;  // THIS IS DUE TO CONFUSE BUG!
    long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20;  // THIS IS DUE TO CONFUSE BUG!
    char *buf=malloc(255*sizeof(char));
    long int brezveze0=1;
    long int brezveze1=1;
    long int brezveze2=1;
    ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
    ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0;
    long int iter=1000, init=1000, mcsw=1000;


    cfg_opt_t opts[] = {
        CFG_SIMPLE_INT("nshell", &nshell),
        CFG_SIMPLE_INT("npoly", &npoly),
        CFG_SIMPLE_INT("nmono", &nmono),
        CFG_SIMPLE_FLOAT("dmax", &dmax),
        CFG_SIMPLE_FLOAT("xk0",&xk0),
        CFG_SIMPLE_FLOAT("k_spring",&kspring),
        CFG_SIMPLE_FLOAT("stepsize",&stepsize),
        CFG_SIMPLE_INT("nxmax", &ncxmax),
        CFG_SIMPLE_INT("nymax", &ncymax),
@@ -335,6 +703,11 @@
    *inititer=init;
    *mcsweeps=mcsw;
    vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
    vesicle->poly_list=init_poly_list(npoly,nmono, vesicle->vlist);
    vesicle->spring_constant=kspring;
    poly_assign_spring_const(vesicle);
	

    vesicle->nshell=nshell;
    vesicle->dmax=dmax*dmax;
    vesicle->bending_rigidity=xk0;
@@ -347,6 +720,9 @@
    cfg_free(cfg);
    free(buf);
  //  fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);

			

    return vesicle;

}

 src/io.h

@@ -50,4 +50,6 @@
ts_bool write_master_xml_file(ts_char *filename);
ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations);

ts_bool dump_state(ts_vesicle *vesicle);
ts_vesicle *restore_state();
#endif

 src/main.c

@@ -10,6 +10,7 @@
#include "initial_distribution.h"
#include "frame.h"
#include "timestep.h"
#include "poly.h"

/** Entrance function to the program
  * @param argv is a number of parameters used in program call (including the program name
@@ -64,6 +65,24 @@
printf("Tests complete.\n");
*/
vesicle=parsetape(&mcsweeps, &inititer, &iterations);

/*Testing */
//vesicle->poly_list=init_poly_list(1400,20,vesicle->vlist);

//poly_list_free(vesicle->poly_list);
/*End testing*/
printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
dump_state(vesicle);
vesicle->vlist->vtx[123]->x=123.0;
printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
write_vertex_xml_file(vesicle,0);
vesicle=restore_state();
printf("Stevilo vertexov je %d\n",vesicle->vlist->n);
printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
write_vertex_xml_file(vesicle,1);
write_master_xml_file("test.pvd");
return 0;

run_simulation(vesicle, mcsweeps, inititer, iterations);

write_master_xml_file("test.pvd");

 src/poly.c

New file
@@ -0,0 +1,120 @@
#include"general.h"
#include"poly.h"
#include<stdlib.h>
#include"vertex.h"
#include"bond.h"
#include<math.h>
#include"energy.h"

ts_bool poly_assign_spring_const(ts_vesicle *vesicle){
    ts_uint i;

    for(i=0;i<vesicle->poly_list->n;i++){
     vesicle->poly_list->poly[i]->k = vesicle->spring_constant;
        }
	
    return TS_SUCCESS;
}

ts_poly    *init_poly(ts_uint n, ts_vertex *grafted_vtx){
    ts_poly    *poly=(ts_poly *)calloc(1,sizeof(ts_poly));
    poly->vlist = init_vertex_list(n);
    poly->blist = init_bond_list();
    poly->grafted_vtx = grafted_vtx;
    grafted_vtx->grafted_poly = poly;

    ts_uint i;
    for(i=0;i<n-1;i++){
        vtx_add_cneighbour(poly->blist, poly->vlist->vtx[i], poly->vlist->vtx[i+1]);
        vtx_add_neighbour(poly->vlist->vtx[i+1], poly->vlist->vtx[i]);
    }

    for(i=0;i<poly->blist->n;i++){
    poly->blist->bond[i]->bond_length=sqrt(vtx_distance_sq(poly->blist->bond[i]->vtx1,poly->blist->bond[i]->vtx2));
    bond_energy(poly->blist->bond[i],poly);
    }

    return poly;
}


ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist){
    ts_poly_list *poly_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
    poly_list->poly    = (ts_poly **)calloc(n_poly,sizeof(ts_poly *));
    ts_uint i=0,j=0, idx;
    ts_uint gvtxi;
    ts_double xnorm,ynorm,znorm,normlength;

    if (n_poly > vlist->n){fatal("Number of polymers larger then numbero f vertices on a vesicle.",310);}
	
    while(i<n_poly){
        gvtxi = rand() % vlist->n;
        if (vlist->vtx[gvtxi]->grafted_poly == NULL){
        poly_list->poly[i] = init_poly(n_mono, vlist->vtx[gvtxi]);
        i++;
        }
    }
	
    poly_list->n = n_poly;

/* Make straight poylmers normal to membrane. Dist. between poly vertices put to 1*/
    for (i=0;i<poly_list->n;i++){

        xnorm=0.0;
        ynorm=0.0;
        znorm=0.0;
        for (j=0;j<poly_list->poly[i]->grafted_vtx->tristar_no;j++){
            xnorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->xnorm;
            ynorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->ynorm;
            znorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->znorm;	
        }
        normlength=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
        xnorm=xnorm/normlength;
        ynorm=ynorm/normlength;
        znorm=znorm/normlength;

        for (j=0;j<poly_list->poly[i]->vlist->n;j++){
            poly_list->poly[i]->vlist->vtx[j]->x = poly_list->poly[i]->grafted_vtx->x + xnorm*(ts_double)(j+1);
            poly_list->poly[i]->vlist->vtx[j]->y = poly_list->poly[i]->grafted_vtx->y + ynorm*(ts_double)(j+1);
            poly_list->poly[i]->vlist->vtx[j]->z = poly_list->poly[i]->grafted_vtx->z + znorm*(ts_double)(j+1);
        }
    }

        //index correction for polymeres. Important, since each vtx has to have unique id
    idx=vlist->n;
    for(i=0;i<n_poly;i++){
        for(j=0;j<n_mono;j++,idx++){

            poly_list->poly[i]->vlist->vtx[j]->idx=idx;

        }
    }


    return poly_list;
}


ts_bool poly_free(ts_poly *poly){

    if (poly->grafted_vtx!=NULL){
        poly->grafted_vtx->grafted_poly=NULL;
    }
    vtx_list_free(poly->vlist);
    bond_list_free(poly->blist);
    free(poly);

    return TS_SUCCESS;
}

ts_bool poly_list_free(ts_poly_list *poly_list){
    ts_uint i;

    for(i=0;i<poly_list->n;i++){
        poly_free(poly_list->poly[i]);
    }
    free(poly_list->poly);
    free(poly_list);
	
    return TS_SUCCESS;
}

 src/poly.h

New file
@@ -0,0 +1,16 @@
#ifndef _POLY_H
#define _POLY_H

#include"general.h"

ts_poly    *init_poly(ts_uint n, ts_vertex *grafted_vtx);

ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist);

ts_bool poly_free(ts_poly *poly);

ts_bool poly_list_free(ts_poly_list *poly_list);

ts_bool poly_assign_spring_const(ts_vesicle *vesicle);

#endif

 src/tape

@@ -1,13 +1,20 @@
####### Vesicle definitions ###########
# nshell is a number of divisions of dipyramid
nshell=17
# dmax is the square of the bond length 
dmax=1.67
# dmax is the max. bond length 
dmax=1.7
# bending rigidity of the membrane 
xk0=25.0
# max step size 
stepsize=0.15

####### Polymer definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
npoly=30
# nmono is a number of monomers in each polymer
nmono=10
# Spring constant between monomers of the polymer
k_spring=800

#######  Cell definitions ############
nxmax=60
@@ -17,11 +24,15 @@

####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
mcsweeps=100
mcsweeps=5000
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=100
inititer=1
#how many records do you want on the disk iteration are there in a run?
<<<<<<< HEAD
iterations=10
=======
iterations=10000
>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768


#shut up if we are using cluster!!!

 src/timestep.c

@@ -8,6 +8,7 @@
#include "bondflip.h"
#include "frame.h"
#include "io.h"

ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){
    ts_uint i, j;

@@ -20,6 +21,7 @@
        }
        centermass(vesicle);
        cell_occupation(vesicle);
        ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps);
        if(i>inititer){
            write_vertex_xml_file(vesicle,i-inititer);
        }
@@ -30,7 +32,7 @@
ts_bool single_timestep(ts_vesicle *vesicle){
    ts_bool retval;
    ts_double rnvec[3];
    ts_uint i, b;
    ts_uint i,j,b;
    for(i=0;i<vesicle->vlist->n;i++){
        rnvec[0]=drand48();
        rnvec[1]=drand48();
@@ -39,7 +41,7 @@
    }

//    ts_int cnt=0;
    for(i=0;i<vesicle->vlist->n;i++){
    for(i=0;i<3*vesicle->vlist->n;i++){
//why is rnvec needed in bondflip?
/*        rnvec[0]=drand48();
        rnvec[1]=drand48();
@@ -51,7 +53,18 @@
        retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
//    if(retval==TS_SUCCESS) cnt++;        
    } 
//    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);

    for(i=0;i<vesicle->poly_list->n;i++){
    for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
        rnvec[0]=drand48();
        rnvec[1]=drand48();
        rnvec[2]=drand48();
        retval=single_poly_vertex_move(vesicle,vesicle->poly_list->poly[i],vesicle->poly_list->poly[i]->vlist->vtx[j],rnvec);	
    }

    }
 
//    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,3*vesicle->blist->n,(double)cnt/(double)vesicle->blist->n/3.0);
    if(retval);
    return TS_SUCCESS;
}

 src/vertexmove.c

@@ -13,8 +13,7 @@
#include "vertexmove.h"
#include <string.h>

ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double
*rn){
ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
    ts_uint i;
    ts_double dist;
    ts_bool retval; 
@@ -58,6 +57,16 @@
        return TS_FAIL;
        }
    }

// Distance with grafted poly-vertex check:	
    if(vtx->grafted_poly!=NULL){
        dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
        if(dist<1.0 || dist>vesicle->dmax) {
        vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
        return TS_FAIL;
        }
    }

    //self avoidance check with distant vertices
     cellidx=vertex_self_avoidance(vesicle, vtx);
    //check occupation number
@@ -88,6 +97,13 @@
        energy_vertex(vtx->neigh->vtx[i]);
        delta_energy+=vtx->neigh->vtx[i]->xk*(vtx->neigh->vtx[i]->energy-oenergy);
    }
/* No poly-bond energy for now!
    if(vtx->grafted_poly!=NULL){
        delta_energy+=
            (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
            pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
    }
*/
//   fprintf(stderr, "DE=%f\n",delta_energy);
    //MONTE CARLOOOOOOOO
    if(delta_energy>=0){
@@ -126,3 +142,107 @@
    return TS_SUCCESS;
}


ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
    ts_uint i;
    ts_bool retval; 
    ts_uint cellidx; 
//    ts_double delta_energy;
    ts_double costheta,sintheta,phi,r;
    ts_double dist;
    //This will hold all the information of vtx and its neighbours
    ts_vertex backupvtx;
//    ts_bond backupbond[2];
    memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));

    //random move in a sphere with radius stepsize:
    r=vesicle->stepsize*rn[0];
    phi=rn[1]*2*M_PI;
    costheta=2*rn[2]-1;
    sintheta=sqrt(1-pow(costheta,2));
    vtx->x=vtx->x+r*sintheta*cos(phi);
    vtx->y=vtx->y+r*sintheta*sin(phi);
    vtx->z=vtx->z+r*costheta;


    //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,sizeof(ts_vertex));
            return TS_FAIL;
        }
    }

// Distance with grafted vesicle-vertex check:	
    if(vtx==poly->vlist->vtx[0]){
        dist=vtx_distance_sq(vtx,poly->grafted_vtx);
        if(dist<1.0 || dist>vesicle->dmax) {
        vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
        return TS_FAIL;
        }
    }


    //self avoidance check with distant vertices
    cellidx=vertex_self_avoidance(vesicle, vtx);
    //check occupation number
    retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	
    if(retval==TS_FAIL){
        vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
        return TS_FAIL;
    } 


    //if all the tests are successful, then energy for vtx and neighbours is calculated
/* Energy ignored for now!
    delta_energy=0;
    for(i=0;i<vtx->bond_no;i++){
        memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));

        vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
        bond_energy(vtx->bond[i],poly);
        delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
    }

    if(vtx==poly->vlist->vtx[0]){
        delta_energy+=
            (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
            pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
		
    }


    if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
        if(exp(-delta_energy)< drand48() )
#endif
#ifdef TS_DOUBLE_FLOAT
        if(expf(-delta_energy)< (ts_float)drand48())
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
        if(expl(-delta_energy)< (ts_ldouble)drand48())
#endif
        {
    //not accepted, reverting changes
    vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
    for(i=0;i<vtx->bond_no;i++){
    vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
    }

    return TS_FAIL; 
    }
    }
*/
		
//    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
    if(vtx->cell!=vesicle->clist->cell[cellidx]){
        retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
//        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
        if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);	
    }
//    if(oldcellidx);
    //END MONTE CARLOOOOOOO
    return TS_SUCCESS;
}

 src/vertexmove.h

@@ -1,3 +1,4 @@

ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn);

ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn);

 src/vesicle.c

@@ -5,10 +5,11 @@
#include "bond.h"
#include "cell.h"
#include "stdlib.h"
#include "poly.h"

ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
    ts_vesicle *vesicle=(ts_vesicle *)malloc(sizeof(ts_vesicle));
    ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
    vesicle->vlist=init_vertex_list(N);
    vesicle->blist=init_bond_list();
    vesicle->tlist=init_triangle_list();
@@ -33,6 +34,7 @@
    bond_list_free(vesicle->blist);
    triangle_list_free(vesicle->tlist);
    cell_list_free(vesicle->clist);
    poly_list_free(vesicle->poly_list);
    free(vesicle);
    return TS_SUCCESS;
}