Merge branch 'master' of bitbucket.org:samop/trisurf-ng

Samo Penic

2014-02-11 4e6285b2342038e3740ff9cea7c31364b4caa5e5

Merge branch 'master' of bitbucket.org:samop/trisurf-ng

 aclocal.m4

@@ -1,4 +1,4 @@
# generated automatically by aclocal 1.11.3 -*- Autoconf -*-
# generated automatically by aclocal 1.11.6 -*- 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.68],,
[m4_warning([this file was generated for autoconf 2.68.
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.69],,
[m4_warning([this file was generated for autoconf 2.69.
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.3], [],
m4_if([$1], [1.11.6], [],
      [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.3])dnl
[AM_AUTOMAKE_VERSION([1.11.6])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])

 src/bond.c

@@ -28,7 +28,7 @@
    //NOW insert vertices into data!    
    blist->bond[blist->n - 1]->vtx1=vtx1;    
    blist->bond[blist->n - 1]->vtx2=vtx2;

    blist->bond[blist->n - 1]->tainted=0;
    //Should we calculate bond length NOW?
    
    return blist->bond[blist->n-1];

 src/cell.c

@@ -53,13 +53,13 @@
    ncy=(ts_uint)((vtx->y-vesicle->cm[1])*clist->dcell+clist->shift);
    ncz=(ts_uint)((vtx->z-vesicle->cm[2])*clist->dcell+clist->shift);

    if(ncx == clist->ncmax[0]-1 || ncx == 2){
    if(ncx >= clist->ncmax[0]-1 || ncx <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500);
    }
    if(ncy == clist->ncmax[1]-1 || ncy == 2){
    if(ncy >= clist->ncmax[1]-1 || ncy <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate y is the problem.",1500);
    }
    if(ncz == clist->ncmax[2]-1 || ncz == 2){
    if(ncz >= clist->ncmax[2]-1 || ncz <= 2){
        fatal("Vesicle is positioned outside the cell covered area. Coordinate z is the problem.",1500);
    }
    cellidx=ncz+(ncy-1)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* 
@@ -69,16 +69,48 @@


//TODO: looks ok, but debug anyway in the future
ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
    ts_uint i;
    for(i=0;i<cell->nvertex;i++){
        if(cell->vertex[i]==vtx){
    //vertex is already in the cell!
            //fprintf(stderr,"VTX in the cell!\n");
            return TS_FAIL;
        }
    }
            //fprintf(stderr,"VTX added to the cell!\n");
    cell->nvertex++;
    cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
        if(cell->vertex == NULL){
            fatal("Reallocation of memory failed during insertion of vertex in cell_add_vertex",3);
        }
    cell->vertex[cell->nvertex-1]=vtx;
    vtx->cell=cell;
    return TS_SUCCESS;
}

inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx){
   ts_uint i,j=0;
    for(i=0;i<cell->nvertex;i++){
        if(cell->vertex[i]!=vtx){
            cell->vertex[j]=cell->vertex[i];
            j++;
        }
    }
    if(j==i){
    fatal("Vertex was not in the cell!",3);
    } 
    //fprintf(stderr, "Vertex deleted from the cell!\n");

/* resize memory. potentionally time consuming */
    cell->nvertex--;
    cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *));
    if(vtx->neigh == NULL && vtx->neigh_no!=0)
        if(cell->vertex == NULL){
            fatal("Reallocation of memory failed during removal of vertex in cell_remove_vertex",3);
        }
    return TS_SUCCESS;
}

ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist){
    ts_uint i;
@@ -93,7 +125,7 @@
}

// TODO: compiles ok, but it is completely untested and undebugged. It was debugged before rewrite, but this was long time ago.
ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx){
ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx){
    ts_uint ncx,ncy,ncz,remainder,cell_occupation;
    ts_uint i,j,k,l,neigh_cidx;
    ts_double dist;
@@ -117,17 +149,18 @@
// cell!
                if(cell_occupation>1){
                    for(l=0;l<cell_occupation;l++){
                        if(clist->cell[neigh_cidx]->vertex[l]!=vtx){

                //carefull with this checks!
                        if(clist->cell[neigh_cidx]->vertex[l]->idx!=vtx->idx){
                    //        fprintf(stderr,"calling dist on vertex %i\n",l);
                           dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],tvtx);
                           dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx);
                    //        fprintf(stderr,"dist was %f\n",dist);
                            if(dist<1) return TS_FAIL;
                            if(dist<=1.0) return TS_FAIL;
                        }
                    }
                }
            }
        }
    }
    
    } 
    return TS_SUCCESS;
}

 src/cell.h

