Merge branch 'trout-rbc'

Samo Penic

2014-04-28 b36841fb1903e9945cb4a1dbc4dd7761556e07e2

Merge branch 'trout-rbc'

 aclocal.m4

@@ -1,4 +1,4 @@
# generated automatically by aclocal 1.13.3 -*- Autoconf -*-
# generated automatically by aclocal 1.14.1 -*- Autoconf -*-

# Copyright (C) 1996-2013 Free Software Foundation, Inc.

@@ -32,10 +32,10 @@
# generated from the m4 files accompanying Automake X.Y.
# (This private macro should not be called outside this file.)
AC_DEFUN([AM_AUTOMAKE_VERSION],
[am__api_version='1.13'
[am__api_version='1.14'
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.13.3], [],
m4_if([$1], [1.14.1], [],
      [AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])

@@ -51,7 +51,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.13.3])dnl
[AM_AUTOMAKE_VERSION([1.14.1])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
@@ -418,6 +418,12 @@
# This macro actually does too much.  Some checks are only needed if
# your package does certain things.  But this isn't really a big deal.

dnl Redefine AC_PROG_CC to automatically invoke _AM_PROG_CC_C_O.
m4_define([AC_PROG_CC],
m4_defn([AC_PROG_CC])
[_AM_PROG_CC_C_O
])

# AM_INIT_AUTOMAKE(PACKAGE, VERSION, [NO-DEFINE])
# AM_INIT_AUTOMAKE([OPTIONS])
# -----------------------------------------------
@@ -526,14 +532,54 @@
AC_CONFIG_COMMANDS_PRE(dnl
[m4_provide_if([_AM_COMPILER_EXEEXT],
  [AM_CONDITIONAL([am__EXEEXT], [test -n "$EXEEXT"])])])dnl
])

# POSIX will say in a future version that running "rm -f" with no argument
# is OK; and we want to be able to make that assumption in our Makefile
# recipes.  So use an aggressive probe to check that the usage we want is
# actually supported "in the wild" to an acceptable degree.
# See automake bug#10828.
# To make any issue more visible, cause the running configure to be aborted
# by default if the 'rm' program in use doesn't match our expectations; the
# user can still override this though.
if rm -f && rm -fr && rm -rf; then : OK; else
  cat >&2 <<'END'
Oops!

Your 'rm' program seems unable to run without file operands specified
on the command line, even when the '-f' option is present.  This is contrary
to the behaviour of most rm programs out there, and not conforming with
the upcoming POSIX standard: <http://austingroupbugs.net/view.php?id=542>

Please tell bug-automake@gnu.org about your system, including the value
of your $PATH and any error possibly output before this message.  This
can help us improve future automake versions.

END
  if test x"$ACCEPT_INFERIOR_RM_PROGRAM" = x"yes"; then
    echo 'Configuration will proceed anyway, since you have set the' >&2
    echo 'ACCEPT_INFERIOR_RM_PROGRAM variable to "yes"' >&2
    echo >&2
  else
    cat >&2 <<'END'
Aborting the configuration process, to ensure you take notice of the issue.

You can download and install GNU coreutils to get an 'rm' implementation
that behaves properly: <http://www.gnu.org/software/coreutils/>.

If you want to complete the configuration process using your problematic
'rm' anyway, export the environment variable ACCEPT_INFERIOR_RM_PROGRAM
to "yes", and re-run configure.

END
    AC_MSG_ERROR([Your 'rm' program is bad, sorry.])
  fi
fi])

dnl Hook into '_AC_COMPILER_EXEEXT' early to learn its expansion.  Do not
dnl add the conditional right here, as _AC_COMPILER_EXEEXT may be further
dnl mangled by Autoconf and run in a shell conditional statement.
m4_define([_AC_COMPILER_EXEEXT],
m4_defn([_AC_COMPILER_EXEEXT])[m4_provide([_AM_COMPILER_EXEEXT])])


# When config.status generates a header, we must update the stamp-h file.
# This file resides in the same directory as the config header
@@ -716,6 +762,70 @@
AC_DEFUN([_AM_IF_OPTION],
[m4_ifset(_AM_MANGLE_OPTION([$1]), [$2], [$3])])

# Copyright (C) 1999-2013 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# _AM_PROG_CC_C_O
# ---------------
# Like AC_PROG_CC_C_O, but changed for automake.  We rewrite AC_PROG_CC
# to automatically call this.
AC_DEFUN([_AM_PROG_CC_C_O],
[AC_REQUIRE([AM_AUX_DIR_EXPAND])dnl
AC_REQUIRE_AUX_FILE([compile])dnl
AC_LANG_PUSH([C])dnl
AC_CACHE_CHECK(
  [whether $CC understands -c and -o together],
  [am_cv_prog_cc_c_o],
  [AC_LANG_CONFTEST([AC_LANG_PROGRAM([])])
  # Make sure it works both with $CC and with simple cc.
  # Following AC_PROG_CC_C_O, we do the test twice because some
  # compilers refuse to overwrite an existing .o file with -o,
  # though they will create one.
  am_cv_prog_cc_c_o=yes
  for am_i in 1 2; do
    if AM_RUN_LOG([$CC -c conftest.$ac_ext -o conftest2.$ac_objext]) \
         && test -f conftest2.$ac_objext; then
      : OK
    else
      am_cv_prog_cc_c_o=no
      break
    fi
  done
  rm -f core conftest*
  unset am_i])
if test "$am_cv_prog_cc_c_o" != yes; then
   # Losing compiler, so override with the script.
   # FIXME: It is wrong to rewrite CC.
   # But if we don't then we get into trouble of one sort or another.
   # A longer-term fix would be to have automake use am__CC in this case,
   # and then we could set am__CC="\$(top_srcdir)/compile \$(CC)"
   CC="$am_aux_dir/compile $CC"
fi
AC_LANG_POP([C])])

# For backward compatibility.
AC_DEFUN_ONCE([AM_PROG_CC_C_O], [AC_REQUIRE([AC_PROG_CC])])

