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