From c2c6364dc7d23f21dabd2179e775e93ba4eeb252 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Tue, 06 Dec 2016 13:45:07 +0000
Subject: [PATCH] Merging again because I was merging in detached state

---
 src/bondflip.h             |    2 
 src/main.c                 |    2 
 src/snapshot.c             |    2 
 src/io.c                   |   95 ++++++++
 src/vesicle.h              |    1 
 src/tape                   |   13 +
 src/vesicle.c              |   12 +
 src/bondflip.c             |    7 
 src/cluster.c              |  166 +++++++++++++++
 src/initial_distribution.c |   29 ++
 src/general.h              |   21 +
 src/cluster.h              |    9 
 src/tspoststat.c           |  156 ++++++++++++++
 src/vertexmove.c           |    3 
 src/io.h                   |    2 
 src/Makefile.am            |    9 
 src/timestep.c             |    2 
 src/initial_distribution.h |    2 
 src/energy.c               |   49 ++++
 src/restore.h              |    1 
 src/restore.c              |   49 ++++
 src/energy.h               |    4 
 src/tsmeasure.c            |    1 
 23 files changed, 620 insertions(+), 17 deletions(-)

diff --git a/src/Makefile.am b/src/Makefile.am
index fb2c46b..19449c4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,5 @@
-bin_PROGRAMS = trisurf tsmeasure
-trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c
+bin_PROGRAMS = trisurf tsmeasure tspoststat
+trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
 GITVERSION:=$(shell git --no-pager describe --tags --always --dirty)
 AM_CFLAGS = -Wall -Werror -DTS_VERSION=\"$(GITVERSION)\" -fgnu89-inline
 AM_CPPFLAGS = ${libxml2_CFLAGS} -fgnu89-inline
@@ -18,7 +18,10 @@
 #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
 
 
-tsmeasure_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c tsmeasure.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c
+tsmeasure_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c tsmeasure.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
 tsmeasure_LDADD = ${libcurl_LIBS} ${libxml2_LIBS}
+
+tspoststat_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c dumpstate.c frame.c energy.c timestep.c vertexmove.c bondflip.c tspoststat.c poly.c stats.c sh.c shcomplex.c constvol.c snapshot.c restore.c cluster.c
+tspoststat_LDADD = ${libcurl_LIBS} ${libxml2_LIBS}
 #gitversion.c: .git/HEAD .git/index
 #    echo "const char *gitversion = \"$(shell git rev-parse HEAD)\";" > $@
diff --git a/src/bondflip.c b/src/bondflip.c
index 858dc74..66c40f7 100644
--- a/src/bondflip.c
+++ b/src/bondflip.c
@@ -168,6 +168,7 @@
   oldenergy+=kp->xk* kp->energy;
   oldenergy+=km->xk* km->energy;
   oldenergy+=it->xk* it->energy;
+  oldenergy+=bond->energy; /* attraction with neighboring vertices, that have spontaneous curvature */
   //Neigbours of k, it, km, kp don't change its energy.
 
 	if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){dvol = -lm->volume - lp->volume;}
@@ -177,7 +178,7 @@
 */
 
 /* fix data structure for flipped bond */
-    ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1);
+    ts_flip_bond(k,it,km,kp, bond,lm, lp, lm2, lp1, vesicle->tape->w);
 
 
 /* Calculating the new energy */
@@ -186,6 +187,7 @@
   delta_energy+=kp->xk* kp->energy;
   delta_energy+=km->xk* km->energy;
   delta_energy+=it->xk* it->energy;
+  delta_energy+=bond->energy; /* attraction with neighboring vertices, that have spontaneous curvature */
   //Neigbours of k, it, km, kp don't change its energy.
 
     delta_energy-=oldenergy;
@@ -369,7 +371,7 @@
 
 
 ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp,
-ts_bond *bond, ts_triangle *lm, ts_triangle *lp, ts_triangle *lm2, ts_triangle *lp1){
+ts_bond *bond, ts_triangle *lm, ts_triangle *lp, ts_triangle *lm2, ts_triangle *lp1, ts_double w_energy){
 
     ts_uint i; //lmidx, lpidx;
 if(k==NULL || it==NULL || km==NULL || kp==NULL){
@@ -435,6 +437,7 @@
   energy_vertex(kp);
   energy_vertex(km);
   energy_vertex(it);
+  attraction_bond_energy(bond, w_energy);
 // END modifications to data structure!
     return TS_SUCCESS;
 }
diff --git a/src/bondflip.h b/src/bondflip.h
index e3bb19d..a0e0f08 100644
--- a/src/bondflip.h
+++ b/src/bondflip.h
@@ -4,5 +4,5 @@
 
 ts_bool single_bondflip_timestep(ts_vesicle *vesicle, ts_bond *bond, ts_double *rn);
 
-ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp, ts_bond *bond, ts_triangle *lm, ts_triangle *lp, ts_triangle *lm2, ts_triangle *lp1);
+ts_bool ts_flip_bond(ts_vertex *k,ts_vertex *it,ts_vertex *km, ts_vertex *kp, ts_bond *bond, ts_triangle *lm, ts_triangle *lp, ts_triangle *lm2, ts_triangle *lp1, ts_double w_energy);
 #endif