# Copyright (C) 2001-2013 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# AM_RUN_LOG(COMMAND)
# -------------------
# Run COMMAND, save the exit status in ac_status, and log it.
# (This has been adapted from Autoconf's _AC_RUN_LOG macro.)
AC_DEFUN([AM_RUN_LOG],
[{ echo "$as_me:$LINENO: $1" >&AS_MESSAGE_LOG_FD
   ($1) >&AS_MESSAGE_LOG_FD 2>&AS_MESSAGE_LOG_FD
   ac_status=$?
   echo "$as_me:$LINENO: \$? = $ac_status" >&AS_MESSAGE_LOG_FD
   (exit $ac_status); }])

# Check to make sure that the build environment is sane.    -*- Autoconf -*-

# Copyright (C) 1996-2013 Free Software Foundation, Inc.

 config.h

@@ -7,6 +7,12 @@
/* Define to 1 if you have the `confuse' library (-lconfuse). */
#define HAVE_LIBCONFUSE 1

/* Define to 1 if you have the `gsl' library (-lgsl). */
#define HAVE_LIBGSL 1

/* Define to 1 if you have the `gslcblas' library (-lgslcblas). */
#define HAVE_LIBGSLCBLAS 1

/* Define to 1 if you have the `m' library (-lm). */
#define HAVE_LIBM 1


 config.h.in

@@ -6,6 +6,12 @@
/* Define to 1 if you have the `confuse' library (-lconfuse). */
#undef HAVE_LIBCONFUSE

/* Define to 1 if you have the `gsl' library (-lgsl). */
#undef HAVE_LIBGSL

/* Define to 1 if you have the `gslcblas' library (-lgslcblas). */
#undef HAVE_LIBGSLCBLAS

/* Define to 1 if you have the `m' library (-lm). */
#undef HAVE_LIBM


 configure.ac

@@ -15,6 +15,9 @@
# FIXME: Replace `main' with a function in `-lm':
AC_CHECK_LIB([m], [pow])

AC_CHECK_LIB([gslcblas],[cblas_dgemm])
AC_CHECK_LIB([gsl],[gsl_blas_dgemm])

# Checks for header files.
AC_CHECK_HEADERS([stdlib.h string.h])


 src/Makefile.am

@@ -1,17 +1,17 @@
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 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 poly.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 stats.c sh.c shcomplex.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 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 poly.c
spherical_trisurf_ffdir=../
spherical_trisurf_ff_PROGRAMS = spherical_trisurf_ff
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
#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 poly.c stats.c shcomplex.c
#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 poly.c stats.c shcomplex.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 poly.c stats.c shcomplex.c
#spherical_trisurf_ffdir=../
#spherical_trisurf_ff_PROGRAMS = spherical_trisurf_ff
#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 stats.c shcomplex.c

 src/cell.c

@@ -44,7 +44,6 @@
}


//TODO: not debugged at all!
inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx){
    ts_uint cellidx;
    ts_uint ncx, ncy,ncz;
@@ -68,7 +67,6 @@
}


//TODO: looks ok, but debug anyway in the future
inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){
    ts_uint i;
    for(i=0;i<cell->nvertex;i++){
@@ -124,7 +122,7 @@
    return TS_SUCCESS;
}

// 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_uint ncx,ncy,ncz,remainder,cell_occupation;
    ts_uint i,j,k,l,neigh_cidx;
@@ -151,11 +149,16 @@
                    for(l=0;l<cell_occupation;l++){

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

//                if(vtx->idx==1)
//                fprintf(stderr,"VTX(0) ima bliznji vertex z indeksom, %d, tipa %d \n", clist->cell[neigh_cidx]->vertex[l]->idx, clist->cell[neigh_cidx]->vertex[l]->id);
//                if(vtx->idx==0 && clist->cell[neigh_cidx]->vertex[l]->idx==0)
//                            fprintf(stderr,"*** dist was %f\n",dist);
				
                            if(dist<=1.0 || (dist<=clist->dmin_interspecies && (clist->cell[neigh_cidx]->vertex[l]->id != vtx->id))) return TS_FAIL;
                        }
                    }
                }

 src/frame.c

@@ -2,6 +2,7 @@
#include "general.h"
#include "cell.h"
#include "frame.h"

