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 ){
11  }
12  }
13  if(p1 <= pNPart()){
14  if( Pm[p1].CId == Pm[p1+1].CId ){
16  }
17  }
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];
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 }
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];
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];
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 }
398  CalcDens(0,pNPart());
399  return SumDens(0,pNPart());
400 }
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();
408 }
409 //------------------Calc-Ch-Nrg----------------------------
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 }
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----------------------------
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 */
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());
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)){
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)){
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 }
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>
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 }
Calculation of the energy from the density functional Hamiltonian for the particle p1...
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
double pkBen()
Bending coupling.
Definition: VarData.h:932
double CheckDomDec(int p)
Check if all the particles are taken in account.
double CalcTotNrgCh()
Calculate and sum up the energy of the chains.
double LJPot(double Dist2, int t1, int t2, double *Pot)
Classical Lennard Jones potential.
int cCurr
Cell where the current particle sits.
Definition: Cubo.h:41
double NrgStepCh(int c, double *Pot)
The energy per chain is constant within the cutoff.
Array of linking between the particles.
Definition: VarData.h:1048
double CalcPairwiseCh(int c, double *Pot)
Calculate the non bonded interaction energy with the neighbouring particles.
double CalcBendingGhost(double *Pos, int pExt)
Calculate the bending energy for a ghost particle.
void Dist2CurrGhost(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:554
int pChain(int p)
Return the chain.
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
Calculate the spring, bending and non bonded interactions and write it in OldNrgPm.
void CheckPairList()
Check the pair list.
double CalcBending(int p)
Calculate the bonded interaction energy with the neighbouring particles.
double DensFuncNrgCh(int c, double *Pot)
Calculate the energy from the density functional Hamiltonian for the chain c.
Calculate and sum up the energy of the part.
double Harmonic(double Dist2, int t1, int t2, double *Pot)
Harmonic potential.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
int pType(int p)
Return the type.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double NrgChBondDens(int c, double *Pot)
Calculate the bond and the density functional energies.
double SumForcesMD()
Iterate all over the particles and calculate the forces.
double Dir[3]
Direction.
Definition: Forces.h:178
int RemDens(int pInit, int pEnd)
Substract the densities connected with the particles between pInit and pEnd.
double pkSpr()
Spring coupling.
Definition: VarData.h:934
double pSprRest()
Rest distance of the harmonic potential.
Definition: VarData.h:942
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
double CosAngle(Vettore *u)
Computes the cosine with respect to.
double CalcSpring(int p)
Calculate the spring interaction energy with the neighbouring particles.
Definition: VarData.h:250
void Shout(const char *s,...)
Internal message.
Definition: Forces.cpp:3
double ElectroPot(double Dist2, int t1, int t2, double *Pot)
Potential for the electrical lines.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int pNPCh()
Number of particle per chain.
int IfCurrGhost()
Tell when the curr loop is over.
Definition: Cubo.cpp:550
void NextCurrGhost()
Increase the iterator to the next couple.
Definition: Cubo.cpp:536
double NrgStep(int p)
The energy is constant within the cutoff.
double pWei3Par()
Parameter of the third order weighting function.
Definition: VarData.h:950
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
double NanoNrg(int p)
Exchange energy with the nano.
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
void NextCurr()
Increase the iterator to the next couple.
Definition: Cubo.cpp:479
double pWei2Par()
Parameter of the second order weighting function.
Definition: VarData.h:948
double DensFuncNrgSys()
Calculation of the energy from the density functional Hamiltonian for the system. ...
double Coeff(int t1, int t2)
Prefactor of the force.
Definition: VarData.cpp:246
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.
void ClearDens()
Set the local densities to zero.
double DensFuncNrgChAv(int c)
Calculate the average energy from the density functional Hamiltonian for the chain c...
void SetCurrGhost(double *Pos)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:512
double CalcNrgCh(int c, double *Pot)
Pointer to the chain energy function.
Definition: Forces.h:513
double LJ
Lennard-Jones.
Definition: Forces.h:155
Definition: VarData.h:252
int pNType()
of types of the particle
double DensFuncNrgChInternal(int c)
Calculate the average energy from the density functional Hamiltonian for the chain c...
double CalcBondedCh(int c, double *Pot)
Calculate the bonded and spring interaction in a cell.
double pPos(int p, int d)
Return back folded position.
int pNChain()
Number of chain.
double CalcPairwise(int p, double *Pot)
Calculate the non bonded interaction energy with the neighbouring particles.
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
double DensFuncNrgGhost(double *Pos, int p1, int t1)
Calculation of the energy from the density functional Hamiltonian for the ghost particle.
double CutOff2
CutOff of the lennard jones.
Definition: Forces.h:165
int p2Curr
Current particle.
Definition: Cubo.h:51
double CalcBonded(int p, double *Pot)
Calculate the spring and the bonded interactions with the other monomers in the chain.
double SumDens(int pInit, int pEnd)
Sum the local density for the particles between pInit and pEnd and multiply the factors by the virial...
void CalcDens()
Calculate the densities.
Definition: ForcesTens.cpp:45
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd.
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)
DomDec * Pc
Pair list.
Definition: Forces.h:785
double DensFuncNrgGhostInternal(double *Pos, int p1, int t1)
Calculation of the energy from the density functional Hamiltonian for the ghost particle.
FORCES * Fm
Array containing the forces for each particle.
Definition: Forces.h:781
double CalcNrgBead(int p, double *Pot)
Pointer to the energy function.
Definition: Forces.h:511
int pNPart()
Number of particle.
void DefForceParam()
Define the parameters for calculating the force.