diff --git a/src/cluster.c b/src/cluster.c
new file mode 100644
index 0000000..bb431f2
--- /dev/null
+++ b/src/cluster.c
@@ -0,0 +1,166 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
+#include<stdlib.h>
+#include "general.h"
+#include "cluster.h"
+#include <math.h>
+
+
+
+
+ts_cluster_list *init_cluster_list(){
+	ts_cluster_list *cstlist=(ts_cluster_list *)malloc(sizeof(ts_cluster_list));
+	cstlist->n=0;
+	cstlist->cluster=NULL;
+	return cstlist;
+}
+
+ts_cluster *new_cluster(ts_cluster_list *cstlist){
+	
+	cstlist->n++;
+	cstlist->cluster=(ts_cluster **)realloc(cstlist->cluster,cstlist->n*sizeof(ts_cluster *));
+	if(cstlist->cluster==NULL) fatal("Cannot reallocate memory for additional **ts_cluster.",100);
+	cstlist->cluster[cstlist->n-1]=(ts_cluster *)calloc(1,sizeof(ts_cluster));
+	if(cstlist->cluster[cstlist->n-1]==NULL) fatal("Cannot allocate memory for additional *ts_cluster.",100);
+	return cstlist->cluster[cstlist->n-1];
+}
+
+ts_bool cluster_add_vertex(ts_cluster *cluster, ts_vertex *vtx){
+	cluster->nvtx++;
+	cluster->vtx=(ts_vertex **)realloc(cluster->vtx, cluster->nvtx*sizeof(ts_vertex *));
+	cluster->vtx[cluster->nvtx-1]=vtx;
+	vtx->cluster=cluster;
+	return TS_SUCCESS;
+}
+
+ts_bool cluster_list_compact(ts_cluster_list *cstlist){
+
+
+	ts_uint i,n=cstlist->n;
+	
+	for(i=0;i<cstlist->n;i++){
+		if(cstlist->cluster[i]==NULL){
+			if(i>=n) break;
+			//printf("Have to do compacting n=%d\n ",n);
+			do{
+				n--;
+			} while(cstlist->cluster[n]==NULL && n>i);
+			cstlist->cluster[i]=cstlist->cluster[n];
+			cstlist->cluster[n]=NULL;
+			//printf("After compacting n=%d\n",n);
+		}
+	}
+	cstlist->cluster=(ts_cluster **)realloc(cstlist->cluster,n*sizeof(ts_cluster *));
+	cstlist->n=n;
+	return TS_SUCCESS;
+}
+
+
+ts_bool cluster_join(ts_cluster_list *cstlist, ts_cluster *cluster1, ts_cluster *cluster2){
+	ts_cluster *master_cluster,*slave_cluster;
+	ts_uint i;
+	if(cluster1->idx<cluster2->idx){
+		master_cluster=cluster1;
+		slave_cluster=cluster2;
+	} else {
+		master_cluster=cluster2;
+		slave_cluster=cluster1;
+	}
+	for(i=0;i<slave_cluster->nvtx;i++){
+		cluster_add_vertex(master_cluster,slave_cluster->vtx[i]);
+	}
+	//find cluster in the list and free the location of the cluster
+	for(i=0;i<cstlist->n;i++){
+		if(cstlist->cluster[i]==slave_cluster){
+			cluster_free(slave_cluster);
+			cstlist->cluster[i]=NULL;
+		}
+	}
+	return TS_SUCCESS;
+}
+
+ts_bool cluster_free(ts_cluster *cluster){
+	if(cluster!=NULL){
+		if(cluster->vtx!=NULL)
+			free(cluster->vtx);
+		free(cluster);
+	}
+	return TS_SUCCESS;
+}
+
+ts_bool cluster_list_free(ts_cluster_list *cstlist){
+	ts_uint i;
+	if(cstlist!=NULL){
+		for(i=0;i<cstlist->n;i++){
+			cluster_free(cstlist->cluster[i]);
+		}
+		free(cstlist->cluster);
+		free(cstlist);
+	}
+	return TS_SUCCESS;
+}
+
+
+/* This is a stub function. User should check whether the vertex is clustering or not. */
+ts_bool is_clusterable(ts_vertex *vtx){
+	if(fabs(vtx->c)>1e-6){
+		return 1;
+	}
+	else {
+		return 0;
+	}
+}
+
+
+ts_cluster *cluster_vertex_neighbor(ts_vertex *vtx){
+	int j;
+	for(j=0;j<vtx->neigh_no;j++){
+		if(vtx->neigh[j]->cluster!=NULL)
+			return vtx->neigh[j]->cluster;
+	}
+	return NULL;
+}
+
+ts_bool cluster_vertex_neighbor_check(ts_cluster_list *cstlist, ts_vertex *vtx){
+
+	int j;
+	for(j=0;j<vtx->neigh_no;j++){
+		if(vtx->neigh[j]->cluster!=NULL){
+			if(vtx->neigh[j]->cluster!=vtx->cluster){
+				cluster_join(cstlist, vtx->cluster, vtx->neigh[j]->cluster);
+			}
+		}
+	}
+	return TS_SUCCESS;
+}
+
+
+ts_bool clusterize_vesicle(ts_vesicle *vesicle, ts_cluster_list *cstlist){
+
+	int i;
+	ts_vertex *vtx;
+	ts_cluster *cst;
+	for(i=0;i<vesicle->vlist->n;i++){
+	//for each vertex
+		vtx=vesicle->vlist->vtx[i];
+		if(is_clusterable(vtx)){
+			if(vtx->cluster==NULL){
+				//find first neigbor with cluster index
+				cst=cluster_vertex_neighbor(vtx);
+				if(cst==NULL){
+					//no clusters are around us, vo we are probably lonely vertex or no surronding vertex has been mapped yet.
+					cst=new_cluster(cstlist);
+					cluster_add_vertex(cst,vtx);
+				} else {
+					//we are added to the first cluster found
+					cluster_add_vertex(cst,vtx);
+					cluster_vertex_neighbor_check(cstlist, vtx);
+					cluster_list_compact(cstlist);
+				}
+
+			}
+		}
+	}
+
+
+	return TS_SUCCESS;
+}
diff --git a/src/cluster.h b/src/cluster.h
new file mode 100644
index 0000000..327fa02
--- /dev/null
+++ b/src/cluster.h
@@ -0,0 +1,9 @@
+#ifndef _H_CLUSTER
+#define _H_CLUSTER
+ts_cluster_list *init_cluster_list();
+ts_cluster *new_cluster(ts_cluster_list *cstlist);
+ts_bool cluster_add_vertex(ts_cluster *cluster, ts_vertex *vtx);
+ts_bool cluster_free(ts_cluster *cluster);
+ts_bool cluster_list_free(ts_cluster_list *cstlist);
+ts_bool clusterize_vesicle(ts_vesicle *vesicle, ts_cluster_list *cstlist);
+#endif
diff --git a/src/energy.c b/src/energy.c
index 4f2b386..996fb16 100644
--- a/src/energy.c
+++ b/src/energy.c
@@ -185,3 +185,52 @@
 
     return TS_SUCCESS;
 }
