From b36841fb1903e9945cb4a1dbc4dd7761556e07e2 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Mon, 28 Apr 2014 10:40:08 +0000
Subject: [PATCH] Merge branch 'trout-rbc'
---
src/io.c | 35 ++
src/stats.c | 83 +++++
src/stats.h | 6
src/tape | 19
src/vesicle.c | 3
config.h.in | 6
src/initial_distribution.c | 16
src/general.h | 5
src/vertexmove.c | 25 +
src/io.h | 4
src/Makefile.am | 28
src/timestep.c | 77 ++++
src/shcomplex.c | 172 +++++++++++
src/timestep.h | 2
config.h | 6
src/shcomplex.h | 9
src/triangle.c | 96 +++---
src/cell.c | 15
src/main.c | 13
configure.ac | 3
src/spherical_trisurf_ff.c | 29 -
src/poly.c | 26 +
src/vertex.c | 8
src/sh.h | 2
src/spherical_trisurf.c | 28 -
aclocal.m4 | 122 +++++++
src/frame.c | 1
src/sh.c | 47 +++
src/vertex.h | 2
29 files changed, 729 insertions(+), 159 deletions(-)
diff --git a/aclocal.m4 b/aclocal.m4
index 72a1419..3a9a040 100644
--- a/aclocal.m4
+++ b/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.
diff --git a/config.h b/config.h
index 9a07c00..ac9c24a 100644
--- a/config.h
+++ b/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
diff --git a/config.h.in b/config.h.in
index c692c55..f865b7a 100644
--- a/config.h.in
+++ b/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
diff --git a/configure.ac b/configure.ac
index 791a930..2ebed26 100644
--- a/configure.ac
+++ b/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])
diff --git a/src/Makefile.am b/src/Makefile.am
index e899cd3..20bae62 100644
--- a/src/Makefile.am
+++ b/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
diff --git a/src/cell.c b/src/cell.c
index 97f2fe7..ed86b16 100644
--- a/src/cell.c
+++ b/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;
}
}
}
diff --git a/src/frame.c b/src/frame.c
index aeaef21..237aa87 100644
--- a/src/frame.c
+++ b/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;
diff --git a/src/general.h b/src/general.h
index 79da2ba..8165505 100644
--- a/src/general.h
+++ b/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;
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 3c27e55..bff9f05 100644
--- a/src/initial_distribution.c
+++ b/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;
}
diff --git a/src/io.c b/src/io.c
index b222081..d001b7e 100644
--- a/src/io.c
+++ b/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;
diff --git a/src/io.h b/src/io.h
index bfc6173..699c73e 100644
--- a/src/io.h
+++ b/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);
diff --git a/src/main.c b/src/main.c
index 588e277..25f1343 100644
--- a/src/main.c
+++ b/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){
diff --git a/src/poly.c b/src/poly.c
index 5cfc6e0..b2b13f9 100644
--- a/src/poly.c
+++ b/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);
}
}
}
diff --git a/src/sh.c b/src/sh.c
index 0c4b06d..f5c310a 100644
--- a/src/sh.c
+++ b/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;
+}
diff --git a/src/sh.h b/src/sh.h
index 7853139..dabece2 100644
--- a/src/sh.h
+++ b/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
diff --git a/src/shcomplex.c b/src/shcomplex.c
new file mode 100644
index 0000000..7d8c9d8
--- /dev/null
+++ b/src/shcomplex.c
@@ -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;
+}
diff --git a/src/shcomplex.h b/src/shcomplex.h
new file mode 100644
index 0000000..fdd0eb0
--- /dev/null
+++ b/src/shcomplex.h
@@ -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
diff --git a/src/spherical_trisurf.c b/src/spherical_trisurf.c
index 7fa40f9..0669981 100644
--- a/src/spherical_trisurf.c
+++ b/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;
-}
diff --git a/src/spherical_trisurf_ff.c b/src/spherical_trisurf_ff.c
index 867f37d..157b0e4 100644
--- a/src/spherical_trisurf_ff.c
+++ b/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;
-}
diff --git a/src/stats.c b/src/stats.c
new file mode 100644
index 0000000..077781e
--- /dev/null
+++ b/src/stats.c
@@ -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;
+}
diff --git a/src/stats.h b/src/stats.h
new file mode 100644
index 0000000..223a99b
--- /dev/null
+++ b/src/stats.h
@@ -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
diff --git a/src/tape b/src/tape
index bd62b6c..0801375 100644
--- a/src/tape
+++ b/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
diff --git a/src/timestep.c b/src/timestep.c
index d631d42..a8f0f94 100644
--- a/src/timestep.c
+++ b/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;
}
diff --git a/src/timestep.h b/src/timestep.h
index dd18481..0d8dc74 100644
--- a/src/timestep.h
+++ b/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
diff --git a/src/triangle.c b/src/triangle.c
index 311e138..7dd5a6e 100644
--- a/src/triangle.c
+++ b/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
diff --git a/src/vertex.c b/src/vertex.c
index 38647aa..f733d49 100644
--- a/src/vertex.c
+++ b/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));
diff --git a/src/vertex.h b/src/vertex.h
index 97775d2..cca12ee 100644
--- a/src/vertex.h
+++ b/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.
diff --git a/src/vertexmove.c b/src/vertexmove.c
index da8e70f..6bc730e 100644
--- a/src/vertexmove.c
+++ b/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;
diff --git a/src/vesicle.c b/src/vesicle.c
index 43b16b5..f63eb12 100644
--- a/src/vesicle.c
+++ b/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;
}
--
Gitblit v1.9.3