ts_bool centermass(ts_vesicle *vesicle){
    ts_uint i,j, n=vesicle->vlist->n;
    ts_vertex **vtx=vesicle->vlist->vtx;

 src/general.h

@@ -3,6 +3,7 @@

#include<stdarg.h>
#include<stdio.h>
#include<gsl/gsl_complex.h>
/* @brief This is a header file, defining general constants and structures.
  * @file header.h
  * @author Samo Penic
@@ -47,6 +48,8 @@
#define TS_FAIL 1

/* CONSTANTS */

#define TS_ID_FILAMENT 1

/* DATA TYPES */
/** @brief Sets the default datatype for ts_double
@@ -213,12 +216,14 @@
    ts_double dcell;
    ts_double shift;
    ts_double max_occupancy;
    ts_double dmin_interspecies;
} ts_cell_list;


typedef struct {
    ts_uint l;
    ts_double **ulm;
    gsl_complex **ulmComplex;
    ts_double **sumUlm2;
    ts_uint N;
    ts_double **co;

 src/initial_distribution.c

@@ -11,6 +11,8 @@
#include "energy.h"
#include "poly.h"
#include "io.h"
#include "sh.h"
#include "shcomplex.h"

ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
    ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
@@ -41,7 +43,9 @@

    vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
    // Nucleus:
    vesicle->R_nucleus=tape->R_nucleus;
    vesicle->R_nucleus=tape->R_nucleus*tape->R_nucleus;

    vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;

    //Initialize grafted polymers (brush):
    vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle);
@@ -69,6 +73,9 @@
        }
    }

    for(i=0;i<vesicle->filament_list->n;i++){
        vertex_list_assign_id(vesicle->filament_list->poly[i]->vlist,TS_ID_FILAMENT);
    }

//    vesicle->spring_constant=tape->kspring;
//    poly_assign_spring_const(vesicle);
@@ -87,7 +94,12 @@

    vesicle->pressure= tape->pressure;
    vesicle->pswitch=tape->pswitch;

    if(tape->shc>0){
        vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
    }
    else {
        vesicle->sphHarmonics=NULL;
    }
    return vesicle;

}

 src/io.c

@@ -918,6 +918,39 @@



ts_bool write_pov_file(ts_vesicle *vesicle, char *filename){
    FILE *fh;
    ts_uint i;
	
    fh=fopen(filename, "w");
    if(fh==NULL){
        err("Cannot open file %s for writing");
        return TS_FAIL;
    }

    for(i=0;i<vesicle->tlist->n;i++){
	
    fprintf(fh,"\ttriangle {");
    fprintf(fh,"\t<%e,%e,%e> <%e,%e,%e> <%e,%e,%e> }\n", 
    vesicle->tlist->tria[i]->vertex[0]->x,
    vesicle->tlist->tria[i]->vertex[0]->y,
    vesicle->tlist->tria[i]->vertex[0]->z,

    vesicle->tlist->tria[i]->vertex[1]->x,
    vesicle->tlist->tria[i]->vertex[1]->y,
    vesicle->tlist->tria[i]->vertex[1]->z,

    vesicle->tlist->tria[i]->vertex[2]->x,
    vesicle->tlist->tria[i]->vertex[2]->y,
    vesicle->tlist->tria[i]->vertex[2]->z
    );
    }
		
    fclose(fh);
    return TS_SUCCESS;
}


ts_tape *parsetape(char *filename){
  //  long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20, pswitch=0;  // THIS IS DUE TO CONFUSE BUG!
    ts_tape *tape=(ts_tape *)calloc(1,sizeof(ts_tape));
@@ -937,6 +970,7 @@
    CFG_SIMPLE_INT("nfono",&tape->nfono),
    CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
    CFG_SIMPLE_FLOAT("dmax", &tape->dmax),
    CFG_SIMPLE_FLOAT("dmin_interspecies", &tape->dmin_interspecies),
        CFG_SIMPLE_FLOAT("xk0",&tape->xk0),
    CFG_SIMPLE_INT("pswitch",&tape->pswitch),
    CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
@@ -954,6 +988,7 @@
        CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
        CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
        CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
    CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
        CFG_END()
    };
    cfg_t *cfg;    

 src/io.h

@@ -26,6 +26,7 @@
        long int brezveze2;
        ts_double xk0;
    ts_double dmax;
    ts_double dmin_interspecies;
    ts_double stepsize;
    ts_double kspring;
    ts_double xi;
@@ -34,6 +35,7 @@
    long int inititer;
    long int mcsweeps;
    long int quiet;
    long int shc;
} ts_tape;

typedef struct{
@@ -98,6 +100,8 @@
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 write_pov_file(ts_vesicle *vesicle, char *filename);

ts_tape *parsetape(char *filename);
ts_bool tape_free(ts_tape *tape);


 src/main.c

@@ -11,6 +11,8 @@
#include "frame.h"
#include "timestep.h"
#include "poly.h"
#include "sh.h"
#include "shcomplex.h"

/** Entrance function to the program
  * @param argv is a number of parameters used in program call (including the program name
@@ -45,11 +47,20 @@
        // nove vrednosti iz tapea...
        vesicle->bending_rigidity=tape->xk0;
        vtx_set_global_values(vesicle);
        vesicle->pswitch =tape->pswitch;
        vesicle->pressure=tape->pressure;
        vesicle->dmax=tape->dmax*tape->dmax;
        poly_assign_filament_xi(vesicle,tape);
        vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
        /* spherical harmonics */
        if(tape->shc>0){
            vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
        }
        else {
            vesicle->sphHarmonics=NULL;
        }

        if(command_line_args.reset_iteration_count) start_iteration=tape->inititer-1;
        if(command_line_args.reset_iteration_count) start_iteration=tape->inititer;
        else start_iteration++;

        if(start_iteration>=tape->iterations){

 src/poly.c

@@ -58,7 +58,6 @@
    ts_uint gvtxi;
    ts_double xnorm,ynorm,znorm,normlength;
    ts_double dphi,dh;
    ts_uint ji;

    // Grafting polymers:
    if (vlist!=NULL){
@@ -108,14 +107,29 @@
    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;
    ts_double a,R,H,tantheta,h,r,phi,A0=1.2;

        a = A0*(ts_double)vesicle->nshell;
        R = A0*((ts_double)vesicle->nshell)/(2.0*sin(M_PI/5.0));
        H = sqrt(a*a - R*R);
        tantheta = sqrt(R*R - a*a/4.0)/H;
		
        h = -H + sqrt(vesicle->clist->dmin_interspecies)*1.5;
        r = (H-fabs(h))*tantheta - sqrt(vesicle->clist->dmin_interspecies)*1.5;
        dphi = 2.0*asin(1.0/2.0/r)*1.001;
        dh = dphi/2.0/M_PI*1.001;
        phi=0.0;
        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);
                h = h + dh;
                r = (H-fabs(h))*tantheta - sqrt(vesicle->clist->dmin_interspecies)*1.5;
                dphi = 2.0*asin(1.0/2.0/r)*1.001;
                dh = dphi/2.0/M_PI*1.001;
                phi+=dphi;
                //ji = j + i*poly_list->poly[i]->vlist->n;
                poly_list->poly[i]->vlist->vtx[j]->x = r*cos(phi);
                poly_list->poly[i]->vlist->vtx[j]->y = r*sin(phi);
                poly_list->poly[i]->vlist->vtx[j]->z = h;// ji*dh - (dh*poly_list->n*poly_list->poly[i]->vlist->n/2.0);
            }
        }
    }

 src/sh.c

@@ -49,6 +49,7 @@

ts_bool sph_free(ts_spharm *sph){
    int i,j;
    if(sph==NULL) return TS_FAIL;
    for(i=0;i<sph->l;i++){
        if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
        if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
@@ -189,7 +190,7 @@
        K=-sqrt(1.0/(M_PI))*cos(m*fi);
    }
    
    return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));	
    return K*sqrt((2.0*l+1.0)/2.0*(ts_double)(fac1/fac2))*plgndr(l,abs(m),cos(theta));	
}


@@ -200,7 +201,7 @@
#ifdef TS_DOUBLE_DOUBLE
    coord->e1=sqrt(x*x+y*y+z*z);
    if(z==0) coord->e3=M_PI/2.0;
    else coord->e3=atan(sqrt(x*x+y*y)/z);
    else coord->e3=atan2(sqrt(x*x+y*y),z);
    coord->e2=atan2(y,x);
