Trisurf Monte Carlo simulator
Samo Penic
2014-03-21 ff8152b46e7957beffeafca580abc90df47d9d28
Merge branch 'trisurf-polyel'
15 files modified
799 ■■■■■ changed files
missing 179 ●●●●● patch | view | raw | blame | history
src/bond.c 11 ●●●●● patch | view | raw | blame | history
src/bond.h 1 ●●●● patch | view | raw | blame | history
src/frame.c 36 ●●●●● patch | view | raw | blame | history
src/general.h 15 ●●●● patch | view | raw | blame | history
src/initial_distribution.c 34 ●●●●● patch | view | raw | blame | history
src/io.c 235 ●●●● patch | view | raw | blame | history
src/io.h 9 ●●●● patch | view | raw | blame | history
src/main.c 15 ●●●● patch | view | raw | blame | history
src/poly.c 96 ●●●● patch | view | raw | blame | history
src/poly.h 5 ●●●● patch | view | raw | blame | history
src/tape 26 ●●●● patch | view | raw | blame | history
src/timestep.c 22 ●●●● patch | view | raw | blame | history
src/vertexmove.c 113 ●●●●● patch | view | raw | blame | history
src/vertexmove.h 2 ●●●●● patch | view | raw | blame | history
missing
@@ -1,12 +1,4 @@
#! /bin/sh
<<<<<<< HEAD
# Common wrapper for a few potentially missing GNU programs.
scriptversion=2012-06-26.16; # UTC
# Copyright (C) 1996-2013 Free Software Foundation, Inc.
# Originally written by Fran,cois Pinard <pinard@iro.umontreal.ca>, 1996.
=======
# Common stub for a few missing GNU programs while installing.
scriptversion=2012-01-06.13; # UTC
@@ -14,7 +6,6 @@
# Copyright (C) 1996, 1997, 1999, 2000, 2002, 2003, 2004, 2005, 2006,
# 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
# Originally by Fran,cois Pinard <pinard@iro.umontreal.ca>, 1996.
>>>>>>> dump-state
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -35,24 +26,6 @@
# the same distribution terms that you use for the rest of that program.
if test $# -eq 0; then
<<<<<<< HEAD
  echo 1>&2 "Try '$0 --help' for more information"
  exit 1
fi
case $1 in
  --is-lightweight)
    # Used by our autoconf macros to check whether the available missing
    # script is modern enough.
    exit 0
    ;;
  --run)
    # Back-compat with the calling convention used by older automake.
    shift
    ;;
=======
  echo 1>&2 "Try \`$0 --help' for more information"
  exit 1
fi
@@ -87,32 +60,17 @@
    msg="probably too old"
  fi
  ;;