@@ -4,8 +4,8 @@
ts_bool cell_free(ts_cell* cell);
ts_bool cell_list_free(ts_cell_list *clist);
inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx);
ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx);
inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx);
inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx);
ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist);

ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx);
inline ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx);
#endif

 src/energy.c

@@ -28,7 +28,7 @@
    ts_uint jjp,jjm;
    ts_vertex *j,*jp, *jm;
    ts_triangle *jt;
    ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0;
    ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0;
    ts_double x1,x2,x3,ctp,ctm,tot,xlen;
    ts_double h,ht;
    for(jj=1; jj<=vtx->neigh_no;jj++){
@@ -72,7 +72,9 @@
#endif
        tot=ctp+ctm;
        tot=0.5*tot;

        xlen=vtx_distance_sq(j,vtx);
/*
#ifdef  TS_DOUBLE_DOUBLE 
        vtx->bond[jj-1]->bond_length=sqrt(xlen); 
#endif
@@ -84,7 +86,7 @@
#endif

        vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;

*/
        s+=tot*xlen;
        xh+=tot*(j->x - vtx->x);
        yh+=tot*(j->y - vtx->y);

 src/frame.c

@@ -23,6 +23,10 @@
        vtx[i]->z-=vesicle->cm[2]; 
    } 

    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;

    return TS_SUCCESS;
}

@@ -37,7 +41,8 @@
    cell_list_cell_occupation_clear(vesicle->clist);
    for(i=0;i<n;i++){
    cellidx=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);
    vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];
//    already done in cell_add_vertex
//    vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];

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


 src/general.h

@@ -169,6 +169,7 @@
    ts_vertex *vtx2;
    ts_double bond_length;
    ts_double bond_length_dual;
    ts_bool tainted;
};
typedef struct ts_bond ts_bond;


 src/io.c

@@ -8,7 +8,7 @@
#include<stdlib.h>
#include <sys/types.h>
#include <dirent.h>

#include "initial_distribution.h"

ts_bool print_vertex_list(ts_vertex_list *vlist){
    ts_uint i;
@@ -294,12 +294,14 @@



ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations){
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!
    char buf[255];
    long int brezveze=1;
    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;
    *iterations=1000;
    long int iter=1000, init=1000, mcsw=1000;
    cfg_opt_t opts[] = {
        CFG_SIMPLE_INT("nshell", &nshell),
        CFG_SIMPLE_FLOAT("dmax", &dmax),
@@ -308,12 +310,14 @@
        CFG_SIMPLE_INT("nxmax", &ncxmax),
        CFG_SIMPLE_INT("nymax", &ncymax),
        CFG_SIMPLE_INT("nzmax", &nczmax),
        CFG_SIMPLE_INT("iterations",iterations),
        CFG_SIMPLE_INT("iterations",&iter),
    CFG_SIMPLE_INT("mcsweeps",&mcsw),
    CFG_SIMPLE_INT("inititer", &init),
        CFG_SIMPLE_BOOL("quiet",&quiet),
        CFG_SIMPLE_STR("multiprocessing",buf),
        CFG_SIMPLE_INT("smp_cores",&brezveze),
        CFG_SIMPLE_INT("cluster_nodes",&brezveze),
        CFG_SIMPLE_INT("distributed_processes",&brezveze),
        CFG_SIMPLE_INT("smp_cores",&brezveze0),
        CFG_SIMPLE_INT("cluster_nodes",&brezveze1),
        CFG_SIMPLE_INT("distributed_processes",&brezveze2),
        CFG_END()
    };
    cfg_t *cfg;    
@@ -326,6 +330,11 @@
    else if(retval==CFG_PARSE_ERROR){
    fatal("Invalid tape!",100);
    }
    ts_vesicle *vesicle;
        *iterations=iter;
    *inititer=init;
    *mcsweeps=mcsw;
    vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
    vesicle->nshell=nshell;
    vesicle->dmax=dmax*dmax;
    vesicle->bending_rigidity=xk0;
@@ -334,9 +343,11 @@
    vesicle->clist->ncmax[1]=ncymax;
    vesicle->clist->ncmax[2]=nczmax;
    vesicle->clist->max_occupancy=8;

    cfg_free(cfg);
//    fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
    return TS_SUCCESS;
    free(buf);
  //  fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
    return vesicle;

}


 src/io.h

@@ -48,7 +48,6 @@
ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text);
ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno);
ts_bool write_master_xml_file(ts_char *filename);
ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations);

ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations);

#endif

 src/main.c

@@ -18,14 +18,15 @@
*/

