Allink  v0.1
ForcesForceField.cpp
1 #include "Forces.h"
2 //-------------------------SPLINES--------------------------------
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 }
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-------------------
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 }
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-------------------------
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 }
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 }
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 }
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;
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 }
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 }
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 }
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 }
void Export(double *x)
Export.
SPLINE Spline4Beg(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline first boundary.
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double LJMin
Lennard-Jones.
Definition: Forces.h:161
double pkBen()
Bending coupling.
Definition: VarData.h:932
double Ext
External force.
Definition: Forces.h:153
double LJPot(double Dist2, int t1, int t2, double *Pot)
Classical Lennard Jones potential.
int ForceFieldLine()
Helfrich Hamiltonian for a line.
double NanoDist2(double *Pos, int n)
Pointer to a generic function.
Definition: VarData.h:764
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
double RigidCoulomb(int nNano, double Dist, double *Pot, double Sign)
Force between to rigid bodies.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
void TabForceAlloc(int NTabExt)
Tabulate the values of the force.
double Gamma
Friction term.
Definition: VarData.h:451
double Hamaker
Strength of the interaction.
Definition: VarData.h:447
double DistThr
Minimum distance threshold.
Definition: VarData.h:465
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
double Lap
Prefactor of the Laplacian.
Definition: Forces.h:147
double RigidLJ(int nNano, double Dist, double *Pot, double Sign)
Force between to rigid bodies.
SPLINE Spline3Beg(double *P1, double *P2, double *P3, int x, int y)
Three order spline first boundary.
Parabolas coefficients.
double Rad
Size.
Definition: VarData.h:445
double Harmonic(double Dist2, int t1, int t2, double *Pot)
Harmonic potential.
double Vel[3]
Velocity.
Definition: VarData.h:431
int pType(int p)
Return the type.
double DerWei2(const double r, const double b)
Derivative of the quadratic weighting function.
double El[3]
Elastic force.
Definition: Forces.h:151
double PotThr
Maximum potential allowed.
Definition: Forces.h:169
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
int IfInterp
If interpolates with the splines.
Definition: Forces.h:755
double Dir[3]
Direction.
Definition: Forces.h:178
double pkSpr()
Spring coupling.
Definition: VarData.h:934
double pSprRest()
Rest distance of the harmonic potential.
Definition: VarData.h:942
void ForceFieldRod()
Bending potential on a rod.
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
int SysAlloc
Which arrays are allocated.
Definition: Forces.h:753
void CalcForcesDensFunc()
Calculate the forces for the density functional.
double Pos[3]
Position.
Definition: VarData.h:427
double ElectroPot(double Dist2, int t1, int t2, double *Pot)
Potential for the electrical lines.
void NanoInteraction()
Calculate all the interaction with the nano.
SPLINE Parab2(double *PA, double *PB, double *PC, int x, int y)
Three points parabolic interpolation.
int pNPCh()
Number of particle per chain.
void ForceFieldRigid()
Interaction between rigid bodies.
double DistThr
Correspondent radial distance for this force.
Definition: Forces.h:173
double Dx
Spatial separation between particles.
Definition: Forces.h:729
double pWei3Par()
Parameter of the third order weighting function.
Definition: VarData.h:950
SPLINE Spline4(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline.
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
int pNBlock()
Number of blocks.
double NanoNrg(int p)
Exchange energy with the nano.
double BaseLine
Baseline of the potential.
Definition: VarData.h:463
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
void TabPot()
Tabulate the values of the potential.
double Casuale()
Random uniform number.
double ProjOnAxis(Vettore *a)
Projects along the axis.
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
void NextCurr()
Increase the iterator to the next couple.
Definition: Cubo.cpp:479
double Zeta
Stochastic term.
Definition: VarData.h:453
double ForThr
Maximum value of the force.
Definition: VarData.h:467
int ForceFieldLeaves()
Helfrich Hamiltonian with an elastic coupling.
void Point2Shape(int iShape)
Point to the shape function.
Definition: VarDataEl.cpp:490
double CutOff
Cut off of the interactions.
Definition: Forces.h:727
void PrintForce()
Print the force and the potential.
double pWei2Par()
Parameter of the second order weighting function.
Definition: VarData.h:948
int pNNano()
Number of nanoparticles.
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.
double Coating
Thickness of the LJ well.
Definition: VarData.h:461
int SumSomeForces(int HowMany)
Boh.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double Coeff(int t1, int t2)
Prefactor of the force.
Definition: VarData.cpp:246
int ForceFieldBulk()
Armonic potential on a lattice.
SPLINE Spline3End(double *P1, double *P2, int x, int y)
Three order spline last boundary.
double StepPot(double Dist2, int t1, int t2, double *Pot)
Step potential.
double Potential(double Dist, int t1, int t2, double *Pot)
Pointer to a potential.
Definition: Forces.h:504
double Wei3(const double r, const double a)
Cubic weighting function.
double RigidDistanceRad(int n, int nn, double *dr)
Distance between two bodies.
double PotThr
Maximum value of the potential.
Definition: VarData.h:469
void ClearDens()
Set the local densities to zero.
double RigidDistanceAxis(int n, int nn, double *dr)
Distance between two bodies.
double LJ
Lennard-Jones.
Definition: Forces.h:155
double Force[3]
Forces.
Definition: VarData.h:433
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int pNType()
of types of the particle
SPLINE Parab(double *P1, double *P2, double *P3, int x, int y)
Three points parabolic interpolation.
double AMom[3]
Angular momentum.
Definition: VarData.h:437
double pPos(int p, int d)
Return back folded position.
double CutOff
Cut off of the potential.
Definition: VarData.h:459
double pTemp()
Temperature.
Definition: VarData.h:912
double DerWei3(const double r, const double a)
Derivative of the cubic weighting function.
int pNChain()
Number of chain.
double AMomTemp[3]
Temporal angular momentum.
Definition: VarData.h:439
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
double Ext[3]
External.
Definition: Forces.h:182
double Elong[3]
Elongation of the springs.
Definition: Forces.h:159
SPLINE Spline3(double *P1, double *P2, double *P3, int x, int y)
Three order spline.
Coefficient of a spline.
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
SPLINE Forth(double *PA, double *PB, double *PC, double *PD, double *PE, int x, int y)
Five points four order interpolation.
int p2Curr
Current particle.
Definition: Cubo.h:51
SPLINE Cubica(double *PA, double *PB, double *PC, double *PD, int x, int y)
Four point cubic interpolation.
double SumDens(int pInit, int pEnd)
Sum the local density for the particles between pInit and pEnd and multiply the factors by the virial...
double BaseLine
Baseline of the potential.
Definition: Forces.h:167
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
double RigidHamaker(int n, double Dist, double *Pot, double Sign)
Hamaker potential.
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd.
double LJ39(double Dist2, int t1, int t2, double *Pot)
Integrated Lennard Jones potential.
void SetCurr(int p)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:453
double Wei2(const double r, const double b)
Quadratic weighting function.
DomDec * Pc
Pair list.
Definition: Forces.h:785
double SLap
Prefactor of the square laplacian.
Definition: Forces.h:149
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
void PerpTo(Vettore *o)
The vector perpendicolar.
int pNPart()
Number of particle.
void DefForceParam()
Define the parameters for calculating the force.