#endif
#ifdef TS_DOUBLE_FLOAT
@@ -218,6 +219,23 @@

    return TS_SUCCESS;
}


ts_bool sph2cart(ts_coord *coord){
    coord->coord_type=TS_COORD_CARTESIAN;
    ts_double x,y,z;

    x=coord->e1*cos(coord->e2)*sin(coord->e3);
    y=coord->e1*sin(coord->e2)*sin(coord->e3);
    z=coord->e1*cos(coord->e3);

    coord->e1=x;
    coord->e2=y;
    coord->e3=z;

    return TS_SUCCESS;
}


/* Function returns radius of the sphere with the same volume as vesicle (r0) */
ts_double getR0(ts_vesicle *vesicle){
@@ -379,3 +397,28 @@
    sph->N++;
return TS_SUCCESS;
}


ts_bool saveAvgUlm2(ts_vesicle *vesicle){

    FILE *fh;
	
    fh=fopen("sph2out.dat", "w");
    if(fh==NULL){
        err("Cannot open file %s for writing");
        return TS_FAIL;
    }

    ts_spharm *sph=vesicle->sphHarmonics;
    ts_int i,j;
    fprintf(fh,"l,\tm,\tulm^2avg\n");
    for(i=0;i<sph->l;i++){
            for(j=0;j<2*i+1;j++){
        fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);

            }
    fprintf(fh,"\n");
    }
    fclose(fh);
    return TS_SUCCESS;
}

 src/sh.h

@@ -6,6 +6,7 @@
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi);

ts_bool *cart2sph(ts_coord *coord, ts_double x, ts_double y, ts_double z);
ts_bool sph2cart(ts_coord *coord);
ts_bool precomputeShCoeff(ts_spharm *sph);
ts_spharm *sph_init(ts_vertex_list *vlist, ts_uint l);
ts_bool sph_free(ts_spharm *sph);
@@ -13,4 +14,5 @@
ts_bool preparationSh(ts_vesicle *vesicle, ts_double r0);
ts_bool calculateYlmi(ts_vesicle *vesicle);
ts_bool calculateUlm(ts_vesicle *vesicle);
ts_bool saveAvgUlm2(ts_vesicle *vesicle);
#endif

 src/shcomplex.c

New file
@@ -0,0 +1,172 @@
#include<math.h>
#include<stdlib.h>
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_sf_legendre.h>

#include<gsl/gsl_matrix.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_linalg.h>
#include "general.h"
#include "sh.h"
#include "shcomplex.h"


ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
    ts_uint j,i;
    ts_spharm *sph=(ts_spharm *)malloc(sizeof(ts_spharm));

    sph->N=0;
    /* lets initialize Ylm for each vertex. */
    sph->Ylmi=(ts_double ***)calloc(l,sizeof(ts_double **));
    for(i=0;i<l;i++){
            sph->Ylmi[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
            for(j=0;j<(2*i+1);j++){
                sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
            }
    }
        
    /* lets initialize ulm */
    sph->ulm=(ts_double **)calloc(l,sizeof(ts_double *));
    sph->ulmComplex=(gsl_complex **)calloc(l,sizeof(gsl_complex *));
    for(j=0;j<l;j++){
        sph->ulm[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
        sph->ulmComplex[j]=(gsl_complex *)calloc(2*j+1,sizeof(gsl_complex));
    }

    /* lets initialize sum of Ulm2 */
    sph->sumUlm2=(ts_double **)calloc(l,sizeof(ts_double *));
    for(j=0;j<l;j++){
        sph->sumUlm2[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
    }

    /* lets initialize co */
//NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2*j+2 instead of 2*j+2. elements starting with 0 are useles and should be ignored!
    sph->co=(ts_double **)calloc(l+1,sizeof(ts_double *));
    for(j=0;j<=l;j++){
        sph->co[j]=(ts_double *)calloc(2*j+2,sizeof(ts_double));
    }

    sph->l=l;   

    /* Calculate coefficients that will remain constant during all the simulation */ 
   precomputeShCoeff(sph);
    
    return sph;
}

ts_bool complex_sph_free(ts_spharm *sph){
    int i,j;
    if(sph==NULL) return TS_FAIL;
    for(i=0;i<sph->l;i++){
        if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
        if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
        if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
        if(sph->co[i]!=NULL) free(sph->co[i]);
    }
        if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
    if(sph->co != NULL) free(sph->co);
    if(sph->ulm !=NULL) free(sph->ulm);
    if(sph->ulmComplex !=NULL) free(sph->ulmComplex);

        if(sph->Ylmi!=NULL) {
            for(i=0;i<sph->l;i++){
                if(sph->Ylmi[i]!=NULL){
                    for(j=0;j<i*2+1;j++){
                        if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
                    }
                    free(sph->Ylmi[i]);
                }
            }
            free(sph->Ylmi);
        }

    free(sph);
    return TS_SUCCESS;
}


ts_bool calculateUlmComplex(ts_vesicle *vesicle){
    ts_int i,j,k,m,l;
    ts_vertex *cvtx;
    ts_coord coord;
/* set all values to zero */
    for(i=0;i<vesicle->sphHarmonics->l;i++){
        for(j=0;j<2*i+1;j++) GSL_SET_COMPLEX(&(vesicle->sphHarmonics->ulmComplex[i][j]),0.0,0.0);
    }

    for(k=0;k<vesicle->vlist->n; k++){
        cvtx=vesicle->vlist->vtx[k];
    cart2sph(&coord,cvtx->x,cvtx->y,cvtx->z);
        for(i=0;i<vesicle->sphHarmonics->l;i++){
            for(j=0;j<2*i+1;j++){
        m=j-i;
        l=i;
        if(m>=0){	
    //    fprintf(stderr, "Racunam za l=%d, m=%d\n", l,m);
                vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*gsl_sf_legendre_sphPlm(l,m,cos(coord.e3)))) );
        } else {
    //    fprintf(stderr, "Racunam za l=%d, abs(m=%d)\n", l,m);
                vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*pow(-1,m)*gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))) );

        }
            }
        }
    }
    return TS_SUCCESS;
}

