Merge branch 'nirgov' into HEAD

Samo Penic

2016-12-06 0634b97e436e480b203a22e1113c473d3ae4d281

Merge branch 'nirgov' into HEAD

 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)\";" > $@

 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;
}

 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

 src/cluster.c

New file
@@ -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;
}

 src/cluster.h

New file
@@ -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

 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;		
	
}

 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

 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;

 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){

 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

 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,11 @@
    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);
<<<<<<< HEAD
    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\">");
>>>>>>> nirgov
       for(i=0;i<vlist->n;i++){
        fprintf(fh,"%u ",vtx[i]->idx);
    }
@@ -887,6 +891,96 @@
    }

        fprintf(fh,"</DataArray>\n");
<<<<<<< HEAD
=======
    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");


>>>>>>> nirgov
    
    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 +1187,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;    

 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);


 src/main.c

@@ -128,6 +128,7 @@
    }
            //printf("nucleus coords: %.17e %.17e %.17e\n",vesicle->nucleus_center[0], vesicle->nucleus_center[1], vesicle->nucleus_center[2]);

            //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");

 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;

 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);

 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);

 src/tape

@@ -1,6 +1,10 @@
####### Vesicle definitions ###########
# nshell is a number of divisions of dipyramid
<<<<<<< HEAD
nshell=17
=======
nshell=10
>>>>>>> nirgov
# dmax is the max. bond length (in units l_min)
dmax=1.7
# dmin_interspecies in the min. dist. between different vertex species (in units l_min)
@@ -56,7 +60,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=10
@@ -84,3 +88,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

 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);

 src/tsmeasure.c

@@ -84,7 +84,6 @@
            tape_free(vesicle->tape);
            vesicle_free(vesicle);
                }
        free(ent);  
        }
    for (n = 0; n < count; n++)
      {

 src/tspoststat.c

New file
@@ -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;
}


 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())

 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;
}

 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