Allink  v0.1
VarDataWrite.cpp
1 /***********************************************************************
2 VarData: This Program reads and writes a specific file format
3 storing all the information relative to a set of equal structure
4 polymers to the CHAIN, PART and GENERAL structures. It provides
5 two different ways to backfold the coordinates, a fanction that
6 creates an initial system with different option and some function
7 for the data analisys. The first calculate the distribution of the
8 monomer in the box, the second the distribution of the bonds.
9 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
10 
11 
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version.
16 
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU General Public License for more details.
21 
22 You should have received a copy of the GNU General Public License
23 along with this program; if not, write to the Free Software
24 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 ***********************************************************************/
26 #include "../include/VarData.h"
27 ;
28 #ifdef USE_BOOST
29 void VarData::OpenComprWrite(char *FName,filtering_ostream ZipFile){
30  ofstream OpenZipFile(FName, ios_base::binary);
31  // create filter:
32  ZipFile.push(bzip2_compressor());
33  ZipFile.push(OpenZipFile);
34 }
35 void VarData::OpenComprRead(char *FName,filtering_istream ZipFile){
36  ifstream OpenZipFile(FName, ios_base::binary);
37  // create filter:
38  ZipFile.push(bzip2_compressor());
39  ZipFile.push(OpenZipFile);
40 }
41 void VarData::CloseCompr(filtering_istream ZipFile){
42  ZipFile.reset();
43 }
44 void VarData::CloseCompr(filtering_ostream ZipFile){
45  ZipFile.reset();
46 }
47 #endif
48 bool VarData::Write(char *OutFile){
49  if( SysFormat == VAR_SYS_TXVL ){
50  if(WriteTxvl(OutFile))
51  return 1;
52  }
53  else if( SysFormat == VAR_SYS_XVT ){
54  if(WriteXvt(OutFile))
55  return 1;
56  }
57  else
58  WriteXyz(OutFile);
59  return 0;
60 }
61 bool VarData::WriteTxvl(char *OutFile){
62  VarMessage("Write");
63  FILE *FileToWrite = fopen(OutFile,"w");
64  SigErr(FileToWrite==NULL,"The output file %s could not be opened\n",OutFile);
65  fprintf(FileToWrite,"# ");
66  if( VAR_IF_TYPE(SysType,VAR_EDGE) ){
67  fprintf(FileToWrite,"l( %.2g %.2g %.2g) ",
68  Gen->Edge[0],Gen->Edge[1],Gen->Edge[2]);
69  }
70  fprintf(FileToWrite,"c[%d] ",Gen->NChain);
71  fprintf(FileToWrite,"s[%d] ",Gen->Step);
72  fprintf(FileToWrite,"n[%d] ",Gen->NPart);
73  fprintf(FileToWrite,"L[%d] ",Gen->NLink);
74  fprintf(FileToWrite,"v[%d] ",NEdge);
75  fprintf(FileToWrite,"d[%s] ",cWhat2Draw);
76  fprintf(FileToWrite,"D(%lf) ",Gen->Deltat);
77  fprintf(FileToWrite,"i(%.2f %.2f %.2f) ",Gen->rho,Gen->chiN,Gen->kappaN);
78  fprintf(FileToWrite,"\n");
79  HeaderNano(FileToWrite);
80  for(int p=0;p<Gen->NPart;p++){
81  fprintf(FileToWrite,"{ t[%d %d %d] ",
82  Pm[p].Idx,Pm[p].CId,Pm[p].Typ);
83  fprintf(FileToWrite,"x(%lf %lf %lf) ",
84  Pm[p].Pos[0],Pm[p].Pos[1],Pm[p].Pos[2]);
85  fprintf(FileToWrite,"v(%lf %lf %lf) ",
86  Pm[p].Vel[0],Pm[p].Vel[1],Pm[p].Vel[2]);
87  // printf("Pm[%d].NLink %d\n",p,Pm[p].NLink);
88  if(pNLink() > 0){
89  for(int i=0;i<Ln[p].NLink;i++)
90  fprintf(FileToWrite,"l[%d] ",Ln[p].Link[i]);
91  }
92  fprintf(FileToWrite,"}\n");//End of particle info
93  }
94  // for(int n=0;n<Gen->NNano;n++){
95  // fprintf(FileToWrite,"{t[%d %d %d] ",Gen->NPart-1+n,-1,2);
96  // fprintf(FileToWrite,"x(%lf %lf %lf) ",Nano[n].Pos[0],Nano[n].Pos[1],Nano[n].Pos[2]);
97  // fprintf(FileToWrite,"v(%lf %lf %lf)}\n",0.,0.,0.);
98  // }
99  fclose(FileToWrite);
100  return 0;
101 }
102 void VarData::HeaderInteraction(FILE *FileToWrite){
103  double Norm2 = CUBE(pReOverCutOff()) / SQR(pNPCh());
104  double Norm3 = CUBE(SQR(pReOverCutOff()) / pNPCh());
105  MInt->Rescale(1./Norm2,2);
106  MInt->Rescale(1./Norm3,3);
107  fprintf(FileToWrite,"# v=%.2f %.2f %.2f %.2f %.2f %.2f ",MInt->Coeff(0,0),MInt->Coeff(0,1),MInt->Coeff(0,2),MInt->Coeff(1,1),MInt->Coeff(1,2),MInt->Coeff(2,2));
108  fprintf(FileToWrite,"w=%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",MInt->Coeff(0,0,0),MInt->Coeff(0,0,1),MInt->Coeff(0,0,2),MInt->Coeff(0,1,1),MInt->Coeff(0,1,2),MInt->Coeff(0,2,2),MInt->Coeff(1,1,1),MInt->Coeff(1,1,2),MInt->Coeff(1,2,2),MInt->Coeff(2,2,2));
109  fprintf(FileToWrite,"# a2=%.1f a3=%.1f Re=%.1f N=%d ks=%.3f kb=%.3f l0=0\n",Gen->WFuncStraight2,Gen->WFuncStraight3,Gen->ReOverCutOff,Gen->NPCh,Gen->kappaSpring,Gen->kappaBend);
110  MInt->Rescale(Norm2,2);
111  MInt->Rescale(Norm3,3);
112 }
113 void VarData::ShapeId(int iShape,char *Shape){
114  sprintf(Shape,"no");
115  if(VAR_IF_TYPE(iShape,SHAPE_SPH)) sprintf(Shape,"sph");
116  else if(VAR_IF_TYPE(iShape,SHAPE_TIP)) sprintf(Shape,"tip");
117  else if(VAR_IF_TYPE(iShape,SHAPE_CYL)) sprintf(Shape,"cyl");
118  else if(VAR_IF_TYPE(iShape,SHAPE_PILL)) sprintf(Shape,"pill");
119  else if(VAR_IF_TYPE(iShape,SHAPE_TILT)) sprintf(Shape,"tilt");
120  else if(VAR_IF_TYPE(iShape,SHAPE_WALL)) sprintf(Shape,"wall");
121  else if(VAR_IF_TYPE(iShape,SHAPE_PORE)) sprintf(Shape,"pore");
122  else if(VAR_IF_TYPE(iShape,SHAPE_EXT)) sprintf(Shape,"ext");
123  else if(VAR_IF_TYPE(iShape,SHAPE_JANUS)) sprintf(Shape,"janus");
124  else if(VAR_IF_TYPE(iShape,SHAPE_STALK)) sprintf(Shape,"stalk");
125  else if(VAR_IF_TYPE(iShape,SHAPE_TORUS)) sprintf(Shape,"torus");
126  else if(VAR_IF_TYPE(iShape,SHAPE_HARM)) sprintf(Shape,"harm");
127  else if(VAR_IF_TYPE(iShape,SHAPE_UMBR)) sprintf(Shape,"umbr");
128  else if(VAR_IF_TYPE(iShape,SHAPE_BOUND)) sprintf(Shape,"bound");
129 }
130 void VarData::StringNano(char *String,int n){
131  char Shape[60];
132  ShapeId(Nano[n].Shape,Shape);
133  if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_CLUSTER)){
134  sprintf(Shape,"Architecture%d.dat",n);
135  sprintf(String,"# Pep x(%.2f %.2f %.2f) a(%.2f %.2f %.2f) g(%.2f %.2f %.2f) i[%d] d[%d %d] fn{%s}",Nano[n].Pos[0],Nano[n].Pos[1],Nano[n].Pos[2],Nano[n].Axis[0],Nano[n].Axis[1],Nano[n].Axis[2],Nano[n].Rad,Nano[n].Height,Nano[n].Hamaker,n,Nano[n].NCircle,Nano[n].NHeight,Shape);
136  }
137  else if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_CLINKS)){
138  sprintf(Shape,"CrossLinks.dat");
139  sprintf(String,"# Ext i[%d] fn{%s}",0,Shape);
140  }
141  else
142  sprintf(String,"# Rigid x(%.2f %.2f %.2f) a(%.2f %.2f %.2f) c(%.2f %.2f %.2f %.4f) s{%s}",Nano[n].Pos[0],Nano[n].Pos[1],Nano[n].Pos[2],Nano[n].Axis[0],Nano[n].Axis[1],Nano[n].Axis[2],Nano[n].Rad,Nano[n].Hamaker,Nano[n].Height,Nano[n].OffSet,Shape);
143 }
144 int VarData::HeaderNano(FILE *FileToWrite){
145  if(Gen->NNano == 0) return 1;
146  for(int n=0;n<Gen->NNano;n++){
147  if(Nano->Shape == SHAPE_NONE)continue;
148  char String[120];
149  StringNano(String,n);
150  fprintf(FileToWrite,"%s\n",String);
151  }
152  return 0;
153 }
154 int VarData::HeaderSoft(char *Line){
155  if(NSoft == 0) return 1;
156  for(int n=0;n<NSoft;n++){
157  char Shape[60];
158  sprintf(Shape,"no");
159  if(Soft[n].Topology == VAR_PLANAR) sprintf(Shape,"planar");
160  if(Soft[n].Topology == VAR_VESICLE) sprintf(Shape,"vesicle");
161  if(Soft[n].Topology == VAR_COATING) sprintf(Shape,"coating");
162  if(Soft[n].Topology == VAR_DISTRIBUTED) sprintf(Shape,"distributed");
163  sprintf(Line,"# Soft x(%.2f %.2f %.2f) c(%.2f %.2f) s{%s} n{%s}\n",Soft[n].Pos[0],Soft[n].Pos[1],Soft[n].Pos[2],Soft[n].Size[0],Soft[n].Size[1],Shape,Soft[n].Name);
164  }
165  return 0;
166 }
167 bool VarData::WriteXvt(char *OutFile){
168  VarMessage("Write");
169  FILE *FileToWrite = NULL;
170  char Line2Put[512];
171  FileToWrite = fopen(OutFile,"w");
172  SigErr(FileToWrite==NULL,"The output file %s could not be opened\n",OutFile);
173  fprintf(FileToWrite,"# L=%lf %lf %lf t=%lf blocks=%d\n",Gen->Edge[0],Gen->Edge[1],Gen->Edge[2],pTime(),Gen->NBlock);
174  HeaderInteraction(FileToWrite);
175  HeaderNano(FileToWrite);
176  for(int b=0,cOff=0,pCurr=0;b<Gen->NBlock;b++,cOff+=Block[b].NChain){
177  if(Block[b].NChain == 0) continue;
178  fprintf(FileToWrite,"# n=%d N=%d name=%s\n",Block[b].NChain,Block[b].NPCh,Block[b].Name);
179  for(int c=cOff;c<Block[b].NChain+cOff;c++){
180  for(int p=0;p<Block[b].NPCh;p++,pCurr++){
181  // if(p > 0){
182  // for(int d=0;d<3;d++){
183  // if(Pm[pCurr].Pos[d] - Pm[pCurr-1].Pos[d] > .5*pEdge(d))
184  // Pm[pCurr].Pos[d] -= pEdge(d);
185  // else if(Pm[pCurr].Pos[d] - Pm[pCurr-1].Pos[d] < -.5*pEdge(d))
186  // Pm[pCurr].Pos[d] += pEdge(d);
187  // }
188  // }
189  fprintf(FileToWrite,"%lf %lf %lf ",
190  pPos(pCurr,0),pPos(pCurr,1),pPos(pCurr,2));
191  fprintf(FileToWrite,"%lf %lf %lf ",
192  Pm[pCurr].Vel[0],Pm[pCurr].Vel[1],Pm[pCurr].Vel[2]);
193  fprintf(FileToWrite," %d ",
194  Pm[pCurr].Typ);
195  fprintf(FileToWrite,"\n");
196  }
197  }
198  }
199  fclose(FileToWrite);
200  return 0;
201 }
202 bool VarData::WriteXyz(char *OutFile){
203  VarMessage("Write");
204  FILE *FileToWrite = NULL;
205  if( (FileToWrite = fopen(OutFile,"w")) == 0){
206  printf("The output file %s could not be opened\n",OutFile);
207  return 0;
208  }
209  for(int p=0;p<Gen->NPart;p++){
210  fprintf(FileToWrite,"%lf %lf %lf\n",
211  Pm[p].Pos[0],Pm[p].Pos[1],Pm[p].Pos[2]);
212  }
213  fclose(FileToWrite);
214  return 0;
215 }
216 void VarData::WriteLinkedSurf(FILE *FWrite,double *Plot,int Values,int NType,double *Bound,int* PId){
217  int p = PId[0];
218  int c = PId[1];
219  int t = PId[2];
220  double InvValues = 1./(double)Values;
221  int link[4] = {0,0,0,0};
222  for(int vv=0;vv<Values-1;vv++){
223  for(int v=0;v<Values-1;v++){
224  link[0] = (v+0)*Values + (vv+0);
225  if(Plot[(link[0])*NType+t] < .1) continue;
226  //------------defines-the-squares---------------------
227  link[1] = (v+0)*Values + (vv+1);
228  link[2] = (v+1)*Values + (vv+0);
229  link[3] = (v+1)*Values + (vv+1);
230  int pRef = p;
231  for(int lx=0;lx<2;lx++){
232  for(int ly=0;ly<2;ly++){
233  int l = 2*lx+ly;
234  double NanoAdded = Plot[(link[l])*NType+2]/Bound[2]+Plot[(link[l])*NType+3]/Bound[3];
235  double Phob = t == 0 ? Plot[(link[l])*NType+0]/Bound[0] : 0.;
236  double Phil = t == 1 ? Plot[(link[l])*NType+1]/Bound[1] : 0.;
237  double r = (v +lx)*InvValues*pEdge(3);
238  double z = (vv+ly)*InvValues*(pEdge(CNorm));
239  double dens = (Plot[(link[l])*NType+t]);//+Plot[(v*Values+vv)*NType+1]+Plot[(v*Values+vv)*NType+2]);
240  int l1 = pRef + (p+1)%4;
241  int l2 = pRef + (p+2)%4;
242  int l3 = pRef + (p+3)%4;
243  fprintf(FWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",p,c,t,r,z,dens,NanoAdded,Phob,Phil,l1,l2,l3);
244  p++;
245  }
246  }
247  c++;
248  }
249  }
250  PId[0] = p;
251  PId[1] = c;
252  PId[2] = t;
253  return ;
254  double Threshold = .1;
255  int *IdSerial = (int *)calloc(Values*Values*2,sizeof(int));
256  for(int vv=0,p=0;vv<Values-1;vv++){
257  for(int v=0;v<Values-1;v++){
258  int l = (vv+0)*Values + (v+0);
259  IdSerial[l] = p;
260  printf("%d %d %d %d\n",l,p,v,vv);
261  if(Plot[l*NType+t] < Threshold) continue;
262  p++;
263  }
264  }
265  for(int vv=0,p=0,c=0;vv<Values-1;vv++){
266  for(int v=0;v<Values-1;v++){
267  link[0] = (v+0)*Values + (vv+0);
268  if(Plot[(link[0])*NType+t] < Threshold) continue;
269  //------------defines-the-squares---------------------
270  link[1] = (vv+0)*Values + (v+1);
271  link[2] = (vv+1)*Values + (v+0);
272  link[3] = (vv+1)*Values + (v+1);
273  printf("%d %d %d %d\n",link[0],link[1],link[2],link[3]);
274  for(int l=0;l<4;l++){
275  link[l] = IdSerial[link[l]];
276  }
277  printf("%d %d %d %d\n",link[0],link[1],link[2],link[3]);
278  double NanoAdded = Plot[(link[0])*NType+2]/Bound[2]+Plot[(link[0])*NType+3]/Bound[3];
279  double Phob = t == 0 ? Plot[(link[0])*NType+0]/Bound[0] : 0.;
280  double Phil = t == 1 ? Plot[(link[0])*NType+1]/Bound[1] : 0.;
281  double r = (v)*InvValues*pEdge(3);
282  double z = (vv)*InvValues*(pEdge(CNorm));
283  double dens = (Plot[(link[0])*NType+t]);//+Plot[(v*Values+vv)*NType+1]+Plot[(v*Values+vv)*NType+2]);
284  fprintf(FWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",link[0],c,t,r,z,dens,NanoAdded,Phob,Phil,link[1],link[2],link[3]);
285  p++;
286  c++;
287  }
288  }
289  free(IdSerial);
290 }
291 void VarData::WriteSurf(FILE *F2Write,double **Plot,int NSample,int OffSet){
292  double InvNSample = 1./(double)NSample;
293  for(int sx=0;sx<NSample;sx++){
294  for(int sy=0;sy<NSample;sy++){
295  double x = (sx+.5)*pEdge(CLat1)*InvNSample;
296  double y = (sy+.5)*pEdge(CLat2)*InvNSample;
297  int l01 = (sx+0)*NSample+(sy+1)+OffSet;
298  if(sy == NSample - 1)
299  l01 = (sx+0)*NSample+(0)+OffSet;
300  int l10 = (sx+1)*NSample+(sy+0)+OffSet;
301  if(sx == NSample - 1)
302  l10 = (0)*NSample+(sy+0)+OffSet;
303  int lm0 = (sx-1)*NSample+(sy+0)+OffSet;
304  if(sx == 0)
305  lm0 = (NSample - 1)*NSample+(sy+0)+OffSet;
306  int l0m = (sx+0)*NSample+(sy-1)+OffSet;
307  if(sy == 0)
308  l0m = (sx)*NSample+(NSample-1)+OffSet;
309  fprintf(F2Write,"{x(%lf %lf %lf) l[%d] l[%d] l[%d] l[%d]}\n",x,y,Plot[sx][sy],l01,l10,lm0,l0m);
310  }
311  }
312 }
int pNLink()
Maximum number of bonds.
double pTime()
Total time.
Definition: VarData.h:908
double Deltat
Timestep.
Definition: VarData.h:337
double kappaN
Incompressibility.
Definition: VarData.h:325
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
void ShapeId(int iShape, char *Shape)
Identifier of the shape.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double chiN
Incompatibility.
Definition: VarData.h:323
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
void StringNano(char *NString, int n)
String for the rigid inclusion in the header file.
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
char cWhat2Draw[STRSIZE]
What to draw.
Definition: VarData.h:1042
double WFuncStraight2
Weighting function straight length.
Definition: VarData.h:339
int CLat2
lateral coordinate
Definition: VarData.h:1078
int NLink
Maximum number of bonds.
Definition: VarData.h:357
void WriteLinkedSurf(FILE *FWrite, double *Plot, int NSample, int NType, double *Bound, int *PId)
Write the positions of the egdes of the rectangles.
int Shape
0 none, 1 spherical, 2 cylindrical 3 wall
Definition: VarData.h:473
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int pNPCh()
Number of particle per chain.
bool WriteXvt(char *OutFile)
Writes a "system-file" or a "x y z" file".
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
int NChain
Number of chain.
Definition: VarData.h:347
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double kappaSpring
Prefactor of the bond spring.
Definition: VarData.h:331
bool WriteTxvl(char *OutFile)
Writes a "system-file" or a "x y z" file".
int HeaderNano(FILE *FileToWrite)
Write the nano section of the header to the file.
int NPCh
Number of particle per chain.
Definition: VarData.h:349
int SysType
Contains the definition of the system.
Definition: VarData.h:1086
double Coeff(int t1, int t2)
Prefactor of the force.
Definition: VarData.cpp:246
void WriteSurf(FILE *F2Write, double **Plot, int NSample, int OffSet)
Write the particle position as linked edges of squares.
int NPCh
particles per chain
Definition: VarData.h:259
bool WriteXyz(char *OutFile)
Writes a "system-file" or a "x y z" file".
int SysFormat
Contains the definition of the file format.
Definition: VarData.h:1088
void HeaderInteraction(FILE *FileToWrite)
Header interactions.
double pReOverCutOff()
Re/CutOff.
Definition: VarData.h:946
SOFT * Soft
Soft bodies.
Definition: VarData.h:1052
int NChain
chains
Definition: VarData.h:263
double WFuncStraight3
Weighting function straight length.
Definition: VarData.h:341
double pPos(int p, int d)
Return back folded position.
unsigned long Step
Courrent step.
Definition: VarData.h:343
int HeaderSoft(char *Line)
Header soft.
int CLat1
lateral coordinate
Definition: VarData.h:1076
int NSoft
Number of soft bodies.
Definition: VarData.h:1062
void Rescale(double SFactor, int Order)
Rescale entries.
Definition: VarData.cpp:312
int NPart
Number of particle.
Definition: VarData.h:345
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
double ReOverCutOff
Convertion unit R_e over CutOff.
Definition: VarData.h:335
int NBlock
Number of blocks.
Definition: VarData.h:359
void VarMessage(const char *s,...)
If enabled call the function position.
Definition: VarData.cpp:29
double rho
Density coexistence.
Definition: VarData.h:327
int CNorm
Normal coordinate.
Definition: VarData.h:1074
double kappaBend
Prefactor of the bending potential.
Definition: VarData.h:329