ts_bool storeUlmComplex2(ts_vesicle *vesicle){

    ts_spharm *sph=vesicle->sphHarmonics;
    ts_int i,j;
    for(i=0;i<sph->l;i++){
            for(j=0;j<2*i+1;j++){
                sph->sumUlm2[i][j]+=gsl_complex_abs2(sph->ulmComplex[i][j]);
            }
    }
    sph->N++;
    return TS_SUCCESS;
}


ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax){
    ts_int min=lmin;
    ts_int max=lmax; //vesicle->sphHarmonics->l-3;
    ts_long i,j;
    ts_double retval, bval;
    gsl_matrix *A=gsl_matrix_alloc(max-min,2);
    gsl_vector *tau=gsl_vector_alloc(2);
    gsl_vector *b=gsl_vector_alloc(max-min);
    gsl_vector *x=gsl_vector_alloc(2);
    gsl_vector *res=gsl_vector_alloc(max-min);

    //solving (A^T*A)*x=A^T*b
    //fill the data for matrix A and vector b
    for(i=min;i<max;i++){
            gsl_matrix_set(A, i-min,0,(ts_double)((i-1)*(i+2)));
            gsl_matrix_set(A, i-min,1,(ts_double)((i-1)*(i+2)*(i+1)*i));
//            fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
            bval=0.0;
            //average for m from 0..l (only positive m's)
            for(j=0;j<=i;j++){
                bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)];
            }
                bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1);

            gsl_vector_set(b,i-min,1.0/bval);
//            fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
    }
//    fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
    gsl_linalg_QR_decomp(A,tau);
    gsl_linalg_QR_lssolve(A,tau,b,x,res);
//    fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
    retval=gsl_vector_get(x,1);
    gsl_matrix_free(A);
    gsl_vector_free(tau);
    gsl_vector_free(b);
    gsl_vector_free(x);
    gsl_vector_free(res);
    
    return retval;
}

 src/shcomplex.h

New file
@@ -0,0 +1,9 @@
#ifndef _H_SH_COMPLEX
#define _H_SH_COMPLEX
#include "general.h"
ts_bool storeUlmComplex2(ts_vesicle *vesicle);
ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l);
ts_bool complex_sph_free(ts_spharm *sph);
ts_bool calculateUlmComplex(ts_vesicle *vesicle);
ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax);
#endif

 src/spherical_trisurf.c

@@ -61,10 +61,11 @@
calculateUlm(vesicle);

//preloop:
ts_double vmsr, bfsr;
for(i=0;i<1000;i++){
    cell_occupation(vesicle);
    for(j=0;j<1000;j++){
        single_timestep(vesicle);
        single_timestep(vesicle, &vmsr, &bfsr);
    }    
    centermass(vesicle);
    fprintf(stderr, "Preloop %d completed.\n",i+1);
@@ -75,7 +76,7 @@
    for(j=0;j<200;j++){
        cell_occupation(vesicle);
        for(k=0;k<5;k++){
        single_timestep(vesicle);
        single_timestep(vesicle, &vmsr, &bfsr);
        }
        centermass(vesicle);
    }    
@@ -103,26 +104,3 @@



ts_bool saveAvgUlm2(ts_vesicle *vesicle){

    FILE *fh;
	
    fh=fopen("sph2out.dat", "w");
    if(fh==NULL){
        err("Cannot open file %s for writing");
        return TS_FAIL;
    }

    ts_spharm *sph=vesicle->sphHarmonics;
    ts_int i,j;
    fprintf(fh,"l,\tm,\tulm^2avg\n");
    for(i=0;i<sph->l;i++){
            for(j=0;j<2*i+1;j++){
        fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);

            }
    fprintf(fh,"\n");
    }
    fclose(fh);
    return TS_SUCCESS;
}

 src/spherical_trisurf_ff.c

@@ -44,10 +44,11 @@
preparationSh(vesicle,r0);
calculateYlmi(vesicle);
calculateUlm(vesicle);
ts_double vmsr,bfsr;
for(i=0;i<500;i++){
    cell_occupation(vesicle);
    for(j=0;j<1000;j++){
        single_timestep(vesicle);
        single_timestep(vesicle,&vmsr,&bfsr);
    }    
    centermass(vesicle);
    fprintf(stderr, "Preloop %d completed.\n",i+1);
@@ -62,7 +63,7 @@
for(i=0;i<10000;i++){
    cell_occupation(vesicle);
    for(j=0;j<1000;j++){
        single_timestep(vesicle);
        single_timestep(vesicle,&vmsr,&bfsr);
    }    
    centermass(vesicle);
    vesicle_volume(vesicle);
@@ -86,27 +87,3 @@
}



ts_bool saveAvgUlm2(ts_vesicle *vesicle){

    FILE *fh;
	
    fh=fopen("sph2out.dat", "w");
    if(fh==NULL){
        err("Cannot open file %s for writing");
        return TS_FAIL;
    }

    ts_spharm *sph=vesicle->sphHarmonics;
    ts_int i,j;
    fprintf(fh,"l,\tm,\tulm^2avg\n");
    for(i=0;i<sph->l;i++){
            for(j=0;j<2*i+1;j++){
        fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);

            }
    fprintf(fh,"\n");
    }
    fclose(fh);
    return TS_SUCCESS;
}

 src/stats.c

New file
@@ -0,0 +1,83 @@
#include <stdio.h>
#include <time.h>
//#include <gsl/gsl_math.h>
//#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sort_vector.h>
#include "general.h"
 