+
+
+
+ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){
+	int i;
+	for(i=0;i<vesicle->blist->n;i++){
+		attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w);
+	}
+	return TS_SUCCESS;
+}
+
+
+inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){
+
+	if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){
+		bond->energy=-w;
+	}
+	else {
+		bond->energy=0.0;
+	}
+	return TS_SUCCESS;
+}
+
+ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){
+	if(fabs(vtx->c)<1e-15) return 0.0;
+//	printf("was here");
+	if(fabs(vesicle->tape->F)<1e-15) return 0.0;
+
+	ts_double norml,ddp=0.0;
+	ts_uint i;
+	ts_double xnorm=0.0,ynorm=0.0,znorm=0.0;
+	/*find normal of the vertex as sum of all the normals of the triangles surrounding it. */
+	for(i=0;i<vtx->tristar_no;i++){
+			xnorm+=vtx->tristar[i]->xnorm;
+			ynorm+=vtx->tristar[i]->ynorm;
+			znorm+=vtx->tristar[i]->znorm;
+	}
+	/*normalize*/
+	norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
+	xnorm/=norml;
+	ynorm/=norml;
+	znorm/=norml;
+	/*calculate ddp, perpendicular displacement*/
+	ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z);
+	/*calculate dE*/
+//	printf("ddp=%e",ddp);
+	return vesicle->tape->F*ddp;		
+	
+}
diff --git a/src/energy.h b/src/energy.h
index e463815..1b97fa9 100644
--- a/src/energy.h
+++ b/src/energy.h
@@ -4,4 +4,8 @@
 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);
+
+ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle);
+inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w);
+ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old);
 #endif