>>>>>>> dump-state
  -h|--h|--he|--hel|--help)
    echo "\
$0 [OPTION]... PROGRAM [ARGUMENT]...
<<<<<<< HEAD
Run 'PROGRAM [ARGUMENT]...', returning a proper advice when this fails due
to PROGRAM being missing or too old.
=======
Handle \`PROGRAM [ARGUMENT]...' for when PROGRAM is missing, or return an
error status if there is no known handling for PROGRAM.
>>>>>>> dump-state
Options:
  -h, --help      display this help and exit
  -v, --version   output version information and exit
<<<<<<< HEAD
Supported PROGRAM values:
  aclocal   autoconf  autoheader   autom4te  automake  makeinfo
  bison     yacc      flex         lex       help2man
Version suffixes to PROGRAM as well as the prefixes 'gnu-', 'gnu', and
'g' are ignored when checking the name.
=======
  --run           try to run the given command, and emulate it if it fails
Supported PROGRAM values:
@@ -130,7 +88,6 @@
Version suffixes to PROGRAM as well as the prefixes \`gnu-', \`gnu', and
\`g' are ignored when checking the name.
>>>>>>> dump-state
Send bug reports to <bug-automake@gnu.org>."
    exit $?
@@ -142,148 +99,13 @@
    ;;
  -*)
<<<<<<< HEAD
    echo 1>&2 "$0: unknown '$1' option"
    echo 1>&2 "Try '$0 --help' for more information"
=======
    echo 1>&2 "$0: Unknown \`$1' option"
    echo 1>&2 "Try \`$0 --help' for more information"
>>>>>>> dump-state
    exit 1
    ;;
esac
<<<<<<< HEAD
# Run the given program, remember its exit status.
"$@"; st=$?
# If it succeeded, we are done.
test $st -eq 0 && exit 0
# Also exit now if we it failed (or wasn't found), and '--version' was
# passed; such an option is passed most likely to detect whether the
# program is present and works.
case $2 in --version|--help) exit $st;; esac
# Exit code 63 means version mismatch.  This often happens when the user
# tries to use an ancient version of a tool on a file that requires a
# minimum version.
if test $st -eq 63; then
  msg="probably too old"
elif test $st -eq 127; then
  # Program was missing.
  msg="missing on your system"
else
  # Program was found and executed, but failed.  Give up.
  exit $st
fi
perl_URL=http://www.perl.org/
flex_URL=http://flex.sourceforge.net/
gnu_software_URL=http://www.gnu.org/software
program_details ()
{
  case $1 in
    aclocal|automake)
      echo "The '$1' program is part of the GNU Automake package:"
      echo "<$gnu_software_URL/automake>"
      echo "It also requires GNU Autoconf, GNU m4 and Perl in order to run:"
      echo "<$gnu_software_URL/autoconf>"
      echo "<$gnu_software_URL/m4/>"
      echo "<$perl_URL>"
      ;;
    autoconf|autom4te|autoheader)
      echo "The '$1' program is part of the GNU Autoconf package:"
      echo "<$gnu_software_URL/autoconf/>"
      echo "It also requires GNU m4 and Perl in order to run:"
      echo "<$gnu_software_URL/m4/>"
      echo "<$perl_URL>"
      ;;
  esac
}
give_advice ()
{
  # Normalize program name to check for.
  normalized_program=`echo "$1" | sed '
    s/^gnu-//; t
    s/^gnu//; t
    s/^g//; t'`
  printf '%s\n' "'$1' is $msg."
  configure_deps="'configure.ac' or m4 files included by 'configure.ac'"
  case $normalized_program in
    autoconf*)
      echo "You should only need it if you modified 'configure.ac',"
      echo "or m4 files included by it."
      program_details 'autoconf'
      ;;
    autoheader*)
      echo "You should only need it if you modified 'acconfig.h' or"
      echo "$configure_deps."
      program_details 'autoheader'
      ;;
    automake*)
      echo "You should only need it if you modified 'Makefile.am' or"
      echo "$configure_deps."
      program_details 'automake'
      ;;
    aclocal*)
      echo "You should only need it if you modified 'acinclude.m4' or"
      echo "$configure_deps."
      program_details 'aclocal'
      ;;
   autom4te*)
      echo "You might have modified some maintainer files that require"
      echo "the 'automa4te' program to be rebuilt."
      program_details 'autom4te'
      ;;
    bison*|yacc*)
      echo "You should only need it if you modified a '.y' file."
      echo "You may want to install the GNU Bison package:"
      echo "<$gnu_software_URL/bison/>"
      ;;
    lex*|flex*)
      echo "You should only need it if you modified a '.l' file."
      echo "You may want to install the Fast Lexical Analyzer package:"
      echo "<$flex_URL>"
      ;;
    help2man*)
      echo "You should only need it if you modified a dependency" \
           "of a man page."
      echo "You may want to install the GNU Help2man package:"
      echo "<$gnu_software_URL/help2man/>"
    ;;
    makeinfo*)
      echo "You should only need it if you modified a '.texi' file, or"
      echo "any other file indirectly affecting the aspect of the manual."
      echo "You might want to install the Texinfo package:"
      echo "<$gnu_software_URL/texinfo/>"
      echo "The spurious makeinfo call might also be the consequence of"
      echo "using a buggy 'make' (AIX, DU, IRIX), in which case you might"
      echo "want to install GNU make:"
      echo "<$gnu_software_URL/make/>"
      ;;
    *)
      echo "You might have modified some files without having the proper"
      echo "tools for further handling them.  Check the 'README' file, it"
      echo "often tells you about the needed prerequisites for installing"
      echo "this package.  You may also peek at any GNU archive site, in"
      echo "case some other package contains this missing '$1' program."
      ;;
  esac
}
give_advice "$1" | sed -e '1s/^/WARNING: /' \
                       -e '2,$s/^/         /' >&2