void gyration_eigen(ts_vesicle *vesicle, ts_double *l1, ts_double *l2, ts_double *l3){
    gsl_matrix *matrix = gsl_matrix_alloc(3,3);
    gsl_vector *eig = gsl_vector_alloc(3);
    gsl_eigen_symm_workspace *workspace = gsl_eigen_symm_alloc(3);
    ts_uint i,j,k;
    ts_double mat[3][3];
    ts_double vec[3];

    for(i = 0; i < 3; i++)
            for(j = 0; j < 3; j++)
                    mat[i][j]=0;
 
    for(k=0;k<vesicle->vlist->n;k++){
        vec[0]=vesicle->vlist->vtx[k]->x;
        vec[1]=vesicle->vlist->vtx[k]->y;
        vec[2]=vesicle->vlist->vtx[k]->z;
        for(i = 0; i < 3; i++)
                for(j = 0; j <= i; j++)
                mat[i][j]+=vec[i]*vec[j];
    }

// Normalize gyration tensor:
    for(i = 0; i < 3; i++)
            for(j = 0; j <= i; j++)
            mat[i][j]=mat[i][j]/(ts_double)vesicle->vlist->n;


// diagonal elements are copied twice!	
    for(i = 0; i < 3; i++)
            for(j = 0; j <= i; j++){
                    gsl_matrix_set(matrix,i,j,mat[i][j]);
                    gsl_matrix_set(matrix,j,i,mat[i][j]);
    }
/*
for(i = 0; i < 3; i++)
    {
        for(j = 0; j < 3; j++)
            printf("%g ", gsl_matrix_get(matrix, i, j));
        printf("\n");
    }
    printf("\n");
*/
    gsl_eigen_symm(matrix, eig, workspace);
    gsl_sort_vector(eig);
/* printf("** eigenvalues    \n");
    for(i = 0; i < 3; i++){
        printf("%3d %25.17e\n", i, gsl_vector_get(eig, i));
} */
    *l1=gsl_vector_get(eig,0);
    *l2=gsl_vector_get(eig,1);
    *l3=gsl_vector_get(eig,2);

    gsl_eigen_symm_free(workspace);
    gsl_matrix_free(matrix);
    gsl_vector_free(eig);
}

ts_ulong get_epoch(){
    time_t currTime;
    currTime = time(NULL);
    return (ts_ulong)currTime;
}

ts_bool get_area_volume(ts_vesicle *vesicle, ts_double *area, ts_double *volume){
    ts_uint i;
    *volume=0.0;
    *area=0.0;
    for(i=0;i<vesicle->tlist->n;i++){
        *volume+=vesicle->tlist->tria[i]->volume;
        *area+=vesicle->tlist->tria[i]->area;
    }
    return TS_SUCCESS;
} 

 src/stats.h

New file
@@ -0,0 +1,6 @@
#ifndef _H_STATS
#define _H_STATS
void gyration_eigen(ts_vesicle *vesicle, ts_double *l1, ts_double *l2, ts_double *l3);
ts_ulong get_epoch();
ts_bool get_area_volume(ts_vesicle *vesicle, ts_double *area, ts_double *volume);
#endif

 src/tape

@@ -3,6 +3,8 @@
nshell=17
# dmax is the max. bond length (in units l_min)
dmax=1.7
# dmin_interspecies in the min. dist. between different vertex species (in units l_min)
dmin_interspecies=1.2
# bending rigidity of the membrane (in units kT)
xk0=10.0
# max step size (in units l_min)
@@ -12,11 +14,11 @@
# (pswitch=1: calc. p*dV energy contribution)
pswitch = 0
# pressure difference: p_inside - p_outside (in units kT/l_min^3):
pressure=0.1
pressure=0.0

####### Polymer (brush) definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
npoly=300
npoly=0
# nmono is a number of monomers in each polymer
nmono=10
# Spring constant between monomers of the polymer
@@ -24,15 +26,15 @@

####### Filament (inside the vesicle) definitions ###########
# nfil is a number of filaments inside the vesicle
nfil=1
nfil=0
# nfono is a number of monomers in each filament
nfono=300
# Persistence lenght of the filaments (in units l_min)
xi=100
xi=0

####### Nucleus (inside the vesicle) ###########
# Radius of an impenetrable hard sphere inside the vesicle
R_nucleus=5
R_nucleus=0

#######  Cell definitions ############
nxmax=60
@@ -42,13 +44,16 @@

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


###### Spherical harmonics ###########
spherical_harmonics_coefficients=21

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


 src/timestep.c

@@ -8,51 +8,103 @@
#include "bondflip.h"
#include "frame.h"
#include "io.h"
#include "stats.h"
#include "sh.h"
#include "shcomplex.h"
#include "vesicle.h"

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

    ts_uint i, j,k;
    ts_double r0,kc1,kc2,kc3,kc4;
    ts_double l1,l2,l3,volume=0.0,area=0.0,vmsr,bfsr, vmsrt, bfsrt;
    ts_ulong epochtime;
    FILE *fd1;
//     char filename[255];
    FILE *fd=fopen("statistics.csv","w");
    if(fd==NULL){
        fatal("Cannot open statistics.csv file for writing",1);
    }
    fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lambda3 Kc(2-9) Kc(6-9) Kc(2-end) Kc(3-6)\n");
    centermass(vesicle);
    cell_occupation(vesicle);
    if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
    for(i=start_iteration;i<inititer+iterations;i++){
        vmsr=0.0;
        bfsr=0.0;
        for(j=0;j<mcsweeps;j++){
            single_timestep(vesicle);
            single_timestep(vesicle, &vmsrt, &bfsrt);
            vmsr+=vmsrt;
            bfsr+=bfsrt;
        }
        vmsr/=(ts_double)mcsweeps;
        bfsr/=(ts_double)mcsweeps;
        centermass(vesicle);
        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){
            write_vertex_xml_file(vesicle,i-inititer);
            write_master_xml_file("test.pvd");
            epochtime=get_epoch();			
            gyration_eigen(vesicle, &l1, &l2, &l3);
            vesicle_volume(vesicle); //calculates just volume. Area is not added to ts_vesicle yet!
            get_area_volume(vesicle, &area,&volume); //that's why I must recalculate area (and volume for no particular reason).
            r0=getR0(vesicle);
            if(vesicle->sphHarmonics!=NULL){
                preparationSh(vesicle,r0);
                //calculateYlmi(vesicle);
                calculateUlmComplex(vesicle);
                storeUlmComplex2(vesicle);
                saveAvgUlm2(vesicle);
                kc1=calculateKc(vesicle, 2,9);
                kc2=calculateKc(vesicle, 6,9);
                kc3=calculateKc(vesicle, 2,vesicle->sphHarmonics->l);
                kc4=calculateKc(vesicle, 3,6);
            
                fd1=fopen("state.dat","w");
                fprintf(fd1,"%e %e\n",vesicle->volume, getR0(vesicle));
                for(k=0;k<vesicle->vlist->n;k++){
                    fprintf(fd1,"%e %e %e %e %e\n",
                        vesicle->vlist->vtx[k]->x,
                        vesicle->vlist->vtx[k]->y,
                        vesicle->vlist->vtx[k]->z,
                        vesicle->vlist->vtx[k]->solAngle,
                        vesicle->vlist->vtx[k]->relR
                    );
                }
                fclose(fd1);
            }

            fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3,kc1, kc2, kc3,kc4);
            fflush(fd);	
        //    sprintf(filename,"timestep-%05d.pov",i-inititer);
        //    write_pov_file(vesicle,filename);
        }
    }
    fclose(fd);
    return TS_SUCCESS;
}