int main(int argv, char *argc[]){
ts_uint inititer,mcsweeps, iterations;
ts_vesicle *vesicle;
/* THIS SHOULD GO INTO UNIT TEST
ts_bool retval;
ts_uint i;
    ts_vertex_list *vlist=init_vertex_list(5);
ts_vertex_list *vlist1;
ts_bond_list *blist=init_bond_list();
ts_triangle_list *tlist=init_triangle_list();
ts_cell_list *clist=init_cell_list(3,3,3,0.3);
ts_vesicle *vesicle;

retval=vtx_add_cneighbour(blist,vlist->vtx[1],vlist->vtx[0]);
if(retval==TS_FAIL) printf("1. already a member or vertex is null!\n");
@@ -61,24 +62,10 @@

vtx_list_free(vlist1);
printf("Tests complete.\n");
vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
//parsetape(vesicle,&i);
*/
vesicle=parsetape(&mcsweeps, &inititer, &iterations);
run_simulation(vesicle, mcsweeps, inititer, iterations);

//these four must come from parsetype!
vesicle->dmax=1.67*1.67;
vesicle->stepsize=0.15;
vesicle->clist->max_occupancy=8;
vesicle->bending_rigidity=25.0;
fprintf(stderr,"xk=%f",vesicle->bending_rigidity);

centermass(vesicle);
cell_occupation(vesicle);
for(i=0;i<100;i++){
single_timestep(vesicle);
if(i%100==0){
write_vertex_xml_file(vesicle,i/100);
}
}
write_master_xml_file("test.pvd");
write_dout_fcompat_file(vesicle,"dout");
vesicle_free(vesicle);

 src/spherical_trisurf.c

@@ -19,11 +19,14 @@
*/
ts_bool saveAvgUlm2(ts_vesicle *vesicle);
int main(int argv, char *argc[]){
ts_uint i,j;
ts_uint i,j,k;
ts_vesicle *vesicle;
ts_double r0;
vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
//parsetape(vesicle,&i);

//similar to nmax in fortran code
ts_uint nmax;

//these four must come from parsetype!
vesicle->dmax=1.67*1.67;
@@ -33,6 +36,21 @@
//fprintf(stderr,"xk=%f",vesicle->bending_rigidity);

    centermass(vesicle);
cell_occupation(vesicle);

//test if the structure is internally organized into cells correctly 
ts_uint cind;
for(i=0;i<vesicle->vlist->n;i++){
    cind=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);

    if(vesicle->clist->cell[cind]==vesicle->vlist->vtx[i]->cell){
        //fprintf(stdout,"(T) Idx match!\n");
    } else {
        fprintf(stderr,"(T) ***** Idx don't match!\n");

    }
}
//end test
vesicle->sphHarmonics=sph_init(vesicle->vlist, 21);

vesicle_volume(vesicle);
@@ -42,14 +60,25 @@
calculateYlmi(vesicle);
calculateUlm(vesicle);



for(i=0;i<10;i++){
//preloop:
for(i=0;i<1000;i++){
    cell_occupation(vesicle);
    for(j=0;j<1000;j++){
        single_timestep(vesicle);
    }    
    centermass(vesicle);
    fprintf(stderr, "Preloop %d completed.\n",i+1);
}

nmax=1000;
for(i=0;i<nmax;i++){
    for(j=0;j<200;j++){
        cell_occupation(vesicle);
        for(k=0;k<5;k++){
        single_timestep(vesicle);
        }
        centermass(vesicle);
    }	
    vesicle_volume(vesicle);
    r0=getR0(vesicle);

@@ -61,8 +90,10 @@
    saveAvgUlm2(vesicle);

    write_vertex_xml_file(vesicle,i);
    fprintf(stderr, "Loop %d completed.\n",i+1);
    fprintf(stderr, "Loop %d out of %d completed.\n",i+1,nmax);

}

write_master_xml_file("test.pvd");
write_dout_fcompat_file(vesicle,"dout");
vesicle_free(vesicle);

 src/tape

@@ -16,8 +16,13 @@


####### Program Control ############
#how many iteration are there in a run?
iterations=20000
#how many MC sweeps between subsequent records of states to disk
mcsweeps=100
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=100
#how many records do you want on the disk iteration are there in a run?
iterations=1000


#shut up if we are using cluster!!!
quiet=false

 src/timestep.c

@@ -6,11 +6,31 @@
#include "timestep.h"
#include "vertexmove.h"
#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;

    centermass(vesicle);
    cell_occupation(vesicle);
    ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
    for(i=0;i<inititer+iterations;i++){
        for(j=0;j<mcsweeps;j++){
            single_timestep(vesicle);
        }
        centermass(vesicle);
        cell_occupation(vesicle);
        if(i>inititer){
            write_vertex_xml_file(vesicle,i-inititer);
        }
    }
    return TS_SUCCESS;
}

ts_bool single_timestep(ts_vesicle *vesicle){
    ts_bool retval;
    ts_double rnvec[3];
    ts_uint i;
    ts_uint i, b;
    for(i=0;i<vesicle->vlist->n;i++){
        rnvec[0]=drand48();
        rnvec[1]=drand48();
@@ -19,13 +39,16 @@
    }

//    ts_int cnt=0;
    for(i=0;i<vesicle->blist->n;i++){
        rnvec[0]=drand48();
    for(i=0;i<vesicle->vlist->n;i++){
//why is rnvec needed in bondflip?
/*        rnvec[0]=drand48();
        rnvec[1]=drand48();
        rnvec[2]=drand48();
*/ 
    b=rand() % vesicle->blist->n;
        //find a bond and return a pointer to a bond...
        //call single_bondflip_timestep...
        retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[i],rnvec);
        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);

 src/timestep.h

@@ -1,4 +1,5 @@
#ifndef _TIMESTEP_H
#define _TIMESTEP_H
ts_bool single_timestep(ts_vesicle *vesicle);
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations);
#endif

 src/vertexmove.c

@@ -11,6 +11,7 @@
//#include "io.h"
#include<stdio.h>
#include "vertexmove.h"
#include <string.h>

ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double
*rn){
@@ -18,50 +19,69 @@
    ts_double dist;
    ts_bool retval; 
    ts_uint cellidx; 
    ts_double xold,yold,zold;
    ts_double delta_energy,oenergy;
    ts_vertex *ovtx;
    ts_vertex *tvtx=(ts_vertex *)calloc(1,sizeof(ts_vertex));
    ts_double costheta,sintheta,phi,r;
    //This will hold all the information of vtx and its neighbours
    ts_vertex backupvtx[20];
    memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));

    //randomly we move the temporary vertex
    tvtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
    tvtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
    tvtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
    //check we if some length to neighbours are too much
    //Some stupid tests for debugging cell occupation!
