Allink  v0.1
ElPolyMeasure.cpp
1 #include "ElPoly.h"
2 
4  int NDim = 2;
5  NFileTot = NFile[1] - NFile[0];
6  int NChain = 0;
7  for(int b=0;b<pNBlock();b++){
8  if(!strcmp(Block[b].Name,"ADDED")) continue;
9  if(!strcmp(Block[b].Name,"CHOL")) continue;
10  NChain += Block[b].NChain;
11  }
12  double *FileAverage = (double *)calloc(NDim*NChain*NFileTot,sizeof(double));
13  double *ChainDiff = (double *)calloc(NFileTot,sizeof(*ChainDiff));
14  double *ChainOrigin = (double *)calloc(NFileTot,sizeof(*ChainDiff));
15  double *Count = (double *)calloc(NFileTot,sizeof(*Count));
16  double *Time = (double *)calloc(NFileTot,sizeof(*Time));
17  double CmInit[3] = {pCm(0),pCm(1),pCm(2)};
18  double InitTime = pTime();
19  for(int f=NFile[0];f<NFile[1];f++){
20  Processing(f);
21  if(OpenRisk(cFile[f],BF_NO))return 1;
22  BfDefChain();
23  Time[f] = pTime() - InitTime;
24  if(Time[f]<=0.) Time[f] = (double)f;
25  for(int c=0;c<NChain;c++){
26  int ff = f - NFile[0];
27  for(int d=0;d<NDim;d++){
28  FileAverage[(c*NFileTot+ff)*NDim+d] = Ch[c].Pos[d];
29  }
30  }
31  }
32  printf("\n");
33 #ifdef OMPI_MPI_H
34  MPI_Allreduce(MPI_IN_PLACE,FileAverage,NDim*NChain*NFileTot, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
35  int Rank=0;
36  MPI_Comm_rank(Proc->CommGrid, &Rank);
37  if(Rank==0){
38 #endif
39  for(int f=0;f<NFileTot;f++){
40  for(int c=0;c<NChain;c++){
41  for(int d=0;d<NDim;d++){
42  ChainOrigin[f] +=
43  SQR((FileAverage[(c*NFileTot)*NDim+d]-FileAverage[(c*NFileTot+f)*NDim+d]));//-(pCm(d)-CmInit[d]) );
44  }
45  }
46  }
47  for(int f=0;f<NFileTot;f++){
48  for(int ff = f; ff<NFileTot;ff++){
49  int fff = ff - f;
50  for(int c=0;c<NChain;c++){
51  for(int d=0;d<NDim;d++){
52  ChainDiff[fff] +=
53  SQR((FileAverage[(c*NFileTot+ff)*NDim+d]-FileAverage[(c*NFileTot+f)*NDim+d]));//-(pCm(d)-CmInit[d]) );
54  }
55  Count[fff] += 1.;
56  }
57  //printf("%d %d %d %lf\n",f,ff,fff,ChainDiff[fff]);
58  }
59  }
60  char FileName[60];
61  sprintf(FileName,"Diffusivity.dat");
62  FILE *FileToWrite = fopen(FileName,"w");
63  char cSystem[STRSIZE];
64  SysDef(cSystem);
65  fprintf(FileToWrite,"%s",cSystem);
66  fprintf(FileToWrite,"#Time ChainOrigin OriginSum DiffSum ChainDiff Count\n");
67  for(int f=0;f<NFileTot;f++){
68  double Inv = Count[f] > 0. ? Count[f] : 1.;
69  fprintf(FileToWrite,"%lf %.5g %.5g %.5g\n",Time[f],ChainOrigin[f]/(double)(NChain),ChainDiff[f]/Inv,ChainDiff[f]/(Inv*Time[f]) );
70  }
71 #ifdef OMPI_MPI_H
72  }
73 #endif
74  free(FileAverage);
75  free(ChainDiff);
76  free(ChainOrigin);
77  free(Count);
78  free(Time);
79  return 0;
80 }
81 void ElPoly::DiffSlab(int NSlab){
82  int NBin = NFile[1] - NFile[0];
83  double Volume1 = 0.;
84  double HeightCyl = pEdge(CNorm);
85  double LenMin = MIN(pEdge(CLat1),pEdge(CLat2));
86  double LenMax = MAX(pEdge(CLat1),pEdge(CLat2));
87  double BoxRad = .5*MIN(pEdge(CLat1),pEdge(CLat2));
88  double InvBoxRad = 1./BoxRad;
89  double *OldPos = (double *)calloc(3*pNPart(),sizeof(double));
90  double *Diff = (double *)calloc(NBin*pNPart(),sizeof(double));
91  double *Count = (double *)calloc(NSlab,sizeof(double));
92  int *PartOr = (int *)calloc(pNPart(),sizeof(int));
93  double dr[3];
94  for(int p = 0;p<pNPart();p++){
95  double dr2 = 0.;
96  for(int d=0;d<3;d++){
97  OldPos[3*p+d] = Pm[p].Pos[d];
98  dr[d] = Nano->Pos[d] - Pm[p].Pos[d];
99  dr[d] -= floor(dr[d]*pInvEdge(d))*pEdge(d);
100  }
101  dr2 = sqrt(SQR(dr[CLat1]) + SQR(dr[CLat2]));
102  int vr = (int)(dr2/BoxRad*NSlab);
103  PartOr[p] = vr;
104  if(vr < 0 || vr >= NSlab) continue;
105  Count[vr] += 1.;
106  }
107  char FName[60];
108  if(1==1){//particles position and slab circles
109  FILE *FRepr = fopen("SystemSlab.dat","w");
110  for(int p=0;p<pNPart();p++){
111  fprintf(FRepr,"%lf %lf %lf 0\n",Pm[p].Pos[CLat1],Pm[p].Pos[CLat2],.5);
112  }
113  for(int v=0;v<NSlab;v++){
114  double Dist = v*BoxRad/(double)NSlab;
115  for(int i=0;i<100;i++){
116  double ang = i*.01*2.*M_PI;
117  double x = Dist*cos(ang) + Nano->Pos[CLat1];
118  double y = Dist*sin(ang) + Nano->Pos[CLat2];
119  fprintf(FRepr,"%lf %lf %lf 2\n",x,y,.51);
120  }
121  }
122  fclose(FRepr);
123  double Norm = 1./(double)pNPart();
124  for(int v=0;v<NSlab;v++){
125  double Dist = v*BoxRad/(double)NSlab;
126  sprintf(FName,"MeanSqDispl%.3f.dat",Dist);
127  FILE *FWrite = fopen(FName,"w");
128  fclose(FWrite);
129  }
130  }
131  int NOut = 0;
132  for(int f=NFile[0];f<NFile[1];f++){
133  Processing(f);
134  if(OpenRisk(cFile[f],BF_PART));
135  double OutRatio = 100.*NOut/(double)(f*pNPart());
136  fprintf(stderr,"step %d)/%d fraction of number of beads out of the slab %lf %%\r",f,NBin,OutRatio);
137  double Rad2Mean = 0.;
138  for(int p = 0;p<pNPart();p++){
139  int vr1 = PartOr[p];
140  if(vr1 < 0 || vr1 >= NSlab) continue;
141  double dr2 = 0.;
142  for(int d=0;d<3;d++){
143  dr[d] = Nano->Pos[d] - Pm[p].Pos[d];
144  dr[d] -= floor(dr[d]*pInvEdge(d))*pEdge(d);
145  }
146  dr2 = sqrt(SQR(dr[CLat1]) + SQR(dr[CLat2]));
147  int vr = (int)(dr2*InvBoxRad*NSlab);
148  if(vr != vr1){
149  NOut++;
150  }
151  double Rad2 = 0.;
152  for(int d=0;d<2;d++){
153  double Dist = Pm[p].Pos[d] - OldPos[3*p+d];
154  Dist -= floor(Dist*pInvEdge(d))*pEdge(d);
155  Rad2 += SQR(Dist);
156  }
157  // vr or vr1?
158  Diff[f*NSlab+vr1] += Rad2;
159  Rad2Mean += Rad2;
160  }
161  if(1==1){//particle position and mean displacement
162  Rad2Mean = pNPart()/Rad2Mean;
163  sprintf(FName,"DisplPos%05d.dat",f);
164  FILE *FWrite = fopen(FName,"w");
165  fprintf(FWrite,"# l(%lf %lf %lf) d[part]\n",pEdge(CLat1),pEdge(CLat2),.5);
166  for(int p=0;p<pNPart();p++){
167  double Rad2 = 0.;
168  for(int d=0;d<2;d++){
169  double Dist = Pm[p].Pos[d] - OldPos[3*p+d];
170  Dist -= floor(Dist*pInvEdge(d))*pEdge(d);
171  Rad2 += SQR(Dist);
172  }
173  int vr1 = PartOr[p];
174  double Vel = Rad2*Rad2Mean;
175  fprintf(FWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)\n",Pm[p].Pos[CLat1],Pm[p].Pos[CLat2],Vel,Vel,Vel,Vel);
176  }
177  fclose(FWrite);
178  }
179  }
180  printf("\n");
181  double DisplMean = 0.;
182  //normalize
183  for(int v=0;v<NSlab;v++){
184  Count[v] = Count[v] > 0. ? 1./Count[v] : 1.;
185  for(int b=0;b<NBin;b++){
186  Diff[b*NSlab+v] *= Count[v];
187  }
188  }
189  for(int v=0;v<NSlab;v++){
190  DisplMean += Diff[(NBin-1)*NSlab+v] ;
191  double Dist = v*BoxRad/(double)NSlab;
192  sprintf(FName,"MeanSqDispl%.3f.dat",Dist);
193  FILE *FWrite = fopen(FName,"w");
194  for(int b=1;b<NBin;b++){
195  fprintf(FWrite,"%d %lf\n",b,Diff[b*NSlab+v]);
196  }
197  fclose(FWrite);
198  }
199  DisplMean /= (double) NSlab;
200  if(1==1){//slab mean displ
201  sprintf(FName,"SlabMeadDispl.dat");
202  FILE *FWrite = fopen(FName,"w");
203  fprintf(FWrite,"# l(%lf %lf %lf) d[part]\n",pEdge(CLat1),pEdge(CLat2),.5);
204  for(int sx=0;sx<100;sx++){
205  for(int sy=0;sy<100;sy++){
206  double x = sx*.01*pEdge(CLat1);
207  double y = sy*.01*pEdge(CLat2);
208  double Rad2 = sqrt(SQR(x-Nano->Pos[CLat1])+SQR(y-Nano->Pos[CLat2]));
209  int vr = (int)(Rad2*InvBoxRad*NSlab);
210  if(vr < 0 || vr >= NSlab) continue;
211  double Vel = Diff[(NBin-1)*NSlab+vr];
212  fprintf(FWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)\n",x,y,Vel,Vel,Vel,Vel);
213  }
214  }
215  fclose(FWrite);
216  }
217  free(OldPos);
218  free(Diff);
219  free(Count);
220  free(PartOr);
221 }
222 #include <Cubo.h>
223 typedef DdLinkedList DomDec;
224 void ElPoly::PairCorr(int NBin,int NDim){
225  int NAlloc = NBin;
226  double CutOff = 4.;
227  if(NDim == 2) NAlloc = NBin*NBin;
228  double *Prob = (double *)calloc(NAlloc,sizeof(double));
229  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
230  double Norm = 0.;
231  for(int f=NFile[0];f<NFile[1];f++){
232  Processing(f);
233  // DomDec *Pc = new DomDec(Edge,pNPart(),1.);
234  // for(int p=0;p<pNPart();p++) Pc->AddPart(p,Pm[p].Pos);
235  // int NeiList[27];
236  double DistRel[4];
237  for(int p1=0;p1<pNPart();p1++){
238  for(int p2=p1+1;p2<pNPart();p2++){
239  // int NNei = Pc->GetNei(Pm[p1].Pos,NeiList);
240  // for(int i=0;i<NNei;i++){
241  // int c1 = NeiList[i];
242  // for(Pc->SetCounters(c1);Pc->IfItCell(c1);Pc->IncrCurr(c1)){
243  // int p2 = Pc->ItCell(c1);
244  // if(p2 < p1) continue;
245  if(!TwoPartDist(p1,p2,DistRel,CutOff)) continue;
246  int r = (int)(DistRel[3]/CutOff*NBin);
247  if(r < 0 || r >= NBin) continue;
248  Prob[r] += 1.;
249  Norm += 1.;
250  }
251  }
252  // }
253  // delete Pc;
254  }
255  printf("\n");
256  Norm *= NFile[1] - NFile[0];
257  double rInv = 1.*CutOff/(double)NBin;
258  FILE *FPair = fopen("PairCorrelation.dat","w");
259  fprintf(FPair,"%lf %lf \n",0.,Prob[0]/Norm);
260  for(int r=1;r<NBin;r++){
261  double rr = r*CutOff/(double)NBin;
262  Prob[r] /= Norm*(CUBE(rr)-CUBE(rr-rInv));
263  fprintf(FPair,"%lf %lf \n",rr,Prob[r]);
264  }
265  fclose(FPair);
266  free(Prob);
267 }
268 int ElPoly::PairCorrelationF(int NBin,int How){
269  double **Plot;
270  Plot = (double **)calloc(NBin,sizeof(double));
271  for(int i=0;i<NBin;i++){
272  *(Plot+i) = (double *)calloc(NBin,sizeof(double));
273  }
274  double *dDensity;
275  dDensity = (double *)calloc(NBin,sizeof(double));
276  double *dDensity1;
277  dDensity1 = (double *)calloc(NBin,sizeof(double));
278  for(int f=NFile[0];f<NFile[1];f++){
279  Processing(f);
280  //printf("Opening: %s\n",cFile[f]);
281  if(OpenRisk(cFile[f],BF_CHAIN))return 1;
282  if(How == 0 || How == 1)
283  PairCorrelation(dDensity,NBin,How,NChType);
284  else if (How == 2)
285  PairCorrelationRound(Plot,NBin,NChType);
286  else if (How == 3)
287  PairCorrelationSquare(Plot,NBin,NChType);
288  else if (How == 4)
289  PairCorrelationPep(Plot,NBin,NChType);
290  for(int v=0;v<NBin;v++){
291  //printf("%d %lf\n",v,dDensity[v]);
292  dDensity1[v] += dDensity[v];
293  }
294  }
295  printf("\n");
296  if(How == 0 || How == 1){
297  FILE *FileToWrite = fopen("PairCorrelation.dat","w");
298  fprintf(FileToWrite,"# RadDist Density \n");
299  for(int v=0;v<NBin;v++){
300  dDensity1[v] /= NFile[1] - NFile[0];
301  fprintf(FileToWrite,"%lf %lf\n",(double)v/(double)NBin*pEdge(3),dDensity1[v]);
302  }
303  fclose(FileToWrite);
304  }
305  else {
306  double Max=0.;
307  double InvNBin=1./(double)NBin;
308  for(int v=0;v<NBin;v++){
309  for(int vv=0;vv<NBin;vv++){
310  // Plot[v][vv] /= DUE_PI*( QUAD((Gen->Edge[3]/InvNBin))*(QUAD((v+1)) - QUAD(v)) );
311  if(Max < Plot[v][vv])
312  Max = Plot[v][vv];
313  }
314  }
315  FILE *FileToWrite = fopen("PairCorrelation.xvl","w");
316  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",2.*pEdge(3),2.*pEdge(3),pEdge(3),NBin,ChooseDraw(EL_QUAD1));
317  int vRef = (int)(NBin/2.);
318  for(int v=0;v<NBin;v++){
319  for(int vv=0;vv<NBin;vv++){
320  if(Plot[v][vv] > 0.){
321  if(How == 2){
322  double PosX = pEdge(3) + InvNBin*(double)v*pEdge(3)*cos((double)vv*InvNBin*DUE_PI);
323  double PosY = pEdge(3) + InvNBin*(double)v*pEdge(3)*sin((double)vv*InvNBin*DUE_PI);
324  fprintf(FileToWrite,"{x(%lf %lf %lf)}\n",
325  PosX,PosY,Plot[v][vv]/Max*pEdge(3));
326  }
327  else{
328  double PosX = pEdge(CLat1)*.5 + (v - vRef)*InvNBin*pEdge(CLat1);
329  double PosY = pEdge(CLat2)*.5 + (vv - vRef)*InvNBin*pEdge(CLat2);
330  fprintf(FileToWrite,"{x(%lf %lf %lf)}\n",
331  PosX,PosY,Plot[v][vv]/Max*pEdge(3));
332  }
333  //printf("%d %d %lf %lf %lf\n",v,vv,PosX,PosY,Plot[v][vv]);
334  }
335  }
336  }
337  fclose(FileToWrite);
338  }
339  free(Plot);
340  free(dDensity);
341  free(dDensity1);
342  return 0;
343 }
344 int ElPoly::ScatteringF(int NBin,int How){
345  double **Plot;
346  Plot = (double **)calloc(NBin,sizeof(double));
347  for(int i=0;i<NBin;i++){
348  *(Plot+i) = (double *)calloc(NBin,sizeof(double));
349  }
350  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
351  for(int f=NFile[0];f<NFile[1];f++){
352  Processing(f);
353  //printf("Opening: %s\n",cFile[f]);
354  if(OpenRisk(cFile[f],BF_CHAIN))return 1;
355  Scattering2d(Plot,NBin,NChType);
356  }
357  printf("\n");
358  double Max=0.;
359  double InvNBin=1./(double)NBin;
360  for(int v=0;v<NBin;v++){
361  for(int vv=0;vv<NBin;vv++){
362  Plot[v][vv] = POS(Plot[v][vv]);
363  if(Max < Plot[v][vv])
364  Max = Plot[v][vv];
365  }
366  }
367  Max = 1./Max;
368  FILE *FileToWrite = NULL;
369  if(How == 0)
370  FileToWrite = fopen("Scattering.xvl","w");
371  else
372  FileToWrite = fopen("Scattering2.xvl","w");
373  fprintf(FileToWrite,"#l(%lf %lf %lf) v(%d) d[%s]\n",2.*pEdge(3),2.*pEdge(3),pEdge(3),NBin,ChooseDraw(EL_QUAD));
374  int vRef = NBin/2;
375  for(int q=0;q<4;q++){
376  double Signx = !(q%2) ? 1. : -1.;
377  double Signy = q > 1 ? 1. : -1.;
378  for(int vx=0;vx<NBin;vx++){
379  for(int vy=0;vy<NBin;vy++){
380  //if(Plot[v][vv] > 0.)
381  {
382  double PosX = Signx*vx*InvNBin*pEdge(CLat1);
383  double PosY = Signy*vy*InvNBin*pEdge(CLat2);
384  //printf("%d %d %lf %lf %lf\n",v,vv,PosX,PosY,Plot[v][vv]);
385  fprintf(FileToWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",
386  PosX,PosY,Plot[vx][vy]*.001,Plot[vx][vy],0.,0.);
387  }
388  }
389  }
390  }
391  fclose(FileToWrite);
392  free(Plot);
393  return 0;
394 }
395 //obsolete?
397  int NTimes = 20;
398  double *ChainPArea = (double *)calloc(3*NTimes,sizeof(double));
399  double MeanSigma = 0.;
400  int SubDiv[3] = {NTimes,NTimes,1};
401  for(int f=NFile[0];f<NFile[1];f++){
402  Processing(f);
403  if(OpenRisk(cFile[f],BF_PART))return 1;
404  // OrderPos();
405  SubDiv[CLat1] = 1;
406  SubDiv[CLat2] = 1;
407  for(int t=0;t<NTimes;t++){
408  if( (t%2)==0)
409  SubDiv[CLat1] += 1;
410  if( (t%2-1)==0)
411  SubDiv[CLat2] += 1;
412  double NEdge = (double)(SubDiv[CLat1]*SubDiv[CLat2]);
413  double *Plot = (double *)calloc((int)NEdge,sizeof(double));
414  double Area = pEdge(CLat1)*pEdge(CLat2)/NEdge;
415  double SumX = 0.;
416  double SumX2 = 0.;
417  for(int p=0;p<pNPart();p++){
418  int vx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*SubDiv[CLat1]);
419  if(vx > SubDiv[CLat1]) continue;
420  int vy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*SubDiv[CLat2]);
421  if(vy > SubDiv[CLat2]) continue;
422  Plot[vx*SubDiv[CLat2] + vy] += 1.;
423  }
424  for(int n=0;n<NEdge;n++){
425  SumX += Plot[n]/Area;
426  SumX2 += QUAD((Plot[n]/Area));
427  Plot[n] /= Area;
428  }
429  int Valori = 20;
430  double *Intervalli = (double *)calloc(Valori,sizeof(double));
431  int IfNorm = 0;
432  MOMENTI m1 = Mat->Distribuzione(Plot,(int)NEdge,Intervalli,Valori,IfNorm);
433  ChainPArea[3*t+0] = NEdge;
434  MeanSigma += m1.Uno;
435  double Media = SumX/NEdge;
436  ChainPArea[3*t+1] = sqrt((SumX2 - Media*Media*NEdge)/(NEdge-1));
437  printf("%lf %lf-%lf\n",m1.Uno,ChainPArea[3*t+1],m1.Due);
438  ChainPArea[3*t+2] = Area*QUAD(m1.Due)/QUAD(m1.Uno);
439  //QUAD(Area)*ChainPArea[4*t+2]/QUAD(ChainPArea[4*t+1]);
440  free(Plot);
441  free(Intervalli);
442  }
443  }
444  char *FileName = (char *)calloc(60,sizeof(char));
445  sprintf(FileName,"ChainPArea%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho());
446  FILE *FileToWrite = fopen(FileName,"w");
447  char cSystem[STRSIZE];
448  SysDef(cSystem);
449  fprintf(FileToWrite,"# %s",cSystem);
450  fprintf(FileToWrite,"#%lf\n",MeanSigma/(double)(NFile[1]-NFile[0]));
451  fprintf(FileToWrite,"# NEdge SDeviation Area*Dev/NEdge\n");
452  for(int t=0;t<NTimes;t++)
453  fprintf(FileToWrite,"%lf %lf %lf\n",ChainPArea[3*t+0],ChainPArea[3*t+1],ChainPArea[3*t+2]);
454  fclose(FileToWrite);
455  free(FileName);
456  free(ChainPArea);
457  return 0;
458 }
459 int ElPoly::SpectrumF(int NSample){
460  int NHalf = NSample/2;
461  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
462  double InvNSample = 1./(double)NSample;
463  double InvNSample2 = 1./(double)SQR(NSample);
464  double *Plot = (double *) calloc(SQR(NSample),sizeof(double));
465  double *PlotUp = (double *) calloc(SQR(NSample),sizeof(double));
466  double *PlotDown = (double *) calloc(SQR(NSample),sizeof(double));
467  double *CountUp = (double *) calloc(SQR(NSample),sizeof(double));
468  double *CountDown = (double *) calloc(SQR(NSample),sizeof(double));
469  double *PlotBS = (double *) calloc(SQR(NSample),sizeof(double));
470  double *PlotA = (double *)calloc(SQR(NSample), sizeof(double));
471  double *PlotAUp = (double *)calloc(SQR(NSample), sizeof(double));
472  double *PlotADown = (double *)calloc(SQR(NSample), sizeof(double));
473  double *PlotThin = (double *)calloc(SQR(NSample), sizeof(double));
474  double *Count = (double *) calloc(NSample*NSample,sizeof(double));
475 #ifdef USE_FFTW
476  fftw_complex *out = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
477  fftw_complex *in = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
478  fftw_plan plan = fftw_plan_dft_2d(NSample,NSample,
479  in, out,FFTW_FORWARD,FFTW_PATIENT);
480  fftw_complex *outUp = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
481  fftw_complex *inUp = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
482  fftw_plan planUp = fftw_plan_dft_2d(NSample,NSample,
483  inUp, outUp,FFTW_FORWARD,FFTW_PATIENT);
484  fftw_complex *outDown = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
485  fftw_complex *inDown = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
486  fftw_plan planDown = fftw_plan_dft_2d(NSample,NSample,
487  inDown, outDown,FFTW_FORWARD,FFTW_PATIENT);
488  //FFTW_MEASURE);
489  fftw_complex *outThin = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
490  fftw_complex *inThin = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
491  fftw_plan planThin = fftw_plan_dft_2d(NSample,NSample,
492  inThin, outThin,FFTW_FORWARD,FFTW_PATIENT);
493  //FFTW_MEASURE);
494 #endif
495  for(int f=NFile[0];f<NFile[1];f++){
496  Processing(f);
497  // memset(Plot,0.,NSample*NSample*sizeof(double));
498  // memset(Count,0.,NSample*NSample*sizeof(double));
499  for(int s=0;s<SQR(NSample);s++){
500  Plot[s] = 0.;
501  Count[s] = 0.;
502  PlotUp[s] = 0.;
503  CountUp[s] = 0.;
504  PlotDown[s] = 0.;
505  CountDown[s] = 0.;
506  }
507  if(OpenRisk(cFile[f],BF_CHAIN))return 0;
508  for(int p=0;p<pNPart();p++){
509  int sx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*NSample);
510  int sy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*NSample);
511  if(sx < 0 || sx >= NSample) continue;
512  if(sy < 0 || sy >= NSample) continue;
513  Plot[sx*NSample+sy] += pPos(p,CNorm);// - pCm(CNorm);
514  Count[sx*NSample+sy] += 1.;
515  //if(Pm[p].Typ != 1) continue;
516  if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){
517  PlotUp[sx*NSample+sy] += pPos(p,CNorm);
518  CountUp[sx*NSample+sy] += 1.;
519  }
520  else {
521  PlotDown[sx*NSample+sy] += pPos(p,CNorm);
522  CountDown[sx*NSample+sy] += 1.;
523  }
524  }
525  double MeanThick = 0.;
526  for(int s=0;s<SQR(NSample);s++){
527  Plot[s] /= Count[s] > 0. ? Count[s] : 1.;
528  PlotUp[s] /= CountUp[s] > 0. ? CountUp[s] : 1.;
529  PlotDown[s] /= CountDown[s] > 0. ? CountDown[s] : 1.;
530  MeanThick += CountUp[s] - CountDown[s];
531  }
532  MeanThick /= (double)SQR(NSample);
533  //InterBSpline2D(Plot,PlotBS,NSample,NSample);
534  //InterBSpline2D(PlotBS,Plot,NSample,NSample);
535 #ifdef USE_FFTW
536  if(1==1){
537  for(int s=0;s<SQR(NSample);s++){
538  in[s][0] = Plot[s];
539  in[s][1] = 0.;
540  inThin[s][0] = PlotUp[s]-PlotDown[s];
541  inThin[s][1] = 0.;
542  inUp[s][0] = PlotUp[s];
543  inUp[s][1] = 0.;
544  inDown[s][0] = PlotDown[s];
545  inDown[s][1] = 0.;
546  }
547  fftw_execute(plan);
548  fftw_execute(planThin);
549  fftw_execute(planUp);
550  fftw_execute(planDown);
551  for(int sx=0;sx<NSample;sx++){
552  int ssx = sx < NHalf ? sx + NHalf : sx - NHalf;
553  for(int sy=0;sy<NSample;sy++){
554  int ssy = sy < NHalf ? sy + NHalf : sy - NHalf;
555  PlotA[ssx*NSample+ssy] += (SQR(out[sx*NSample+sy][0])+SQR(out[sx*NSample+sy][1]))*SQR(InvNSample2);
556  PlotThin[ssx*NSample+ssy] += (SQR(outThin[sx*NSample+sy][0])+SQR(outThin[sx*NSample+sy][1]))*SQR(InvNSample2);
557  PlotAUp[ssx*NSample+ssy] += (SQR(outUp[sx*NSample+sy][0])+SQR(outUp[sx*NSample+sy][1]))*SQR(InvNSample2);
558  PlotADown[ssx*NSample+ssy] += (SQR(outDown[sx*NSample+sy][0])+SQR(outDown[sx*NSample+sy][1]))*SQR(InvNSample2);
559  }
560  }
561  }
562  else {
563  int NMax = NSample;
564  double *st = Plot;
565  double *sw = PlotA;
566  double dNMax = 1./(double)NMax;
567  int NHalf = (int)(NMax/2.);
568  for(int kx=-NHalf;kx<NHalf;kx++){
569  double qx = kx*dNMax;
570  for(int ky=-NHalf;ky<NHalf;ky++){
571  double qy = ky*dNMax;
572  double Re2=0.,Im2=0.;
573  double Re1=0.,Im1=0.;
574  for(int lx=0;lx<NMax;lx++){
575  for(int ly=0;ly<NMax;ly++){
576  double Arg = 2.*M_PI*(lx*qx + ly*qy);
577  double cy = cos(Arg);
578  double sy = sin(Arg);
579  Re1 += st[lx*NMax+ly]*cy;
580  Im1 += st[lx*NMax+ly]*sy;
581  }
582  }
583  int kkx = kx + NHalf;
584  int kky = ky + NHalf;
585  sw[kkx*NMax+kky] += SQR(Re1*dNMax) + SQR(Im1*dNMax);
586  }
587  }
588  }
589  if(1==0){
590  FILE *FMidplane = fopen("Midplane.dat","w");
591  FILE *FThickness = fopen("Thickness.dat","w");
592  for(int sx=0;sx<NSample;sx++){
593  double x = sx*pEdge(CLat1)*InvNSample;
594  for(int sy=0;sy<NSample;sy++){
595  double y = sy*pEdge(CLat2)*InvNSample;
596  fprintf(FMidplane,"%lf %lf %lf \n",x,y, Plot[sx*NSample+sy]);
597  fprintf(FThickness,"%lf %lf %lf \n",x,y,PlotUp[sx*NSample+sy]-PlotDown[sx*NSample+sy]);
598  }
599  }
600  fclose(FMidplane);
601  fclose(FThickness);
602  }
603 #else
604  Mat->Spettro2d(Plot,PlotA,NSample);
605 #endif //FFTW3_H
606  }
607  printf("\n");
608 #ifdef OMPI_MPI_H
609  MPI_Allreduce(MPI_IN_PLACE,PlotA,NSample*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
610  int Rank=0;
611  MPI_Comm_rank(Proc->CommGrid, &Rank);
612  if(Rank==0){
613 #endif
614  char *NomeFile = (char*)calloc(60,sizeof(char));
615  sprintf(NomeFile,"Spectrum.xvl",pchiN(),pkappaN(),prho());
616  FILE *FileToWrite = fopen(NomeFile,"w");
617  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[squaremesh]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NSample);
618  for(int sx=0;sx<NSample;sx++){
619  for(int sy=0;sy<NSample;sy++){
620  double x = pEdge(CLat1)*sx*InvNSample;
621  double y = pEdge(CLat2)*sy*InvNSample;
622  fprintf(FileToWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",x,y,PlotA[sx*NSample+sy]*InvNFile,fabs(log10(PlotA[sx*NSample+sy]*InvNFile)),0.,0.);
623  }
624  }
625  printf("Write\n");
626  fclose(FileToWrite);
627  sprintf(NomeFile,"Spectrum.dat",pchiN(),pkappaN(),prho());
628  FileToWrite = fopen(NomeFile,"w");
629  fprintf(FileToWrite,"# 1) q^2 2) h(q)^2*L_xL_y | q_x=s2pi/(L_x) NSample %d\n",NSample);
630  for(int kx=0;kx<NSample;kx++){
631  int kx1 = kx - NSample/2;
632  double qx = 2.*M_PI*kx*pInvEdge(CLat1);
633  double qxt = 2.*M_PI*kx1*pInvEdge(CLat1);
634  for(int ky=0;ky<NSample;ky++){
635  int ky1 = ky - NSample/2;
636  double qy = 2.*M_PI*ky*pInvEdge(CLat2);
637  double qyt = 2.*M_PI*ky1*pInvEdge(CLat2);
638  double qq = SQR(qxt) + SQR(qyt);
639  double Spe = PlotA[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2);
640  double SpeUp = PlotAUp[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2);
641  double SpeDown = PlotADown[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2);
642  double SpeThin = PlotThin[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2);
643  fprintf(FileToWrite,"%lf %g %g %g %g\n",qq,Spe,SpeUp,SpeDown,SpeThin);
644  }
645  }
646  fclose(FileToWrite);
647  free(NomeFile);
648 #ifdef OMPI_MPI_H
649  }
650 #endif
651  free(Plot);
652  free(PlotA);
653  free(PlotBS);
654  free(Count);
655 #ifdef USE_FFTW
656  fftw_destroy_plan(plan);
657 #endif
658  return 0;
659 }
660 void ElPoly::Midplane(int NSample){
661  int NHalf = NSample/2;
662  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
663  double InvNSample = 1./(double)NSample;
664  double *PlotBS = (double *) calloc(SQR(NSample),sizeof(double));
665  double *Count = (double *) calloc(NSample*NSample,sizeof(double));
666  double *Plot = (double *) calloc(SQR(NSample),sizeof(double));
667  double *PlotA = (double *) calloc(SQR(NSample),sizeof(double));
668  for(int f=NFile[0];f<NFile[1];f++){
669  Processing(f);
670  // memset(Plot,0.,NSample*NSample*sizeof(double));
671  // memset(Count,0.,NSample*NSample*sizeof(double));
672  for(int sx=0;sx<NSample;sx++){
673  for(int sy=0;sy<NSample;sy++){
674  Plot[sx*NSample+sy] = 0.;
675  }
676  }
677  if(OpenRisk(cFile[f],NBackFold))return;
678  for(int p=0;p<pNPart();p++){
679  int sx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*NSample);
680  int sy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*NSample);
681  if(sx < 0 || sx >= NSample) continue;
682  if(sy < 0 || sy >= NSample) continue;
683  Plot[sx*NSample+sy] += pPos(p,CNorm);// - pCm(CNorm);
684  Count[sx*NSample+sy] += 1.;
685  }
686  for(int sx=0;sx<NSample;sx++){
687  for(int sy=0;sy<NSample;sy++){
688  Plot[sx*NSample+sy] /= Count[sx*NSample+sy] > 0. ? Count[sx*NSample+sy] : 1.;
689  PlotA[sx*NSample+sy] += Plot[sx*NSample+sy];
690  }
691  }
692  // InterBSpline2D(Plot,PlotBS,NSample,NSample);
693  // InterBSpline2D(PlotBS,Plot,NSample,NSample);
694  char FName[60];
695  sprintf(FName,"Midplane%05d.dat",f);
696  FILE *FMid = fopen(FName,"w");
697  for(int sy=0;sy<NSample;sy++){
698  for(int sx=0;sx<NSample;sx++){
699  double x = pEdge(CLat1)*sx*InvNSample;
700  double y = pEdge(CLat2)*sy*InvNSample;
701  fprintf(FMid,"%lf %lf %lf\n",x,y,Plot[sx*NSample+sy]);
702  }
703  }
704  fclose(FMid);
705  }
706  printf("\n");
707  FILE *FMid = fopen("average_middplane.dat","w");
708  for(int sy=0;sy<NSample;sy++){
709  for(int sx=0;sx<NSample;sx++){
710  double x = pEdge(CLat1)*sx*InvNSample;
711  double y = pEdge(CLat2)*sy*InvNSample;
712  fprintf(FMid,"%lf %lf %lf\n",x,y,PlotA[sx*NSample+sy]*InvNFile);
713  }
714  }
715  fclose(FMid);
716  free(Plot);
717  free(Count);
718  free(PlotBS);
719 }
720 void ElPoly::SpectrumMidplane(int NSample){
721  // FILE *Ciccia = fopen("Ciccia.dat","w");
722  // for(int p=0;p<pNPart();p++){
723  // fprintf(Ciccia,"%lf %lf\n",SQR(pPos(p,0)-.5*pEdge(0))+SQR(pPos(p,1)-.5*pEdge(1)),pPos(p,2));
724  // }
725  // fclose(Ciccia);
726  // return;
727  int NHalf = NSample/2;
728  double InvNSample = 1./(double)NSample;
729  double *Plot = (double *) calloc(NSample*NSample,sizeof(double));
730  double *Count = (double *) calloc(NSample*NSample,sizeof(double));
731  double *PlotA = (double *) calloc(NSample*NSample,sizeof(double));
732  int NPoint = (int)(NSample*sqrt(2.));
733  double *Density = (double *) calloc(NPoint,sizeof(double));
734 #ifdef USE_FFTW
735  fftw_complex *out = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
736  fftw_complex *in = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex));
737  fftw_plan plan1 = fftw_plan_dft_2d(NSample,NSample,
738  in, out,FFTW_FORWARD,FFTW_PATIENT);
739  for(int f=NFile[0];f<NFile[1];f++){
740  Processing(f);
741  for(int sx=0;sx<NSample;sx++){
742  for(int sy=0;sy<NSample;sy++){
743  Plot[sx*NSample+sy] = 0.;
744  Count[sx*NSample+sy] = 0.;
745  }
746  }
747  FILE *Ciccia = fopen(cFile[f],"r");
748  char cLine[256];
749  for(int p=0;!(fgets(cLine,STRSIZE,Ciccia)==NULL);p++){
750  double x = 0., y = 0., z = 0.;
751  sscanf(cLine,"%lf %lf %lf\n",&x,&y,&z);
752  int v = (int)(x*NSample*pInvEdge(CLat1));//(pEdge(0)+.4));
753  int vv = (int)(y*NSample*pInvEdge(CLat2));//(pEdge(1)+.4));
754  if(v < 0 || v >= NSample) continue;
755  if(vv < 0 || vv >= NSample) continue;
756  Plot[v*NSample+vv] += z;
757  Count[v*NSample+vv] += 1.;
758  }
759  for(int s=0;s<NSample;s++){
760  for(int ss=0;ss<NSample;ss++){
761  Plot[s*NSample+ss] /= Count[s*NSample+ss] > 0. ? Count[s*NSample+ss] : 1.;
762  }
763  }
764  fclose(Ciccia);
765  if(1==0){
766  for(int sx=0;sx<NSample;sx++){
767  for(int sy=0;sy<NSample;sy++){
768  in[sx*NSample+sy][0] = Plot[sx*NSample+sy];
769  in[sx*NSample+sy][1] = 0.;
770  }
771  }
772  fftw_execute(plan1);
773  for(int sx=0;sx<NSample;sx++){
774  int sx1 = sx + NHalf;
775  if(sx1 >= NSample) sx1 -= NSample;
776  for(int sy=0;sy<NSample;sy++){
777  int sy1 = sy + NHalf;
778  if(sy1 >= NSample) sy1 -= NSample;
779  PlotA[sx*NSample+sy] += (SQR(out[sx*NSample+sy][0])+SQR(out[sx*NSample+sy][1]))*SQR(InvNSample);
780  //PlotA[sx*NSample+sy] += Plot[sx][sy];
781  }
782  }
783  }
784  else {
785  int NMax = NSample;
786  double *st = Plot;
787  double *sw = PlotA;
788  double dNMax = 1./(double)NMax;
789  int NHalf = (int)(NMax/2.);
790  for(int kx=-NHalf;kx<NHalf;kx++){
791  double qx = kx*dNMax;
792  for(int ky=-NHalf;ky<NHalf;ky++){
793  double qy = ky*dNMax;
794  double Re2=0.,Im2=0.;
795  double Re1=0.,Im1=0.;
796  for(int lx=0;lx<NMax;lx++){
797  for(int ly=0;ly<NMax;ly++){
798  double Arg = 2.*M_PI*(lx*qx + ly*qy);
799  double cy = cos(Arg);
800  double sy = sin(Arg);
801  Re1 += st[lx*NMax+ly]*cy;
802  Im1 += st[lx*NMax+ly]*sy;
803  }
804  }
805  int kkx = kx + NMax/2;
806  int kky = ky + NMax/2;
807  sw[kkx*NMax+kky] += SQR(Re1*dNMax) + SQR(Im1*dNMax);
808  }
809  }
810  }
811  }
812 #endif // USE_FFTW
813  printf("\n");
814 #ifdef OMPI_MPI_H
815  MPI_Allreduce(MPI_IN_PLACE,PlotA,NSample*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
816  int Rank=0;
817  MPI_Comm_rank(Proc->CommGrid, &Rank);
818  if(Rank==0){
819 #endif
820  char *NomeFile = (char*)calloc(60,sizeof(char));
821  sprintf(NomeFile,"SpectrumChi%.0fKappa%.0fRho%.0f.xvl",pchiN(),pkappaN(),prho());
822  FILE *FileToWrite = fopen(NomeFile,"w");
823  // fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[squaremesh]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NSample);
824  double DivInv = 1./(double)(NFile[1]-NFile[0]);
825  for(int s=0;s<NSample;s++){
826  for(int ss=0;ss<NSample;ss++){
827  fprintf(FileToWrite,"%lf %lf %lf\n",pEdge(CLat1)*s*InvNSample,pEdge(CLat2)*ss*InvNSample,PlotA[s*NSample+ss]*DivInv);
828  // fprintf(FileToWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",pEdge(CLat1)*s*InvNSample,pEdge(CLat2)*ss*InvNSample,PlotA[s*NSample+ss]*DivInv,0.,0.,log10(PlotA[s*NSample+ss]*DivInv));
829  }
830  }
831  fclose(FileToWrite);
832  sprintf(NomeFile,"SpectrumChi%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho());
833  FileToWrite = fopen(NomeFile,"w");
834  double qHalf = QUAD(pEdge(CLat1)*(NHalf)*InvNSample)+QUAD(pEdge(CLat2)*(NHalf)*InvNSample);
835  double dx = pEdge(CLat1)*InvNSample;
836  double dy = pEdge(CLat2)*InvNSample;
837  for(int kx=0;kx<NSample;kx++){
838  int kx1 = kx - NSample/2;
839  double qx = kx*InvNSample;
840  double qxt = kx1*InvNSample;
841  for(int ky=0;ky<NSample;ky++){
842  int ky1 = ky - NSample/2;
843  double qy = ky*InvNSample;
844  double qyt = ky1*InvNSample;
845  double qq = SQR(qx/dx - .5/dx) + SQR(qy/dy - .5/dy);
846  //double qq = SQR(qx) + SQR(qy);
847  fprintf(FileToWrite,"%lf %g\n",qq,PlotA[kx*NSample+ky]*DivInv);
848  //fprintf(FileToWrite,"%lf %g\n",QUAD(pEdge(CLat1)*(s-NHalf)*InvSample)+QUAD(pEdge(CLat2)*(ss-NHalf)*InvSample),PlotA[s*NSample+ss]*DivInv);
849  }
850  }
851  // for(int p=1;p<NPoint;p++)
852  // fprintf(FileToWrite,"%d %lf\n",p,Density[p]*DivInv);
853  fclose(FileToWrite);
854  free(NomeFile);
855 #ifdef OMPI_MPI_H
856  }
857 #endif
858  free(Density);
859  free(PlotA);
860  free(Plot);
861 }
862 void ElPoly::HeaderAverage(int nNano){
863  if(nNano >= pNNano()){
864  printf("The specified nanoparticle doesn't exist\n");
865  nNano = 0;
866  // return;
867  }
868  int NDim = 2;
869  int NumFile = NFile[1]-NFile[0];
870  double *Area = (double *)calloc(NumFile,sizeof(double));
871  double *NanoPos = (double *)calloc(6*NumFile,sizeof(double));
872  double *NanoDiff = (double *)calloc(2*NumFile,sizeof(double));
873  double *NanoCount = (double *)calloc(NumFile,sizeof(double));
874  double *NanoDist = (double *)calloc(2*NumFile,sizeof(double));
875  double *Time = (double *)calloc(NumFile,sizeof(double));
876  double NanoInit[3] = {Nano->Pos[0],Nano->Pos[1],Nano->Pos[2]};
877  double AreaAv = 0.;
878  double AreaErr = 0.;
879  double InitTime = pTime();
880  //indentation
881  double InitPos = pPos(0,CNorm) - pPos(1,CNorm) - MAX(Nano[0].Rad,Nano[0].Height);
882  for(int f=NFile[0];f<NFile[1];f++){
883  Processing(f);
884  FILE *File2Read;
885  if((File2Read = fopen(cFile[f],"r"))==0){
886  printf("The file is missing\n");
887  return ;
888  }
889  ReadHeader(File2Read);
890  fclose(File2Read);
891  AreaAv += pEdge(CLat1)*pEdge(CLat2);
892  AreaErr += SQR(pEdge(CLat1)*pEdge(CLat2));
893  Area[f] = pEdge(CLat1)*pEdge(CLat2);
894  Time[f] = pTime() - InitTime;
895  for(int d=0;d<3;d++) NanoPos[f*6+d] = Nano[nNano].Pos[d];
896  NanoPos[f*6+3] = Nano[nNano].Rad;
897  NanoPos[f*6+4] = Nano[nNano].Height;
898  NanoPos[f*6+5] = Nano[nNano].Area;
899  for(int ff=f;ff>=0;ff--){
900  int fff = f-ff;
901  NanoDiff[fff*2] += SQR(NanoPos[f*6+0]-NanoPos[ff*6+0]) + SQR(NanoPos[f*6+1]-NanoPos[ff*6+1]);
902  NanoCount[fff] += 1.;
903  }
904  NanoDiff[f*2+1] = SQR(NanoPos[f*6+0]-NanoInit[0]) + SQR(NanoPos[f*6+1]-NanoInit[1]);
905  if(pNNano() > 1){
906  //indentation
907  NanoDist[f*2] = InitPos - pNanoPos(0,CNorm) + pNanoPos(1,CNorm) + MAX(Nano[nNano].Rad,Nano[nNano].Height);
908  double Dist = 0.;
909  for(int d=0;d<3;d++){
910  double Pos1 = pNanoPos(0,d) - floor(pNanoPos(0,d)*pInvEdge(d))*pEdge(d);
911  Dist += SQR(pNanoPos(0,d)-pNanoPos(1,d));
912  }
913  Dist = sqrt(Dist);
914  //two np dist
915  NanoDist[f*2+1] = Dist;
916  }
917  }
918  AreaAv /= (double)(NFile[1] - NFile[0]);
919  AreaErr = sqrt( AreaErr - (NFile[1] - NFile[0])*SQR(AreaAv))/(double)(NFile[1] - NFile[0]);
920  printf("\n");
921  printf("Area %lf \\pm %lf\n",AreaAv,AreaErr);
922  FILE *File2Write = fopen("AreaTime.dat","w");
923  fprintf(File2Write,"# Time Area %lf pm %lf\n",AreaAv,AreaErr);
924  char FName[120];
925  sprintf(FName,"NanoDiffR%.1fS%.1fH%.1f.dat",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height);
926  FILE *FNano = fopen(FName,"w");
927  fprintf(FNano,"# t <msd> msd Dz NpDist rad hei\n");
928  for(int f=0;f<NumFile;f++){
929  NanoDiff[f*2] /= NanoCount[f] > 0. ? NanoCount[f] : 1.;
930  fprintf(File2Write,"%lf %lf\n",Time[f],Area[f]);
931  fprintf(FNano,"%lf %lf %lf %lf %lf %lf %lf\n",Time[f],NanoDiff[f*2 ],NanoDiff[f*2+1],NanoDist[f*2],NanoDist[f*2+1],NanoPos[f*6+3],NanoPos[f*6+4]);
932  }
933  fclose(File2Write);
934  fclose(FNano);
935  free(NanoDiff);
936  free(NanoDist);
937  free(Area);
938  free(Time);
939 }
940 //only one block
941 void ElPoly::WidomOut(char *NrgFile,int NBin){
942  int NChains = 1;
943  Block[0].NChain -= NChains;
944  WriteXvt("SnapOut.dat");
945  FILE *FileLast = fopen("LastChain.dat","w");
946  for(int p=pNPart()-(NChains)*pNPCh();p<pNPart();p++){
947  fprintf(FileLast,"%lf %lf %lf %lf %lf %lf %d\n",pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0),pVel(p,1),pVel(p,2),pType(p));
948  }
949  fclose(FileLast);
950 }
951 //only one block
953  int NChains = 1;
954  Block[0].NChain -= NChains;
955  char FileName[60];
956  sprintf(FileName,"SnapOut%05d.dat",pNChain()-1);
957  WriteXvt(FileName);
958  for(int c=0;c<pNChain()-1;c++){
959  SwapChain(c,pNChain()-1);
960  sprintf(FileName,"SnapOut%05d.dat",c);
961  WriteXvt(FileName);
962  }
963 }
964 //the system has an additional fake chain
966  int NFile = 1000;
967  char FileName[60];
968  SetNPart(pNPart()+pNPCh());
969  SetNChain(pNChain()+1);
970  Block[0].NChain += 1;
971  int p1 = pNPCh()*(pNChain()-1);
972  double Var = sqrt(SQR(pReOverCutOff())/(double)(pNPCh()-1)/3.)/2.;
973  for(int f=0;f<NFile;f++){
974  for(int d=0;d<3;d++){
975  SetPos(p1,d,Mat->Casuale()*pEdge(d));
976  }
977  for(int p=p1+1;p<pNPart();p++){
978  for(int d=0;d<3;d++){
979  SetPos(p1,d,pPos(p-1,d) + Mat->Gaussiano(0.,Var));
980  }
981  if( p%pNPCh() >= Block[0].Asym )
982  SetType(p,1);
983  }
984  sprintf(FileName,"SnapIn%05d.dat",f);
985  Write(FileName);
986  }
987 }
988 //the system has an additional fake chain
989 void ElPoly::WidomIn(char *NrgFile,int NBin){
990  int NChains = 1;
991  Block[0].NChain += NChains;
992  WriteXvt("SnapIn.dat");
993  int NAddPart = NChains*pNPCh();
994  FILE *FileLast = fopen("SnapIn.dat","a");
995  PART *Pn = (PART *)calloc(NAddPart,sizeof(PART));
996  double GaussVar = sqrt(SQR(pReOverCutOff())/(double)(pNPCh()-1)/3.)*2./3.;
997  for(int c=0;c<NChains;c++){
998  int p1 = c*pNPCh();
999  for(int d=0;d<3;d++){
1000  SetPos(p1,d,Mat->Casuale()*pEdge(d));
1001  }
1002  for(int p=1;p<pNPCh();p++){
1003  for(int d=0;d<3;d++){
1004  Pn[p+p1].Pos[d] = Pn[p-1+p1].Pos[d] + Mat->Gaussiano(0,GaussVar);
1005  }
1006  Pn[p+p1].Typ = 0;
1007  if(p >= Block[0].Asym){
1008  Pn[p+p1].Typ = 1;
1009  }
1010  }
1011  }
1012  for(int p=0;p<NAddPart;p++){
1013  fprintf(FileLast,"%lf %lf %lf %lf %lf %lf %d\n",Pn[p].Pos[0],Pn[p].Pos[1],Pn[p].Pos[2],Pn[p].Vel[0],Pn[p].Vel[1],Pn[p].Vel[2],Pn[p].Typ);
1014  }
1015  free(Pn);
1016  fclose(FileLast);
1017 }
1018 void ElPoly::End2EndDistr(char *OutFile){
1019  char FileName[256];
1020  for(int f=NFile[0];f<NFile[1];f++){
1021  sprintf(FileName,"%s%04d.dat",OutFile,f);
1022  FILE *File2Write = fopen(FileName,"w");
1023  Processing(f);
1024  if(OpenRisk(cFile[f],BF_NO))return;
1025  for(int c=0;c<pNChain();c++){
1026  int p1 = c*pNPCh();
1027  int p2 = c*pNPCh() + pNPCh() - 1;
1028  double End2End = 0.;
1029  for(int d=0;d<3;d++){
1030  End2End += SQR(pPos(p1,d) - pPos(p2,d));
1031  }
1032  fprintf(File2Write,"%lf ",End2End);
1033  fprintf(File2Write,"\n");
1034  }
1035  fclose(File2Write);
1036  }
1037 }
1038 void ElPoly::Decoupling(int What){
1039  NFileTot = NFile[1] - NFile[0];
1040  double *Angles = (double *)calloc((NFile[1]-NFile[0])*pNChain(),sizeof(double));
1041  double *Count = (double *)calloc((NFile[1]-NFile[0]),sizeof(double));
1042  double *AnglesNano = (double *)calloc((NFile[1]-NFile[0])*pNNano(),sizeof(double));
1043  double *Time = (double *)calloc((NFile[1]-NFile[0]),sizeof(double));
1044  double *AngleDiff = (double *)calloc(pNChain(),sizeof(double));
1045  double *AngleDiff2 = (double *)calloc(pNChain(),sizeof(double));
1046  Vettore Ax0(1.,0.,0.);
1047  Vettore Ax2(0.,0.,1.);
1048  double InvPhob = Block[0].Asym > 0 ? 1./(double)Block[0].Asym : 1.;
1049  double InvPhil = Block[0].Asym > 0 ? 1./(double)(pNPCh() - Block[0].Asym) : 1.;
1050  double DistRel[4];
1051  Vettore ChDir(0.,0.,0.);
1052  Vettore NanoAx(0.,0.,0.);
1053  double Average = 0.;
1054  double Variance = 0.;
1055  double InitTime = pTime();
1056  for(int f=NFile[0];f<NFile[1];f++){
1057  Processing(f);
1058  int ff = f - NFile[0];
1059  if(OpenRisk(cFile[f],BF_CHAIN))return ;
1060  Time[f] = pTime() - InitTime;
1061  for(int n=0;n<pNNano();n++){
1062  for(int d=0;d<3;d++){
1063  NanoAx.Set(Nano[n].Axis[d],d);
1064  }
1065  NanoAx.Set(0.,CNorm);
1066  AnglesNano[ff*pNNano()+n] = Ax0.Angle(&NanoAx);
1067  }
1068  for(int c=0,b=0;c<pNChain();c++){
1069  if(What == 0){
1070  ChDir.Set(Ch[c].Dir[0],0);
1071  ChDir.Set(Ch[c].Dir[1],1);
1072  double Angle = Ax0.Angle(&ChDir);
1073  if(isnan(Angle)) Angles[ff*pNChain()+c] = Angles[(ff-1)*pNChain()+c];
1074  else Angles[ff*pNChain()+c] = Ax0.Angle(&ChDir);
1075  }
1076  else if(What == 1){
1077  Angles[ff*pNChain()+c] = TwoPartDist(c*pNPCh(),(c+1)*pNPCh()-1,DistRel);
1078  }
1079  Average += Angles[ff*pNChain()+c];
1080  Variance += SQR(Angles[ff*pNChain()+c]);
1081  }
1082  }
1083  printf("\n");
1084  double NTot = (double)(pNChain()*(NFile[1]-NFile[0]));
1085  Average /= NTot > 0. ? NTot : 1.;
1086  double PreVar = Variance;
1087  Variance = (Variance - SQR(Average)*NTot)/(NTot-1);
1088  if(isnan(Variance)){
1089  printf("Variance is not a number %lf %lf %lf\n",PreVar,Average,NTot);
1090  Variance = 1.;
1091  }
1092  for(int f=0;f<NFileTot;f++){
1093  for(int c=0;c<pNChain();c++){
1094  AngleDiff2[f] += (Angles[f*pNChain()+c]-Average)*(Angles[0*pNChain()+c]-Average);
1095  }
1096  for(int ff = f+1; ff<NFileTot;ff++){
1097  int fff = ff - f;
1098  for(int c=0;c<pNChain();c++){
1099  AngleDiff[fff] += (Angles[ff*pNChain()+c]-Average)*(Angles[f*pNChain()+c]-Average);
1100  Count[fff] += 1.;
1101  }
1102  }
1103  }
1104  char FileName[256];
1105  if(What == 0)
1106  sprintf(FileName,"DirDecoupling.dat");
1107  else if(What == 1)
1108  sprintf(FileName,"E2EDecoupling.dat");
1109  FILE *FileToWrite = fopen(FileName,"w");
1110  char cSystem[STRSIZE];
1111  SysDef(cSystem);
1112  // fprintf(FileToWrite,"%s",cSystem);
1113  fprintf(FileToWrite,"#Time AngleProd\n");
1114  for(int f=1;f<NFileTot;f++){
1115  double Inv = Count[f] > 0. ? Count[f] : 1.;
1116  fprintf(FileToWrite,"%lf %lf %lf\n",Time[f],AngleDiff[f]/(Inv*Variance),AngleDiff2[f]/(Variance*pNChain()) );
1117  }
1118  fclose(FileToWrite);
1119  for(int n=0;n<pNNano();n++){
1120  sprintf(FileName,"NanoAngles%d.dat",n);
1121  FILE *FNanoAng = fopen(FileName,"w");
1122  for(int f=0;f<NFileTot;f++)
1123  fprintf(FNanoAng,"%lf %lf \n",Time[f],AnglesNano[f*pNNano()+n]*360./DUE_PI);
1124  fclose(FNanoAng);
1125  }
1126  // FileToWrite = fopen("RawAngles.dat","w");
1127  // for(int f=0;f<NFileTot;f++){
1128  // for(int c=0;c<pNChain();c++){
1129  // fprintf(FileToWrite,"%lf\n",Angles[f*pNChain()+c]);
1130  // }
1131  // }
1132  // fclose(FileToWrite);
1133  free(Count);
1134  free(Time);
1135  free(Angles);
1136  free(AngleDiff);
1137  free(AngleDiff2);
1138 }
1146 void ElPoly::AreaCompr(int NSample){
1147  int SubDiv[3] = {NSample,NSample,1};
1148  int NEdge = (SubDiv[CLat1]*SubDiv[CLat2]);
1149  //mean StdDev Area*StdDev/Mean^2
1150  double *ChainPArea = (double *)calloc(3*NSample,sizeof(double));
1151  double *Plot = (double *)calloc(NEdge,sizeof(double));
1152  double *Area = (double *)calloc(NSample,sizeof(double));
1153  double *NDiv = (double *)calloc(NSample,sizeof(double));
1154  char FileName[256];
1155  //for the simulation in NPtT a smaller box is used
1156  double Edge[3] = {pEdge(0)-1.,pEdge(1)-1.,pEdge(2)};
1157  double InvEdge[3] = {1./Edge[0],1./Edge[1],1./Edge[2]};
1158  for(int f=NFile[0];f<NFile[1];f++){
1159  Processing(f);
1160  OpenRisk(cFile[f],BF_CHAIN);
1161  for(int t=0;t<NSample;t++){
1162  SubDiv[CLat1] = NSample-t;
1163  SubDiv[CLat2] = NSample-t;
1164  if(t == NSample -1) SubDiv[CLat2] = 2;
1165  NDiv[t] = SubDiv[CLat1]*SubDiv[CLat2];
1166  Area[t] = Edge[CLat1]*Edge[CLat2]/(SubDiv[CLat1]*SubDiv[CLat2]);
1167  for(int n=0;n<NEdge;n++) Plot[n] = 0.;
1168  for(int c=0;c<pNChain();c++){
1169  double Posx = pChPos(c,CLat1) - floor(pChPos(c,CLat1)*pInvEdge(CLat1))*pEdge(CLat1);
1170  double Posy = pChPos(c,CLat2) - floor(pChPos(c,CLat2)*pInvEdge(CLat2))*pEdge(CLat2);
1171  int vx = (int)(Posx*InvEdge[CLat1]*SubDiv[CLat1]);
1172  int vy = (int)(Posy*InvEdge[CLat2]*SubDiv[CLat2]);
1173  if(vx < 0 || vx >= SubDiv[CLat1]) continue;
1174  if(vy < 0 || vy >= SubDiv[CLat2]) continue;
1175  Plot[vx*NSample+vy] += 1.;
1176  }
1177  double SumX = 0.;
1178  double SumX2 = 0.;
1179  double Count = 0.;
1180  for(int sx = 0;sx<SubDiv[CLat1];sx++){
1181  for(int sy = 0;sy<SubDiv[CLat2];sy++){
1182  int n = sx*NSample+sy;
1183  Plot[n] /= Area[t];
1184  SumX += Plot[n];
1185  Count += 1.;
1186  }
1187  }
1188  double NInv = Count > 0. ? 1./Count : 1.;
1189  double MeanVal = SumX*NInv;
1190  for(int sx = 0;sx < SubDiv[CLat1];sx++){
1191  for(int sy = 0;sy < SubDiv[CLat2];sy++){
1192  int n = sx*NSample+sy;
1193  SumX2 += SQR(Plot[n] - MeanVal);
1194  }
1195  }
1196  ChainPArea[3*t ] += MeanVal;
1197  ChainPArea[3*t+1] += SumX2/(Count - 1.);
1198  ChainPArea[3*t+2] += 1.;
1199  if(1==0){//writes the single distribution files
1200  sprintf(FileName,"Distr%.0fKappa%02d.dat",pkappaN(),t);
1201  FILE *FileToWrite = fopen(FileName,"a");
1202  for(int l1=0;l1<SubDiv[CLat1];l1++){
1203  for(int l2=0;l2<SubDiv[CLat2];l2++){
1204  fprintf(FileToWrite,"%lf\n",Plot[l1*NSample+l2]);
1205  }
1206  }
1207  fclose(FileToWrite);
1208  }
1209  }
1210  }
1211 #ifdef OMPI_MPI_H
1212  MPI_Allreduce(MPI_IN_PLACE,ChainPArea,3*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
1213  int Rank=0;
1214  MPI_Comm_rank(Proc->CommGrid, &Rank);
1215  if(Rank==0){
1216 #endif
1217  printf("\n");
1218  free(Plot);
1219  //--------normalizing--------------------
1220  double Norm = 1./(double)(NFile[1]-NFile[0]);
1221  for(int t=0;t<NSample;t++){
1222  ChainPArea[3*t ] *= Norm;
1223  ChainPArea[3*t+1] *= Norm;
1224  ChainPArea[3*t+2] *= Norm;
1225  }
1226  sprintf(FileName,"ChainPArea%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN());
1227  FILE *FileToWrite = fopen(FileName,"w");
1228  char cSystem[STRSIZE];
1229  SysDef(cSystem);
1230  fprintf(FileToWrite,"# %s",cSystem);
1231  fprintf(FileToWrite,"# Area MeanChPArea SDeviation Area*Dev/MeanChPArea\n");
1232  for(int t=0;t<NSample;t++){
1233  double Compr = Area[t]*ChainPArea[3*t+1]/SQR(ChainPArea[3*t ]);
1234  fprintf(FileToWrite,"%lf %lf %lf %lf\n",Area[t],(ChainPArea[3*t ]),ChainPArea[3*t+1],Compr);
1235  //fprintf(FileToWrite,"%lf %lf %lf %lf\n",ChainPArea[3*t+0],Area[t]*ChainPArea[3*t+1],ChainPArea[3*t+1],ChainPArea[3*t+2]);
1236  }
1237  fclose(FileToWrite);
1238 #ifdef OMPI_MPI_H
1239  }
1240 #endif
1241  free(ChainPArea);
1242 }
1248 void ElPoly::ElasticCoupling(int NSample){
1249  double Edge[3] = {pEdge(0)-1.,pEdge(1)-1.,pEdge(2)};
1250  int SubDiv[3] = {NSample,NSample,1};
1251  //mean StdDev Area*StdDev/Mean^2
1252  double *ElCoup = (double *)calloc(3*(NSample+1),sizeof(double));
1253  double *Area = (double *)calloc((NSample+1),sizeof(double));
1254  double *NDiv = (double *)calloc((NSample+1),sizeof(double));
1255  double MeanSigma = 0.;
1256  char FileName[256];
1257  //for the simulation in NPtT a smaller box is used
1258  double InvEdge[3] = {1./Edge[0],1./Edge[1],1./Edge[2]};
1259  double *PosUp = (double *)calloc(SQR(NSample),sizeof(double));
1260  double *PosDown = (double *)calloc(SQR(NSample),sizeof(double));
1261  double *CountUp = (double *)calloc(SQR(NSample),sizeof(double));
1262  double *CountDown= (double *)calloc(SQR(NSample),sizeof(double));
1263  for(int f=NFile[0];f<NFile[1];f++){
1264  Processing(f);
1265  OpenRisk(cFile[f],BF_CHAIN);
1266  for(int t=0;t<NSample+1;t++){
1267  //defining the partitions
1268  // SubDiv[CLat1] = (int)pow(2.,t);
1269  // SubDiv[CLat2] = (int)pow(2.,t);
1270  SubDiv[CLat1] = NSample - t;
1271  SubDiv[CLat2] = NSample - t;
1272  if(t == NSample-1) SubDiv[CLat2] = 2;
1273  if(t == NSample){
1274  SubDiv[CLat1] = 1;
1275  SubDiv[CLat2] = 1;
1276  }
1277  NDiv[t] = SubDiv[CLat1]*SubDiv[CLat2];
1278  Area[t] = Edge[CLat1]*Edge[CLat2]/(double)(SubDiv[CLat1]*SubDiv[CLat2]);
1279  for(int n=0;n<SQR(NSample);n++){
1280  PosUp[n] = 0.;
1281  PosDown[n] = 0.;
1282  CountUp[n] = 0.;
1283  CountDown[n] = 0.;
1284  }
1285  for(int p=0;p<pNPart();p++){
1286  if(pType(p) != 1) continue;
1287  double x = pPos(p,CLat1) - floor(pPos(p,CLat1)*pInvEdge(CLat1))*pEdge(CLat1);
1288  double y = pPos(p,CLat2) - floor(pPos(p,CLat2)*pInvEdge(CLat2))*pEdge(CLat2);
1289  int sx = (int)(x*InvEdge[CLat1]*SubDiv[CLat1]);
1290  if(sx < 0 || sx >= SubDiv[CLat1]) continue;
1291  int sy = (int)(y*InvEdge[CLat2]*SubDiv[CLat2]);
1292  if(sy < 0 || sy >= SubDiv[CLat2]) continue;
1293  if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){
1294  PosUp[sx*NSample+sy] += pPos(p,CNorm);
1295  CountUp[sx*NSample+sy] += 1.;
1296  }
1297  else {
1298  PosDown[sx*NSample+sy] += pPos(p,CNorm);
1299  CountDown[sx*NSample+sy] += 1.;
1300  }
1301  }
1302  double SumX = 0.;
1303  double SumX2 = 0.;
1304  double Count = 0.;
1305  for(int sx = 0;sx<SubDiv[CLat1];sx++){
1306  for(int sy = 0;sy<SubDiv[CLat2];sy++){
1307  int n = sx*NSample+sy;
1308  PosUp[n] /= CountUp[n] > 1. ? CountUp[n] : 1.;
1309  PosDown[n] /= CountDown[n] > 1. ? CountDown[n] : 1.;
1310  double Val = (PosUp[n]-PosDown[n]);
1311  SumX += Val;
1312  Count += 1.;
1313  ElCoup[3*t+2] += 1.;
1314  }
1315  }
1316  double NInv = Count > 0. ? 1./Count : 1.;
1317  double MeanVal = SumX*NInv;
1318  for(int sx = 0;sx < SubDiv[CLat1];sx++){
1319  for(int sy = 0;sy < SubDiv[CLat2];sy++){
1320  int n = sx*NSample+sy;
1321  double Val = (PosUp[n]-PosDown[n]);
1322  SumX2 += SQR(Val - MeanVal);
1323  }
1324  }
1325  ElCoup[3*t ] += MeanVal;
1326  ElCoup[3*t+1] += SumX2/(Count - 1.);
1327  if(1==0){//writes the single distribution files
1328  sprintf(FileName,"Distr%.0fKappa%02d.dat",pkappaN(),t);
1329  FILE *FileToWrite = fopen(FileName,"a");
1330  for(int l1=0;l1<SubDiv[CLat1];l1++){
1331  for(int l2=0;l2<SubDiv[CLat2];l2++){
1332  int n = l1*NSample+l2;
1333  fprintf(FileToWrite,"%lf\n",(PosUp[n]-PosDown[n]));
1334  }
1335  }
1336  fclose(FileToWrite);
1337  }
1338  }
1339  }
1340  #ifdef OMPI_MPI_H
1341  MPI_Allreduce(MPI_IN_PLACE,ElCoup,3*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
1342  int Rank=0;
1343  MPI_Comm_rank(Proc->CommGrid, &Rank);
1344  if(Rank==0){
1345 #endif
1346  printf("\n");
1347  //--------normalizing--------------------
1348  double Norm = 1./(double)(NFile[1]-NFile[0]);
1349  for(int t=0;t<NSample+1;t++){
1350  ElCoup[3*t ] *= Norm;
1351  ElCoup[3*t+1] *= Norm;
1352  ElCoup[3*t+2] *= Norm;
1353  }
1354  ElCoup[3*NSample+1] = 0.;
1355  sprintf(FileName,"ElasticCoupling%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN());
1356  FILE *FileToWrite = fopen(FileName,"w");
1357  char cSystem[STRSIZE];
1358  SysDef(cSystem);
1359  fprintf(FileToWrite,"# %s",cSystem);
1360  fprintf(FileToWrite,"# Area Thickness ThickStdDev Area*ThickStdDev\n");
1361  for(int t=0;t<NSample+1;t++){
1362  fprintf(FileToWrite,"%lf %lf %lf %lf \n",Area[t],ElCoup[3*t ],ElCoup[3*t+1],Area[t]*ElCoup[3*t+1]/SQR(ElCoup[3*t ]));
1363  }
1364  fclose(FileToWrite);
1365 #ifdef OMPI_MPI_H
1366  }
1367 #endif
1368  free(ElCoup);
1369  free(PosUp);
1370  free(PosDown);
1371  free(CountUp);
1372  free(CountDown);
1373 }
1378  //mean StdDev Area*StdDev/Mean^2
1379  double ElCoup[3] = {0.,0.,0.};
1380  char FileName[256];
1381  char cSystem[STRSIZE];
1382  //for the simulation in NPtT a smaller box is used
1383  sprintf(FileName,"ElCoupNVT%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN());
1384  FILE *FileToWrite = fopen(FileName,"w");
1385  SysDef(cSystem);
1386  fprintf(FileToWrite,"# %s",cSystem);
1387  fprintf(FileToWrite,"# step Thickness ThickStdDev\n");
1388  for(int f=NFile[0];f<NFile[1];f++){
1389  Processing(f);
1390  OpenRisk(cFile[f],BF_CHAIN);
1391  double PosUp = 0.;
1392  double PosDown = 0.;
1393  double CountUp = 0.;
1394  double CountDown = 0.;
1395  for(int p=0;p<pNPart();p++){
1396  if(pType(p) != 1) continue;
1397  if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){
1398  PosUp += pPos(p,CNorm);
1399  CountUp += 1.;
1400  }
1401  else {
1402  PosDown += pPos(p,CNorm);
1403  CountDown += 1.;
1404  }
1405  }
1406  PosUp /= CountUp > 1. ? CountUp : 1.;
1407  PosDown /= CountDown > 1. ? CountDown : 1.;
1408  ElCoup[0] += (PosUp-PosDown);
1409  ElCoup[1] += SQR(PosUp-PosDown);
1410  fprintf(FileToWrite,"%d %lf\n",f,PosUp-PosDown);
1411  }
1412  printf("\n");
1413  ElCoup[1] = sqrt((ElCoup[1] - SQR(ElCoup[0])*(NFile[1]-NFile[1]))/((double)(NFile[0]-NFile[1])-1.));
1414  fprintf(FileToWrite,"# %lf %lf %lf\n",pEdge(CLat1)*pEdge(CLat2),ElCoup[0],ElCoup[1]);
1415  fclose(FileToWrite);
1416 }
1417 void ElPoly::BilayerDistance(char *OutFile,int NSample){
1418  FILE *File2Write = fopen(OutFile,"w");
1419  int NLayer = 2;
1420  double PosBf[3];
1421  for(int f=NFile[0];f<NFile[1];f++){
1422  Processing(f);
1423  if(OpenRisk(cFile[f],BF_CHAIN))return;
1424  double *PosLay = (double *)calloc(4*NSample*NSample,sizeof(double));
1425  double *CountLay = (double *)calloc(4*NSample*NSample,sizeof(double));
1426  for(int b=0,cOff=0,pOff=0;b<pNBlock();cOff+=pNChain(b++)){
1427  int Level = 0;
1428  if(!strcmp(Block[b].Name,"LIPID1")) Level = 2;
1429  for(int c=cOff;c<pNChain(b)+cOff;c++,pOff+=pNPCh(b)){
1430  for(int d=0;d<3;d++){
1431  PosBf[d] = Ch[c].Pos[d] - floor(Ch[c].Pos[d]*pInvEdge(d))*pEdge(d);
1432  }
1433  int cLevel = Level;
1434  if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_UP)) cLevel += 1;
1435  int sx = (int)(PosBf[CLat1]*pInvEdge(CLat1)*NSample);
1436  int sy = (int)(PosBf[CLat2]*pInvEdge(CLat2)*NSample);
1437  if(sx < 0 || sx >= NSample) continue;
1438  if(sy < 0 || sy >= NSample) continue;
1439  PosLay[(sx*NSample+sy)*4+cLevel] += Ch[c].Pos[CNorm];//PosBf[CNorm];
1440  CountLay[(sx*NSample+sy)*4+cLevel] += 1.;
1441  }
1442  }
1443  for(int sx = 0;sx < NSample;sx++){
1444  for(int sy = 0;sy < NSample;sy++){
1445  for(int l=0;l<4;l++){
1446  PosLay[(sx*NSample+sy)*4+l] /= CountLay[(sx*NSample+sy)*4+l] > 0. ? CountLay[(sx*NSample+sy)*4+l] : 1.;
1447  }
1448  double h1 = PosLay[(sx*NSample+sy)*4+2]-PosLay[(sx*NSample+sy)*4+1];
1449  double h2 = PosLay[(sx*NSample+sy)*4+3]-PosLay[(sx*NSample+sy)*4+0];
1450  h2 -= floor(h2*pInvEdge(CNorm))*pEdge(CNorm);
1451  //if(PosLay[(sx*NSample+sy)*4+0]*PosLay[(sx*NSample+sy)*4+1]*PosLay[(sx*NSample+sy)*4+2]*PosLay[(sx*NSample+sy)*4+3] <= 0.) continue;
1452  fprintf(File2Write,"%lf %lf %lf %lf %lf\n",MIN(h1,h2),PosLay[(sx*NSample+sy)*4+0],PosLay[(sx*NSample+sy)*4+1],PosLay[(sx*NSample+sy)*4+2],PosLay[(sx*NSample+sy)*4+3]);
1453  }
1454  }
1455  free(PosLay);
1456  free(CountLay);
1457  }
1458  printf("\n");
1459  fclose(File2Write);
1460 }
1461 void ElPoly::BondDistr(char *OutFile,int NBin){
1462  FILE *File2Write = fopen(OutFile,"w");
1463  // double **Histo = (double **)calloc(3,sizeof(double));
1464  // double **Raw = (double **)calloc(3,sizeof(double));
1465  // for(int d=0;d<3;d++){
1466  // Histo[d] = (double *)calloc(NBin,sizeof(double));
1467  // Raw[d] = (double *)calloc(pNPart,sizeof(double));
1468  // }
1469  for(int f=NFile[0];f<NFile[1];f++){
1470  Processing(f);
1471  if(OpenRisk(cFile[f],BF_NO)) return;
1472  double Dist[3];
1473  for(int c=0;c<pNChain();c++){
1474  for(int p=0;p<pNPCh()-1;p++){
1475  int p1 = p + c*pNPCh();
1476  for(int d=0;d<3;d++){
1477  Dist[d] = remainder(pPos(p1,d) - pPos(p1+1,d),pEdge(d));
1478  }
1479  fprintf(File2Write,"%lf %lf %lf %lf\n",Dist[0],Dist[1],Dist[2],sqrt(SQR(Dist[0])+SQR(Dist[1])+SQR(Dist[2])));
1480  }
1481  }
1482  }
1483  // MOMENTI Mom[3];
1484  // int xMin = Raw[0][0];
1485  // int xMax = Raw[0][0];
1486  // for(int d=0;d<3;d++){
1487  // Mom[d] = Mat->Distribuzione(Raw[d],pNPart-pNChain,Histo[d],NBin);
1488  // Mat->Normalize(Histo[d],NBin);
1489  // if(xMin > Mom[d].Minimo) xMin = Mom[d].Minimo;
1490  // if(xMax > Mom[d].Massimo) xMax = Mom[d].Massimo;
1491  // }
1492  // for(int v=0;v<NBin;v++){
1493  // double x = v/(double)(NBin)*(xMax-xMin)+xMin;
1494  // fprintf(File2Write,"%lf %lf %lf %lf\n",x,Histo[0][v],Histo[1][v],Histo[2][v]);
1495  // }
1496  // for(int d=0;d<3;d++){
1497  // free(Histo[d]);
1498  // free(Raw[d]);
1499  // }
1500  // free(Histo);
1501  // free(Raw);
1502  fclose(File2Write);
1503 }
1505  double *EndToEnd = (double *) calloc(pNPart()*(NFile[1]-NFile[0]),sizeof(double));
1506  double *Distr = (double *) calloc((NFile[1]-NFile[0]),sizeof(double));
1507  for(int f=NFile[0];f<NFile[1];f++){
1508  Processing(f);
1509  if(OpenRisk(cFile[f],BF_NO))return;
1510 
1511 
1512  }
1513 
1514 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
void SpectrumMidplane(int NBin)
Calculates the 1d, 2d spectrum spe.
double pTime()
Total time.
Definition: VarData.h:908
int Diffusivity()
Diffusivity coefficient of the chains diff.
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.
double pChPos(int p, int d)
Return back folded position.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void Midplane(int NBin)
Midplanes for a sequence of snapshots.
int NFileTot
Total number of file.
Definition: ElPoly.h:420
NANO * Nano
Extra particle.
Definition: VarData.h:1044
int PairCorrelationPep(double **Point, int NSample, int Type)
2-d pair correlation on a square fererring to the pep position
Definition: VarDataExp.cpp:105
void BilayerDistance(char *OutFile, int NBin)
Measure the elastic coupling between the two sheets.
int SetNChain(int NewNCh)
Set and reallocate the number of chains.
Definition: VarDataComm.cpp:78
int pChain(int p)
Return the chain.
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
int Asym
Diblock limit of the chain.
Definition: VarData.h:273
void SetPos(int p, double *Pos)
Set the particle position.
int BfDefChain()
Definition of the chain.
double Rad
Size.
Definition: VarData.h:445
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
int Scattering2d(double **Point, int NSample, int Type)
2-d Scattering
Definition: VarDataExp.cpp:157
double Due
Second moment.
double Height
Height of the cylinder.
Definition: VarData.h:449
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
void Decoupling(int What)
Decoupling of the direction/end to end distance of the chains.
int ScatteringF(int NBin, int How)
2d Scattering scatt scatt2
void EndToEndDist()
Distribution of the end to end distances.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
void AreaCompr(int NSample)
Calculates the number of chain per area area.
MOMENTI Distribuzione(const double *st, int NMass)
Moments of a signal.
void PairCorr(int NBin, int NDim)
1d pair correlation
double Area
Area of a pore or a stalk.
Definition: VarData.h:471
void HeaderAverage(int nNano)
Read the header and average the information.
int Angle(int NBin)
Calculation of the contact angle angle Boh?
Definition: ElPolyEl.cpp:78
double Pos[3]
Position.
Definition: VarData.h:427
void ReadHeader(FILE *FileToRead)
Reads a header.
void WidomIn()
Create a set of files with a chain more.
double pNanoPos(int n, int d)
Return back folded nano position.
int SetNPart(int NewNPart)
Set and reallocate the number of particles.
Definition: VarDataComm.cpp:2
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.
Moments of a distribution.
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
int NChainPSquareF()
Calculates the number of chain per area area.
int pNBlock()
Number of blocks.
double pchiN()
Incompatibility.
Definition: VarData.h:928
bool WriteXvt(char *OutFile)
Writes a "system-file" or a "x y z" file".
int SpectrumF(int NBin)
Calculates the 1d, 2d spectrum spe.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double Casuale()
Random uniform number.
void SetType(int p, int t)
Set the particle type.
void ElasticCouplingNVT()
Measure the elastic coupling between the two sheets.
void DiffSlab(int NSlab)
Diffusivity of particles starting from an initial slab.
void WidomOut()
Create NCh systems with a chain less.
int pNNano()
Number of nanoparticles.
int NBackFold
Type of backfold.
Definition: ElPoly.h:439
void SwapChain(int c1, int c2, int b)
Swap two chains.
Definition: VarDataEl.cpp:355
void Spettro2d(double *st, double *sw, int NMass)
Compute the 2d spectrum of.
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int PairCorrelation(double *Point, int NSample, int How, int Type)
1-d pair correlation
Definition: VarDataExp.cpp:24
void End2EndDistr(char *OutFile)
Write the end to end distance of the chains.
int PairCorrelationF(int NBin, int How)
1d, 2d pair correlation of the chains pairRound pairSquare (obsolete)
double Gaussiano(double Media, double Scarto)
Gaussian random number.
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
double Uno
First moment.
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
void ElasticCoupling(int NSample)
Measure the elastic coupling between the two sheets.
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
Information of every particle.
Definition: VarData.h:214
int PairCorrelationSquare(double **Point, int NSample, int Type)
2-d pair correlation on a square
Definition: VarDataExp.cpp:73
Domain decomposition as pointer to linked particles.
Definition: Cubo.h:162
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
void Set(double Val, int Col)
Set the N column.
int PairCorrelationRound(double **Point, int NSample, int Type)
Circular 2-d pair correlation.
Definition: VarDataExp.cpp:122
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 pkappaN()
Incompressibility.
Definition: VarData.h:930
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.