Allink  v0.1
ElPolyEl.cpp
1 /***********************************************************************
2 ElPoly:This progam provide a graphical visualisation of the data
3 opend by VarData using openGL glut. The most important option are
4 the possibility of changing the backfold of the polymers with 'c',
5 see the subsequent file in the list with '>', see the bond with 'b'.
6 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
7 
8 
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2 of the License, or
12 (at your option) any later version.
13 
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along with this program; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 ***********************************************************************/
23 #include "ElPoly.h"
24 ;
25 
26 ElPoly::ElPoly(int argc,char **argv,int *FilePos){
27  cFile = (char **)calloc(argc,4);//sizeof(**cFile));
28  if(cFile == NULL){printf("cFile not allocated\n");return;}
29  for(int f=0;f<argc;f++){
30  cFile[f] = (char *)calloc(STRSIZE,sizeof(*cFile));
31  if(cFile[f] == NULL){printf("cFile not allocated\n");return;}
32  cFile[f] = argv[FilePos[f]];
33  }
34  NFileTot = argc;
35  NFile[0] = 0; NFile[1] = NFileTot;
36  NPro =0;//?
37  quando=0;
38  if(!pNLink()) IfLine = 1;
39  else IfLine = 0;
40  // IfIntorno=0;//?
41  // IfColour=0;//?
42  // IfChains=0;//?
43  // IfChType=0;//?
44  // Diameter=0.1;//?
45  // StepDiameter=0.;//?
46  // ExtraDiam=0.;//?
47  // Vicinanze=1.;//?
48  InvScaleUn=1.;
49  ScaleFact = 1.;
50  LineSize = 1;
51  NVisSkip = 1;
52  NChType = CHAIN_EVERY;
53  NBackFold = BF_PART;
54  What2Draw = EL_PART;
55  for(int d=0;d<3;d++) ShiftPos[d] = 0.;
56  for(int d=0;d<3;d++) ScaleF[d] = 1.;
57 #ifdef OMPI_MPI_H
58  MPI_Init(&argc,&argv);
59  int Rank=0,Size=0;
60  MPI_Comm_rank(MPI_COMM_WORLD, &Rank);
61  MPI_Comm_size(MPI_COMM_WORLD, &Size);
62  int Partition = (int)(NFileTot/(double)Size);
63  NFile[0] = Partition*Rank;
64  NFile[1] = Partition*(Rank+1);
65  if(Rank==Size-1) NFile[1] += NFileTot%Size;
66  Proc = new SingProc(Size,Rank);
67 #endif
68  //TrialSys();
69 }
71  // for(int f=0;f<NFileTot;f++)
72  // free(cFile[f]);
73  // free(cFile);
74 #ifdef OMPI_MPI_H
75  MPI_Finalize();
76 #endif
77 }
78 int ElPoly::Angle(int Values){
79  int NRadici = 4;
80  int NZeri = NRadici;
81  double *Radici;
82  Radici = (double *)malloc(NRadici*sizeof(double));
83  double *Average;
84  Average = (double *)malloc(NRadici*sizeof(double));
85  for(int i=0;i<NRadici;i++){
86  Radici[i] = 0.;
87  Average[i] = 0.;
88  }
89  for(int f=NFile[0];f<NFile[1];f++){
90  Processing(f);
91  if(OpenRisk(cFile[f],BF_PART)==0)return 1;
92  Mat->Ypsilon = pEdge(2)-pCm(2)-1.;
93  Mat->PreFact = pow((2.*pNPart()/10.)/DUE_PI,.33);
95  printf("Volume %lf Cm %lf # to Invert %lf\n",(double)pNPart()/10.,pCm(2),pow(DUE_PI/(2.*pNPart()/10.),.25));
96  NZeri = Mat->Zeri(0.,DUE_PI/2.,Radici,NRadici);
97  for(int i=0;i<NZeri;i++){
98  Radici[i] *= 360./DUE_PI;
99  Average[i] += Radici[i];
100  }
101  if(NZeri == 0){
102  printf("The function has no real roots\n");
103  }
104  }
105  for(int i=0;i<NZeri;i++){
106  Average[i] /= NFile[1]-NFile[0];
107  printf("Angle %lf\n",Average[i]);
108  }
109  return 0;
110 }
111 int ElPoly::CenterOfMass(int Coord){
112  double SysCm;
113  FILE *FileToWrite = NULL;
114  FileToWrite = fopen("CenterOfMass.dat","w");
115  for(int f=NFile[0];f<NFile[1];f++){
116  Processing(f);
117  if(OpenRisk(cFile[f],1)==0)return 1;
118  SysCm=0.;
119  if(Coord == 3){
120  for(int p=0;p<pNPart();p++){
121  SysCm += pVel(p,2);
122  }
123  }
124  else{
125  for(int p=0;p<pNPart();p++){
126  SysCm += pPos(p,Coord);
127  }
128  }
129  fprintf(FileToWrite,"%lf \n",SysCm/pNPart());
130  }
131  printf("\n");
132  fclose(FileToWrite);
133  return 0;
134 }
136  printf("Introduce the first of %d file from the list: ",NFileTot);
137  scanf("%d",&NFile[0]);
138  if(NFile[0] < 0 || NFile[0] > NFile[1])
139  NFile[0] = 0;
140  quando = NFile[0];
141  printf("Introduce the second of %d file from the list: ",NFileTot);
142  scanf("%d",&NFile[1]);
143  if(NFile[1] < NFile[0] || NFile[1] > NFileTot)
144  NFile[1] = NFileTot;
145  return 0;
146 }
147 int ElPoly::OpenFile(int f){
148  if(f >= NFileTot || f < 0)
149  return 0;
150  if(Open(cFile[f],NBackFold)) return 1;
151  NFile[0] = f;
152  quando = f;
153  return 0;
154 }
155 int ElPoly::OpenFile(char *FileName){
156  if(Open(FileName,NBackFold))return 1;
157  return 0;
158 }
159 double ElPoly::ContactAngle(double x){
160  double Ris;
161  Ris = QUAD( 2*x - QUAD(sin(2*x)) ) / (4./3. * CUB(sin(x))-2.*x*cos(x)+cos(x)*sin(2*x) ) - 5.;
162  return Ris;
163 }
164 // OpenGl----------------------------------
166  Properties Pr;
167  for(int f=NFile[0];f<NFile[1];f++){
168  Processing(f);
169  if(OpenRisk(cFile[f],BF_PART))return 0;
170  Pr = Pr + SysProperties();
171  }
172  Pr = Pr * (1./(NFile[1] - NFile[0]));
173  printf("\n");
174  char *NomeFile = (char *) malloc(60*sizeof(char));
175  sprintf(NomeFile,"PropertiesChi%.0fRho%.0fKappa%.0f.txt",pchiN(),prho(),pkappaN());
176  FILE *FileToWrite = NULL;
177  FileToWrite = fopen(NomeFile,"w");
178  fprintf(FileToWrite,"RePhob %lf RePhil %lf\n",Pr.RePhob,Pr.RePhil);
179  fprintf(FileToWrite,"RadGyrPhob %lf RadGyrPhil %lf\n",Pr.GyrPhob,Pr.GyrPhil);
180  fprintf(FileToWrite,"StructFactPhob %lf StructFactPhil %lf\n",Pr.FactPhob,Pr.FactPhil);
181  fprintf(FileToWrite,"VolPhob %lf VolPhil %lf\n",Pr.VolPhob,Pr.VolPhil);
182  fprintf(FileToWrite,"ChainDiff %lf\n",Pr.ChDiff);
183  Pr.Print();
184  fclose(FileToWrite);
185  return 0;
186 }
187 void ElPoly::DivideLayers(int How){
188  int NLayer = 2;
189  if(How == VAR_OPPOSED) NLayer = 4;
190  int *NChStep = (int *)calloc(NLayer,sizeof(int));
191  int *ChainLevel = (int *)calloc(pNChain(),sizeof(int));
192  if(How == VAR_OPPOSED){
193  for(int c=0,b=0;c<pNChain();c++){
194  if(c*Block[b].NPCh > Block[b].EndIdx) b++;
195  int Level = 0;
196  if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_UP))
197  Level = 1;
198  if(!strcmp(Block[b].Name,"LIPID1")) Level += 2;
199  NChStep[Level]++;
200  ChainLevel[c] = Level;
201  }
202  }
203  else if(How == VAR_TUBE || How == VAR_VESICLE){
204  for(int c=0,b=0;c<pNChain();c++){
205  if(c*Block[b].NPCh > Block[b].EndIdx) b++;
206  int Level = 0;
207  if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_OUTER))
208  Level = 1;
209  NChStep[Level]++;
210  ChainLevel[c] = Level;
211  }
212  }
213  //for(int c=0;c<pNChain();c++)printf("%lf %d %d\n",Ch[c].Pos[CNorm],c,ChainLevel[c]);
214  for(int f=NFile[0];f<NFile[1];f++){
215  Processing(f);
216  if(OpenRisk(cFile[f],BF_NO))return ;
217  PART *Pn = (PART *)calloc(pNPart(),sizeof(PART));
218  int pp=0;
219  for(int p=0;p<pNPart();p++){
220  for(int d=0;d<3;d++){
221  Pn[p].Pos[d] = pPos(p,d);
222  Pn[p].Vel[d] = pVel(p,d);
223  }
224  Pn[p].Typ = pType(p);
225  }
226  for(int l=0;l<NLayer;l++){
227  for(int p=0;p<pNPart();p++){
228  if(ChainLevel[pChain(p)] == l){
229  for(int d=0;d<3;d++){
230  SetPos(pp,Pn[p].Pos);
231  SetVel(pp,Pn[p].Vel);
232  }
233  SetType(pp,Pn[p].Typ);
234  pp++;
235  }
236  }
237  }
238  DefBlock(NChStep,How);
239  Write(cFile[f]);
240  free(Pn);
241  }
242  printf("\n");
243  free(ChainLevel);
244  free(NChStep);
245 }
247  int nChain = pNChain();
248  for(int b=0;b<pNBlock();b++){
249  for(int c=0;c<Block[b].NChain;c++){
250  if(Ch[c].Pos[CLat1] > .4*pEdge(CLat1) && Ch[c].Pos[CLat1] < .6*pEdge(CLat1) ){
251  SwapChain(c,nChain-1,b);
252  nChain--;
253  SetNPart(pNPart()-pNPCh());
254  }
255  }
256  Block[b].NChain = nChain;
257  }
258  SetNChain(nChain);
259  Block[0].EndIdx = Block[0].NChain*Block[0].NPCh;
260  for(int b=1;b<pNBlock();b++){
261  Block[b].InitIdx = Block[0].NChain*Block[b].NPCh + Block[b-1].EndIdx;
262  }
263  Write("cleft.dat");
264 }
265 // int ElPoly::DensPatch(int NPatch,int Values){
266 // int NType = 2;
267 // double *DensProf = (double *)calloc(NPatch*Values,sizeof(double));
268 // for(int f=NFile[0];f<NFile[1];f++){
269 // Processing(f);
270 // fprintf(stderr,"Processing: %d/%d %.1f %%\r",f,(NFile[1]-NFile[0]),(f-NFile[0])/(double)(NFile[1]-NFile[0])*100.);
271 // if(OpenRisk(cFile[f],BF_PARTICLE))return 0;
272 
273 // }
274 // char *FileName = (char *)calloc(60,sizeof(char));
275 // sprintf(FileName,"ChainPArea%.0fRho%.0fChiN%.0fKappa.dat",Gen->rho,Gen->chiN,Gen->kappaN);
276 // FileToWrite = fopen(FileName,"w");
277 // fprintf(FileToWrite,"#%s\n",SysInfo());
278 // }
279 char *ElPoly::ChooseDraw(int ExtWhat2Draw){
280  char *String = (char *)calloc(16,sizeof(char));
281  if(ExtWhat2Draw == EL_PART)
282  sprintf(String,"part");
283  else if(ExtWhat2Draw == EL_PART)
284  sprintf(String,"part");
285  else if(ExtWhat2Draw == EL_CHAIN)
286  sprintf(String,"chain");
287  else if(ExtWhat2Draw == EL_DENS)
288  sprintf(String,"dens");
289  else if(ExtWhat2Draw == EL_QUAD)
290  sprintf(String,"quad");
291  else if(ExtWhat2Draw == EL_QUAD1)
292  sprintf(String,"quad1");
293  else if(ExtWhat2Draw == EL_POTENTIAL)
294  sprintf(String,"pot");
295  else if(ExtWhat2Draw == EL_SURF)
296  sprintf(String,"surf");
297  else if(ExtWhat2Draw == EL_CROSS)
298  sprintf(String,"cross");
299  else if(ExtWhat2Draw == EL_ISOIPSE)
300  sprintf(String,"isoipse");
301  else if(ExtWhat2Draw == EL_SAMPLE)
302  sprintf(String,"sample");
303  else if(ExtWhat2Draw == EL_COLOR)
304  sprintf(String,"color");
305  else if(ExtWhat2Draw == EL_SKIN)
306  sprintf(String,"skin");
307  else if(ExtWhat2Draw == EL_POLYGON)
308  sprintf(String,"polygon");
309  else if(ExtWhat2Draw == EL_ISOLEVEL)
310  sprintf(String,"isolevel");
311  else if(ExtWhat2Draw == EL_SQUAREMESH)
312  sprintf(String,"squaremesh");
313  return String;
314 }
315 void ElPoly::ChooseDraw(char *String){
316  if(!strcmp(String,"part"))
317  What2Draw = EL_PART;
318  else if(!strcmp(String,"quad"))
319  What2Draw = EL_QUAD;
320  else if(!strcmp(String,"chain"))
321  What2Draw = EL_CHAIN;
322  else if(!strcmp(String,"dens"))
323  What2Draw = EL_DENS;
324  else if(!strcmp(String,"quad1"))
325  What2Draw = EL_QUAD1;
326  else if(!strcmp(String,"pot"))
327  What2Draw = EL_POTENTIAL;
328  else if(!strcmp(String,"surf"))
329  What2Draw = EL_SURF;
330  else if(!strcmp(String,"cross"))
331  What2Draw = EL_CROSS;
332  else if(!strcmp(String,"isoipse"))
333  What2Draw = EL_ISOIPSE;
334  else if(!strcmp(String,"sample"))
335  What2Draw = EL_SAMPLE;
336  else if(!strcmp(String,"color"))
337  What2Draw = EL_COLOR;
338  else if(!strcmp(String,"stalk"))
339  What2Draw = EL_STALK;
340  else if(!strcmp(String,"skin"))
341  What2Draw = EL_SKIN;
342  else if(!strcmp(String,"polygon"))
343  What2Draw = EL_POLYGON;
344  else if(!strcmp(String,"isolevel"))
345  What2Draw = EL_ISOLEVEL;
346  else if(!strcmp(String,"squaremesh"))
347  What2Draw = EL_SQUAREMESH;
348  else
349  What2Draw = 0;
350 }
351 void ElPoly::SetBoundFile(int InitFile,int EndFile){
352  if(InitFile < 0)
353  NFile[0] = 0;
354  else if(InitFile < EndFile){
355  NFile[0] = InitFile;
356  quando = InitFile;
357  }
358  else {
359  printf("Wrong init file %d !< %d\n",InitFile,EndFile);
360  }
361  if(EndFile >= NFileTot)
362  NFile[1] = NFileTot-1;
363  else if(EndFile > InitFile)
364  NFile[1] = EndFile;
365  else {
366  printf("Wrong end file %d !< %d\n",InitFile,EndFile);
367  }
368 }
369 void ElPoly::Processing(int f){
370  fprintf(stderr,"Processing: %s %d/%d %.1f %%\r",cFile[f],f,(NFile[1]-NFile[0]),(f-NFile[0])/(double)(NFile[1]-NFile[0])*100.);
371 }
373  double sigma =
374  sqrt(QUAD(pReOverCutOff())/(pNPCh()-1)/3.);
375  for(int f=NFile[0];f<NFile[1];f++){
376  Processing(f);
377  Open(cFile[f],BF_NO);
378  // for(int d=0;d<3;d++){
379  // ShiftPos[d] = .5 - pCm(d)*pInvEdge(d);
380  // }
381  // ShiftRef(BF_CHAIN);
382  double Temp = pEdge(CLat1);
384  SetEdge(Temp,CNorm);
385 
386  // FILE *FWrite = fopen("Response.dat","w");
387  // fprintf(FWrite,"# L=%lf %lf %lf t=%lf blocks=%d\n",pEdge(0),pEdge(1),pEdge(2),pTime(),pNBlock());
388  // HeaderInteraction(FWrite);
389  // HeaderNano(FWrite);
390  // fprintf(FWrite,"# n=%d N=%d name=%s\n",Block[0].NChain,Block[0].NPCh,Block[0].Name);
391  // int NChain = 0;
392  // for(int c=0;c<pNChain();c++){
393  // double Dist = 0.;
394  // for(int d=0;d<3;d++){
395  // Dist += SQR(Nano->Pos[d] - Ch[c].Pos[d]);
396  // }
397  // if(Dist > SQR(2.5*Nano->Rad)) continue;
398  // for(int p=c*pNPCh();p<(c+1)*pNPCh();p++){
399  // fprintf(FWrite,"%lf %lf %lf %lf %lf %lf %d\n",
400  // pPos(p,0),pPos(p,1),pPos(p,2),
401  // pVel(p,0),pVel(p,1),pVel(p,2),pType(p));
402  // }
403  // NChain++;
404  // }
405  // printf("Number of chains %d\n",NChain);
406  // return;
407 
408  for(int p=0;p<pNPart();p++){
409  double Pos[3] = {pPos(p,CNorm),pPos(p,CLat1),pPos(p,CLat2)};
410  SetPos(p,Pos);
411  // for(int d=0;d<3;d++){
412  // Pm[p].Vel[d] = Mat->Gaussiano(0.,sigma);
413  // }
414  }
415  Write(cFile[f]);
416  }
417 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int pNLink()
Maximum number of bonds.
void SetBoundFile(int InitFile, int EndFile)
Set the initial and final number of files.
Definition: ElPolyEl.cpp:351
double Ypsilon
External parameter to calculate the contact angle.
Definition: Matematica.h:330
double InvScaleUn
Define the shrink factor between the box edges.
Definition: ElPoly.h:457
int LineSize
Line size of the gl.
Definition: ElPoly.h:427
bool OpenRisk(char *InFile, int BF)
Opens a file without reallocationg.
Definition: VarData.cpp:126
double prho()
Density coexistence.
Definition: VarData.h:944
double pVel(int p, int d)
Return the velocity.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
int NFileTot
Total number of file.
Definition: ElPoly.h:420
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
double ShiftPos[3]
Reference position.
Definition: VarData.h:1058
int SetNChain(int NewNCh)
Set and reallocate the number of chains.
Definition: VarDataComm.cpp:78
int OpenFile(int f)
Opens the f file of the list open.
Definition: ElPolyEl.cpp:147
int pChain(int p)
Return the chain.
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
int CenterOfMass(int Coord)
Averaged center of mass Boh.
Definition: ElPolyEl.cpp:111
int Zeri(double a, double b, double *Radici, int NRadici)
Find the.
void RemoveChains()
Remove chains satisfying a condition.
Definition: ElPolyEl.cpp:246
void SetPos(int p, double *Pos)
Set the particle position.
int quando
Current number of the file list.
Definition: ElPoly.h:424
Properties SysProperties()
Calculate some basis properties.
Definition: VarDataEl.cpp:234
int NChType
Type of chain selected.
Definition: VarData.h:1080
int pType(int p)
Return the type.
double pCm(int d)
Center of mass of the system.
Definition: VarData.h:924
int CLat2
lateral coordinate
Definition: VarData.h:1078
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
Some calculated properties of the system.
Definition: VarData.h:364
double ScaleFact
Normal scaling factor (z zooming)
Definition: ElPoly.h:461
void DivideLayers(int How)
Reorder the LIPID block in four different layers.
Definition: ElPolyEl.cpp:187
bool Open(char *InFile, int BF)
Open the.
Definition: VarData.cpp:105
FUNC Func
Pointer to a function.
Definition: Matematica.h:111
int InitIdx
Initial particle position.
Definition: VarData.h:265
int Angle(int NBin)
Calculation of the contact angle angle Boh?
Definition: ElPolyEl.cpp:78
int SetNPart(int NewNPart)
Set and reallocate the number of particles.
Definition: VarDataComm.cpp:2
~ElPoly()
Desctructor.
Definition: ElPolyEl.cpp:70
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int pNPCh()
Number of particle per chain.
int NPro
Boh.
Definition: ElPoly.h:422
int pNBlock()
Number of blocks.
double pchiN()
Incompatibility.
Definition: VarData.h:928
int ChangeFile()
Defines the first and the last file to be elaborated file.
Definition: ElPolyEl.cpp:135
double ContactAngle(double x)
Definition of the contact angle.
int PropertiesF()
Calculates some properties of the system prop.
Definition: ElPolyEl.cpp:165
void SetType(int p, int t)
Set the particle type.
Basics class to start a MPI grid.
Definition: SingProc.h:6
double PreFact
External parameter in the definition of the contact angle.
Definition: Matematica.h:332
int NBackFold
Type of backfold.
Definition: ElPoly.h:439
void SwapChain(int c1, int c2, int b)
Swap two chains.
Definition: VarDataEl.cpp:355
int NVisSkip
How many lipids are skipped in the visualization.
Definition: ElPoly.h:443
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int NPCh
particles per chain
Definition: VarData.h:259
void DefBlock(int *NChStep, int How)
Define four different blocks.
double ScaleF[3]
Scale factor.
Definition: VarData.h:1060
ElPoly(int argc, char **argv, int *FilePos)
Constructor.
Definition: ElPolyEl.cpp:26
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
int IfLine
Boh.
Definition: ElPoly.h:435
double pReOverCutOff()
Re/CutOff.
Definition: VarData.h:946
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int NChain
chains
Definition: VarData.h:263
int What2Draw
Draw the content in the appropriate visulization.
Definition: ElPoly.h:441
int EndIdx
End particle position.
Definition: VarData.h:267
double pPos(int p, int d)
Return back folded position.
double RePhob
End2End distance of the.
Definition: VarData.h:373
int pNChain()
Number of chain.
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
int CLat1
lateral coordinate
Definition: VarData.h:1076
Information of every particle.
Definition: VarData.h:214
void Print()
Print the current values.
Definition: VarData.h:369
int Typ
Type.
Definition: VarData.h:226
void SetVel(int p, double *Vel)
Set the particle velocity.
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
double ContactAngle(double x)
Definition of the contact angle function Boh.
Definition: ElPolyEl.cpp:159
double pkappaN()
Incompressibility.
Definition: VarData.h:930
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.
void Shift2Center()
Shift the system to the center.
Definition: ElPolyEl.cpp:372