# Propagate the correct exit status (expected to be 127 for a program
# not found, 63 for a program that failed due to version mismatch).
exit $st
=======
# normalize program name to check for.
program=`echo "$1" | sed '
  s/^gnu-//; t
@@ -499,7 +321,6 @@
esac
exit 0
>>>>>>> dump-state
# Local variables:
# eval: (add-hook 'write-file-hooks 'time-stamp)
src/bond.c
@@ -35,6 +35,17 @@
    return blist->bond[blist->n-1];
}
ts_bool bond_vector(ts_bond *bond){
    bond->x = bond->vtx1->x - bond->vtx2->x;
    bond->y = bond->vtx1->y - bond->vtx2->y;
    bond->z = bond->vtx1->z - bond->vtx2->z;
    return TS_SUCCESS;
}
ts_bool bond_list_free(ts_bond_list *blist){
    ts_uint i;
    for(i=0;i<blist->n;i++){
src/bond.h
@@ -20,6 +20,7 @@
 */
ts_bond *bond_add(ts_bond_list *blist, ts_vertex *vtx1, ts_vertex *vtx2);
ts_bool bond_vector(ts_bond *bond);
ts_bool bond_list_free(ts_bond_list *blist);
src/frame.c
@@ -22,7 +22,7 @@
        vtx[i]->y-=vesicle->cm[1];
        vtx[i]->z-=vesicle->cm[2]; 
    } 
//move polymers for the same vector as we moved vesicle
    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];
@@ -30,22 +30,24 @@
            vesicle->poly_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
        }
    }
//move filaments for the same vector as we moved vesicle
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            vesicle->filament_list->poly[i]->vlist->vtx[j]->x-=vesicle->cm[0];
            vesicle->filament_list->poly[i]->vlist->vtx[j]->y-=vesicle->cm[1];
            vesicle->filament_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
        }
    }
    vesicle->cm[0]=0;
    vesicle->cm[1]=0;
    vesicle->cm[2]=0;
    vesicle->cm[0]=0.0;
    vesicle->cm[1]=0.0;
    vesicle->cm[2]=0.0;
    return TS_SUCCESS;
}
ts_bool cell_occupation(ts_vesicle *vesicle){
    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);
    for(i=0;i<n;i++){
@@ -56,17 +58,21 @@
    cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);
    }
//Add all polymers to cells
    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]);
    }
    }
//Add all filaments to cells
     for(i=0;i<vesicle->filament_list->n;i++){
    for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
        cellidx=vertex_self_avoidance(vesicle, vesicle->filament_list->poly[i]->vlist->vtx[j]);
        cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->filament_list->poly[i]->vlist->vtx[j]);
    }
    }
    //fprintf(stderr, "Bil sem tu\n");
    //if(dcell);
    //if(shift);
    return TS_SUCCESS;
}
src/general.h
@@ -166,10 +166,11 @@
        ts_uint idx;
    ts_vertex *vtx1;
    ts_vertex *vtx2;
    ts_double bond_length;
    ts_double bond_length_dual;
    ts_bool tainted;
        ts_double bond_length;
        ts_double bond_length_dual;
    ts_bool tainted; //TODO: remove
    ts_double energy;
    ts_double x,y,z;
};
typedef struct ts_bond ts_bond;
@@ -256,11 +257,17 @@
       ts_double cm[3];
    ts_double volume;
    ts_spharm *sphHarmonics;
// Polymers outside the vesicle and attached to the vesicle membrane (polymer brush):
    ts_poly_list *poly_list;
// Filaments inside the vesicle (not attached to the vesicel membrane:
    ts_poly_list *filament_list;
    ts_double spring_constant;
    ts_double pressure;
    ts_int pswitch;
    ts_double R_nucleus;
} ts_vesicle;
src/initial_distribution.c
@@ -37,10 +37,42 @@
ts_vesicle *create_vesicle_from_tape(ts_tape *tape){
    ts_vesicle *vesicle;
    ts_vertex *vtx;
    vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
    vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist);
    // Nucleus:
    vesicle->R_nucleus=tape->R_nucleus;
    //Initialize grafted polymers (brush):
    vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle);
    vesicle->spring_constant=tape->kspring;
    poly_assign_spring_const(vesicle);
    //Initialize filaments (polymers inside the vesicle):
    vesicle->filament_list=init_poly_list(tape->nfil,tape->nfono, NULL, vesicle);
    poly_assign_filament_xi(vesicle,tape);
    ts_uint i,j;
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
            bond_vector(vesicle->filament_list->poly[i]->blist->bond[j]);
            vesicle->filament_list->poly[i]->blist->bond[j]->bond_length = sqrt(vtx_distance_sq(vesicle->filament_list->poly[i]->blist->bond[j]->vtx1,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2));
        }
    }
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            vtx = vesicle->filament_list->poly[i]->vlist->vtx[j];
            if(vtx->bond_no == 2){
            vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
            }
        }
    }
//    vesicle->spring_constant=tape->kspring;
//    poly_assign_spring_const(vesicle);
    
    vesicle->nshell=tape->nshell;
    vesicle->dmax=tape->dmax*tape->dmax; /* dmax^2 in the vesicle dmax variable */
src/io.c
@@ -36,6 +36,8 @@
    fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
    /* dump poly list */
    fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
    /* dump filament list */
    fwrite(vesicle->filament_list, sizeof(ts_poly_list),1,fh);
    /* level 1 complete */
    /*dump vertices*/
@@ -118,6 +120,43 @@
        }
    }
  /*dump filamentes grandes svinjas */
    for(i=0;i<vesicle->filament_list->n;i++){
        fwrite(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh);
        fwrite(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
        fwrite(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
    }
    /* dump filamentes vertex(monomer) list*/
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            fwrite(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
            /* dump offset for neigh and bond */
            for(k=0;k<vesicle->filament_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(&vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]->idx,sizeof(ts_uint),1,fh);
            }
            for(k=0;k<vesicle->filament_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(&vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]->idx,sizeof(ts_uint),1,fh);
            }
        }
    }
    /* dump poly bonds between monomers list*/
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
            fwrite(vesicle->filament_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(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx,sizeof(ts_uint),1,fh);
//            off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
            fwrite(&vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx,sizeof(ts_uint),1,fh);
        }
    }
/* pointer offsets for fixing the restored pointers */
/* need pointers for 
    vlist->vtx
@@ -142,6 +181,20 @@
ts_vesicle *restore_state(ts_uint *iteration){
    ts_uint i,j,k;
    FILE *fh=fopen("dump.bin","rb");
    struct stat sb;
    if (stat("dump.bin", &sb) == -1) {
        //dump file does not exist.
        return NULL;
    }
    //check if it is regular file
    if((sb.st_mode & S_IFMT) != S_IFREG) {
        //dump file is not a regular file.
        ts_fprintf(stderr,"Dump file is not a regular file!\n");
        return NULL;
    }
    ts_uint retval;
    ts_uint idx;
@@ -166,6 +219,9 @@
    /* 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);
    /* restore filament list */
    vesicle->filament_list=(ts_poly_list *)calloc(1,sizeof(ts_poly_list));
    retval=fread(vesicle->filament_list, sizeof(ts_poly_list),1,fh);
    /* level 1 complete */
/* prerequisity. Bonds must be malloced before vertexes are recreated */
@@ -308,6 +364,61 @@
        }
    }
    /*restore filaments */
    vesicle->filament_list->poly = (ts_poly **)calloc(vesicle->filament_list->n,sizeof(ts_poly *));
    for(i=0;i<vesicle->filament_list->n;i++){
        vesicle->filament_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
        retval=fread(vesicle->filament_list->poly[i],sizeof(ts_poly),1,fh);
        vesicle->filament_list->poly[i]->vlist=(ts_vertex_list *)calloc(1,sizeof(ts_vertex_list));
        retval=fread(vesicle->filament_list->poly[i]->vlist,sizeof(ts_vertex_list),1,fh);
        vesicle->filament_list->poly[i]->blist=(ts_bond_list *)calloc(1,sizeof(ts_bond_list));
        retval=fread(vesicle->filament_list->poly[i]->blist,sizeof(ts_bond_list),1,fh);
    /* initialize adress space for pointers that will hold specific vertices (monomers) and bonds */
        vesicle->filament_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->n,sizeof(ts_vertex *));
        vesicle->filament_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->blist->n,sizeof(ts_bond *));
     for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            vesicle->filament_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
    }
    for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
            vesicle->filament_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
    }
    }
    /* restore poly vertex(monomer) list*/
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            retval=fread(vesicle->filament_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
            /* restore neigh and bonds */
            vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *));
            for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
                retval=fread(&idx,sizeof(ts_uint),1,fh);
                vesicle->filament_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->filament_list->poly[i]->vlist->vtx[idx];
            }
            vesicle->filament_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *));
            for(k=0;k<vesicle->filament_list->poly[i]->vlist->vtx[j]->bond_no;k++){
                retval=fread(&idx,sizeof(ts_uint),1,fh);
                vesicle->filament_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->filament_list->poly[i]->blist->bond[idx];
            }
        }
    }
    /* restore poly bonds between monomers list*/
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
       //     vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
            retval=fread(vesicle->filament_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
            /* restore vtx1 and vtx2 */
                retval=fread(&idx,sizeof(ts_uint),1,fh);
                vesicle->filament_list->poly[i]->blist->bond[j]->vtx1=vesicle->filament_list->poly[i]->vlist->vtx[idx];
                retval=fread(&idx,sizeof(ts_uint),1,fh);
                vesicle->filament_list->poly[i]->blist->bond[j]->vtx2=vesicle->filament_list->poly[i]->vlist->vtx[idx];
        }
    }
// recreating space for cells // 
    vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
    retval=fread(vesicle->clist, sizeof(ts_cell_list), 1,fh); 
@@ -326,9 +437,14 @@
ts_bool parse_args(int argc, char **argv){
 int c, retval;
 DIR *dir;
    path[0]=0;
    int c, retval;
    struct stat sb;
    sprintf(command_line_args.path, "./"); //clear string;
    sprintf(command_line_args.output_fullfilename,"./output.pvd");
    sprintf(command_line_args.dump_fullfilename,"./dump.bin");
    sprintf(command_line_args.tape_fullfilename,"./tape");
            FILE *file;
while (1)
     {
       static struct option long_options[] =
@@ -338,13 +454,13 @@
           {"tape",     no_argument,       0, 't'},
           {"output-file",  required_argument, 0, 'o'},
           {"directory",  required_argument, 0, 'd'},
           {"dump-file", required_argument,0, 'f'},
           {"dump-filename", required_argument,0, 'f'},
           {0, 0, 0, 0}
         };
       /* getopt_long stores the option index here. */
       int option_index = 0;
       c = getopt_long (argc, argv, "d:fot",
       c = getopt_long (argc, argv, "d:f:o:t:",
                        long_options, &option_index);
       /* Detect the end of the options. */
@@ -365,46 +481,56 @@
         case 't':
            //check if tape exists. If not, fail immediately.
           puts ("option -t\n");
            if (stat(optarg, &sb) == -1) {
                ts_fprintf(stderr,"Tape '%s' does not exist!\n",optarg);
                fatal("Please select correct tape",1);
            } else {
                strcpy(command_line_args.tape_fullfilename,optarg);
            }
           break;
         case 'o':
            //set filename of master output file
           printf ("option -o with value `%s'\n", optarg);
           break;
            if ((file = fopen(optarg, "w")) == NULL) {
                fprintf(stderr,"Could not create output file!\n");
                fatal("Please specify correct output file",1);
            } else {
                fclose(file);
            }
            strcpy(command_line_args.output_fullfilename, optarg);
            break;
         case 'd':
            //check if directory exists. If not create one. If creation is
            //successful, set directory for output files.
            //printf ("option -d with value `%s'\n", optarg);
            dir = opendir(optarg);
            if (dir)
            {
                /* Directory exists. */
                closedir(dir);
            }
            else if (ENOENT == errno)
            {
                /* Directory does not exist. */
            if (stat(optarg, &sb) == -1) {
                //directory does not exist
                retval=mkdir(optarg, 0700);
                if(retval){
                fatal("Could not create requested directory. Check if you have permissions",1);
                    fatal("Could not create requested directory. Check if you have permissions",1);
                }
            }
            else
            {
                /* opendir() failed for some other reason. */
                fatal("Could not check if directory exists. Reason unknown",1);
            //check if is a proper directory
            else if((sb.st_mode & S_IFMT) != S_IFDIR) {
                //it is not a directory. fatal error.
                ts_fprintf(stderr,"%s is not a directory!\n",optarg);
                fatal("Cannot continue",1);
            }
            ts_fprintf(stdout,"\n*** Using output directory: %s\n\n", optarg);
//            sprintf(path,"%s", optarg);
            strcpy(path, optarg);
           // ts_fprintf(stdout,"ok!\n");
            strcpy(command_line_args.path, optarg);
           break;
        case 'f':
            //check if dump file specified exists. Defaults to dump.bin
            if ((file = fopen(optarg, "w")) == NULL) {
                fprintf(stderr,"Could not create dump file!\n");
                fatal("Please specify correct dump file",1);
            } else {
                fclose(file);
            }
            strcpy(command_line_args.dump_fullfilename, optarg);
            break;
         case '?':
@@ -635,8 +761,8 @@
    /* Here comes header of the file */
    //find number of extra vtxs and bonds of polymeres
    ts_uint monono=0, polyno=0, poly_idx=0;
    ts_bool poly=0;
    ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
    ts_bool poly=0, fil=0;
    if(vesicle->poly_list!=NULL){
        if(vesicle->poly_list->poly[0]!=NULL){
        polyno=vesicle->poly_list->n;
@@ -644,8 +770,17 @@
        poly=1;
        }
    }
    if(vesicle->filament_list!=NULL){
        if(vesicle->filament_list->poly[0]!=NULL){
        filno=vesicle->filament_list->n;
        fonono=vesicle->filament_list->poly[0]->vlist->n;
        fil=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+monono*polyno, blist->n+monono*polyno);
    fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1));
    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);
