Allink  v0.1
ForcesForceField.cpp
1 #include "Forces.h"
2 //-------------------------SPLINES--------------------------------
4 int Forces::ForceFieldLine(){
5  double Sign=1.;
6  SPLINE Sp1;
7  SPLINE Sp2;
8  SPLINE Cub1;
9  SPLINE Cub2;
10  SPLINE Par1;
11  SPLINE Par2;
12  double fSLap=0.;
13  double fLap=0.;
14  //Fm[Part2Move].fExt[2] = -Fm[Part2Move].fHel[2];
15  //Fm[Part2Move].fExt[2] = 0.;
16  for(int p=0;p<pNPart();p++){
17  Fm[p].Ext[2] = -Fm[p].Dir[2]/1.1;
18  if(IfInterp==FIT_SPLINE3){
19  if( p == 0)
20  Sp2 = Mat->Spline3Beg(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
21  else if( p < pNPart() - 2 )
22  Sp2 = Mat->Spline3(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
23  else if( p == pNPart() - 2 )
24  Sp2 = Mat->Spline3End(Pm[p].Pos,Pm[p+1].Pos,0,2);
25  fSLap = 0.;//atan(Sp2.a1);
26  fLap = 2.*Sp2.a2;
27  }
28  else if(IfInterp==FIT_SPLINE4){
29  if( p == 0){
30  Sp2 = Mat->Spline4Beg(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,2);
31  //Pm[p].Typ = 2;
32  }
33  else if( p < pNPart() - 3 )
34  Sp2 = Mat->Spline4(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,2);
35  else if( p == pNPart() - 3 ){
36  Sp2 = Mat->Spline3(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
37  //Sp2 = Mate->Spline4PreEnd(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,1);
38  //Pm[p].Typ = 2;
39  }
40  else if( p == pNPart() - 2 ){
41  Sp2 = Mat->Spline3End(Pm[p].Pos,Pm[p+1].Pos,0,2);
42  //Sp2 = Mate->Spline4End(Pm[p].Pos,Pm[p+1].Pos,0,1);
43  //Pm[p].Typ = 2;
44  }
45  fSLap = (24.*Sp2.a4);
46  fLap = (2.*Sp2.a2);
47  }
48 // else if(IfInterp==FIT_PARAB){
49 // if(p < Gen->NPart - 2){
50 // Par1 = Mate->Parab(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
51 // }
52 // fSLap = 0.;//(2.*Par2.A*Pm[p].Pos[0]+Par2.B);
53 // fLap = 2.*Par1.A/100.;
54 // }
55  else if(IfInterp==FIT_PARAB2){
56  if(p == 0){
57  Sp2 = Mat->Parab(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
58  }
59  if(p < pNPart() - 1){
60  Sp2 = Mat->Parab2(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,0,2);
61  }
62  fSLap = 0.;//(Par2.B);
63  fLap = 2.*Sp2.a2;
64  }
65  else if(IfInterp==FIT_CUBIC){
66  if(p == 0){
67  Cub2 = Mat->Parab(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
68  }
69  if(p < pNPart() - 2){
70  Cub2 = Mat->Cubica(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
71  }
72  if(p == pNPart() - 2){
73  Cub2 = Mat->Parab2(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,0,2);
74  }
75  fSLap = 0.;
76  fLap = 2.*Cub2.a2;
77  }
78  else if(IfInterp==FIT_FORTH){
79  double Ref=0.;
80  if(p < 1){
81  //Sp1 = Mat->Forth(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,Pm[p+4].Pos,0,1);
82  Sp2 = Mat->Forth(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,Pm[p+4].Pos,0,2);
83  Ref = Pm[p].Pos[0] - Pm[p+2].Pos[0];
84  }
85  else if(p < 2){
86  //Sp1 = Mat->Forth(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,1);
87  Sp2 = Mat->Forth(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,2);
88  Ref = Pm[p].Pos[0] - Pm[p+1].Pos[0];
89  }
90  else if(p < pNPart()-2){
91  //Sp1 = Mat->Forth(Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,1);
92  Sp2 = Mat->Forth(Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
93  Ref = 0.;
94  }
95  else if(p == pNPart()-2){
96  //Sp1 = Mat->Forth(Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,0,1);
97  Sp2 = Mat->Forth(Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,0,2);
98  Ref = Pm[p].Pos[0] - Pm[p-1].Pos[0];
99  }
100  else if(p <= pNPart()-1){
101  //Sp1 = Mat->Forth(Pm[p-4].Pos,Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,0,1);
102  Sp2 = Mat->Forth(Pm[p-4].Pos,Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,0,2);
103  Ref = Pm[p].Pos[0] - Pm[p-2].Pos[0];
104  }
105  fSLap = -(24.*Sp2.a4)*QUAD(QUAD(Dx));
106  fLap = (2.*Sp2.a2)*QUAD(Dx);
107  //printf("%d) %lf %lf\n",p,Kf.kLap*fLap,Kf.kSLap*fSLap);
108  }
109  if(fLap > 100.)
110  fLap = 100.;
111  if(fSLap < -100.)
112  fSLap = -100.;
113  if(fSLap > 100.)
114  fSLap = 100.;
115  if(fSLap < -100.)
116  fSLap = -100.;
117  Fm[p].Dir[2] = Kf.SLap*fSLap + Kf.Lap*fLap;
118  //if(Pm[p].Pos[2] <0) Sign = -1.;
119  //else Sign = 1;
120  //Fm[p].fHel[2] *= Sign;
121  //if(p == Part2Move)
122  //if(p > Part2Move -4 && p < Part2Move +4)
123  {
124  //printf("%d in (%lf,%lf) ha %lf|%lf Sp %lf %lf\n",p,Pm[p].Pos[2],Pm[p].Vel[2],Fm[p].fHel[2],Fm[p].fExt[2],fSLap,fLap);
125  }
126  }
127  return 0;
128 }
130 int Forces::ForceFieldLeaves(){
131  double Sign=1.;
132  SPLINE Sp1;
133  SPLINE Sp2;
134  PARABOLA Par1;
135  PARABOLA Par2;
136  double fSLap=0.;
137  double fLap=0.;
138  for(int c=0;c<pNChain();c++){
139  for(int p=c*NEdge;p<NEdge*(c+1);p++){
140  Fm[p].Ext[2] = -(Fm[p].Dir[2])/1.1;
141  // Fm[p].El[2] = - Kf.El[2]*( c*Kf.Elong[2] + .5 - Pm[p].Pos[2]);
142  if( p == c*NEdge){
143  Sp2 = Mat->Spline4Beg(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,2);
144  //Pm[p].Typ = 2;
145  }
146  else if( p < NEdge*(c+1) - 3 )
147  Sp2 = Mat->Forth(Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
148  //Sp2 = Mat->Spline4(Pm[p].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
149  // else Pm[p].Typ = 2;
150  else if( p == NEdge*(c+1) - 3 ){
151  Sp2 = Mat->Spline3(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
152  //Pm[p].Typ = 2;
153  }
154  else if( p == pNPart() - 2 ){
155  Sp2 = Mat->Spline3End(Pm[p].Pos,Pm[p+1].Pos,0,2);
156  //Pm[p].Typ = 2;
157  }
158  fSLap = -(24.*Sp2.a4)*QUAD(QUAD(Dx));
159  fLap = (2.*Sp2.a2)*QUAD(Dx);
160  if(fLap > 100.)fLap = 100.;
161  if(fLap < -100.)fLap = -100.;
162  if(fSLap > 100.)fSLap = 100.;
163  if(fSLap < -100.)fSLap = -100.;
164  Fm[p].Dir[2] = +Kf.SLap*fSLap + Kf.Lap*fLap;
165  continue;
166  {
167  double Rad2=0.;
168  for(int d=0;d<3;d++){
169  Rad2 += QUAD((Pm[p].Pos[d] - Nano->Pos[d]));
170  }
171  if( Rad2 > QUAD((3.*Nano->Rad))) continue;
172  double CosTh = (Nano->Pos[0] - Pm[p].Pos[0])/Nano->Rad;
173  double SenTh = sqrt(1 - QUAD(CosTh));
174  double Sign = 1.;
175  double Rad = sqrt(Rad2);
176  double Strength = 12.*pow( (Nano->Rad/Rad) , 13) + 6.*pow( (Nano->Rad/Rad) , 7);
177  if(p <= NEdge)
178  Sign *= -1.;
179  Fm[p].Dir[0] = Kf.LJ * Strength * CosTh;
180  Fm[p].Dir[2] = Kf.LJ * Strength * SenTh * Sign;
181  }
182  }
183  }
184 }
185 //----------------------HARMONIC-POTENTIAL-------------------
187 int Forces::ForceFieldBulk(){
188  for(int p=0;p<pNPart();p++){
189  for(int d=0;d<3;d++){
190  Fm[p].Ext[d] = -(Fm[p].Dir[d])/2.;
191  }
192  for(int l=0;l<Ln[p].NLink;l++){
193  int link = Ln[p].Link[l];
194  for(int d=0;d<3;d++){
195  double Delta = ASS((Pm[link].Pos[d] - Pm[p].Pos[d]));
196  if(Delta < 0.001) continue;
197  // printf("%d-%d (%d) %lf ",p,link,d,Delta);
198  double Sign = 1.;
199  if(Pm[p].Pos[d] < Pm[link].Pos[d] )
200  Sign = -1.;
201  Fm[p].Dir[d] += - Sign*Kf.El[d]*(Delta - Kf.Elong[d])/1.1;
202  }
203  // printf("\n");
204  }
205  //if(p > Part2Move -4 && p < Part2Move +4)
206  //printf("%d) in (%lf,%lf) ha %lf|%lf|%lf Sp %lf %lf\n",p,Pm[p].Pos[2],Pm[p].Vel[2],Fm[p].fHel[2],Fm[p].fExt[2],Fm[p].fEl[2]);
207  }
208 }
210 void Forces::ForceFieldRod(){
211  double DistRelBA[4];
212  double DistRelBC[4];
213  double Pre[12];
214  for(int b=0;b<pNBlock();b++){
215  for(int c=0;c<pNChain(b);c++){
216  for(int p=c*pNPCh(b);p<(c+1)*pNPCh(b)-1;p++){
217  TwoPartDist(p+1,p,DistRelBA);
218  double ForceSp = pkSpr()*(1. - pSprRest()/DistRelBA[3]);
219  for(int d=0;d<3;d++){
220  Fm[p].Dir[d] += ForceSp*DistRelBA[d];
221  Fm[p+1].Dir[d] -= ForceSp*DistRelBA[d];
222  }
223  if(p < (c+1)*pNPCh(b)-2){
224  TwoPartDist(p+2,p+1,DistRelBC);
225  double CosAngle = 0.;
226  for(int d=0;d<3;d++){
227  DistRelBA[d] /= DistRelBA[3];
228  DistRelBC[d] /= DistRelBC[3];
229  CosAngle += DistRelBA[d]*DistRelBC[d];
230  }
231  double PreFactBA = pkBen()/DistRelBA[3];
232  double PreFactBC = pkBen()/DistRelBC[3];
233  for(int d=0;d<3;d++){
234  Fm[p+0].Dir[d] += PreFactBA*(DistRelBC[d]-DistRelBA[d]*CosAngle);
235  Fm[p+1].Dir[d] -= PreFactBA*(DistRelBC[d]-DistRelBA[d]*CosAngle);
236  Fm[p+1].Dir[d] += PreFactBC*(DistRelBA[d]-DistRelBC[d]*CosAngle);
237  Fm[p+2].Dir[d] -= PreFactBC*(DistRelBA[d]-DistRelBC[d]*CosAngle);
238  }
239  }
240  }
241  }
242  }
243 }
244 /* obsolete? */
245 int Forces::SumSomeForces(int HowMany){
246  double Sign=1.;
247  SPLINE Par1;
248  SPLINE Par2;
249  double Slope1;
250  double Slope2;
251  if(Bead2Move - HowMany < 1 || Bead2Move + HowMany > pNPart() - 2)
252  return 1;
253  //Fm[Part2Move].fExt[2] = -Fm[Part2Move].fHel[2];
254  //Fm[Part2Move].fExt[2] = 0.;
255  int pp = Bead2Move;
256  for(int h = 0;h<HowMany;h++){
257  int link1 = Ln[pp].Link[0];
258  int link2 = Ln[pp].Link[1];
259 // printf("(%lf %lf %lf), (%lf %lf %lf), (%lf %lf %lf)\n",
260 // Pm[link1].Pos[0],Pm[link1].Pos[1],Pm[link1].Pos[2],
261 // Pm[pp].Pos[0],Pm[pp].Pos[1],Pm[pp].Pos[2],
262 // Pm[link2].Pos[0],Pm[link2].Pos[1],Pm[link2].Pos[2]);
263 // Par1 = Mat->Parab2(Pm[link1].Pos,Pm[pp].Pos,Pm[link2].Pos,0,2);
264  Par2 = Mat->Parab2(Pm[link1].Pos,Pm[pp].Pos,Pm[link2].Pos,1,2);
265  Fm[pp].Dir[2] = - Kf.Lap*QUAD(2.*Par1.a2);
266  Fm[pp].Dir[2] += - Kf.Lap*QUAD(2.*Par2.a2);
267  //printf("%d-%d|%d) %lf %lf %lf %lf ",pp,link1,link2,Slope1,Par1.A,Slope2,Par2.A);
268  link1 = Ln[pp].Link[2];
269  link2 = Ln[pp].Link[3];
270  Par1 = Mat->Parab2(Pm[link1].Pos,Pm[pp].Pos,Pm[link2].Pos,0,2);
271  //Par2 = Mat->Parab2(Pm[link1].Pos,Pm[pp].Pos,Pm[link2].Pos,1,2);
272  Fm[pp].Dir[2] += - Kf.Lap*QUAD(2.*Par1.a2);
273  Fm[pp].Dir[2] += - Kf.Lap*QUAD(2.*Par2.a2);
274  //printf("%d-%d|%d) %lf %lf %lf %lf\n",pp,link1,link2,Slope1,Par1.A,Slope2,Par2.A);
275  if(Pm[pp].Pos[2] <0) Sign = -1.;
276  else Sign = 1;
277  Fm[pp].Dir[2] *= Sign;
278  //printf("%d) Pos %lf Force %lf %lf\n",pp,Pm[pp].Pos[2],Fm[pp].fHel[2],Fm[pp].fExt[2]);
279  }
280  return 0;
281 }
282 //-----------------------------RIGID-------------------------
283 void Forces::ForceFieldRigid(){
284  double Pot[2];
285  for(int n=0;n<pNNano();n++){
286  for(int d=0;d<3;d++){
287  Nano[n].Force[d] = 0.;
288  Nano[n].AMom[d] = 0.;
289  }
290  }
291  for(int n=0;n<pNNano();n++){
292  for(int nn=n+1;nn<pNNano();nn++){
293  double dr[3] = {0.,0.,0.};
294  double Dist = RigidDistanceAxis(n,nn,dr);
295  //if( Dist > CutOff ) continue;
296  double Cons = RigidLJ(n,Dist,Pot,0);
297  for(int d=0;d<3;d++){
298  Nano[n].Force[d] -= Cons*dr[d];
299  Nano[nn].Force[d] += Cons*dr[d];
300  Nano[n].AMom[d] -= Nano[n].AMomTemp[d]*Cons;
301  Nano[nn].AMom[d] += Nano[n].AMomTemp[d]*Cons;
302  //printf("%d %lf %lf\n",d,Cons*dr[d],Nano[nn].AMom[d]);
303  }
304  }
305  for(int d=0;d<3;d++){
306  double Ran = Nano[n].Zeta * (2.*Mat->Casuale() - 1.);
307  double Dis = - Nano[n].Gamma*Nano[n].Vel[d];
308  //Nano->Force[d] += Ran + Dis;
309  }
310  //printf("%d %lf %lf %lf\n",n,Nano[n].AMom[0],Nano[n].AMom[1],Nano[n].AMom[2]);
311  }
312 }
313 double Forces::RigidDistanceAxis(int n,int nn,double *dr){
314  Vettore Pos1(Nano[n].Pos,3);
315  Vettore Pos2(Nano[nn].Pos,3);
316  Vettore PosRel(Pos2[0]-Pos1[0],Pos2[1]-Pos1[1],Pos2[2]-Pos1[2]);
317  Vettore Axis(Nano[n].Axis,3);
318  Vettore Perp(3);
319  Vettore AMomTemp(3);
320  double Dist = Perp.PerpTo(&PosRel,&Axis);
321  Perp.Export(dr);
322  AMomTemp.VetV(&Perp,&PosRel);
323  AMomTemp.Export(Nano[n].AMomTemp);
324  //AMomTemp.Export(Nano[nn].AMomTemp);
325  double HeiOnAxis = Perp.ProjOnAxis(&PosRel,&Axis);
326  if( fabs(HeiOnAxis) > Nano[n].Height*.5){
327  double Sign = HeiOnAxis > 0 ? 1. : -1.;
328  Dist = 0.;
329  for(int d=0;d<3;d++){
330  dr[d] = Pos2[d] - Pos1[d];
331  Dist = SQR(dr[d]);
332  }
333  Dist = sqrt(Dist);
334  }
335  return Dist;
336 }
337 void Forces::NanoInteraction(){
338  double Pot[2];
339  for(int n=0;n<pNNano();n++){
340  for(int d=0;d<3;d++){
341  Nano[n].Force[d] = 0.;
342  Nano[n].AMom[d] = 0.;
343  }
344  }
345  for(int n=0;n<pNNano();n++){
346  Point2Shape(Nano[n].Shape);
347  for(int p=0;p<pNPart();p++){
348  double dr[3] = {0.,0.,0.};
349  double Dist2 = NanoDist2(Pm[p].Pos,n);
350  if( Dist2 > Kf.CutOff2 ) continue;
351  double Cons = RigidLJ(n,Dist2,Pot,0);
352  for(int d=0;d<3;d++){
353  Nano[n].Force[d] -= Cons*dr[d];
354  Fm[p].Dir[d] += Cons*dr[d];
355  Nano[n].AMom[d] -= Nano[n].AMomTemp[d]*Cons;
356  }
357  }
358  for(int d=0;d<3;d++){
359  double Ran = Nano[n].Zeta * (2.*Mat->Casuale() - 1.);
360  double Dis = - Nano[n].Gamma*Nano[n].Vel[d];
361  //Nano->Force[d] += Ran + Dis;
362  }
363  //printf("%d %lf %lf %lf\n",n,Nano[n].AMom[0],Nano[n].AMom[1],Nano[n].AMom[2]);
364  }
365 }
366 double Forces::NanoNrg(int p){
367  double Nrg = NanoNrg(pPos(p),pType(p));
368  return Nrg;
369 }
370 double Forces::NanoNrg(double *Pos,int t){
371  double Pot[2];
372  double Cons = 0.;
373  double Nrg = 0.;
374  for(int n=0;n<pNNano();n++){
375  Point2Shape(Nano[n].Shape);
376  double Dist2 = NanoDist2(Pos,n);
377  if( Dist2 > SQR(Nano[n].CutOff) ) continue;
378  double Dist = sqrt(Dist2);
379  double Sign = -1.;
380  if(t == 1) Sign = 1.;
381  Cons += RigidHamaker(n,Dist,Pot,Sign);
382  //Cons += RigidLJ(n,Dist2,Pot,0);
383  Nrg += Pot[0];
384  }
385  Nrg *= Nano->Coating;
386  return Nrg;
387 }
388 void Forces::DefForceParam(){
389  double Pot[2];
390  Kf.ForThr = 50.;
391  Kf.DistThr = 0.;
392  Potential(Kf.CutOff2,0,0,Pot);
393  Kf.BaseLine = -Pot[0];
394  double CutOff = sqrt(Kf.CutOff2);
395  for(int i=10000;i>0;i--){
396  double x = CutOff/10000.*i;
397  double For = Potential(SQR(x),0,0,Pot)/x;
398  if(For >= Kf.ForThr){
399  Kf.DistThr = x;
400  Kf.PotThr = Pot[0];
401  break;
402  }
403  }
404 }
405 void Forces::DefNanoForceParam(){
406  double Pot[2];
407  for(int n=0;n<pNNano();n++){
408  Nano[n].ForThr = 50.;
409  Point2Shape(Nano[n].Shape);
410  Nano[n].CutOff = Nano[n].Rad*3.;
411  Nano[n].BaseLine = 0.;
412  double For = RigidLJ(n,Nano[n].CutOff,Pot,0);
413  Nano[n].BaseLine = -Pot[0];
414  Nano[n].DistThr = 0.;
415  for(int i=10000;i>0;i--){
416  double x = (Nano[n].CutOff)/10000.*i;
417  For = RigidHamaker(n,x,Pot,0);
418  //For = RigidLJ(n,x,Pot,0);
419  if(For >= Nano[n].ForThr){
420  Nano[n].DistThr = x;
421  Nano[n].PotThr = Pot[0];
422  break;
423  }
424  }
425  //printf("DefForce %lf %lf %lf\n",Nano[n].DistThr,Nano[n].PotThr,Nano[n].BaseLine);
426  }
427 }
428 double Forces::RigidDistanceRad(int n,int nn,double *dr){
429  Vettore Pos1(Nano[n].Pos,3);
430  Vettore Pos2(Nano[nn].Pos,3);
431  double Dist = 0.;
432  for(int d=0;d<3;d++){
433  dr[d] = Pos2[d] - Pos1[d];
434  Dist += SQR(dr[d]);
435  }
436  return sqrt(Dist);
437 }
438 double Forces::RigidLJ(int n,double Dist,double *Pot,double Sign){
439  double Cons = 0.;
440  double Strength = pow(Nano[n].Coating,9.)*Nano[n].Hamaker;
441  double Hamaker = pow(Nano[n].Coating,3.)*Nano[n].Hamaker;
442  if( Dist < Nano[n].DistThr){
443  Pot[0] = Nano[n].PotThr - Nano[n].ForThr*(Dist - Nano[n].DistThr);
444  Pot[1] = Nano[n].ForThr;
445  Cons = Nano[n].ForThr;
446  return Cons;
447  }
448  double idr = 1./(Dist-Nano[n].Rad+Nano[n].Coating);
449  double idr2 = idr*idr;
450  double idr4 = idr2*idr2;
451  double idr6 = idr2*idr4;
452  Pot[0] = idr*idr2*(Strength*idr6 +Sign*Hamaker)+Nano[n].BaseLine;
453  Cons = idr4 *(Strength*9.*idr6+Sign*3.*Hamaker);
454  Pot[1] = idr4 *(Strength*9.*idr6+Sign*3.*Hamaker);
455  return Cons;
456 }
458 double Forces::RigidHamaker(int n,double Dr,double *Pot,double Sign){
459  double SurfTens = 3.*Nano[n].Hamaker;
460  double Cons = 0.;
461  double Sigma = 1.0;//0.04;
462  double Slope = 1.00203;
463  double Intercept = 0.31739;
464  double Rnp = Nano[n].Rad/Slope-Intercept;
465  if( Dr < Nano[n].DistThr){
466  Pot[0] = Nano[n].PotThr - Nano[n].ForThr*(Dr - Nano[n].DistThr);
467  Pot[1] = Nano[n].PotThr;
468  return Nano[n].ForThr;
469  }
470  double DrInv = 1./Dr;
471  double DrP1 = 1./(Dr+Rnp);
472  double DrP3 = CUBE(DrP1);
473  double DrP6 = SQR(DrP3);
474  double DrP9 = DrP6*DrP3;
475  double DrM1 = 1./(-Dr+Rnp);
476  double DrM3 = CUBE(DrM1);
477  double DrM6 = SQR(DrM3);
478  double DrM9 = DrM6*DrM3;
479  double PreRep = Sigma*1./360.*DrInv*SurfTens;
480  double PreAttr = 1./12.*DrInv*SurfTens;
481  double Rep = (9.0*Rnp+Dr)*DrP9 - (9.0*Rnp-Dr)*DrM9;
482  double Attr= (3.0*Rnp+Dr)*DrP3 - (3.0*Rnp-Dr)*DrM3;
483  double RepP = -DrP9+9.*(9.*Rnp+Dr)*DrP9*DrP1
484  -(9.*Rnp+Dr)*DrP9*DrInv;
485  double RepM = -DrM9+9.*(9.*Rnp-Dr)*DrM9*DrM1
486  -(9.*Rnp-Dr)*DrM9*DrInv;
487  double AttrP = -DrP3+3.*(3.*Rnp+Dr)*DrP3*DrP1
488  -(3.*Rnp+Dr)*DrP3*DrInv;
489  double AttrM = -DrM3+3.*(3.*Rnp-Dr)*DrM3*DrM1
490  -(3.*Rnp-Dr)*DrM3*DrInv;
491  double RepChem = 9.*DrP9 - (9.*Rnp+Dr)*9.*DrP9*DrP1
492  - 9.*DrM9 + 9.*(9.*Rnp-Dr)*DrM9*DrM1;
493  double AttrChem = 3.*DrP3 - (3.*Rnp+Dr)*3.*DrP3*DrP1
494  - 3.*DrM9 + 3.*(3.*Rnp-Dr)*DrM3*DrM1;
495  Pot[0] = PreRep*Rep + Sign*PreAttr*Attr + Nano[n].BaseLine;
496  Pot[1] = PreRep*RepChem + Sign*PreAttr*AttrChem;
497  Cons = PreRep*(RepP+RepM) + Sign*PreAttr*(AttrP+AttrM);
498  return Cons;
499 }
500 double Forces::RigidCoulomb(int n,double Dist,double *Pot,double Sign){
501  if(Dist < 2.*Nano[n].Rad)
502  return 300.;
503  Pot[0] = Sign*Nano[n].Hamaker/(Dist);
504  return -Sign*Nano[n].Hamaker/SQR(Dist);
505 }
506 //----------------------------DEF-FORCES-POTENTIAL-------------
507 void Forces::PointShape(int iShape){
508  Point2Shape(iShape);
509  if(VAR_IF_TYPE(iShape,SHAPE_SPH))
510  CalcPot = &Forces::LJ39;
511  if(VAR_IF_TYPE(iShape,SHAPE_CYL))
512  CalcPot = &Forces::LJ39;
513  if(VAR_IF_TYPE(iShape,SHAPE_PILL))
514  CalcPot = &Forces::LJ39;
515  if(VAR_IF_TYPE(iShape,SHAPE_TILT))
516  CalcPot = &Forces::LJ39;
517  if(VAR_IF_TYPE(iShape,SHAPE_WALL))
518  CalcPot = &Forces::LJ39;
519  if(VAR_IF_TYPE(iShape,SHAPE_PORE))
520  CalcPot = &Forces::LJ39;
521  if(VAR_IF_TYPE(iShape,SHAPE_EXT))
522  CalcPot = &Forces::LJ39;
523  if(VAR_IF_TYPE(iShape,SHAPE_JANUS))
524  CalcPot = &Forces::LJ39;
525  if(VAR_IF_TYPE(iShape,SHAPE_STALK))
526  CalcPot = &Forces::LJ39;
527  if(VAR_IF_TYPE(iShape,SHAPE_TIP))
528  CalcPot = &Forces::LJ39;
529  if(VAR_IF_TYPE(iShape,SHAPE_TORUS))
530  CalcPot = &Forces::LJ39;
531  if(VAR_IF_TYPE(iShape,SHAPE_HARM))
532  CalcPot = &Forces::LJ39;
533  if(VAR_IF_TYPE(iShape,SHAPE_UMBR))
534  CalcPot = &Forces::LJ39;
535  DefNanoForceParam();
536 }
537 double Forces::Harmonic(double Dist2,int t1,int t2,double *Pot){
538  double Dist = sqrt(Dist2);
539  double a = Kf.LJ;
540  double b = -2.*a*Kf.LJMin;
541  double c = -a*SQR(Kf.LJMin) - b*Kf.LJMin;
542  Pot[2] = a*Dist2 + b*Dist + c;
543  double For = - 2.*a*Dist2 - b*Dist;
544  if(t1 != t2){
545  Pot[2] *= -1.;
546  return -For;
547  }
548  //Pot[2] = 0.;
549  return For;
550 }
551 double Forces::LJ39(double Dist2,int t1,int t2,double *Pot){
552  double idr6 = CUBE(SQR(Kf.LJMin)/Dist2);
553  double idr3 = CUBE(Kf.LJMin)/(Dist2*sqrt(Dist2));
554  // if( Dist2 < Kf.DistThr){
555  // Pot[0] = Kf.PotThr - Kf.ForThr*(Dist2 - Kf.DistThr);
556  // return Cons = Kf.ForThr;
557  // }
558  if(t1 == t2){
559  Pot[2] = Kf.LJ*idr3*(idr6-1.)+Kf.BaseLine;
560  return Kf.LJ*idr3*(9.*idr6-3.);
561  }
562  else{
563  Pot[2] = Kf.LJ*idr3*(idr6+1.)+Kf.BaseLine;
564  return Kf.LJ*idr3*(9.*idr6+3.);
565  }
566 }
567 double Forces::LJPot(double Dist2,int t1,int t2,double *Pot){
568  double idr6 = CUBE(SQR(Kf.LJMin)/Dist2);
569  // if( Dist2 < Kf.DistThr){
570  // Pot[0] = Kf.PotThr - Kf.ForThr*(Dist2 - Kf.DistThr);
571  // return Kf.ForThr;
572  // }
573  if(t1 == t2){
574  Pot[2] = Kf.LJ*4.*(idr6*idr6-idr6) + Kf.BaseLine;
575  return Kf.LJ*48.*(idr6*idr6-.5*idr6);
576  }
577  else{
578  Pot[2] = Kf.LJ*4.*(idr6*idr6+idr6) + Kf.BaseLine;
579  return Kf.LJ*48.*(idr6*idr6+.5*idr6);
580  }
581 }
582 double Forces::StepPot(double Dist2,int t1,int t2,double *Pot){
583  if(t1 != t2){
584  Pot[2] = -Kf.LJ;
585  return -Kf.LJ*sqrt(Dist2/Kf.CutOff2);
586  }
587  Pot[2] = Kf.LJ;
588  return Kf.LJ*sqrt(Dist2/Kf.CutOff2);
589 }
590 double Forces::ElectroPot(double Dist2,int t1,int t2,double *Pot){
591  if(t1 == t2)
592  Pot[2] = Kf.Ext/Dist2;
593  else
594  Pot[2] = -Kf.LJ/Dist2;
595  return Pot[2];
596 }
597 void Forces::CalcForcesDensFunc(){
598  if(!VAR_IF_TYPE(SysAlloc,ALL_FORCES)){
599  printf("Forces not allocated\n");
600  return;
601  }
602  int NTabF = 300;
603  TabForceAlloc(NTabF);
604  double *Count = (double *)calloc(NTab*pNType()*pNType(),sizeof(double));
605  ClearDens();
606  AddDens(0,pNPart());
607  SumDens(0,pNPart());
608  double Dist = 0.;
609  double DistRelBA[4];
610  double InvCutOff = 1./sqrt(Kf.CutOff2);
611  for(int p=0;p<pNPart();p++){
612  for(int d=0;d<3;d++){
613  Fm[p].Dir[d] = 0.;
614  }
615  }
616  // non bonded
617  for(int p1=0;p1<pNPart();p1++){
618  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
619  int p2 = Pc->p2Curr;
620  if(p2 <= p1) continue;
621  Pc->Dist2Curr(DistRelBA);
622  if(DistRelBA[3] > Kf.CutOff2) continue;
623  double Force = 0.;
624  double Dist = sqrt(DistRelBA[3]);
625  for(int t=0;t<pNType();t++){
626  Force += MInt->Coeff(Pm[p1].Typ,Pm[p2].Typ,t)*(Dens3[p1*pNType()+t]+Dens3[p2*pNType()+t]);
627  }
628  Force *= DerWei3(Dist,pWei3Par())*2./3.;
629  Force += DerWei2(Dist,pWei2Par())*MInt->Coeff(Pm[p1].Typ,Pm[p2].Typ);
630  Force /= -Dist;
631  int x = (int)(Dist*InvCutOff*NTab);
632  if( x < 0 || x >= NTab) continue;
633  int t = MIN(Pm[p1].Typ,Pm[p2].Typ)*pNType()+MAX(Pm[p1].Typ,Pm[p2].Typ);
634  FTab[x*pNType()*pNType()+t] += Force;
635  Count[x*pNType()*pNType()+t] += 1.;
636  for(int d=0;d<3;d++){
637  Fm[p1].Dir[d] += Force*DistRelBA[d];
638  Fm[p2].Dir[d] -= Force*DistRelBA[d];
639  }
640  }
641  }
642  for(int t=0;t<NTab;t++){
643  for(int t1=0;t1<pNType();t1++){
644  for(int t2=0;t2<pNType();t2++){
645  int tType = MIN(t1,t2)*pNType()+MAX(t1,t2);
646  double Norm = Count[t*pNType()*pNType()+tType] > 0. ? 1./Count[t*pNType()*pNType()+tType] : 1.;
647  FTab[t*pNType()*pNType()+tType] *= Norm;
648  }
649  }
650  }
651  // bonded
652  double DistRelBC[4];
653  double Pre[12];
654  for(int b=0;b<pNBlock();b++){
655  for(int c=0;c<pNChain(b);c++){
656  for(int p=c*pNPCh(b);p<(c+1)*pNPCh(b)-1;p++){
657  TwoPartDist(p+1,p,DistRelBA);
658  double ForceSp = pkSpr()*(1. - pSprRest()/DistRelBA[3]);
659  for(int d=0;d<3;d++){
660  Fm[p].Dir[d] += ForceSp*DistRelBA[d];
661  Fm[p+1].Dir[d] -= ForceSp*DistRelBA[d];
662  }
663  if(p < (c+1)*pNPCh(b)-2){
664  TwoPartDist(p+2,p+1,DistRelBC);
665  double CosAngle = 0.;
666  for(int d=0;d<3;d++){
667  DistRelBA[d] /= DistRelBA[3];
668  DistRelBC[d] /= DistRelBC[3];
669  CosAngle += DistRelBA[d]*DistRelBC[d];
670  }
671  double PreFactBA = pkBen()/DistRelBA[3];
672  double PreFactBC = pkBen()/DistRelBC[3];
673  for(int d=0;d<3;d++){
674  Fm[p+0].Dir[d] -= PreFactBA*(DistRelBC[d]-DistRelBA[d]*CosAngle);
675  Fm[p+1].Dir[d] += PreFactBA*(DistRelBC[d]-DistRelBA[d]*CosAngle);
676  Fm[p+1].Dir[d] -= PreFactBC*(DistRelBA[d]-DistRelBC[d]*CosAngle);
677  Fm[p+2].Dir[d] += PreFactBC*(DistRelBA[d]-DistRelBC[d]*CosAngle);
678  }
679  }
680  }
681  }
682  }
683 }
684 void Forces::PrintForce(){
685  double NPoint = 1000.;
686  FILE *FWrite = fopen("Force.dat","w");
687  double Pot[2];
688  double Delta = sqrt(Kf.CutOff2)/NPoint;
689  double CutOff = sqrt(Kf.CutOff2);
690  for(double x=Delta;x<CutOff;x+=Delta){
691  double Cons = LJPot(x,0,0,Pot);
692  fprintf(FWrite,"%lf %lf %lf\n",x,Cons,*Pot);
693  }
694  fclose(FWrite);
695 }
696 void Forces::TabForceAlloc(int NTabExt){
697  NTab = NTabExt;
698  double Pot[2];
699  double Delta = sqrt(Kf.CutOff2)/(double)NTab;
700  double NrgLim = 20.;
701  FTab = (double *)calloc(NTab*pNType()*pNType(),sizeof(double));
702  MTab = (double *)calloc(NTab,sizeof(double));
703  PTab = (double *)calloc(NTab*pNType()*pNType(),sizeof(double));
704  VAR_ADD_TYPE(DynFlag,ALL_FORCES);
705 }
706 void Forces::TabPot(){
707  NTab = 10000;
708  double NrgLim = 10.;
709  double Delta = 2.*NrgLim/(double)NTab;
710  int i = 0;
711  for(double Nrg=-NrgLim;Nrg<NrgLim;Nrg+=Delta,i++){
712  PTab[i] = Nrg;
713  MTab[i] = exp(-pTemp()*Nrg);
714  }
715  VAR_ADD_TYPE(DynFlag,ALL_METR);
716  VAR_ADD_TYPE(DynFlag,ALL_POT);
717 }
718 double Forces::Wei3(const double r, const double b){
719  if(r < 0.0001) return 0.;
720  return (r < b) ? 15. / (2. * M_PI * CUBE(b)) * SQR(1. - r / b) : 0.;
721 }
722 double Forces::DerWei3(const double r, const double b){
723  if(r < 0.0001) return 0.;
724  return (r < b) ? 15. / (M_PI * SQR(SQR(b))) * (r / b - 1.) : 0.;
725 }
726 double Forces::Wei2(const double r, const double a){
727  if(r < 0.0001) return 0.;
728  const double b = 1.;
729  double Ai = .5 * 15. / M_PI / (CUBE(a) + CUBE(a + b) + CUBE(b));
730  return (r < a) ? Ai : (Ai / CUBE(a - b)) *
731  ((((-2. * r) + 3. * (a + b)) * r - 6. * a * b) * r +
732  3. * a * b * b - b * b * b);
733 }
734 double Forces::DerWei2(const double r, const double a){
735  if(r < 0.0001) return 0.;
736  const double b = 1.;
737  double Ai = 15. / (M_PI * (((((4. * a) + 6.) * a) + 6.) * a + 4.));
738  return (r > a && r < b) ? (6. * Ai / CUBE(a - b) * (-r * r + (a + b) *
739  r - a * b)) : 0.;
740 }
Vettore::Export
void Export(double *x)
Export.
Definition: MatematicaVect.cpp:376
Matematica::Spline4Beg
SPLINE Spline4Beg(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline first boundary.
Definition: MatematicaInterp.cpp:68
SPLINE::a4
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Definition: MatematicaStruct.h:114
KFORCES::LJMin
double LJMin
Lennard-Jones.
Definition: Forces.h:161
VarData::pkBen
double pkBen()
Bending coupling.
Definition: VarData.h:932
KFORCES::Ext
double Ext
External force.
Definition: Forces.h:153
Forces::LJPot
double LJPot(double Dist2, int t1, int t2, double *Pot)
Classical Lennard Jones potential.
Definition: ForcesForceField.cpp:567
Forces::ForceFieldLine
int ForceFieldLine()
Helfrich Hamiltonian for a line.
Definition: ForcesForceField.cpp:4
VarData::NanoDist2
double NanoDist2(double *Pos, int n)
Pointer to a generic function.
Definition: VarData.h:764
VarData::Ln
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
Forces::RigidCoulomb
double RigidCoulomb(int nNano, double Dist, double *Pot, double Sign)
Force between to rigid bodies.
Definition: ForcesForceField.cpp:500
VarData::Nano
NANO * Nano
Extra particle.
Definition: VarData.h:1044
Forces::TabForceAlloc
void TabForceAlloc(int NTabExt)
Tabulate the values of the force.
Definition: ForcesForceField.cpp:696
NANO::Gamma
double Gamma
Friction term.
Definition: VarData.h:451
NANO::Hamaker
double Hamaker
Strength of the interaction.
Definition: VarData.h:447
NANO::DistThr
double DistThr
Minimum distance threshold.
Definition: VarData.h:465
Vettore
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
KFORCES::Lap
double Lap
Prefactor of the Laplacian.
Definition: Forces.h:147
Forces::RigidLJ
double RigidLJ(int nNano, double Dist, double *Pot, double Sign)
Force between to rigid bodies.
Definition: ForcesForceField.cpp:438
Matematica::Spline3Beg
SPLINE Spline3Beg(double *P1, double *P2, double *P3, int x, int y)
Three order spline first boundary.
Definition: MatematicaInterp.cpp:15
PARABOLA
Parabolas coefficients.
Definition: MatematicaStruct.h:7
NANO::Rad
double Rad
Size.
Definition: VarData.h:445
Forces::Harmonic
double Harmonic(double Dist2, int t1, int t2, double *Pot)
Harmonic potential.
Definition: ForcesForceField.cpp:537
NANO::Vel
double Vel[3]
Velocity.
Definition: VarData.h:431
VarData::pType
int pType(int p)
Return the type.
Definition: VarDataComm.cpp:219
Forces::DerWei2
double DerWei2(const double r, const double b)
Derivative of the quadratic weighting function.
Definition: ForcesForceField.cpp:734
KFORCES::El
double El[3]
Elastic force.
Definition: Forces.h:151
KFORCES::PotThr
double PotThr
Maximum potential allowed.
Definition: Forces.h:169
PART::Pos
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
Forces::IfInterp
int IfInterp
If interpolates with the splines.
Definition: Forces.h:755
FORCES::Dir
double Dir[3]
Direction.
Definition: Forces.h:178
VarData::pkSpr
double pkSpr()
Spring coupling.
Definition: VarData.h:934
VarData::pSprRest
double pSprRest()
Rest distance of the harmonic potential.
Definition: VarData.h:942
Forces::ForceFieldRod
void ForceFieldRod()
Bending potential on a rod.
Definition: ForcesForceField.cpp:210
DdLinkedList::Dist2Curr
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
Forces::SysAlloc
int SysAlloc
Which arrays are allocated.
Definition: Forces.h:753
Forces::CalcForcesDensFunc
void CalcForcesDensFunc()
Calculate the forces for the density functional.
Definition: ForcesForceField.cpp:597
NANO::Pos
double Pos[3]
Position.
Definition: VarData.h:427
Forces::ElectroPot
double ElectroPot(double Dist2, int t1, int t2, double *Pot)
Potential for the electrical lines.
Definition: ForcesForceField.cpp:590
Forces::NanoInteraction
void NanoInteraction()
Calculate all the interaction with the nano.
Definition: ForcesForceField.cpp:337
Matematica::Parab2
SPLINE Parab2(double *PA, double *PB, double *PC, int x, int y)
Three points parabolic interpolation.
Definition: MatematicaInterp.cpp:199
VarData::pNPCh
int pNPCh()
Number of particle per chain.
Definition: VarDataComm.cpp:249
Forces::ForceFieldRigid
void ForceFieldRigid()
Interaction between rigid bodies.
Definition: ForcesForceField.cpp:283
KFORCES::DistThr
double DistThr
Correspondent radial distance for this force.
Definition: Forces.h:173
Forces::Dx
double Dx
Spatial separation between particles.
Definition: Forces.h:729
VarData::pWei3Par
double pWei3Par()
Parameter of the third order weighting function.
Definition: VarData.h:950
Matematica::Spline4
SPLINE Spline4(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline.
Definition: MatematicaInterp.cpp:87
VarData::TwoPartDist
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
VarData::pNBlock
int pNBlock()
Number of blocks.
Definition: VarDataComm.cpp:254
Forces::NanoNrg
double NanoNrg(int p)
Exchange energy with the nano.
Definition: ForcesForceField.cpp:366
NANO::BaseLine
double BaseLine
Baseline of the potential.
Definition: VarData.h:463
VarData::MInt
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
VarData::NEdge
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
Forces::TabPot
void TabPot()
Tabulate the values of the potential.
Definition: ForcesForceField.cpp:706
Matematica::Casuale
double Casuale()
Random uniform number.
Definition: MatematicaSign.cpp:31
Vettore::ProjOnAxis
double ProjOnAxis(Vettore *a)
Projects along the axis.
Definition: MatematicaVect.cpp:351
Forces::Kf
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
DdLinkedList::NextCurr
void NextCurr()
Increase the iterator to the next couple.
Definition: Cubo.cpp:479
NANO::Zeta
double Zeta
Stochastic term.
Definition: VarData.h:453
NANO::ForThr
double ForThr
Maximum value of the force.
Definition: VarData.h:467
Forces::ForceFieldLeaves
int ForceFieldLeaves()
Helfrich Hamiltonian with an elastic coupling.
Definition: ForcesForceField.cpp:130
VarData::Point2Shape
void Point2Shape(int iShape)
Point to the shape function.
Definition: VarDataEl.cpp:490
Forces::CutOff
double CutOff
Cut off of the interactions.
Definition: Forces.h:727
Forces::PrintForce
void PrintForce()
Print the force and the potential.
Definition: ForcesForceField.cpp:684
VarData::pWei2Par
double pWei2Par()
Parameter of the second order weighting function.
Definition: VarData.h:948
VarData::pNNano
int pNNano()
Number of nanoparticles.
Definition: VarDataComm.cpp:253
Vettore::VetV
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.
Definition: MatematicaVect.cpp:310
NANO::Coating
double Coating
Thickness of the LJ well.
Definition: VarData.h:461
Forces::SumSomeForces
int SumSomeForces(int HowMany)
Boh.
Definition: ForcesForceField.cpp:245
SPLINE::a2
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Definition: MatematicaStruct.h:110
MatInt::Coeff
double Coeff(int t1, int t2)
Prefactor of the force.
Definition: VarData.cpp:246
Forces::ForceFieldBulk
int ForceFieldBulk()
Armonic potential on a lattice.
Definition: ForcesForceField.cpp:187
Matematica::Spline3End
SPLINE Spline3End(double *P1, double *P2, int x, int y)
Three order spline last boundary.
Definition: MatematicaInterp.cpp:52
Forces::StepPot
double StepPot(double Dist2, int t1, int t2, double *Pot)
Step potential.
Definition: ForcesForceField.cpp:582
Forces::Potential
double Potential(double Dist, int t1, int t2, double *Pot)
Pointer to a potential.
Definition: Forces.h:504
Forces::Wei3
double Wei3(const double r, const double a)
Cubic weighting function.
Definition: ForcesForceField.cpp:718
Forces::RigidDistanceRad
double RigidDistanceRad(int n, int nn, double *dr)
Distance between two bodies.
Definition: ForcesForceField.cpp:428
NANO::PotThr
double PotThr
Maximum value of the potential.
Definition: VarData.h:469
Forces::ClearDens
void ClearDens()
Set the local densities to zero.
Definition: ForcesCalcEnergies.cpp:401
Forces::RigidDistanceAxis
double RigidDistanceAxis(int n, int nn, double *dr)
Distance between two bodies.
Definition: ForcesForceField.cpp:313
KFORCES::LJ
double LJ
Lennard-Jones.
Definition: Forces.h:155
NANO::Force
double Force[3]
Forces.
Definition: VarData.h:433
VarData::Mat
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
VarData::pNType
int pNType()
of types of the particle
Definition: VarDataComm.cpp:251
Matematica::Parab
SPLINE Parab(double *P1, double *P2, double *P3, int x, int y)
Three points parabolic interpolation.
Definition: MatematicaInterp.cpp:178
NANO::AMom
double AMom[3]
Angular momentum.
Definition: VarData.h:437
VarData::pPos
double pPos(int p, int d)
Return back folded position.
Definition: VarDataComm.cpp:192
NANO::CutOff
double CutOff
Cut off of the potential.
Definition: VarData.h:459
VarData::pTemp
double pTemp()
Temperature.
Definition: VarData.h:912
Forces::DerWei3
double DerWei3(const double r, const double a)
Derivative of the cubic weighting function.
Definition: ForcesForceField.cpp:722
VarData::pNChain
int pNChain()
Number of chain.
Definition: VarDataComm.cpp:247
NANO::AMomTemp
double AMomTemp[3]
Temporal angular momentum.
Definition: VarData.h:439
DdLinkedList::IfCurr
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
FORCES::Ext
double Ext[3]
External.
Definition: Forces.h:182
KFORCES::Elong
double Elong[3]
Elongation of the springs.
Definition: Forces.h:159
Matematica::Spline3
SPLINE Spline3(double *P1, double *P2, double *P3, int x, int y)
Three order spline.
Definition: MatematicaInterp.cpp:34
SPLINE
Coefficient of a spline.
Definition: MatematicaStruct.h:104
KFORCES::CutOff2
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
Matematica::Forth
SPLINE Forth(double *PA, double *PB, double *PC, double *PD, double *PE, int x, int y)
Five points four order interpolation.
Definition: MatematicaInterp.cpp:262
DomDecBasics::p2Curr
int p2Curr
Current particle.
Definition: Cubo.h:51
Matematica::Cubica
SPLINE Cubica(double *PA, double *PB, double *PC, double *PD, int x, int y)
Four point cubic interpolation.
Definition: MatematicaInterp.cpp:239
Forces::SumDens
double SumDens(int pInit, int pEnd)
Sum the local density for the particles between pInit and pEnd and multiply the factors by the virial...
Definition: ForcesCalcEnergies.cpp:243
KFORCES::BaseLine
double BaseLine
Baseline of the potential.
Definition: Forces.h:167
VarData::Pm
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
Forces::RigidHamaker
double RigidHamaker(int n, double Dist, double *Pot, double Sign)
Hamaker potential.
Definition: ForcesForceField.cpp:458
Forces::AddDens
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd.
Definition: ForcesCalcEnergies.cpp:218
Forces::LJ39
double LJ39(double Dist2, int t1, int t2, double *Pot)
Integrated Lennard Jones potential.
Definition: ForcesForceField.cpp:551
DdLinkedList::SetCurr
void SetCurr(int p)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:453
Forces::Wei2
double Wei2(const double r, const double b)
Quadratic weighting function.
Definition: ForcesForceField.cpp:726
Forces::Pc
DomDec * Pc
Pair list.
Definition: Forces.h:785
KFORCES::SLap
double SLap
Prefactor of the square laplacian.
Definition: Forces.h:149
Forces::DefNanoForceParam
void DefNanoForceParam()
Define the parameters for calculating the force.
Definition: ForcesForceField.cpp:405
KFORCES::ForThr
double ForThr
Maximum force allowed.
Definition: Forces.h:171
Forces::Bead2Move
int Bead2Move
Bead to move.
Definition: Forces.h:737
Forces::Fm
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
Vettore::PerpTo
void PerpTo(Vettore *o)
The vector perpendicolar.
Definition: MatematicaVect.cpp:380
VarData::pNPart
int pNPart()
Number of particle.
Definition: VarDataComm.cpp:246
Forces::DefForceParam
void DefForceParam()
Define the parameters for calculating the force.
Definition: ForcesForceField.cpp:388
Generated by   doxygen 1.8.11