diff --git a/src/general.h b/src/general.h
index 4eab674..983361b 100644
--- a/src/general.h
+++ b/src/general.h
@@ -156,7 +156,8 @@
         ts_double projArea;
         ts_double relR;
         ts_double solAngle;
-	struct ts_poly *grafted_poly;
+		struct ts_poly *grafted_poly;
+		struct ts_cluster *cluster;
 };
 typedef struct ts_vertex ts_vertex;
 
@@ -284,6 +285,10 @@
 	long int mcsweeps;
 	long int quiet;
 	long int shc;
+	long int number_of_vertices_with_c0;
+	ts_double c0;
+	ts_double w;
+	ts_double F;
 } ts_tape;
 
 
@@ -320,6 +325,20 @@
 
 
 
+struct ts_cluster{
+	ts_uint nvtx;
+	ts_uint idx;
+	ts_vertex **vtx;
+};
+
+typedef struct ts_cluster ts_cluster;
+
+typedef struct{
+	ts_uint n;
+	ts_cluster **cluster;
+} ts_cluster_list;
+
+
 /* GLOBAL VARIABLES */
 
 int quiet;
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 71a2c4f..143d9d3 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -44,6 +44,7 @@
 	vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize);
     	vesicle->tape=tape;
 	set_vesicle_values_from_tape(vesicle);
+		initial_population_with_c0(vesicle,tape);
 	return vesicle;
 }
 
@@ -110,12 +111,38 @@
     else {
         vesicle->sphHarmonics=NULL;
     }
+
+    
     return TS_SUCCESS;
 
 }
 
 
-
+ts_bool initial_population_with_c0(ts_vesicle *vesicle, ts_tape *tape){
+	int rndvtx,i,j;
+	if(tape->number_of_vertices_with_c0>0){
+		ts_fprintf(stderr,"Setting values for spontaneous curvature as defined in tape\n");
+		j=0;
+		for(i=0;i<tape->number_of_vertices_with_c0;i++){
+			rndvtx=rand() % vesicle->vlist->n;
+			if(fabs(vesicle->vlist->vtx[rndvtx]->c-tape->c0)<1e-15){
+				j++;
+				i--;
+				if(j>10*vesicle->vlist->n){
+					fatal("cannot populate vesicle with vertices with spontaneous curvature. Too many spontaneous curvature vertices?",100);
+				}
+				continue;
+			}
+			vesicle->vlist->vtx[rndvtx]->c=tape->c0;
+		}
+		mean_curvature_and_energy(vesicle);
+		if(fabs(tape->w)>1e-16){ //if nonzero energy
+			ts_fprintf(stderr,"Setting attraction between vertices with spontaneous curvature\n");
+			sweep_attraction_bond_energy(vesicle);
+		}
+	}
+	return TS_SUCCESS;
+}
 
 
 ts_bool pentagonal_dipyramid_vertex_distribution(ts_vertex_list *vlist){
diff --git a/src/initial_distribution.h b/src/initial_distribution.h
index b1a20e1..e54ee08 100644
--- a/src/initial_distribution.h
+++ b/src/initial_distribution.h
@@ -15,6 +15,8 @@
 
 ts_vesicle *create_vesicle_from_tape(ts_tape *tape);
 ts_bool set_vesicle_values_from_tape(ts_vesicle *vesicle);
+
+ts_bool initial_population_with_c0(ts_vesicle *vesicle, ts_tape *tape);
 /** Sets the initial position of the vertexes to dipyramid
  *
  *      @param *vlist is a pointer to list of vertices
diff --git a/src/io.c b/src/io.c
index a478ca0..21c8f1d 100644
--- a/src/io.c
+++ b/src/io.c
@@ -820,7 +820,7 @@
 	return TS_SUCCESS;
 }
 
-ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno){
+ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno, ts_cluster_list *cstlist){
 	ts_vertex_list *vlist=vesicle->vlist;
 	ts_bond_list *blist=vesicle->blist;
 	ts_vertex **vtx=vlist->vtx;
@@ -862,7 +862,7 @@
 	xml_trisurf_data(fh,vesicle);
 	fprintf(fh, " <UnstructuredGrid>\n");
     fprintf(fh, "<Piece NumberOfPoints=\"%u\" NumberOfCells=\"%u\">\n",vlist->n+monono*polyno+fonono*filno, blist->n+monono*polyno+filno*(fonono-1)+vesicle->tlist->n);
-    fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
+    fprintf(fh,"<PointData Scalars=\"vertices_idx\">\n<DataArray type=\"Int64\" Name=\"vertices_idx\" format=\"ascii\">");
    	for(i=0;i<vlist->n;i++){
 		fprintf(fh,"%u ",vtx[i]->idx);
     }
@@ -887,6 +887,93 @@
 	}
 
     	fprintf(fh,"</DataArray>\n");
+	if(cstlist!=NULL){
+		fprintf(fh,"<DataArray type=\"Int64\" Name=\"vertices_in_cluster\" format=\"ascii\">");
+		for(i=0;i<vlist->n;i++){
+			if(vtx[i]->cluster!=NULL){
+				fprintf(fh,"%u ",vtx[i]->cluster->nvtx);
+			} else {
+				fprintf(fh,"-1 ");
+			}
+	    	}
+		//polymeres
+		if(poly){
+			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,"-1 ");
+				}
+			}
+		}
+		//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,"-1 ");
+				}
+			}
+		}
+
+		fprintf(fh,"</DataArray>\n");
+
+
+	}
+
+	//here comes additional data as needed. Currently only spontaneous curvature
+	fprintf(fh,"<DataArray type=\"Float64\" Name=\"spontaneous_curvature\" format=\"ascii\">");
+	for(i=0;i<vlist->n;i++){
+		fprintf(fh,"%.17e ",vtx[i]->c);
+	}
+		//polymeres
+		if(poly){
+			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,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->c);
+				}
+			}
+		}
+		//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,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->c);
+				}
+			}
+		}
+    fprintf(fh,"</DataArray>\n");
+
+	//here comes additional data. Energy!
+	fprintf(fh,"<DataArray type=\"Float64\" Name=\"bending_energy\" format=\"ascii\">");
+	for(i=0;i<vlist->n;i++){
+		fprintf(fh,"%.17e ",vtx[i]->energy*vtx[i]->xk);
+	}
+		//polymeres
+		if(poly){
+			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,"%.17e ", vesicle->poly_list->poly[i]->vlist->vtx[j]->energy* vesicle->poly_list->poly[i]->k);
+				}
+			}
+		}
+		//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,"%.17e ",  vesicle->filament_list->poly[i]->vlist->vtx[j]->energy*  vesicle->filament_list->poly[i]->k);
+				}
+			}
+		}
+    fprintf(fh,"</DataArray>\n");
+
+
 	
 	fprintf(fh,"</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++){
@@ -1093,6 +1180,10 @@
         CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
         CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
 	CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
+	CFG_SIMPLE_INT("number_of_vertices_with_c0", &tape->number_of_vertices_with_c0),
+	CFG_SIMPLE_FLOAT("c0",&tape->c0),
+	CFG_SIMPLE_FLOAT("w",&tape->w),
+	CFG_SIMPLE_FLOAT("F",&tape->F),
         CFG_END()
     };
     cfg_t *cfg;    
diff --git a/src/io.h b/src/io.h
index 2243555..6ca8331 100644
--- a/src/io.h
+++ b/src/io.h
@@ -72,7 +72,7 @@
  *	@param *text is a description line (max. 255 characters) to be included in the file
  */
 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_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno, ts_cluster_list *cstlist);
 ts_bool write_master_xml_file(ts_char *filename);
 ts_bool write_pov_file(ts_vesicle *vesicle, char *filename);
 
