Allink  v0.1
ForcesCalcEnergies.cpp
1 #include "Forces.h"
2 //-----------------------ENERGY-BONDED-----------------------
3 double Forces::CalcSpring(int p1){
4  double Nrg = 0.;
5  double Dist[4];
6  int NLink = 0;
7  int Link[2] = {-1,-1};
8  if(p1 > 0){
9  if( Pm[p1].CId == Pm[p1-1].CId ){
10  Link[NLink++] = p1 - 1;
11  }
12  }
13  if(p1 <= pNPart()){
14  if( Pm[p1].CId == Pm[p1+1].CId ){
15  Link[NLink++] = p1 + 1;
16  }
17  }
18  for(int l=0;l<NLink;l++){
19  TwoPartDist(p1,Link[l],Dist);
20  Nrg += .5*pkSpr()*SQR(Dist[3] - pSprRest());
21  }
22  return Nrg;
23 }
29 double Forces::CalcBonded(int p,double *Pot){
30  double NrgSpr = 0.;
31  double NrgBen = 0.;
32  int pA = p-1;
33  int pB = p;
34  int pC = p+1;
35  int IfBefore = 1;
36  int IfAfter = 1;
37  if(p > 0){
38  if(Pm[p].CId -1 == Pm[p-1].CId){
39  pA = p;
40  pB = p+1;
41  pC = p+2;
42  IfBefore = 0;
43  }
44  }
45  if(p <= pNPart() ){
46  if(Pm[p].CId + 1 == Pm[p+1].CId || p == pNPart()-1){
47  pC = p;
48  pB = p-1;
49  pA = p-2;
50  IfAfter = 0;
51  }
52  }
53  double DistBA[4];
54  double DistCB[4];
55  double CosAngle = 0.;
56  for(int d=0;d<3;d++){
57  DistBA[d] = Pm[pB].Pos[d] - Pm[pA].Pos[d];
58  DistCB[d] = Pm[pC].Pos[d] - Pm[pB].Pos[d];
59  }
60  DistBA[3] = (SQR(DistBA[0])+SQR(DistBA[1])+SQR(DistBA[2]));
61  DistCB[3] = (SQR(DistCB[0])+SQR(DistCB[1])+SQR(DistCB[2]));
62  for(int d=0;d<3;d++)
63  CosAngle += DistBA[d]*DistCB[d];
64  CosAngle /= (DistBA[3]*DistCB[3]);
65  if(IfAfter)
66  NrgSpr += .5*pkSpr()*SQR(DistBA[3] - pSprRest());
67  if(IfBefore)
68  NrgSpr += .5*pkSpr()*SQR(DistCB[3] - pSprRest());
69  NrgBen += pkBen()*(1.-CosAngle);
70  Pot[0] = NrgSpr;
71  Pot[1] = NrgBen;
72  return NrgSpr + NrgBen;
73 }
80 double Forces::CalcBendingGhost(double *Pos,int p){
81  if(pkBen() <= 0.) return 0.;
82  double NrgBen = 0.;
83  int pA = p;
84  int pB = p-1;
85  int pC = p-2;
86  if(pC < pChain(p)*pNPCh()){
87  return 0.;
88  }
89  double DistBA[3];
90  double DistCB[3];
91  double DistBA2 = 0.;
92  double DistCB2 = 0.;
93  double CosAngle = 0.;
94  for(int d=0;d<3;d++){
95  DistBA[d] = pPos(pB,d) - pPos(pA,d);
96  DistBA[d] -= floor(DistBA[d]*pInvEdge(d) + .5)*pEdge(d);
97  DistCB[d] = pPos(pC,d) - pPos(pB,d);
98  DistCB[d] -= floor(DistCB[d]*pInvEdge(d) + .5)*pEdge(d);
99  DistBA2 += SQR(DistBA[d]);
100  DistCB2 += SQR(DistCB[d]);
101  CosAngle += DistBA[d]*DistCB[d];
102  }
103  DistCB2 = sqrt(DistCB2);
104  DistBA2 = sqrt(DistBA2);
105  CosAngle /= (DistBA2*DistCB2);
106  NrgBen += pkBen()*(1.-CosAngle);
107  return NrgBen;
108 }
109 double Forces::CalcBondedCh(int c,double *Pot){
110  double DistBA[3];
111  double DistCB[3];
112  double NrgBend = 0.;
113  double NrgSpr = 0.;
114  for(int p=c*pNPCh();p<(c+1)*pNPCh()-1;p++){
115  double DistBA2 = 0.;
116  double DistCB2 = 0.;
117  double CosAngle = 0.;
118  for(int d=0;d<3;d++){
119  DistCB[d] = pPos(p+1,d) + pPos(p,d);
120  DistCB[d] -= floor(DistCB[d]*pInvEdge(d) + .5)*pEdge(d);
121  DistCB2 += SQR(DistCB[d]);
122  }
123  DistCB2 = sqrt(DistCB2);
124  NrgSpr += .5*pkSpr()*SQR(DistCB2 - pSprRest());
125  if(p == c*pNPCh()) continue;
126  for(int d=0;d<3;d++){
127  DistBA[d] = pPos(p,d) - pPos(p-1,d);
128  DistBA[d] -= floor(DistBA[d]*pInvEdge(d) + .5)*pEdge(d);
129  DistBA2 += SQR(DistBA[d]);
130  CosAngle += DistBA[d]*DistCB[d];
131  }
132  DistBA2 = sqrt(DistBA2);
133  CosAngle /= (DistBA2*DistCB2);
134  NrgBend += pkBen()*(1.-CosAngle);
135  }
136  Pot[0] = NrgSpr;
137  Pot[1] = NrgBend;
138  return NrgBend + NrgSpr;
139 }
140 double Forces::CalcBending(int p1){
141  if(Ln[p1].NLink < 1) return 0.;
142  if( !(p1%pNPCh()) ) return 0.;
143  int l1 = p1+1;
144  int l2 = p1-1;
145  Vettore v1(pPos(l1,0) - pPos(p1,0),pPos(l1,1) - pPos(p1,1),pPos(l1,2) - pPos(p1,2));
146  Vettore v2(pPos(p1,0) - pPos(l2,0),pPos(p1,1) - pPos(l2,1),pPos(p1,2) - pPos(l2,2));
147  return pkBen()*(1. - v1.CosAngle(&v2));
148 }
149 //%%%%%%%%%%%%%%%%%%%%%%%%%MC%NON%BONDED%%%%%%%%%%%%%%%%%%%%%
153 double Forces::NrgStep(int p1){
154  double DistRel[4];
155  double Nrg = 0.;
156  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
157  int p2 = Pc->p2Curr;
158  Pc->Dist2Curr(DistRel);
159  int t2 = pType(p2);
160  if(DistRel[3] > Kf.CutOff2) continue;
161  if(Pm[p1].Typ == Pm[p2].Typ)
162  Nrg += Kf.LJ;
163  else
164  Nrg += -Kf.LJ;
165  }
166  return Nrg;
167 }
168 double Forces::NrgStepCh(int c,double *Pot){
169  Pot[0] = 0.;Pot[1] = 0.;Pot[2] = 0.;
170  int p1 = c*pNPCh();
171  double DistRel[4];
172  double Pos[3];
173  //loop in the chain
174  for(int p=p1;p<p1+pNPCh();p++){
175  int t1 = pType(p);
176  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
177  int p2 = Pc->p2Curr;
178  int t2 = pType(p2);
179  // avoid self interactions in the chains
180  if(p2 > p1 && p2 < p1 + pNPCh() && p2 < p) continue;
181  Pc->Dist2Curr(DistRel);
182  if(DistRel[3] > Kf.CutOff2) continue;
183  // if(DistRel[3] < .1) Pot[2] += 55.;
184  // else Pot[2] += -.1;
185  Pot[2] += 1.;
186  }
187  }
188  return Pot[2];
189 }
190 //--------------------------Dens----------------------
191 int Forces::RemDens(int pInit,int pEnd){
192  double DistRel[4] = {0.,0.,0.,0.};
193  int NCutOff = 0;
194  double Count = 0.;
195  double Pos[3];
196  for(int p1=pInit;p1<pEnd;p1++){
197  int t1 = pType(p1);
198  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
199  int p2 = Pc->p2Curr;
200  int t2 = pType(p2);
201  // avoid self interactions in the chains
202  if(p2>=pInit && p2<pEnd && p2<=p1) continue;
203  Pc->Dist2Curr(DistRel);
204  if(DistRel[3] > Kf.CutOff2) continue;
205  double Dist = sqrt(DistRel[3]);
206  double w2 = Wei2(Dist,pWei2Par());
207  double w3 = Wei3(Dist,pWei3Par());
208  Dens2[p2*pNType()+t1] -= w2;
209  Dens2[p1*pNType()+t2] -= w2;
210  Dens3[p2*pNType()+t1] -= w3;
211  Dens3[p1*pNType()+t2] -= w3;
212  NCutOff++;
213  }
214  }
215  return NCutOff;
216 }
217 //Bottleneck!!!
218 int Forces::AddDens(int pInit,int pEnd){
219  double DistRel[4] = {0.,0.,0.,0.};
220  int NCutOff = 0;
221  double Pos[3];
222  for(int p1=pInit;p1<pEnd;p1++){
223  int t1 = pType(p1);
224  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
225  int p2 = Pc->p2Curr;
226  int t2 = pType(p2);
227  // avoid self interactions in the chains
228  if(p2>=pInit && p2<pEnd && p2<=p1)continue;
229  Pc->Dist2Curr(DistRel);
230  if(DistRel[3] > Kf.CutOff2) continue;
231  double Dist = sqrt(DistRel[3]);
232  double w2 = Wei2(Dist,pWei2Par());
233  double w3 = Wei3(Dist,pWei3Par());
234  Dens2[p2*pNType()+t1] += w2;
235  Dens2[p1*pNType()+t2] += w2;
236  Dens3[p2*pNType()+t1] += w3;
237  Dens3[p1*pNType()+t2] += w3;
238  NCutOff++;
239  }
240  }
241  return NCutOff;
242 }
243 double Forces::SumDens(int pInit,int pEnd){
244  double Nrg = 0.;
245  double OneThird = 1./3.;
246  for(int p=pInit;p<pEnd;p++){
247  int t1 = pType(p);
248  for(int t2=0;t2<pNType();t2++){
249  Nrg += .5*Dens2[p*pNType()+t2]*MInt->Coeff(t1,t2);
250  for(int t3=0;t3<pNType();t3++){
251  Nrg += Dens3[p*pNType()+t2]*Dens3[p*pNType()+t3]*OneThird*MInt->Coeff(t1,t2,t3);
252  }
253  }
254  }
255  return Nrg;
256 }
257 double Forces::DensFuncNrgGhost(double *Pos,int p1,int t1){
258  double DistRel[4] = {0.,0.,0.,0.};
259  double Nrg = 0.;
260  int NCutOff = 0;
261  double OneThird = 1./3.;
262  memset(LocDens2,0,pNType()*sizeof(double));
263  memset(LocDens3,0,pNType()*sizeof(double));
264  for(Pc->SetCurrGhost(Pos);Pc->IfCurrGhost();Pc->NextCurrGhost()){
265  int p2 = Pc->p2Curr;
266  int t2 = pType(p2);
267  if(p1 == p2) continue;
268  Pc->Dist2CurrGhost(DistRel);
269  if(DistRel[3] > Kf.CutOff2) continue;
270  double Dist = sqrt(DistRel[3]);
271  double w2 = Wei2(Dist,pWei2Par());
272  double w3 = Wei3(Dist,pWei3Par());
273  LocDens2[t2] += w2*2.;
274  LocDens3[t2] += w3;
275  double W3Add[3] = {0.,0.,0.};
276  W3Add[t1] = w3;
277  for(int t3=0;t3<pNType();t3++){
278  for(int t4=0;t4<pNType();t4++){
279  double Fact = Dens3[p2*pNType()+t3]*W3Add[t4];
280  Fact += W3Add[t3]*Dens3[p2*pNType()+t4];
281  Fact += W3Add[t3]*W3Add[t4];
282  //Why t2?
283  Nrg += Fact*OneThird*MInt->Coeff(t2,t3,t4);
284  }
285  }
286  NCutOff++;
287  }
288  for(int t2=0;t2<pNType();t2++){
289  Nrg += .5*LocDens2[t2]*MInt->Coeff(t1,t2);
290  for(int t3=0;t3<pNType();t3++){
291  Nrg += LocDens3[t2]*LocDens3[t3]*OneThird*MInt->Coeff(t1,t2,t3);
292  }
293  }
294  return Nrg;
295 }
296 double Forces::DensFuncNrgChInternal(int c){
297  double DistRel[4] = {0.,0.,0.,0.};
298  double Nrg = 0.;
299  int NCutOff = 0;
300  double OneThird = 1./3.;
301  memset(LocDens2,0,pNPCh()*pNType()*sizeof(double));
302  memset(LocDens3,0,pNPCh()*pNType()*sizeof(double));
303  // for(int pt=0;pt<pNPCh()*pNType();pt++){
304  // LocDens2[pt] = 0.;
305  // LocDens3[pt] = 0.;
306  // }
307  int pInit = Ch[c].InitBead;
308  int pEnd = Ch[c].EndBead;
309  double Pos[3];
310  for(int p1=pInit,pc=0;p1<pEnd;p1++,pc++){
311  int t1 = pType(p1);
312  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
313  int p2 = Pc->p2Curr;
314  int t2 = pType(p2);
315  if(p1 == p2) continue;
316  Pc->Dist2Curr(DistRel);
317  if(DistRel[3] > Kf.CutOff2) continue;
318  double Dist = sqrt(DistRel[3]);
319  double w2 = Wei2(Dist,pWei2Par());
320  double w3 = Wei3(Dist,pWei3Par());
321  LocDens2[pc*pNType()+t2] += w2;
322  LocDens3[pc*pNType()+t2] += w3;
323  //energy change for the neighbouring chains
324  if(p2 < p1 || p2 >= pEnd){
325  LocDens2[pc*pNType()+t2] += w2;
326  double W3Add[3] = {0.,0.,0.};
327  W3Add[t1] = w3;
328  //2FIX: multiple change in the 3 order density
329  for(int t3=0;t3<pNType();t3++){
330  for(int t4=0;t4<pNType();t4++){
331  double Fact = Dens3[p2*pNType()+t3]*W3Add[t4];
332  Fact += W3Add[t3]*Dens3[p2*pNType()+t4];
333  Fact += W3Add[t3]*W3Add[t4];
334  //Why t2?
335  Nrg += Fact*OneThird*MInt->Coeff(t2,t3,t4);
336  }
337  }
338  NCutOff++;
339  }
340  }
341  }
342  for(int pt=0;pt<pNPCh();pt++){
343  int t1 = pType(pInit+pt);
344  for(int t2=0;t2<pNType();t2++){
345  Nrg += .5*LocDens2[pt*pNType()+t2]*MInt->Coeff(t1,t2);
346  for(int t3=0;t3<pNType();t3++){
347  Nrg += LocDens3[pt*pNType()+t2]*LocDens3[pt*pNType()+t3]*OneThird*MInt->Coeff(t1,t2,t3);
348  }
349  }
350  }
351  return Nrg;
352 }
353 double Forces::DensFuncNrgGhostInternal(double *Pos,int p1,int t1){
354  double DistRel[4] = {0.,0.,0.,0.};
355  double Nrg = 0.;
356  int NCutOff = 0;
357  double OneThird = 1./3.;
358  memset(LocDens2,0,pNType()*sizeof(double));
359  memset(LocDens3,0,pNType()*sizeof(double));
360  int c = (int)(p1/pNPCh());
361  int pInit = Ch[c].InitBead;
362  int pEnd = Ch[c].EndBead;
363  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
364  int p2 = Pc->p2Curr;
365  int t2 = pType(p2);
366  if(p1 == p2) continue;
367  if(p2 < pInit) continue;
368  if(p2 >= pEnd) continue;
369  Pc->Dist2Curr(DistRel);
370  if(DistRel[3] > Kf.CutOff2) continue;
371  double Dist = sqrt(DistRel[3]);
372  double w2 = Wei2(Dist,pWei2Par());
373  double w3 = Wei3(Dist,pWei3Par());
374  LocDens2[t2] += w2*2.;
375  LocDens3[t2] += w3;
376  double W3Add[3] = {0.,0.,0.};
377  W3Add[t1] = w3;
378  for(int t3=0;t3<pNType();t3++){
379  for(int t4=0;t4<pNType();t4++){
380  double Fact = Dens3[p2*pNType()+t3]*W3Add[t4];
381  Fact += W3Add[t3]*Dens3[p2*pNType()+t4];
382  Fact += W3Add[t3]*W3Add[t4];
383  //Why t2?
384  Nrg += Fact*OneThird*MInt->Coeff(t2,t3,t4);
385  }
386  }
387  NCutOff++;
388  }
389  for(int t2=0;t2<pNType();t2++){
390  Nrg += .5*LocDens2[t2]*MInt->Coeff(t1,t2);
391  for(int t3=0;t3<pNType();t3++){
392  Nrg += LocDens3[t2]*LocDens3[t3]*OneThird*MInt->Coeff(t1,t2,t3);
393  }
394  }
395  return Nrg;
396 }
397 double Forces::DensFuncNrgSys(){
398  CalcDens(0,pNPart());
399  return SumDens(0,pNPart());
400 }
401 void Forces::ClearDens(){
402  memset(Dens2,0,pNPart()*pNType()*sizeof(double));
403  memset(Dens3,0,pNPart()*pNType()*sizeof(double));
404 }
405 void Forces::CalcDens(int pInit,int pEnd){
406  ClearDens();
407  AddDens(pInit,pEnd);
408 }
409 //------------------Calc-Ch-Nrg----------------------------
411 double Forces::CalcTotNrgCh(){
412  Shout("Calculate chains energy");
413  double Pot[3] = {0.,0.,0.};
414  double Nrg = 0.;
415  for(int c=0;c<pNChain();c++){
416  /* sustract to the densities the values of the c chain */
417  int p1 = c*pNPCh();
418  CalcNrgCh(c,Pot);
419  OldNrgCh[c*3 ] = Pot[0];
420  OldNrgCh[c*3+1] = Pot[1];
421  OldNrgCh[c*3+2] = Pot[2];
422  Nrg += Pot[0] + Pot[1] + Pot[2];
423  //fprintf(StatFile1,"%d %d %lf\n",c,pNPCh(c),Pot[0]);
424  }
425  return Nrg;
426 }
428 double Forces::DensFuncNrgCh(int c,double *Pot){
429  int p1 = c*pNPCh();
430  Pot[2] = DensFuncNrgChInternal(c);
431  for(int p=0;p<pNPCh();p++) Pot[2] += NanoNrg(p+p1);
432  return Pot[2];
433 }
435 double Forces::NrgChBondDens(int c,double *Pot){
436  CalcBondedCh(c,Pot);
437  int p1 = c*pNPCh();
438  Pot[2] = DensFuncNrgChInternal(c);
439  return Pot[0] + Pot[1] + Pot[2];
440 }
441 double Forces::DensFuncNrgChAv(int Ch){
442  CalcDens(Ch*pNPCh(),(Ch+1)*pNPCh());
443  return SumDens(Ch*pNPCh(),(Ch+1)*pNPCh());
444 }
445 double Forces::CalcPairwiseCh(int c,double *Pot){
446  Pot[0] = 0.;Pot[1] = 0.;Pot[2] = 0.;
447  int p1 = c*pNPCh();
448  double DistRel[4] = {0.,0.,0.,0.};;
449  double Pot1[3] = {0.,0.,0.};
450  //loop in the chain
451  double Pos[3];
452  for(int p=p1;p<p1+pNPCh();p++){
453  int t1 = pType(p);
454  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
455  int p2 = Pc->p2Curr;
456  Pc->Dist2Curr(DistRel);
457  if(p1 >= p2) continue;
458  if(p2 > p1 && p2 < p1 + pNPCh() && p2 < p) continue;
459  int t2 = pType(p2);
460  // avoid self interactions in the chains
461  if(DistRel[3] > Kf.CutOff2) continue;
462  Potential(DistRel[3],pType(p1),pType(p2),Pot1);
463  Pot[2] += Pot1[0];
464  }
465  }
466  return Pot[0] + Pot[1] + Pot[2];
467 }
468 //------------------Calc-Part-Nrg----------------------------
470 double Forces::CalcTotNrgBead(){
471  Shout("Calculate chains energy");
472  double Pot[3] = {0.,0.,0.};
473  double Nrg = 0.;
474  for(int p=0;p<pNChain();p++){
475  /* sustract to the densities the values of the c chain */
476  CalcNrgBead(p,Pot);
477  OldNrgBead[p*3 ] = Pot[0];
478  OldNrgBead[p*3+1] = Pot[1];
479  OldNrgBead[p*3+2] = Pot[2];
480  Nrg += Pot[0] + Pot[1] + Pot[2];
481  }
482  return Nrg;
483 }
484 double Forces::DensFuncNrgBead(int p1){
485  double Pos[3];
486  pPos(p1,Pos);
487  return DensFuncNrgGhost(Pos,p1,pType(p1));
488  RemDens(p1,p1+1);
489  double Nrg = OldNrgSys - SumDens(0,pNPart());
490  AddDens(p1,p1+1);
491  return Nrg;
492 }
494 double Forces::CalcNrgBeadDensFunc(int p,double *Pot){
495  CalcBonded(p,Pot);
496  Pot[2] = DensFuncNrgBead(p);
497  return Pot[0] + Pot[1] + Pot[2];
498 }
499 double Forces::CalcPairwise(int p1,double *Pot){
500  //return CheckDomDec(p1);
501  double NrgSum = 0.;
502  int NCutOff = 0;
503  double Pos[3];
504  double DistRel[4] = {0.,0.,0.,0.};
505  Pot[0] = 0.;Pot[1] = 0.;Pot[2] = 0.;
506  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
507  int p2 = Pc->p2Curr;
508  // for(int p2 = 0;p2<pNPart();p2++){
509  if(p2 == p1) continue;
510  //if(p2 >= p1) continue;
511  // TwoPartDist(p1,p2,DistRel);
512  Pc->Dist2Curr(DistRel);
513  if(DistRel[3] > Kf.CutOff2) continue;
514  Potential(DistRel[3],pType(p1),pType(p2),Pot);
515  NrgSum += Pot[0] + Pot[1] + Pot[2];
516  NCutOff++;
517  }
518  return NrgSum;
519 }
520 //----------------------MOLECULAR-DYNAMICS----------------------
521 void Forces::ChooseCalcMode(int Mode){
522  printf("Calculation mode: ");
523  if(VAR_IF_TYPE(Mode,CALC_PAIR)){
524  NrgBead = &Forces::CalcPairwise;
525  NrgCh = &Forces::CalcPairwiseCh;
526  printf("Pairwise interactions\n");
527  }
528  else if(VAR_IF_TYPE(Mode,CALC_DENS)){
529  NrgBead = &Forces::CalcNrgBeadDensFunc;
530  NrgCh = &Forces::DensFuncNrgCh;
531  Kf.CutOff2 = 1.;
532  printf("Lipid model particle energy\n");
533  }
534  else if(VAR_IF_TYPE(Mode,CALC_DENS_CH)){
535  NrgBead = &Forces::CalcNrgBeadDensFunc;
536  NrgCh = &Forces::DensFuncNrgCh;
537  Kf.CutOff2 = 1.;
538  printf("Lipid model chain energy\n");
539  }
540  else{
541  printf("mode not present\n");
542  exit(1);
543  }
544 }
545 void Forces::ChoosePot(int Mode){
546  printf("Potential: ");
547  if(VAR_IF_TYPE(Mode,CALC_LJ)){
548  CalcPot = &Forces::LJPot;
549  printf("Lennard Jones\n");
550  }
551  else if(VAR_IF_TYPE(Mode,CALC_LJ39)){
552  CalcPot = &Forces::LJ39;
553  printf("Lennard Jones 3-9\n");
554  }
555  else if(VAR_IF_TYPE(Mode,CALC_HARM)){
556  CalcPot = &Forces::Harmonic;
557  printf("Harmonic\n");
558  }
559  else if(VAR_IF_TYPE(Mode,CALC_STEP)){
560  CalcPot = &Forces::StepPot;
561  printf("Step\n");
562  }
563  else if(VAR_IF_TYPE(Mode,CALC_ELECTRO)){
564  CalcPot = &Forces::ElectroPot;
565  printf("Step\n");
566  }
567  else{
568  printf("not present\n");
569  exit(1);
570  }
571  DefForceParam();
572 }
573 double Forces::SumForcesMD(){
574  double Pot[3];
575  double DistRel[4] = {0.,0.,0.,0.};
576  double Nrg = 0.;
577  for(int p1=0;p1<pNPart();p1++){
578  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
579  int p2 = Pc->p2Curr;
580  if(p1 >= p2) continue;
581  Pc->Dist2Curr(DistRel);
582  if(DistRel[3] > Kf.CutOff2) continue;
583  double InvDist = 1./sqrt(DistRel[3]);
584  double Cons = Potential(DistRel[3],pType(p1),pType(p2),Pot);
585  for(int d=0;d<3;d++){
586  Fm[p1].Dir[d] += Cons*DistRel[d]*InvDist;
587  Fm[p2].Dir[d] -= Cons*DistRel[d]*InvDist;
588  }
589  Nrg += Pot[0];
590  }
591  }
592  return Nrg;
593 }
594 #include <sys/time.h>
595 void Forces::CheckPairList(){
596  int NAll = pNPart();
597  int *PDom = (int *)calloc(NAll*NAll,sizeof(int));
598  if(PDom == NULL){printf("Could not alloc PDom\n");return;};
599  int *PLoop = (int *)calloc(NAll*NAll,sizeof(int));
600  if(PLoop == NULL){printf("Could not alloc PLoop\n");return;};
601  int NDom = 0;
602  int NLoop = 0;
603  double NPairDom = 0.;
604  double NPairLoop = 0.;
605  double Dist[3];
606  double Nrg = 0.;
607  double Pot[2];
608  double DistRel[4] = {0.,0.,0.,0.};
609  double TimeDomDec;
610  double TimeLoop;
611  timespec TimeInit;
612  timespec TimeEnd;
613 // clock_gettime(CLOCK_REALTIME, &TimeInit);
614  double Pos[3];
615  for(int p1=0;p1<pNPart();p1++){
616  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
617  int p2 = Pc->p2Curr;
618  if(p2 <= p1) continue;
619  Pc->Dist2Curr(DistRel);
620  NPairDom += 1.;
621  if(DistRel[3] > Kf.CutOff2) continue;
622  PLoop[NDom*2+0] = p2;
623  PLoop[NDom*2+1] = p2;
624  NDom++;
625  }
626  }
627 // clock_gettime(CLOCK_REALTIME, &TimeEnd);
628  TimeDomDec = (double)(TimeInit.tv_nsec - TimeEnd.tv_nsec);
629 // clock_gettime(CLOCK_REALTIME, &TimeInit);
630  for(int p1=0;p1<pNPart();p1++){
631  for(int p2=p1+1;p2<pNPart();p2++){
632  double Dist2 = 0.;
633  for(int d=0;d<3;d++){
634  Dist[d] = pPos(p1,d) - pPos(p2,d);
635  Dist[d] -= floor(Dist[d]*pInvEdge(d) + .5)*pEdge(d);
636  Dist2 += SQR(Dist[d]);
637  }
638  NPairLoop += 1.;
639  if(Dist2 > Kf.CutOff2) continue;
640  PLoop[NLoop*2+0] = p1;
641  PLoop[NLoop*2+1] = p2;
642  NLoop++;
643  }
644  }
645 // clock_gettime(CLOCK_REALTIME, &TimeEnd);
646  TimeLoop = (double)(TimeInit.tv_nsec - TimeEnd.tv_nsec);
647  //Mat->Sort(PDom,NDom);
648  printf("[Pairs] %d=%d Gain (pair): %lf (time): %lf\n",NDom,NLoop,NPairLoop/NPairDom,TimeLoop/TimeDomDec);
649  for(int p=0;p<MAX(NDom,NLoop);p++){
650  //printf("%d %d %d\n",p,PDom[p],PLoop[p]);
651  //if(p >= NDom)
652  //printf(" %d %d\n",PLoop[p],Pc->pCell(PLoop[p]));
653  }
654  free(PDom);
655  free(PLoop);
656 }
658 double Forces::CheckDomDec(int p1){
659  const int NAll = pNPart();
660  int *PDom = new int[NAll];
661  int *PCell = new int[NAll];
662  int NDom = 0;
663  int *PLoop = new int[NAll];
664  int NLoop = 0;
665  double NCell = 0.;
666  double Pos[3];
667  double DistRel[4] = {0.,0.,0.,0.};
668  double Nrg=0.;
669  double Pot[2];
670  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
671  int p2 = Pc->p2Curr;
672  if(p2 == p1) continue;
673  Pc->Dist2Curr(DistRel);
674  if(DistRel[3] > Kf.CutOff2) continue;
675  PDom[NDom++] = p2;
676  PCell[NDom] = Pc->cCurr;
677  }
678  for(int p2=0;p2<pNPart();p2++){
679  if(p1 == p2) continue;
680  double Dist2 = 0.;
681  for(int d=0;d<3;d++){
682  DistRel[d] = pPos(p1,d) - pPos(p2,d);
683  DistRel[d] -= floor(DistRel[d]*pInvEdge(d) + .5)*pEdge(d);
684  Dist2 += SQR(DistRel[d]);
685  }
686  if(Dist2 > Kf.CutOff2) continue;
687  PLoop[NLoop++] = p2;
688  }
689  //Mat->Sort(PDom,NDom);
690  printf("-------------%d=%d Ratio %lf\n",NDom,NLoop,pNPart()/NCell);
691  for(int p=0;p<MAX(NDom,NLoop);p++){
692  printf("%d %d %d %d\n",p,PDom[p],PLoop[p],PCell[p]);
693  //if(p >= NDom)
694  //printf(" %d %d\n",PLoop[p],Pc->pCell(PLoop[p]));
695  }
696  printf("nei \n");
697  delete [] PDom;
698  delete [] PCell;
699  delete [] PLoop;
700  return 0.;
701 }
Forces::DensFuncNrgBead
double DensFuncNrgBead(int p1)
Calculation of the energy from the density functional Hamiltonian for the particle p1...
Definition: ForcesCalcEnergies.cpp:484
VarData::Ch
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
VarData::pkBen
double pkBen()
Bending coupling.
Definition: VarData.h:932
Forces::CheckDomDec
double CheckDomDec(int p)
Check if all the particles are taken in account.
Definition: ForcesCalcEnergies.cpp:658
Forces::CalcTotNrgCh
double CalcTotNrgCh()
Calculate and sum up the energy of the chains.
Definition: ForcesCalcEnergies.cpp:411
Forces::LJPot
double LJPot(double Dist2, int t1, int t2, double *Pot)
Classical Lennard Jones potential.
Definition: ForcesForceField.cpp:567
DomDecBasics::cCurr
int cCurr
Cell where the current particle sits.
Definition: Cubo.h:41
Forces::NrgStepCh
double NrgStepCh(int c, double *Pot)
The energy per chain is constant within the cutoff.
Definition: ForcesCalcEnergies.cpp:168
VarData::Ln
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
Forces::CalcPairwiseCh
double CalcPairwiseCh(int c, double *Pot)
Calculate the non bonded interaction energy with the neighbouring particles.
Definition: ForcesCalcEnergies.cpp:445
Forces::CalcBendingGhost
double CalcBendingGhost(double *Pos, int pExt)
Calculate the bending energy for a ghost particle.
Definition: ForcesCalcEnergies.cpp:80
DdLinkedList::Dist2CurrGhost
void Dist2CurrGhost(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:554
VarData::pChain
int pChain(int p)
Return the chain.
Definition: VarDataComm.cpp:222
Vettore
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
Forces::CalcNrgBeadDensFunc
void CalcNrgBeadDensFunc()
Calculate the spring, bending and non bonded interactions and write it in OldNrgPm.
Forces::CheckPairList
void CheckPairList()
Check the pair list.
Definition: ForcesCalcEnergies.cpp:595
Forces::CalcBending
double CalcBending(int p)
Calculate the bonded interaction energy with the neighbouring particles.
Definition: ForcesCalcEnergies.cpp:140
Forces::DensFuncNrgCh
double DensFuncNrgCh(int c, double *Pot)
Calculate the energy from the density functional Hamiltonian for the chain c.
Definition: ForcesCalcEnergies.cpp:428
Forces::CalcTotNrgBead
double CalcTotNrgBead()
Calculate and sum up the energy of the part.
Definition: ForcesCalcEnergies.cpp:470
Forces::Harmonic
double Harmonic(double Dist2, int t1, int t2, double *Pot)
Harmonic potential.
Definition: ForcesForceField.cpp:537
VarData::pInvEdge
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
VarData::pType
int pType(int p)
Return the type.
Definition: VarDataComm.cpp:219
PART::Pos
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
Forces::NrgChBondDens
double NrgChBondDens(int c, double *Pot)
Calculate the bond and the density functional energies.
Definition: ForcesCalcEnergies.cpp:435
Forces::SumForcesMD
double SumForcesMD()
Iterate all over the particles and calculate the forces.
Definition: ForcesCalcEnergies.cpp:573
FORCES::Dir
double Dir[3]
Direction.
Definition: Forces.h:178
Forces::RemDens
int RemDens(int pInit, int pEnd)
Substract the densities connected with the particles between pInit and pEnd.
Definition: ForcesCalcEnergies.cpp:191
VarData::pkSpr
double pkSpr()
Spring coupling.
Definition: VarData.h:934
VarData::pSprRest
double pSprRest()
Rest distance of the harmonic potential.
Definition: VarData.h:942
DdLinkedList::Dist2Curr
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
Vettore::CosAngle
double CosAngle(Vettore *u)
Computes the cosine with respect to.
Definition: MatematicaVect.cpp:229
Forces::CalcSpring
double CalcSpring(int p)
Calculate the spring interaction energy with the neighbouring particles.
Definition: ForcesCalcEnergies.cpp:3
CHAIN::InitBead
int InitBead
Initial bead.
Definition: VarData.h:250
Forces::Shout
void Shout(const char *s,...)
Internal message.
Definition: Forces.cpp:3
Forces::ElectroPot
double ElectroPot(double Dist2, int t1, int t2, double *Pot)
Potential for the electrical lines.
Definition: ForcesForceField.cpp:590
VarData::pEdge
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
VarData::pNPCh
int pNPCh()
Number of particle per chain.
Definition: VarDataComm.cpp:249
DdLinkedList::IfCurrGhost
int IfCurrGhost()
Tell when the curr loop is over.
Definition: Cubo.cpp:550
DdLinkedList::NextCurrGhost
void NextCurrGhost()
Increase the iterator to the next couple.
Definition: Cubo.cpp:536
Forces::NrgStep
double NrgStep(int p)
The energy is constant within the cutoff.
Definition: ForcesCalcEnergies.cpp:153
VarData::pWei3Par
double pWei3Par()
Parameter of the third order weighting function.
Definition: VarData.h:950
VarData::TwoPartDist
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
Forces::NanoNrg
double NanoNrg(int p)
Exchange energy with the nano.
Definition: ForcesForceField.cpp:366
VarData::MInt
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
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
VarData::pWei2Par
double pWei2Par()
Parameter of the second order weighting function.
Definition: VarData.h:948
Forces::DensFuncNrgSys
double DensFuncNrgSys()
Calculation of the energy from the density functional Hamiltonian for the system. ...
Definition: ForcesCalcEnergies.cpp:397
MatInt::Coeff
double Coeff(int t1, int t2)
Prefactor of the force.
Definition: VarData.cpp:246
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::ClearDens
void ClearDens()
Set the local densities to zero.
Definition: ForcesCalcEnergies.cpp:401
Forces::DensFuncNrgChAv
double DensFuncNrgChAv(int c)
Calculate the average energy from the density functional Hamiltonian for the chain c...
Definition: ForcesCalcEnergies.cpp:441
DdLinkedList::SetCurrGhost
void SetCurrGhost(double *Pos)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:512
Forces::CalcNrgCh
double CalcNrgCh(int c, double *Pot)
Pointer to the chain energy function.
Definition: Forces.h:513
KFORCES::LJ
double LJ
Lennard-Jones.
Definition: Forces.h:155
CHAIN::EndBead
int EndBead
End bead.
Definition: VarData.h:252
VarData::pNType
int pNType()
of types of the particle
Definition: VarDataComm.cpp:251
Forces::DensFuncNrgChInternal
double DensFuncNrgChInternal(int c)
Calculate the average energy from the density functional Hamiltonian for the chain c...
Definition: ForcesCalcEnergies.cpp:296
Forces::CalcBondedCh
double CalcBondedCh(int c, double *Pot)
Calculate the bonded and spring interaction in a cell.
Definition: ForcesCalcEnergies.cpp:109
VarData::pPos
double pPos(int p, int d)
Return back folded position.
Definition: VarDataComm.cpp:192
VarData::pNChain
int pNChain()
Number of chain.
Definition: VarDataComm.cpp:247
Forces::CalcPairwise
double CalcPairwise(int p, double *Pot)
Calculate the non bonded interaction energy with the neighbouring particles.
Definition: ForcesCalcEnergies.cpp:499
DdLinkedList::IfCurr
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
Forces::DensFuncNrgGhost
double DensFuncNrgGhost(double *Pos, int p1, int t1)
Calculation of the energy from the density functional Hamiltonian for the ghost particle.
Definition: ForcesCalcEnergies.cpp:257
KFORCES::CutOff2
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
DomDecBasics::p2Curr
int p2Curr
Current particle.
Definition: Cubo.h:51
Forces::CalcBonded
double CalcBonded(int p, double *Pot)
Calculate the spring and the bonded interactions with the other monomers in the chain.
Definition: ForcesCalcEnergies.cpp:29
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
Forces::CalcDens
void CalcDens()
Calculate the densities.
Definition: ForcesTens.cpp:45
VarData::Pm
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
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
Forces::DensFuncNrgGhostInternal
double DensFuncNrgGhostInternal(double *Pos, int p1, int t1)
Calculation of the energy from the density functional Hamiltonian for the ghost particle.
Definition: ForcesCalcEnergies.cpp:353
Forces::Fm
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
Forces::CalcNrgBead
double CalcNrgBead(int p, double *Pot)
Pointer to the energy function.
Definition: Forces.h:511
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