From 719c9febac2eaff9483fda487b57684afbb59bb2 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Wed, 05 Mar 2014 16:47:50 +0000
Subject: [PATCH] dont know
---
src/cell.c | 2
src/main.c | 19 +
src/vertexmove.h | 1
src/io.c | 388 +++++++++++++++++++++++
src/tape | 19
src/vesicle.c | 4
src/bondflip.c | 73 +--
src/general.h | 45 ++
src/vertexmove.c | 124 +++++++
src/io.h | 2
src/poly.h | 16 +
src/Makefile.am | 12
src/timestep.c | 21 +
src/poly.c | 120 +++++++
src/energy.c | 8
aclocal.m4 | 10
src/frame.c | 36 +
src/energy.h | 1
18 files changed, 809 insertions(+), 92 deletions(-)
diff --git a/aclocal.m4 b/aclocal.m4
index 3094178..da06da5 100644
--- a/aclocal.m4
+++ b/aclocal.m4
@@ -1,4 +1,4 @@
-# generated automatically by aclocal 1.11.6 -*- Autoconf -*-
+# generated automatically by aclocal 1.11.3 -*- Autoconf -*-
# Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
# 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation,
@@ -14,8 +14,8 @@
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
-m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.69],,
-[m4_warning([this file was generated for autoconf 2.69.
+m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],,
+[m4_warning([this file was generated for autoconf 2.68.
You have another version of autoconf. It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])
@@ -38,7 +38,7 @@
[am__api_version='1.11'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version. Point them to the right macro.
-m4_if([$1], [1.11.6], [],
+m4_if([$1], [1.11.3], [],
[AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])
@@ -54,7 +54,7 @@
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
-[AM_AUTOMAKE_VERSION([1.11.6])dnl
+[AM_AUTOMAKE_VERSION([1.11.3])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
diff --git a/src/Makefile.am b/src/Makefile.am
index 10ae3c2..30d71d3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,21 +1,25 @@
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
+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
+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
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
+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
+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
+<<<<<<< HEAD
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
neigh_testdir=../
neigh_test_PROGRAMS = neigh_test
neigh_test_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 neigh_test.c
+=======
+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
+>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768
diff --git a/src/bondflip.c b/src/bondflip.c
index cc29021..4ceaf2b 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -15,10 +15,10 @@
ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn){
/*c Vertex and triangle (lm and lp) indexing for bond flip:
c +----- k-------+ +----- k ------+
-c |lm1 / | \ lp1 | |lm1 / \ lp1 |
+c |lm1 / ^ \ lp1 | |lm1 / \ lp1 |
c | / | \ | | / \ |
c |/ | \ | FLIP |/ lm \ |
-c km lm | lp kp ---> km ---------- kp
+c km lm | lp kp ---> km ---------> kp
c |\ | / | |\ lp / |
c | \ | / | | \ / |
c |lm2 \ | / lp2 | |lm2 \ / lp2 |
@@ -113,7 +113,7 @@
#endif
{
//not accepted, reverting changes
- // fprintf(stderr,"Failed to move, due to MC\n");
+ fprintf(stderr,"Failed to move, due to MC\n");
// ts_flip_bond(km,kp,it,k, bond);
ts_flip_bond(bond);
@@ -142,7 +142,7 @@
return TS_FAIL;
}
}
- // fprintf(stderr,"Success\n");
+ fprintf(stderr,"Success\n");
/* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
@@ -153,49 +153,8 @@
ts_bool ts_flip_bond(ts_bond *bond){
fprintf(stderr,"Called!\n");
ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL;
- ts_uint i; //lmidx, lpidx;
-/*This is no longer necessary! if(k==NULL || it==NULL || km==NULL || kp==NULL){
- fatal("ts_flip_bond: You called me with invalid pointers to vertices",999);
-}*/
-// 1. step. We find lm and lp from k->tristar !
-/* for(i=0;i<it->tristar_no;i++){
- for(j=0;j<k->tristar_no;j++){
- if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
- if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
-== km || it->tristar[i]->vertex[2]== km )){
- lm=it->tristar[i];
- // lmidx=i;
- }
- else
- {
- lp=it->tristar[i];
- // lpidx=i;
- }
+ ts_uint i,j; //lmidx, lpidx;
- }
- }
- }
-if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999);
-
-//we look for important triangles lp1 and lm2.
-
- for(i=0;i<k->tristar_no;i++){
- for(j=0;j<kp->tristar_no;j++){
- if((k->tristar[i] == kp->tristar[j]) && k->tristar[i]!=lp){ //ce gre za skupen trikotnik
- lp1=k->tristar[i];
- }
- }
-}
-
- for(i=0;i<it->tristar_no;i++){
- for(j=0;j<km->tristar_no;j++){
- if((it->tristar[i] == km->tristar[j]) && it->tristar[i]!=lm){ //ce gre za skupen trikotnik
- lm2=it->tristar[i];
- }
- }
- }
-*/
-//Upper block no longer neccessary.
ts_vertex *it=bond->vtx1;
ts_vertex *k=bond->vtx2;
ts_vertex *km=bond->adjvtx[0];
@@ -219,6 +178,22 @@
*/
//find bonds between k and kp (adj[1] vtx)
//find bonds between it and km (adj[0] vtx)
+
+for(i=0;i<3;i++){
+if(lm->neigh[i]==lp) break;
+ for(j=0;j<3;j++){
+ if(lm->neigh[i]->vertex[j]==it) lm2=lm->neigh[i];
+ }
+}
+
+for(i=0;i<3;i++){
+if(lp->neigh[i]==lm) break;
+ for(j=0;j<3;j++){
+ if(lp->neigh[i]->vertex[j]==k) lp1=lp->neigh[i];
+ }
+}
+
+/*
for(i=0;i<it->tristar_no;i++){
if((it->tristar[i]->vertex[0]==km || it->tristar[i]->vertex[1]==km || it->tristar[i]->vertex[2]==km) && (it->tristar[i]->vertex[0]==it || it->tristar[i]->vertex[1]==it || it->tristar[i]->vertex[2]==it)){
lm2=it->tristar[i];
@@ -227,7 +202,11 @@
lp1=it->tristar[i];
}
}
+*/
+if(lm2==NULL ) fatal("ts_flip_bond: Cannot find triangles lm2!",999);
+if(lp1==NULL ) fatal("ts_flip_bond: Cannot find triangles lp1!",999);
if(lm2==NULL || lp1==NULL) fatal("ts_flip_bond: Cannot find triangles lm2 and lp1!",999);
+
/*
//DEBUG TESTING
fprintf(stderr,"1. step: lm, lm2, lp1 and lp found!\n");
@@ -347,9 +326,9 @@
energy_vertex(it);
//Extra steps for new bond data structure
+
bond->adjvtx[0]=k;
bond->adjvtx[1]=it;
-
// END modifications to data structure!
diff --git a/src/cell.c b/src/cell.c
index 32fb9ea..fcd974f 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -147,7 +147,7 @@
}
// Now we check whether we didn't come close to some other vertices in the same
// cell!
- if(cell_occupation>1){
+ if(cell_occupation>0){
for(l=0;l<cell_occupation;l++){
//carefull with this checks!
diff --git a/src/energy.c b/src/energy.c
index a11dce8..abe95d3 100644
--- a/src/energy.c
+++ b/src/energy.c
@@ -20,6 +20,14 @@
}
+inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){
+//TODO: This value to be changed and implemented in data structure:
+ ts_double d_relaxed=1.0;
+ bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2);
+ return TS_SUCCESS;
+};
+
+
inline ts_bool energy_vertex(ts_vertex *vtx){
// ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
// ts_triangle *tristar=vtx->tristar-1;
diff --git a/src/energy.h b/src/energy.h
index 90463d8..8a37b85 100644
--- a/src/energy.h
+++ b/src/energy.h
@@ -2,4 +2,5 @@
#define _ENERGY_H
ts_bool mean_curvature_and_energy(ts_vesicle *vesicle);
inline ts_bool energy_vertex(ts_vertex *vtx);
+inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly);
#endif
diff --git a/src/frame.c b/src/frame.c
index 5cfede9..d2d7e10 100644
--- a/src/frame.c
+++ b/src/frame.c
@@ -3,7 +3,7 @@
#include "cell.h"
#include "frame.h"
ts_bool centermass(ts_vesicle *vesicle){
- ts_uint i, n=vesicle->vlist->n;
+ ts_uint i,j, n=vesicle->vlist->n;
ts_vertex **vtx=vesicle->vlist->vtx;
vesicle->cm[0]=0;
vesicle->cm[1]=0;
@@ -23,6 +23,15 @@
vtx[i]->z-=vesicle->cm[2];
}
+ 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];
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->y-=vesicle->cm[1];
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->z-=vesicle->cm[2];
+ }
+ }
+
+
vesicle->cm[0]=0;
vesicle->cm[1]=0;
vesicle->cm[2]=0;
@@ -31,11 +40,11 @@
}
ts_bool cell_occupation(ts_vesicle *vesicle){
- ts_uint i,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);
+ 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);
@@ -45,10 +54,19 @@
// vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx];
cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]);
-
}
+
+ 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]);
+ }
+ }
+
+
+
//fprintf(stderr, "Bil sem tu\n");
- if(dcell);
- if(shift);
+ //if(dcell);
+ //if(shift);
return TS_SUCCESS;
}
diff --git a/src/general.h b/src/general.h
index 304da95..643aab1 100644
--- a/src/general.h
+++ b/src/general.h
@@ -8,10 +8,11 @@
* @file header.h
* @author Samo Penic
* @date 5.3.2001
- *
+ *
* Header file for general inclusion in all the code, defining data structures
* and general constans. All datatypes used in the code is also defined here.
*
+ * Miha: branch trisurf-polyel
*/
/* Defines */
@@ -130,9 +131,12 @@
ts_uint coord_type;
} ts_coord;
+<<<<<<< HEAD
+=======
+>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768
/** @brief Data structure of all data connected to a vertex
*
* ts_vertex holds the data for one single point (bead, vertex). To understand how to use it
@@ -161,6 +165,7 @@
ts_double projArea;
ts_double relR;
ts_double solAngle;
+ struct ts_poly *grafted_poly;
};
typedef struct ts_vertex ts_vertex;
@@ -173,7 +178,7 @@
typedef struct ts_vertex_list ts_vertex_list;
struct ts_bond {
- ts_uint idx;
+ ts_uint idx;
ts_vertex *vtx1;
ts_vertex *vtx2;
ts_vertex *adjvtx[2];
@@ -181,6 +186,7 @@
ts_double bond_length;
ts_double bond_length_dual;
ts_bool tainted;
+ ts_double energy;
};
typedef struct ts_bond ts_bond;
@@ -237,18 +243,39 @@
+struct ts_poly {
+ ts_vertex_list *vlist;
+ ts_bond_list *blist;
+ ts_vertex *grafted_vtx;
+ ts_double k;
+};
+typedef struct ts_poly ts_poly;
+
+
+struct ts_poly_list {
+ ts_uint n;
+ ts_poly **poly;
+};
+typedef struct ts_poly_list ts_poly_list;
+
+
+
+
typedef struct {
ts_vertex_list *vlist;
ts_bond_list *blist;
ts_triangle_list *tlist;
- ts_cell_list *clist;
+ ts_cell_list *clist;
ts_uint nshell;
- ts_double bending_rigidity;
- ts_double dmax;
- ts_double stepsize;
- ts_double cm[3];
- ts_double volume;
- ts_spharm *sphHarmonics;
+ ts_double bending_rigidity;
+ ts_double dmax;
+ ts_double stepsize;
+ ts_double cm[3];
+ ts_double volume;
+ ts_spharm *sphHarmonics;
+
+ ts_poly_list *poly_list;
+ ts_double spring_constant;
} ts_vesicle;
diff --git a/src/io.c b/src/io.c
index f35b57d..b2659db 100644
--- a/src/io.c
+++ b/src/io.c
@@ -9,6 +9,329 @@
#include <sys/types.h>
#include <dirent.h>
#include "initial_distribution.h"
+#include "poly.h"
+
+
+
+/** DUMP STATE TO DISK DRIVE **/
+
+ts_bool dump_state(ts_vesicle *vesicle){
+
+ /* save current state with wrong pointers. Will fix that later */
+ ts_uint i,j,k;
+ ts_ulong off;
+ FILE *fh=fopen("dump.bin","wb");
+
+ /* dump vesicle */
+ fwrite(vesicle, sizeof(vesicle),1,fh);
+ /* dump vertex list */
+ fwrite(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
+ /* dump bond list */
+ fwrite(vesicle->blist, sizeof(ts_bond_list),1,fh);
+ /* dump triangle list */
+ fwrite(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
+ /* dump cell list */
+ fwrite(vesicle->clist, sizeof(ts_cell_list),1,fh);
+ /* dump poly list */
+ fwrite(vesicle->poly_list, sizeof(ts_poly_list),1,fh);
+ /* level 1 complete */
+
+ /*dump vertices*/
+ for(i=0;i<vesicle->vlist->n;i++){
+ fwrite(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
+ /* dump pointer offsets for:
+ neigh
+ bond
+ tria
+ cell is ignored
+ */
+ for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
+ off=(ts_ulong)(vesicle->vlist->vtx[i]->neigh[j]-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+ for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
+ off=(ts_ulong)(vesicle->vlist->vtx[i]->bond[j]-vesicle->blist->bond[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+ for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
+ off=(ts_ulong)(vesicle->vlist->vtx[i]->tristar[j]-vesicle->tlist->tria[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+// off=(ts_ulong)(vesicle->vlist->vtx[i]->cell-vesicle->clist->cell[0]);
+// fwrite(&off,sizeof(ts_ulong),1,fh);
+
+ }
+
+ /*dump bonds*/
+ for(i=0;i<vesicle->blist->n;i++){
+ fwrite(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
+ /* dump pointer offsets for vtx1 and vtx2 */
+ off=(ts_ulong)(vesicle->blist->bond[i]->vtx1-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ off=(ts_ulong)(vesicle->blist->bond[i]->vtx2-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+
+ /*dump triangles*/
+ for(i=0;i<vesicle->tlist->n;i++){
+ fwrite(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
+ /* dump pointer offsets for vertex */
+ off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[0]-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[1]-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ off=(ts_ulong)(vesicle->tlist->tria[i]->vertex[2]-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ /* dump pointer offsets for neigh */
+ for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
+ off=(ts_ulong)(vesicle->tlist->tria[i]->neigh[j]-vesicle->tlist->tria[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+ }
+
+ /*dump cells and its contents is not dumped but recalculated on restore*/
+
+/* for(i=0;i<vesicle->clist->cellno;i++){
+ fwrite(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
+ for(j=0;j<vesicle->clist->cell[i]->nvertex;j++){
+ off=(ts_ulong)(vesicle->clist->cell[i]->vertex[j]-vesicle->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+ }
+*/
+
+
+ /*dump polymeres */
+ for(i=0;i<vesicle->poly_list->n;i++){
+ fwrite(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
+ }
+
+ /* dump poly vertex(monomer) list*/
+ for(i=0;i<vesicle->poly_list->n;i++){
+ for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
+ fwrite(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
+ /* dump offset for neigh and bond */
+ for(k=0;k<vesicle->poly_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(&off,sizeof(ts_ulong),1,fh);
+ }
+ for(k=0;k<vesicle->poly_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(&off,sizeof(ts_ulong),1,fh);
+ }
+ }
+ }
+ /* dump poly bonds between monomers list*/
+ for(i=0;i<vesicle->poly_list->n;i++){
+ for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
+ fwrite(vesicle->poly_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(&off,sizeof(ts_ulong),1,fh);
+ off=(ts_ulong)(vesicle->poly_list->poly[i]->blist->bond[j]->vtx2-vesicle->poly_list->poly[i]->vlist->vtx[0]);
+ fwrite(&off,sizeof(ts_ulong),1,fh);
+ }
+ }
+
+/* pointer offsets for fixing the restored pointers */
+/* need pointers for
+ vlist->vtx
+ blist->bond
+ tlist->tria
+ clist->cell
+ poly_list->poly
+ and for each poly:
+ poly_list->poly->vtx
+ poly_list->poly->bond
+*/
+
+ fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
+ fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
+ fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
+ fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
+ fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
+ for(i=0;i<vesicle->poly_list->n;i++){
+ fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
+ fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
+ }
+ fclose(fh);
+ return TS_SUCCESS;
+}
+
+
+/** RESTORE DUMP FROM DISK **/
+ts_vesicle *restore_state(){
+ ts_uint i,j,k;
+ FILE *fh=fopen("dump.bin","rb");
+ ts_uint retval;
+ ts_ulong off;
+
+/* we restore all the data from the dump */
+ /* restore vesicle */
+ ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
+ retval=fread(vesicle, sizeof(vesicle),1,fh);
+ /* restore vertex list */
+ vesicle->vlist=(ts_vertex_list *)malloc(sizeof(ts_vertex_list));
+ retval=fread(vesicle->vlist, sizeof(ts_vertex_list),1,fh);
+ /* restore bond list */
+ vesicle->blist=(ts_bond_list *)malloc(sizeof(ts_bond_list));
+ retval=fread(vesicle->blist, sizeof(ts_bond_list),1,fh);
+ /* restore triangle list */
+ vesicle->tlist=(ts_triangle_list *)malloc(sizeof(ts_triangle_list));
+ retval=fread(vesicle->tlist, sizeof(ts_triangle_list),1,fh);
+ /* restore cell list */
+ vesicle->clist=(ts_cell_list *)malloc(sizeof(ts_cell_list));
+ retval=fread(vesicle->clist, sizeof(ts_cell_list),1,fh);
+ /* 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);
+ /* level 1 complete */
+
+ /*restore vertices*/
+ vesicle->vlist->vtx=(ts_vertex **)calloc(vesicle->vlist->n,sizeof(ts_vertex *));
+ for(i=0;i<vesicle->vlist->n;i++){
+ vesicle->vlist->vtx[i]=(ts_vertex *)malloc(sizeof(ts_vertex));
+ retval=fread(vesicle->vlist->vtx[i],sizeof(ts_vertex),1,fh);
+ /*restore neigh, bond, tristar. Ignoring cell */
+ vesicle->vlist->vtx[i]->neigh=(ts_vertex **)calloc(vesicle->vlist->vtx[i]->neigh_no, sizeof(ts_vertex *));
+ for(j=0;j<vesicle->vlist->vtx[i]->neigh_no;j++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->vlist->vtx[i]->neigh[j]=vesicle->vlist->vtx[0]+off;
+ }
+ vesicle->vlist->vtx[i]->bond=(ts_bond **)calloc(vesicle->vlist->vtx[i]->bond_no, sizeof(ts_bond *));
+ for(j=0;j<vesicle->vlist->vtx[i]->bond_no;j++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->vlist->vtx[i]->bond[j]=vesicle->blist->bond[0]+off;
+ }
+
+ vesicle->vlist->vtx[i]->tristar=(ts_triangle **)calloc(vesicle->vlist->vtx[i]->tristar_no, sizeof(ts_triangle *));
+ for(j=0;j<vesicle->vlist->vtx[i]->tristar_no;j++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->vlist->vtx[i]->tristar[j]=vesicle->tlist->tria[0]+off;
+ }
+
+ }
+
+ /*restore bonds*/
+ vesicle->blist->bond=(ts_bond **)calloc(vesicle->blist->n,sizeof(ts_bond *));
+ for(i=0;i<vesicle->blist->n;i++){
+ vesicle->blist->bond[i]=(ts_bond *)malloc(sizeof(ts_bond));
+ retval=fread(vesicle->blist->bond[i],sizeof(ts_bond),1,fh);
+ /* restore vtx1 and vtx2 */
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->blist->bond[i]->vtx1=vesicle->vlist->vtx[0]+off;
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->blist->bond[i]->vtx2=vesicle->vlist->vtx[0]+off;
+ }
+
+ /*restore triangles*/
+ vesicle->tlist->tria=(ts_triangle **)calloc(vesicle->tlist->n,sizeof(ts_triangle *));
+ for(i=0;i<vesicle->tlist->n;i++){
+ vesicle->tlist->tria[i]=(ts_triangle *)malloc(sizeof(ts_triangle));
+ retval=fread(vesicle->tlist->tria[i],sizeof(ts_triangle),1,fh);
+ /* restore pointers for vertices */
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->tlist->tria[i]->vertex[0]=vesicle->vlist->vtx[0]+off;
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->tlist->tria[i]->vertex[1]=vesicle->vlist->vtx[0]+off;
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->tlist->tria[i]->vertex[2]=vesicle->vlist->vtx[0]+off;
+ /* restore pointers for neigh */
+ for(j=0;j<vesicle->tlist->tria[i]->neigh_no;j++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->tlist->tria[i]->neigh[j]=vesicle->tlist->tria[0]+off;
+ }
+
+ }
+
+ /*restore cells */
+/*TODO: do we need to recalculate cells here? */
+/* vesicle->clist->cell=(ts_cell **)malloc(vesicle->clist->cellno*sizeof(ts_cell *));
+ for(i=0;i<vesicle->clist->cellno;i++){
+ vesicle->clist->cell[i]=(ts_cell *)malloc(sizeof(ts_cell));
+ retval=fread(vesicle->clist->cell[i],sizeof(ts_cell),1,fh);
+ }
+*/
+ /*restore polymeres */
+ vesicle->poly_list->poly = (ts_poly **)calloc(vesicle->poly_list->n,sizeof(ts_poly *));
+ for(i=0;i<vesicle->poly_list->n;i++){
+ vesicle->poly_list->poly[i]=(ts_poly *)calloc(1,sizeof(ts_poly));
+ retval=fread(vesicle->poly_list->poly[i],sizeof(ts_poly),1,fh);
+ }
+
+ /* restore poly vertex(monomer) list*/
+ for(i=0;i<vesicle->poly_list->n;i++){
+ vesicle->poly_list->poly[i]->vlist->vtx=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->n,sizeof(ts_vertex *));
+ for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
+ vesicle->poly_list->poly[i]->vlist->vtx[j]=(ts_vertex *)malloc(sizeof(ts_vertex));
+ retval=fread(vesicle->poly_list->poly[i]->vlist->vtx[j],sizeof(ts_vertex),1,fh);
+
+ /* restore neigh and bonds */
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh=(ts_vertex **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no, sizeof(ts_vertex *));
+ for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh_no;k++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->neigh[k]=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
+ }
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no, sizeof(ts_bond *));
+ for(k=0;k<vesicle->poly_list->poly[i]->vlist->vtx[j]->bond_no;k++){
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->poly_list->poly[i]->vlist->vtx[j]->bond[k]=vesicle->poly_list->poly[i]->blist->bond[0]+off;
+ }
+
+ }
+ }
+ /* restore poly bonds between monomers list*/
+ for(i=0;i<vesicle->poly_list->n;i++){
+ vesicle->poly_list->poly[i]->blist->bond=(ts_bond **)calloc(vesicle->poly_list->poly[i]->blist->n,sizeof(ts_bond *));
+ for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
+ vesicle->poly_list->poly[i]->blist->bond[j]=(ts_bond *)malloc(sizeof(ts_bond));
+ retval=fread(vesicle->poly_list->poly[i]->blist->bond[j],sizeof(ts_bond),1,fh);
+ /* restore vtx1 and vtx2 */
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->poly_list->poly[i]->blist->bond[j]->vtx1=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
+ retval=fread(&off,sizeof(ts_ulong),1,fh);
+ vesicle->poly_list->poly[i]->blist->bond[j]->vtx2=vesicle->poly_list->poly[i]->vlist->vtx[0]+off;
+ }
+ }
+
+
+ /* recalculate pointers using saved information on pointer addresses */
+ /* following pointers need to be recalculated:
+ all vtx->neigh
+ all vtx->tristar
+ all vtx->bond
+ vtx->cell
+ vtx->grafted_poly
+ bond->vtx1
+ bond->vtx2
+ tria->vertex[*]
+ tria->neigh
+ cell->vertex
+ poly->grafted_vtx
+ poly->vlist->vtx->neigh
+ poly->vlist->vtx->bond
+ poly->vlist->vtx->cell
+ */
+
+ /* read pointers from disk */
+/* fwrite(vesicle->vlist->vtx, sizeof(ts_vertex *),1,fh);
+ fwrite(vesicle->blist->bond, sizeof(ts_bond *),1,fh);
+ fwrite(vesicle->tlist->tria, sizeof(ts_triangle *),1,fh);
+ fwrite(vesicle->clist->cell, sizeof(ts_cell *),1,fh);
+ fwrite(vesicle->poly_list->poly, sizeof(ts_poly *),1,fh);
+ for(i=0;i<vesicle->poly_list->n;i++){
+ fwrite(vesicle->poly_list->poly[i]->vlist->vtx, sizeof(ts_vertex *),1,fh);
+ fwrite(vesicle->poly_list->poly[i]->blist->bond, sizeof(ts_bond *),1,fh);
+ }
+*/
+
+
+ if(retval);
+ fclose(fh);
+ return vesicle;
+}
+
+
ts_bool print_vertex_list(ts_vertex_list *vlist){
ts_uint i;
@@ -206,7 +529,7 @@
ts_vertex_list *vlist=vesicle->vlist;
ts_bond_list *blist=vesicle->blist;
ts_vertex **vtx=vlist->vtx;
- ts_uint i;
+ ts_uint i,j;
char filename[255];
FILE *fh;
@@ -217,29 +540,69 @@
return TS_FAIL;
}
/* Here comes header of the file */
+
+ //find number of extra vtxs and bonds of polymeres
+ ts_uint monono=0, polyno=0;
+ ts_bool poly=0;
+ if(vesicle->poly_list!=NULL){
+ if(vesicle->poly_list->poly[0]!=NULL){
+ polyno=vesicle->poly_list->n;
+ monono=vesicle->poly_list->poly[0]->vlist->n;
+ poly=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, blist->n);
+ fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno, blist->n+monono*polyno);
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);
}
+ //polymeres
+ if(poly){
+ for(i=0;i<vesicle->poly_list->n;i++){
+ for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
+ fprintf(fh,"%u ", vesicle->poly_list->poly[i]->vlist->vtx[j]->idx);
+ }
+ }
+ }
fprintf(fh,"</DataArray>\n</PointData>\n<CellData>\n</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
for(i=0;i<vlist->n;i++){
fprintf(fh,"%e %e %e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
+ }
+ //polymeres
+ if(poly){
+ 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 );
+ }
+ }
}
fprintf(fh,"</DataArray>\n</Points>\n<Cells>\n<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">");
for(i=0;i<blist->n;i++){
fprintf(fh,"%u %u\n",blist->bond[i]->vtx1->idx,blist->bond[i]->vtx2->idx);
}
+ //polymeres
+ if(poly){
+ for(i=0;i<vesicle->poly_list->n;i++){
+ for(j=0;j<vesicle->poly_list->poly[i]->blist->n;j++){
+ fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->blist->bond[j]->vtx1->idx,vesicle->poly_list->poly[i]->blist->bond[j]->vtx2->idx);
+ }
+ //grafted bonds
+ fprintf(fh,"%u %u\n", vesicle->poly_list->poly[i]->grafted_vtx->idx, vesicle->poly_list->poly[i]->vlist->vtx[0]->idx);
+ }
+
+ }
+
+
fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
- for (i=2;i<blist->n*2+1;i+=2){
+ for (i=2;i<(blist->n+monono*polyno)*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;i++){
+ for (i=0;i<blist->n+monono*polyno;i++){
fprintf(fh,"3 ");
}
@@ -295,17 +658,22 @@
ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){
- long int nshell=17,ncxmax=60, ncymax=60, nczmax=60; // THIS IS DUE TO CONFUSE BUG!
+ long int nshell=17,ncxmax=60, ncymax=60, nczmax=60, npoly=10, nmono=20; // THIS IS DUE TO CONFUSE BUG!
char *buf=malloc(255*sizeof(char));
long int brezveze0=1;
long int brezveze1=1;
long int brezveze2=1;
- ts_double xk0=25.0, dmax=1.67,stepsize=0.15;
+ ts_double xk0=25.0, dmax=1.67,stepsize=0.15,kspring=800.0;
long int iter=1000, init=1000, mcsw=1000;
+
+
cfg_opt_t opts[] = {
CFG_SIMPLE_INT("nshell", &nshell),
+ CFG_SIMPLE_INT("npoly", &npoly),
+ CFG_SIMPLE_INT("nmono", &nmono),
CFG_SIMPLE_FLOAT("dmax", &dmax),
CFG_SIMPLE_FLOAT("xk0",&xk0),
+ CFG_SIMPLE_FLOAT("k_spring",&kspring),
CFG_SIMPLE_FLOAT("stepsize",&stepsize),
CFG_SIMPLE_INT("nxmax", &ncxmax),
CFG_SIMPLE_INT("nymax", &ncymax),
@@ -335,6 +703,11 @@
*inititer=init;
*mcsweeps=mcsw;
vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize);
+ vesicle->poly_list=init_poly_list(npoly,nmono, vesicle->vlist);
+ vesicle->spring_constant=kspring;
+ poly_assign_spring_const(vesicle);
+
+
vesicle->nshell=nshell;
vesicle->dmax=dmax*dmax;
vesicle->bending_rigidity=xk0;
@@ -347,6 +720,9 @@
cfg_free(cfg);
free(buf);
// fprintf(stderr,"NSHELL=%u\n",vesicle->nshell);
+
+
+
return vesicle;
}
diff --git a/src/io.h b/src/io.h
index 20eab38..7ac2160 100644
--- a/src/io.h
+++ b/src/io.h
@@ -50,4 +50,6 @@
ts_bool write_master_xml_file(ts_char *filename);
ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations);
+ts_bool dump_state(ts_vesicle *vesicle);
+ts_vesicle *restore_state();
#endif
diff --git a/src/main.c b/src/main.c
index 7e12c1c..941daaf 100644
--- a/src/main.c
+++ b/src/main.c
@@ -10,6 +10,7 @@
#include "initial_distribution.h"
#include "frame.h"
#include "timestep.h"
+#include "poly.h"
/** Entrance function to the program
* @param argv is a number of parameters used in program call (including the program name
@@ -64,6 +65,24 @@
printf("Tests complete.\n");
*/
vesicle=parsetape(&mcsweeps, &inititer, &iterations);
+
+/*Testing */
+//vesicle->poly_list=init_poly_list(1400,20,vesicle->vlist);
+
+//poly_list_free(vesicle->poly_list);
+/*End testing*/
+printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
+dump_state(vesicle);
+vesicle->vlist->vtx[123]->x=123.0;
+printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
+write_vertex_xml_file(vesicle,0);
+vesicle=restore_state();
+printf("Stevilo vertexov je %d\n",vesicle->vlist->n);
+printf("Vertex %d ima x komponento %e\n",123,vesicle->vlist->vtx[123]->x);
+write_vertex_xml_file(vesicle,1);
+write_master_xml_file("test.pvd");
+return 0;
+
run_simulation(vesicle, mcsweeps, inititer, iterations);
write_master_xml_file("test.pvd");
diff --git a/src/poly.c b/src/poly.c
new file mode 100644
index 0000000..1ddcec7
--- /dev/null
+++ b/src/poly.c
@@ -0,0 +1,120 @@
+#include"general.h"
+#include"poly.h"
+#include<stdlib.h>
+#include"vertex.h"
+#include"bond.h"
+#include<math.h>
+#include"energy.h"
+
+ts_bool poly_assign_spring_const(ts_vesicle *vesicle){
+ ts_uint i;
+
+ for(i=0;i<vesicle->poly_list->n;i++){
+ vesicle->poly_list->poly[i]->k = vesicle->spring_constant;
+ }
+
+ return TS_SUCCESS;
+}
+
+ts_poly *init_poly(ts_uint n, ts_vertex *grafted_vtx){
+ 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;
+
+ ts_uint i;
+ for(i=0;i<n-1;i++){
+ vtx_add_cneighbour(poly->blist, poly->vlist->vtx[i], poly->vlist->vtx[i+1]);
+ vtx_add_neighbour(poly->vlist->vtx[i+1], poly->vlist->vtx[i]);
+ }
+
+ for(i=0;i<poly->blist->n;i++){
+ poly->blist->bond[i]->bond_length=sqrt(vtx_distance_sq(poly->blist->bond[i]->vtx1,poly->blist->bond[i]->vtx2));
+ bond_energy(poly->blist->bond[i],poly);
+ }
+
+ return poly;
+}
+
+
+ts_poly_list *init_poly_list(ts_uint n_poly, ts_uint n_mono, ts_vertex_list *vlist){
+ 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;
+
+ 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++;
+ }
+ }
+
+ 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++){
+
+ 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;
+
+ 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);
+ }
+ }
+
+ //index correction for polymeres. Important, since each vtx has to have unique id
+ idx=vlist->n;
+ for(i=0;i<n_poly;i++){
+ for(j=0;j<n_mono;j++,idx++){
+
+ poly_list->poly[i]->vlist->vtx[j]->idx=idx;
+
+ }
+ }
+
+
+ return poly_list;
+}
+
+
+ts_bool poly_free(ts_poly *poly){
+
+ if (poly->grafted_vtx!=NULL){
+ poly->grafted_vtx->grafted_poly=NULL;
+ }
+ vtx_list_free(poly->vlist);
+ bond_list_free(poly->blist);
+ free(poly);
+
+ return TS_SUCCESS;
+}
+
+ts_bool poly_list_free(ts_poly_list *poly_list){
+ ts_uint i;
+
+ for(i=0;i<poly_list->n;i++){
+ poly_free(poly_list->poly[i]);
+ }
+ free(poly_list->poly);
+ free(poly_list);
+
+ return TS_SUCCESS;
+}
diff --git a/src/poly.h b/src/poly.h
new file mode 100644
index 0000000..9bdec52
--- /dev/null
+++ b/src/poly.h
@@ -0,0 +1,16 @@
+#ifndef _POLY_H
+#define _POLY_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_bool poly_free(ts_poly *poly);
+
+ts_bool poly_list_free(ts_poly_list *poly_list);
+
+ts_bool poly_assign_spring_const(ts_vesicle *vesicle);
+
+#endif
diff --git a/src/tape b/src/tape
index a879c23..914bf81 100644
--- a/src/tape
+++ b/src/tape
@@ -1,13 +1,20 @@
####### Vesicle definitions ###########
# nshell is a number of divisions of dipyramid
nshell=17
-# dmax is the square of the bond length
-dmax=1.67
+# dmax is the max. bond length
+dmax=1.7
# bending rigidity of the membrane
xk0=25.0
# max step size
stepsize=0.15
+####### Polymer definitions ###########
+# npoly is a number of polymers attached to npoly distinct vertices on vesicle
+npoly=30
+# nmono is a number of monomers in each polymer
+nmono=10
+# Spring constant between monomers of the polymer
+k_spring=800
####### Cell definitions ############
nxmax=60
@@ -17,11 +24,15 @@
####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
-mcsweeps=100
+mcsweeps=5000
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
-inititer=100
+inititer=1
#how many records do you want on the disk iteration are there in a run?
+<<<<<<< HEAD
iterations=10
+=======
+iterations=10000
+>>>>>>> e9eab4fb2f72383a1a2adbe5f09f7bbd1fd45768
#shut up if we are using cluster!!!
diff --git a/src/timestep.c b/src/timestep.c
index 9e5e61c..d717abb 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -8,6 +8,7 @@
#include "bondflip.h"
#include "frame.h"
#include "io.h"
+
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){
ts_uint i, j;
@@ -20,6 +21,7 @@
}
centermass(vesicle);
cell_occupation(vesicle);
+ ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps);
if(i>inititer){
write_vertex_xml_file(vesicle,i-inititer);
}
@@ -30,7 +32,7 @@
ts_bool single_timestep(ts_vesicle *vesicle){
ts_bool retval;
ts_double rnvec[3];
- ts_uint i, b;
+ ts_uint i,j,b;
for(i=0;i<vesicle->vlist->n;i++){
rnvec[0]=drand48();
rnvec[1]=drand48();
@@ -39,7 +41,7 @@
}
// ts_int cnt=0;
- for(i=0;i<vesicle->vlist->n;i++){
+ for(i=0;i<3*vesicle->vlist->n;i++){
//why is rnvec needed in bondflip?
/* rnvec[0]=drand48();
rnvec[1]=drand48();
@@ -50,8 +52,19 @@
//call single_bondflip_timestep...
retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
// if(retval==TS_SUCCESS) cnt++;
- }
-// printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);
+ }
+
+ 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);
+ }
+
+ }
+
+// 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 7047450..2a26bbe 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -13,8 +13,7 @@
#include "vertexmove.h"
#include <string.h>
-ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double
-*rn){
+ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
ts_uint i;
ts_double dist;
ts_bool retval;
@@ -58,6 +57,16 @@
return TS_FAIL;
}
}
+
+// Distance with grafted poly-vertex check:
+ if(vtx->grafted_poly!=NULL){
+ dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
+ if(dist<1.0 || dist>vesicle->dmax) {
+ vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
+ return TS_FAIL;
+ }
+ }
+
//self avoidance check with distant vertices
cellidx=vertex_self_avoidance(vesicle, vtx);
//check occupation number
@@ -88,6 +97,13 @@
energy_vertex(vtx->neigh->vtx[i]);
delta_energy+=vtx->neigh->vtx[i]->xk*(vtx->neigh->vtx[i]->energy-oenergy);
}
+/* No poly-bond energy for now!
+ if(vtx->grafted_poly!=NULL){
+ delta_energy+=
+ (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
+ pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
+ }
+*/
// fprintf(stderr, "DE=%f\n",delta_energy);
//MONTE CARLOOOOOOOO
if(delta_energy>=0){
@@ -126,3 +142,107 @@
return TS_SUCCESS;
}
+
+ts_bool single_poly_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;
+ //This will hold all the information of vtx and its neighbours
+ ts_vertex backupvtx;
+// ts_bond backupbond[2];
+ 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->neigh_no;i++){
+ dist=vtx_distance_sq(vtx,vtx->neigh[i]);
+ if(dist<1.0 || dist>vesicle->dmax) {
+ vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
+ return TS_FAIL;
+ }
+ }
+
+// Distance with grafted vesicle-vertex check:
+ if(vtx==poly->vlist->vtx[0]){
+ dist=vtx_distance_sq(vtx,poly->grafted_vtx);
+ if(dist<1.0 || dist>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;
+ }
+
+
+ //if all the tests are successful, then energy for vtx and neighbours is calculated
+/* Energy ignored for now!
+ delta_energy=0;
+ for(i=0;i<vtx->bond_no;i++){
+ memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
+
+ vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
+ bond_energy(vtx->bond[i],poly);
+ delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
+ }
+
+ if(vtx==poly->vlist->vtx[0]){
+ delta_energy+=
+ (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
+ pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *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->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 0c24d54..6387ab8 100644
--- a/src/vertexmove.h
+++ b/src/vertexmove.h
@@ -1,3 +1,4 @@
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);
diff --git a/src/vesicle.c b/src/vesicle.c
index d2caf4b..43b16b5 100644
--- a/src/vesicle.c
+++ b/src/vesicle.c
@@ -5,10 +5,11 @@
#include "bond.h"
#include "cell.h"
#include "stdlib.h"
+#include "poly.h"
ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
- ts_vesicle *vesicle=(ts_vesicle *)malloc(sizeof(ts_vesicle));
+ ts_vesicle *vesicle=(ts_vesicle *)calloc(1,sizeof(ts_vesicle));
vesicle->vlist=init_vertex_list(N);
vesicle->blist=init_bond_list();
vesicle->tlist=init_triangle_list();
@@ -33,6 +34,7 @@
bond_list_free(vesicle->blist);
triangle_list_free(vesicle->tlist);
cell_list_free(vesicle->clist);
+ poly_list_free(vesicle->poly_list);
free(vesicle);
return TS_SUCCESS;
}
--
Gitblit v1.9.3