ts_bool single_timestep(ts_vesicle *vesicle){
ts_bool single_timestep(ts_vesicle *vesicle,ts_double *vmsr, ts_double *bfsr){
    ts_bool retval;
    ts_double rnvec[3];
    ts_uint i,j,b;
    ts_uint vmsrcnt=0;
    for(i=0;i<vesicle->vlist->n;i++){
        rnvec[0]=drand48();
        rnvec[1]=drand48();
        rnvec[2]=drand48();
        retval=single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec);
    if(retval==TS_SUCCESS) vmsrcnt++;        
    }

//    ts_int cnt=0;
    ts_int bfsrcnt=0;
    for(i=0;i<3*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[b],rnvec);
//    if(retval==TS_SUCCESS) cnt++;        
    if(retval==TS_SUCCESS) bfsrcnt++;        
    }

    for(i=0;i<vesicle->poly_list->n;i++){
@@ -76,7 +128,8 @@
 

//    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);
    *vmsr=(ts_double)vmsrcnt/(ts_double)vesicle->vlist->n;
    *bfsr=(ts_double)bfsrcnt/(ts_double)vesicle->vlist->n/3.0;
    return TS_SUCCESS;
}


 src/timestep.h

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

 src/triangle.c

@@ -5,6 +5,7 @@
#include<math.h>

/** @brief Prepares the list for triangles.
  * @returns pointer to empty data structure for maintaining triangle list.
  *
  * Create empty list for holding the information on triangles. Triangles are
  * added later on with triangle_add().
@@ -27,14 +28,17 @@
}

/** @brief Add the triangle to the triangle list and create necessary data
 * structures.
  * structures.
  * @param *tlist is a pointer to triangle list where triangle should be created
  * @param *vtx1, *vtx2, *vtx3 are the three vertices defining the triangle
  * @returns pointer to the newly created triangle on success and NULL if
  * triangle could not be created. It breaks program execution if memory
  * allocation of triangle list can't be done.
  *
  * Add the triangle ts_triangle with ts_triangle_data to the ts_triangle_list.
  * Add the triangle ts_triangle to the ts_triangle_list.
  * The triangle list is resized, the ts_triangle is allocated and
  * ts_triangle_data is allocated and zeroed. The function receives 4 arguments:
  * ts_triangle_list *tlist as list of triangles and 3 ts_vertex *vtx as
  * vertices that are used to form a triangle. Returns a pointer to newly
  * created triangle. This pointer doesn't need assigning, since it is
  * triangle data is zeroed. Returned pointer to newly
  * created triangle doesn't need assigning, since it is
  * referenced by triangle list.
  *
  * WARNING: Function can be accelerated a bit by removing the NULL checks.
@@ -57,7 +61,6 @@

        tlist->tria[tlist->n-1]=(ts_triangle *)calloc(1,sizeof(ts_triangle));
        if(tlist->tria[tlist->n-1]==NULL) fatal("Cannot reallocate memory for additional ts_triangle.",5);
  //      tlist->tria[tlist->n-1]->data=(ts_triangle_data *)calloc(1,sizeof(ts_triangle_data));

        //NOW insert vertices!
        tlist->tria[tlist->n - 1]->idx=tlist->n-1;
@@ -68,9 +71,14 @@
}

/** @brief Add the neigbour to triangles.
  * @param *tria is a first triangle.
  * @param *ntria is a second triangle.
  * @returns TS_SUCCES on sucessful adition to the list, TS_FAIL if triangles
  * are NULL and breaks execution FATALY if memory allocation error occurs.
  *
  * Add the neigbour to the list of neighbouring triangles. The
  * neighbouring triangles are those, who share two vertices. Function resizes
  * neighbouring triangles are those, who share two vertices and corresponding
  * bond. Function resizes
  * the list and adds the pointer to neighbour. It receives two arguments of
  * ts_triangle type. It then adds second triangle to the list of first
  * triangle, but not the opposite. Upon
@@ -84,42 +92,40 @@
  * debugging stupid NULL pointers.
  *
  * Example of usage:
  *        triangle_remove_neighbour(tlist->tria[3], tlist->tria[4]);
  *        triangle_add_neighbour(tlist->tria[3], tlist->tria[4]);
  *
  *        Triangles 3 and 4 are not neighbours anymore.
  *        Triangle 4 is a neighbour of triangle 3, but (strangely) not the
  *        oposite. The function should be called again with the changed order of
  *        triangles to make neighbourship mutual.
  *        
  */

ts_bool triangle_add_neighbour(ts_triangle *tria, ts_triangle *ntria){
    if(tria==NULL || ntria==NULL) return TS_FAIL;
/*TODO: check if the neighbour already exists! Now there is no such check
 * because of the performance issue. */
    tria->neigh_no++;
    tria->neigh=realloc(tria->neigh,tria->neigh_no*sizeof(ts_triangle *));
    if(tria->neigh == NULL)
            fatal("Reallocation of memory failed during insertion of triangle neighbour in triangle_add_neighbour",3);
    tria->neigh[tria->neigh_no-1]=ntria;
  
 
/* we repeat the procedure for the neighbour */  
/*    ntria->data->neigh_no++;
    ntria->data->neigh=realloc(ntria->data->neigh,ntria->data->neigh_no*sizeof(ts_triangle *));
    if(ntria->data->neigh == NULL)
            fatal("Reallocation of memory failed during insertion of triangle neighbour in triangle_add_neighbour",3);
    ntria->data->neigh[ntria->data->neigh_no-1]=tria;
*/
    tria->neigh[tria->neigh_no-1]=ntria; 
    return TS_SUCCESS;
}

