Allink  v0.1
Forces.cpp
1 #include "Forces.h"
2 #include <stdarg.h>
3 void Forces::Shout(const char * s, ...){
4 #ifndef DIN_DEBUG
5  va_list args;
6  va_start(args, s);
7  fprintf(stderr, ">");
8  vfprintf(stderr, s, args);
9  fprintf(stderr, "\n");
10  va_end(args);
11 #else
12  return;
13 #endif
14 }
15 Forces::Forces(int argc,char **argv,int NInEdge,char *ConfFileExt){
16  Shout("Constructor/no starting configuration");
17  InitConst();
18  sprintf(ConfFile,"%s",ConfFileExt);
19  char ConfF[40];
20  sprintf(ConfF,ConfFile);
21  if(ReadConfDinamica(ConfF)){
22  Dx = 1./(double)(NEdge-1);
23  }
24  if(VAR_IF_TYPE(SysShape,SYS_2D)){
25  nEdge[0] = NEdge;
26  double Ratio = pEdge(1)*pInvEdge(0);
27  nEdge[1] = (int)(nEdge[0]*Ratio + 0.0001);
28  for(int i=0;;i++){
29  if((int)(nEdge[1]/Ratio) == nEdge[0]) break;
30  nEdge[0]++;
31  nEdge[1] = (int)(nEdge[0]*Ratio + 0.0001);
32  printf("using: nEdge[0] %d nEdge[1] %d\n",nEdge[0],nEdge[1]);
33  if(i>=10){
34  printf("Could not find the appropriate border partition %d %d\n",nEdge[0],nEdge[1]);
35  exit(0);
36  }
37  }
38  SetNLink(4);
39  ReSetNPart(nEdge[0]*nEdge[1]);
40  ReSetNChain(nEdge[1]);
41  SetNPCh(nEdge[0]);
42  }
43  else if(VAR_IF_TYPE(SysShape,SYS_3D)){
44  SetNLink(6);
46  ReSetNPart(NEdge*NEdge*NEdge);
47  }
48  else if(VAR_IF_TYPE(SysShape,SYS_ROD)){
49  SetNLink(1);
50  ReSetNChain(1);
52  }
53  else if(VAR_IF_TYPE(SysShape,SYS_STALK)){
54  SetNLink(3);
55  ReSetNPart(4*NEdge);
56  ReSetNChain(4);
57  // Kf.El[0] *= Gen->NPart;
58  }
59  else if(VAR_IF_TYPE(SysShape,SYS_LEAVES)){
60  SetNLink(3);
62  ReSetNChain(1);
63  // Kf.El[0] *= Gen->NPart;
64  }
65  else if(VAR_IF_TYPE(SysShape,SYS_PORE)){
66  SetNLink(3);
68  ReSetNChain(1);
69  }
70  else if(VAR_IF_TYPE(SysShape,SYS_1D)){
71  SetNLink(2);
73  ReSetNChain(1);
74  }
75  else if(VAR_IF_TYPE(SysShape,SYS_RIGID)){
76  ReSetNPart(0);
77  SetNLink(0);
78  }
79  else if(VAR_IF_TYPE(SysShape,SYS_TRIAL)){
80  ReSetNPart(3*3*3);
81  SetNLink(0);
82  ReSetNChain(1);
83  }
84  else if(VAR_IF_TYPE(SysShape,SYS_MD)){
85  SetNLink(0);
87  ReSetNChain(1);
88  }
89  else if(VAR_IF_TYPE(SysShape,SYS_MC)){
90  SetNLink(0);
92  ReSetNChain(1);
93  }
94  else if(VAR_IF_TYPE(SysShape,SYS_ELECTRO)){
95  SetNLink(4);
97  ReSetNChain(1);
98  }
100  //SetNPCh(NEdge);
101  SetDeltat(Deltat);
102  SetStep(0);
103  SetNType(2);
104  CreateInitial();
105  {
106  MInt = new MatInt(3,3);
107  MInt->SetCoeff(-24.33,0,0);
108  MInt->SetCoeff(-7.22,0,1);
109  MInt->SetCoeff(-24.33,0,2);
110  MInt->SetCoeff(-0.1,1,1);
111  MInt->SetCoeff(-7.22,1,2);
112  MInt->SetCoeff(0.,2,2);
113  MInt->SetCoeff(3.,0,0,0);
114  MInt->SetCoeff(3.,0,0,1);
115  MInt->SetCoeff(3.,0,0,2);
116  MInt->SetCoeff(3.,0,1,1);
117  MInt->SetCoeff(3.,0,1,2);
118  MInt->SetCoeff(3.,0,2,2);
119  MInt->SetCoeff(0.,1,1,1);
120  MInt->SetCoeff(3.,1,1,2);
121  MInt->SetCoeff(3.,1,2,2);
122  MInt->SetCoeff(0.,2,2,2);
123  }
124  AllocMethod();
125  PrepareSys();
126  PrepareParallel(argc,argv);
127  if(VAR_IF_TYPE(SysShape,SYS_LEAVES)) IfNano = 2;
128  if(VAR_IF_TYPE(SysShape,SYS_PORE)) IfNano = 2;
129  if(VAR_IF_TYPE(SysShape,SYS_ELECTRO)){
130  Pm[NEdge].Typ = 2;
131  Pm[NEdge+1].Typ = 2;
132  }
133  //Pc->PrintCells();
134  //for(int p=0;p<pNPart();p++) CheckDomDec(p);
135  //Interp();
136 }
138  Shout("Freeing");
139  fclose(StatFile1);
140  fclose(StatFile2);
141  Shout("Freeing/md: forces domain decomposition");
142  if(VAR_IF_TYPE(SysAlloc,ALL_FORCES)){
143  delete [] Fm;
144  }
145  if(VAR_IF_TYPE(SysAlloc,ALL_MD)){
146  delete Pc;
147  }
148  if(VAR_IF_TYPE(SysAlloc,ALL_MC)){
149  delete Pc;
150  delete [] OldNrgBead;
151  delete [] OldNrgCh;
152  for(int p=0;p<pNPCh();p++){
153  delete [] OldPos[p];
154  }
155  delete [] OldPos;
156  }
157  if(VAR_IF_TYPE(SysAlloc,ALL_BIAS)){
158  for(int t=0;t<NTrialBias;t++){
159  delete [] BondPosBias[t];
160  }
161  delete [] BondPosBias;
162  delete [] CumProbBias;
163  }
164  if(VAR_IF_TYPE(SysAlloc,ALL_SPLINE)){
165  free(Pl);
166  }
167  if(VAR_IF_TYPE(SysAlloc,ALL_DENS)){
168  free(LocDens2);
169  free(LocDens3);
170  free(Dens2);
171  free(Dens3);
172  }
173  if(VAR_IF_TYPE(SysShape,SYS_ROD)){
174  delete IntMatrix;
175  }
176  if(VAR_IF_TYPE(SysShape,SYS_2D)){
177  delete IntMatrix;
178  }
179  if(VAR_IF_TYPE(SysShape,SYS_LEAVES)){
180  delete IntMatrix;
181  }
182 #ifdef __glut_h__
183  free(Cylinder);
184 #endif
185 }
186 Forces::Forces(int argc,char **argv,char *ConfF,char *Snapshot){
187  Shout("constructor with intial configuration");
188  InitConst();
189  if(ReadConfDinamica(ConfF)){
190  Dx = 1./(double)(NEdge-1);
191  }
192  Open(Snapshot,BF_PART);
193  if(VAR_IF_TYPE(SysShape,SYS_LEAVES)) IfNano = 2;
194  if(VAR_IF_TYPE(SysShape,SYS_PORE)) IfNano = 2;
195  AllocMethod();
196  PrepareParallel(argc,argv);
197  PrepareSys();
198  //Interp();
199 }
201  Shout("Preparing system");
202  OldNrgSys = 0.;
204  /*calculate the energies per chain or per particle for
205  the calculation of the MC.
206  */
207  if(VAR_IF_TYPE(SysShape,SYS_MC)){
208  if(VAR_IF_TYPE(SysAlloc,ALL_DENS)){
209  Shout("Preparing system/adding densities\n");
210  OldNrgSys = DensFuncNrgSys();
211  ChemPotEx += NrgPBead;
212  if(VAR_IF_TYPE(CalcMode,CALC_NcVT)){
213  //CalcTotNrgCh();
214  }
215  if(VAR_IF_TYPE(CalcMode,CALC_NVT)){
216  ;//CalcTotNrgBead();
217  }
218  if(VAR_IF_TYPE(CalcMode,CALC_mcVT)){
219  if(!VAR_IF_TYPE(CalcMode,CALC_CONF_BIAS))
220  ;//CalcTotNrgCh();
221  }
222  }
223  //defining all the allocated chains
224  for(int p=0;p<pNAllocP();p++){
225  int c = (int)(p/(double)pNPCh());
226  Pm[p].CId = c;
227  Pm[p].Typ = 0;
228  if( p%pNPCh() >= Block[0].Asym ) Pm[p].Typ = 1;
229  if( p%pNPCh() == pNPCh() - 1) continue;
230  Ln[p].NLink = 1;
231  Ln[p].Link[0] = p+1;
232  }
233  }
234  if(VAR_IF_TYPE(SysShape,SYS_ELECTRO)){
235  for(int p=0;p<NEdge;p++){
236  Pm[p].Typ = 2;
237  Pm[p].Idx = p;
238  Pm[p].CId = p;
239  Ln[p].NLink = 0;
240  SetBkf(p);
241  }
242  for(int p=NEdge;p<NEdge+NSpline;p++){
243  Pm[p].Typ = 0;
244  Pm[p].Idx = p;
245  Pm[p].CId = NEdge;
246  Ln[p].NLink = 1;
247  Ln[p].Link[0] = p+1;
248  if(p==NEdge+NSpline-1) Ln[p].NLink = 0;
249  SetBkf(p);
250  }
251  }
252  if(VAR_IF_TYPE(SysShape,SYS_MD)){
253  Shout("Preparing system/calculating forces\n");
254  OldNrgSys = SumForcesMD();
255  }
256  else if(VAR_IF_TYPE(SysShape,SYS_ROD)){
257  IntMatrix = new Matrice(pNPart(),pNPart());
258  }
259  else if(VAR_IF_TYPE(SysShape,SYS_LEAVES)){
260  IntMatrix = new Matrice(pNPCh(),pNPCh());
261  }
262  else if(VAR_IF_TYPE(SysShape,SYS_2D)){
263  IntMatrix = new Matrice(pNPart(),pNPart());
264  }
265  /*for the chemical potential, to avoid the calculation of
266  large exponential */
267  NrgPBead = 0.;//2.*OldNrgSys/(double)pNPart();
268  Shout("Prepared");
269 }
271  if(!IfFillMatrix) return;
272  IntMatrix->Clear();
273  double Inter = pEdge(CLat1)/(double)NEdge;//fabs(pPos(1,0) - pPos(0,0));
274  SPLINE Weight;
275  Weight.a0 = Kf.El[2];
276  Weight.a1 = 0./Inter;
277  Weight.a2 = Kf.Lap/SQR(Inter);
278  Weight.a3 = 0./(Inter*SQR(Inter));
279  Weight.a4 = Kf.SLap/(SQR(Inter)*SQR(Inter));
280  if(VAR_IF_TYPE(SysShape,SYS_ROD)){
281  int NDim = 1;
282  Matrice *CoeffMatrix = new Matrice(Weight,NDim);
283  CoeffMatrix->Print();
284  for(int r=0;r<pNPart();r++){
285  if(Pm[r].Typ != 0){ IntMatrix->Set(r,r,1.);continue;}
286  if(r >= 2) IntMatrix->Set(r,r-2,CoeffMatrix->Val(2,0));
287  if(r >= 1) IntMatrix->Set(r,r-1,CoeffMatrix->Val(2,1));
288  if(r < pNPart()-1) IntMatrix->Set(r,r+1,CoeffMatrix->Val(2,3));
289  if(r < pNPart()-2) IntMatrix->Set(r,r+2,CoeffMatrix->Val(2,4));
290  IntMatrix->Set(r,r,CoeffMatrix->Val(2,2));
291  }
292  IntMatrix->Invert();
293  //IntMatrix->Print();
294  delete CoeffMatrix;
295  }
296  else if(VAR_IF_TYPE(SysShape,SYS_LEAVES)){
297  int NDim = 1;
298  Matrice *CoeffMatrix = new Matrice(Weight,NDim);
299  CoeffMatrix->Print();
300  for(int r=0;r<pNPCh();r++){
301  if(Pm[r].Typ != 0){ IntMatrix->Set(r,r,1.);continue;}
302  if(r >= 2) IntMatrix->Set(r,r-2,CoeffMatrix->Val(2,0));
303  if(r >= 1) IntMatrix->Set(r,r-1,CoeffMatrix->Val(2,1));
304  if(r < pNPart()-1) IntMatrix->Set(r,r+1,CoeffMatrix->Val(2,3));
305  if(r < pNPart()-2) IntMatrix->Set(r,r+2,CoeffMatrix->Val(2,4));
306  IntMatrix->Set(r,r,CoeffMatrix->Val(2,2));
307  }
308  //IntMatrix->Invert();
309  IntMatrix->Print();
310  delete CoeffMatrix;
311  }
312  else if(VAR_IF_TYPE(SysShape,SYS_2D)){
313  int NDim = 2;
314  Matrice *CoeffMatrix = new Matrice(Weight,NDim);
315  CoeffMatrix->Print();
316  for(int p=0;p<pNPart();p++){
317  if(Pm[p].Typ != 0){
318  IntMatrix->Set(p,p,1.);
319  continue;
320  }
321  int pym1 = Ln[p].Link[0];
322  int pyp1 = Ln[p].Link[1];
323  int pym2 = Ln[pym1].Link[0];
324  int pyp2 = Ln[pyp1].Link[1];
325  int pxm1 = Ln[p].Link[2];
326  int pxp1 = Ln[p].Link[3];
327  int pxm2 = Ln[pxm1].Link[2];
328  int pxp2 = Ln[pxp1].Link[3];
329  // printf("%d)\n",p);
330  //printf("%d %d %d %d\n",pym2,pym1,pyp1,pyp2);
331  // printf("%d %d %d %d\n",pxm2,pxm1,pxp1,pxp2);
332  if(PeriodicImage[0]){
333  IntMatrix->Set(p,pxm2,CoeffMatrix->Val(2,0));
334  IntMatrix->Set(p,pxm1,CoeffMatrix->Val(2,1));
335  IntMatrix->Set(p,pxp1,CoeffMatrix->Val(2,3));
336  IntMatrix->Set(p,pxp2,CoeffMatrix->Val(2,4));
337  }
338  else{
339  if(pxm2 == p-2*nEdge[1])
340  IntMatrix->Set(p,pxm2,CoeffMatrix->Val(2,0));
341  if(pxm1 == p-nEdge[1])
342  IntMatrix->Set(p,pxm1,CoeffMatrix->Val(2,1));
343  if(pxp1 == p+nEdge[1])
344  IntMatrix->Set(p,pxp1,CoeffMatrix->Val(2,3));
345  if(pxp2 == p+2*nEdge[1])
346  IntMatrix->Set(p,pxp2,CoeffMatrix->Val(2,4));
347  }
348  IntMatrix->Add(p,p,CoeffMatrix->Val(2,2));
349  if(PeriodicImage[1]){
350  IntMatrix->Set(p,pym2,CoeffMatrix->Val(2,0));
351  IntMatrix->Set(p,pym1,CoeffMatrix->Val(2,1));
352  IntMatrix->Set(p,pyp1,CoeffMatrix->Val(2,3));
353  IntMatrix->Set(p,pyp2,CoeffMatrix->Val(2,4));
354  }
355  else{
356  if(pym2 == p-2)
357  IntMatrix->Set(p,pym2,CoeffMatrix->Val(2,0));
358  if(pym1 == p-1)
359  IntMatrix->Set(p,pym1,CoeffMatrix->Val(2,1));
360  if(pyp1 == p+1)
361  IntMatrix->Set(p,pyp1,CoeffMatrix->Val(2,3));
362  if(pyp2 == p+2)
363  IntMatrix->Set(p,pyp2,CoeffMatrix->Val(2,4));
364  }
365  }
366  IntMatrix->Invert();
367  //IntMatrix->Print();
368  delete CoeffMatrix;
369  }
370  IfFillMatrix = 0;
371 }
372 void Forces::PrepareParallel(int argc,char **argv){
373 #ifdef OMPI_MPI_H
374  Shout("Preparing parallelisation");
375  MPI_Init(&argc,&argv);
376  int Rank=0,Size=0;
377  MPI_Comm_rank(MPI_COMM_WORLD, &Rank);
378  MPI_Comm_size(MPI_COMM_WORLD, &Size);
379  int Partition = (int)(argc/(double)Size);
380  NFile[0] = Partition*Rank;
381  NFile[1] = Partition*(Rank+1);
382  if(Rank==Size-1) NFile[1] += argc%Size;
383  Proc = new SingProc(Size,Rank);
384 #endif
385 }
387  Shout("Allocation");
388  StatFile1 = fopen("StatDyn1.dat","w");
389  StatFile2 = fopen("StatDyn2.dat","w");
390  ChooseCalcMode(CalcMode);
391  ChoosePot(CalcMode);
393  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
394  double Pos[3];
395  //Pc->PrintCells();
396  //CheckPairList();
397  SetDeltat(Deltat);
398  if(VAR_IF_TYPE(SysShape,SYS_MD)){
399  Shout("Allocating/md: forces, domain decomposition");
400  Fm = (FORCES *)calloc(pNAllocP(),sizeof(FORCES));
401  VAR_ADD_TYPE(SysAlloc,ALL_MD);
402  VAR_ADD_TYPE(SysAlloc,ALL_FORCES);
403  Pc = new DomDec(Edge,pNPart(),sqrt(Kf.CutOff2));
404  for(int p=0;p<pNPart();p++){pPos(p,Pos);Pc->AddPart(p,Pos);}
405  }
406  if(VAR_IF_TYPE(SysShape,SYS_ROD)){
407  Shout("Allocating/rod: forces");
408  Fm = (FORCES *)calloc(pNAllocP(),sizeof(FORCES));
409  VAR_ADD_TYPE(SysAlloc,ALL_FORCES);
410  }
411  else if(VAR_IF_TYPE(SysShape,SYS_MC)){
412  Shout("Allocating/mc: domain decomposition old chain positions, old energies for particles and chains, first bead distribution, bias (cumulative probabilities and bead positions)");
413  Pc = new DomDec(Edge,pNPart(),sqrt(Kf.CutOff2));
414  double Pos[3];
415  for(int p=0;p<pNPart();p++){pPos(p,Pos);Pc->AddPart(p,Pos);}
416  //Pc->PrintCells();
417  if(VAR_IF_TYPE(CalcMode,CALC_NVT)){
418  OldNrgBead = new double[3*pNAllocP()];
419  if(!OldNrgBead){
420  printf("Could not allocate OldNrgPm \n");
421  exit(1);
422  }
423  }
424  OldNrgCh = new double[3*pNAllocC()];
425  FirstBeadDistr = new double[NBin];
426  if(!OldNrgCh){
427  printf("Could not allocate OldNrgCh\n");
428  exit(1);
429  }
430  OldPos = (double **)calloc(pNPCh(),sizeof(double));
431  if( !OldPos){
432  printf("Could not allocate OldPos\n");
433  exit(1);
434  }
435  for(int p=0;p<pNPCh();p++){
436  OldPos[p] = (double *)calloc(3,sizeof(double));
437  }
438  CumProbBias = new double[NTrialBias];
439  BondPosBias = new double*[NTrialBias];
440  for(int t=0;t<NTrialBias;t++){
441  BondPosBias[t] = new double[3];
442  }
443  VAR_ADD_TYPE(SysAlloc,ALL_MC);
444  VAR_ADD_TYPE(SysAlloc,ALL_BIAS);
445  GaussVar = sqrt(pReOverCutOff()/(3.*(pNPCh()-1)))/2.;
446  GaussVar = sqrt(1./pkSpr());
447  }
448  else if(VAR_IF_TYPE(SysShape,SYS_ELECTRO)){
449  Shout("Allocating/mc: domain decomposition old chain positions, old energies for particles and chains, first bead distribution, bias (cumulative probabilities and bead positions)");
450  Pc = new DomDec(Edge,pNPart(),sqrt(Kf.CutOff2));
451  double Pos[3];
452  for(int p=0;p<pNPart();p++){pPos(p,Pos);Pc->AddPart(p,Pos);}
453  //Pc->PrintCells();
454  if(VAR_IF_TYPE(CalcMode,CALC_NVT)){
455  OldNrgBead = new double[3*pNAllocP()];
456  if(!OldNrgBead){
457  printf("Could not allocate OldNrgPm \n");
458  exit(1);
459  }
460  }
461  OldNrgCh = new double[3*pNAllocC()];
462  FirstBeadDistr = new double[NBin];
463  if(!OldNrgCh){
464  printf("Could not allocate OldNrgCh\n");
465  exit(1);
466  }
467  OldPos = (double **)calloc(pNPCh(),sizeof(double));
468  if( !OldPos){
469  printf("Could not allocate OldPos\n");
470  exit(1);
471  }
472  for(int p=0;p<pNPCh();p++){
473  OldPos[p] = (double *)calloc(3,sizeof(double));
474  }
475  CumProbBias = new double[NTrialBias];
476  BondPosBias = new double*[NTrialBias];
477  for(int t=0;t<NTrialBias;t++){
478  BondPosBias[t] = new double[3];
479  }
480  VAR_ADD_TYPE(SysAlloc,ALL_MC);
481  VAR_ADD_TYPE(SysAlloc,ALL_BIAS);
482  }
483  if(VAR_IF_TYPE(SysShape,SYS_1D) || VAR_IF_TYPE(SysShape,SYS_TRIAL)){
484  Shout("Allocating/splines");
485  Pl = (PART *)calloc(NSpline,sizeof(PART));
486  VAR_ADD_TYPE(SysAlloc,ALL_SPLINE);
487  }
488  if(VAR_IF_TYPE(CalcMode,CALC_DENS)){
489  Shout("Allocating/particle densities, sum of local densities");
490  LocDens2 = (double *)calloc(pNPCh()*pNType(),sizeof(double));
491  LocDens3 = (double *)calloc(pNPCh()*pNType(),sizeof(double));
492  Dens2 = (double *)calloc(pNAllocP()*pNType(),sizeof(double));
493  Dens3 = (double *)calloc(pNAllocP()*pNType(),sizeof(double));
494  if( !Dens2 || !Dens3){
495  printf("Could not allocate Dens2 Dens3\n");
496  exit(1);
497  }
498  VAR_ADD_TYPE(SysAlloc,ALL_DENS);
499  }
500 #ifdef __glut_h__
501  Cylinder = (GLuint *)calloc(pNNano(),sizeof(GLuint));
502 #endif
503  ChemPotId = log(NChemPotId/pVol());
504 }
506  Shout("Init constants");
507  NEdge = 80+1;//NInEdge;
508  Dx = 1./(double)(NEdge-1);
509  Kf.El[0] = 11.;//2.*NEdge/10.;//Elastic coupling
510  Kf.El[1] = 11.;//Elastic coupling
511  Kf.El[2] = 11.;//Elastic coupling
512  Kf.Lap = 20.*pow(Dx,2.);//bending rigidity
513  Kf.SLap = 0.;//.015*pow(Dx,2.);//Surface tension
514  Kf.Ext = 1.;
515  Kf.LJ = 1.;
516  Kf.LJMin = .5;
517  Kf.CutOff2 = 1.;
518  Kf.Cont = 0.;//.001;
519  Kf.Elong[0] = 0.1;
520  Kf.Elong[1] = 0.1;
521  Kf.Elong[2] = 0.1;
522  Kf.ForThr = 100.;
523  IntMax = 100;
524  IncrDist = 0.01;
525  IfInterp=FIT_FORTH;
526  IfFillMatrix=1;
527  Nano->Rad = .05;
528  Nano->Height = .3;
529  ChemPotId = 0.;
530  ChemPotEx = 0.;
531  Nano->Pos[0] = 0.;
532  Nano->Pos[1] = 0.;
533  Nano->Pos[2] = .5;
534  BoundCond[0] = 1;BoundCond[1] = 1;
535  BoundCond[2] = 1;BoundCond[3] = 1;
536  BoundCond[4] = 1;BoundCond[5] = 1;
537  PeriodicImage[0] = 1;PeriodicImage[1] = 1;PeriodicImage[2] = 1;
538  Time = 0.;
539  NUpdate = 100;
540  NWrite = 1000;
541  NBin = 100;
542  NTrialBias = 10;
543  // GaussVar = sqrt(1./pkSpr());
544  ThermMode = THERM_LANG;
545  CalcMode = CALC_LJ39;
546  // Part2Move = NEdge/2+NEdge;
547 #ifdef __glut_h__
548  NSpline = 100;
549  IfSphere = 0;
550  IfLine = 1;
551  IfSpline = 0;
552  NShow = 1;
553  IfMovie=0;
554  Frame = 0;
555  IfExt = 0;
556  IfRot = 0;
557  BeadType = 0;
558 #endif
559  IfExit = 0;
560  IfNano = 0;
561  DynFlag = 0;
562  SetNBlock(1);
563 }
565  printf("-----INFO-------\n");
566  if(VAR_IF_TYPE(SysShape,SYS_MC)){
567  // Pc->Erase();
568  // ClearDens();
569  // for(int p=0;p<pNPart();p++)
570  // Pc->AddPart(p,Pm[p].Pos);
571  // AddDens(0,pNPart());
572  printf("Energy: Sys %lf Ch Part \n",OldNrgSys);
573  // CalcTotNrgCh();
574  // for(int c=0;c<pNChain();c++)
575  // printf("%d %lf %lf %lf\n",c,OldNrgCh[c*3],OldNrgCh[c*3+1],OldNrgCh[c*3+2]);
576  }
577  else if(VAR_IF_TYPE(SysShape,SYS_MD)){
578  for(int p=0;p<pNPart();p++){
579  if(p == Bead2Move) printf("p ");
580  printf("%d) (%lf,%lf,%lf) (%lf %lf %lf) %lf-%lf-%lf\n",p,pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0),pVel(p,1),pVel(p,2),Fm[p].Dir[2],Fm[p].Dir[2],Fm[p].Dir[2]);
581  }
582  }
583  printf("------------\n");
584 }
585 int Forces::ReadConfDinamica(char *InFile){
586  SysType = 0;
587  SysShape = 0;
588  CalcMode = 0;
589  FILE *FileToRead = fopen(InFile,"r");
590  if(FileToRead == NULL){
591  printf("The conf file %s is missing\n",InFile);
592  return 1;
593  }
594  double buff[12];
595  char SysInit[20];
596  char *Line = (char *)malloc(256*sizeof(char));
597  //fgets(Line,256,FileToRead);
598  int NNano = 0;
599  for(int k=0;!(fgets(Line,256,FileToRead)==NULL);k++){
600  if(strstr(Line, "Rigid") == Line) NNano++;
601  }
602  SetNNano(NNano);
603  NNano = 0;
604  rewind(FileToRead);
605  for(int k=0;!(fgets(Line,256,FileToRead)==NULL);k++){
606  //printf("%s",Line);
607  if(1 == sscanf(Line,"NEdge %lf",buff) )
608  NEdge = (int)*buff;
609  else if(1 == sscanf(Line,"SysShape %s",SysInit) ){
610  if(!strcmp(SysInit,"leaves") )
611  VAR_ADD_TYPE(SysShape,SYS_LEAVES);
612  else if(!strcmp(SysInit,"pore") )
613  VAR_ADD_TYPE(SysShape,SYS_PORE);
614  else if(!strcmp(SysInit,"1d") )
615  VAR_ADD_TYPE(SysShape,SYS_1D);
616  else if(!strcmp(SysInit,"2d") )
617  VAR_ADD_TYPE(SysShape,SYS_2D);
618  else if(!strcmp(SysInit,"3d") )
619  VAR_ADD_TYPE(SysShape,SYS_3D);
620  else if(!strcmp(SysInit,"rod") )
621  VAR_ADD_TYPE(SysShape,SYS_ROD);
622  else if(!strcmp(SysInit,"trial") )
623  VAR_ADD_TYPE(SysShape,SYS_TRIAL);
624  else if(!strcmp(SysInit,"rigid") )
625  VAR_ADD_TYPE(SysShape,SYS_RIGID);
626  else if(!strcmp(SysInit,"stalk") )
627  VAR_ADD_TYPE(SysShape,SYS_STALK);
628  else if(!strcmp(SysInit,"md") )
629  VAR_ADD_TYPE(SysShape,SYS_MD);
630  else if(!strcmp(SysInit,"mc") )
631  VAR_ADD_TYPE(SysShape,SYS_MC);
632  else if(!strcmp(SysInit,"electro") ){
633  VAR_ADD_TYPE(SysShape,SYS_ELECTRO);
634  //VAR_ADD_TYPE(SysShape,SYS_MC);
635  }
636  else{
637  printf("system type not recognized\n");
638  exit(1);
639  }
640  }
641  else if(1 == sscanf(Line,"CalcMode %s",SysInit) ){
642  if(!strcmp(SysInit,"NVT") )
643  VAR_ADD_TYPE(CalcMode,CALC_NVT);
644  else if(!strcmp(SysInit,"NcVT") )
645  VAR_ADD_TYPE(CalcMode,CALC_NcVT);
646  else if(!strcmp(SysInit,"mcVT") )
647  VAR_ADD_TYPE(CalcMode,CALC_mcVT);
648  else if(!strcmp(SysInit,"mVT") )
649  VAR_ADD_TYPE(CalcMode,CALC_mVT);
650  else{
651  printf("calculation type not recognized\n");
652  exit(1);
653  }
654  }
655  else if(1 == sscanf(Line,"Thermostat %s",SysInit) ){
656  if(!strcmp(SysInit,"Langevin") )
657  ThermMode = THERM_LANG;
658  else if(!strcmp(SysInit,"Andersen") )
659  ThermMode = THERM_AND;
660  else if(!strcmp(SysInit,"Berendsen") )
661  ThermMode = THERM_BERE;
662  else if(!strcmp(SysInit,"no") )
663  ThermMode = THERM_NO;
664  else{
665  printf("thermostat not recognized\n");
666  exit(1);
667  }
668  }
669  else if(1 == sscanf(Line,"PotentialMode %s",SysInit) ){
670  if(!strcmp(SysInit,"Pair") )
671  VAR_ADD_TYPE(CalcMode,CALC_PAIR);
672  else if(!strcmp(SysInit,"DensFunc") )
673  VAR_ADD_TYPE(CalcMode,CALC_DENS);
674  else if(!strcmp(SysInit,"DensFuncCh") )
675  VAR_ADD_TYPE(CalcMode,CALC_DENS);
676  }
677  else if(1 == sscanf(Line,"Potential %s",SysInit) ){
678  if(!strcmp(SysInit,"LJ") )
679  VAR_ADD_TYPE(CalcMode,CALC_LJ);
680  else if(!strcmp(SysInit,"LJ39") )
681  VAR_ADD_TYPE(CalcMode,CALC_LJ39);
682  else if(!strcmp(SysInit,"Harmonic") )
683  VAR_ADD_TYPE(CalcMode,CALC_HARM);
684  else if(!strcmp(SysInit,"Step") )
685  VAR_ADD_TYPE(CalcMode,CALC_STEP);
686  else if(!strcmp(SysInit,"Electro") )
687  VAR_ADD_TYPE(CalcMode,CALC_ELECTRO);
688  else{
689  printf("interaction potential not recognized\n");
690  exit(1);
691  }
692  }
693  else if(1 == sscanf(Line,"IfInterp %lf",buff) )
694  IfInterp = (int)*buff;
695  else if(1 == sscanf(Line,"Lap %lf",buff) )
696  Kf.Lap = *buff;
697  else if(1 == sscanf(Line,"SLap %lf",buff) )
698  Kf.SLap = *buff;
699  else if(1 == sscanf(Line,"Ext %lf",buff) )
700  Kf.Ext = *buff;
701  else if(1 == sscanf(Line,"LJ %lf",buff) )
702  Kf.LJ = *buff;
703  else if(1 == sscanf(Line,"LJMin %lf",buff) )
704  Kf.LJMin = *buff;
705  else if(1 == sscanf(Line,"kSpr %lf",buff) )
706  SetkSpr(*buff);
707  else if(1 == sscanf(Line,"SprRest %lf",buff) )
708  SetSprRest(*buff);
709  else if(1 == sscanf(Line,"kBen %lf",buff) )
710  SetkBen(*buff);
711  else if(1 == sscanf(Line,"SimLimit %lf",buff) )
712  SimLimit = (int)*buff;
713  else if(1 == sscanf(Line,"CutOff %lf",buff) )
714  Kf.CutOff2 = SQR(*buff);
715  else if(1 == sscanf(Line,"Cont %lf",buff) )
716  Kf.Cont = *buff;
717  else if(1 == sscanf(Line,"NWrite %lf",buff) )
718  NWrite = (int)*buff;
719  else if(1 == sscanf(Line,"NBin %lf",buff) )
720  NBin = (int)*buff;
721  else if(1 == sscanf(Line,"NGrid %lf",buff) )
722  NGrid = (int)*buff;
723  else if(1 == sscanf(Line,"NUpdate %lf",buff) )
724  NUpdate = (int)*buff;
725  else if(strstr(Line, "Rigid") == Line){
726  NanoString(Line,NNano++);
727  }
728 #ifdef __glut_h__
729  else if(1 == sscanf(Line,"IfMovie %lf",buff) )
730  IfMovie = (int)*buff;
731  else if(1 == sscanf(Line,"IfLine %lf",buff) )
732  IfLine = (int)*buff;
733  else if(1 == sscanf(Line,"NSpline %lf",buff) )
734  NSpline = *buff;
735  else if(1 == sscanf(Line,"IfSphere %lf",buff) )
736  IfSphere = (int)*buff;
737 #endif
738  else if(3 == sscanf(Line,"Edge %lf %lf %lf",buff,buff+1,buff+2) ){
739  SetEdge(buff[0],0);
740  SetEdge(buff[1],1);
741  SetEdge(buff[2],2);
742  }
743  else if(3 == sscanf(Line,"El %lf %lf %lf",buff,buff+1,buff+2) ){
744  Kf.El[0] = *buff;
745  Kf.El[1] = *(buff+1);
746  Kf.El[2] = *(buff+2);
747  }
748  else if(3 == sscanf(Line,"Elong %lf %lf %lf",buff,buff+1,buff+2) ){
749  Kf.Elong[0] = *buff;
750  Kf.Elong[1] = *(buff+1);
751  Kf.Elong[2] = *(buff+2);
752  }
753  else if(6 == sscanf(Line,"Boundary %lf %lf %lf %lf %lf %lf",buff,buff+1,buff+2,buff+3,buff+4,buff+5)){
754  BoundCond[0] = (int)buff[0];
755  BoundCond[1] = (int)buff[1];
756  BoundCond[2] = (int)buff[2];
757  BoundCond[3] = (int)buff[3];
758  BoundCond[4] = (int)buff[4];
759  BoundCond[5] = (int)buff[5];
760  }
761  else if(3 == sscanf(Line,"Periodic %lf %lf %lf",buff,buff+1,buff+2)){
762  PeriodicImage[0] = (int)*buff;
763  PeriodicImage[1] = (int)*(buff+1);
764  PeriodicImage[2] = (int)*(buff+2);
765  }
766  else if(1 == sscanf(Line,"Deltat %lf",buff) )
767  Deltat = *buff;
768  else if(1 == sscanf(Line,"Temp %lf",buff) )
769  SetTemp(buff[0]);
770  else if(1 == sscanf(Line,"NChemPotId %lf",buff) )
771  NChemPotId = *buff;
772  else if(1 == sscanf(Line,"NTrialBias %lf",buff) ){
773  NTrialBias =(int) *buff;
774  }
775  else if(1 == sscanf(Line,"ChemPotEx %lf",buff) )
776  ChemPotEx = *buff;
777  else if(1 == sscanf(Line,"IfConfBias %lf",buff) ){
778  if((int)*buff == 1)
779  VAR_ADD_TYPE(CalcMode,CALC_CONF_BIAS);
780  }
781  else if(1 == sscanf(Line,"IfBilBias %lf",buff) ){
782  if((int)*buff == 1){
783  VAR_ADD_TYPE(CalcMode,CALC_BIL_BIAS);
784  StudySys();
785  }
786  }
787  else if(1 == sscanf(Line,"IfSphBias %lf",buff) ){
788  if((int)*buff == 1)
789  VAR_ADD_TYPE(CalcMode,CALC_SPH_BIAS);
790  }
791  else if(1 == sscanf(Line,"Viscosity %lf",buff) )
792  Viscosity = *buff;
793  else if(1 == sscanf(Line,"TNSlab %lf",buff) )
794  Tens.NSlab = (int)*buff;
795  else if(1 == sscanf(Line,"TNComp %lf",buff) )
796  Tens.NComp = (int)*buff;
797  else if(1 == sscanf(Line,"TCalcMode %s",SysInit) ){
798  if(!strcmp(SysInit,"2d") ){
799  VAR_ADD_TYPE(Tens.CalcMode,CALC_2d);
800  Tens.NDim = 2;
801  }
802  else if(!strcmp(SysInit,"3d") ){
803  VAR_ADD_TYPE(Tens.CalcMode,CALC_3d);
804  Tens.NDim = 3;
805  }
806  else{
807  printf("Pressure summation not recognized\n");
808  exit(1);
809  }
810  }
811  }
812  //for(int n=0;n<pNNano();n++) for(int d=0;d<3;d++) Nano[n].Pos[d] *= pEdge(d);
813  ChemPotId = log(NChemPotId/pVol());
814  printf("Sys] NEdge %d SysShape %d %s Interp %d \n",NEdge,SysShape,SysInit,IfInterp);
815  printf("Forces] Lap %lf SLap %lf Ext %lf LJ %lf Cont %lf El %lf %lf %lf Elong %lf %lf %lf\n",Kf.Lap,Kf.SLap,Kf.Ext,Kf.LJ,Kf.Cont,Kf.El[0],Kf.El[1],Kf.El[2],Kf.Elong[0],Kf.Elong[1],Kf.Elong[2]);
816  printf("External] Center %lf %lf %lf Rad %lf Hei %lf\n",Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],Nano->Rad,Nano->Height);
817  fclose(FileToRead);
818  return 0;
819 }
820 int Forces::ReSetNPart(int NewNPart){
821  int OldNPart = pNPart();
822  //if the number of allocated particle is less then the new particles returns
823  Block[0].NPart = NewNPart;
824  if(SetNPart(NewNPart)) return 1;
825  if(VAR_IF_TYPE(SysAlloc,ALL_MC)){
826  double *Dens2T = (double *)calloc(pNPart()*pNType(),sizeof(double));
827  double *Dens3T = (double *)calloc(pNPart()*pNType(),sizeof(double));
828  for(int p=0;p<OldNPart*pNType();p++){
829  Dens2T[p] = Dens2[p];
830  Dens3T[p] = Dens3[p];
831  }
832  //memcpy(Dens2,Dens2T,OldNPart*pNType());
833  //memcpy(Dens3,Dens3T,OldNPart*pNType());
834  free(Dens2);
835  free(Dens3);
836  Dens2 = (double *)calloc(pNAllocP()*pNType(),sizeof(double));
837  Dens3 = (double *)calloc(pNAllocP()*pNType(),sizeof(double));
838  // memcpy(Dens2T,Dens2,OldNPart*pNType());
839  // memcpy(Dens3T,Dens3,OldNPart*pNType());
840  for(int p=0;p<OldNPart*pNType();p++){
841  Dens2[p] = Dens2T[p];
842  Dens3[p] = Dens3T[p];
843  }
844  free(Dens2T);
845  free(Dens3T);
846  double *TempNrgPart = new double[OldNPart];
847  for(int p=0;p<OldNPart;p++){
848  TempNrgPart[p] = OldNrgBead[p];
849  }
850  delete[] OldNrgBead;
851  OldNrgBead = new double[pNAllocP()];
852  for(int p=0;p<OldNPart;p++){
853  OldNrgBead[p] = TempNrgPart[p];
854  }
855  delete []TempNrgPart;
856  }
857  if(VAR_IF_TYPE(SysAlloc,ALL_MD)){
858  FORCES *TmpFm = new FORCES [OldNPart];
859  for(int p=0;p<OldNPart;p++){
860  for(int d=0;d<3;d++){
861  TmpFm[p].Dir[d] = Fm[p].Dir[d];
862  TmpFm[p].Ext[d] = Fm[p].Ext[d];
863  }
864  }
865  free(Fm);
866  Fm = (FORCES *)calloc(pNAllocP(),sizeof(FORCES));
867  for(int p=0;p<OldNPart;p++){
868  for(int d=0;d<3;d++){
869  Fm[p].Dir[d] = TmpFm[p].Dir[d];
870  Fm[p].Ext[d] = TmpFm[p].Ext[d];
871  }
872  }
873  free(TmpFm);
874  }
875  return 0;
876 }
877 int Forces::ReSetNChain(int NewNChain){
878  int OldNChain = pNChain();
879  Block[0].NChain = NewNChain;
880  if(SetNChain(NewNChain))return 1;
881  if(VAR_IF_TYPE(SysAlloc,ALL_MC)){
882  double *TempNrgCh = new double[OldNChain];
883  for(int c=0;c<3*OldNChain;c++){
884  TempNrgCh[c] = OldNrgCh[c];
885  }
886  delete[] OldNrgCh;
887  OldNrgCh = new double[pNAllocC()];
888  for(int c=0;c<3*OldNChain;c++){
889  OldNrgCh[c] = TempNrgCh[c];
890  }
891  delete []TempNrgCh;
892  }
893  if(VAR_IF_TYPE(SysAlloc,ALL_MD)){
894  }
895  return 0;
896 }
897 void Forces::ReSetNPCh(int NewNPCh){
898  Block[0].NPCh = NewNPCh;
899  SetNPCh(NewNPCh);
900 }
901 void Forces::ReOpen(char *FName,int Bf){
902  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
903  Pc->Erase();
904  delete Pc;
905  Open(FName,Bf);
906  Pc = new DomDec(Edge,pNPart(),sqrt(Kf.CutOff2));
907  double Pos[3];
908  for(int p=0;p<pNPart();p++){pPos(p,Pos);Pc->AddPart(p,Pos);}
909  PrepareSys();
910 }
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
int SetNLink(int NewNCh)
Set and reallocate the number of links.
Definition: VarDataComm.cpp:41
double NChemPotId
Equilibrium number of particles/chains.
Definition: Forces.h:719
void SetNPCh(int NewNCh)
Set and reallocate the number of particles per chains.
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
int CId
Chain Identifier.
Definition: VarData.h:224
double LJMin
Lennard-Jones.
Definition: Forces.h:161
double Ext
External force.
Definition: Forces.h:153
void ReSetNPCh(int NewNPCh)
Set a new number of particle per chain.
Definition: Forces.cpp:897
double Cont
Contact/Friction.
Definition: Forces.h:157
double pVel(int p, int d)
Return the velocity.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void Erase()
Erase the pairlist.
Definition: Cubo.cpp:258
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
NANO * Nano
Extra particle.
Definition: VarData.h:1044
void PrepareSys()
Calculate some initial quantities for the succesive calculations.
Definition: Forces.cpp:200
int IntMax
Boh.
Definition: Forces.h:741
Energy contribution from the density functional hamiltonian.
Definition: VarData.h:388
int SetNChain(int NewNCh)
Set and reallocate the number of chains.
Definition: VarDataComm.cpp:78
double Lap
Prefactor of the Laplacian.
Definition: Forces.h:147
int NDim
of dimensions
Definition: Forces.h:193
void SetCoeff(double Co, int t1, int t2)
Set the prefactor.
Definition: VarData.cpp:252
int pNAllocC()
Allocated number of chains.
int Asym
Diblock limit of the chain.
Definition: VarData.h:273
int SetNBlock(int Val)
Set NBlock.
double Rad
Size.
Definition: VarData.h:445
double Viscosity
Viscosity of the medium.
Definition: Forces.h:731
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
double GaussVar
Standard deviation of the gaussian chain.
Definition: Forces.h:725
double Height
Height of the cylinder.
Definition: VarData.h:449
double El[3]
Elastic force.
Definition: Forces.h:151
int IfInterp
If interpolates with the splines.
Definition: Forces.h:755
void FillMatrix()
Fill the entries of the interaction matrix.
Definition: Forces.cpp:270
double SumForcesMD()
Iterate all over the particles and calculate the forces.
int Frame
Current number of frame.
Definition: Forces.h:823
double Dir[3]
Direction.
Definition: Forces.h:178
int NWrite
How many timesteps before write the snapshot.
Definition: Forces.h:775
double Deltat
Time step.
Definition: Forces.h:717
double pkSpr()
Spring coupling.
Definition: VarData.h:934
Forces(int argc, char **argv, int NPart, char *ConfFile)
Create an initial system and choose the simulation thecnique.
Definition: Forces.cpp:15
int IfExt
Boh.
Definition: Forces.h:825
~Forces()
Frees the memory.
Definition: Forces.cpp:137
int IfExit
Exit from the loop.
Definition: Forces.h:765
void SetBkf(int p)
Set the back folded array for the particle p.
int NShow
Boh.
Definition: Forces.h:819
bool Open(char *InFile, int BF)
Open the.
Definition: VarData.cpp:105
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 Invert()
Invert.
int IfFillMatrix
If the matrix has to be changed.
Definition: Forces.h:763
double pVol()
xyzr edges of the simulation box
Definition: VarData.h:922
void PrepareParallel(int argc, char **argv)
Create the grid for the parallelisation.
Definition: Forces.cpp:372
double ChemPotEx
Chemical potential of the particles.
Definition: Forces.h:723
void CreateInitial()
Create an initial configuration and an appropriate force field.
Definition: ForcesCreate.cpp:2
int SimLimit
Maximum number of time steps.
Definition: Forces.h:751
int CalcMode
Calculation mode.
Definition: Forces.h:747
double Pos[3]
Position.
Definition: VarData.h:427
void ChooseThermostat(int Mode)
Choose a calculation mode.
void Shout(const char *s,...)
Internal message.
Definition: Forces.cpp:3
int pNAllocP()
Allocated number of particles.
double Time
Total time.
Definition: Forces.h:733
void SetDeltat(double Val)
Set DeltaT.
Definition: VarData.h:1000
int IfSpline
Boh.
Definition: Forces.h:827
void Info()
System&#39;s info.
Definition: Forces.cpp:564
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 Clear()
Set all the entries to zero.
int pNPCh()
Number of particle per chain.
int Idx
Particle identifier.
Definition: VarData.h:222
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
GLuint * Cylinder
List referring the cylinder.
Definition: Forces.h:815
double Dx
Spatial separation between particles.
Definition: Forces.h:729
int ReSetNPart(int NewNPart)
Realloc the number of particles.
Definition: Forces.cpp:820
void ReOpen(char *FName, int Bf)
Open a new file.
Definition: Forces.cpp:901
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double Val(int row)
Returns a value in 1d.
void InitConst()
Read the config file.
Definition: Forces.cpp:505
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
int NUpdate
How many timesteps before redrawing.
Definition: Forces.h:773
int NPart
particles
Definition: VarData.h:261
void SetkSpr(double Val)
Spring coupling.
Definition: VarData.h:938
double IncrDist
Step to move a particle.
Definition: Forces.h:715
void SetTemp(double Val)
Set Temperature.
Definition: VarData.h:1004
void StudySys()
Obtain informations for a better performance in inserting the chains.
Basics class to start a MPI grid.
Definition: SingProc.h:6
int SysShape
Shape of system.
Definition: Forces.h:745
int SysType
Contains the definition of the system.
Definition: VarData.h:1086
bool Add(int row, int col, double Val)
Add the value of the coefficient to the previous one.
int pNNano()
Number of nanoparticles.
double DensFuncNrgSys()
Calculation of the energy from the density functional Hamiltonian for the system. ...
void AllocMethod()
Allocated the structures needed for the corresponding simulation method.
Definition: Forces.cpp:386
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Single contribution of the forces.
Definition: Forces.h:176
int nEdge[3]
Number of particle per edge.
Definition: Forces.h:743
int NPCh
particles per chain
Definition: VarData.h:259
Matrice computes the algebric operations on matrices.
int CalcMode
Calculation mode.
Definition: Forces.h:197
bool Set(int row, int column, double Val)
Set a coefficient.
double ChemPotId
Chemical potential of the particles.
Definition: Forces.h:721
void SetStep(int Val)
Set Step.
Definition: VarData.h:1002
int IfSphere
Visualize spheres or points.
Definition: Forces.h:831
int NanoString(char *cLine, int n)
Reads and set the specifics of the nano.
int ReSetNChain(int NewNChain)
Realloc the number of chains.
Definition: Forces.cpp:877
double LJ
Lennard-Jones.
Definition: Forces.h:155
int IfLine
Boh.
Definition: Forces.h:837
double pReOverCutOff()
Re/CutOff.
Definition: VarData.h:946
void SetNType(int NewNType)
Set the number of species.
int IfMovie
Produce the images for the video.
Definition: Forces.h:821
int NChain
chains
Definition: VarData.h:263
void SetSprRest(double Val)
Rest distance of the harmonic potential.
Definition: VarData.h:940
int pNType()
of types of the particle
int NComp
of components
Definition: Forces.h:191
double pPos(int p, int d)
Return back folded position.
int pNChain()
Number of chain.
int CLat1
lateral coordinate
Definition: VarData.h:1076
Information of every particle.
Definition: VarData.h:214
double Ext[3]
External.
Definition: Forces.h:182
double Elong[3]
Elongation of the springs.
Definition: Forces.h:159
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Coefficient of a spline.
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
int ReadConfDinamica(char *File)
Read the config file.
Definition: Forces.cpp:585
int IfRot
Boh.
Definition: Forces.h:829
int NSlab
of slabs
Definition: Forces.h:195
int SetNNano(int Val)
Set NNano.
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 Print()
Print the entries.
int ThermMode
Thermostat mode.
Definition: Forces.h:749
TENS Tens
Structure for the pressure calculation.
Definition: Forces.h:783
PART * Pl
Splines.
Definition: Forces.h:787
int IfNano
Boh.
Definition: Forces.h:761
DomDec * Pc
Pair list.
Definition: Forces.h:785
double NrgPBead
Average energy per particle.
Definition: Forces.h:735
double SLap
Prefactor of the square laplacian.
Definition: Forces.h:149
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
void DefNanoForceParam()
Define the parameters for calculating the force.
double ForThr
Maximum force allowed.
Definition: Forces.h:171
int Bead2Move
Bead to move.
Definition: Forces.h:737
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
int pNPart()
Number of particle.
int NSpline
Total number of points for drawing a spline.
Definition: Forces.h:777
void SetkBen(double Val)
Bending coupling.
Definition: VarData.h:936