@@ -655,6 +790,16 @@
        poly_idx=vlist->n;
        for(i=0;i<vesicle->poly_list->n;i++){
            for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++,poly_idx++){
                fprintf(fh,"%u ", poly_idx);
            }
        }
    }
    //filaments
    if(fil){
        poly_idx=vlist->n+monono*polyno;
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++,poly_idx++){
    //    fprintf(stderr,"was here\n");
                fprintf(fh,"%u ", poly_idx);
            }
        }
@@ -669,6 +814,14 @@
        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 );
            }
        }
    }
    //filaments
    if(fil){
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
                fprintf(fh,"%e %e %e\n", vesicle->filament_list->poly[i]->vlist->vtx[j]->x,vesicle->filament_list->poly[i]->vlist->vtx[j]->y, vesicle->filament_list->poly[i]->vlist->vtx[j]->z );
            }
        }
    }
@@ -691,14 +844,26 @@
    }
    
    //filaments
    if(fil){
        poly_idx=vlist->n+monono*polyno;
        for(i=0;i<vesicle->filament_list->n;i++){
            for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){
                fprintf(fh,"%u %u\n", vesicle->filament_list->poly[i]->blist->bond[j]->vtx1->idx+vlist->n+monono*polyno+i*fonono,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2->idx+vlist->n+monono*polyno+i*fonono);
//        fprintf(stderr,"was here\n");
            }
        }
    }
    fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
    for (i=2;i<(blist->n+monono*polyno)*2+1;i+=2){
    for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*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+monono*polyno;i++){
     for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){
        fprintf(fh,"3 ");
    }