/** @brief Remove the neigbours from triangle.
  * @param *tria is a first triangle.
  * @param *ntria is neighbouring triangle.
  * @returns TS_SUCCESS on successful removal, TS_FAIL if triangles are not
  * neighbours and it breaks program execution FATALY if memory allocation
  * problem occurs.
  *
  * Removes the neigbour from the list of neighbouring triangles. The
  * neighbouring triangles are those, who share two vertices. Function resizes
  * neighbouring triangles are those, who share two vertices and corresponding
  * bond. Function resizes
  * the list and deletes the pointer to neighbour. It receives two arguments of
  * ts_triangle type. It then removes eachother form eachother's list. Upon
  * ts_triangle type. It then mutually removes triangles from eachouther
  * neighbour list. Upon
  * success it returns TS_SUCCESS, upon failure to find the triangle in the
  * neighbour list returns TS_FAIL and it FATALY ends when the datastructure
  * cannot be resized.
  * neighbour list returns TS_FAIL. It FATALY breaks program execution when the datastructure
  * cannot be resized due to memory constrain problems.
  *
  * WARNING: The function doesn't check whether the pointer is NULL or invalid. It is the
  * job of programmer to make sure the pointer is valid.
@@ -144,10 +150,8 @@
    }
    if(j==i) {
        return TS_FAIL; 
        //fatal("In triangle_remove_neighbour: Specified neighbour does not exist for given triangle",3);
    }
    tria->neigh_no--;
//    fprintf(stderr,"*** tria_number=%d\n",tria->neigh_no);
    tria->neigh=(ts_triangle **)realloc(tria->neigh,tria->neigh_no*sizeof(ts_triangle *));
    if(tria->neigh == NULL){
        fprintf(stderr,"Ooops: tria->neigh_no=%d\n",tria->neigh_no);
@@ -163,10 +167,8 @@
    }
    if(j==i) {
        return TS_FAIL; 
        //fatal("In triangle_remove_neighbour: Specified neighbour does not exist for given triangle",3);
    }
    ntria->neigh_no--;
//    fprintf(stderr,"*** ntria_number=%d\n",ntria->neigh_no);
    ntria->neigh=(ts_triangle **)realloc(ntria->neigh,ntria->neigh_no*sizeof(ts_triangle *));
    if(ntria->neigh == NULL){
        fprintf(stderr,"Ooops: ntria->neigh_no=%d\n",ntria->neigh_no);
@@ -176,25 +178,33 @@
}


/** @brief Calculates normal vector of the triangle.

/** @brief Calculates normal vector of the triangle, its corresponding area and volume.
  * @param *tria is a triangle pointer for which normal, area and volume is
  * to be calculated.
  * @returns TS_SUCCESS on success. (always)
  *
  * Calculate normal vector of the triangle (xnorm, ynorm and znorm) and stores
  * information in underlying ts_triangle_data data_structure.
  * information. At the same time
  * triangle area is determined, since we already have the normal and volume of
  * triangular pyramid with given triangle as a base and vesicle centroid as a
  * tip. 
  *
  * Function receives one argument of type ts_triangle. It should be corectly
  * initialized with underlying data structure of type ts_triangle_data. the
  * result is stored in triangle->data->xnorm, triangle->data->ynorm,
  * triangle->data->znorm. Returns TS_SUCCESS on completion. 
  * initialized. The
  * result is stored in triangle->xnorm, triangle->ynorm, triangle->znorm.
  * Area and volume are stored into triangle->area and triangle->volume.
  * Returns TS_SUCCESS on completion. 
  *
  * NOTE: Function uses math.h library. pow function implementation is selected
  * accordind to the setting in genreal.h
  * NOTE: Function uses math.h library. Function pow implementation is selected
  * accordind to the used TS_DOUBLE_* definition set in general.h, so it should
  * be compatible with any type of floating point precision.
  *
  * Example of usage:
  *        triangle_normal_vector(tlist->tria[3]);
  *
  *        Computes normals and stores information into tlist->tria[3]->xnorm,
  *        tlist->tria[3]->ynorm, tlist->tria[3]->znorm.
  *        tlist->tria[3]->ynorm, tlist->tria[3]->znorm tlist->tria[3]->area and
  *        tlist->tria[3]->volume.
  *        
  */
ts_bool triangle_normal_vector(ts_triangle *tria){
@@ -239,13 +249,9 @@
    return TS_SUCCESS;
}






/** @brief Frees the memory allocated for data structure of triangle list
 * (ts_triangle_list)
  * @param *tlist is a pointer to datastructure triangle list to be freed.
  * @returns TS_SUCCESS on success (always).
  *
  * Function frees the memory of ts_triangle_list previously allocated. It
  * accepts one argument, the address of data structure. It destroys all

 src/vertex.c

@@ -6,6 +6,14 @@
#include "bond.h"
#include<stdio.h>

ts_bool vertex_list_assign_id(ts_vertex_list *vlist, ts_uint id){
    ts_uint i;	
    for(i=0;i<vlist->n;i++){
        vlist->vtx[i]->id = id;
    }
    return TS_SUCCESS;
}

ts_vertex_list *init_vertex_list(ts_uint N){    
    ts_int i;
    ts_vertex_list *vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));

 src/vertex.h

@@ -1,6 +1,8 @@
#ifndef _VERTEX_H
#define _VERTEX_H

ts_bool vertex_list_assign_id(ts_vertex_list *vlist, ts_uint id);

/** @brief Creates initial vertex list
 *  
 *  Allocates memory and initializes the vertices.

 src/vertexmove.c

@@ -67,11 +67,18 @@
        }
    }

    //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);
	
// TODO: Maybe faster if checks only nucleus-neighboring cells
// Nucleus penetration check:
    if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
        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
    retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);

    if(retval==TS_FAIL){
        vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
        return TS_FAIL;
@@ -291,12 +298,18 @@
        }
    }

// TODO: Maybe faster if checks only nucleus-neighboring cells
// Nucleus penetration check:
    if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
        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;

 src/vesicle.c

@@ -6,6 +6,8 @@
#include "cell.h"
#include "stdlib.h"
#include "poly.h"
#include "sh.h"
#include "shcomplex.h"

ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
@@ -35,6 +37,7 @@
    triangle_list_free(vesicle->tlist);
    cell_list_free(vesicle->clist);
    poly_list_free(vesicle->poly_list);
    complex_sph_free(vesicle->sphHarmonics);
    free(vesicle);
    return TS_SUCCESS;
}