From 6487a03362d35589c286d0660eba123a7ab00535 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Mon, 10 May 2021 11:14:04 +0000
Subject: [PATCH] Raj and Yoav session in making multiple protein species
---
src/io.c | 201 +++++++++++++++++++++++++++++++++++++++++++-------
1 files changed, 172 insertions(+), 29 deletions(-)
diff --git a/src/io.c b/src/io.c
index f3029db..034207b 100644
--- a/src/io.c
+++ b/src/io.c
@@ -1,3 +1,4 @@
+
/* vim: set ts=4 sts=4 sw=4 noet : */
#include "general.h"
#include<stdio.h>
@@ -423,7 +424,8 @@
vesicle->tape=parsetape(command_line_args.tape_fullfilename);
// recreating space for cells //
vesicle->clist=init_cell_list(vesicle->tape->ncxmax, vesicle->tape->ncymax, vesicle->tape->nczmax, vesicle->tape->stepsize);
- vesicle->clist->max_occupancy=8;
+//THIS IS HARDCODED IN CELL.C NOW
+// vesicle->clist->max_occupancy=16;
// vesicle->tape=(ts_tape *)malloc(sizeof(ts_tape));
// retval=fread(vesicle->tape, sizeof(ts_tape),1,fh);
retval=fread(iteration,sizeof(ts_uint),1,fh);
@@ -749,22 +751,6 @@
}
-ts_bool write_dout_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
- FILE *fh;
- fh=fopen(filename, "w");
- if(fh==NULL){
- err("Cannot open file %s for writing");
- return TS_FAIL;
- }
- fprintf(fh,"%.17E\n%.17E\n",vesicle->stepsize,vesicle->dmax);
- fprint_vertex_list(fh,vesicle->vlist);
- fprint_tristar(fh,vesicle);
- fprint_triangle_list(fh,vesicle);
- fprint_vertex_data(fh,vesicle->vlist);
- fprint_bonds(fh,vesicle);
- fclose(fh);
- return TS_SUCCESS;
-}
ts_bool read_tape_fcompat_file(ts_vesicle *vesicle, ts_char *filename){
FILE *fh;
@@ -820,11 +806,12 @@
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;
ts_uint i,j;
+ //ts_double senergy=0.0;
char filename[10000];
char just_name[255];
FILE *fh;
@@ -843,7 +830,8 @@
ts_uint monono=0, polyno=0, poly_idx=0, filno=0, fonono=0;
ts_bool poly=0, fil=0;
if(vesicle->poly_list!=NULL){
- if(vesicle->poly_list->poly[0]!=NULL){
+ if(vesicle->poly_list->n!=0){
+ //if(vesicle->poly_list->poly[0]!=NULL){
polyno=vesicle->poly_list->n;
monono=vesicle->poly_list->poly[0]->vlist->n;
poly=1;
@@ -851,7 +839,8 @@
}
if(vesicle->filament_list!=NULL){
- if(vesicle->filament_list->poly[0]!=NULL){
+ if(vesicle->filament_list->n!=0){
+ //if(vesicle->filament_list->poly[0]!=NULL){
filno=vesicle->filament_list->n;
fonono=vesicle->filament_list->poly[0]->vlist->n;
fil=1;
@@ -861,8 +850,8 @@
fprintf(fh, "<?xml version=\"1.0\"?>\n<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
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));
- fprintf(fh,"<PointData Scalars=\"scalars\">\n<DataArray type=\"Int64\" Name=\"scalars\" format=\"ascii\">");
+ 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=\"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,9 +876,136 @@
}
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");
+ fprintf(fh,"<DataArray type=\"Float64\" Name=\"direct_interaction_force\" format=\"ascii\">");
+ for(i=0;i<vlist->n;i++){
+ fprintf(fh,"%.17e ",vtx[i]->direct_interaction_force);
+ }
+ //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]->direct_interaction_force);
+ }
+ }
+ }
+ //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]->direct_interaction_force);
+ }
+ }
+ }
+ 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");
+ fprintf(fh,"</PointData>\n<CellData>\n");
+
+ if(vesicle->tape->stretchswitch==1){
+ fprintf(fh,"<DataArray type=\"Float64\" Name=\"stretching_energy\" format=\"ascii\">");
+ for(i=0;i<blist->n;i++){
+ fprintf(fh, "0.0 ");
+ }
+ for(i=0;i<monono*polyno+filno*(fonono-1);i++){
+ fprintf(fh,"0.0 ");
+ }
+ for(i=0;i<vesicle->tlist->n;i++){
+ fprintf(fh,"%.17e ",vesicle->tlist->tria[i]->energy);
+ }
+ fprintf(fh,"</DataArray>\n");
+ }
+
+
+
+ fprintf(fh,"</CellData>\n<Points>\n<DataArray type=\"Float64\" Name=\"Koordinate tock\" NumberOfComponents=\"3\" format=\"ascii\">\n");
for(i=0;i<vlist->n;i++){
fprintf(fh,"%.17e %.17e %.17e\n",vtx[i]->x,vtx[i]->y, vtx[i]->z);
}
@@ -940,17 +1056,24 @@
}
}
-
+ for(i=0;i<vesicle->tlist->n;i++){
+ fprintf(fh,"%u %u %u\n", vesicle->tlist->tria[i]->vertex[0]->idx, vesicle->tlist->tria[i]->vertex[1]->idx, vesicle->tlist->tria[i]->vertex[2]->idx);
+ }
fprintf(fh,"</DataArray>\n<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">");
for (i=2;i<(blist->n+monono*polyno+(fonono-1)*filno)*2+1;i+=2){
fprintf(fh,"%u ",i);
}
+ for(j=i+1;j<i+3*(vesicle->tlist->n);j+=3){ //let's continue counting from where we left of
+ fprintf(fh,"%u ", j);
+ }
fprintf(fh,"\n");
fprintf(fh,"</DataArray>\n<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
- for (i=0;i<blist->n+monono*polyno+fonono*filno;i++){
+ for (i=0;i<blist->n+monono*polyno+(fonono-1)*filno;i++){
fprintf(fh,"3 ");
}
-
+ for(i=0;i<vesicle->tlist->n;i++){
+ fprintf(fh,"5 ");
+ }
fprintf(fh,"</DataArray>\n</Cells>\n</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
fclose(fh);
return TS_SUCCESS;
@@ -1059,6 +1182,7 @@
CFG_SIMPLE_INT("nmono", &tape->nmono),
CFG_SIMPLE_INT("nfil",&tape->nfil),
CFG_SIMPLE_INT("nfono",&tape->nfono),
+ CFG_SIMPLE_INT("internal_poly",&tape->internal_poly),
CFG_SIMPLE_INT("R_nucleus",&tape->R_nucleus),
CFG_SIMPLE_FLOAT("R_nucleusX",&tape->R_nucleusX),
CFG_SIMPLE_FLOAT("R_nucleusY",&tape->R_nucleusY),
@@ -1070,6 +1194,8 @@
CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch),
CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch),
CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision),
+ CFG_SIMPLE_INT("stretchswitch",&tape->stretchswitch),
+ CFG_SIMPLE_FLOAT("xkA0",&tape->xkA0),
CFG_SIMPLE_FLOAT("pressure",&tape->pressure),
CFG_SIMPLE_FLOAT("k_spring",&tape->kspring),
CFG_SIMPLE_FLOAT("xi",&tape->xi),
@@ -1080,18 +1206,33 @@
CFG_SIMPLE_INT("iterations",&tape->iterations),
CFG_SIMPLE_INT("mcsweeps",&tape->mcsweeps),
CFG_SIMPLE_INT("inititer", &tape->inititer),
- CFG_SIMPLE_BOOL("quiet",&tape->quiet),
- CFG_SIMPLE_STR("multiprocessing",tape->multiprocessing),
+ CFG_SIMPLE_BOOL("quiet",(cfg_bool_t *)&tape->quiet),
+ CFG_SIMPLE_STR("multiprocessing",&tape->multiprocessing),
CFG_SIMPLE_INT("smp_cores",&tape->brezveze0),
CFG_SIMPLE_INT("cluster_nodes",&tape->brezveze1),
CFG_SIMPLE_INT("distributed_processes",&tape->brezveze2),
CFG_SIMPLE_INT("spherical_harmonics_coefficients",&tape->shc),
+ CFG_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),
+/* Variables related to plane confinement */
+ CFG_INT("plane_confinement_switch", 0, CFGF_NONE),
+ CFG_FLOAT("plane_d", 15, CFGF_NONE),
+ CFG_FLOAT("plane_F", 1000, CFGF_NONE),
+/* Variables related to stretching */
+// CFG_FLOAT("stretchswitch", 0, CFGF_NONE),
+// CFG_FLOAT("xkA0",0,CFGF_NONE),
CFG_END()
};
cfg_t *cfg;
ts_uint retval;
cfg = cfg_init(opts, 0);
retval=cfg_parse_buf(cfg, buffer);
+ tape->plane_confinement_switch=cfg_getint(cfg,"plane_confinement_switch");
+ tape->plane_d=cfg_getfloat(cfg,"plane_d");
+ tape->plane_F=cfg_getfloat(cfg,"plane_F");
+
if(retval==CFG_FILE_ERROR){
fatal("No tape file.",100);
}
@@ -1120,6 +1261,8 @@
ts_bool getcmdline_tape(cfg_t *cfg, char *opts){
char *commands, *backup, *saveptr, *saveopptr, *command, *operator[2];
+ operator[0]=0;
+ operator[1]=0;
ts_uint i,j;
commands=(char *)malloc(10000*sizeof(char));
backup=commands; //since the pointer to commands will be lost, we acquire a pointer that will serve as backup.
--
Gitblit v1.9.3