/*         cellidx=vertex_self_avoidance(vesicle, vtx);
    if(vesicle->clist->cell[cellidx]==vtx->cell){
        fprintf(stderr,"Idx match!\n");
    } else {
        fprintf(stderr,"***** Idx don't match!\n");
        fatal("ENding.",1);
    }
*/

        //temporarly moving the vertex
//    vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
//        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:
    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(tvtx,vtx->neigh[i]);
        dist=vtx_distance_sq(vtx,vtx->neigh[i]);
        if(dist<1.0 || dist>vesicle->dmax) {
        vtx_free(tvtx);
//    fprintf(stderr,"Fail 1, dist=%f, vesicle->dmax=%f\n", 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, tvtx);
     cellidx=vertex_self_avoidance(vesicle, vtx);
    //check occupation number
     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx);
     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	
    if(retval==TS_FAIL){
    vtx_free(tvtx);
//    fprintf(stderr,"Fail 2\n");
        vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
        return TS_FAIL;
    } 
    
    //if all the tests are successful, then we update the vertex position
    xold=vtx->x;
    yold=vtx->y;
    zold=vtx->z;
    ovtx=malloc(sizeof(ts_vertex));
    vtx_copy(ovtx,vtx);
    vtx->x=tvtx->x;
    vtx->y=tvtx->y;
    vtx->z=tvtx->z;
   
 
    //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));
    }



    delta_energy=0;
    //update the normals of triangles that share bead i.
    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
    //energy and curvature
    oenergy=vtx->energy;
    energy_vertex(vtx);
    delta_energy=vtx->xk*(vtx->energy - ovtx->energy);
    delta_energy=vtx->xk*(vtx->energy - oenergy);
    //the same is done for neighbouring vertices
    for(i=0;i<vtx->neigh_no;i++){
        oenergy=vtx->neigh[i]->energy;
@@ -82,30 +102,27 @@
#endif
    {
    //not accepted, reverting changes
    vtx->x=xold;
    vtx->y=yold;
    vtx->z=zold;
    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));
    }
	
    //update the normals of triangles that share bead i.
    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
    //energy and curvature
    energy_vertex(vtx);
    //the same is done for neighbouring vertices
    for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]);
    free(ovtx->bond_length);
    free(ovtx->bond_length_dual);
    free(ovtx);
    vtx_free(tvtx);
   for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);

    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[0].cell,vtx);
		
    }
//    if(oldcellidx);
    //END MONTE CARLOOOOOOO

    //TODO: change cell occupation if necessary!
//    fprintf(stderr,"Success!!\n");
    free(ovtx->bond_length);
    free(ovtx->bond_length_dual);
    free(ovtx);
    vtx_free(tvtx);
    return TS_SUCCESS;
}