Allink  v0.1
ForcesLoop.cpp
1 #include "Forces.h"
4  if(VAR_IF_TYPE(SysShape,SYS_2D)){
5  //CalcUpdate = &Forces::Sim2d;
6  }
7 }
9 void Forces::Solve(){
10  if(VAR_IF_TYPE(SysShape,SYS_2D))
12  //SolveLinks();
13  else if(VAR_IF_TYPE(SysShape,SYS_LEAVES))
14  SolveLeaves();
15  else if(VAR_IF_TYPE(SysShape,SYS_PORE))
16  SolveLeaves();
17  else if(VAR_IF_TYPE(SysShape,SYS_ROD))
18  SolveRod();
19 }
21  if(VAR_IF_TYPE(SysShape,SYS_MD)){
22  VelVerlet1();
23  }
24  if(VAR_IF_TYPE(SysShape,SYS_2D)){
25  Wave();
26  }
27  else if(VAR_IF_TYPE(SysShape,SYS_LEAVES)){
29  }
30  else if(VAR_IF_TYPE(SysShape,SYS_1D)){
32  }
33  else if(VAR_IF_TYPE(SysShape,SYS_3D)){
35  }
36  else if(VAR_IF_TYPE(SysShape,SYS_ROD)){
37  ForceFieldRod();
38  }
39  else if(VAR_IF_TYPE(SysShape,SYS_RIGID)){
42  }
43  else if(VAR_IF_TYPE(SysShape,SYS_MC)){
44  if(VAR_IF_TYPE(CalcMode,CALC_NVT)){
45  //for(int p=0;p<2*pNPCh();p++)
46  NInsertion += TryMove();
47  }
48  else if(VAR_IF_TYPE(CalcMode,CALC_NcVT)){
49  NInsertion += TryMoveCh();
50  }
51  else if(VAR_IF_TYPE(CalcMode,CALC_mVT)){
52  if(Mat->Casuale() < .5)
53  NRemoval += TryRemove();
54  else
55  NInsertion += TryInsert();
56  }
57  else if(VAR_IF_TYPE(CalcMode,CALC_mcVT)){
58  if(!VAR_IF_TYPE(CalcMode,CALC_CONF_BIAS)){
59  if(Mat->Casuale() < .5) NRemoval += TryRemoveCh();
60  else NInsertion += TryInsertCh();
61  }
62  else {
63  if(Mat->Casuale() < .5)
65  else
67  }
68  }
69  }
70  else if(VAR_IF_TYPE(SysShape,SYS_ELECTRO)){
71  if(VAR_IF_TYPE(CalcMode,CALC_NVT)){
72  NInsertion += TryMove();
73  }
74  else if(VAR_IF_TYPE(CalcMode,CALC_mVT)){
75  if(Mat->Casuale() < .5)
76  NRemoval += TryRemove();
77  else
78  NInsertion += TryInsert();
79  }
80  }
81  else if(VAR_IF_TYPE(SysShape,SYS_MD)){
82  OldNrgSys = SumForcesMD();
83  Pc->Erase();
84  for(int p=0;p<pNPart();p++){
85  Pc->AddPart(p,Pm[p].Pos);
86  }
87  }
88  if(VAR_IF_TYPE(SysShape,SYS_MD)){
89  ApplyTherm();
90  VelVerlet2();
91  }
92  Time += pDeltat();
93  IncrStep();
94  SetTime(Time);
95  Task();
96 }
98  Shout("Explore pep size\n");
99  if(SysShape != SYS_LEAVES) return ;
100  char FName[120];
101  double HeiMin = .15;
102  double HeiMax = .4;
103  double HeiStep = (HeiMax-HeiMin)/(double)NGrid;
104  double AngMin = 5.;
105  double AngMax = 35.;
106  double AngStep = (AngMax-AngMin)/(double)NGrid;
107  double SLapMin = 0.001;
108  double SLapMax = 10.0;
109  double SLapStep = 10.;
110  for(double SLap=SLapMin;SLap<SLapMax;SLap*=SLapStep){
111  sprintf(FName,"PepHeiElThickSLap%lf.dat",SLap);
112  FILE *FWrite = fopen(FName,"w");
113  for(double Hei=HeiMin;Hei<HeiMax;Hei+=HeiStep){
114  for(double Angle=AngMin;Angle<AngMax;Angle+=AngStep){
115  fprintf(stderr,"%lf %lf %lf\r",SLap,Hei,Angle);
116  Nano->Height = Hei;
117  Kf.SLap = SLap;
118  Nano->Hamaker = Angle;
119  Solve();
120  int c = 0;
121  double Min = pEdge(2);
122  double Max = 0.;
123  double xMin = 0.;
124  double xMax = 0.;
125  for(int p = c*pNPCh();p<(c+1)*pNPCh();p++){
126  if(Pm[p].Pos[2] > Max && Pm[p].Typ == 0){
127  Max = Pm[p].Pos[2];
128  xMax = Pm[p].Pos[0];
129  }
130  }
131  c = 1;
132  for(int p = c*pNPCh();p<(c+1)*pNPCh();p++){
133  if(Pm[p].Pos[2] < Min && Pm[p].Typ == 0){
134  Min = Pm[p].Pos[2];
135  xMin = Pm[p].Pos[0];
136  }
137  }
138  fprintf(FWrite,"%lf %lf %lf %lf %lf\n",Hei/.2,Angle,(Min-Max)/.2,SLap,.5*(xMin+xMax)/.2);
139  }
140  }
141  fclose(FWrite);
142  }
143  printf("\n");
144 }
146  Shout("Explore pep size 2d\n");
147  if(SysShape != SYS_2D) return ;
148  SetEdge(.5*MIN(pEdge(0),pEdge(1)),3);
149  char FName[120];
150  double HeiMin = 0.2;
151  double HeiMax = 4.2;
152  double HeiStep = (HeiMax-HeiMin)/(double)NGrid;
153  double AngMin = 0.;
154  double AngMax = 80.;
155  double AngStep = (AngMax-AngMin)/(double)NGrid;
156  double InvNBin = 1./(double)NBin;
157  SysFormat = VAR_SYS_TXVL;
158  for(double Angle=AngMin;Angle<AngMax;Angle+=AngStep){
159  for(double Hei=HeiMin;Hei<HeiMax;Hei+=HeiStep){
160  fprintf(stderr,"%lf %lf\r",Hei,Angle);
161  Nano->Height = Hei;
162  Nano->Hamaker = Angle;
163  Solve();
164  sprintf(FName,"Sol2dHei%.1fAng%02d.dat",Hei,(int)Angle);
165  Write(FName);
166  }
167  }
168  printf("\n");
169 }
171  Shout("Explore pep size\n");
172  char String[120];
173  if(SysShape != SYS_2D) return;
174  if(pNNano() < 2) return;
175  double DistMin = pNanoPos(1,0);
176  double DistDelta = (pEdge(0)-pNanoPos(1,0))/(double)NGrid;
177  double DeltaGrid = pEdge(0)/(double)nEdge[0];
178  // DistDelta = floor(2.*DistDelta/DeltaGrid)*DeltaGrid;
179  // SigErr(DistDelta <= 0.,"The grid step is too short for %d distances ",NGrid);
180  char FName[60];
181  FILE *FProf = fopen("DistThinMin.dat","w");
182  for(int i=0;i<NGrid;i++){
183  fprintf(stderr,"done %.2f %%\r",100.*i/(double)NGrid);
184  IfFillMatrix = 1;
185  Solve();
186  sprintf(FName,"MinProf%03d.dat",i);
187  FILE *FWrite = fopen(FName,"w");
188  StringNano(String,0);
189  fprintf(FWrite,String);
190  StringNano(String,1);
191  fprintf(FWrite,String);
192  double Min = 100000.;
193  double Dist = 0.;
194  double zOld = 100000.;
195  int IfContinue = 0;
196  int nBin = 0;
197  for(int p=nEdge[1]/2+1;p<pNPart();p+=nEdge[1]){
198  double x = pPos(p,0) - pNanoPos(0,0);
199  double z = pPos(p,2);
200  if(pType(p) == 1) z = pNanoPos(0,2) + .5*Nano[0].Height;
201  if( x > (pNanoPos(1,0) - pNanoPos(0,0) )) continue;
202  fprintf(FWrite,"%lf %lf\n",x,z);
203  if( x > .6*(pNanoPos(1,0) - pNanoPos(0,0) )) continue;
204  if(Min > z){// && IfContinue){
205  Min = z;
206  Dist = x;
207  nBin = p;
208  }
209  // if(z > zOld){IfContinue = 1;}
210  // zOld = z;
211  }
212  fprintf(FProf,"%lf %lf %lf\n",Nano[1].Pos[0]-Nano[0].Pos[0],Dist,Min);
213  fflush(FProf);
214  fclose(FWrite);
215  sprintf(FName,"Sol2dDist%.3f.dat",Nano[1].Pos[0]-Nano[0].Pos[0]);
216  WriteTxvl(FName);
217  Nano[1].Pos[0] += DistDelta;
218  if(Nano[1].Pos[0] > pEdge(0)) return;
219  }
220  printf("\n");
221  fclose(FProf);
222 }
223 void Forces::RunWidom(char *File2Open,int f){
224  Shout("Widom on particles\n");
225  int NInt = 10000;
226  double NrgDiff = 0.;
227  char File2Write[60];
228  sprintf(File2Write,"WidomOut%05d.dat",f);
229  ReOpen(File2Open,BF_PART);
230  FILE *WidomOut = fopen(File2Write,"w");
231  SigErr(WidomOut==NULL,"Cannot allocate %s\n",File2Write);
232  for(int p=0;p<pNPart();p++){
233  WidomRemove(&NrgDiff,p);
234  fprintf(WidomOut,"%lf\n",NrgDiff);
235  }
236  fclose(WidomOut);
237  sprintf(File2Write,"WidomIn%05d.dat",f);
238  FILE *WidomIn = fopen(File2Write,"w");
239  for(int p=0;p<NInt;p++){
240  WidomInsert(&NrgDiff);
241  fprintf(WidomIn,"%lf\n",NrgDiff);
242  }
243  fclose(WidomIn);
244 }
245 void Forces::RosenIn(FILE *WidomIn){
246  Shout("Widom insertion biased/Rosenbluth\n");
247  int NInt = 10000;
248  double Weight = 0.;
249  for(int p=0;p<NInt;p++){
250  WidomBiasChIn(&Weight);
251  fprintf(WidomIn,"%g\n",Weight);
252  }
253 }
254 void Forces::RosenOut(FILE *WidomIn){
255  Shout("Widom deletion biased/Rosenbluth\n");
256  int NInt = 10000;
257  double Weight = 0.;
258  for(int c=0;c<pNChain();c++){
259  WidomBiasChOut(&Weight,c);
260  fprintf(WidomIn,"%g\n",Weight);
261  }
262 }
263 void Forces::RunWidomChIn(char *File2Open,int f){
264  Shout("Widom chain in\n");
265  int NInt = 5*pNChain();
266  double NrgDiff[3];
267  char File2Write[60];
268  ReOpen(File2Open,BF_PART);
269  sprintf(File2Write,"WidomIn%05d.dat",f);
270  FILE *WidomIn = fopen(File2Write,"w");
271  //StudySys();
272  for(int p=0;p<NInt;p++){
273  IncrStep();
274  WidomInsertCh(NrgDiff);
275  fprintf(WidomIn ,"%lf\n",NrgDiff[2]);
276  }
277  fclose(WidomIn);
278 }
279 void Forces::RunWidomChOut(char *File2Open,int f){
280  Shout("Widom chain out\n");
281  double NrgDiff[3];
282  char File2Write[60];
283  sprintf(File2Write,"WidomOut%05d.dat",f);
284  ReOpen(File2Open,BF_PART);
285  FILE *WidomOut = fopen(File2Write,"w");
286  CalcTotNrgCh();
287  for(int c=0;c<pNChain();c++){
288  WidomRemoveCh(NrgDiff,c);
289  fprintf(WidomOut,"%lf\n",NrgDiff[2]);
290  }
291  fclose(WidomOut);
292 }
293 void Forces::CalcTotNrg(char *File2Open,int f){
294  Shout("Calculating total energy\n");
295  ReOpen(File2Open,BF_PART);
296  double NrgNb = 0.;
297  double NrgSpr = 0.;
298  double NrgSpr2 = 0.;
299  double NrgBend = 0.;
300  for(int c=0;c<pNChain();c++){
301  for(int p=c*pNPCh();p<(c+1)*pNPCh()-1;p++){
302  double DistBA[3];
303  double DistCB[3];
304  double DistBA2 = 0.;
305  double DistCB2 = 0.;
306  double CosAngle = 0.;
307  for(int d=0;d<3;d++){
308  DistBA[d] = Pm[p].Pos[d] - Pm[p-1].Pos[d];
309  DistBA[d] -= floor(DistBA[d]/(pEdge(d)) + .5)*pEdge(d);
310  DistCB[d] = Pm[p+1].Pos[d] - Pm[p].Pos[d];
311  DistCB[d] -= floor(DistCB[d]/(pEdge(d)) + .5)*pEdge(d);
312  DistBA2 += SQR(DistBA[d]);
313  DistCB2 += SQR(DistCB[d]);
314  CosAngle += DistBA[d]*DistCB[d];
315  }
316  DistCB2 = sqrt(DistCB2);
317  NrgSpr += .5*pkSpr()*SQR(DistCB2 - pSprRest());
318  if(p == c*pNPCh()) continue;
319  DistBA2 = sqrt(DistBA2);
320  CosAngle /= (DistBA2*DistCB2);
321  NrgBend += pkBen()*(1.-CosAngle);
322  }
323  }
324  double NrgNano = 0.;
325  for(int p=0;p<pNPart();p++){
326  NrgNano += NanoNrg(p);
327  }
328  CalcTotNrgCh();
329  printf("Nb %lf + Nano %lf = %lf Spr %lf + Ben %lf = %lf\n",OldNrgSys,NrgNano,OldNrgSys+NrgNano,NrgSpr,NrgBend,NrgSpr+NrgBend);
330  fprintf(StatFile1,"%d %lf %lf %lf %lf\n",f,OldNrgSys,NrgSpr,NrgBend,NrgNano);
331  fflush(StatFile1);
332 }
333 void Forces::CalcNrgPep(char *File2Open,int f){
334  Shout("Calculating total energy\n");
335  ReOpen(File2Open,BF_PART);
336  AddDens(0,pNPart());
337  OldNrgSys = SumDens(0,pNPart());
338  ClearDens();
339  for(int b=0,NPep=0,cOff=0,pOff=0;b<pNBlock();cOff+=pNChain(b++)){
340  if(!strncmp(Block[b].Name,"PEP",3)){
341  continue;
342  }
343  for(int c=cOff;c<cOff+pNChain(b);c++,pOff+=pNPCh(b)){
344  for(int p=pOff,link=0;p<MIN(pOff+pNPCh(b),pNPart());p++){
345  AddDens(p,p+1);
346  }
347  }
348  }
349  double NrgPep = OldNrgSys - SumDens(0,pNPart());
350  printf("Nrg pep %lf\n",NrgPep);
351  fprintf(StatFile1,"%d %lf %lf\n",f,NrgPep,OldNrgSys);
352  fflush(StatFile1);
353 }
354 void Forces::CalcTens(char **argv,int *FilePos,int NFile){
355  Shout("Calculating tension\n");
356  AllocTens();
357  char FName[60];
358  for(int f=0;f<NFile;f++){
359  fprintf(stderr,"Elaborating file %s %.3f %%\r",argv[FilePos[f]],f/(double)NFile*100.);
360  ReOpen(argv[FilePos[f]],BF_PART);
361  CalcTens();
362  CalcDens();
363  if( !(f%(NWrite)) ){
364  for(int c=0;c<Tens.NComp;c++){
365  if(VAR_IF_TYPE(Tens.CalcMode,CALC_2d))
366  sprintf(FName,"Tension2d%05dL%d.dat",f,c);
367  if(VAR_IF_TYPE(Tens.CalcMode,CALC_3d))
368  sprintf(FName,"Tension3d%05dL%d.dat",f,c);
369  WriteTens(FName,c,1./(double)NFile);
370  }
371  }
372  }
373  printf("\n");
374 }
375 void Forces::AvForces(char **argv,int *FilePos,int NFile){
376  Shout("Calculating force field/average force\n");
377  AllocTens();
378  char FName[60];
379  const int NBin = 120;
380  double *Profile = (double *)calloc(3*NBin,sizeof(double));
381  for(int f=0;f<NFile;f++){
382  fprintf(stderr,"Elaborating file %s %.3f %%\r",argv[FilePos[f]],f/(double)NFile*100.);
383  ReOpen(argv[FilePos[f]],BF_PART);
385  FILE *FForce = fopen("ForceProfile.dat","w");
386  double Force[3] = {0.,0.,0.};
387  int cBin[3];
388  for(int p=0;p<pNPart();p++){
389  for(int d=0;d<3;d++){
390  cBin[d] = (int)(Pm[p].Pos[d]*pInvEdge(d)*NBin);
391  if(cBin[d] < 0 || cBin[d] >= NBin) continue;
392  Profile[cBin[d]*3+d] += Fm[p].Dir[d];
393  Force[d] += Fm[p].Dir[d];
394  fprintf(FForce,"%d %g %g %g\n",p,Force[0],Force[1],Force[2]);
395  }
396  }
397  printf("%g %g %g\n",Force[0],Force[1],Force[2]);
398  fclose(FForce);
399  }
400  FILE *FProf = fopen("ForceProfile.dat","w");
401  for(int b=0;b<NBin;b++){
402  fprintf(FProf,"%d %lf %lf %lf\n",b,Profile[b*3],Profile[b*3+1],Profile[b*3+2]);
403  }
404  FILE *FField = fopen("ForceField.dat","w");
405  double Delta = sqrt(Kf.CutOff2)/(double)NTab;
406  for(int b=0;b<NTab;b++){
407  double x = b*Delta;
408  fprintf(FField,"%lf ",x);
409  for(int t1=0;t1<pNType();t1++){
410  for(int t2=0;t2<pNType();t2++){
411  fprintf(FField,"%lf ",FTab[(b*pNType()+t1)*pNType()+t2]);
412  }
413  }
414  fprintf(FField,"\n");
415  }
416  free(Profile);
417  fclose(FProf);
418  fclose(FField);
419  printf("\n");
420 }
422  if(VAR_IF_TYPE(SysAlloc,ALL_MC)){
423  return;
424  if((pStep()%NUpdate)!=0) return;
425  OldNrgSys = 0.;
426  double Pot[3];
427  for(int c=0;c<pNChain();c++) OldNrgSys += CalcNrgCh(c,Pot);
428  //fprintf(TempFile,"%d %d %lf \n",pStep(),pNChain(),NrgSys/(double)pNChain());
429  //fprintf(StatFile1,"%d %d %lf \n",pStep(),pNChain(),OldNrgSys/(double)pNChain());
430  //fflush(StatFile1);
431  // int c = 0;//p%pNPCh();
432  // double Pot[3];
433  // CalcNrgCh(c,Pot);
434  // for(int e=0;e<3;e++) OldNrgCh[c*3+e] = Pot[e];
435  // double Dist[4];
436  // TwoPartDist(0,pNPCh()-1,Dist);
437  // fprintf(TempFile,"%d %lf %lf %lf %lf\n",pStep(),OldNrgCh[3*0],OldNrgCh[3*0+1],OldNrgCh[3*0+2],Dist[3]);
438  }
439  else if(VAR_IF_TYPE(SysAlloc,ALL_MD)){
440  if( !(pStep()%NUpdate) ){
441  double v2 = 0.;
442  for(int p=0;p<pNPart();p++){
443  for(int d=0;d<3;d++){
444  v2 += SQR(Pm[p].Vel[d]);
445  }
446  }
447  //fprintf(StatFile1,"%d %lf %lf\n",pStep(),v2/(double)(3*pNPart()),OldNrgSys);
448  //fflush(StatFile1);
449  }
450  }
451 }
453  CheckPairList();
454 }
456  Shout("run dynamics\n");
457  double NLoopSec = 0.;
458  char FileName[60];
459  DefForceParam();
460  //MinimalNrg();
461  for(int s=0;s<SimLimit;s++){
462  Dynamics();
463  //MinimalMD();
464  if( !(s%NUpdate) ){
465  CurrTime = time(NULL);
466  NLoopSec += (s)/(double)(CurrTime-InitTime);
467  double Temp = 0.;
468  for(int p=0;p<pNPart();p++){
469  for(int d=0;d<3;d++){
470  Temp += SQR(Pm[p].Vel[d]);
471  }
472  }
473  Temp = Temp/(double)(3*pNPart());
474  fprintf(stderr,"NPart %d loop/ms %.3g acc/step %.3f in/out %.4f T %.4f accomplished %.3f %% Nrg %lf\n",pNPart(),NLoopSec,(NRemoval+NInsertion)/(double)s,NInsertion/(double)NRemoval,Temp,s/(double)(SimLimit)*100.,OldNrgSys);
475  fprintf(StatFile1,"%d %lf %lf\n",s,OldNrgSys,Temp);
476  }
477  if( !(s%(NWrite)) ){
478  sprintf(FileName,"Trajectory%09d.dat",pStep());
479  Write(FileName);
480  }
481  }
482  printf("\n");
483  CurrTime = time(NULL);
484  NLoopSec += (SimLimit)/(double)(CurrTime-InitTime);
485  double v2 = 0.;
486  for(int p=0;p<pNPart();p++){
487  for(int d=0;d<3;d++){
488  v2 += SQR(Pm[p].Vel[d]);
489  }
490  }
491  fprintf(StatFile1,"##NPart %d loop/ms %.3g acc/step %.3f in/out %.4f T %.4f accomplished Nrg %lf\r",pNPart(),NLoopSec,(NRemoval+NInsertion)/(double)pStep(),NInsertion/(double)NRemoval,v2/(double)(3*pNPart()),OldNrgSys);
492 }
494  double Pos[3];
495  double Pot[3];
496  double DistRel[4];
497  OldNrgSys = 0.;
498  Pc->Erase();
499  for(int p=0;p<pNPart();p++){pPos(p,Pos);Pc->AddPart(p,Pos);}
500  for(int p1=0;p1<pNPart();p1++){
501  if(VAR_IF_TYPE(SysAlloc,ALL_FORCES)){
502  for(int d=0;d<3;d++) Fm[p1].Dir[d] = 0.;
503  }
504  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
505  int p2 = Pc->p2Curr;
506  if(p2 <= p1) continue;
507  Pc->Dist2Curr(DistRel);
508  if(DistRel[3] > Kf.CutOff2) continue;
509  double Cons = Potential(DistRel[3],pType(p1),pType(p2),Pot);
510  double InvDist = DistRel[3];
511  if(VAR_IF_TYPE(SysAlloc,ALL_FORCES)){
512  for(int d=0;d<3;d++){
513  Fm[p1].Dir[d] += Cons*DistRel[d]*InvDist;
514  Fm[p2].Dir[d] -= Cons*DistRel[d]*InvDist;
515  }
516  }
517  OldNrgSys += Pot[0];
518  }
519  }
520 }
522  double Sigma = sqrt(pTemp());
523  double Pos[3];
524  for(int p=0;p<pNPart();p++){
525  for(int d=0;d<3;d++){
526  Pm[p].Vel[d] += .5*Fm[p].Dir[d]*pDeltat();
527  Pm[p].Pos[d] += Pm[p].Vel[d]*pDeltat();
528  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
529  }
530  }
531  MinimalNrg();
532  //vv2
533  double Temp = 0.;
534  for(int p=0;p<pNPart();p++){
535  for(int d=0;d<3;d++){
536  Pm[p].Vel[d] += .5*Fm[p].Dir[d]*pDeltat();
537  Temp += SQR(Pm[p].Vel[d]);
538  }
539  }
540  //Andersen
541  for(int p=0;p<pNPart();p++){
542  if(Mat->Casuale() < pDeltat()){
543  for(int d=0;d<3;d++){
544  Pm[p].Vel[d] = Mat->Gaussiano(0.,Sigma);//Mat->Gaussiano(0.,Sigma);
545  }
546  }
547  }
548  Temp = Temp/(3*pNPart());
549  double Norm = 1./sqrt(Temp);
550  // for(int p=0;p<pNPart();p++){
551  // for(int d=0;d<3;d++){
552  // Pm[p].Vel[d] *= Norm;
553  // }
554  // }
555 }
557  double *Sol2d = (double *)calloc(SQR(NEdge),sizeof(double));
558  double *Count = (double *)calloc(SQR(NEdge),sizeof(double));
559  double NInvEdge = 1./(double)NEdge;
560  for(int p=0;p<pNPart();p++){
561  if(pType(p) != 0) continue;
562  int vx = (Pm[p].Pos[0]*NEdge*pInvEdge(0));
563  if(vx < 0 || vx >= NEdge) continue;
564  int vy = (Pm[p].Pos[1]*NEdge*pInvEdge(1));
565  if(vy < 0 || vy >= NEdge) continue;
566  Sol2d[vx*NEdge+vx] = Pm[p].Pos[2];
567  Count[vx*NEdge+vx] = 1.;
568  }
569  for(int vx=0;vx<SQR(NEdge);vx++){
570  double Norm = Count[vx] > 0. ? 1./Count[vx] : 1.;
571  Sol2d[vx] *= Norm;
572  }
573  Create2d();
574  int NInt = 1000;
575  for(int n=0;n<NInt;n++){
576 
577  }
578 }
int TryMoveCh()
Trial desplacement of a chain.
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
void RunDynamics()
Run a step further.
Definition: ForcesLoop.cpp:455
void CalcTens()
Calculate the forces for the tension profile.
Definition: ForcesTens.cpp:85
void WidomBiasChIn(double *Weight)
Widom with Rosenbluth weight.
double pkBen()
Bending coupling.
Definition: VarData.h:932
void ExplorePepSize2d()
Find the minimun bilayer thickness for different peptide sizes.
Definition: ForcesLoop.cpp:145
void CalcTotNrg(char *FName, int nFile)
Total energy of the system.
Definition: ForcesLoop.cpp:293
double CalcTotNrgCh()
Calculate and sum up the energy of the chains.
void Trial()
Trial loop.
Definition: ForcesLoop.cpp:452
void VelVerlet1()
First step of the velocity Verlet.
int NInsertion
Count accepted moves.
Definition: Forces.h:767
int ForceFieldLine()
Helfrich Hamiltonian for a line.
void WriteTens(char *TFile, int Comp, double InvNFile)
Write the pressure and density profile.
Definition: ForcesTens.cpp:324
void Create2d()
Create a plane of connected beads.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void Erase()
Erase the pairlist.
Definition: Cubo.cpp:258
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
int TryRemoveChBias()
Trial biased removal of a chain.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double Hamaker
Strength of the interaction.
Definition: VarData.h:447
void RunWidom(char *File2Read, int f)
Build the widom histograms.
Definition: ForcesLoop.cpp:223
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
void StringNano(char *NString, int n)
String for the rigid inclusion in the header file.
void CheckPairList()
Check the pair list.
void SetTime(double Val)
Set Time.
Definition: VarData.h:1006
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
double Height
Height of the cylinder.
Definition: VarData.h:449
int pType(int p)
Return the type.
void Dynamics()
Sum up all the forces and update the positions.
Definition: ForcesLoop.cpp:20
int TryInsertCh()
Trial insertion of a chain.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
void RunWidomChIn(char *File2Read, int f)
Build the widom histograms.
Definition: ForcesLoop.cpp:263
double SumForcesMD()
Iterate all over the particles and calculate the forces.
void WidomInsert(double *NrgDiff)
Widom insertion.
double Dir[3]
Direction.
Definition: Forces.h:178
int NWrite
How many timesteps before write the snapshot.
Definition: Forces.h:775
double pkSpr()
Spring coupling.
Definition: VarData.h:934
double pSprRest()
Rest distance of the harmonic potential.
Definition: VarData.h:942
void CalcNrgPep(char *File2Open, int f)
Exchange energy of the protein.
Definition: ForcesLoop.cpp:333
void ForceFieldRod()
Bending potential on a rod.
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
void SolveLinksIterative()
Solve a system of four oder differential equation of particles connected by links.
int TryInsertChBias()
Trial biased insertion of a chain.
void WidomInsertCh(double *NrgDiff)
Widom insertion.
int NFile[2]
First and last file of the list.
Definition: Forces.h:771
int SysAlloc
Which arrays are allocated.
Definition: Forces.h:753
void WidomRemoveCh(double *NrgDiff, int c)
Widom removal.
int IfFillMatrix
If the matrix has to be changed.
Definition: Forces.h:763
void CalcForcesDensFunc()
Calculate the forces for the density functional.
int NRemoval
Count accepted moves.
Definition: Forces.h:769
int SimLimit
Maximum number of time steps.
Definition: Forces.h:751
void ExploreDoubleMin()
Find the minimum for different interpeptide distances.
Definition: ForcesLoop.cpp:170
int CalcMode
Calculation mode.
Definition: Forces.h:747
void ApplyTherm()
Pointer to the energy function.
Definition: Forces.h:441
double Pos[3]
Position.
Definition: VarData.h:427
void Shout(const char *s,...)
Internal message.
Definition: Forces.cpp:3
double Time
Total time.
Definition: Forces.h:733
double pNanoPos(int n, int d)
Return back folded nano position.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
void VelVerletRigid()
Velocity Verlet for a rigid body, first step.
int pNPCh()
Number of particle per chain.
void ForceFieldRigid()
Interaction between rigid bodies.
int pNBlock()
Number of blocks.
void RosenOut(FILE *WidomIn)
Rosenbluth histograms for deletion.
Definition: ForcesLoop.cpp:254
double NanoNrg(int p)
Exchange energy with the nano.
int TryMove()
Trial movement.
void ReOpen(char *FName, int Bf)
Open a new file.
Definition: Forces.cpp:901
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double Casuale()
Random uniform number.
bool WriteTxvl(char *OutFile)
Writes a "system-file" or a "x y z" file".
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
int NUpdate
How many timesteps before redrawing.
Definition: Forces.h:773
void NextCurr()
Increase the iterator to the next couple.
Definition: Cubo.cpp:479
double MinimalNrg()
Minmal nrg.
Definition: ForcesLoop.cpp:493
int ForceFieldLeaves()
Helfrich Hamiltonian with an elastic coupling.
int TryRemove()
Trial removal.
void ExplorePepSize()
Find the minimun bilayer thickness for different peptide sizes.
Definition: ForcesLoop.cpp:97
int SysShape
Shape of system.
Definition: Forces.h:745
void RosenIn(FILE *WidomIn)
Rosenbluth weights for insertion.
Definition: ForcesLoop.cpp:245
int pNNano()
Number of nanoparticles.
void RunWidomChOut(char *File2Read, int f)
Build the widom histograms.
Definition: ForcesLoop.cpp:279
void MinimalMD()
Minmal md.
Definition: ForcesLoop.cpp:521
int TryRemoveCh()
Trial removal of a chain.
int nEdge[3]
Number of particle per edge.
Definition: Forces.h:743
void IncrStep()
Increment Step.
Definition: VarData.h:1008
int ForceFieldBulk()
Armonic potential on a lattice.
void MinimizeSol()
Iterative process to approach to the solution.
Definition: ForcesLoop.cpp:556
double Potential(double Dist, int t1, int t2, double *Pot)
Pointer to a potential.
Definition: Forces.h:504
void SolveLeaves()
Solve a system of four oder differential equation of particles in a line.
int CalcMode
Calculation mode.
Definition: Forces.h:197
void ClearDens()
Set the local densities to zero.
double Gaussiano(double Media, double Scarto)
Gaussian random number.
int SysFormat
Contains the definition of the file format.
Definition: VarData.h:1088
void Wave()
Sinusoidal surface wave.
double CalcNrgCh(int c, double *Pot)
Pointer to the chain energy function.
Definition: Forces.h:513
void AllocTens()
Alloc the pressure profile.
Definition: ForcesTens.cpp:2
void VelVerlet2()
Second step of the velocity Verlet.
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int pNType()
of types of the particle
int NComp
of components
Definition: Forces.h:191
void Task()
Perform a operation every time step.
Definition: ForcesLoop.cpp:421
double pPos(int p, int d)
Return back folded position.
double pTemp()
Temperature.
Definition: VarData.h:912
int pNChain()
Number of chain.
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
void WidomRemove(double *NrgDiff, int p)
Widom removal.
void SolveRod()
Solve a system of four oder differential equation of particles connected by links.
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
int p2Curr
Current particle.
Definition: Cubo.h:51
double SumDens(int pInit, int pEnd)
Sum the local density for the particles between pInit and pEnd and multiply the factors by the virial...
void ChooseSimMode()
Choose the simulation method.
Definition: ForcesLoop.cpp:3
void CalcDens()
Calculate the densities.
Definition: ForcesTens.cpp:45
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd.
void SetCurr(int p)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:453
double pDeltat()
Delta t.
Definition: VarData.h:910
TENS Tens
Structure for the pressure calculation.
Definition: Forces.h:783
DomDec * Pc
Pair list.
Definition: Forces.h:785
double SLap
Prefactor of the square laplacian.
Definition: Forces.h:149
void AvForces(char **argv, int *FilePos, int NFile)
Average of the forces.
Definition: ForcesLoop.cpp:375
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
int TryInsert()
Trial insertion.
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
int pStep()
Number of steps.
int pNPart()
Number of particle.
void DefForceParam()
Define the parameters for calculating the force.
void WidomBiasChOut(double *Weight, int c)
Widom with Rosenbluth weight.
void Solve()
Solve a system of four oder differential equation.
Definition: ForcesLoop.cpp:9