Allink  v0.1
ForcesIntegration.cpp
1 #include "Forces.h"
6  AddRigid();
7  double *Known = (double *) calloc(pNPCh(),sizeof(double));
8  double *UnKnown = (double *) calloc(pNPCh(),sizeof(double));
9  for(int c=0;c<pNChain();c++){
10  int pInit = c*pNPCh();
11  int pEnd = (c+1)*pNPCh();
12  double Offset = 0.;
13  for(int p=pInit;p<pEnd;p++)
14  Known[p-pInit] = Pm[p].Pos[CNorm] - Offset;
15  IntMatrix->Apply(Known,UnKnown);
16  //IntMatrix->Solve(Known,UnKnown);
17  for(int p=pInit;p<pEnd;p++){
18  int p1 = p-pInit;
19  if(Pm[p].Typ != 0){continue;}
20  Pm[p].Pos[CNorm] = UnKnown[p1] + Offset;
21  }
22  }
23  free(Known);
24  free(UnKnown);
25 }
30  double *Known = (double *) calloc(pNPart(),sizeof(double));
31  double *UnKnown = (double *) calloc(pNPart(),sizeof(double));
32  double Offset = .5*pEdge(2);
33  for(int p=0;p<pNPart();p++)
34  Known[p] = Pm[p].Pos[CNorm] - Offset;
35  IntMatrix->Apply(Known,UnKnown);
36  //IntMatrix->Solve(Known,UnKnown);
37  for(int p=0;p<pNPart();p++){
38  if(Pm[p].Typ != 0){continue;}
39  Pm[p].Pos[CNorm] = UnKnown[p] + Offset;
40  }
41  free(Known);
42  free(UnKnown);
43 }
48  AddRigid();
49  //SmoothGrid(NEdge);
50  double Inter = pEdge(CLat1)/(double)NEdge;
51  SPLINE Weight;
52  int NIter = 2000;
53  Weight.a0 = Kf.El[2];
54  Weight.a1 = 0./Inter;
55  Weight.a2 = Kf.Lap/SQR(Inter);
56  Weight.a3 = 0./(Inter*SQR(Inter));
57  Weight.a4 = Kf.SLap/(SQR(Inter)*SQR(Inter));
58  int NDim = 2;
59  Matrice *CoeffMatrix = new Matrice(Weight,NDim);
60  double Offset = 0.;
61  double *Known = (double *) calloc(pNPart(),sizeof(double));
62  double *UnKnown = (double *) calloc(pNPart(),sizeof(double));
63  double Sol = 0.;
64  double Zed = 0.;
65  int Link[2][5];
66  CoeffMatrix->Print();
67  for(int i=0;i<NIter;i++){
68  for(int p=0;p<pNPart();p++)
69  Known[p] = Pm[p].Pos[CNorm] - Offset;
70  for(int p=0;p<pNPart();p++){
71  if(Pm[p].Typ != 0){
72  UnKnown[p] = Known[p];
73  continue;
74  }
75  Sol = 0.;
76  for(int d=0;d<2;d++){
77  Link[d][1] = Ln[p].Link[d*2];
78  Link[d][3] = Ln[p].Link[d*2+1];
79  Link[d][0] = Ln[Link[d][1]].Link[d*2];
80  Link[d][4] = Ln[Link[d][3]].Link[d*2+1];
81  for(int l=0,l1=-2;l<5;l++,l1++){
82  if(l==2){
83  l1 = 0;
84  continue;
85  }
86  // if(!PeriodicImage[d]){
87  // if(d == 0){
88  // if(Link[d][l] != p+l1*nEdge[1]) continue;
89  // if(Link[d][l] != p+l1) continue;
90  // }
91  // }
92  Sol += CoeffMatrix->Val(2,l)*Pm[Link[d][l]].Pos[2];
93  }
94  }
95  UnKnown[p] = .4*Known[p] - .6*Sol/CoeffMatrix->Val(2,2);
96  }
97  for(int p=0;p<pNPart();p++){
98  if(Pm[p].Typ != 0){continue;}
99  Pm[p].Pos[CNorm] = UnKnown[p] + Offset;
100  }
101  }
102  free(Known);
103  free(UnKnown);
104  delete CoeffMatrix;
105 }
110  AddRigid();
111  FillMatrix();
112  double Offset = 0.;
113  double *Known = (double *) calloc(pNPart(),sizeof(double));
114  double *UnKnown = (double *) calloc(pNPart(),sizeof(double));
115  for(int p=0;p<pNPart();p++)
116  Known[p] = Pm[p].Pos[CNorm] - Offset;
117  IntMatrix->Apply(Known,UnKnown);
118  for(int p=0;p<pNPart();p++){
119  if(Pm[p].Typ != 0){continue;}
120  Pm[p].Pos[CNorm] = UnKnown[p] + Offset;
121  }
122  free(Known);
123  free(UnKnown);
124 }
125 // #include "compcol_double.h"
126 // #include "comprow_double.h"
127 // #include "ilupre_double.h"
128 // #include "icpre_double.h"
129 // #include "iotext_double.h"
130 // void Forces::SolveLinksSparse(){
131 // int NEntries = 0;
132 // int NColMax = pNPart();
133 // int NRowMax = pNPart();
134 // for(int p=0;p<pNPart();p++){
135 // if(Pm[p].Typ != 0){
136 // IntMatrix->Set(p,p,1.);Entries++;
137 // continue;
138 // }
139 // int pym1 = Ln[p].Link[0];
140 // int pyp1 = Ln[p].Link[1];
141 // int pym2 = Ln[pym1].Link[0];
142 // int pyp2 = Ln[pyp1].Link[1];
143 // int pxm1 = Ln[p].Link[2];
144 // int pxp1 = Ln[p].Link[3];
145 // int pxm2 = Ln[pxm1].Link[2];
146 // int pxp2 = Ln[pxp1].Link[3];
147 // if(PeriodicImage[0]){ Entries += 4;
148 // }
149 // else{
150 // if(pxm2 == p-2*nEdge[1]) Entries++;
151 // if(pxm1 == p-nEdge[1]) Entries++;
152 // if(pxp1 == p+nEdge[1]) Entries++;
153 // if(pxp2 == p+2*nEdge[1]) Entries++;
154 // }
155 // Entries++;
156 // if(PeriodicImage[1]){Entries += 4;
157 // }
158 // else{
159 // if(pym2 == p-2) Entries++;
160 // if(pym1 == p-1) Entries++;
161 // if(pyp1 == p+1) Entries++;
162 // if(pyp2 == p+2) Entries++;
163 // }
164 // }
165 // double *MatVal = (double *)calloc(Entries,sizeof(double));
166 // int *ColIdx = (int *)calloc(Entries,sizeof(int));
167 // int *RowIdx = (int *)calloc(Entries,sizeof(int));
168 // for(int p=0,i=0;p<pNPart();p++){
169 // if(Pm[p].Typ != 0){
170 // ColIdx[i] = p;
171 // RowIdx[i] = p;
172 // MatVal[i] = 1.;
173 // i++;
174 // continue;
175 // }
176 // int pym1 = Ln[p].Link[0];
177 // int pyp1 = Ln[p].Link[1];
178 // int pym2 = Ln[pym1].Link[0];
179 // int pyp2 = Ln[pyp1].Link[1];
180 // int pxm1 = Ln[p].Link[2];
181 // int pxp1 = Ln[p].Link[3];
182 // int pxm2 = Ln[pxm1].Link[2];
183 // int pxp2 = Ln[pxp1].Link[3];
184 // // printf("%d)\n",p);
185 // //printf("%d %d %d %d\n",pym2,pym1,pyp1,pyp2);
186 // // printf("%d %d %d %d\n",pxm2,pxm1,pxp1,pxp2);
187 // if(PeriodicImage[0]){
188 // RowIdx[i]=p;ColIdx[i]=pxm2;ValMat[i]=CoeffMatrix->Val(2,0);i++;
189 // RowIdx[i]=p;ColIdx[i]=pxm1;ValMat[i]=CoeffMatrix->Val(2,1);i++;
190 // RowIdx[i]=p;ColIdx[i]=pxp1;ValMat[i]=CoeffMatrix->Val(2,3);i++;
191 // RowIdx[i]=p;ColIdx[i]=pxp2;ValMat[i]=CoeffMatrix->Val(2,4);i++;
192 // }
193 // else{
194 // if(pxm2 == p-2*nEdge[1])
195 // RowIdx[i]=p;ColIdx[i]=pxm2;ValMat[i]=CoeffMatrix->Val(2,0);i++;
196 // if(pxm1 == p-nEdge[1])
197 // RowIdx[i]=p;ColIdx[i]=pxm1;ValMat[i]=CoeffMatrix->Val(2,1);i++;
198 // if(pxp1 == p+nEdge[1])
199 // RowIdx[i]=p;ColIdx[i]=pxp1;ValMat[i]=CoeffMatrix->Val(2,3);i++;
200 // if(pxp2 == p+2*nEdge[1])
201 // RowIdx[i]=p;ColIdx[i]=pxp2;ValMat[i]=CoeffMatrix->Val(2,4);i++;
202 // }
203 // IntMatrix->Add(p,p,CoeffMatrix->Val(2,2));Entries++;
204 // if(PeriodicImage[1]){
205 // RowIdx[i]=p;ColIdx[i]=pym2;ValMat[i]=CoeffMatrix->Val(2,0);i++;
206 // RowIdx[i]=p;ColIdx[i]=pym1;ValMat[i]=CoeffMatrix->Val(2,1);i++;
207 // RowIdx[i]=p;ColIdx[i]=pyp1;ValMat[i]=CoeffMatrix->Val(2,3);i++;
208 // RowIdx[i]=p;ColIdx[i]=pyp2;ValMat[i]=CoeffMatrix->Val(2,4);i++;
209 // }
210 // else{
211 // if(pym2 == p-2)
212 // RowIdx[i]=p;ColIdx[i]=pym2;ValMat[i]=CoeffMatrix->Val(2,0);i++;
213 // if(pym1 == p-1)
214 // RowIdx[i]=p;ColIdx[i]=pym1;ValMat[i]=CoeffMatrix->Val(2,1);i++;
215 // if(pyp1 == p+1)
216 // RowIdx[i]=p;ColIdx[i]=pyp1;ValMat[i]=CoeffMatrix->Val(2,3);i++;
217 // if(pyp2 == p+2)
218 // RowIdx[i]=p;ColIdx[i]=pyp2;ValMat[i]=CoeffMatrix->Val(2,4);i++;
219 // }
220 // }
221 // Coord_Mat_double C(NRowMax,NColMax,Entries,MatVal,RowIdx,ColIdx);
222 // free(ColIdx);
223 // free(RowIdx);
224 // free(MatVal);
225 // }
226 // obsolete?
227 // int Forces::Update(){
228 // for(int i=0;i<IntMax;i++){
229 // double Diff = Fm[Part2Move].Dir[2] - Fm[Part2Move].Ext[2];
230 // if(QUAD(Diff) < DiffForce ){
231 // printf("%d,%d) Force %lf= %lf - %lf Pos %lf \n",Part2Move,i,Diff,Fm[Part2Move].Dir[2],Fm[Part2Move].Ext[2],Pm[Part2Move].Pos[2]);
232 // break;
233 // }
234 // double Sign=1.;
235 // if( Diff < 0. ){
236 // Sign = -1.;
237 // Pm[Part2Move].Pos[2] -= IncrDist;
238 // }
239 // else
240 // Pm[Part2Move].Pos[2] += IncrDist;
241 // for(int p=0,pp=Part2Move-1,ppp=Part2Move+1;p<pNPart();p++)
242 // ;
243 // if( Diff < 0. ){
244 // Pm[Part2Move].Pos[2] -= IncrDist;
245 // }
246 // else if(Diff > 0.){
247 // Pm[Part2Move].Pos[2] += IncrDist;
248 // }
249 // ForceFieldLine();
250 // Diff = Fm[Part2Move].Dir[2] - Fm[Part2Move].Ext[2];
251 // if(QUAD(Diff) < DiffForce ){
252 // printf("%d,%d) Force %lf= %lf - %lf Pos %lf \n",Part2Move,i,Diff,Fm[Part2Move].Dir[2],Fm[Part2Move].Ext[2],Pm[Part2Move].Pos[2]);
253 // break;
254 // }
255 // if(Part2Move>0){
256 // if( Diff < 0. ){
257 // Pm[Part2Move-1].Pos[2] += IncrDist;
258 // }
259 // else if(Diff > 0.){
260 // Pm[Part2Move-1].Pos[2] -= IncrDist;
261 // }
262 // }
263 // if(Part2Move<pNPart()-1){
264 // if( Diff < 0. ){
265 // Pm[Part2Move+1].Pos[2] += IncrDist;
266 // }
267 // else if(Diff > 0.){
268 // Pm[Part2Move+1].Pos[2] -= IncrDist;
269 // }
270 // }
271 // ForceFieldLine();
272 // }
273 // }
274 // // obsolete?
275 // int Forces::MinHelfrich(){
276 // SPLINE Par;
277 // SPLINE Sp2;
278 // double a1=0.;
279 // double a2=0.;
280 // double a3=0.;
281 // double a4=0.;
282 // double Ref=0.;
283 // double Dx = pEdge(0)/(double)(NEdge-1);
284 // for(int p=0;p<pNPart()-1;p++){
285 // if(p < 1){
286 // Sp2 = Mat->Forth(Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,Pm[p+4].Pos,0,2);
287 // Ref = Pm[p].Pos[0] - Pm[p+2].Pos[0];
288 // }
289 // else if(p < 2){
290 // Sp2 = Mat->Forth(Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,Pm[p+3].Pos,0,2);
291 // Ref = Pm[p].Pos[0] - Pm[p+1].Pos[0];
292 // }
293 // else if(p < pNPart()-2){
294 // Sp2 = Mat->Forth(Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,Pm[p+2].Pos,0,2);
295 // Ref = 0.;
296 // }
297 // else if(p == pNPart()-2){
298 // Sp2 = Mat->Forth(Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,Pm[p+1].Pos,0,2);
299 // Ref = Pm[p].Pos[0] - Pm[p-1].Pos[0];
300 // }
301 // else if(p <= pNPart()-1){
302 // Sp2 = Mat->Forth(Pm[p-4].Pos,Pm[p-3].Pos,Pm[p-2].Pos,Pm[p-1].Pos,Pm[p].Pos,0,2);
303 // Ref = Pm[p].Pos[0] - Pm[p-2].Pos[0];
304 // }
305 // if(Pm[p].Typ == 1 || Pm[p].Typ == 2)
306 // continue;
307 // a2 = 2.*Sp2.a2 + 6.*Sp2.a3*Ref + 12.*Sp2.a4*Ref*Ref;
308 // a4 = 24.*Sp2.a4;
309 // a3 = a4*(12.*Dx*Dx - 24.*Kf.Lap/Kf.SLap) - 2.*a2;
310 // //a3 = a4*(12.*Dx*Dx - 24.*.1) - 2.*a2;
311 // a3 /= Dx;
312 // a1 = Pm[p+1].Pos[2] - Pm[p].Pos[2] - a2*Dx*Dx - a3*Dx*Dx*Dx - a4*Dx*Dx*Dx*Dx;
313 // a1 /= Dx;
314 // double x = Pm[p+1].Pos[0] - Pm[p].Pos[0];
315 // printf("%d) %lf %lf %lf %lf -> %lf\n",p,a1,a2,a3,a4,Pm[p].Pos[2] + a1*x + a2*x*x + a3*x*x*x + a4*x*x*x*x);
316 // Pm[p+1].Pos[2] = Pm[p].Pos[2] + a1*x + a2*x*x + a3*x*x*x + a4*x*x*x*x;
317 // }
318 // }
319 //---------------------------------RIGID--------------------------
322  for(int n=0;n<pNNano();n++){
323  Vettore Axis(Nano[n].Axis,3);
324  Vettore AMom(Nano[n].AMom,3);
325  Vettore BackBone(3);
326  Vettore VelTang(3);
327  for(int d=0;d<3;d++){
328  BackBone.x[d] = Nano[n].Axis[d]*Nano->Height*.5;
329  double Ran = Nano[n].Zeta * (2.*Mat->Casuale() - 1.);
330  double Dis = - Nano[n].Gamma*Nano[n].AVel[d];
331  Nano[n].AMom[d] += Ran + Dis;
332  }
333  Matrice InTensor(3,3);
334  Matrice RotT(3,3);
335  Matrice Resp(3,3);
336  Matrice Resp2(3,3);
337  Vettore AVel(3);
338  Vettore Normal(0.,0.,1.);
339  Vettore Ax(3);
340  Axis.VetV(&Axis,&Normal);
341  double Angle = Ax.Angle(&Axis,&Normal);
342  Quadri q(Ax.x,Angle);
343  Matrice Rot(q,3);
344  Rot.CopyOn(&RotT);
345  RotT.Transpose();
346  InTensor.Set(0,0,1./(Nano->Mass*(.25*SQR(Nano->Rad)+1./12.*SQR(Nano->Height) ) ) );
347  InTensor.Set(1,1,1./(Nano->Mass*(.25*SQR(Nano->Rad)+1./12.*SQR(Nano->Height) ) ) );
348  InTensor.Set(2,2,1./(.5*Nano->Mass*SQR(Nano->Rad)) );
349  // Calculating r^t I^-1 r L = w
350  Resp.Mult(RotT,InTensor);
351  Resp2.Mult(Resp,Rot);
352  // MatrVect(Resp2,Nano->AMom,AVel);
353 
354  //printf("%d) %d\n",n,Gen->Step);
355  for(int d=0;d<3;d++){
356  double tmp = .5*Nano[n].Force[d]*pDeltat()/Nano[n].Mass;
357  Nano[n].Vel[d] += tmp;
358  Nano[n].AVel[d] += Nano[n].AMom[d]/Nano[n].Mass*.01;
359  Nano[n].Pos[d] += Nano[n].Vel[d]*pDeltat();
360  Nano[n].Pos[d] -= floor(Nano[n].Pos[d]/pEdge(d))*pEdge(d);
361  }
362  if(Nano[n].Shape != 2) continue;
363  Vettore Omega(Nano[n].AVel,3);
364  for(int d=0;d<3;d++)
365  BackBone.x[d] = Nano[n].Height*.5*Nano[n].Axis[d];
366  VelTang.VetV(&Omega,&BackBone);
367  for(int d=0;d<3;d++){
368  BackBone.x[d] += VelTang.x[d]*pDeltat();
369  Axis.x[d] = BackBone.x[d];
370  }
371  Axis.Normalize();
372  Axis.Export(Nano[n].Axis);
373  //printf("%lf %lf %lf\n",Nano[n].Vel[0],Nano[n].Vel[1],Nano[n].Vel[2]);
374  //printf("%lf %lf %lf\n",Nano[n].Force[0],Nano[n].Force[1],Nano[n].Force[2]);
375  //printf("%lf %lf %lf\n",Nano[n].Pos[0],Nano[n].Pos[1],Nano[n].Pos[2]);
376  //printf("%lf %lf %lf\n",Nano[n].AMom[0],Nano[n].AMom[1],Nano[n].AMom[2]);
377  }
378 }
383  for(int n=0;n<pNNano();n++) {
384  for(int d=0;d<3;d++){
385  Nano[n].AVel[d] += .5*Nano[n].AMom[d];
386  Nano[n].Vel[d] += .5*Nano[n].Force[d]*pDeltat()/Nano[n].Mass;
387  }
388  }
389 }
390 //------------------------------------------------------------
391 //-------------------------MONTE-CARLO------------------------
392 //------------------------------------------------------------
394 int Forces::IfMetropolis(double Arg,double Weight){
395  if(exp(pBeta()*Arg)*Weight > 1. ) return 1;
396  double Ran = Mat->Casuale();
397  if(exp(pBeta()*Arg)*Weight > Ran ){
398  return 1;
399  }
400  return 0;
401 }
402 //-----------------Single-particle----------------------------
405  Pm[p].Pos[0] = pEdge(0)*Mat->Casuale();
406  Pm[p].Pos[1] = pEdge(1)*Mat->Casuale();
407  Pm[p].Pos[2] = pEdge(2)*Mat->Casuale();
408  Pc->AddPart(p,Pm[p].Pos);
409 }
414  int p = pNPart();
415  ReSetNPart(pNPart()+1);
416  InsertBead(p);
417  double Pot[3];
418  double NrgDiff = CalcNrgBead(p,Pot);
419  double Arg = ChemPotId + ChemPotEx - NrgDiff;
420  double Weight = pVol()/(double)(pNPart()+1);
421  if( !IfMetropolis(Arg,Weight) ){
422  Pc->RemPart(p,Pm[p].Pos);
423  ReSetNPart(pNPart()-1);
424  return 0;
425  }
426  //printf("Added a particle %d\n",p);
427  return 1;
428 }
433  printf("remove\n");
434  int p = (int)(Mat->Casuale()*pNPart());
435  double Pot[3];
436  double NrgDiff = CalcNrgBead(p,Pot);
437  double Arg = - ChemPotId - ChemPotEx - NrgDiff;
438  if( !IfMetropolis(Arg,pNPart()/pVol()) ){
439  return 0;
440  }
441  SwapPart(p,pNPart()-1);
442  Pc->SwapPart(p,Pm[p].Pos,pNPart()-1,Pm[pNPart()-1].Pos);
443  Pc->RemPart(pNPart()-1,Pm[pNPart()-1].Pos);
444  ReSetNPart(pNPart()-1);
445  return 1;
446 }
448 int Forces::MoveBead(int p){
449  for(int d=0;d<3;d++){
450  Pm[p].Pos[d] += Mat->Gaussiano(0.,GaussVar);
451  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
452  }
453 }
458  int p = (int)(Mat->Casuale()*pNPart());
459  if(Pm[p].Typ >= 2) return 1;
460  //old situa
461  double Pos[3] = {Pm[p].Pos[0],Pm[p].Pos[1],Pm[p].Pos[2]};
462  double Pot[3] = {0.,0.,0.};
463  double Nrg0 = CalcNrgBead(p,Pot);
464  //RemDens(p,p+1);
465  //double Nrg0 = OldNrgPm[p*3] + OldNrgPm[p*3+1] + OldNrgPm[p*3+2];
466  //new situa
467  MoveBead(p);
468  Pc->MovePart(p,Pos,Pm[p].Pos);
469  //AddDens(p,p+1);
470  double Nrg1 = CalcNrgBead(p,Pot);
471  //like the new situa?
472  double NrgDiff = Nrg1 - Nrg0;
473  if(!IfMetropolis(-NrgDiff,1.)){
474  //RemDens(p,p+1);
475  Pc->MovePart(p,Pm[p].Pos,Pos);
476  for(int d=0;d<3;d++) Pm[p].Pos[d] = Pos[d];
477  //AddDens(p,p+1);
478  return 0;
479  }
480  //OldNrgSys = SumDens(0,pNPart());
481  OldNrgSys -= Pot[2];//DensFuncNrgCh(c);
482  return 1;
483 }
484 //-----------------------------Chains---------------------
489  //getting the distribution
490  double Norm = 0.;
491  double *LineS = new double[NBin];
492  for(int i=0;i<NBin;i++){
493  FirstBeadDistr[i] = 0.;
494  }
495  for(int c=0;c<pNChain();c++){
496  int p1 = c*pNPCh();
497  int i = (int)(Pm[p1].Pos[CNorm]*NBin*pInvEdge(CNorm));
498  if(i < 0 || i >= NBin) continue;
499  FirstBeadDistr[i] += 1.;
500  Norm += 1.;
501  }
502  //normalizing and smoothing
503  for(int v=0;v<NBin;v++) LineS[v] = FirstBeadDistr[v];
504  InterBSpline1D(LineS,FirstBeadDistr,NBin,NBin);
505  for(int v=0;v<NBin;v++) LineS[v] = FirstBeadDistr[v];
506  InterBSpline1D(LineS,FirstBeadDistr,NBin,NBin);
507  Norm = Norm >= 0. ? 1./Norm : 1.;
508  for(int i=0;i<NBin;i++){
509  FirstBeadDistr[i] *= Norm;
510  }
511  //cumulative
512  for(int i=1;i<NBin;i++){
513  FirstBeadDistr[i] += FirstBeadDistr[i-1];
514  if(FirstBeadDistr[i] > .25 && FirstBeadDistr[i-1] < .25)
515  BorderBias[0] = i/NBin*pEdge(CNorm);
516  else if(FirstBeadDistr[i] > .75 && FirstBeadDistr[i-1] < .75)
517  BorderBias[1] = i/NBin*pEdge(CNorm);
518  }
519  delete [] LineS;
520 }
524 void Forces::ConsiderCh(int c){
525  int p1 = c*pNPCh();
526  for(int p=0;p<pNPCh();p++) Pc->AddPart(p+p1,Pm[p+p1].Pos);
527  AddDens(p1,p1+pNPCh());
528 }
532 void Forces::IgnoreCh(int c){
533  int p1 = c*pNPCh();
534  RemDens(p1,p1+pNPCh());
535  for(int p=0;p<pNPCh();p++) Pc->RemPart(p+p1,Pm[p+p1].Pos);
536 }
541  int c2 = pNChain()-1;
542  int p1 = c*pNPCh();
543  int p2 = c2*pNPCh();
544  RemDens(p1,p1+pNPCh());
545  if(c != c2){
546  for(int e=0;e<3;e++)
547  OldNrgCh[c*3+e] = OldNrgCh[c2*3+e];
548  RemDens(p2,p2+pNPCh());
549  // for(int p=0;p<pNPCh();p++)
550  // Pc->RemPart(p2+p,Pm[p2+p].Pos);
551  // for(int p=0;p<pNPCh();p++)
552  // Pc->RemPart(p1+p,Pm[p1+p].Pos);
553  // SwapChain(c,pNChain()-1);
554  // for(int p=0;p<pNPCh();p++)
555  // Pc->AddPart(p1+p,Pm[p1+p].Pos);
556  for(int p=0;p<pNPCh();p++)
557  Pc->SwapPart(p+p1,Pm[p+p1].Pos,p2+p,Pm[p2+p].Pos);
558  SwapChain(c,pNChain()-1);
559  AddDens(p1,p1+pNPCh());
560  }
561  for(int p=0;p<pNPCh();p++)
562  Pc->RemPart(p2+p,Pm[p2+p].Pos);
563  ReSetNPart(pNPart()-pNPCh());
564  ReSetNChain(pNChain()-1);
565 }
569 void Forces::SaveCh(int c){
570  int p1 = c*pNPCh();
571  for(int p=0;p<pNPCh();p++){
572  for(int d=0;d<3;d++){
573  OldPos[p][d] = Pm[p+p1].Pos[d];
574  }
575  }
576 }
580 void Forces::ReInsertCh(int c){
581  IgnoreCh(c);
582  int p1 = c*pNPCh();
583  //old situa
584  for(int p=0;p<pNPCh();p++){
585  for(int d=0;d<3;d++){
586  Pm[p+p1].Pos[d] = OldPos[p][d];
587  }
588  }
589  ConsiderCh(c);
590 }
592 double Forces::InsertCh(int c){
593  int p1 = c*pNPCh();
594  double Weight = 1.;
595  ReSetNPart(pNPart()+pNPCh());
596  ReSetNChain(pNChain()+1);
597  Pm[p1].Pos[CLat1] = Mat->Casuale()*pEdge(CLat1);
598  Pm[p1].Pos[CLat2] = Mat->Casuale()*pEdge(CLat2);
599  Pm[p1].Pos[CNorm] = Mat->Casuale()*pEdge(CNorm);
600  if(VAR_IF_TYPE(CalcMode,CALC_BIL_BIAS)){
601  Pm[p1].Pos[CNorm] = Mat->Casuale()*(BorderBias[1] - BorderBias[0]) + BorderBias[0];
602  //Mat->RandDiscrProb(FirstBeadDistr,NBin)*pEdge(CNorm);
603  Weight *= (BorderBias[1]-BorderBias[0])*pEdge(CLat1)*pEdge(CLat2)/pVol();
604  }
605  else if(VAR_IF_TYPE(CalcMode,CALC_SPH_BIAS)){
606  double Incr = 0.2;
607  double z = 2.0 * Mat->Casuale() - 1.0;
608  double t = 2.0 * M_PI * Mat->Casuale();
609  double w = sqrt( 1 - z*z );
610  double x = w * cos( t );
611  double y = w * sin( t );
612  double Rad = Mat->Casuale()*Incr + Nano->Rad;
613  Pm[p1].Pos[CLat1] = x*Rad + Nano->Pos[0];
614  Pm[p1].Pos[CLat2] = y*Rad + Nano->Pos[1];
615  Pm[p1].Pos[CNorm] = z*Rad + Nano->Pos[2];
616  double Vol = 4./3.*M_PI*( CUBE(Nano->Rad + Incr) - CUBE(Nano->Rad));
617  Weight *= Vol/pVol();
618  }
619  Weight *= InsertRest(p1,0);
620  // for(int p=0;p<pNPCh();p++){
621  // fprintf(StatFile1,"%lf %lf %lf 0. 0. 0. %d\n",Pm[p1+p].Pos[0],Pm[p1+p].Pos[1],Pm[p1+p].Pos[2],pType(p+p1));
622  // }
623  // fflush(StatFile1);
624  ConsiderCh(c);
625  return Weight;
626 }
628 double Forces::InsertRest(int pCurr,int StartPos){
629  double Weight = 1.;
630  for(int p=StartPos+1;p<pNPCh();p++){
631  for(int d=0;d<3;d++){
632  double Ran = Mat->Gaussiano(0.,GaussVar);
633  Pm[p+pCurr].Pos[d] = Pm[p-1+pCurr].Pos[d] + Ran;
634  Pm[p+pCurr].Pos[d] -= floor(Pm[p+pCurr].Pos[d]*pInvEdge(d))*pEdge(d);
635  }
636  }
637  for(int p=StartPos-1;p>=0;p--){
638  for(int d=0;d<3;d++){
639  double Ran = Mat->Gaussiano(0.,GaussVar);
640  Pm[p+pCurr].Pos[d] = Pm[p+1+pCurr].Pos[d] + Ran;
641  Pm[p+pCurr].Pos[d] -= floor(Pm[p+pCurr].Pos[d]*pInvEdge(d))*pEdge(d);
642  }
643  }
644  return Weight;
645 }
650  int c = (int)(Mat->Casuale()*pNChain());
651  int p1 = c*pNPCh();
652  double Pot[3];
653  // old situa
654  SaveCh(c);
655  IgnoreCh(c);
656  double NrgOld = CalcNrgCh(c,Pot);
657  // new situa
658  ReSetNPart(pNPart()-pNPCh());
659  ReSetNChain(pNChain()-1);
660  InsertCh(c);
661  double NrgNew = CalcNrgCh(c,Pot);
662  //like the new situa?
663  double NrgDiff = NrgNew - NrgOld;
664  if(!IfMetropolis(-NrgDiff,1.)){
665  ReInsertCh(c);
666  return 0;
667  }
668  OldNrgCh[c*3+2] = NrgNew;
669  OldNrgSys += NrgDiff;
670  return 1;
671 }
676  if(pNPart() <= 0) return 0;
677  int c = (int)(Mat->Casuale()*pNChain());
678  int p1 = c*pNPCh();
679  double Pot[3];
680  //like the new situa?
681  double NrgDiff = CalcNrgCh(c,Pot);
682  //double NrgDiff = OldNrgCh[c*3+2];
683  //fprintf(StatFile1,"%lf\n",NrgDiff);
684  //fflush(StatFile1);
685  double Arg = -ChemPotId -ChemPotEx +NrgDiff;
686  double Weight = pNChain()/pVol();
687  if(!IfMetropolis(Arg,Weight) ){
688  return 0;
689  }
690  RemChFromSys(c);
691  OldNrgSys -= NrgDiff;
692  return 1;
693 }
698  //old situa
699  int c = pNChain();
700  int p1 = c*pNPCh();
701  double Pot[3];
702  //new situa
703  InsertCh(c);
704  //like the new situa?
705  double NrgDiff = CalcNrgCh(c,Pot);
706  //fprintf(StatFile2,"%lf\n",NrgDiff);
707  fflush(StatFile2);
708  //printf("add %lf sys %lf\n",NrgDiff,OldNrgSys);
709  double Arg = ChemPotId + ChemPotEx - NrgDiff;
710  double Weight = pVol()/(double)(pNChain()+1);
711  //printf("Ins %lf\n",NrgDiff);
712  if(!IfMetropolis(Arg,Weight)){
713  RemChFromSys(c);
714  return 0;
715  }
716  //printf("accepted\n");
717  //printf("Inserted %lf\n",NrgDiff);
718  OldNrgCh[c*3+2] = NrgDiff;
719  OldNrgSys += NrgDiff;
720  return 1;
721 }
722 //---------------------------Bias-------------------------------
724 double Forces::RemoveChBias(int c){
725  int p1 = c*pNPCh();
726  RemDens(p1,p1+pNPCh());
727  for(int p=0;p<pNPCh();p++) Pc->RemPart(p+p1,Pm[p+p1].Pos);
728  double NrgDiff = DensFuncNrgGhost(Pm[p1].Pos,p1,Pm[p1].Typ);
729  //fprintf(StatFile1,"%lf %lf\n",NrgDiff,NanoNrg(p1));
730  NrgDiff += CalcBendingGhost(Pm[p1].Pos,p1);
731  NrgDiff += NanoNrg(p1);
732  double Weight = exp(-NrgDiff+NrgPBead);
733  Pc->AddPart(p1,Pm[p1].Pos);
734  AddDens(p1,p1+1);
735  for(int p=p1+1;p<p1+pNPCh();p++){
736  Weight *= WeightSetBond(p,Pm[p].Typ);
737  Pc->AddPart(p,Pm[p].Pos);
738  AddDens(p,p+1);
739  }
740  return Weight;
741 }
743 double Forces::WeightSetBond(int p,int Type){
744  double NrgDiff = DensFuncNrgGhost(Pm[p].Pos,p,Pm[p].Typ);
745  NrgDiff += CalcBendingGhost(Pm[p].Pos,p);
746  NrgDiff += NanoNrg(p);
747  CumProbBias[0] = exp(-NrgDiff+NrgPBead);
748  for(int t=1;t<NTrialBias;t++){
749  for(int d=0;d<3;d++){
750  BondPosBias[t][d] = Pm[p-1].Pos[d] + Mat->Gaussiano(0.,GaussVar);
751  BondPosBias[t][d] -= floor(BondPosBias[t][d]*pInvEdge(d))*pEdge(d);
752  }
753  NrgDiff = DensFuncNrgGhost(BondPosBias[t],p,Type);
754  //fprintf(StatFile1,"%lf %lf\n",NrgDiff,NanoNrg(BondPosBias[t],Type));
755  NrgDiff += CalcBendingGhost(BondPosBias[t],p);
756  NrgDiff += NanoNrg(BondPosBias[t],Type);
757  CumProbBias[t] = exp(-NrgDiff+NrgPBead);
758  }
759  double Weight = 0.;
760  for(int t=0;t<NTrialBias;t++){
761  Weight += CumProbBias[t];
762  }
763  return Weight/(double)NTrialBias;
764 }
766 double Forces::InsertChBias(int c){
767  int p1 = c*pNPCh();
768  double Weight = 1.;
769  ReSetNPart(pNPart()+pNPCh());
770  ReSetNChain(pNChain()+1);
771  Pm[p1].Pos[CLat1] = Mat->Casuale()*pEdge(CLat1);
772  Pm[p1].Pos[CLat2] = Mat->Casuale()*pEdge(CLat2);
773  Pm[p1].Pos[CNorm] = Mat->Casuale()*pEdge(CNorm);
774  if(VAR_IF_TYPE(CalcMode,CALC_BIL_BIAS)){
775  Pm[p1].Pos[CNorm] = Mat->Casuale()*(BorderBias[1] - BorderBias[0]) + BorderBias[0];
776  //Mat->RandDiscrProb(FirstBeadDistr,NBin)*pEdge(CNorm);
777  Weight *= (BorderBias[1]-BorderBias[0])*pEdge(CLat1)*pEdge(CLat2)/pVol();
778  }
779  else if(VAR_IF_TYPE(CalcMode,CALC_SPH_BIAS)){
780  double Incr = 0.2;
781  double z = 2.0 * Mat->Casuale() - 1.0;
782  double t = 2.0 * M_PI * Mat->Casuale();
783  double w = sqrt( 1 - z*z );
784  double x = w * cos( t );
785  double y = w * sin( t );
786  double Rad = Mat->Casuale()*Incr + Nano->Rad;
787  Pm[p1].Pos[CLat1] = x*Rad + Nano->Pos[0];
788  Pm[p1].Pos[CLat2] = y*Rad + Nano->Pos[1];
789  Pm[p1].Pos[CNorm] = z*Rad + Nano->Pos[2];
790  double Vol = 4./3.*M_PI*( CUBE(Nano->Rad + Incr) - CUBE(Nano->Rad));
791  Weight *= Vol/pVol();
792  //fprintf(StatFile2,"%lf %lf %lf %d\n",Pm[p1].Pos[0],Pm[p1].Pos[1],Pm[p1].Pos[2],pType(p1));
793  //fflush(StatFile2);
794  }
795  double NrgDiff = DensFuncNrgGhost(Pm[p1].Pos,p1,Pm[p1].Typ);
796  fprintf(StatFile2,"%lf %lf\n",NrgDiff,NanoNrg(p1));
797  NrgDiff += CalcBendingGhost(Pm[p1].Pos,p1);
798  NrgDiff += NanoNrg(p1);
799  Weight *= exp(-NrgDiff+NrgPBead);
800  // printf("Weight1 %lf %lf\n",NrgDiff,Weight);
801  Pc->AddPart(p1,Pm[p1].Pos);
802  AddDens(p1,p1+1);
803  for(int p=p1+1;p<(c+1)*pNPCh();p++){
804  double Weight1 = CreateSetBond(p,Pm[p].Typ);
805  Weight *= Weight1;
806  Pc->AddPart(p,Pm[p].Pos);
807  AddDens(p,p+1);
808  }
809  // for(int p=0;p<pNPCh();p++){
810  // fprintf(StatFile1,"%lf %lf %lf %d\n",Pm[p1+p].Pos[0],Pm[p1+p].Pos[1],Pm[p1+p].Pos[2],pType(p+p1));
811  // }
812  // fflush(StatFile1);
813  return Weight;
814 }
816 double Forces::CreateSetBond(int p,int Type){
817  int tChoosen = 0;
818  for(int t=0;t<NTrialBias;t++){
819  for(int d=0;d<3;d++){
820  BondPosBias[t][d] = Pm[p-1].Pos[d] + Mat->Gaussiano(0.,GaussVar);
821  BondPosBias[t][d] -= floor(BondPosBias[t][d]*pInvEdge(d))*pEdge(d);
822  }
823  double NrgDiff = DensFuncNrgGhost(BondPosBias[t],p,Type);
824  fprintf(StatFile2,"%lf %lf\n",NrgDiff,NanoNrg(BondPosBias[t],Type));
825  NrgDiff += CalcBendingGhost(BondPosBias[t],p);
826  NrgDiff += NanoNrg(BondPosBias[t],Type);
827  CumProbBias[t] = exp(-NrgDiff+NrgPBead);
828  }
829  double Weight = 0.;
830  for(int t=0;t<NTrialBias;t++){
831  Weight += CumProbBias[t];
832  }
833  CumProbBias[0] /= Weight;
834  double Ran = Mat->Casuale();
835  for(int t=1;t<NTrialBias;t++){
836  CumProbBias[t] /= Weight;
837  CumProbBias[t] += CumProbBias[t-1];
838  }
839  for(int t=0;t<NTrialBias;t++){
840  if(Ran < CumProbBias[t]){
841  tChoosen = t;
842  break;
843  }
844  }
845  for(int d=0;d<3;d++){
846  Pm[p].Pos[d] = BondPosBias[tChoosen][d];
847  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
848  }
849  return Weight/(double)NTrialBias;
850 }
855  if(pNPart() <= 0) return 0;
856  int c = (int)(Mat->Casuale()*pNChain());
857  int p1 = c*pNPCh();
858  double Weight1 = RemoveChBias(c);
859  double Arg = - ChemPotId - ChemPotEx ;
860  double Weight = pNChain()/(Weight1*pVol());
861  //printf("rem %lf\n",Weight1);
862  //fprintf(TempFile,"1 %lf\n",NrgDiff);
863  if(!IfMetropolis(Arg,Weight)){
864  return 0;
865  }
866  //fprintf(TempFile,"2 %lf\n",NrgDiff);
867  RemChFromSys(c);
868  return 1;
869 }
874  //add one chain
875  int c = pNChain();
876  int p1 = c*pNPCh();
877  double Weight1 = InsertChBias(c);
878  double Arg = ChemPotId + ChemPotEx;
879  double Weight = Weight1*pVol()/(double)(pNChain()+1);
880  if(!IfMetropolis(Arg,Weight)){
881  RemChFromSys(c);
882  return 0;
883  }
884  return 1;
885 }
886 //-------------------WIDOM--------------------------
890 void Forces::WidomInsert(double *NrgDiff){
891  int p1 = pNPart();
892  InsertBead(p1);
893  *NrgDiff = DensFuncNrgGhost(Pm[p1].Pos,p1,Pm[p1].Typ);
894  *NrgDiff += CalcBendingGhost(Pm[p1].Pos,p1);
895  *NrgDiff += NanoNrg(p1);
896  Pc->RemPart(p1,Pm[p1].Pos);
897 }
901 void Forces::WidomRemove(double *NrgDiff,int p){
902  double Pot[3];
903  *NrgDiff = CalcNrgBead(p,Pot);
904 }
908 void Forces::WidomInsertCh(double *NrgDiff){
909  int c = pNChain();
910  int p1 = c*pNPCh();
911  InsertCh(c);
912  CalcNrgCh(c,NrgDiff);
913  RemChFromSys(c);
914 }
918 void Forces::WidomRemoveCh(double *NrgDiff,int c){
919  NrgDiff[0] = OldNrgCh[c*3 ];
920  NrgDiff[1] = OldNrgCh[c*3+1];
921  NrgDiff[2] = OldNrgCh[c*3+2];
922 }
926 void Forces::WidomBiasChIn(double *Weight){
927  int c = pNChain();
928  int p1 = c*pNPCh();
929  *Weight = InsertChBias(c);
930  RemChFromSys(c);
931 }
935 void Forces::WidomBiasChOut(double *Weight,int c){
936  *Weight = RemoveChBias(c);
937 }
938 //----------------------MOLECULAR-DYNAMICS----------------------
943  for(int p=0;p<pNPart();p++){
944  if(Pm[p].Typ != 0) continue;
945  for(int d=0;d<3;d++){
946  //if(fabs(Fm[p].Dir[d]) > 500.) continue;
947  Pm[p].Vel[d] += .5*(Fm[p].Dir[d] + Fm[p].Ext[d])*pDeltat();
948  Pm[p].Pos[d] += Pm[p].Vel[d]*pDeltat();
949  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
950  //printf("%d %d) %lf %g %lf\n",p,d,Pm[p].Pos[d],Fm[p].Dir[d],Pm[p].Vel[d]);
951  }
952  }
953 }
958  double Gamma = Viscosity;//3.*3.14*10.*Viscosity;
959  double Zeta = sqrt(12.*2.*1.*Gamma/pDeltat());
960  for(int p=0;p<pNPart();p++){
961  for(int d=0;d<3;d++){
962  Fm[p].Dir[d] += Zeta*(2.*Mat->Casuale() - 1.);
963  //Fm[p].Dir[d] += Zeta*Mat->Gaussiano(0.,1.);
964  Fm[p].Dir[d] -= Gamma*Pm[p].Vel[d];
965  }
966  }
967 }
972  double Sigma = sqrt(pTemp());
973  for(int p=0;p<pNPart();p++){
974  if(Mat->Casuale() < pDeltat()){
975  for(int d=0;d<3;d++){
976  Pm[p].Vel[d] = Mat->Gaussiano(0.,Sigma);
977  }
978  }
979  }
980 }
985  double Temp = 0.;
986  for(int p=0;p<pNPart();p++){
987  for(int d=0;d<3;d++){
988  Pm[p].Vel[d] += .5*Fm[p].Dir[d]*pDeltat();
989  Temp += SQR(Pm[p].Vel[d]);
990  }
991  }
992  Temp = Temp/(3*pNPart());
993  double Norm = 1./sqrt(Temp);
994  for(int p=0;p<pNPart();p++){
995  for(int d=0;d<3;d++){
996  Pm[p].Vel[d] *= Norm;
997  }
998  }
999 }
1004  printf("Thermostat: ");
1005  if(Mode==THERM_LANG){
1007  printf("Langevin\n");
1008  }
1009  else if((Mode==THERM_AND)){
1011  printf("Andersen\n");
1012  }
1013  else if((Mode==THERM_BERE)){
1015  printf("Berendsen\n");
1016  }
1017  else if((Mode==THERM_NO)){
1019  printf("no thermostat\n");
1020  }
1021  else{
1022  printf("mode not present\n");
1023  exit(1);
1024  }
1025 }
1030  for(int p=0;p<pNPart();p++){
1031  for(int d=0;d<3;d++){
1032  Pm[p].Vel[d] += .5*(Fm[p].Dir[d] + Fm[p].Ext[d])*pDeltat();
1033  Fm[p].Dir[d] = 0.;
1034  }
1035  }
1036 }
int TryMoveCh()
Trial desplacement of a chain.
void LangevinTherm()
Langevin thermostat.
void Export(double *x)
Export.
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
void WidomBiasChIn(double *Weight)
Widom with Rosenbluth weight.
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double InsertChBias(int c)
Put a chain in the box.
int InterBSpline1D(double *PlIn, double *PmOut, int NIn, int NOut)
1-d BSpline
double InsertRest(int pCurr, int StartPos)
Build the rest of the chain.
void VelVerlet1()
First step of the velocity Verlet.
double RemoveChBias(int c)
Remove a chain in the box.
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
int TryRemoveChBias()
Trial biased removal of a chain.
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
void AddRigid()
Add all rigid bodies as a boundary condition.
CALC_THERM CalcTherm
Pointer to the energy calculation function.
Definition: Forces.h:439
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double CalcBendingGhost(double *Pos, int pExt)
Calculate the bending energy for a ghost particle.
double Normalize()
Normlizes a Vettore.
double Gamma
Friction term.
Definition: VarData.h:451
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
double Lap
Prefactor of the Laplacian.
Definition: Forces.h:147
double CreateSetBond(int p, int t)
Create a set of bonds for the configurational bias.
void RemPart(const int p, double *Pos)
Remove a particle form the cell c.
Definition: Cubo.cpp:307
double Mass
Mass.
Definition: VarData.h:443
double Rad
Size.
Definition: VarData.h:445
int MoveBead(int p)
Move a particle in the box.
double Viscosity
Viscosity of the medium.
Definition: Forces.h:731
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
double GaussVar
Standard deviation of the gaussian chain.
Definition: Forces.h:725
double Height
Height of the cylinder.
Definition: VarData.h:449
double Vel[3]
Velocity.
Definition: VarData.h:431
int CLat2
lateral coordinate
Definition: VarData.h:1078
double El[3]
Elastic force.
Definition: Forces.h:151
void CopyOn(Matrice *B)
Copy on a matrix.
int TryInsertCh()
Trial insertion of a chain.
void Mult(Matrice &A, Matrice &B)
Multiplication between two matrices.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
void FillMatrix()
Fill the entries of the interaction matrix.
Definition: Forces.cpp:270
int IfMetropolis(double ExpArg, double Weight)
Metropolis acceptance criterium.
void WidomInsert(double *NrgDiff)
Widom insertion.
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 * x
Where the data are stored.
void SolveLinksIterative()
Solve a system of four oder differential equation of particles connected by links.
void NoTherm()
No thermostat.
Definition: Forces.h:430
int TryInsertChBias()
Trial biased insertion of a chain.
void WidomInsertCh(double *NrgDiff)
Widom insertion.
void WidomRemoveCh(double *NrgDiff, int c)
Widom removal.
double pVol()
xyzr edges of the simulation box
Definition: VarData.h:922
double ChemPotEx
Chemical potential of the particles.
Definition: Forces.h:723
void SwapPart(int p1, int p2)
Swap two particle.
Definition: VarDataEl.cpp:378
int CalcMode
Calculation mode.
Definition: Forces.h:747
double Pos[3]
Position.
Definition: VarData.h:427
void ChooseThermostat(int Mode)
Choose a calculation mode.
void ConsiderCh(int c)
Consider a chain in the system (densities, pairlist)
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
void IgnoreCh(int c)
Ignore a chain in the system (densities, pairlist)
void VelVerletRigid()
Velocity Verlet for a rigid body, first step.
int pNPCh()
Number of particle per chain.
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
int InsertBead(int p)
Insert a particle in the box.
void Transpose()
Transpose the matrix.
int ReSetNPart(int NewNPart)
Realloc the number of particles.
Definition: Forces.cpp:820
double NanoNrg(int p)
Exchange energy with the nano.
int TryMove()
Trial movement.
void BackBone(double *Line, int NBin)
Describe the backbone of a filament.
void AndersenTherm()
Andersen thermostat.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double Casuale()
Random uniform number.
double Val(int row)
Returns a value in 1d.
KFORCES Kf
Prefactor of the forces.
Definition: Forces.h:779
double Zeta
Stochastic term.
Definition: VarData.h:453
void StudySys()
Obtain informations for a better performance in inserting the chains.
int TryRemove()
Trial removal.
int pNNano()
Number of nanoparticles.
double VetV(const Vettore *u, const Vettore *v)
Vectorial product between two Vettore returns the area.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
void RemChFromSys(int c)
Delete a chain in the system.
void VelVerletRigid2()
Velocity Verlet for a rigid body, second step.
void SwapChain(int c1, int c2, int b)
Swap two chains.
Definition: VarDataEl.cpp:355
int TryRemoveCh()
Trial removal of a chain.
double pBeta()
Beta factor 1/kTB.
Definition: VarData.h:914
Matrice computes the algebric operations on matrices.
void SolveLeaves()
Solve a system of four oder differential equation of particles in a line.
bool Set(int row, int column, double Val)
Set a coefficient.
double ChemPotId
Chemical potential of the particles.
Definition: Forces.h:721
void Apply(double *Known, double *UnKnown)
Apply Ax = y.
double Gaussiano(double Media, double Scarto)
Gaussian random number.
double CalcNrgCh(int c, double *Pot)
Pointer to the chain energy function.
Definition: Forces.h:513
int ReSetNChain(int NewNChain)
Realloc the number of chains.
Definition: Forces.cpp:877
void VelVerlet2()
Second step of the velocity Verlet.
int SwapPart(int p1, double *Pos1, int p2, double *Pos2)
Swap two particles.
Definition: Cubo.cpp:342
double Force[3]
Forces.
Definition: VarData.h:433
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
Quaternion class.
double AMom[3]
Angular momentum.
Definition: VarData.h:437
double pTemp()
Temperature.
Definition: VarData.h:912
int pNChain()
Number of chain.
double Angle(Vettore *u, Vettore *v)
Computes the angle between two Vetttore.
int CLat1
lateral coordinate
Definition: VarData.h:1076
double DensFuncNrgGhost(double *Pos, int p1, int t1)
Calculation of the energy from the density functional Hamiltonian for the ghost particle.
void WidomRemove(double *NrgDiff, int p)
Widom removal.
double Ext[3]
External.
Definition: Forces.h:182
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
void SolveRod()
Solve a system of four oder differential equation of particles connected by links.
Coefficient of a spline.
void BerendsenTherm()
Berendsen thermostat.
double AVel[3]
Angular velocity.
Definition: VarData.h:441
double Axis[3]
Rotation axis.
Definition: VarData.h:435
void ReInsertCh(int c)
Reinsert the chain in the previous position.
void MovePart(const int p, double *OldPos, double *NewPos)
Move a particle form the cell c1 to the cell c2.
Definition: Cubo.cpp:403
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
void SaveCh(int c)
Save the chain configuration.
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd.
void Print()
Print the entries.
double pDeltat()
Delta t.
Definition: VarData.h:910
DomDec * Pc
Pair list.
Definition: Forces.h:785
double NrgPBead
Average energy per particle.
Definition: Forces.h:735
double WeightSetBond(int p, int t)
Weight a set of bonds for the configurational bias.
double InsertCh(int c)
Displace a chain in the box.
double SLap
Prefactor of the square laplacian.
Definition: Forces.h:149
void SolveLinks()
Solve a system of four oder differential equation of particles connected by links.
int TryInsert()
Trial insertion.
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 CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.
void WidomBiasChOut(double *Weight, int c)
Widom with Rosenbluth weight.