Allink  v0.1
ElPolyProfDens.cpp
1 /***********************************************************************
2 ElPolyProfDens:
3 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
4 
5 
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10 
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 ***********************************************************************/
20 #include "ElPoly.h"
24 int ElPoly::DensProf(int NBin,int NSample,int Coord){
25  int NType = 3;
26  double *DensityAv = (double *)calloc(NType*NBin,sizeof(double));
27  double FreeE = 0.;
28  double ChDensMean = 0.;
29  double ChDensSqr = 0.;
30  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
31  double *VolEl = (double *)calloc(NBin,sizeof(double));
32  VolumeCircSlab(VolEl,NBin);
33  if(Coord > 3 || Coord < 0) return 1;
34  for(int f=NFile[0];f<NFile[1];f++){
35  Processing(f);
36  //printf("Opening: %s\n",cFile[f]);
37  if(Coord < 3)
38  if(OpenRisk(cFile[f],BF_PART))return 0;
39  else if(Coord == 3){
40  if(OpenRisk(cFile[f],BF_NANO))return 0;
41  }
42  double RefZed = .5*pCm(CNorm);
43  if(pNNano() > 0) RefZed = pNanoPos(0,CNorm);
44  // if(Coord == 2 && VAR_IF_TYPE(SysType,VAR_MEMBRANE) )
45  // ShiftSys(SHIFT_CM);
47  if(Coord < 3){
48  //DensityProfile(Coord,NBin,NType,DensityAv);
49  BfDefChain();
50  int NPartition = NSample;
51  double *Density = (double *)calloc(NPartition*NPartition*NType*NBin,sizeof(double));
52  double Pos[4];
53  char cLine[STRSIZE];
54  for(int b=0;b<pNBlock();b++){
55  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
56  if(!strcmp(Block[b].Name,"WATER") ) continue;
57  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
58  int t = pType(p);
59  for(int d=0;d<3;d++) Pos[d] = pPos(p,d);
60  Pos[CNorm] = pPos(p,CNorm) - RefZed + .5*pEdge(CNorm);
61  int vx = (int)(NPartition*Pos[CLat1]*pInvEdge(CLat2));
62  if(vx < 0 || vx >= NBin) continue;
63  int vy = (int)(NPartition*Pos[CLat2]*pInvEdge(CLat2));
64  if(vy < 0 || vy >= NBin) continue;
65  int vz = (int)(NBin*Pos[CNorm]*pInvEdge(CNorm));
66  if(vz < 0 || vz >= NBin) continue;
67  if(CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED)) t = 2;
68  if(!strncmp(Block[b].Name,"ADDED",5)) t = 2;
69  if(!strncmp(Block[b].Name,"CHOL",4)) t = 2;
70  Density[((vz*NType+t)*NPartition+vy)*NPartition+vx] += 1.;
71  }
72  }
73  for(int vx=0;vx<NPartition;vx++){
74  for(int vy=0;vy<NPartition;vy++){
75  double Media = 0.;
76  double Weight = 0.;
77  for(int v=0;v<NBin;v++){
78  Weight += Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
79  Media += v*Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
80  }
81  int vMedia = Weight > 0 ? (int)(Media/Weight) : NBin/2;
82  int vOffset = vMedia - NBin/2;
83  for(int t=0;t<NType;t++){
84  for(int v=0;v<NBin;v++){
85  int vOther = v + vOffset;
86  if(vOther < 0)
87  vOther = NBin - vOther;
88  if(vOther >= NBin)
89  vOther = vOther - NBin;
90  DensityAv[v*NType+t] += Density[((vOther*NType+t)*NPartition+vy)*NPartition+vx];
91  }
92  }
93  }
94  }
95  free(Density);
96  }
97  else{
98  Vettore NanoAxis(Nano[0].Axis[0],Nano[0].Axis[1],Nano[0].Axis[2]);
99  Vettore PosRel(3);
100  Vettore Dist(3);
101  for(int b=0;b<pNBlock();b++){
102  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
103  for(int p=0;p<pNPart();p++){
104  for(int d=0;d<3;d++){
105  PosRel.Set(remainder(pNanoPos(0,d)-pPos(p,d),pEdge(d)),d);
106  }
107  double RadDist = ABS(Dist.PerpTo3(&PosRel,&NanoAxis));
108  int v = (int)(RadDist/pEdge(3)*NBin);
109  if(v < 0 || v >= NBin) continue;
110  int t = pType(p);
111  if(CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED)) t = 2;
112  DensityAv[v*NType+t] += 1.;
113  }
114  }
115  }
116  ChDensMean += pNChain()/(pEdge(CLat1)*pEdge(CLat2));
117  ChDensSqr += SQR(pNChain()/(pEdge(CLat1)*pEdge(CLat2)));
118  }
119 #ifdef OMPI_MPI_H
120  MPI_Allreduce(MPI_IN_PLACE,DensityAv,NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
121  int Rank=0;
122  MPI_Comm_rank(Proc->CommGrid, &Rank);
123  if(Rank==0){
124 #endif
125  printf("\n");
126  double NormaF = (NFile[1] - NFile[0]);
127  double VolumeSlab = (pEdge(0)*pEdge(1)*pEdge(2))/(double)NBin;
128  for(int v=0;v<NBin;v++){
129  for(int t=0;t<NType;t++){
130  if(Coord < 3)
131  DensityAv[v*NType+t] /= VolumeSlab*NormaF;
132  else
133  DensityAv[v*NType+t] /= NormaF*VolEl[v];
134  }
135  }
136  char FileName[256];
137  sprintf(FileName,"DensityChi%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho());
138  FILE *FileToWrite = fopen(FileName,"w");
139  char cSystem[STRSIZE];
140  SysDef(cSystem);
141  double ChDensErr = sqrt(ChDensSqr - SQR(ChDensMean)/NormaF )/NormaF;
142  fprintf(FileToWrite,"# %s\n",cSystem);
143  fprintf(FileToWrite,"# ChainDens %lf %lf\n",ChDensMean/NormaF,ChDensErr);
144  for(int v=0;v<NBin;v++){
145  fprintf(FileToWrite,"%lf %lf %lf %lf\n",(double)v/(double)NBin*pEdge(Coord),DensityAv[v*NType],DensityAv[v*NType+1],DensityAv[v*NType+2]);
146  }
147  fclose(FileToWrite);
148 #ifdef OMPI_MPI_H
149  }
150 #endif
151  free(DensityAv);
152  return 0;
153 }
157 void ElPoly::DensProfNormalSlab(int NBin,int NSample,int Coord){
158  int NType = 3;
159  double *DensityAv = (double *)calloc(NType*NBin,sizeof(double));
160  double FreeE = 0.;
161  double ChDensMean = 0.;
162  double ChDensSqr = 0.;
163  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
164  double *VolEl = (double *)calloc(NBin,sizeof(double));
165  VolumeCircSlab(VolEl,NBin);
166  if(Coord > 3 || Coord < 0) return;
167  for(int f=NFile[0];f<NFile[1];f++){
168  Processing(f);
169  //printf("Opening: %s\n",cFile[f]);
170  if(Coord < 3)
171  if(OpenRisk(cFile[f],BF_PART))return;
172  else if(Coord == 3){
173  if(OpenRisk(cFile[f],BF_NANO))return;
174  }
175  double RefZed = .5*pCm(CNorm);
176  if(pNNano() > 0) RefZed = pNanoPos(0,CNorm);
177  // if(Coord == 2 && VAR_IF_TYPE(SysType,VAR_MEMBRANE) )
178  // ShiftSys(SHIFT_CM);
180  if(Coord < 3){
181  //DensityProfile(Coord,NBin,NType,DensityAv);
182  BfDefChain();
183  int NPartition = NSample;
184  double *Density = (double *)calloc(NPartition*NPartition*NType*NBin,sizeof(double));
185  double Pos[4];
186  char cLine[STRSIZE];
187  for(int b=0;b<pNBlock();b++){
188  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
189  if(!strcmp(Block[b].Name,"WATER") ) continue;
190  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
191  int t = pType(p);
192  for(int d=0;d<3;d++)
193  Pos[d] = pPos(p,d);
194  Pos[CNorm] = pPos(p,CNorm) - RefZed + .5*pEdge(CNorm);
195  int vx = (int)(NPartition*Pos[CLat1]*pInvEdge(CLat2));
196  if(vx < 0 || vx >= NBin) continue;
197  int vy = (int)(NPartition*Pos[CLat2]*pInvEdge(CLat2));
198  if(vy < 0 || vy >= NBin) continue;
199  int vz = (int)(NBin*Pos[CNorm]*pInvEdge(CNorm));
200  if(vz < 0 || vz >= NBin) continue;
201  if(CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED)) t = 2;
202  if(!strncmp(Block[b].Name,"ADDED",5)) t = 2;
203  if(!strncmp(Block[b].Name,"CHOL",4)) t = 2;
204  Density[((vz*NType+t)*NPartition+vy)*NPartition+vx] += 1.;
205  }
206  }
207  for(int vx=0;vx<NPartition;vx++){
208  for(int vy=0;vy<NPartition;vy++){
209  double Media = 0.;
210  double Weight = 0.;
211  for(int v=0;v<NBin;v++){
212  Weight += Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
213  Media += v*Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
214  }
215  int vMedia = Weight > 0 ? (int)(Media/Weight) : NBin/2;
216  int vOffset = vMedia - NBin/2;
217  for(int t=0;t<NType;t++){
218  for(int v=0;v<NBin;v++){
219  int vOther = v + vOffset;
220  if(vOther < 0)
221  vOther = NBin - vOther;
222  if(vOther >= NBin)
223  vOther = vOther - NBin;
224  DensityAv[v*NType+t] += Density[((vOther*NType+t)*NPartition+vy)*NPartition+vx];
225  }
226  }
227  }
228  }
229  free(Density);
230  }
231  else{
232  Vettore NanoAxis(Nano[0].Axis[0],Nano[0].Axis[1],Nano[0].Axis[2]);
233  Vettore PosRel(3);
234  Vettore Dist(3);
235  for(int b=0;b<pNBlock();b++){
236  if(!strncmp(Block[b].Name,"PEP",3)) continue;
237  for(int p=0;p<pNPart();p++){
238  for(int d=0;d<3;d++){
239  PosRel.Set(remainder(pNanoPos(0,d)-pPos(p,d),pEdge(d)),d);
240  }
241  double RadDist = ABS(Dist.PerpTo3(&PosRel,&NanoAxis));
242  int v = (int)(RadDist/pEdge(3)*NBin);
243  if(v < 0 || v >= NBin) continue;
244  int t = pType(p);
245  if(CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED))
246  t = 2;
247  DensityAv[v*NType+t] += 1.;
248  }
249  }
250  }
251  ChDensMean += pNChain()/(pEdge(CLat1)*pEdge(CLat2));
252  ChDensSqr += SQR(pNChain()/(pEdge(CLat1)*pEdge(CLat2)));
253  }
254 #ifdef OMPI_MPI_H
255  MPI_Allreduce(MPI_IN_PLACE,DensityAv,NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
256  int Rank=0;
257  MPI_Comm_rank(Proc->CommGrid, &Rank);
258  if(Rank==0){
259 #endif
260  printf("\n");
261  double NormaF = (NFile[1] - NFile[0]);
262  double VolumeSlab = (pEdge(0)*pEdge(1)*pEdge(2))/(double)NBin;
263  for(int v=0;v<NBin;v++){
264  for(int t=0;t<NType;t++){
265  if(Coord < 3)
266  DensityAv[v*NType+t] /= VolumeSlab*NormaF;
267  else
268  DensityAv[v*NType+t] /= NormaF*VolEl[v];
269  }
270  }
271  char FileName[256];
272  sprintf(FileName,"DensityChi%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho());
273  FILE *FileToWrite = fopen(FileName,"w");
274  char cSystem[STRSIZE];
275  SysDef(cSystem);
276  double ChDensErr = sqrt(ChDensSqr - SQR(ChDensMean)/NormaF )/NormaF;
277  fprintf(FileToWrite,"# %s\n",cSystem);
278  fprintf(FileToWrite,"# ChainDens %lf %lf\n",ChDensMean/NormaF,ChDensErr);
279  for(int v=0;v<NBin;v++){
280  fprintf(FileToWrite,"%lf %lf %lf %lf\n",(double)v/(double)NBin*pEdge(Coord),DensityAv[v*NType],DensityAv[v*NType+1],DensityAv[v*NType+2]);
281  }
282  fclose(FileToWrite);
283 #ifdef OMPI_MPI_H
284  }
285 #endif
286  free(DensityAv);
287  return;
288 }
292 void ElPoly::SlabProf(int NBin,int nNano,int Coord1){
293  if(Coord1 != CLat1 && Coord1 != CLat2){
294  printf("%d is not a lateral coordinate\n",Coord1,CLat1,CLat2);
295  }
296  //SigErr(pNNano() < 2,"The number of nano should be at least two");
297  int Coord2 = Coord1 == CLat1 ? CLat2 : CLat1;
298  double *Dens = (double *) calloc(NBin,sizeof(double));
299  double *Temp = (double *) calloc(NBin,sizeof(double));
300  double *Prof = (double *) calloc(6*NBin,sizeof(double));
301  double **Plot = (double **)calloc(3,sizeof(double));
302  for(int t=0;t<3;t++){
303  Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
304  }
305  double *Count = (double *) calloc(6*NBin,sizeof(double));
306  Vettore Pos(3);
307  Vettore PosRel(3);
308  //FILE *Ciccia = fopen("Ciccia.dat","w");
309  for(int f=NFile[0];f<NFile[1];f++){
310  Processing(f);
311  //--------------open,-back-fold----------------------
312  if(OpenRisk(cFile[f],BF_NO)) return;
313  if(pNNano() == 1){
314  SetNNano(2);
315  Nano[1].Pos[0] = 0.;
316  Nano[1].Pos[1] = 0.;
317  Nano[1].Pos[2] = .5*pEdge(CNorm);
318  }
319  Vettore Axis(3);
320  for(int d=0;d<3;d++){
321  Axis.Set(Nano[1].Pos[d]-Nano[0].Pos[d],d);
322  }
323  Axis.Set(0.,CNorm);
324  Axis.Normalize();
325  double InvEdge = (pEdge(CLat1)*Axis.Val(0) + pEdge(CLat2)*Axis.Val(1))/(Axis.Val(0) + Axis.Val(1));
326  InvEdge = 1./InvEdge;
327  //Nano[nNano].Pos[CLat1] = 0.;
328  //Nano[nNano].Pos[CLat2] = 0.;
329  for(int b=0;b<pNBlock();b++){
330  int TypePlus = 0;
331  if(!strncmp(Block[b].Name,"PEP",3) ) TypePlus = 1;
332  if(!strcmp(Block[b].Name,"WATER") ) continue;
333  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
334  int t = 2*pType(p);
335  for(int d=0;d<3;d++){
336  double Dist = Pm[p].Pos[d] - Nano[nNano].Pos[d];
337  Dist -= floor(Dist*pInvEdge(d)+.5)*pEdge(d);
338  Pos.Set(Dist,d);
339  }
340  Pos.Set(0.,CNorm);
341  double LatDist = PosRel.PerpTo(&Pos,&Axis);
342  if(LatDist > Nano[nNano].Rad) continue;
343  double Dist = Pos.ProjOnAxis(&Axis);
344  Pos.Set(pPos(p,CNorm)-Nano[nNano].Pos[CNorm],CNorm);
345  int t1 = TypePlus ? 2 : pType(p);
346  //fprintf(Ciccia,"%lf %lf %lf %d\n",pPos(p,0),pPos(p,1),pPos(p,2),t1);
347  //fprintf(Ciccia,"%lf %lf %lf %d\n",Pos.Val(0)+Nano[nNano].Pos[0]-5.,Pos.Val(1)+Nano[nNano].Pos[1]+5.,Pos.Val(2)+Nano[nNano].Pos[2],t1);
348  if(pPos(p,CNorm) > pNanoPos(nNano,CNorm)) t += 1;
349  int v = (int)(Dist*InvEdge*NBin);
350  if(v < 0 || v >= NBin) continue;
351  int z = (int)(pPos(p,CNorm)*pInvEdge(CNorm)*NBin);
352  if(z < 0 || z >= NBin) continue;
353  Prof[v*6+t] += pPos(p,CNorm) - pCm(CNorm);
354  Count[v*6+t] += 1.;
355  Plot[t1][v*NBin+z] += 1.;
356  if(!TypePlus) Dens[v] += 1.;
357  }
358  }
359  }
360  //fclose(Ciccia);
361  printf("\n");
362  //normalize
363  for(int b=0;b<NBin;b++){
364  for(int t=0;t<6;t++){
365  Prof[b*6+t] /= Count[b*6+t] > 0. ? Count[b*6+t] : 1.;
366  }
367  }
368  double Norm = 1./(NFile[1]-NFile[0]);
369  for(int v=0;v<SQR(NBin);v++){
370  for(int t=0;t<3;t++){
371  Plot[t][v] *= Norm;
372  }
373  }
374  Matrice Mask1(5);
375  Mask1.FillGaussian(.5,3.);
376  Mask1.Print();
377  int NDim = 1;
378  int IfMinImConv = 1;
379  Mask1.ConvoluteMatrix(Dens,NBin,NDim,IfMinImConv);
380  //write profile
381  FILE *FSlab = fopen("SlabProfile.dat","w");
382  FILE *FDens = fopen("SlabDensProfile.dat","w");
383  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
384  for(int b=0;b<NBin;b++){
385  double x = b/(double)NBin*pEdge(Coord1);
386  //fprintf(FSlab,"%lf %lf %lf %lf %lf %lf %lf\n",x,Prof[b*6 ],Prof[b*6+1],Prof[b*6+2],Prof[b*6+3],Prof[b*6+4],Prof[b*6+5],Prof[b*6+6]);
387  fprintf(FSlab,"%lf %lf %lf %lf %lf \n",x,Prof[b*6 ],Prof[b*6+1],Prof[b*6+2],Prof[b*6+3],Prof[b*6+4]);
388  fprintf(FDens,"%lf %lf\n",x,Dens[b]*InvNFile);
389  }
390  fclose(FSlab);
391  fclose(FDens);
392  //smooth and write dens
393  Matrice Mask(5,5);
394  Mask.FillGaussian(.5,3.);
395  Mask.Print();
396  NDim = 2;
397  IfMinImConv = 1;
398  for(int t=0;t<3;t++){
399  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
400  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
401  }
402  FILE *FSlabProf = fopen("SlabDens.xzd","w");
403  double LatDim[3] = {pEdge(Coord1),pEdge(CNorm),30.};
404  PrintDens(FSlabProf,Plot,LatDim,NBin);
405  fclose(FSlabProf);
406  free(Prof);
407  free(Count);
408  free(Temp);
409  free(Dens);
410 }
414 void ElPoly::PrintDens(FILE *FileToWrite,double **Plot,double *LatDim,int NBin){
415  int link[4] = {0,0,0,0};
416  double InvNBin = 1./(double)NBin;
417  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",LatDim[0],LatDim[1],LatDim[2],NBin,ChooseDraw(EL_DENS));
418  double Dr = InvNBin*LatDim[0];
419  double Dz = InvNBin*LatDim[1];
420  for(int t=0,p=0,c=0;t<3;t++){
421  for(int vv=0;vv<NBin-1;vv++){
422  for(int v=0;v<NBin-1;v++){
423  link[0] = (v+0)*NBin + (vv+0);
424  if(Plot[t][link[0]] < .02) continue;
425  //fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf)}\n",p,c,t, v*InvNBin*pEdge(3),(vv)*InvNBin*(Border[1]-Border[0]),(Plot[(v*NBin+vv)*NType+t]));
426  //continue;
427  //------------defines-the-squares---------------------
428  link[1] = (v+0)*NBin + (vv+1);
429  link[2] = (v+1)*NBin + (vv+0);
430  link[3] = (v+1)*NBin + (vv+1);
431  int pRef = p;
432  for(int lx=0;lx<2;lx++){
433  for(int ly=0;ly<2;ly++){
434  int l = 2*lx+ly;
435  double NanoAdded = Plot[2][(link[l])];//+Plot[3][(link[l])];
436  double Phob = t == 0 ? Plot[0][(link[l])] : 0.;
437  double Phil = t == 1 ? Plot[1][(link[l])] : 0.;
438  double r = (v +lx)*Dr + .5*Dr;
439  double z = (vv+ly)*Dz + .5*Dz;
440  double dens = Plot[t][(link[l])];//+Plot[(v*NBin+vv)*NType+1]+Plot[(v*NBin+vv)*NType+2]);
441  int l1 = pRef + (p+1)%4;
442  int l2 = pRef + (p+2)%4;
443  int l3 = pRef + (p+3)%4;
444  fprintf(FileToWrite,"{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);
445  p++;
446  }
447  }
448  c++;
449  }
450  }
451  }
452 }
456 void ElPoly::CartDens(int NBin,int nNano){
457  int NType = 5;
458  double **Plot = (double **)calloc(NType,sizeof(double));
459  for(int t=0;t<NType;t++){
460  Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
461  }
462  double NBinInv = 1./(double)NBin;
463  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
464  double LatDim[3] = {pEdge(CLat1),pEdge(CNorm),30.};
465  char FName[20];
466  double Width = 2.;
467  for(int f=NFile[0];f<NFile[1];f++){
468  Processing(f);
469  //--------------open,-back-fold----------------------
470  if(OpenRisk(cFile[f],NBackFold)) return;
471  for(int p=0;p<pNPart();p++){
472  double Dx = fabs(pCm(CLat1) - pPos(p,CLat1));
473  double Dy = fabs(pCm(CLat2) - pPos(p,CLat2));
474  if(Dy > Width) continue;
475  int vx = (int)(2.*Dx*pInvEdge(CLat1)*NBin);
476  if(vx<0 || vx>=NBin) continue;
477  int vz = (int)(pPos(p,CNorm)*pInvEdge(CNorm)*NBin);
478  if(vz<0 || vz>=NBin) continue;
479  int t = pType(p);
480  //if((p%Block[0].Asym)!=0)continue;
481  // if(VAR_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_INNER)) t = 0;
482  // else if(VAR_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_OUTER)) t = 1;
483  // else continue;
484  Plot[t][vx*NBin+vz] += 1.;
485  }
486  }
487  Matrice Mask(5,5);
488  Mask.FillGaussian(.5,3.);
489  Mask.Print();
490  for(int v=0;v<NBin*NBin;v++) for(int t=0;t<NType;t++) Plot[t][v] *= InvNFile;
491  int NDim = 2;
492  int IfMinImConv = 1;
493  for(int t=0;t<NType;t++){
494  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
495  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
496  }
497  sprintf(FName,"CartDens.dat");
498  FILE *FileToWrite = fopen(FName,"w");
499  PrintDens(FileToWrite,Plot,LatDim,NBin);
500  int NAvx = 4;
501  int NAvz = 4;
502  double NAv = 1./(double)(NAvx*NAvz);
503  double Avx = 0.;
504  double Avz = 0.;
505  double Av = 0.;
506  for(int vx=0;vx<NBin*NBin;vx++) Plot[3][vx] = 1.;
507  for(int t=0;t<2;t++){
508  int p = 0;
509  sprintf(FName,"LineLayer%d.dat",t);
510  FILE *FLine = fopen(FName,"w");
511  if(1==1){
512  fprintf(FLine,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(0),pEdge(1),pEdge(2),NBin);
513  for(int vx=0;vx<NBin;vx++){
514  for(int vz=0;vz<NBin;vz++){
515  fprintf(FLine,"{t[%d %d %d] x(%lf %lf %lf)}\n",p++,0,t,(double)(vx*pEdge(CLat1)/(double)NBin),(double)(vz*pEdge(CNorm)/(double)NBin),Plot[t][vx*NBin+vz],t);
516  }
517  }
518  }
519  for(int vx=0;vx<NBin;vx+=NAvx){
520  for(int vz=0;vz<NBin;vz+=NAvz){
521  for(int vvx=0;vvx<NAvx;vvx++){
522  for(int vvz=0;vvz<NAvz;vvz++){
523  int vTot = (vx+vvx)*NBin + vz+vvz;
524  if(vTot >= SQR(NBin)) break;
525  Plot[3][vTot] *= Plot[t][vTot];
526  double Weight = Plot[t][vTot];
527  if(Weight < 1.0) continue;
528  Avx += (double)(vx*pEdge(CLat1)/(double)NBin)*Weight;
529  Avz += (double)(vz*pEdge(CNorm)/(double)NBin)*Weight;
530  Av += Weight;
531  }
532  }
533  //fprintf(FLine,"%lf %lf %lf %d\n",Avx/Av,Avz/Av,Av/(double)NAverage,2);
534  fprintf(FLine,"{t[%d %d %d] x(%lf %lf %lf) l[%d]}\n",p++,0,2,Avx/Av,Avz/Av,Av*NAv,p+1);
535  Avx = 0.;
536  Avz = 0.;
537  Av = 0.;
538  }
539  }
540  fclose(FLine);
541  }
542  sprintf(FName,"LineLayer2.dat");
543  FILE *FLine = fopen(FName,"w");
544  for(int vx=0;vx<NBin;vx++){
545  for(int vz=0;vz<NBin;vz++){
546  double Weight = Plot[3][vx*NBin+vz];
547  if(Weight < 1.1) continue;
548  fprintf(FLine,"%lf %lf %lf\n",(double)(vx*pEdge(CLat1)/(double)NBin),(double)(vz*pEdge(CNorm)/(double)NBin),Weight);
549  }
550  }
551  for(int t=0;t<NType;t++) free(Plot[t]);
552  free(Plot);
553 }
557 void ElPoly::RadDistrF(int NBin,int How,int nNano){
558  if(nNano > pNNano()){
559  printf("The specified nanoparticle doesn't exist\n");
560  return ;
561  }
562  int NType = 5;
563  double *Dens = (double *)calloc(NBin,sizeof(double));
564  double **Plot = (double **)calloc(NType,sizeof(double));
565  for(int t=0;t<NType;t++){
566  Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
567  }
568  //--------The-first-file-defines-the-reference-distance--------
569  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
570  double RefPos[3] = {pCm(CLat1),pCm(CLat2),pCm(CNorm)};
571  double LatDim[3] = {pEdge(3),pEdge(CNorm),30.};
572  if(How >= 6){
573  LatDim[0] = pEdge(CLat1);
574  LatDim[1] = pEdge(CLat2);
575  }
576  double InvNBin = 1./(double)NBin;
577  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
578  double OldPos[3] = {pNanoPos(nNano,0),pNanoPos(nNano,1),pNanoPos(nNano,2)};
579  for(int f=NFile[0];f<NFile[1];f++){
580  Processing(f);
581  //--------------open,-back-fold----------------------
582  if(OpenRisk(cFile[f],NBackFold)) return;
583  BfDefChain();
584  if(How == 5) PorePos();
585  RefPos[CLat1] = pNanoPos(nNano,CLat1);
586  RefPos[CLat2] = pNanoPos(nNano,CLat2);
587  RefPos[CNorm] = pNanoPos(nNano,CNorm);
588  if(How == 2){
589  RefPos[CLat1] = pCm(CLat1);
590  RefPos[CLat2] = pCm(CLat2);
591  RefPos[CNorm] = pCm(CNorm);
592  }
593  if(How == 0){
594  for(int d=0;d<3;d++){
595  RefPos[d] = pCm(d);
596  }
597  }
598  if(How < 6){
599  //SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
600  //----------define-the-distance,-fill-Plot---------------
601  double TangInv = ABS(Nano[nNano].Axis[CNorm]) > 0. ? sqrt(QUAD(Nano[nNano].Axis[CLat1])+QUAD(Nano[nNano].Axis[CLat2]))/Nano[nNano].Axis[CNorm] : 0.;
602  if(How == 0) TangInv = 0.;
603  Vettore Rotated(3);
604  Vettore Dist(3);
605  Vettore PosRel(3);
606  Vettore Zed(3);
607  Zed.Set(0.,CLat1);Zed.Set(0.,CLat2);Zed.Set(1.,CNorm);
608  //Vettore NanoAxis(1.,0.,0.);
609  Vettore NanoAxis(Nano[nNano].Axis[0],Nano[nNano].Axis[1],Nano[nNano].Axis[2]);
610  if(How == 0){
611  NanoAxis.Set(0.,CLat1);
612  NanoAxis.Set(0.,CLat2);
613  NanoAxis.Set(1.,CNorm);
614  }
615  Vettore RotAxis(3);
616  RotAxis.VetV(&NanoAxis,&Zed);
617  double Angle = Zed.Angle(&NanoAxis,&Zed);
618  Matrice M(RotAxis.x,Angle,3);
619  int vzNano = (int)(pNanoPos(nNano,CNorm)*InvNBin*pInvEdge(CNorm));
620  for(int b=0;b<pNBlock();b++){
621  int IfAdded = 0;
622  if(!strncmp(Block[b].Name,"PEP",3) ) IfAdded = 2;
623  if(!strcmp(Block[b].Name,"STUFFING") ) continue;
624  if(!strcmp(Block[b].Name,"WATER") ) continue;
625  if(!strcmp(Block[b].Name,"CHOL") ) IfAdded = 1;
626  if(!strcmp(Block[b].Name,"OIL") ) continue;//IfAdded = 1;
627  if(!strcmp(Block[b].Name,"ADDED") ) continue;//IfAdded = 1;
628  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
629  // double RadCh = SQR(Ch[Pm[pPos(.CId,0)-Nano[nNano].PosBf[0])+SQR(Ch[Pm[pPos(.CId,1)-Nano[nNano].PosBf[1]);
630  // if( RadCh < SQR(Nano[nNano].Rad)){continue;}
631  //finding the position on the axis of the nanoparticle
632 
633  int c = pChain(p);
634  int p1 = c*pNPCh(b);
635  //if(fabs(Pm[p1].Pos[CNorm] - Nano[0].Pos[CNorm]) < 2.5) continue;
636  if(fabs(Pm[p1].Pos[CNorm] - Nano[0].Pos[CNorm]) < 3.5) continue;
637 
638  // if(!strcmp(Block[b].Name,"TT0"))
639  // if(VAR_IF_TYPE(Ch[c].Type,CHAIN_INNER)) continue;
640  // if(!strcmp(Block[b].Name,"TT1"))
641  // if(VAR_IF_TYPE(Ch[c].Type,CHAIN_OUTER)) continue;
642  // if(VAR_IF_TYPE(Ch[c].Type,CHAIN_OUTER)) continue;
643 
644 
645  for(int d=0;d<3;d++){
646  double Dist = pPos(p,d) - RefPos[d];
647  Dist -= floor( Dist*pInvEdge(d) + .5)*pEdge(d);
648  PosRel.Set(Dist,d);
649  }
650  double RadDist = fabs(Dist.PerpTo3(&PosRel,&NanoAxis));
651  // double RadDist = hypot(PosRel.Val(CLat1),PosRel.Val(CLat2));
652  //correspondent bin on the radial distance
653  int vr = (int)(NBin*RadDist*pInvEdge(3));
654  //Projecting the position on the axis and rotating to the normal
655  PosRel.ProjOnAxis(&NanoAxis);
656  M.Mult(PosRel,Rotated);
657  //double PosNorm = Rotated[CNorm] + .5*pEdge(CNorm) - floor(PosNorm*pInvEdge(CNorm))*pEdge(CNorm);
658  double PosNorm = PosRel.Val(CNorm) + .5*pEdge(CNorm);
659  int vz = (int)(NBin*(PosNorm)*pInvEdge(CNorm));
660  if( vz < 0 || vz >= NBin){
661  printf("Axis %lf %lf %lf \n",Nano[nNano].Axis[0],Nano[nNano].Axis[1],Nano[nNano].Axis[2]);
662  printf("Normal %lf > %lf %d > %d\n",PosNorm,pEdge(CNorm),vz,NBin);
663  continue;
664  }
665  if( vr < 0 || vr >= NBin){
666  continue;
667  }
668  //shifting wrt the inclination of the peptide
669  // v += (int)( (vv-NBin/2)*TangInv);
670  // if( v < 0) v += NBin;
671  // if( v >= NBin) v -= NBin;
672  int t = pType(p);
673  if(IfAdded==1){
674  t = 2;
675  Plot[0][vr*NBin+vz] += 1.;
676  }
677  if(IfAdded==2) t = 3;
678  else Dens[vr] += 1.;
679  Plot[t][vr*NBin+vz] += 1.;
680  }
681  }
682  }
683  else if(How >= 6){
684  for(int b=0;b<pNBlock();b++){
685  int IfAdded = 0;
686  if(!strncmp(Block[b].Name,"PEP",3) ) continue;//IfAdded = 1;
687  //if(!strncmp(Block[b].Name,"ADDED",3) ) IfAdded = 1;
688  if(!strcmp(Block[b].Name,"STUFFING") ) continue;
689  if(!strcmp(Block[b].Name,"WATER") ) continue;
690  if(!strcmp(Block[b].Name,"CHOL") ) IfAdded = 1;
691  if(!strcmp(Block[b].Name,"OIL") ) IfAdded = 1;
692  if(!strcmp(Block[b].Name,"ADDED") ) IfAdded = 1;
693  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
694  double Dist = pPos(p,CLat1) - RefPos[CLat1];
695  Dist -= floor(Dist*pInvEdge(CLat1)+.5)*pEdge(CLat1) - .5*pEdge(CLat1);
696  int vx = (int)(NBin*Dist*pInvEdge(CLat1));
697  Dist = pPos(p,CLat2) - RefPos[CLat2];
698  Dist -= floor(Dist*pInvEdge(CLat2)+.5)*pEdge(CLat2) - .5*pEdge(CLat2);
699  int vy = (int)(NBin*Dist*pInvEdge(CLat2));
700  int t = pType(p);
701  if( IfAdded ) Plot[2][vx*NBin+vy] += 1.;
702  else Plot[t][(vx*NBin+vy)] += 1.;
703  }
704  }
705  }
706  }
707 #ifdef OMPI_MPI_H
708  MPI_Allreduce(MPI_IN_PLACE,Plot, NType*NBin*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
709  int Rank=0;
710  MPI_Comm_rank(Proc->CommGrid, &Rank);
711  if(Rank==0){
712 #endif
713  printf("\n");
714  //------------Normalize----------------------
715  if(How < 6){
716  double *VolContr = (double *)calloc(NBin,sizeof(double));
717  VolumeCircSlab(VolContr,NBin);
718  for(int v=0;v<NBin;v++){
719  Dens[v] /= VolContr[v]*(NFile[1] - NFile[0]);
720  for(int vv=0;vv<NBin;vv++){
721  for(int t=0,n=0;t<NType;t++){
722  Plot[t][(v*NBin+vv)] /= VolContr[v]*(NFile[1]-NFile[0]);
723  }
724  }
725  }
726  free(VolContr);
727  }
728  if(How >= 6){
729  for(int v=0;v<NBin;v++){
730  Dens[v] /= (NFile[1] - NFile[0]);
731  for(int vv=0;vv<NBin;vv++){
732  for(int t=0,n=0;t<NType;t++){
733  Plot[t][(v*NBin+vv)] *= InvNFile;
734  }
735  }
736  }
737  }
738  Matrice Mask(5,5);
739  Mask.FillGaussian(.5,3.);
740  Mask.Print();
741  int NDim = 2;
742  int IfMinImConv = 0;
743  if(How == 6) IfMinImConv = 1;
744  for(int t=0;t<NType;t++){
745  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
746  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
747  }
748  //--------------------writes-the-files---------------------
749  char FileName[256];
750  char FileName2[256];
751  if(How == 0){//Cm Reference
752  sprintf(FileName,"RadDistrCmR%.1fH%.1fh%.1f",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
753  }
754  else if(How==1 || How == 4 || How == 5){//Nano reference
755  sprintf(FileName,"RadDistrNanoR%.1fH%.1fh%.1f",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
756  }
757  else if(How == 2){//Cm Nano Reference
758  sprintf(FileName,"RadDistrCmNR%.1fH%.1fh%.1f",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
759  }
760  else if(How == 6){//Cm Nano Reference
761  sprintf(FileName,"CartDistrNanoR%.1fH%.1fh%.1f",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
762  }
763  sprintf(FileName2,"%s.rzd",FileName);
764  FILE *FileToWrite = fopen(FileName2,"w");
765  PrintDens(FileToWrite,Plot,LatDim,NBin);
766  fclose(FileToWrite);
767  sprintf(FileName,"DensProfR%.1fH%.1fh%.1f.dat",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
768  FILE *FToWrite = fopen(FileName,"w");
769  for(int v=0;v<NBin;v++){
770  double r = v*InvNBin*pEdge(3);
771  fprintf(FToWrite,"%lf %lf\n",r,Dens[v]);
772  }
773 #ifdef OMPI_MPI_H
774  }
775 #endif
776  for(int t=0;t<NType;t++)
777  free(Plot[t]);
778  free(Plot);
779 }
783 void ElPoly::BondDistr(int NSample){
784  int NType = 3;
785  double **Plot = (double **)calloc(NType,sizeof(double));
786  double **Count = (double **)calloc(NType,sizeof(double));
787  for(int t=0;t<3;t++){
788  Plot[t] = (double *)calloc(SQR(NSample),sizeof(double));
789  Count[t] = (double *)calloc(SQR(NSample),sizeof(double));
790  }
791  double InvNSample = 1./(double)NSample;
792  double Cm[4] = {0.,0.,0.,0.};
793  double Dist[4] = {0.,0.,0.,0.};
794  double Bond[4] = {0.,0.,0.,0.};
795  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
796  for(int f=NFile[0];f<NFile[1];f++){
797  Processing(f);
798  OpenRisk(cFile[f],BF_CHAIN);
799  for(int b=0;b<pNBlock();b++){
800  int IfAdded = 0;
801  if(!strncmp(Block[b].Name,"PEP",3) ) continue;// IfAdded = 1;
802  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx-1;p++){
803  if(pChain(p) != pChain(p+1)) continue;
804  Dist[3] = 0.;
805  Bond[3] = 0.;
806  for(int d=0;d<3;d++){
807  Cm[d] = .5*(pPos(p,d) + pPos(p+1,d));
808  Bond[d] = Pm[p].Pos[d] - Pm[p+1].Pos[d];
809  Bond[3] += SQR(Bond[d]);
810  Dist[d] = Cm[d] - Nano[0].Pos[d];
811  Dist[d] -= floor( Dist[d]*pInvEdge(d) + .5)*pEdge(d);
812  }
813  double DistR = sqrt(SQR(Dist[CLat1])+SQR(Dist[CLat2]));
814  int t = pType(p);
815  if(t != pType(p+1)) continue;
816  if(IfAdded) t = 2;
817  double z = .5*pEdge(CNorm) + Dist[CNorm];
818  double r = DistR;
819  int vz = (int)(z*pInvEdge(CNorm)*NSample);
820  if(vz < 0 || vz >= NSample) continue;
821  int vr = (int)(r*pInvEdge(3)*NSample);
822  if(vr < 0 || vr >= NSample) continue;
823  Plot[t][vr*NSample+vz] += sqrt(Bond[3]);
824  Count[t][vr*NSample+vz] += 1.;
825  }
826  }
827  }
828  // for(int t=0;t<NType;t++){
829  // for(int v=0;v<SQR(NSample);v++){
830  // Plot[t][v] /= Count[t][v] > 0. ? 1./Count[t][v] : 1.;
831  // }
832  // }
833  double *VolContr = (double *)calloc(NSample,sizeof(double));
834  VolumeCircSlab(VolContr,NSample);
835  for(int v=0;v<NSample;v++){
836  for(int vv=0;vv<NSample;vv++){
837  for(int t=0,n=0;t<NType;t++){
838  Plot[t][(v*NSample+vv)] /= VolContr[v]*(NFile[1]-NFile[0]);
839  //double Norm = Count[t][v*NSample+vv] > 1. ? 1./Count[t][v*NSample+vv] : 1.;
840  //Plot[t][(v*NSample+vv)] *= Norm;
841  }
842  }
843  }
844  free(VolContr);
845  Matrice Mask(5,5);
846  Mask.FillGaussian(.5,3.);
847  Mask.Print();
848  int NDim = 2;
849  int IfMinImConv = 0;
850  for(int t=0;t<NType;t++){
851  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
852  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
853  }
854  char FileName2[60];
855  sprintf(FileName2,"BondDistr.rzd");
856  FILE *FileToWrite = fopen(FileName2,"w");
857  double LatDim[3] = {pEdge(3),pEdge(CNorm),15.};
858  PrintDens(FileToWrite,Plot,LatDim,NSample);
859  fclose(FileToWrite);
860  for(int t=0;t<NType;t++){
861  free(Plot[t]);
862  free(Count[t]);
863  }
864  free(Plot);
865  free(Count);
866 }
870 void ElPoly::SplayDistr(int NSample){
871  int NType = 3;
872  double **Plot = (double **)calloc(NType,sizeof(double));
873  double **Count = (double **)calloc(NType,sizeof(double));
874  for(int t=0;t<3;t++){
875  Plot[t] = (double *)calloc(SQR(NSample),sizeof(double));
876  Count[t] = (double *)calloc(SQR(NSample),sizeof(double));
877  }
878  double InvNSample = 1./(double)NSample;
879  double Cm[4] = {0.,0.,0.,0.};
880  double CmP[6] = {0.,0.,0.,0.,0.,0.};
881  double E2EP[6] = {0.,0.,0.,0.,0.,0.};
882  double GyrP[2] = {0.,0.};
883  double DistBA[4] = {0.,0.,0.,0.};
884  double DistCB[4] = {0.,0.,0.,0.};
885  int IfMartini = 0;
886  int HalfLim = (int)(Block[0].Asym*.5);
887  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
888  double Direct[6] = {0.,0.,0.,0.,0.,0.};
889  //int pLimit[4] = {4,8,8,13};//Martini
890  //double TailLength[2] = {1./4.,1./5.};//Martini
891  int pLimit[4] = {0,4,4,9};//DFT
892  double TailLength[2] = {1./4.,1./5.};//DFT
893  Vettore Ax0(1.,0.,0.);
894  Vettore Ax2(0.,0.,1.);
895  for(int f=NFile[0];f<NFile[1];f++){
896  Processing(f);
897  OpenRisk(cFile[f],BF_CHAIN);
898  for(int b=0;b<pNBlock();b++){
899  double *xc = (double *)calloc(pNPCh(b),sizeof(double));
900  double *yc = (double *)calloc(pNPCh(b),sizeof(double));
901  double *zc = (double *)calloc(pNPCh(b),sizeof(double));
902  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
903  if(!strncmp(Block[b].Name,"ADDED",3) ) continue;
904  for(int c=0;c<Block[b].NChain;c++){
905  for(int i=0;i<2;i++){
906  GyrP[i] = 0.;
907  int p1 = Block[b].InitIdx + c*Block[b].NPCh + pLimit[i*2];
908  int p2 = Block[b].InitIdx + c*Block[b].NPCh + pLimit[i*2+1];
909  int NBead = 0;
910  for(int d=0;d<3;d++){
911  CmP[i*3+d] = 0.;
912  E2EP[i*3+d] = Pm[p1].Pos[d] - Pm[p2].Pos[d];
913  }
914  for(int p=p1;p<p2;p++){
915  int ppc = p-p1;
916  xc[ppc] = Pm[p].Pos[0];
917  yc[ppc] = Pm[p].Pos[1];
918  zc[ppc] = Pm[p].Pos[2];
919  NBead += 1;
920  for(int d=0;d<3;d++){
921  CmP[i*3+d] += Pm[p].Pos[d];
922  }
923  }
924  for(int d=0;d<3;d++){
925  CmP[i*3+d] /= (double) NBead;
926  }
927  for(int p=p1;p<p2;p++){
928  for(int d=0;d<3;d++){
929  GyrP[i] += SQR(Pm[p].Pos[d]-CmP[i*3+d]);
930  }
931  }
932  RETTA rzx = Mat->InterRett(zc,xc,NBead);
933  RETTA rzy = Mat->InterRett(zc,yc,NBead);
934  double x1 = Pm[p1].Pos[2]*rzx.m + rzx.q;
935  double x2 = Pm[p2].Pos[2]*rzx.m + rzx.q;
936  double y1 = Pm[p1].Pos[2]*rzy.m + rzy.q;
937  double y2 = Pm[p2].Pos[2]*rzy.m + rzy.q;
938  Direct[2*i+2] = Pm[p1].Pos[2] - Pm[p2].Pos[2];
939  Direct[2*i+1] = y1 - y2;
940  Direct[2*i+0] = x1 - x2;
941  }
942  Vettore Dir1(Direct[0],Direct[1],Direct[2]);
943  Vettore Dir2(Direct[3],Direct[4],Direct[5]);
944  double Angle = Ax0.Angle(&Dir1,&Dir2);
945  if(isnan(Angle)) continue;
946  //int pA = Block[b].InitIdx + c*Block[b].NPCh + 3;//Martini
947  int pA = Block[b].InitIdx + c*Block[b].NPCh + 3;//DFT
948  for(int d=0;d<3;d++){
949  Cm[d] = Nano[0].Pos[d] - Pm[pA].Pos[d];
950  CmP[0+d] = CmP[0+d] - Nano[0].Pos[d];
951  CmP[0+d] -= floor( CmP[0+d]*pInvEdge(d) + .5)*pEdge(d);
952  CmP[3+d] = CmP[3+d] - Nano[0].Pos[d];
953  CmP[3+d] -= floor( CmP[3+d]*pInvEdge(d) + .5)*pEdge(d);
954  }
955  double E2E1 = sqrt(SQR(E2EP[0+0])+SQR(E2EP[0+1])+SQR(E2EP[0+2]));
956  double E2E2 = sqrt(SQR(E2EP[3+0])+SQR(E2EP[3+1])+SQR(E2EP[3+2]));
957  double DistR = sqrt(SQR(Cm[CLat1])+SQR(Cm[CLat2]));
958  double z = .5*pEdge(CNorm) + Cm[CNorm];
959  double r = DistR;
960  int vz = (int)(z*pInvEdge(CNorm)*NSample);
961  if(vz < 0 || vz >= NSample) continue;
962  int vr = (int)(r*pInvEdge(3)*NSample);
963  if(vr < 0 || vr >= NSample) continue;
964  Plot[0][vr*NSample+vz] += Angle;
965  Count[0][vr*NSample+vz] += 1.;
966  r = sqrt(SQR(CmP[0+CLat1])+SQR(CmP[0+CLat2]));
967  z = .5*pEdge(CNorm) + CmP[0+CNorm];
968  vz = (int)(z*pInvEdge(CNorm)*NSample);
969  if(vz < 0 || vz >= NSample) continue;
970  vr = (int)(r*pInvEdge(3)*NSample);
971  if(vr < 0 || vr >= NSample) continue;
972  // Plot[1][vr*NSample+vz] += Dir1.Norm()*TailLength[0];
973  Plot[1][vr*NSample+vz] += E2E1*TailLength[0];
974  Count[1][vr*NSample+vz] += 1.;
975  Plot[2][vr*NSample+vz] += sqrt(GyrP[0]*TailLength[0]);
976  Count[2][vr*NSample+vz] += 1.;
977  r = sqrt(SQR(CmP[3+CLat1])+SQR(CmP[3+CLat2]));
978  z = .5*pEdge(CNorm) + CmP[3+CNorm];
979  vz = (int)(z*pInvEdge(CNorm)*NSample);
980  if(vz < 0 || vz >= NSample) continue;
981  vr = (int)(r*pInvEdge(3)*NSample);
982  if(vr < 0 || vr >= NSample) continue;
983  //Plot[1][vr*NSample+vz] += Dir2.Norm()*TailLength[1];
984  Plot[1][vr*NSample+vz] += E2E2*TailLength[1];
985  Count[1][vr*NSample+vz] += 1.;
986  Plot[2][vr*NSample+vz] += sqrt(GyrP[1]*TailLength[1]);
987  Count[2][vr*NSample+vz] += 1.;
988  }
989  free(xc);
990  free(yc);
991  free(zc);
992  }
993  }
994  double *VolContr = (double *)calloc(NSample,sizeof(double));
995  VolumeCircSlab(VolContr,NSample);
996  for(int v=0;v<NSample;v++){
997  for(int vv=0;vv<NSample;vv++){
998  for(int t=0,n=0;t<NType;t++){
999  Plot[t][(v*NSample+vv)] /= VolContr[v]*(NFile[1]-NFile[0]);
1000  //double Norm = Count[t][v*NSample+vv] > 1. ? 1./Count[t][v*NSample+vv] : 1.;
1001  //Plot[t][(v*NSample+vv)] *= Norm;
1002  }
1003  }
1004  }
1005  free(VolContr);
1006  Matrice Mask(5,5);
1007  Mask.FillGaussian(.5,3.);
1008  Mask.Print();
1009  int NDim = 2;
1010  int IfMinImConv = 0;
1011  for(int t=0;t<NType;t++){
1012  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
1013  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
1014  }
1015  char FileName2[60];
1016  sprintf(FileName2,"SplayDistr.rzd");
1017  FILE *FileToWrite = fopen(FileName2,"w");
1018  double LatDim[3] = {pEdge(3),pEdge(CNorm),15.};
1019  PrintDens(FileToWrite,Plot,LatDim,NSample);
1020  fclose(FileToWrite);
1021  for(int t=0;t<NType;t++){
1022  free(Plot[t]);
1023  free(Count[t]);
1024  }
1025  free(Plot);
1026  free(Count);
1027 }
1031 void ElPoly::RadDens2Thick(int NBin){
1032  double *ThickProf = (double *)calloc(3*NBin,sizeof(double));
1033  double *Count = (double *)calloc(3*NBin,sizeof(double));
1034  double **Plot = (double **)calloc(3,sizeof(double));
1035  for(int t=0;t<3;t++){
1036  Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
1037  }
1038  LoadDensFile(Plot,NBin);
1039  double Cm[2] = {0.,0.};
1040  double CmCount[2] = {0.,0.};
1041  double InvNBin = 1./(double)NBin;
1042  for(int sr=0;sr<NBin;sr++){
1043  for(int sz=0;sz<NBin/2;sz++){
1044  double Pos = (sz-NBin/2)*pEdge(1)*InvNBin;
1045  double Weight = Plot[0][sr*NBin+sz]*Plot[1][sr*NBin+sz];
1046  ThickProf[sr*3 ] += Pos*Weight;
1047  Count[sr*3 ] += Weight;
1048  }
1049  for(int sz=NBin/2;sz<NBin;sz++){
1050  double Pos = (sz-NBin/2)*pEdge(1)*InvNBin;
1051  double Weight = Plot[0][sr*NBin+sz]*Plot[1][sr*NBin+sz];
1052  ThickProf[sr*3+1] += Pos*Weight;
1053  Count[sr*3+1] += Weight;
1054  }
1055  for(int sz=0;sz<NBin;sz++){
1056  ThickProf[sr*3+2] += Plot[2][sr*NBin+sz];
1057  Count[sr*3+2] += 1.;
1058  }
1059  }
1060  //normalize
1061  double Max = 0.;
1062  for(int sr=0;sr<NBin;sr++){
1063  ThickProf[sr*3 ] /= Count[sr*3 ] > 0. ? Count[sr*3 ] : 1.;
1064  ThickProf[sr*3+1] /= Count[sr*3+1] > 0. ? Count[sr*3+1] : 1.;
1065  ThickProf[sr*3+2] /= Count[sr*3+2] > 0. ? Count[sr*3+2] : 1.;
1066  if(Max < ThickProf[sr*3+2]) Max = ThickProf[sr*3+2];
1067  }
1068  for(int sr=0;sr<NBin;sr++){
1069  if(ThickProf[sr*3+2] > .5*Max){
1070  ThickProf[sr*3 ] = 0.;
1071  ThickProf[sr*3+1] = 0.;
1072  }
1073  }
1074  char FName[256];
1075  sprintf(FName,"ThickDensProf.dat",cFile[NFile[0]]);
1076  FILE *F2Write = fopen(FName,"w");
1077  for(int sr=0;sr<NBin;sr++){
1078  double x = sr*pEdge(0)/(double)NBin;
1079  fprintf(F2Write,"%lf %lf %lf %lf %lf\n",x,.5*(-ThickProf[sr*3 ]+ThickProf[sr*3+1]),-ThickProf[sr*3 ],ThickProf[sr*3+1],ThickProf[sr*3+2]);
1080  }
1081  for(int t=0;t<3;t++)
1082  free(Plot[t]);
1083  free(Plot);
1084  free(ThickProf);
1085  free(Count);
1086 }
1091  double *ThickProf = (double *)calloc(NBin*NBin,sizeof(double));
1092  double *Count = (double *)calloc(NBin*NBin,sizeof(double));
1093  double Round = 0.001;
1094  for(int p=0;p<pNPart();p++){
1095  if(pType(p) != 1) continue;
1096  int sr = (int)((pPos(p,0)+Round)*pInvEdge(0)*NBin);
1097  int sz = (int)((pPos(p,1)+Round)*pInvEdge(1)*NBin);
1098  ThickProf[sr*NBin+sz] += pPos(p,2);
1099  Count[sr*NBin+sz] += 1.;
1100  }
1101  double Norm = 0.;
1102  for(int sr=0;sr<NBin;sr++){
1103  for(int sz=0;sz<NBin;sz++){
1104  if(Norm < ThickProf[sr*NBin+sz]) Norm = ThickProf[sr*NBin+sz];
1105  }
1106  }
1107  Norm = 1./Norm;
1108  char FName[256];
1109  sprintf(FName,"ThickDensProf2d.dat",cFile[NFile[0]]);
1110  FILE *F2Write = fopen(FName,"w");
1111  fprintf(F2Write,"%lf %lf %lf\n",0.,0.,0.);
1112  fprintf(F2Write,"%lf %lf %lf\n",pEdge(0),pEdge(1),1.);
1113  for(int sr=0;sr<NBin;sr++){
1114  double r = sr*pEdge(0)/(double)NBin;
1115  for(int sz=0;sz<NBin;sz++){
1116  double z = sz*pEdge(1)/(double)NBin;
1117  //if(Count[sr*2] <= 1. || Count[sr*2+1] <= 1.)continue;
1118  if(ThickProf[sr*NBin+sz]*Norm < .3) continue;
1119  fprintf(F2Write,"%lf %lf %lf\n",r,z,ThickProf[sr*NBin+sz]*Norm);
1120  }
1121  }
1122  fclose(F2Write);
1123  free(ThickProf);
1124 }
1128 void ElPoly::ThickFromDens(int NBin){
1129  int NType = 3;
1130  double Round = 0.001;
1131  double InvNBin = 1./(double)NBin;
1132  double *Plot = (double *)calloc(NType*NBin*NBin,sizeof(double));
1133  double *Line = (double *)calloc(2*NBin,sizeof(double));
1134  double *Count = (double *)calloc(2*NBin,sizeof(double));
1135  for(int p=0;p<pNPart();p++){
1136  int vr = (int)((pPos(p,CLat1)+Round)*pInvEdge(0)*NBin);
1137  if (vr<0 || vr>NBin) continue;
1138  int vz = (int)((pPos(p,CLat2)+Round)*pInvEdge(1)*NBin);
1139  if (vz<0 || vz>NBin) continue;
1140  //int t = Pm[p].Typ;
1141  for(int t=0;t<3;t++){
1142  Plot[(vr*NBin+vz)*NType+t] += pVel(p,t);
1143  }
1144  }
1145  for(int vr=0;vr<NBin;vr++){
1146  for(int vz=0;vz<NBin;vz++){
1147  double z = vz*InvNBin*pEdge(1);
1148  double Weight = Plot[(vr*NBin+vz)*NType+1]*Plot[(vr*NBin+vz)*NType+2];
1149  int t = 0;
1150  if(z>pCm(1)) t = 1;
1151  Line[vr*2+t] += z*Weight;
1152  Count[vr*2+t] += Weight;
1153  }
1154  }
1155  FILE *FLine = fopen("ThickProf.dat","w");
1156  for(int vr=0;vr<NBin;vr++){
1157  double r = vr*InvNBin*pEdge(0);
1158  Line[vr*2] /= Count[vr*2] > 0. ? Count[vr*2] : 1.;
1159  Line[vr*2+1] /= Count[vr*2+1] > 0. ? Count[vr*2+1] : 1.;
1160  fprintf(FLine,"%lf %lf %lf\n",r,Line[vr*2],Line[vr*2+1]);
1161  }
1162  fclose(FLine);
1163  free(Plot);
1164  free(Count);
1165  free(Line);
1166 }
1170 void ElPoly::SlabAngleProfs(int NBin,int NAngle,int Coord){
1171  NAngle = 5;
1172  double *ProfAngle = (double *)calloc(NAngle*NBin,sizeof(double));
1173  double *Count = (double *)calloc(NAngle*NBin,sizeof(double));
1174  double Ref[3] = {.5*pEdge(1),.5*pEdge(1),pEdge(2)};
1175  double BoxRad = sqrt( SQR(.5*pEdge(0)) + SQR(.5*pEdge(1)) ) ;
1176  double InvBoxRad = 1./BoxRad;
1177  int Coord1 = (Coord)%2;
1178  int Coord2 = (Coord+1)%2;
1179  FILE *Ciccia = fopen("SlabAngleSectors.dat","w");
1180  fprintf(Ciccia,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(0),pEdge(1),10.,NBin);
1181  for(int p=0,p1=0;p<pNPart();p++){
1182  if(pPos(p,CNorm) < 1.) continue;
1183  double RelDist[2] = {pPos(p,CLat1) - Ref[CLat1],pPos(p,CLat2) - Ref[CLat2]};
1184  double r = hypot(RelDist[0],RelDist[1]);
1185  int rb = (int)(r*InvBoxRad*NBin);
1186  if(rb < 0 || rb >= NBin) continue;
1187  if(Coord == 2) RelDist[CLat2] = - RelDist[CLat2];
1188  if(Coord == 3) {RelDist[CLat1] = - RelDist[CLat1];RelDist[CLat2] = - RelDist[CLat2];}
1189  double phi = fabs(atan2(RelDist[Coord1],RelDist[Coord2]));
1190  int phib = (int)(phi*NAngle/(M_PI*.5));
1191  if(phib < 0 || phib >= NAngle) continue;
1192  ProfAngle[rb*NAngle+phib] += pPos(p,CNorm);
1193  Count[rb*NAngle+phib] += 1.;
1194  fprintf(Ciccia,"{t[%d 0 %d] x(%lf %lf %lf)}\n",p1++,phib,pPos(p,0),pPos(p,1),pPos(p,2));
1195  }
1196  fclose(Ciccia);
1197  FILE *FProf = fopen("SlabRadProfs.dat","w");
1198  for(int rb=0;rb<NBin;rb++){
1199  double r = rb*BoxRad/(double)NBin;
1200  fprintf(FProf,"%lf ",r);
1201  for(int phib=0;phib<NAngle;phib++){
1202  double Norm = Count[rb*NAngle+phib] > 0. ? 1./Count[rb*NAngle+phib] : 1.;
1203  fprintf(FProf,"%lf ",ProfAngle[rb*NAngle+phib]*Norm);
1204  }
1205  fprintf(FProf,"\n");
1206  }
1207  free(ProfAngle);
1208  free(Count);
1209 }
1213 int ElPoly::CoreF(int NSample,int How){
1214  double Cm=0.;
1215  double Volume=0.;
1216  int NVolume = NSample*NSample*NSample;
1217  int NType = 3;//Gen->NType;
1218  double *Plot = (double *) calloc(NVolume*NType,sizeof(double));
1219  double Border[3][2];
1220  Border[0][0]= 0.;
1221  Border[0][1]= pEdge(0);
1222  Border[1][0]= 0.;
1223  Border[1][1]= pEdge(1);
1224  Border[2][0]= 0.;
1225  Border[2][1]= pEdge(2);
1226  for(int f=NFile[0];f<NFile[1];f++){
1227  Processing(f);
1228  if(How == 0){
1229  if(OpenRisk(cFile[f],BF_CHAIN))return 1;
1230  }
1231  else if(How == 1){
1232  if(OpenRisk(cFile[f],BF_NANO))return 1;
1233  BfDefChain();
1234  }
1235  for(int b=0;b<pNBlock();b++){
1236  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
1237  //if(strcmp(Block[b].Name,"LIPID") ) continue;
1238  int vx = (int)(NSample*(pPos(p,0)-Border[0][0])/(Border[0][1]-Border[0][0]));
1239  if( vx < 0 || vx >= NSample) continue;
1240  int vy = (int)(NSample*(pPos(p,1)-Border[1][0])/(Border[1][1]-Border[1][0]));
1241  if( vy < 0 || vy >= NSample) continue;
1242  int vz = (int)(NSample*(pPos(p,2)-Border[2][0])/(Border[2][1]-Border[2][0]));
1243  if( vz < 0 || vz >= NSample) continue;
1244  int t = pType(p);
1245  if( CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED) ) t = 2;
1246  if( t < 0 || t > 3) continue;
1247  Plot[((vx*NSample+vy)*NSample+vz)*NType+t] += 1.;
1248  }
1249  Cm += pCm(2);
1250  }
1251  }
1252  printf("\n");
1253  Cm /= NFile[1] - NFile[0];
1254  FILE *FileToWrite = NULL;
1255  FileToWrite = fopen("Core.xvl","w");
1256  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(0),pEdge(1),pEdge(2),NSample,ChooseDraw(EL_COLOR));
1257  double *Norm = (double *)calloc(NType,sizeof(double));
1258  for(int t=0;t<NType;t++){
1259  for(int v=0;v<NVolume;v++){
1260  if(Norm[t] < Plot[v*NType+t])
1261  Norm[t] = Plot[v*NType+t];
1262  }
1263  Norm[t] = Norm[t] <= 0. ? 1. : Norm[t];
1264  }
1265  //for(int t=0;t<NType;t++){
1266  for(int vx=0;vx<NSample;vx++){
1267  for(int vy=0;vy<NSample;vy++){
1268  for(int vz=NSample-1;vz>0;vz--){
1269  double Dens0 = Plot[((vx*NSample+vy)*NSample+vz)*NType+2]/Norm[2];
1270  double Dens1 = Plot[((vx*NSample+vy)*NSample+vz)*NType+0]/Norm[0];
1271  double Dens2 = Plot[((vx*NSample+vy)*NSample+vz)*NType+1]/Norm[1];
1272  if(Dens0+Dens1+Dens2 <= 0.0)continue;
1273  fprintf(FileToWrite,"{x(%lf %lf %lf) ",
1274  Border[0][0]+vx/(double)NSample*(Border[0][1]-Border[0][0]),
1275  Border[1][0]+vy/(double)NSample*(Border[1][1] - Border[1][0]),
1276  Border[2][0]+vz/(double)NSample*(Border[2][1] - Border[2][0]));
1277  fprintf(FileToWrite,"v(%lf %lf %lf)}\n",Dens0,Dens1,Dens2);
1278  Volume += 1.;
1279  }
1280  }
1281  }
1282  //}
1283  Volume /= NSample*NSample*NSample;
1284  Volume *= (Border[0][1] - Border[0][0])*(Border[1][1] - Border[1][0])*(Border[2][1] - Border[2][0]);
1285  fprintf(FileToWrite,"#Volume %lf\n",Volume);
1286  fclose(FileToWrite);
1287  free(Plot);
1288  free(Norm);
1289  return 0;
1290 }
void RadDens2Thick2d(int NBin)
Calculate the thickness from the radial density profile.
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
void SysDef(char *cSystem)
Print a string with the system definitions.
Definition: VarData.cpp:167
int Type
Type of the chain (see list CHAIN_)
Definition: VarData.h:246
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.
void ThickFromDens(int NBin)
Calculate the thickness profile from the density plot.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double Normalize()
Normlizes a Vettore.
int pChain(int p)
Return the chain.
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
int BfDefChain()
Definition of the chain.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
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
void Mult(Matrice &A, Matrice &B)
Multiplication between two matrices.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double PerpTo3(Vettore *Pos, Vettore *Axis)
The vector perpendicolar in three dimension (faster)
double Val(int N)
Value of the N column.
double * x
Where the data are stored.
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
double m
y = m*x + q
int InitIdx
Initial particle position.
Definition: VarData.h:265
void RadDistrF(int NBin, int How, int nNano)
3d (rad,norma,dens) density profile with respect to a initial position radNano radCm radCmN ...
int Angle(int NBin)
Calculation of the contact angle angle Boh?
Definition: ElPolyEl.cpp:78
double Pos[3]
Position.
Definition: VarData.h:427
Linear interpolation.
double pNanoPos(int n, int d)
Return back folded nano position.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
void BondDistr(int NSample)
Distribution of the bond lengths.
int pNPCh()
Number of particle per chain.
double PorePos()
Find the position of the pore.
void SlabAngleProfs(int NBin, int NAngle, int Coord)
Radial profiles of the slab density between 0⁰ and 90⁰ subdivided in NAngle angles.
int CoreF(int NBin, int How)
Sampling of the three dimentional space.
int pNBlock()
Number of blocks.
double pchiN()
Incompatibility.
Definition: VarData.h:928
double ProjOnAxis(Vettore *a)
Projects along the axis.
int DensProf(int NBin, int NSample, int Coord)
Calculate the density profile dens.
void RadDens2Thick(int NBin)
Calculate the thickness from the radial density profile.
int pNNano()
Number of nanoparticles.
void CartDens(int NBin, int nNano)
3d (rad,norma,dens) density profile in cartesian coordinates
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
int NBackFold
Type of backfold.
Definition: ElPoly.h:439
void PrintDens(FILE *FileToWrite, double **Plot, double *LatDim, int NBin)
Print the density profile in the surfaces representation.
int NPCh
particles per chain
Definition: VarData.h:259
Matrice computes the algebric operations on matrices.
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int NChain
chains
Definition: VarData.h:263
int EndIdx
End particle position.
Definition: VarData.h:267
double pPos(int p, int d)
Return back folded position.
int pNChain()
Number of chain.
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
double Angle(Vettore *u, Vettore *v)
Computes the angle between two Vetttore.
int CLat1
lateral coordinate
Definition: VarData.h:1076
void DensProfNormalSlab(int NBin, int NSample, int Coord)
Calculate the density profile projecting the particles on the normal density.
void SplayDistr(int NSample)
Distribution of end to end distances.
double Axis[3]
Rotation axis.
Definition: VarData.h:435
int SetNNano(int Val)
Set NNano.
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
void SlabProf(int NBin, int nNano, int Coord1)
Thickness and density profile along Coord1.
void Print()
Print the entries.
void Set(double Val, int Col)
Set the N column.
double q
y = m*x + q
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
void VolumeCircSlab(double *VolContr, int NSample)
Fill an array of.
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
double pkappaN()
Incompressibility.
Definition: VarData.h:930
void PerpTo(Vettore *o)
The vector perpendicolar.
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.
RETTA InterRett(double *Px, double *Py, int NMass)
Linear interpolation.
void LoadDensFile(double **Plot, int NBin)
Load in the array Plot the density of the system.
Definition: VarDataEl.cpp:150