diff --git a/src/main.c b/src/main.c
index cb6dd05..3f91f43 100644
--- a/src/main.c
+++ b/src/main.c
@@ -129,6 +129,8 @@
 			//printf("nucleus coords: %.17e %.17e %.17e\n",vesicle->nucleus_center[0], vesicle->nucleus_center[1], vesicle->nucleus_center[2]);
 //	write_vertex_xml_file(vesicle,0);
 //	exit(1);
+
+			//write_vertex_xml_file(vesicle,1000);
 	run_simulation(vesicle, tape->mcsweeps, tape->inititer, tape->iterations, start_iteration);
 	write_master_xml_file(command_line_args.output_fullfilename);
 	write_dout_fcompat_file(vesicle,"dout");
diff --git a/src/restore.c b/src/restore.c
index 086bf1c..c89ff6d 100644
--- a/src/restore.c
+++ b/src/restore.c
@@ -54,6 +54,10 @@
 				if ((!xmlStrcmp(cur1->name, (const xmlChar *)"Piece"))){
 					cur2=cur1->xmlChildrenNode;
 					while(cur2!=NULL){
+						if ((!xmlStrcmp(cur2->name, (const xmlChar *)"PointData"))){
+							if(vesicle!=NULL)
+								parseXMLPointData(vesicle,doc,cur2);
+						}
 						if ((!xmlStrcmp(cur2->name, (const xmlChar *)"Points"))){
 							//fprintf(stderr,"Found point data\n");
 							if(vesicle!=NULL)
@@ -81,6 +85,7 @@
 
 	init_normal_vectors(vesicle->tlist);
 	mean_curvature_and_energy(vesicle);
+	sweep_attraction_bond_energy(vesicle);
 
 /* TODO: filaments */
 
@@ -182,7 +187,6 @@
 
 	vesicle->tape=tape;
 	set_vesicle_values_from_tape(vesicle);
-
 	return vesicle;
 }
 
@@ -318,7 +322,48 @@
 	return TS_SUCCESS;
 }
 
-
+/* this parses the data for vertices (like spontaneous curvature, etc.) */
+ts_bool parseXMLPointData(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
+	xmlNodePtr child = cur->xmlChildrenNode;
+	xmlChar *property_name;
+	xmlChar *values;
+	char *vals;
+	char *token;
+	int idx, polyidx, monoidx, filidx, fonoidx;
+	while (child != NULL) {
+		if ((!xmlStrcmp(child->name, (const xmlChar *)"DataArray"))){
+			property_name=xmlGetProp(child, (xmlChar *)"Name");
+		//	fprintf(stderr,"Name: %s\n", property_name);
+			if(!xmlStrcmp(property_name,(const xmlChar *)"spontaneous_curvature")){
+				values=xmlNodeListGetString(doc,child->xmlChildrenNode,1);
+				vals=(char *)values;
+				token=strtok(vals," ");
+				idx=0;
+				while(token!=NULL){
+					if(idx<vesicle->vlist->n){
+						vesicle->vlist->vtx[idx]->c=atof(token);
+					} else if(vesicle->tape->nmono && vesicle->tape->npoly && idx<vesicle->vlist->n+vesicle->tape->nmono*vesicle->tape->npoly) {
+						polyidx=(idx-vesicle->vlist->n)/vesicle->tape->nmono;
+						monoidx=(idx-vesicle->vlist->n)%vesicle->tape->nmono;
+						vesicle->poly_list->poly[polyidx]->vlist->vtx[monoidx]->c=atof(token);
+					} else {
+						filidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)/vesicle->tape->nfono;
+						fonoidx=(idx-vesicle->vlist->n-vesicle->tape->nmono*vesicle->tape->npoly)%vesicle->tape->nfono;
+						//fprintf(stderr,"filidx=%d, fonoidx=%d, coord=%s,%s,%s\n",filidx,fonoidx,token[0],token[1],token[2]);
+						vesicle->filament_list->poly[filidx]->vlist->vtx[fonoidx]->c=atof(token);
+					}
+					idx++;
+					token=strtok(NULL," ");
+				}
+				xmlFree(values);		
+			}
+			xmlFree(property_name);
+		}
+		
+		child=child->next;
+	}
+	return TS_SUCCESS;
+}
 /* this is a parser of vertex positions and bonds from main xml data */
 ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur){
 	xmlNodePtr child = cur->xmlChildrenNode;
diff --git a/src/restore.h b/src/restore.h
index cfc6b66..b05c6a8 100644
--- a/src/restore.h
+++ b/src/restore.h
@@ -9,6 +9,7 @@
 ts_bool parseTrisurfTria(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
 ts_bool parseTrisurfTriaNeigh(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
 ts_bool parseTrisurfTristar(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
+ts_bool parseXMLPointData(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
 ts_bool parseXMLVertexPosition(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
 ts_bool parseXMLBonds(ts_vesicle *vesicle,xmlDocPtr doc, xmlNodePtr cur);
 ts_bool parseTrisurfNucleus(ts_vesicle *vesicle, xmlDocPtr doc, xmlNodePtr cur);
diff --git a/src/snapshot.c b/src/snapshot.c
index fbb6d4b..da21498 100644
--- a/src/snapshot.c
+++ b/src/snapshot.c
@@ -24,7 +24,7 @@
 ts_bool xml_trisurf_data(FILE *fh, ts_vesicle *vesicle){
 
 	ts_string *data=(ts_string *)malloc(sizeof(ts_sprintf));
-	data->string=(char *)malloc(512000*sizeof(char)); /*TODO: warning, can break if the string is to long */
+	data->string=(char *)malloc(5120000*sizeof(char)); /*TODO: warning, can break if the string is to long */
 	data->beg=0;
 	
 	xml_trisurf_header(fh, vesicle);
diff --git a/src/tape b/src/tape
index 93e42d7..f2edf02 100644
--- a/src/tape
+++ b/src/tape
@@ -56,7 +56,7 @@
 ####### Program Control ############
 #how many MC sweeps between subsequent records of states to disk
 #200000
-mcsweeps=200
+mcsweeps=2000
 #how many initial mcsweeps*inititer MC sweeps before recording to disk?
 #2
 inititer=1
@@ -84,3 +84,14 @@
 #max number of processes in distributed (voluntary) environment
 distributed_processes=50
 #cuda options???
+
+
+#NirGov branch only!
+#number of vertices with spontaneous curvature (integer)
+number_of_vertices_with_c0=100
+#spontaneous curvature (float)
+c0=0.5
+#energy of attraction of vertices with spontaneous curvature (float, positive value for attraction)
+w=10.0
+#direct force on vesicles with spontaneous curvature (float)
+F=2.0
diff --git a/src/timestep.c b/src/timestep.c
index be14619..f7cb5dd 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -93,7 +93,7 @@
 		cell_occupation(vesicle);
             dump_state(vesicle,i);
 		if(i>=inititer){
-			write_vertex_xml_file(vesicle,i-inititer);
+			write_vertex_xml_file(vesicle,i-inititer,NULL);
 			write_master_xml_file(command_line_args.output_fullfilename);
 			epochtime=get_epoch();			
 			gyration_eigen(vesicle, &l1, &l2, &l3);
diff --git a/src/tsmeasure.c b/src/tsmeasure.c
index 4ed39eb..beabc61 100644
--- a/src/tsmeasure.c
+++ b/src/tsmeasure.c
@@ -84,7 +84,6 @@
 			tape_free(vesicle->tape);
 			vesicle_free(vesicle);
             	}
-		free(ent);  
 		}
 	for (n = 0; n < count; n++)
   	{
diff --git a/src/tspoststat.c b/src/tspoststat.c
new file mode 100644
index 0000000..f518000
--- /dev/null
+++ b/src/tspoststat.c
@@ -0,0 +1,156 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include "general.h"
+//#include "vertex.h"
+//#include "bond.h"
+//#include "triangle.h"
+//#include "cell.h"
+#include "vesicle.h"
+#include "io.h"
+//#include "initial_distribution.h"
+//#include "frame.h"
+//#include "timestep.h"
+//#include "poly.h"
+#include "stats.h"
+#include "sh.h"
+#include "shcomplex.h"
+#include "dumpstate.h"
+#include "restore.h"
+#include "cluster.h"
+#include <string.h>
+#include <getopt.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <dirent.h>
+#include <errno.h>
+#include <snapshot.h>
+#include<gsl/gsl_complex.h>
+#include<gsl/gsl_complex_math.h>
+#include<stdio.h>
+
+ts_vesicle *restoreVesicle(char *filename){
+	ts_vesicle *vesicle = parseDump(filename);
+	return vesicle;
+}
+
+void vesicle_calculate_ulm2(ts_vesicle *vesicle){
+	//complex_sph_free(vesicle->sphHarmonics);
+
+	//vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,21);
+	vesicle_volume(vesicle);
+	preparationSh(vesicle,getR0(vesicle));
+	calculateUlmComplex(vesicle);
+	ts_int i,j;
+	for(i=0;i<vesicle->sphHarmonics->l;i++){
+    		for(j=i;j<2*i+1;j++){
+			printf("%e ", gsl_complex_abs2(vesicle->sphHarmonics->ulmComplex[i][j]));
+    		}
+	}
+		printf("\n");
+
+}
+
+
+
+int count_bonds_with_energy(ts_bond_list *blist){
+
+	unsigned int i, cnt;
+	cnt=0;
+	for(i=0;i<blist->n;i++){
+		if(fabs(blist->bond[i]->energy)>1e-16) cnt++;
+	}
+	return cnt;
+}
+
+
+
+ts_bool write_histogram_data(ts_uint timestep_no, ts_vesicle *vesicle){
+	ts_cluster_list *cstlist=init_cluster_list();
+	clusterize_vesicle(vesicle,cstlist);
+	//printf("No clusters=%d\n",cstlist->n);
+	int k,i,cnt, test=0;
+	int max_nvtx=0;
+	char filename[255];
+	sprintf(filename,"histogram_%.6u.csv",timestep_no);
+	FILE *fd=fopen(filename,"w");
+	fprintf(fd,"Number_of_vertices_in cluster Number_of_clusters\n");
+	for(k=0;k<cstlist->n;k++)
+		if(cstlist->cluster[k]->nvtx>max_nvtx) max_nvtx=cstlist->cluster[k]->nvtx;
+	//printf("Max. number of vertices in cluster: %d\n",max_nvtx);
+	for(i=1;i<=max_nvtx;i++){
+		cnt=0;
+		for(k=0;k<cstlist->n;k++)
+			if(cstlist->cluster[k]->nvtx==i) cnt++;
+		fprintf(fd,"%d %d\n",i,cnt);
+		test+=cnt*i;
+	}
+	//for(k=0;k<cstlist->n;k++){
+//		printf("*Cluster %d has %d vertices\n",k,cstlist->cluster[k]->nvtx);
+//	}
+
+	fclose(fd);
+//	printf("*Sum of all vertices in clusters: %d\n", test);
+//	write_vertex_xml_file(vesicle,timestep_no,cstlist);
+	cluster_list_free(cstlist);
+	
+	return TS_SUCCESS;
+}
+
+
+int main(){
+	ts_vesicle *vesicle;
+	ts_char *i,*j;
+	ts_uint tstep,n;
+    	ts_char *number;
+	struct dirent **list;
+	ts_double l1,l2,l3,hbar;
+	int count;
+	ts_fprintf(stderr,"TRISURF-NG v. %s, compiled on: %s %s.\n", TS_VERSION, __DATE__, __TIME__);
+
+	fprintf(stdout, "OuterLoop Volume Area lamdba1 lambda2 lambda3 Nbw/Nb hbar\n");
+
+
+	count=scandir(".",&list,0,alphasort);
+	if(count<0){
+		fatal("Error, cannot open directory.",1);
+	}
+        tstep=0;
+	for(n=0;n<count;n++){
+		struct dirent *ent;
+		ent=list[n];	
+            	i=rindex(ent->d_name,'.');
+            	if(i==NULL) {
+				continue;
+		}
+            	if(strcmp(i+1,"vtu")==0){
+                    j=rindex(ent->d_name,'_');
+                    if(j==NULL) continue;
+                    number=strndup(j+1,j-i); 
+			quiet=1;
+                    ts_fprintf(stdout,"timestep: %u filename: %s\n",atoi(number),ent->d_name);
+//			printf("%u ",atoi(number));
+			vesicle=restoreVesicle(ent->d_name);
+//			vesicle_calculate_ulm2(vesicle);
+			vesicle_volume(vesicle);
+			vesicle_area(vesicle);
+			gyration_eigen(vesicle,&l1,&l2,&l3);
+			hbar=vesicle_meancurvature(vesicle)/vesicle->area;			
+			fprintf(stdout,"%d %.17e %.17e %.17e %.17e %.17e %.17e %.17e\n",atoi(number),vesicle->volume, vesicle->area,l1,l2,l3, (ts_double)count_bonds_with_energy(vesicle->blist)/(ts_double)vesicle->blist->n,hbar);
+                    	tstep++;
+			write_histogram_data(atoi(number), vesicle);
+                    free(number);
+			tape_free(vesicle->tape);
+			vesicle_free(vesicle);
+            	}
+		}
+	for (n = 0; n < count; n++)
+  	{
+  		free(list[n]);
+  	}
+	
+	free(list);
+	return 0;
+}
+
diff --git a/src/vertexmove.c b/src/vertexmove.c
index fb74a06..2517f69 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -189,6 +189,8 @@
     delta_energy+=delta_energy_cv;
 //    fprintf(stderr,"Denergy after=%e\n",delta_energy);
     }
+/* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */
+	delta_energy+=direct_force_energy(vesicle,vtx,backupvtx);
 /* No poly-bond energy for now!
 	if(vtx->grafted_poly!=NULL){
 		delta_energy+=
@@ -198,6 +200,7 @@
 */
 //   fprintf(stderr, "DE=%f\n",delta_energy);
     //MONTE CARLOOOOOOOO
+//	if(vtx->c!=0.0) printf("DE=%f\n",delta_energy);
     if(delta_energy>=0){
 #ifdef TS_DOUBLE_DOUBLE
         if(exp(-delta_energy)< drand48())
diff --git a/src/vesicle.c b/src/vesicle.c
index bc69e62..fe4c4b1 100644
--- a/src/vesicle.c
+++ b/src/vesicle.c
@@ -79,3 +79,15 @@
     vesicle->area=area;
     return TS_SUCCESS;
 }
+
+ts_double vesicle_meancurvature(ts_vesicle *vesicle){
+// Integrates (H dA) over vesicle area A, where H=(C1+C2)/2.
+// (To be devided by A outside of function)
+	ts_double mc;
+	ts_uint i;
+	mc=0;
+	for(i=0;i<vesicle->vlist->n;i++){
+		mc=mc+vesicle->vlist->vtx[i]->curvature;
+	}
+	return mc/2.0;
+}
diff --git a/src/vesicle.h b/src/vesicle.h
index e8f09fb..953e0a2 100644
--- a/src/vesicle.h
+++ b/src/vesicle.h
@@ -7,4 +7,5 @@
 ts_bool vesicle_free(ts_vesicle *vesicle);
 ts_bool vesicle_volume(ts_vesicle *vesicle);
 ts_bool vesicle_area(ts_vesicle *vesicle);
+ts_double vesicle_meancurvature(ts_vesicle *vesicle);
 #endif

--
Gitblit v1.9.3