|
Allink
v0.1
|
00001 #include "ElPoly.h" 00002 00003 int ElPoly::Diffusivity(){ 00004 int NDim = 2; 00005 NFileTot = NFile[1] - NFile[0]; 00006 int NChain = 0; 00007 for(int b=0;b<pNBlock();b++){ 00008 if(!strcmp(Block[b].Name,"ADDED")) continue; 00009 if(!strcmp(Block[b].Name,"CHOL")) continue; 00010 NChain += Block[b].NChain; 00011 } 00012 double *FileAverage = (double *)calloc(NDim*NChain*NFileTot,sizeof(double)); 00013 double *ChainDiff = (double *)calloc(NFileTot,sizeof(*ChainDiff)); 00014 double *ChainOrigin = (double *)calloc(NFileTot,sizeof(*ChainDiff)); 00015 double *Count = (double *)calloc(NFileTot,sizeof(*Count)); 00016 double *Time = (double *)calloc(NFileTot,sizeof(*Time)); 00017 double CmInit[3] = {pCm(0),pCm(1),pCm(2)}; 00018 double InitTime = pTime(); 00019 for(int f=NFile[0];f<NFile[1];f++){ 00020 Processing(f); 00021 if(OpenRisk(cFile[f],BF_NO))return 1; 00022 BfDefChain(); 00023 Time[f] = pTime() - InitTime; 00024 if(Time[f]<=0.) Time[f] = (double)f; 00025 for(int c=0;c<NChain;c++){ 00026 int ff = f - NFile[0]; 00027 for(int d=0;d<NDim;d++){ 00028 FileAverage[(c*NFileTot+ff)*NDim+d] = Ch[c].Pos[d]; 00029 } 00030 } 00031 } 00032 printf("\n"); 00033 #ifdef OMPI_MPI_H 00034 MPI_Allreduce(MPI_IN_PLACE,FileAverage,NDim*NChain*NFileTot, MPI_DOUBLE, MPI_SUM, Proc->CommGrid); 00035 int Rank=0; 00036 MPI_Comm_rank(Proc->CommGrid, &Rank); 00037 if(Rank==0){ 00038 #endif 00039 for(int f=0;f<NFileTot;f++){ 00040 for(int c=0;c<NChain;c++){ 00041 for(int d=0;d<NDim;d++){ 00042 ChainOrigin[f] += 00043 SQR((FileAverage[(c*NFileTot)*NDim+d]-FileAverage[(c*NFileTot+f)*NDim+d]));//-(pCm(d)-CmInit[d]) ); 00044 } 00045 } 00046 } 00047 for(int f=0;f<NFileTot;f++){ 00048 for(int ff = f; ff<NFileTot;ff++){ 00049 int fff = ff - f; 00050 for(int c=0;c<NChain;c++){ 00051 for(int d=0;d<NDim;d++){ 00052 ChainDiff[fff] += 00053 SQR((FileAverage[(c*NFileTot+ff)*NDim+d]-FileAverage[(c*NFileTot+f)*NDim+d]));//-(pCm(d)-CmInit[d]) ); 00054 } 00055 Count[fff] += 1.; 00056 } 00057 //printf("%d %d %d %lf\n",f,ff,fff,ChainDiff[fff]); 00058 } 00059 } 00060 char FileName[60]; 00061 sprintf(FileName,"Diffusivity.dat"); 00062 FILE *FileToWrite = fopen(FileName,"w"); 00063 char cSystem[STRSIZE]; 00064 SysDef(cSystem); 00065 fprintf(FileToWrite,"%s",cSystem); 00066 fprintf(FileToWrite,"#Time ChainOrigin OriginSum DiffSum ChainDiff Count\n"); 00067 for(int f=0;f<NFileTot;f++){ 00068 double Inv = Count[f] > 0. ? Count[f] : 1.; 00069 fprintf(FileToWrite,"%lf %.5g %.5g %.5g\n",Time[f],ChainOrigin[f]/(double)(NChain),ChainDiff[f]/Inv,ChainDiff[f]/(Inv*Time[f]) ); 00070 } 00071 #ifdef OMPI_MPI_H 00072 } 00073 #endif 00074 free(FileAverage); 00075 free(ChainDiff); 00076 free(ChainOrigin); 00077 free(Count); 00078 free(Time); 00079 return 0; 00080 } 00081 void ElPoly::DiffSlab(int NSlab){ 00082 int NBin = NFile[1] - NFile[0]; 00083 double Volume1 = 0.; 00084 double HeightCyl = pEdge(CNorm); 00085 double LenMin = MIN(pEdge(CLat1),pEdge(CLat2)); 00086 double LenMax = MAX(pEdge(CLat1),pEdge(CLat2)); 00087 double BoxRad = .5*MIN(pEdge(CLat1),pEdge(CLat2)); 00088 double InvBoxRad = 1./BoxRad; 00089 double *OldPos = (double *)calloc(3*pNPart(),sizeof(double)); 00090 double *Diff = (double *)calloc(NBin*pNPart(),sizeof(double)); 00091 double *Count = (double *)calloc(NSlab,sizeof(double)); 00092 int *PartOr = (int *)calloc(pNPart(),sizeof(int)); 00093 double dr[3]; 00094 for(int p = 0;p<pNPart();p++){ 00095 double dr2 = 0.; 00096 for(int d=0;d<3;d++){ 00097 OldPos[3*p+d] = Pm[p].Pos[d]; 00098 dr[d] = Nano->Pos[d] - Pm[p].Pos[d]; 00099 dr[d] -= floor(dr[d]*pInvEdge(d))*pEdge(d); 00100 } 00101 dr2 = sqrt(SQR(dr[CLat1]) + SQR(dr[CLat2])); 00102 int vr = (int)(dr2/BoxRad*NSlab); 00103 PartOr[p] = vr; 00104 if(vr < 0 || vr >= NSlab) continue; 00105 Count[vr] += 1.; 00106 } 00107 char FName[60]; 00108 if(1==1){//particles position and slab circles 00109 FILE *FRepr = fopen("SystemSlab.dat","w"); 00110 for(int p=0;p<pNPart();p++){ 00111 fprintf(FRepr,"%lf %lf %lf 0\n",Pm[p].Pos[CLat1],Pm[p].Pos[CLat2],.5); 00112 } 00113 for(int v=0;v<NSlab;v++){ 00114 double Dist = v*BoxRad/(double)NSlab; 00115 for(int i=0;i<100;i++){ 00116 double ang = i*.01*2.*M_PI; 00117 double x = Dist*cos(ang) + Nano->Pos[CLat1]; 00118 double y = Dist*sin(ang) + Nano->Pos[CLat2]; 00119 fprintf(FRepr,"%lf %lf %lf 2\n",x,y,.51); 00120 } 00121 } 00122 fclose(FRepr); 00123 double Norm = 1./(double)pNPart(); 00124 for(int v=0;v<NSlab;v++){ 00125 double Dist = v*BoxRad/(double)NSlab; 00126 sprintf(FName,"MeanSqDispl%.3f.dat",Dist); 00127 FILE *FWrite = fopen(FName,"w"); 00128 fclose(FWrite); 00129 } 00130 } 00131 int NOut = 0; 00132 for(int f=NFile[0];f<NFile[1];f++){ 00133 Processing(f); 00134 if(OpenRisk(cFile[f],BF_PART)); 00135 double OutRatio = 100.*NOut/(double)(f*pNPart()); 00136 fprintf(stderr,"step %d)/%d fraction of number of beads out of the slab %lf %%\r",f,NBin,OutRatio); 00137 double Rad2Mean = 0.; 00138 for(int p = 0;p<pNPart();p++){ 00139 int vr1 = PartOr[p]; 00140 if(vr1 < 0 || vr1 >= NSlab) continue; 00141 double dr2 = 0.; 00142 for(int d=0;d<3;d++){ 00143 dr[d] = Nano->Pos[d] - Pm[p].Pos[d]; 00144 dr[d] -= floor(dr[d]*pInvEdge(d))*pEdge(d); 00145 } 00146 dr2 = sqrt(SQR(dr[CLat1]) + SQR(dr[CLat2])); 00147 int vr = (int)(dr2*InvBoxRad*NSlab); 00148 if(vr != vr1){ 00149 NOut++; 00150 } 00151 double Rad2 = 0.; 00152 for(int d=0;d<2;d++){ 00153 double Dist = Pm[p].Pos[d] - OldPos[3*p+d]; 00154 Dist -= floor(Dist*pInvEdge(d))*pEdge(d); 00155 Rad2 += SQR(Dist); 00156 } 00157 // vr or vr1? 00158 Diff[f*NSlab+vr1] += Rad2; 00159 Rad2Mean += Rad2; 00160 } 00161 if(1==1){//particle position and mean displacement 00162 Rad2Mean = pNPart()/Rad2Mean; 00163 sprintf(FName,"DisplPos%05d.dat",f); 00164 FILE *FWrite = fopen(FName,"w"); 00165 fprintf(FWrite,"# l(%lf %lf %lf) d[part]\n",pEdge(CLat1),pEdge(CLat2),.5); 00166 for(int p=0;p<pNPart();p++){ 00167 double Rad2 = 0.; 00168 for(int d=0;d<2;d++){ 00169 double Dist = Pm[p].Pos[d] - OldPos[3*p+d]; 00170 Dist -= floor(Dist*pInvEdge(d))*pEdge(d); 00171 Rad2 += SQR(Dist); 00172 } 00173 int vr1 = PartOr[p]; 00174 double Vel = Rad2*Rad2Mean; 00175 fprintf(FWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)\n",Pm[p].Pos[CLat1],Pm[p].Pos[CLat2],Vel,Vel,Vel,Vel); 00176 } 00177 fclose(FWrite); 00178 } 00179 } 00180 printf("\n"); 00181 double DisplMean = 0.; 00182 //normalize 00183 for(int v=0;v<NSlab;v++){ 00184 Count[v] = Count[v] > 0. ? 1./Count[v] : 1.; 00185 for(int b=0;b<NBin;b++){ 00186 Diff[b*NSlab+v] *= Count[v]; 00187 } 00188 } 00189 for(int v=0;v<NSlab;v++){ 00190 DisplMean += Diff[(NBin-1)*NSlab+v] ; 00191 double Dist = v*BoxRad/(double)NSlab; 00192 sprintf(FName,"MeanSqDispl%.3f.dat",Dist); 00193 FILE *FWrite = fopen(FName,"w"); 00194 for(int b=1;b<NBin;b++){ 00195 fprintf(FWrite,"%d %lf\n",b,Diff[b*NSlab+v]); 00196 } 00197 fclose(FWrite); 00198 } 00199 DisplMean /= (double) NSlab; 00200 if(1==1){//slab mean displ 00201 sprintf(FName,"SlabMeadDispl.dat"); 00202 FILE *FWrite = fopen(FName,"w"); 00203 fprintf(FWrite,"# l(%lf %lf %lf) d[part]\n",pEdge(CLat1),pEdge(CLat2),.5); 00204 for(int sx=0;sx<100;sx++){ 00205 for(int sy=0;sy<100;sy++){ 00206 double x = sx*.01*pEdge(CLat1); 00207 double y = sy*.01*pEdge(CLat2); 00208 double Rad2 = sqrt(SQR(x-Nano->Pos[CLat1])+SQR(y-Nano->Pos[CLat2])); 00209 int vr = (int)(Rad2*InvBoxRad*NSlab); 00210 if(vr < 0 || vr >= NSlab) continue; 00211 double Vel = Diff[(NBin-1)*NSlab+vr]; 00212 fprintf(FWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)\n",x,y,Vel,Vel,Vel,Vel); 00213 } 00214 } 00215 fclose(FWrite); 00216 } 00217 free(OldPos); 00218 free(Diff); 00219 free(Count); 00220 free(PartOr); 00221 } 00222 #include <Cubo.h> 00223 typedef DdLinkedList DomDec; 00224 void ElPoly::PairCorr(int NBin,int NDim){ 00225 int NAlloc = NBin; 00226 double CutOff = 4.; 00227 if(NDim == 2) NAlloc = NBin*NBin; 00228 double *Prob = (double *)calloc(NAlloc,sizeof(double)); 00229 double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)}; 00230 double Norm = 0.; 00231 for(int f=NFile[0];f<NFile[1];f++){ 00232 Processing(f); 00233 // DomDec *Pc = new DomDec(Edge,pNPart(),1.); 00234 // for(int p=0;p<pNPart();p++) Pc->AddPart(p,Pm[p].Pos); 00235 // int NeiList[27]; 00236 double DistRel[4]; 00237 for(int p1=0;p1<pNPart();p1++){ 00238 for(int p2=p1+1;p2<pNPart();p2++){ 00239 // int NNei = Pc->GetNei(Pm[p1].Pos,NeiList); 00240 // for(int i=0;i<NNei;i++){ 00241 // int c1 = NeiList[i]; 00242 // for(Pc->SetCounters(c1);Pc->IfItCell(c1);Pc->IncrCurr(c1)){ 00243 // int p2 = Pc->ItCell(c1); 00244 // if(p2 < p1) continue; 00245 if(!TwoPartDist(p1,p2,DistRel,CutOff)) continue; 00246 int r = (int)(DistRel[3]/CutOff*NBin); 00247 if(r < 0 || r >= NBin) continue; 00248 Prob[r] += 1.; 00249 Norm += 1.; 00250 } 00251 } 00252 // } 00253 // delete Pc; 00254 } 00255 printf("\n"); 00256 Norm *= NFile[1] - NFile[0]; 00257 double rInv = 1.*CutOff/(double)NBin; 00258 FILE *FPair = fopen("PairCorrelation.dat","w"); 00259 fprintf(FPair,"%lf %lf \n",0.,Prob[0]/Norm); 00260 for(int r=1;r<NBin;r++){ 00261 double rr = r*CutOff/(double)NBin; 00262 Prob[r] /= Norm*(CUBE(rr)-CUBE(rr-rInv)); 00263 fprintf(FPair,"%lf %lf \n",rr,Prob[r]); 00264 } 00265 fclose(FPair); 00266 free(Prob); 00267 } 00268 int ElPoly::PairCorrelationF(int NBin,int How){ 00269 double **Plot; 00270 Plot = (double **)calloc(NBin,sizeof(double)); 00271 for(int i=0;i<NBin;i++){ 00272 *(Plot+i) = (double *)calloc(NBin,sizeof(double)); 00273 } 00274 double *dDensity; 00275 dDensity = (double *)calloc(NBin,sizeof(double)); 00276 double *dDensity1; 00277 dDensity1 = (double *)calloc(NBin,sizeof(double)); 00278 for(int f=NFile[0];f<NFile[1];f++){ 00279 Processing(f); 00280 //printf("Opening: %s\n",cFile[f]); 00281 if(OpenRisk(cFile[f],BF_CHAIN))return 1; 00282 if(How == 0 || How == 1) 00283 PairCorrelation(dDensity,NBin,How,NChType); 00284 else if (How == 2) 00285 PairCorrelationRound(Plot,NBin,NChType); 00286 else if (How == 3) 00287 PairCorrelationSquare(Plot,NBin,NChType); 00288 else if (How == 4) 00289 PairCorrelationPep(Plot,NBin,NChType); 00290 for(int v=0;v<NBin;v++){ 00291 //printf("%d %lf\n",v,dDensity[v]); 00292 dDensity1[v] += dDensity[v]; 00293 } 00294 } 00295 printf("\n"); 00296 if(How == 0 || How == 1){ 00297 FILE *FileToWrite = fopen("PairCorrelation.dat","w"); 00298 fprintf(FileToWrite,"# RadDist Density \n"); 00299 for(int v=0;v<NBin;v++){ 00300 dDensity1[v] /= NFile[1] - NFile[0]; 00301 fprintf(FileToWrite,"%lf %lf\n",(double)v/(double)NBin*pEdge(3),dDensity1[v]); 00302 } 00303 fclose(FileToWrite); 00304 } 00305 else { 00306 double Max=0.; 00307 double InvNBin=1./(double)NBin; 00308 for(int v=0;v<NBin;v++){ 00309 for(int vv=0;vv<NBin;vv++){ 00310 // Plot[v][vv] /= DUE_PI*( QUAD((Gen->Edge[3]/InvNBin))*(QUAD((v+1)) - QUAD(v)) ); 00311 if(Max < Plot[v][vv]) 00312 Max = Plot[v][vv]; 00313 } 00314 } 00315 FILE *FileToWrite = fopen("PairCorrelation.xvl","w"); 00316 fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",2.*pEdge(3),2.*pEdge(3),pEdge(3),NBin,ChooseDraw(EL_QUAD1)); 00317 int vRef = (int)(NBin/2.); 00318 for(int v=0;v<NBin;v++){ 00319 for(int vv=0;vv<NBin;vv++){ 00320 if(Plot[v][vv] > 0.){ 00321 if(How == 2){ 00322 double PosX = pEdge(3) + InvNBin*(double)v*pEdge(3)*cos((double)vv*InvNBin*DUE_PI); 00323 double PosY = pEdge(3) + InvNBin*(double)v*pEdge(3)*sin((double)vv*InvNBin*DUE_PI); 00324 fprintf(FileToWrite,"{x(%lf %lf %lf)}\n", 00325 PosX,PosY,Plot[v][vv]/Max*pEdge(3)); 00326 } 00327 else{ 00328 double PosX = pEdge(CLat1)*.5 + (v - vRef)*InvNBin*pEdge(CLat1); 00329 double PosY = pEdge(CLat2)*.5 + (vv - vRef)*InvNBin*pEdge(CLat2); 00330 fprintf(FileToWrite,"{x(%lf %lf %lf)}\n", 00331 PosX,PosY,Plot[v][vv]/Max*pEdge(3)); 00332 } 00333 //printf("%d %d %lf %lf %lf\n",v,vv,PosX,PosY,Plot[v][vv]); 00334 } 00335 } 00336 } 00337 fclose(FileToWrite); 00338 } 00339 free(Plot); 00340 free(dDensity); 00341 free(dDensity1); 00342 return 0; 00343 } 00344 int ElPoly::ScatteringF(int NBin,int How){ 00345 double **Plot; 00346 Plot = (double **)calloc(NBin,sizeof(double)); 00347 for(int i=0;i<NBin;i++){ 00348 *(Plot+i) = (double *)calloc(NBin,sizeof(double)); 00349 } 00350 SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3); 00351 for(int f=NFile[0];f<NFile[1];f++){ 00352 Processing(f); 00353 //printf("Opening: %s\n",cFile[f]); 00354 if(OpenRisk(cFile[f],BF_CHAIN))return 1; 00355 Scattering2d(Plot,NBin,NChType); 00356 } 00357 printf("\n"); 00358 double Max=0.; 00359 double InvNBin=1./(double)NBin; 00360 for(int v=0;v<NBin;v++){ 00361 for(int vv=0;vv<NBin;vv++){ 00362 Plot[v][vv] = POS(Plot[v][vv]); 00363 if(Max < Plot[v][vv]) 00364 Max = Plot[v][vv]; 00365 } 00366 } 00367 Max = 1./Max; 00368 FILE *FileToWrite = NULL; 00369 if(How == 0) 00370 FileToWrite = fopen("Scattering.xvl","w"); 00371 else 00372 FileToWrite = fopen("Scattering2.xvl","w"); 00373 fprintf(FileToWrite,"#l(%lf %lf %lf) v(%d) d[%s]\n",2.*pEdge(3),2.*pEdge(3),pEdge(3),NBin,ChooseDraw(EL_QUAD)); 00374 int vRef = NBin/2; 00375 for(int q=0;q<4;q++){ 00376 double Signx = !(q%2) ? 1. : -1.; 00377 double Signy = q > 1 ? 1. : -1.; 00378 for(int vx=0;vx<NBin;vx++){ 00379 for(int vy=0;vy<NBin;vy++){ 00380 //if(Plot[v][vv] > 0.) 00381 { 00382 double PosX = Signx*vx*InvNBin*pEdge(CLat1); 00383 double PosY = Signy*vy*InvNBin*pEdge(CLat2); 00384 //printf("%d %d %lf %lf %lf\n",v,vv,PosX,PosY,Plot[v][vv]); 00385 fprintf(FileToWrite,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n", 00386 PosX,PosY,Plot[vx][vy]*.001,Plot[vx][vy],0.,0.); 00387 } 00388 } 00389 } 00390 } 00391 fclose(FileToWrite); 00392 free(Plot); 00393 return 0; 00394 } 00395 //obsolete? 00396 int ElPoly::NChainPSquareF(){ 00397 int NTimes = 20; 00398 double *ChainPArea = (double *)calloc(3*NTimes,sizeof(double)); 00399 double MeanSigma = 0.; 00400 int SubDiv[3] = {NTimes,NTimes,1}; 00401 for(int f=NFile[0];f<NFile[1];f++){ 00402 Processing(f); 00403 if(OpenRisk(cFile[f],BF_PART))return 1; 00404 // OrderPos(); 00405 SubDiv[CLat1] = 1; 00406 SubDiv[CLat2] = 1; 00407 for(int t=0;t<NTimes;t++){ 00408 if( (t%2)==0) 00409 SubDiv[CLat1] += 1; 00410 if( (t%2-1)==0) 00411 SubDiv[CLat2] += 1; 00412 double NEdge = (double)(SubDiv[CLat1]*SubDiv[CLat2]); 00413 double *Plot = (double *)calloc((int)NEdge,sizeof(double)); 00414 double Area = pEdge(CLat1)*pEdge(CLat2)/NEdge; 00415 double SumX = 0.; 00416 double SumX2 = 0.; 00417 for(int p=0;p<pNPart();p++){ 00418 int vx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*SubDiv[CLat1]); 00419 if(vx > SubDiv[CLat1]) continue; 00420 int vy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*SubDiv[CLat2]); 00421 if(vy > SubDiv[CLat2]) continue; 00422 Plot[vx*SubDiv[CLat2] + vy] += 1.; 00423 } 00424 for(int n=0;n<NEdge;n++){ 00425 SumX += Plot[n]/Area; 00426 SumX2 += QUAD((Plot[n]/Area)); 00427 Plot[n] /= Area; 00428 } 00429 int Valori = 20; 00430 double *Intervalli = (double *)calloc(Valori,sizeof(double)); 00431 int IfNorm = 0; 00432 MOMENTI m1 = Mat->Distribuzione(Plot,(int)NEdge,Intervalli,Valori,IfNorm); 00433 ChainPArea[3*t+0] = NEdge; 00434 MeanSigma += m1.Uno; 00435 double Media = SumX/NEdge; 00436 ChainPArea[3*t+1] = sqrt((SumX2 - Media*Media*NEdge)/(NEdge-1)); 00437 printf("%lf %lf-%lf\n",m1.Uno,ChainPArea[3*t+1],m1.Due); 00438 ChainPArea[3*t+2] = Area*QUAD(m1.Due)/QUAD(m1.Uno); 00439 //QUAD(Area)*ChainPArea[4*t+2]/QUAD(ChainPArea[4*t+1]); 00440 free(Plot); 00441 free(Intervalli); 00442 } 00443 } 00444 char *FileName = (char *)calloc(60,sizeof(char)); 00445 sprintf(FileName,"ChainPArea%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho()); 00446 FILE *FileToWrite = fopen(FileName,"w"); 00447 char cSystem[STRSIZE]; 00448 SysDef(cSystem); 00449 fprintf(FileToWrite,"# %s",cSystem); 00450 fprintf(FileToWrite,"#%lf\n",MeanSigma/(double)(NFile[1]-NFile[0])); 00451 fprintf(FileToWrite,"# NEdge SDeviation Area*Dev/NEdge\n"); 00452 for(int t=0;t<NTimes;t++) 00453 fprintf(FileToWrite,"%lf %lf %lf\n",ChainPArea[3*t+0],ChainPArea[3*t+1],ChainPArea[3*t+2]); 00454 fclose(FileToWrite); 00455 free(FileName); 00456 free(ChainPArea); 00457 return 0; 00458 } 00459 int ElPoly::SpectrumF(int NSample){ 00460 int NHalf = NSample/2; 00461 double InvNFile = 1./(double)(NFile[1]-NFile[0]); 00462 double InvNSample = 1./(double)NSample; 00463 double InvNSample2 = 1./(double)SQR(NSample); 00464 double *Plot = (double *) calloc(SQR(NSample),sizeof(double)); 00465 double *PlotUp = (double *) calloc(SQR(NSample),sizeof(double)); 00466 double *PlotDown = (double *) calloc(SQR(NSample),sizeof(double)); 00467 double *CountUp = (double *) calloc(SQR(NSample),sizeof(double)); 00468 double *CountDown = (double *) calloc(SQR(NSample),sizeof(double)); 00469 double *PlotBS = (double *) calloc(SQR(NSample),sizeof(double)); 00470 double *PlotA = (double *)calloc(SQR(NSample), sizeof(double)); 00471 double *PlotAUp = (double *)calloc(SQR(NSample), sizeof(double)); 00472 double *PlotADown = (double *)calloc(SQR(NSample), sizeof(double)); 00473 double *PlotThin = (double *)calloc(SQR(NSample), sizeof(double)); 00474 double *Count = (double *) calloc(NSample*NSample,sizeof(double)); 00475 #ifdef USE_FFTW 00476 fftw_complex *out = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00477 fftw_complex *in = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00478 fftw_plan plan = fftw_plan_dft_2d(NSample,NSample, 00479 in, out,FFTW_FORWARD,FFTW_PATIENT); 00480 fftw_complex *outUp = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00481 fftw_complex *inUp = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00482 fftw_plan planUp = fftw_plan_dft_2d(NSample,NSample, 00483 inUp, outUp,FFTW_FORWARD,FFTW_PATIENT); 00484 fftw_complex *outDown = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00485 fftw_complex *inDown = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00486 fftw_plan planDown = fftw_plan_dft_2d(NSample,NSample, 00487 inDown, outDown,FFTW_FORWARD,FFTW_PATIENT); 00488 //FFTW_MEASURE); 00489 fftw_complex *outThin = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00490 fftw_complex *inThin = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00491 fftw_plan planThin = fftw_plan_dft_2d(NSample,NSample, 00492 inThin, outThin,FFTW_FORWARD,FFTW_PATIENT); 00493 //FFTW_MEASURE); 00494 #endif 00495 for(int f=NFile[0];f<NFile[1];f++){ 00496 Processing(f); 00497 // memset(Plot,0.,NSample*NSample*sizeof(double)); 00498 // memset(Count,0.,NSample*NSample*sizeof(double)); 00499 for(int s=0;s<SQR(NSample);s++){ 00500 Plot[s] = 0.; 00501 Count[s] = 0.; 00502 PlotUp[s] = 0.; 00503 CountUp[s] = 0.; 00504 PlotDown[s] = 0.; 00505 CountDown[s] = 0.; 00506 } 00507 if(OpenRisk(cFile[f],BF_CHAIN))return 0; 00508 for(int p=0;p<pNPart();p++){ 00509 int sx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*NSample); 00510 int sy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*NSample); 00511 if(sx < 0 || sx >= NSample) continue; 00512 if(sy < 0 || sy >= NSample) continue; 00513 Plot[sx*NSample+sy] += pPos(p,CNorm);// - pCm(CNorm); 00514 Count[sx*NSample+sy] += 1.; 00515 //if(Pm[p].Typ != 1) continue; 00516 if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){ 00517 PlotUp[sx*NSample+sy] += pPos(p,CNorm); 00518 CountUp[sx*NSample+sy] += 1.; 00519 } 00520 else { 00521 PlotDown[sx*NSample+sy] += pPos(p,CNorm); 00522 CountDown[sx*NSample+sy] += 1.; 00523 } 00524 } 00525 double MeanThick = 0.; 00526 for(int s=0;s<SQR(NSample);s++){ 00527 Plot[s] /= Count[s] > 0. ? Count[s] : 1.; 00528 PlotUp[s] /= CountUp[s] > 0. ? CountUp[s] : 1.; 00529 PlotDown[s] /= CountDown[s] > 0. ? CountDown[s] : 1.; 00530 MeanThick += CountUp[s] - CountDown[s]; 00531 } 00532 MeanThick /= (double)SQR(NSample); 00533 //InterBSpline2D(Plot,PlotBS,NSample,NSample); 00534 //InterBSpline2D(PlotBS,Plot,NSample,NSample); 00535 #ifdef USE_FFTW 00536 if(1==1){ 00537 for(int s=0;s<SQR(NSample);s++){ 00538 in[s][0] = Plot[s]; 00539 in[s][1] = 0.; 00540 inThin[s][0] = PlotUp[s]-PlotDown[s]; 00541 inThin[s][1] = 0.; 00542 inUp[s][0] = PlotUp[s]; 00543 inUp[s][1] = 0.; 00544 inDown[s][0] = PlotDown[s]; 00545 inDown[s][1] = 0.; 00546 } 00547 fftw_execute(plan); 00548 fftw_execute(planThin); 00549 fftw_execute(planUp); 00550 fftw_execute(planDown); 00551 for(int sx=0;sx<NSample;sx++){ 00552 int ssx = sx < NHalf ? sx + NHalf : sx - NHalf; 00553 for(int sy=0;sy<NSample;sy++){ 00554 int ssy = sy < NHalf ? sy + NHalf : sy - NHalf; 00555 PlotA[ssx*NSample+ssy] += (SQR(out[sx*NSample+sy][0])+SQR(out[sx*NSample+sy][1]))*SQR(InvNSample2); 00556 PlotThin[ssx*NSample+ssy] += (SQR(outThin[sx*NSample+sy][0])+SQR(outThin[sx*NSample+sy][1]))*SQR(InvNSample2); 00557 PlotAUp[ssx*NSample+ssy] += (SQR(outUp[sx*NSample+sy][0])+SQR(outUp[sx*NSample+sy][1]))*SQR(InvNSample2); 00558 PlotADown[ssx*NSample+ssy] += (SQR(outDown[sx*NSample+sy][0])+SQR(outDown[sx*NSample+sy][1]))*SQR(InvNSample2); 00559 } 00560 } 00561 } 00562 else { 00563 int NMax = NSample; 00564 double *st = Plot; 00565 double *sw = PlotA; 00566 double dNMax = 1./(double)NMax; 00567 int NHalf = (int)(NMax/2.); 00568 for(int kx=-NHalf;kx<NHalf;kx++){ 00569 double qx = kx*dNMax; 00570 for(int ky=-NHalf;ky<NHalf;ky++){ 00571 double qy = ky*dNMax; 00572 double Re2=0.,Im2=0.; 00573 double Re1=0.,Im1=0.; 00574 for(int lx=0;lx<NMax;lx++){ 00575 for(int ly=0;ly<NMax;ly++){ 00576 double Arg = 2.*M_PI*(lx*qx + ly*qy); 00577 double cy = cos(Arg); 00578 double sy = sin(Arg); 00579 Re1 += st[lx*NMax+ly]*cy; 00580 Im1 += st[lx*NMax+ly]*sy; 00581 } 00582 } 00583 int kkx = kx + NHalf; 00584 int kky = ky + NHalf; 00585 sw[kkx*NMax+kky] += SQR(Re1*dNMax) + SQR(Im1*dNMax); 00586 } 00587 } 00588 } 00589 if(1==0){ 00590 FILE *FMidplane = fopen("Midplane.dat","w"); 00591 FILE *FThickness = fopen("Thickness.dat","w"); 00592 for(int sx=0;sx<NSample;sx++){ 00593 double x = sx*pEdge(CLat1)*InvNSample; 00594 for(int sy=0;sy<NSample;sy++){ 00595 double y = sy*pEdge(CLat2)*InvNSample; 00596 fprintf(FMidplane,"%lf %lf %lf \n",x,y, Plot[sx*NSample+sy]); 00597 fprintf(FThickness,"%lf %lf %lf \n",x,y,PlotUp[sx*NSample+sy]-PlotDown[sx*NSample+sy]); 00598 } 00599 } 00600 fclose(FMidplane); 00601 fclose(FThickness); 00602 } 00603 #else 00604 Mat->Spettro2d(Plot,PlotA,NSample); 00605 #endif //FFTW3_H 00606 } 00607 printf("\n"); 00608 #ifdef OMPI_MPI_H 00609 MPI_Allreduce(MPI_IN_PLACE,PlotA,NSample*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid); 00610 int Rank=0; 00611 MPI_Comm_rank(Proc->CommGrid, &Rank); 00612 if(Rank==0){ 00613 #endif 00614 char *NomeFile = (char*)calloc(60,sizeof(char)); 00615 sprintf(NomeFile,"Spectrum.xvl",pchiN(),pkappaN(),prho()); 00616 FILE *FileToWrite = fopen(NomeFile,"w"); 00617 fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[squaremesh]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NSample); 00618 for(int sx=0;sx<NSample;sx++){ 00619 for(int sy=0;sy<NSample;sy++){ 00620 double x = pEdge(CLat1)*sx*InvNSample; 00621 double y = pEdge(CLat2)*sy*InvNSample; 00622 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.); 00623 } 00624 } 00625 printf("Write\n"); 00626 fclose(FileToWrite); 00627 sprintf(NomeFile,"Spectrum.dat",pchiN(),pkappaN(),prho()); 00628 FileToWrite = fopen(NomeFile,"w"); 00629 fprintf(FileToWrite,"# 1) q^2 2) h(q)^2*L_xL_y | q_x=s2pi/(L_x) NSample %d\n",NSample); 00630 for(int kx=0;kx<NSample;kx++){ 00631 int kx1 = kx - NSample/2; 00632 double qx = 2.*M_PI*kx*pInvEdge(CLat1); 00633 double qxt = 2.*M_PI*kx1*pInvEdge(CLat1); 00634 for(int ky=0;ky<NSample;ky++){ 00635 int ky1 = ky - NSample/2; 00636 double qy = 2.*M_PI*ky*pInvEdge(CLat2); 00637 double qyt = 2.*M_PI*ky1*pInvEdge(CLat2); 00638 double qq = SQR(qxt) + SQR(qyt); 00639 double Spe = PlotA[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2); 00640 double SpeUp = PlotAUp[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2); 00641 double SpeDown = PlotADown[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2); 00642 double SpeThin = PlotThin[kx*NSample+ky]*InvNFile*pEdge(CLat1)*pEdge(CLat2); 00643 fprintf(FileToWrite,"%lf %g %g %g %g\n",qq,Spe,SpeUp,SpeDown,SpeThin); 00644 } 00645 } 00646 fclose(FileToWrite); 00647 free(NomeFile); 00648 #ifdef OMPI_MPI_H 00649 } 00650 #endif 00651 free(Plot); 00652 free(PlotA); 00653 free(PlotBS); 00654 free(Count); 00655 #ifdef USE_FFTW 00656 fftw_destroy_plan(plan); 00657 #endif 00658 return 0; 00659 } 00660 void ElPoly::Midplane(int NSample){ 00661 int NHalf = NSample/2; 00662 double InvNFile = 1./(double)(NFile[1]-NFile[0]); 00663 double InvNSample = 1./(double)NSample; 00664 double *PlotBS = (double *) calloc(SQR(NSample),sizeof(double)); 00665 double *Count = (double *) calloc(NSample*NSample,sizeof(double)); 00666 double *Plot = (double *) calloc(SQR(NSample),sizeof(double)); 00667 double *PlotA = (double *) calloc(SQR(NSample),sizeof(double)); 00668 for(int f=NFile[0];f<NFile[1];f++){ 00669 Processing(f); 00670 // memset(Plot,0.,NSample*NSample*sizeof(double)); 00671 // memset(Count,0.,NSample*NSample*sizeof(double)); 00672 for(int sx=0;sx<NSample;sx++){ 00673 for(int sy=0;sy<NSample;sy++){ 00674 Plot[sx*NSample+sy] = 0.; 00675 } 00676 } 00677 if(OpenRisk(cFile[f],NBackFold))return; 00678 for(int p=0;p<pNPart();p++){ 00679 int sx = (int)(pPos(p,CLat1)*pInvEdge(CLat1)*NSample); 00680 int sy = (int)(pPos(p,CLat2)*pInvEdge(CLat2)*NSample); 00681 if(sx < 0 || sx >= NSample) continue; 00682 if(sy < 0 || sy >= NSample) continue; 00683 Plot[sx*NSample+sy] += pPos(p,CNorm);// - pCm(CNorm); 00684 Count[sx*NSample+sy] += 1.; 00685 } 00686 for(int sx=0;sx<NSample;sx++){ 00687 for(int sy=0;sy<NSample;sy++){ 00688 Plot[sx*NSample+sy] /= Count[sx*NSample+sy] > 0. ? Count[sx*NSample+sy] : 1.; 00689 PlotA[sx*NSample+sy] += Plot[sx*NSample+sy]; 00690 } 00691 } 00692 // InterBSpline2D(Plot,PlotBS,NSample,NSample); 00693 // InterBSpline2D(PlotBS,Plot,NSample,NSample); 00694 char FName[60]; 00695 sprintf(FName,"Midplane%05d.dat",f); 00696 FILE *FMid = fopen(FName,"w"); 00697 for(int sy=0;sy<NSample;sy++){ 00698 for(int sx=0;sx<NSample;sx++){ 00699 double x = pEdge(CLat1)*sx*InvNSample; 00700 double y = pEdge(CLat2)*sy*InvNSample; 00701 fprintf(FMid,"%lf %lf %lf\n",x,y,Plot[sx*NSample+sy]); 00702 } 00703 } 00704 fclose(FMid); 00705 } 00706 printf("\n"); 00707 FILE *FMid = fopen("average_middplane.dat","w"); 00708 for(int sy=0;sy<NSample;sy++){ 00709 for(int sx=0;sx<NSample;sx++){ 00710 double x = pEdge(CLat1)*sx*InvNSample; 00711 double y = pEdge(CLat2)*sy*InvNSample; 00712 fprintf(FMid,"%lf %lf %lf\n",x,y,PlotA[sx*NSample+sy]*InvNFile); 00713 } 00714 } 00715 fclose(FMid); 00716 free(Plot); 00717 free(Count); 00718 free(PlotBS); 00719 } 00720 void ElPoly::SpectrumMidplane(int NSample){ 00721 // FILE *Ciccia = fopen("Ciccia.dat","w"); 00722 // for(int p=0;p<pNPart();p++){ 00723 // fprintf(Ciccia,"%lf %lf\n",SQR(pPos(p,0)-.5*pEdge(0))+SQR(pPos(p,1)-.5*pEdge(1)),pPos(p,2)); 00724 // } 00725 // fclose(Ciccia); 00726 // return; 00727 int NHalf = NSample/2; 00728 double InvNSample = 1./(double)NSample; 00729 double *Plot = (double *) calloc(NSample*NSample,sizeof(double)); 00730 double *Count = (double *) calloc(NSample*NSample,sizeof(double)); 00731 double *PlotA = (double *) calloc(NSample*NSample,sizeof(double)); 00732 int NPoint = (int)(NSample*sqrt(2.)); 00733 double *Density = (double *) calloc(NPoint,sizeof(double)); 00734 #ifdef USE_FFTW 00735 fftw_complex *out = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00736 fftw_complex *in = (fftw_complex *)fftw_malloc(SQR(NSample)*sizeof(fftw_complex)); 00737 fftw_plan plan1 = fftw_plan_dft_2d(NSample,NSample, 00738 in, out,FFTW_FORWARD,FFTW_PATIENT); 00739 for(int f=NFile[0];f<NFile[1];f++){ 00740 Processing(f); 00741 for(int sx=0;sx<NSample;sx++){ 00742 for(int sy=0;sy<NSample;sy++){ 00743 Plot[sx*NSample+sy] = 0.; 00744 Count[sx*NSample+sy] = 0.; 00745 } 00746 } 00747 FILE *Ciccia = fopen(cFile[f],"r"); 00748 char cLine[256]; 00749 for(int p=0;!(fgets(cLine,STRSIZE,Ciccia)==NULL);p++){ 00750 double x = 0., y = 0., z = 0.; 00751 sscanf(cLine,"%lf %lf %lf\n",&x,&y,&z); 00752 int v = (int)(x*NSample*pInvEdge(CLat1));//(pEdge(0)+.4)); 00753 int vv = (int)(y*NSample*pInvEdge(CLat2));//(pEdge(1)+.4)); 00754 if(v < 0 || v >= NSample) continue; 00755 if(vv < 0 || vv >= NSample) continue; 00756 Plot[v*NSample+vv] += z; 00757 Count[v*NSample+vv] += 1.; 00758 } 00759 for(int s=0;s<NSample;s++){ 00760 for(int ss=0;ss<NSample;ss++){ 00761 Plot[s*NSample+ss] /= Count[s*NSample+ss] > 0. ? Count[s*NSample+ss] : 1.; 00762 } 00763 } 00764 fclose(Ciccia); 00765 if(1==0){ 00766 for(int sx=0;sx<NSample;sx++){ 00767 for(int sy=0;sy<NSample;sy++){ 00768 in[sx*NSample+sy][0] = Plot[sx*NSample+sy]; 00769 in[sx*NSample+sy][1] = 0.; 00770 } 00771 } 00772 fftw_execute(plan1); 00773 for(int sx=0;sx<NSample;sx++){ 00774 int sx1 = sx + NHalf; 00775 if(sx1 >= NSample) sx1 -= NSample; 00776 for(int sy=0;sy<NSample;sy++){ 00777 int sy1 = sy + NHalf; 00778 if(sy1 >= NSample) sy1 -= NSample; 00779 PlotA[sx*NSample+sy] += (SQR(out[sx*NSample+sy][0])+SQR(out[sx*NSample+sy][1]))*SQR(InvNSample); 00780 //PlotA[sx*NSample+sy] += Plot[sx][sy]; 00781 } 00782 } 00783 } 00784 else { 00785 int NMax = NSample; 00786 double *st = Plot; 00787 double *sw = PlotA; 00788 double dNMax = 1./(double)NMax; 00789 int NHalf = (int)(NMax/2.); 00790 for(int kx=-NHalf;kx<NHalf;kx++){ 00791 double qx = kx*dNMax; 00792 for(int ky=-NHalf;ky<NHalf;ky++){ 00793 double qy = ky*dNMax; 00794 double Re2=0.,Im2=0.; 00795 double Re1=0.,Im1=0.; 00796 for(int lx=0;lx<NMax;lx++){ 00797 for(int ly=0;ly<NMax;ly++){ 00798 double Arg = 2.*M_PI*(lx*qx + ly*qy); 00799 double cy = cos(Arg); 00800 double sy = sin(Arg); 00801 Re1 += st[lx*NMax+ly]*cy; 00802 Im1 += st[lx*NMax+ly]*sy; 00803 } 00804 } 00805 int kkx = kx + NMax/2; 00806 int kky = ky + NMax/2; 00807 sw[kkx*NMax+kky] += SQR(Re1*dNMax) + SQR(Im1*dNMax); 00808 } 00809 } 00810 } 00811 } 00812 #endif // USE_FFTW 00813 printf("\n"); 00814 #ifdef OMPI_MPI_H 00815 MPI_Allreduce(MPI_IN_PLACE,PlotA,NSample*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid); 00816 int Rank=0; 00817 MPI_Comm_rank(Proc->CommGrid, &Rank); 00818 if(Rank==0){ 00819 #endif 00820 char *NomeFile = (char*)calloc(60,sizeof(char)); 00821 sprintf(NomeFile,"SpectrumChi%.0fKappa%.0fRho%.0f.xvl",pchiN(),pkappaN(),prho()); 00822 FILE *FileToWrite = fopen(NomeFile,"w"); 00823 // fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[squaremesh]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NSample); 00824 double DivInv = 1./(double)(NFile[1]-NFile[0]); 00825 for(int s=0;s<NSample;s++){ 00826 for(int ss=0;ss<NSample;ss++){ 00827 fprintf(FileToWrite,"%lf %lf %lf\n",pEdge(CLat1)*s*InvNSample,pEdge(CLat2)*ss*InvNSample,PlotA[s*NSample+ss]*DivInv); 00828 // 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)); 00829 } 00830 } 00831 fclose(FileToWrite); 00832 sprintf(NomeFile,"SpectrumChi%.0fKappa%.0fRho%.0f.dat",pchiN(),pkappaN(),prho()); 00833 FileToWrite = fopen(NomeFile,"w"); 00834 double qHalf = QUAD(pEdge(CLat1)*(NHalf)*InvNSample)+QUAD(pEdge(CLat2)*(NHalf)*InvNSample); 00835 double dx = pEdge(CLat1)*InvNSample; 00836 double dy = pEdge(CLat2)*InvNSample; 00837 for(int kx=0;kx<NSample;kx++){ 00838 int kx1 = kx - NSample/2; 00839 double qx = kx*InvNSample; 00840 double qxt = kx1*InvNSample; 00841 for(int ky=0;ky<NSample;ky++){ 00842 int ky1 = ky - NSample/2; 00843 double qy = ky*InvNSample; 00844 double qyt = ky1*InvNSample; 00845 double qq = SQR(qx/dx - .5/dx) + SQR(qy/dy - .5/dy); 00846 //double qq = SQR(qx) + SQR(qy); 00847 fprintf(FileToWrite,"%lf %g\n",qq,PlotA[kx*NSample+ky]*DivInv); 00848 //fprintf(FileToWrite,"%lf %g\n",QUAD(pEdge(CLat1)*(s-NHalf)*InvSample)+QUAD(pEdge(CLat2)*(ss-NHalf)*InvSample),PlotA[s*NSample+ss]*DivInv); 00849 } 00850 } 00851 // for(int p=1;p<NPoint;p++) 00852 // fprintf(FileToWrite,"%d %lf\n",p,Density[p]*DivInv); 00853 fclose(FileToWrite); 00854 free(NomeFile); 00855 #ifdef OMPI_MPI_H 00856 } 00857 #endif 00858 free(Density); 00859 free(PlotA); 00860 free(Plot); 00861 } 00862 void ElPoly::HeaderAverage(int nNano){ 00863 if(nNano >= pNNano()){ 00864 printf("The specified nanoparticle doesn't exist\n"); 00865 nNano = 0; 00866 // return; 00867 } 00868 int NDim = 2; 00869 int NumFile = NFile[1]-NFile[0]; 00870 double *Area = (double *)calloc(NumFile,sizeof(double)); 00871 double *NanoPos = (double *)calloc(6*NumFile,sizeof(double)); 00872 double *NanoDiff = (double *)calloc(2*NumFile,sizeof(double)); 00873 double *NanoCount = (double *)calloc(NumFile,sizeof(double)); 00874 double *NanoDist = (double *)calloc(2*NumFile,sizeof(double)); 00875 double *Time = (double *)calloc(NumFile,sizeof(double)); 00876 double NanoInit[3] = {Nano->Pos[0],Nano->Pos[1],Nano->Pos[2]}; 00877 double AreaAv = 0.; 00878 double AreaErr = 0.; 00879 double InitTime = pTime(); 00880 //indentation 00881 double InitPos = pPos(0,CNorm) - pPos(1,CNorm) - MAX(Nano[0].Rad,Nano[0].Height); 00882 for(int f=NFile[0];f<NFile[1];f++){ 00883 Processing(f); 00884 FILE *File2Read; 00885 if((File2Read = fopen(cFile[f],"r"))==0){ 00886 printf("The file is missing\n"); 00887 return ; 00888 } 00889 ReadHeader(File2Read); 00890 fclose(File2Read); 00891 AreaAv += pEdge(CLat1)*pEdge(CLat2); 00892 AreaErr += SQR(pEdge(CLat1)*pEdge(CLat2)); 00893 Area[f] = pEdge(CLat1)*pEdge(CLat2); 00894 Time[f] = pTime() - InitTime; 00895 for(int d=0;d<3;d++) NanoPos[f*6+d] = Nano[nNano].Pos[d]; 00896 NanoPos[f*6+3] = Nano[nNano].Rad; 00897 NanoPos[f*6+4] = Nano[nNano].Height; 00898 NanoPos[f*6+5] = Nano[nNano].Area; 00899 for(int ff=f;ff>=0;ff--){ 00900 int fff = f-ff; 00901 NanoDiff[fff*2] += SQR(NanoPos[f*6+0]-NanoPos[ff*6+0]) + SQR(NanoPos[f*6+1]-NanoPos[ff*6+1]); 00902 NanoCount[fff] += 1.; 00903 } 00904 NanoDiff[f*2+1] = SQR(NanoPos[f*6+0]-NanoInit[0]) + SQR(NanoPos[f*6+1]-NanoInit[1]); 00905 if(pNNano() > 1){ 00906 //indentation 00907 NanoDist[f*2] = InitPos - pNanoPos(0,CNorm) + pNanoPos(1,CNorm) + MAX(Nano[nNano].Rad,Nano[nNano].Height); 00908 double Dist = 0.; 00909 for(int d=0;d<3;d++){ 00910 double Pos1 = pNanoPos(0,d) - floor(pNanoPos(0,d)*pInvEdge(d))*pEdge(d); 00911 Dist += SQR(pNanoPos(0,d)-pNanoPos(1,d)); 00912 } 00913 Dist = sqrt(Dist); 00914 //two np dist 00915 NanoDist[f*2+1] = Dist; 00916 } 00917 } 00918 AreaAv /= (double)(NFile[1] - NFile[0]); 00919 AreaErr = sqrt( AreaErr - (NFile[1] - NFile[0])*SQR(AreaAv))/(double)(NFile[1] - NFile[0]); 00920 printf("\n"); 00921 printf("Area %lf \\pm %lf\n",AreaAv,AreaErr); 00922 FILE *File2Write = fopen("AreaTime.dat","w"); 00923 fprintf(File2Write,"# Time Area %lf pm %lf\n",AreaAv,AreaErr); 00924 char FName[120]; 00925 sprintf(FName,"NanoDiffR%.1fS%.1fH%.1f.dat",Nano[nNano].Rad,Nano[nNano].Hamaker,Nano[nNano].Height); 00926 FILE *FNano = fopen(FName,"w"); 00927 fprintf(FNano,"# t <msd> msd Dz NpDist rad hei\n"); 00928 for(int f=0;f<NumFile;f++){ 00929 NanoDiff[f*2] /= NanoCount[f] > 0. ? NanoCount[f] : 1.; 00930 fprintf(File2Write,"%lf %lf\n",Time[f],Area[f]); 00931 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]); 00932 } 00933 fclose(File2Write); 00934 fclose(FNano); 00935 free(NanoDiff); 00936 free(NanoDist); 00937 free(Area); 00938 free(Time); 00939 } 00940 //only one block 00941 void ElPoly::WidomOut(char *NrgFile,int NBin){ 00942 int NChains = 1; 00943 Block[0].NChain -= NChains; 00944 WriteXvt("SnapOut.dat"); 00945 FILE *FileLast = fopen("LastChain.dat","w"); 00946 for(int p=pNPart()-(NChains)*pNPCh();p<pNPart();p++){ 00947 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)); 00948 } 00949 fclose(FileLast); 00950 } 00951 //only one block 00952 void ElPoly::WidomOut(){ 00953 int NChains = 1; 00954 Block[0].NChain -= NChains; 00955 char FileName[60]; 00956 sprintf(FileName,"SnapOut%05d.dat",pNChain()-1); 00957 WriteXvt(FileName); 00958 for(int c=0;c<pNChain()-1;c++){ 00959 SwapChain(c,pNChain()-1); 00960 sprintf(FileName,"SnapOut%05d.dat",c); 00961 WriteXvt(FileName); 00962 } 00963 } 00964 //the system has an additional fake chain 00965 void ElPoly::WidomIn(){ 00966 int NFile = 1000; 00967 char FileName[60]; 00968 SetNPart(pNPart()+pNPCh()); 00969 SetNChain(pNChain()+1); 00970 Block[0].NChain += 1; 00971 int p1 = pNPCh()*(pNChain()-1); 00972 double Var = sqrt(SQR(pReOverCutOff())/(double)(pNPCh()-1)/3.)/2.; 00973 for(int f=0;f<NFile;f++){ 00974 for(int d=0;d<3;d++){ 00975 SetPos(p1,d,Mat->Casuale()*pEdge(d)); 00976 } 00977 for(int p=p1+1;p<pNPart();p++){ 00978 for(int d=0;d<3;d++){ 00979 SetPos(p1,d,pPos(p-1,d) + Mat->Gaussiano(0.,Var)); 00980 } 00981 if( p%pNPCh() >= Block[0].Asym ) 00982 SetType(p,1); 00983 } 00984 sprintf(FileName,"SnapIn%05d.dat",f); 00985 Write(FileName); 00986 } 00987 } 00988 //the system has an additional fake chain 00989 void ElPoly::WidomIn(char *NrgFile,int NBin){ 00990 int NChains = 1; 00991 Block[0].NChain += NChains; 00992 WriteXvt("SnapIn.dat"); 00993 int NAddPart = NChains*pNPCh(); 00994 FILE *FileLast = fopen("SnapIn.dat","a"); 00995 PART *Pn = (PART *)calloc(NAddPart,sizeof(PART)); 00996 double GaussVar = sqrt(SQR(pReOverCutOff())/(double)(pNPCh()-1)/3.)*2./3.; 00997 for(int c=0;c<NChains;c++){ 00998 int p1 = c*pNPCh(); 00999 for(int d=0;d<3;d++){ 01000 SetPos(p1,d,Mat->Casuale()*pEdge(d)); 01001 } 01002 for(int p=1;p<pNPCh();p++){ 01003 for(int d=0;d<3;d++){ 01004 Pn[p+p1].Pos[d] = Pn[p-1+p1].Pos[d] + Mat->Gaussiano(0,GaussVar); 01005 } 01006 Pn[p+p1].Typ = 0; 01007 if(p >= Block[0].Asym){ 01008 Pn[p+p1].Typ = 1; 01009 } 01010 } 01011 } 01012 for(int p=0;p<NAddPart;p++){ 01013 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); 01014 } 01015 free(Pn); 01016 fclose(FileLast); 01017 } 01018 void ElPoly::End2EndDistr(char *OutFile){ 01019 char FileName[256]; 01020 for(int f=NFile[0];f<NFile[1];f++){ 01021 sprintf(FileName,"%s%04d.dat",OutFile,f); 01022 FILE *File2Write = fopen(FileName,"w"); 01023 Processing(f); 01024 if(OpenRisk(cFile[f],BF_NO))return; 01025 for(int c=0;c<pNChain();c++){ 01026 int p1 = c*pNPCh(); 01027 int p2 = c*pNPCh() + pNPCh() - 1; 01028 double End2End = 0.; 01029 for(int d=0;d<3;d++){ 01030 End2End += SQR(pPos(p1,d) - pPos(p2,d)); 01031 } 01032 fprintf(File2Write,"%lf ",End2End); 01033 fprintf(File2Write,"\n"); 01034 } 01035 fclose(File2Write); 01036 } 01037 } 01038 void ElPoly::Decoupling(int What){ 01039 NFileTot = NFile[1] - NFile[0]; 01040 double *Angles = (double *)calloc((NFile[1]-NFile[0])*pNChain(),sizeof(double)); 01041 double *Count = (double *)calloc((NFile[1]-NFile[0]),sizeof(double)); 01042 double *AnglesNano = (double *)calloc((NFile[1]-NFile[0])*pNNano(),sizeof(double)); 01043 double *Time = (double *)calloc((NFile[1]-NFile[0]),sizeof(double)); 01044 double *AngleDiff = (double *)calloc(pNChain(),sizeof(double)); 01045 double *AngleDiff2 = (double *)calloc(pNChain(),sizeof(double)); 01046 Vettore Ax0(1.,0.,0.); 01047 Vettore Ax2(0.,0.,1.); 01048 double InvPhob = Block[0].Asym > 0 ? 1./(double)Block[0].Asym : 1.; 01049 double InvPhil = Block[0].Asym > 0 ? 1./(double)(pNPCh() - Block[0].Asym) : 1.; 01050 double DistRel[4]; 01051 Vettore ChDir(0.,0.,0.); 01052 Vettore NanoAx(0.,0.,0.); 01053 double Average = 0.; 01054 double Variance = 0.; 01055 double InitTime = pTime(); 01056 for(int f=NFile[0];f<NFile[1];f++){ 01057 Processing(f); 01058 int ff = f - NFile[0]; 01059 if(OpenRisk(cFile[f],BF_CHAIN))return ; 01060 Time[f] = pTime() - InitTime; 01061 for(int n=0;n<pNNano();n++){ 01062 for(int d=0;d<3;d++){ 01063 NanoAx.Set(Nano[n].Axis[d],d); 01064 } 01065 NanoAx.Set(0.,CNorm); 01066 AnglesNano[ff*pNNano()+n] = Ax0.Angle(&NanoAx); 01067 } 01068 for(int c=0,b=0;c<pNChain();c++){ 01069 if(What == 0){ 01070 ChDir.Set(Ch[c].Dir[0],0); 01071 ChDir.Set(Ch[c].Dir[1],1); 01072 double Angle = Ax0.Angle(&ChDir); 01073 if(isnan(Angle)) Angles[ff*pNChain()+c] = Angles[(ff-1)*pNChain()+c]; 01074 else Angles[ff*pNChain()+c] = Ax0.Angle(&ChDir); 01075 } 01076 else if(What == 1){ 01077 Angles[ff*pNChain()+c] = TwoPartDist(c*pNPCh(),(c+1)*pNPCh()-1,DistRel); 01078 } 01079 Average += Angles[ff*pNChain()+c]; 01080 Variance += SQR(Angles[ff*pNChain()+c]); 01081 } 01082 } 01083 printf("\n"); 01084 double NTot = (double)(pNChain()*(NFile[1]-NFile[0])); 01085 Average /= NTot > 0. ? NTot : 1.; 01086 double PreVar = Variance; 01087 Variance = (Variance - SQR(Average)*NTot)/(NTot-1); 01088 if(isnan(Variance)){ 01089 printf("Variance is not a number %lf %lf %lf\n",PreVar,Average,NTot); 01090 Variance = 1.; 01091 } 01092 for(int f=0;f<NFileTot;f++){ 01093 for(int c=0;c<pNChain();c++){ 01094 AngleDiff2[f] += (Angles[f*pNChain()+c]-Average)*(Angles[0*pNChain()+c]-Average); 01095 } 01096 for(int ff = f+1; ff<NFileTot;ff++){ 01097 int fff = ff - f; 01098 for(int c=0;c<pNChain();c++){ 01099 AngleDiff[fff] += (Angles[ff*pNChain()+c]-Average)*(Angles[f*pNChain()+c]-Average); 01100 Count[fff] += 1.; 01101 } 01102 } 01103 } 01104 char FileName[256]; 01105 if(What == 0) 01106 sprintf(FileName,"DirDecoupling.dat"); 01107 else if(What == 1) 01108 sprintf(FileName,"E2EDecoupling.dat"); 01109 FILE *FileToWrite = fopen(FileName,"w"); 01110 char cSystem[STRSIZE]; 01111 SysDef(cSystem); 01112 // fprintf(FileToWrite,"%s",cSystem); 01113 fprintf(FileToWrite,"#Time AngleProd\n"); 01114 for(int f=1;f<NFileTot;f++){ 01115 double Inv = Count[f] > 0. ? Count[f] : 1.; 01116 fprintf(FileToWrite,"%lf %lf %lf\n",Time[f],AngleDiff[f]/(Inv*Variance),AngleDiff2[f]/(Variance*pNChain()) ); 01117 } 01118 fclose(FileToWrite); 01119 for(int n=0;n<pNNano();n++){ 01120 sprintf(FileName,"NanoAngles%d.dat",n); 01121 FILE *FNanoAng = fopen(FileName,"w"); 01122 for(int f=0;f<NFileTot;f++) 01123 fprintf(FNanoAng,"%lf %lf \n",Time[f],AnglesNano[f*pNNano()+n]*360./DUE_PI); 01124 fclose(FNanoAng); 01125 } 01126 // FileToWrite = fopen("RawAngles.dat","w"); 01127 // for(int f=0;f<NFileTot;f++){ 01128 // for(int c=0;c<pNChain();c++){ 01129 // fprintf(FileToWrite,"%lf\n",Angles[f*pNChain()+c]); 01130 // } 01131 // } 01132 // fclose(FileToWrite); 01133 free(Count); 01134 free(Time); 01135 free(Angles); 01136 free(AngleDiff); 01137 free(AngleDiff2); 01138 } 01146 void ElPoly::AreaCompr(int NSample){ 01147 int SubDiv[3] = {NSample,NSample,1}; 01148 int NEdge = (SubDiv[CLat1]*SubDiv[CLat2]); 01149 //mean StdDev Area*StdDev/Mean^2 01150 double *ChainPArea = (double *)calloc(3*NSample,sizeof(double)); 01151 double *Plot = (double *)calloc(NEdge,sizeof(double)); 01152 double *Area = (double *)calloc(NSample,sizeof(double)); 01153 double *NDiv = (double *)calloc(NSample,sizeof(double)); 01154 char FileName[256]; 01155 //for the simulation in NPtT a smaller box is used 01156 double Edge[3] = {pEdge(0)-1.,pEdge(1)-1.,pEdge(2)}; 01157 double InvEdge[3] = {1./Edge[0],1./Edge[1],1./Edge[2]}; 01158 for(int f=NFile[0];f<NFile[1];f++){ 01159 Processing(f); 01160 OpenRisk(cFile[f],BF_CHAIN); 01161 for(int t=0;t<NSample;t++){ 01162 SubDiv[CLat1] = NSample-t; 01163 SubDiv[CLat2] = NSample-t; 01164 if(t == NSample -1) SubDiv[CLat2] = 2; 01165 NDiv[t] = SubDiv[CLat1]*SubDiv[CLat2]; 01166 Area[t] = Edge[CLat1]*Edge[CLat2]/(SubDiv[CLat1]*SubDiv[CLat2]); 01167 for(int n=0;n<NEdge;n++) Plot[n] = 0.; 01168 for(int c=0;c<pNChain();c++){ 01169 double Posx = pChPos(c,CLat1) - floor(pChPos(c,CLat1)*pInvEdge(CLat1))*pEdge(CLat1); 01170 double Posy = pChPos(c,CLat2) - floor(pChPos(c,CLat2)*pInvEdge(CLat2))*pEdge(CLat2); 01171 int vx = (int)(Posx*InvEdge[CLat1]*SubDiv[CLat1]); 01172 int vy = (int)(Posy*InvEdge[CLat2]*SubDiv[CLat2]); 01173 if(vx < 0 || vx >= SubDiv[CLat1]) continue; 01174 if(vy < 0 || vy >= SubDiv[CLat2]) continue; 01175 Plot[vx*NSample+vy] += 1.; 01176 } 01177 double SumX = 0.; 01178 double SumX2 = 0.; 01179 double Count = 0.; 01180 for(int sx = 0;sx<SubDiv[CLat1];sx++){ 01181 for(int sy = 0;sy<SubDiv[CLat2];sy++){ 01182 int n = sx*NSample+sy; 01183 Plot[n] /= Area[t]; 01184 SumX += Plot[n]; 01185 Count += 1.; 01186 } 01187 } 01188 double NInv = Count > 0. ? 1./Count : 1.; 01189 double MeanVal = SumX*NInv; 01190 for(int sx = 0;sx < SubDiv[CLat1];sx++){ 01191 for(int sy = 0;sy < SubDiv[CLat2];sy++){ 01192 int n = sx*NSample+sy; 01193 SumX2 += SQR(Plot[n] - MeanVal); 01194 } 01195 } 01196 ChainPArea[3*t ] += MeanVal; 01197 ChainPArea[3*t+1] += SumX2/(Count - 1.); 01198 ChainPArea[3*t+2] += 1.; 01199 if(1==0){//writes the single distribution files 01200 sprintf(FileName,"Distr%.0fKappa%02d.dat",pkappaN(),t); 01201 FILE *FileToWrite = fopen(FileName,"a"); 01202 for(int l1=0;l1<SubDiv[CLat1];l1++){ 01203 for(int l2=0;l2<SubDiv[CLat2];l2++){ 01204 fprintf(FileToWrite,"%lf\n",Plot[l1*NSample+l2]); 01205 } 01206 } 01207 fclose(FileToWrite); 01208 } 01209 } 01210 } 01211 #ifdef OMPI_MPI_H 01212 MPI_Allreduce(MPI_IN_PLACE,ChainPArea,3*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid); 01213 int Rank=0; 01214 MPI_Comm_rank(Proc->CommGrid, &Rank); 01215 if(Rank==0){ 01216 #endif 01217 printf("\n"); 01218 free(Plot); 01219 //--------normalizing-------------------- 01220 double Norm = 1./(double)(NFile[1]-NFile[0]); 01221 for(int t=0;t<NSample;t++){ 01222 ChainPArea[3*t ] *= Norm; 01223 ChainPArea[3*t+1] *= Norm; 01224 ChainPArea[3*t+2] *= Norm; 01225 } 01226 sprintf(FileName,"ChainPArea%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN()); 01227 FILE *FileToWrite = fopen(FileName,"w"); 01228 char cSystem[STRSIZE]; 01229 SysDef(cSystem); 01230 fprintf(FileToWrite,"# %s",cSystem); 01231 fprintf(FileToWrite,"# Area MeanChPArea SDeviation Area*Dev/MeanChPArea\n"); 01232 for(int t=0;t<NSample;t++){ 01233 double Compr = Area[t]*ChainPArea[3*t+1]/SQR(ChainPArea[3*t ]); 01234 fprintf(FileToWrite,"%lf %lf %lf %lf\n",Area[t],(ChainPArea[3*t ]),ChainPArea[3*t+1],Compr); 01235 //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]); 01236 } 01237 fclose(FileToWrite); 01238 #ifdef OMPI_MPI_H 01239 } 01240 #endif 01241 free(ChainPArea); 01242 } 01248 void ElPoly::ElasticCoupling(int NSample){ 01249 double Edge[3] = {pEdge(0)-1.,pEdge(1)-1.,pEdge(2)}; 01250 int SubDiv[3] = {NSample,NSample,1}; 01251 //mean StdDev Area*StdDev/Mean^2 01252 double *ElCoup = (double *)calloc(3*(NSample+1),sizeof(double)); 01253 double *Area = (double *)calloc((NSample+1),sizeof(double)); 01254 double *NDiv = (double *)calloc((NSample+1),sizeof(double)); 01255 double MeanSigma = 0.; 01256 char FileName[256]; 01257 //for the simulation in NPtT a smaller box is used 01258 double InvEdge[3] = {1./Edge[0],1./Edge[1],1./Edge[2]}; 01259 double *PosUp = (double *)calloc(SQR(NSample),sizeof(double)); 01260 double *PosDown = (double *)calloc(SQR(NSample),sizeof(double)); 01261 double *CountUp = (double *)calloc(SQR(NSample),sizeof(double)); 01262 double *CountDown= (double *)calloc(SQR(NSample),sizeof(double)); 01263 for(int f=NFile[0];f<NFile[1];f++){ 01264 Processing(f); 01265 OpenRisk(cFile[f],BF_CHAIN); 01266 for(int t=0;t<NSample+1;t++){ 01267 //defining the partitions 01268 // SubDiv[CLat1] = (int)pow(2.,t); 01269 // SubDiv[CLat2] = (int)pow(2.,t); 01270 SubDiv[CLat1] = NSample - t; 01271 SubDiv[CLat2] = NSample - t; 01272 if(t == NSample-1) SubDiv[CLat2] = 2; 01273 if(t == NSample){ 01274 SubDiv[CLat1] = 1; 01275 SubDiv[CLat2] = 1; 01276 } 01277 NDiv[t] = SubDiv[CLat1]*SubDiv[CLat2]; 01278 Area[t] = Edge[CLat1]*Edge[CLat2]/(double)(SubDiv[CLat1]*SubDiv[CLat2]); 01279 for(int n=0;n<SQR(NSample);n++){ 01280 PosUp[n] = 0.; 01281 PosDown[n] = 0.; 01282 CountUp[n] = 0.; 01283 CountDown[n] = 0.; 01284 } 01285 for(int p=0;p<pNPart();p++){ 01286 if(pType(p) != 1) continue; 01287 double x = pPos(p,CLat1) - floor(pPos(p,CLat1)*pInvEdge(CLat1))*pEdge(CLat1); 01288 double y = pPos(p,CLat2) - floor(pPos(p,CLat2)*pInvEdge(CLat2))*pEdge(CLat2); 01289 int sx = (int)(x*InvEdge[CLat1]*SubDiv[CLat1]); 01290 if(sx < 0 || sx >= SubDiv[CLat1]) continue; 01291 int sy = (int)(y*InvEdge[CLat2]*SubDiv[CLat2]); 01292 if(sy < 0 || sy >= SubDiv[CLat2]) continue; 01293 if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){ 01294 PosUp[sx*NSample+sy] += pPos(p,CNorm); 01295 CountUp[sx*NSample+sy] += 1.; 01296 } 01297 else { 01298 PosDown[sx*NSample+sy] += pPos(p,CNorm); 01299 CountDown[sx*NSample+sy] += 1.; 01300 } 01301 } 01302 double SumX = 0.; 01303 double SumX2 = 0.; 01304 double Count = 0.; 01305 for(int sx = 0;sx<SubDiv[CLat1];sx++){ 01306 for(int sy = 0;sy<SubDiv[CLat2];sy++){ 01307 int n = sx*NSample+sy; 01308 PosUp[n] /= CountUp[n] > 1. ? CountUp[n] : 1.; 01309 PosDown[n] /= CountDown[n] > 1. ? CountDown[n] : 1.; 01310 double Val = (PosUp[n]-PosDown[n]); 01311 SumX += Val; 01312 Count += 1.; 01313 ElCoup[3*t+2] += 1.; 01314 } 01315 } 01316 double NInv = Count > 0. ? 1./Count : 1.; 01317 double MeanVal = SumX*NInv; 01318 for(int sx = 0;sx < SubDiv[CLat1];sx++){ 01319 for(int sy = 0;sy < SubDiv[CLat2];sy++){ 01320 int n = sx*NSample+sy; 01321 double Val = (PosUp[n]-PosDown[n]); 01322 SumX2 += SQR(Val - MeanVal); 01323 } 01324 } 01325 ElCoup[3*t ] += MeanVal; 01326 ElCoup[3*t+1] += SumX2/(Count - 1.); 01327 if(1==0){//writes the single distribution files 01328 sprintf(FileName,"Distr%.0fKappa%02d.dat",pkappaN(),t); 01329 FILE *FileToWrite = fopen(FileName,"a"); 01330 for(int l1=0;l1<SubDiv[CLat1];l1++){ 01331 for(int l2=0;l2<SubDiv[CLat2];l2++){ 01332 int n = l1*NSample+l2; 01333 fprintf(FileToWrite,"%lf\n",(PosUp[n]-PosDown[n])); 01334 } 01335 } 01336 fclose(FileToWrite); 01337 } 01338 } 01339 } 01340 #ifdef OMPI_MPI_H 01341 MPI_Allreduce(MPI_IN_PLACE,ElCoup,3*NSample, MPI_DOUBLE, MPI_SUM, Proc->CommGrid); 01342 int Rank=0; 01343 MPI_Comm_rank(Proc->CommGrid, &Rank); 01344 if(Rank==0){ 01345 #endif 01346 printf("\n"); 01347 //--------normalizing-------------------- 01348 double Norm = 1./(double)(NFile[1]-NFile[0]); 01349 for(int t=0;t<NSample+1;t++){ 01350 ElCoup[3*t ] *= Norm; 01351 ElCoup[3*t+1] *= Norm; 01352 ElCoup[3*t+2] *= Norm; 01353 } 01354 ElCoup[3*NSample+1] = 0.; 01355 sprintf(FileName,"ElasticCoupling%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN()); 01356 FILE *FileToWrite = fopen(FileName,"w"); 01357 char cSystem[STRSIZE]; 01358 SysDef(cSystem); 01359 fprintf(FileToWrite,"# %s",cSystem); 01360 fprintf(FileToWrite,"# Area Thickness ThickStdDev Area*ThickStdDev\n"); 01361 for(int t=0;t<NSample+1;t++){ 01362 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 ])); 01363 } 01364 fclose(FileToWrite); 01365 #ifdef OMPI_MPI_H 01366 } 01367 #endif 01368 free(ElCoup); 01369 free(PosUp); 01370 free(PosDown); 01371 free(CountUp); 01372 free(CountDown); 01373 } 01377 void ElPoly::ElasticCouplingNVT(){ 01378 //mean StdDev Area*StdDev/Mean^2 01379 double ElCoup[3] = {0.,0.,0.}; 01380 char FileName[256]; 01381 char cSystem[STRSIZE]; 01382 //for the simulation in NPtT a smaller box is used 01383 sprintf(FileName,"ElCoupNVT%.0fRho%.0fChiN%.0fKappa.dat",prho(),pchiN(),pkappaN()); 01384 FILE *FileToWrite = fopen(FileName,"w"); 01385 SysDef(cSystem); 01386 fprintf(FileToWrite,"# %s",cSystem); 01387 fprintf(FileToWrite,"# step Thickness ThickStdDev\n"); 01388 for(int f=NFile[0];f<NFile[1];f++){ 01389 Processing(f); 01390 OpenRisk(cFile[f],BF_CHAIN); 01391 double PosUp = 0.; 01392 double PosDown = 0.; 01393 double CountUp = 0.; 01394 double CountDown = 0.; 01395 for(int p=0;p<pNPart();p++){ 01396 if(pType(p) != 1) continue; 01397 if(VAR_IF_TYPE(Ch[pChain(p)].Type,CHAIN_UP)){ 01398 PosUp += pPos(p,CNorm); 01399 CountUp += 1.; 01400 } 01401 else { 01402 PosDown += pPos(p,CNorm); 01403 CountDown += 1.; 01404 } 01405 } 01406 PosUp /= CountUp > 1. ? CountUp : 1.; 01407 PosDown /= CountDown > 1. ? CountDown : 1.; 01408 ElCoup[0] += (PosUp-PosDown); 01409 ElCoup[1] += SQR(PosUp-PosDown); 01410 fprintf(FileToWrite,"%d %lf\n",f,PosUp-PosDown); 01411 } 01412 printf("\n"); 01413 ElCoup[1] = sqrt((ElCoup[1] - SQR(ElCoup[0])*(NFile[1]-NFile[1]))/((double)(NFile[0]-NFile[1])-1.)); 01414 fprintf(FileToWrite,"# %lf %lf %lf\n",pEdge(CLat1)*pEdge(CLat2),ElCoup[0],ElCoup[1]); 01415 fclose(FileToWrite); 01416 } 01417 void ElPoly::BilayerDistance(char *OutFile,int NSample){ 01418 FILE *File2Write = fopen(OutFile,"w"); 01419 int NLayer = 2; 01420 double PosBf[3]; 01421 for(int f=NFile[0];f<NFile[1];f++){ 01422 Processing(f); 01423 if(OpenRisk(cFile[f],BF_CHAIN))return; 01424 double *PosLay = (double *)calloc(4*NSample*NSample,sizeof(double)); 01425 double *CountLay = (double *)calloc(4*NSample*NSample,sizeof(double)); 01426 for(int b=0,cOff=0,pOff=0;b<pNBlock();cOff+=pNChain(b++)){ 01427 int Level = 0; 01428 if(!strcmp(Block[b].Name,"LIPID1")) Level = 2; 01429 for(int c=cOff;c<pNChain(b)+cOff;c++,pOff+=pNPCh(b)){ 01430 for(int d=0;d<3;d++){ 01431 PosBf[d] = Ch[c].Pos[d] - floor(Ch[c].Pos[d]*pInvEdge(d))*pEdge(d); 01432 } 01433 int cLevel = Level; 01434 if(CHAIN_IF_TYPE(Ch[c].Type,CHAIN_UP)) cLevel += 1; 01435 int sx = (int)(PosBf[CLat1]*pInvEdge(CLat1)*NSample); 01436 int sy = (int)(PosBf[CLat2]*pInvEdge(CLat2)*NSample); 01437 if(sx < 0 || sx >= NSample) continue; 01438 if(sy < 0 || sy >= NSample) continue; 01439 PosLay[(sx*NSample+sy)*4+cLevel] += Ch[c].Pos[CNorm];//PosBf[CNorm]; 01440 CountLay[(sx*NSample+sy)*4+cLevel] += 1.; 01441 } 01442 } 01443 for(int sx = 0;sx < NSample;sx++){ 01444 for(int sy = 0;sy < NSample;sy++){ 01445 for(int l=0;l<4;l++){ 01446 PosLay[(sx*NSample+sy)*4+l] /= CountLay[(sx*NSample+sy)*4+l] > 0. ? CountLay[(sx*NSample+sy)*4+l] : 1.; 01447 } 01448 double h1 = PosLay[(sx*NSample+sy)*4+2]-PosLay[(sx*NSample+sy)*4+1]; 01449 double h2 = PosLay[(sx*NSample+sy)*4+3]-PosLay[(sx*NSample+sy)*4+0]; 01450 h2 -= floor(h2*pInvEdge(CNorm))*pEdge(CNorm); 01451 //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; 01452 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]); 01453 } 01454 } 01455 free(PosLay); 01456 free(CountLay); 01457 } 01458 printf("\n"); 01459 fclose(File2Write); 01460 } 01461 void ElPoly::BondDistr(char *OutFile,int NBin){ 01462 FILE *File2Write = fopen(OutFile,"w"); 01463 // double **Histo = (double **)calloc(3,sizeof(double)); 01464 // double **Raw = (double **)calloc(3,sizeof(double)); 01465 // for(int d=0;d<3;d++){ 01466 // Histo[d] = (double *)calloc(NBin,sizeof(double)); 01467 // Raw[d] = (double *)calloc(pNPart,sizeof(double)); 01468 // } 01469 for(int f=NFile[0];f<NFile[1];f++){ 01470 Processing(f); 01471 if(OpenRisk(cFile[f],BF_NO)) return; 01472 double Dist[3]; 01473 for(int c=0;c<pNChain();c++){ 01474 for(int p=0;p<pNPCh()-1;p++){ 01475 int p1 = p + c*pNPCh(); 01476 for(int d=0;d<3;d++){ 01477 Dist[d] = remainder(pPos(p1,d) - pPos(p1+1,d),pEdge(d)); 01478 } 01479 fprintf(File2Write,"%lf %lf %lf %lf\n",Dist[0],Dist[1],Dist[2],sqrt(SQR(Dist[0])+SQR(Dist[1])+SQR(Dist[2]))); 01480 } 01481 } 01482 } 01483 // MOMENTI Mom[3]; 01484 // int xMin = Raw[0][0]; 01485 // int xMax = Raw[0][0]; 01486 // for(int d=0;d<3;d++){ 01487 // Mom[d] = Mat->Distribuzione(Raw[d],pNPart-pNChain,Histo[d],NBin); 01488 // Mat->Normalize(Histo[d],NBin); 01489 // if(xMin > Mom[d].Minimo) xMin = Mom[d].Minimo; 01490 // if(xMax > Mom[d].Massimo) xMax = Mom[d].Massimo; 01491 // } 01492 // for(int v=0;v<NBin;v++){ 01493 // double x = v/(double)(NBin)*(xMax-xMin)+xMin; 01494 // fprintf(File2Write,"%lf %lf %lf %lf\n",x,Histo[0][v],Histo[1][v],Histo[2][v]); 01495 // } 01496 // for(int d=0;d<3;d++){ 01497 // free(Histo[d]); 01498 // free(Raw[d]); 01499 // } 01500 // free(Histo); 01501 // free(Raw); 01502 fclose(File2Write); 01503 } 01504 void ElPoly::EndToEndDist(){ 01505 double *EndToEnd = (double *) calloc(pNPart()*(NFile[1]-NFile[0]),sizeof(double)); 01506 double *Distr = (double *) calloc((NFile[1]-NFile[0]),sizeof(double)); 01507 for(int f=NFile[0];f<NFile[1];f++){ 01508 Processing(f); 01509 if(OpenRisk(cFile[f],BF_NO))return; 01510 01511 01512 } 01513 01514 }
1.7.6.1