@@ -768,11 +933,15 @@
        CFG_SIMPLE_INT("nshell", &tape->nshell),
        CFG_SIMPLE_INT("npoly", &tape->npoly),
        CFG_SIMPLE_INT("nmono", &tape->nmono),
        CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
    CFG_SIMPLE_INT("nfil",&tape->nfil),
    CFG_SIMPLE_INT("nfono",&tape->nfono),
    CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
    CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
        CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
    CFG_SIMPLE_INT("pswitch",&tape->pswitch),
    CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
    CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
    CFG_SIMPLE_FLOAT("xi",&tape->xi),
        CFG_SIMPLE_FLOAT("stepsize",&tape->stepsize),
        CFG_SIMPLE_INT("nxmax", &tape->ncxmax),
        CFG_SIMPLE_INT("nymax", &tape->ncymax),
src/io.h
@@ -6,7 +6,6 @@
static char prefixname[1024];
static ts_bool restore=0;
static char tape[1024]; */
char path[1024];
int force_from_tape;
@@ -17,6 +16,9 @@
    long int nczmax;
    long int npoly;
    long int nmono;
    long int nfil;
    long int nfono;
    long int R_nucleus;
    long int pswitch;
        char *multiprocessing;
       long int brezveze0;
@@ -26,6 +28,7 @@
    ts_double dmax;
    ts_double stepsize;
    ts_double kspring;
    ts_double xi;
    ts_double pressure;
    long int iterations;
    long int inititer;
@@ -36,6 +39,10 @@
typedef struct{
    ts_int force_from_tape;
    ts_int reset_iteration_count;
    char path[1024]; //path where all files should be added
    char output_fullfilename[1024]; //name of the master file
    char dump_fullfilename[1024]; //name of the dump file
    char tape_fullfilename[1024]; //name of the tape file
} ts_args;
ts_args command_line_args;
src/main.c
@@ -22,9 +22,10 @@
    ts_vesicle *vesicle;
    ts_tape *tape;
    ts_uint start_iteration=0;
    force_from_tape=0;
    parse_args(argv,argc); // sets global variable command_line_args (defined in io.h)
    ts_fprintf(stdout,"Starting program...\n\n");
    if(force_from_tape){
    if(command_line_args.force_from_tape){
        ts_fprintf(stdout,"************************************************\n");
        ts_fprintf(stdout,"**** Generating initial geometry from tape *****\n");
        ts_fprintf(stdout,"************************************************\n\n");
@@ -37,8 +38,18 @@
        ts_fprintf(stdout,"**********************************************************************\n\n");
        tape=parsetape("tape");
        vesicle=restore_state(&start_iteration);
        if(vesicle==NULL){
            ts_fprintf(stderr, "Dump file does not exist or is not a regular file! Did you mean to invoke trisurf with --force-from-tape option?\n\n");
            return 1;
        }
        // nove vrednosti iz tapea...
        vesicle->bending_rigidity=tape->xk0;
        vtx_set_global_values(vesicle);
        vesicle->pressure=tape->pressure;
        vesicle->dmax=tape->dmax*tape->dmax;
        poly_assign_filament_xi(vesicle,tape);
        if(command_line_args.reset_iteration_count) start_iteration=0;
        if(command_line_args.reset_iteration_count) start_iteration=tape->inititer-1;
        else start_iteration++;
        if(start_iteration>=tape->iterations){
src/poly.c
@@ -6,6 +6,17 @@
#include<math.h>
#include"energy.h"
ts_bool poly_assign_filament_xi(ts_vesicle *vesicle, ts_tape *tape){
    ts_uint i;
    for(i=0;i<vesicle->filament_list->n;i++){
     vesicle->filament_list->poly[i]->k = tape->xi;
        }
    return TS_SUCCESS;
}
ts_bool poly_assign_spring_const(ts_vesicle *vesicle){
    ts_uint i;
@@ -20,8 +31,10 @@
    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;
    if (grafted_vtx!=NULL){
        poly->grafted_vtx = grafted_vtx;
        grafted_vtx->grafted_poly = poly;
    }
    ts_uint i;
    for(i=0;i<n-1;i++){
@@ -38,45 +51,72 @@
}
ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist){
ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist, ts_vesicle *vesicle){
    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;
    ts_double dphi,dh;
    ts_uint ji;
    if (n_poly > vlist->n){fatal("Number of polymers larger then numbero f vertices on a vesicle.",310);}
    // Grafting polymers:
    if (vlist!=NULL){
        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++;
        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++;
            }
        }
    }
    else
    {
        for(i=0;i<n_poly;i++){
            poly_list->poly[i] = init_poly(n_mono, NULL);
        }
    }
    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++){
    if (vlist!=NULL){
    /* Make straight grafted poylmers normal to membrane (polymer brush). 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;
        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;
            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);
            }
        }
        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);
    }
    else
    {
    /* Make filaments inside the vesicle. Helix with radius... Dist. between poly vertices put to 1*/
        dphi = 2.0*asin(1.0/2.0/vesicle->R_nucleus)*1.001;
        dh = dphi/2.0/M_PI*1.001;
        for(i=0;i<poly_list->n;i++){
            for (j=0;j<poly_list->poly[i]->vlist->n;j++){
                ji = j + i*poly_list->poly[i]->vlist->n;
                poly_list->poly[i]->vlist->vtx[j]->x = vesicle->R_nucleus*cos(ji*dphi);
                poly_list->poly[i]->vlist->vtx[j]->y = vesicle->R_nucleus*sin(ji*dphi);
                poly_list->poly[i]->vlist->vtx[j]->z = ji*dh - (dh*poly_list->n*poly_list->poly[i]->vlist->n/2.0);
            }
        }
    }
src/poly.h
@@ -1,11 +1,12 @@
#ifndef _POLY_H
#define _POLY_H
#include"io.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_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist, ts_vesicle *vesicle);
ts_bool poly_free(ts_poly *poly);
@@ -13,4 +14,6 @@
ts_bool poly_assign_spring_const(ts_vesicle *vesicle);
ts_bool poly_assign_filament_xi(ts_vesicle *vesicle, ts_tape *tape);
#endif
src/tape
@@ -4,23 +4,35 @@
# 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=10.0
# max step size (in units l_min)
stepsize=0.15
# Pressure calculations
# (pswitch=1: calc. p*dV energy contribution)
pswitch = 0
# pressure difference: p_inside - p_outside (in units l_min^3/kT):
pressure=0.0
# pressure difference: p_inside - p_outside (in units kT/l_min^3):
pressure=0.1
####### Polymer definitions ###########
####### Polymer (brush) definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
npoly=20
npoly=300
# nmono is a number of monomers in each polymer
nmono=10
# Spring constant between monomers of the polymer
k_spring=800
####### Filament (inside the vesicle) definitions ###########
# nfil is a number of filaments inside the vesicle
nfil=1
# nfono is a number of monomers in each filament
nfono=300
# Persistence lenght of the filaments (in units l_min)
xi=100
####### Nucleus (inside the vesicle) ###########
# Radius of an impenetrable hard sphere inside the vesicle
R_nucleus=5
#######  Cell definitions ############
nxmax=60
@@ -32,9 +44,9 @@
#how many MC sweeps between subsequent records of states to disk
mcsweeps=500
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=1
inititer=0
#how many records do you want on the disk iteration are there in a run?
iterations=10
iterations=20
#shut up if we are using cluster!!!
src/timestep.c
@@ -23,7 +23,7 @@
        cell_occupation(vesicle);
        ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps);
            dump_state(vesicle,i);
        if(i>inititer){
        if(i>=inititer){
            write_vertex_xml_file(vesicle,i-inititer);
        }
    }
@@ -56,15 +56,25 @@
    }
    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);
        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);
        }
    }
    for(i=0;i<vesicle->filament_list->n;i++){
        for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
            rnvec[0]=drand48();
            rnvec[1]=drand48();
            rnvec[2]=drand48();
            retval=single_filament_vertex_move(vesicle,vesicle->filament_list->poly[i],vesicle->filament_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
@@ -254,3 +254,116 @@
    //END MONTE CARLOOOOOOO
    return TS_SUCCESS;
}
ts_bool single_filament_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[2];
    //This will hold all the information of vtx and its neighbours
    ts_vertex backupvtx,backupneigh[2];
    ts_bond backupbond[2];
    //backup vertex:
    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->bond_no;i++){
        dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
        if(dist[i]<1.0 || dist[i]>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;
    }
    //backup bonds
    for(i=0;i<vtx->bond_no;i++){
        memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
        vtx->bond[i]->bond_length=sqrt(dist[i]);
        bond_vector(vtx->bond[i]);
    }
    //backup neighboring vertices:
    for(i=0;i<vtx->neigh_no;i++){
        memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
    }
    //if all the tests are successful, then energy for vtx and neighbours is calculated
    delta_energy=0;
    if(vtx->bond_no == 2){
        vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
        delta_energy += vtx->energy - backupvtx.energy;
    }
    for(i=0;i<vtx->neigh_no;i++){
        if(vtx->neigh[i]->bond_no == 2){
            vtx->neigh[i]->energy = -(vtx->neigh[i]->bond[0]->x*vtx->neigh[i]->bond[1]->x + vtx->neigh[i]->bond[0]->y*vtx->neigh[i]->bond[1]->y + vtx->neigh[i]->bond[0]->z*vtx->neigh[i]->bond[1]->z)/vtx->neigh[i]->bond[0]->bond_length/vtx->neigh[i]->bond[1]->bond_length;
            delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
        }
    }
    // poly->k is filament persistence length (in units l_min)
    delta_energy *= 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->neigh_no;i++){
        memcpy(vtx->neigh[i],&backupneigh[i],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
@@ -2,3 +2,5 @@
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);
ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn);