Allink  v0.1
VarDataCreate.cpp
1 /***********************************************************************
2 VarDataCreate: Functions that create an initial configuration system as read from Polymer.conf.
3 Copyright (C) 2008-2010 by Giovanni Marelli <sabeiro@virgilio.it>
4 
5 
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10 
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; If not, see <http://www.gnu.org/licenses/>.
17 ***********************************************************************/
18 #include "../include/VarData.h"
19 
20 static int Tries = 0;
21 
22 int VarData::DefSoft(char *nome2,char *ConfF){
23  /*---------Variables--------------------*/
24  char cSystem[512];
25  double ReOSigma = 5.;
26  int NChain = 1250;
27  // Dens = 5.;
28  double ChainPArea = 1./(0.01974*QUAD(ReOSigma));
29  double Thickness = 0.849*ReOSigma;
30  double Volume = 1.;
31  double Area =0.;
32  int *arch;
33  if(ReadConf(ConfF)) return 1;
34  /*----------Setting----------*/
35  Gen->NChain = 0;
36  Gen->NPart = 0;
37  Gen->Step = 0;
38  Gen->NBlock = 0;
39  //if(Gen->kappaBend <= 0.) Gen->kappaSpring = 3.*(Gen->NPCh-1)/(2.*QUAD(Gen->ReOverCutOff));
40  VAR_ADD_TYPE(SysType,VAR_EDGE);
41  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
42  sigma = 1./sqrt(pkSpr());
43  Mat->InizializzaGaussiano(sigma,100);
44  int AddNSoft = 0;
45  for(int s=0;s<NSoft;s++){
46  //------------------------16---------------------------
47  if (Gen->NPCh == 16){
48  ChainPArea = 2.6*Soft[s].Size[2];
49  Thickness = 0.93*ReOSigma;
50  }
51  //------------------------11------------------------
52  else if (Gen->NPCh == 11){
53  ChainPArea = 3.19*Soft[s].Size[2];
54  Thickness = 2.55;
55  }
56  //------------------------12------------------------
57  else if (Gen->NPCh == 12){
58  ChainPArea = 3.19*Soft[s].Size[2];
59  Thickness = 2.55;
60  }
61  //------------------------13------------------------
62  else if (Gen->NPCh == 13){
63  ChainPArea = 3.19*Soft[s].Size[2];
64  Thickness = 2.55;
65  }
66  //------------------------10------------------------
67  else if (Gen->NPCh == 10){
68  ChainPArea = 3.98*Soft[s].Size[2];
69  Thickness = 0.8*ReOSigma;
70  }
71  //------------------------14------------------------
72  else if (Gen->NPCh == 14){
73  ChainPArea = 3.65*Soft[s].Size[2];
74  Thickness = 0.5*ReOSigma;
75  }
76  //------------------------20------------------------
77  else if (Gen->NPCh==21||Gen->NPCh==20){
78  ChainPArea = 3.65*Soft[s].Size[2];
79  Thickness = 0.6*ReOSigma;
80  }
81  //------------------------32---------------------
82  else if (Gen->NPCh == 32){
83  ChainPArea = 2.6*Soft[s].Size[2];
84  Thickness = 3.5;//1.14*ReOSigma;
85  }
86  if(VAR_IF_TYPE(SysType,VAR_TWOTAILS)){
87  if (Gen->NPCh == 21 ||Gen->NPCh == 20 ){
88  ChainPArea = 3.65*Soft[s].Size[2];
89  Thickness = 0.6*ReOSigma;
90  }
91  if (Gen->NPCh == 11||Gen->NPCh == 12||Gen->NPCh == 13){
92  ChainPArea = 3.19*Soft[s].Size[2];
93  Thickness = 2.55;
94  }
95  }
96  //-------------------arch--------------------
97  arch = (int *)calloc(Gen->NPCh,sizeof(int));
98  for(int b=0;b<Block[0].Asym;b++){
99  arch[b] = 0;
100  }
101  for(int b=Block[0].Asym;b<Gen->NPCh;b++){
102  arch[b] = 1;
103  }
104  //--------------------Area---------------------
105  Volume = Gen->Edge[CLat1]*Gen->Edge[CLat2]*Thickness;
106  if(VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR)){
107  Area = Gen->Edge[CLat1]*Gen->Edge[CLat2];
108  }
109  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_COATING)){
110  Thickness = 4.0;//1.15;
111  Area = .7*DUE_PI*(Nano->Rad+Thickness)*Nano->Height;
112  if(VAR_IF_TYPE(Nano->Shape,SHAPE_SPH)){
113  Area = .3*2.*DUE_PI*SQR(Nano->Rad+Thickness);
114  if (Gen->NPCh == 32){
115  Thickness = 6.0;//1.15;
116  Area = .1*2.*DUE_PI*SQR(Nano->Rad+Thickness);
117  }
118  }
119  }
120  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_TUBE)){
121  Area = DUE_PI*(Soft[s].Size[0]+Thickness)*Soft[s].Size[1];
122  }
123  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_VESICLE)){
124  ChainPArea = 3.6/Soft[s].Size[2];
125  //Thickness = 3.5;
126  Area = DUE_PI*SQR(Soft[s].Size[0]+.5*Thickness);
127  Area += DUE_PI*SQR(Soft[s].Size[0]-.5*Thickness);
128  AddNSoft += 1;
129  }
130  for(int n=0;n<Gen->NNano;n++){
131  if(Nano[n].Shape == SHAPE_NONE) continue;
132  if(Nano[n].Shape == SHAPE_WALL) continue;
133  Area -= PI*SQR(Nano[n].Rad);
134  }
135  //------------------NChain-------------------
136  if(Nano->Shape == SHAPE_SPH){
137  Volume -= 4.*PI*CUB((Nano->Rad))/3.;
138  Nano->Height = 0.;
139  }
140  if(Nano->Shape == SHAPE_CYL || Nano->Shape == SHAPE_TILT)
141  Volume -= PI*QUAD(Nano->Rad)*Nano->Height;
142  if(Nano->Shape == SHAPE_CLUSTER){
143  Volume -= PI*QUAD(Nano->Rad)*Nano->Height;
144  }
145  if(VAR_IF_TYPE(Soft[s].Topology,VAR_DISTRIBUTED)){
146  Volume = Gen->Edge[CLat1]*Gen->Edge[CLat2]*Gen->Edge[CNorm];
147  NChain = (int)( Volume*Soft[s].Size[2]/Gen->NPCh);
148  }
149  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR)){
150  NChain = (int)(Area*ChainPArea);
151  }
152  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_OBSTACLE)){
153  Soft[s].NPCh = 2*Gen->NPCh;
154  NChain = (int)(Area*ChainPArea*Soft[s].Size[2]);
155  NChain = 200;
156  }
157  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_COATING)){
158  NChain = (int)( .53*Area*ChainPArea);
159  }
160  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_TUBE)){
161  NChain = (int)( .69*Area*ChainPArea);
162  }
163  else if(VAR_IF_TYPE(Soft[s].Topology,VAR_VESICLE)){
164  ChainPArea = 3.0/Soft[s].Size[2];
165  //Thickness = 3.5;
166  NChain = (int)(Area*ChainPArea);
167  }
168  //------------DefBlock---------------
169  Soft[s].NChain = NChain;
170  Soft[s].NPCh = Gen->NPCh;
171  if(VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR_PE)){
172  Soft[s].NPCh = Gen->NPCh-1;
173  }
174  if(VAR_IF_TYPE(Soft[s].Topology,VAR_OBSTACLE))
175  Soft[s].NPCh = 2*Gen->NPCh;
176  Soft[s].NPart = Soft[s].NPCh*NChain;
177  //Soft[s].NPart = NChain*Gen->NPCh;
178  Gen->NChain += Soft[s].NChain;
179  Gen->NPart += Soft[s].NPart;
180  }
181  Gen->NBlock = NSoft+AddNSoft;
182  //----------Nano---------------
183  for(int n=0;n<Gen->NNano;n++){
184  if(Nano[n].Shape == SHAPE_CLUSTER){
185  Gen->NBlock++;
186  Gen->NChain += 1;
187  Gen->NPart += Nano[n].NCircle*Nano[n].NHeight;
188  if(Gen->NLink < 6) Gen->NLink = 6;
189  }
190  }
191  if(NAddChain > 0) Gen->NBlock++;
192  for(int s=0;s<NSoft;s++)
193  if( VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR) )
194  if(NAddChol > 0) Gen->NBlock++;
195  if(NStuffing > 0) Gen->NBlock++;
196  if(NSolvent > 0) Gen->NBlock++;
197  //--------------Block---------------
198  int DiblockLim = Block[0].Asym;
199  Block = (BLOCK *)calloc(Gen->NBlock,sizeof(BLOCK));
200  Soft[0].InitIdx = 0;
201  for(int s=1;s<NSoft;s++){
202  Soft[s].InitIdx += Soft[s-1].NPart + Soft[s-1].InitIdx;
203  }
204  for(int s=0,s1=0;s<NSoft;s++,s1++){
205  Soft[s].EndIdx = Soft[s].InitIdx + Soft[s].NPart;
206  Block[s1].Asym = DiblockLim;
207  Block[s1].InitIdx = Soft[s].InitIdx;
208  Block[s1].NChain = Soft[s].NChain;
209  Block[s1].EndIdx = Soft[s].EndIdx;
210  Block[s1].NPCh = Soft[s].NPCh;
211  sprintf(Block[s1].Name,"LIPID%d",s);
212  if(VAR_IF_TYPE(SysType,VAR_TWOTAILS))
213  sprintf(Block[s1].Name,"TT%d",s);
214  if(VAR_IF_TYPE(Soft[s].Topology,VAR_VESICLE)){
215  int NLayerIn = (int)(Soft[s].NChain*SQR(Soft[s].Size[0])/(SQR(Soft[s].Size[0])+SQR(Soft[s].Size[0]+Thickness)));
216  int NLayerOut = Soft[s].NChain - NLayerIn;
217  Block[s1].NChain = NLayerIn;
218  Block[s1].EndIdx = Block[s1].InitIdx+Soft[s].NPCh*NLayerIn;
219  Block[s1+1].InitIdx = Block[s1].EndIdx;
220  Block[s1+1].EndIdx = Soft[s].EndIdx;
221  Block[s1+1].NChain = NLayerOut;
222  Block[s1+1].NPCh = Soft[s].NPCh;
223  sprintf(Block[s1].Name,"INNER%d",s);
224  sprintf(Block[s1+1].Name,"OUTER%d",s);
225  s1++;
226  }
227  }
228  for(int n=0,b=NSoft+AddNSoft;n<Gen->NNano;n++){
229  if(Nano[n].Shape == SHAPE_CLUSTER){
230  Block[b].InitIdx = b != 0 ? Block[b-1].EndIdx : 0 ;
231  Block[b].NChain = 1;
232  Block[b].EndIdx = Block[b].InitIdx + Nano[n].NCircle*Nano[n].NHeight;
233  Block[b].NPCh = Nano[n].NCircle*Nano[n].NHeight;
234  Block[b].NPart = Block[b].NPCh;
235  sprintf(Block[b].Name,"PEP%d",n);
236  Nano[n].nBlock = b;
237  b++;
238  }
239  }
240  //-------------Alloc---------------
241  Pm = (PART *)calloc(Gen->NPart,sizeof(PART));
242  Ln = (LINKS *)calloc(Gen->NPart,sizeof(LINKS));
243  if(Pm == NULL){printf("Non s'alloca\n"); return 1;}
244  for(int p=0;p<Gen->NPart;p++){
245  Ln[p].Link = (int *)calloc(Gen->NLink,sizeof(int));
246  if(Ln[p].Link == NULL){printf("Non s'alloca\n"); return 1;}
247  }
248  SysInfo(cSystem);
249  printf("%s\n",cSystem);
250  SysDef(cSystem);
251  printf("%s",cSystem);
252  if(!HeaderSoft(cSystem))
253  printf("%s",cSystem);
254  /*---------Creating--------*/
255  char ArchFile[60];
256  int NMonCluster = 0;
257  for(int s=0,b=0;s<NSoft;s++,b++){
258  CreateSoft(arch,Thickness,s);
259  if(VAR_IF_TYPE(Soft[s].Topology,VAR_VESICLE))
260  b++;
261  NMonCluster = Block[b].EndIdx;
262  }
263  for(int n=0;n<Gen->NNano;n++){
264  if(Nano[n].Shape == SHAPE_CLINKS){
265  FindNeighbours("CrossLinks.dat");
266  }
267  if(Nano[n].Shape != SHAPE_CLUSTER) continue;
268  sprintf(ArchFile,"Architecture%d.dat",n);
269  CreateProtein(n,NMonCluster);
270  NMonCluster += Nano[n].NCircle*Nano[n].NHeight;
271  }
272  SetCoeff();
273  double Norm2 = CUBE(pReOverCutOff()) / SQR(pNPCh());
274  double Norm3 = CUBE(SQR(pReOverCutOff()) / pNPCh());
275  MInt->Rescale(Norm2,2);
276  MInt->Rescale(Norm3,3);
277  Write(nome2);
278  AddStuffing(nome2,NStuffing,0);
279  AddChains(nome2,Thickness);
280  AddSolvent(nome2,NSolvent);
281  for(int s=0;s<NSoft;s++){
282  if( VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR) ){
283  AddCholesterol(nome2,Thickness,s);
284  }
285  }
286  return 0;
287 }
288 bool VarData::CreateSoft(int *arch,double Thickness,int s){
289  if( VAR_IF_TYPE(Soft[s].Topology,VAR_COATING) )
290  CreateCoating(arch,Thickness,s);
291  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_VESICLE) )
292  CreateVesicle(arch,Thickness,s);
293  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_TUBE) )
294  CreateTube(arch,Thickness,s);
295  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR) )
296  CreatePlanar(arch,Thickness,s);
297  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_PLANAR_PE) )
298  CreatePlanar(arch,Thickness,s);
299  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_OBSTACLE) )
300  CreateObstacle(arch,Thickness,s);
301  else if( VAR_IF_TYPE(Soft[s].Topology,VAR_DISTRIBUTED) ){
302  double sigma = 1./sqrt(pkSpr());
303  //sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
304  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;){
305  for(int d=0;d<3;d++)
306  Pm[p].Pos[d] = Mat->Casuale()*Gen->Edge[d];
307  int IfContinue = 1;
308  for(int i=1;i<Gen->NPCh;i++){
309  for(int d=0;d<3;d++)
310  Pm[p+i].Pos[d] = Pm[p+i-1].Pos[d]+Mat->Gaussiano(0.,sigma);
311  if(CheckNano(Pm[p+i].Pos,s)){IfContinue=0;break;}
312  }
313  if(IfContinue){
314  c++;
315  p += Gen->NPCh;
316  }
317  }
318  }
319  else {
320  printf("System topology not recognized\n");
321  return 1;
322  }
323  DefRest(arch,s);
324  printf("Efficency: %d Tries for %d Particle\n",Tries,Soft[s].NPart);
325  return 0;
326 }
328  int Values=32;
329  double dValues = 1./(double)(Values);
330  double **Plot = (double **)calloc(Values,sizeof(double));
331  Gen->NPart = Values*Values;
332  Gen->NChain = Gen->NPart;
333  Gen->NPCh = 1;
334  for(int d=0;d<3;d++){
335  Gen->Edge[d] = (double)Values;
336  }
337  Pm = (PART *)calloc(Gen->NPart,sizeof(PART));
338  Ch = (CHAIN *)calloc(Gen->NChain,sizeof(CHAIN));
339  for(int i=0;i<Values;i++){
340  Plot[i] = (double *)malloc(Values*sizeof(double));
341  for(int j=0;j<Values;j++){
342  Pm[i*Values+j].Pos[CLat1] = Ch[i*Values+j].Pos[CLat1] = (double) j;
343  Pm[i*Values+j].Pos[CLat2] = Ch[i*Values+j].Pos[CLat2] = (double) i;
344  Pm[i*Values+j].Pos[CNorm] = Ch[i*Values+j].Pos[CNorm] =
345  Gen->Edge[CNorm]*.5 + cos(DUE_PI*(i)*dValues) + cos(DUE_PI*(i)*dValues*2.) + cos(DUE_PI*(i)*dValues*4.) +
346  cos(DUE_PI*(j)*dValues) + cos(DUE_PI*(j)*dValues*2.) + cos(DUE_PI*(j)*dValues*4.);
347  }
348  }
349  return 0;
350 }
351 void VarData::DefRest(int *arch,int s){
352  double sigma = 1./sqrt(pkSpr());
353  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;c++){
354  for(int i=0;i<Soft[s].NPCh;i++){
355  for(int d=0;d <3;d++){
356  Pm[p].Vel[d] = Mat->Gaussiano(0.,sigma) + Soft[s].Vel[d];
357  }
358  if(i == Soft[s].NPCh-1){
359  Ln[p].Link[0] = p-1;
360  Ln[p].NLink = 1;
361  }
362  else if(i == 0){
363  Ln[p].Link[0] = p+1;
364  Ln[p].NLink = 1;
365  }
366  else{
367  Ln[p].NLink = 2;
368  Ln[p].Link[0] = p-1;
369  Ln[p].Link[1] = p+1;
370  }
371  if( VAR_IF_TYPE(Soft[s].Topology,VAR_ADDED) )
372  Pm[p].Typ = 0;
373  else
374  Pm[p].Typ = arch[i];
375  Pm[p].CId = c;
376  Pm[p].Idx = p;
377  p++;
378  }
379  }
380 }
381 int VarData::PutPart(int j,int p,int HalfLim,double sigma){
382  int i=0;
383  if( j <= Block[0].Asym){
384  i = Block[0].Asym - j;
385  for(int d=0;d <3;d++){
386  Pm[p+i].Pos[d] = Pm[p+1+i].Pos[d]+Mat->Gaussiano(0.,sigma);
387  }
388  if( VAR_IF_TYPE(SysType,VAR_TWOTAILS) ){
389  if(i==HalfLim-1){
390  for(int d=0;d <3;d++)
391  Pm[p+i].Pos[d] = Pm[p+Block[0].Asym].Pos[d]+Mat->Gaussiano(0.,sigma);
392  }
393  }
394  }
395  else if(j>Block[0].Asym){
396  i = j;
397  for(int d=0;d <3;d++){
398  Pm[p+i].Pos[d] = Pm[p-1+i].Pos[d]+Mat->Gaussiano(0.,sigma);
399  }
400  }
401  return i;
402 }
403 void VarData::CreateObstacle(int *arch,double Thickness,int s){
404  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
405  int HalfLim = (int)(Block[0].Asym*.5);
406  double Dz = Thickness/16.;
407  double Leaflet = -.5*Thickness - 2.*Dz;
408  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;c++,p+=Soft[s].NPCh){
409  double Cas1 = Mat->Casuale(),Cas2 = Mat->Casuale();
410  Pm[p].Pos[CLat1] = Cas1*Gen->Edge[CLat1]*.5;
411  Pm[p].Pos[CLat2] = Cas2*Gen->Edge[CLat2]*.5;
412  Pm[p].Pos[CNorm] = Soft[s].Pos[CNorm]+Leaflet;
413  Pm[p].Typ = 1;
414  for(int pn=1;pn<Soft[s].NPCh;pn++){
415  Pm[p+pn].Pos[CLat1] = Pm[p].Pos[CLat1];
416  Pm[p+pn].Pos[CLat2] = Pm[p].Pos[CLat2];
417  Pm[p+pn].Pos[CNorm] = Pm[p+pn-1].Pos[CNorm] + Dz;
418  Pm[p+pn].Typ = 0;
419  if( (pn < 2) || (pn >= Soft[s].NPCh - 2)){
420  Pm[p+pn].Typ = 1;
421  }
422  }
423  }
424 }
425 void VarData::CreatePlanar(int *arch,double Thickness,int s){
426  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
427  int HalfLim = (int)(Block[0].Asym*.5);
428  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;){
429  printf("%d %d %d %lf \r",p,c,Tries,p/(double)Soft[s].NPart);
430  double Leaflet = -.5*Thickness;
431  if(c >= Soft[s].NChain/2) Leaflet = .5*Thickness;
432  //first part
433  double Cas1 = Mat->Casuale(),Cas2 = Mat->Casuale();
434  Pm[p+Block[0].Asym].Pos[CLat1] = Cas1*Gen->Edge[CLat1];
435  Pm[p+Block[0].Asym].Pos[CLat2] = Cas2*Gen->Edge[CLat2];
436  Pm[p+Block[0].Asym].Pos[CNorm] = Soft[s].Pos[CNorm]+Leaflet;
437  int IfContinue = 1;
438  //others
439  for(int j=1;j<Soft[s].NPCh;j++){
440  int i = PutPart(j,p,HalfLim,sigma);
441  //absorbing boundary condition
442  if(arch[i] == 0 ){
443  if(Pm[p+i].Pos[CNorm] > Soft[s].Pos[CNorm]+.5*Thickness ||
444  Pm[p+i].Pos[CNorm] < Soft[s].Pos[CNorm]-.5*Thickness ){
445  Tries++;
446  IfContinue = 0;
447  break;
448  }
449  }
450  else if (arch[i] == 1 ){
451  if(Pm[p+i].Pos[CNorm] < Soft[s].Pos[CNorm]+.5*Thickness &&
452  Pm[p+i].Pos[CNorm] > Soft[s].Pos[CNorm]-.5*Thickness ){
453  Tries++;
454  IfContinue = 0;
455  break;
456  }
457  }
458  if(CheckNano(Pm[p+i].Pos,s)){IfContinue=0;break;}
459  }
460  if(IfContinue){
461  c++;
462  p += Soft[s].NPCh;
463  }
464  }
465 }
466 void VarData::CreateVesicle(int *arch,double Thickness,int s){
467  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
468  int NLayerIn = (int)(Soft[s].NChain*SQR(Soft[s].Size[0])/(SQR(Soft[s].Size[0])+SQR(Soft[s].Size[0]+Thickness)) );
469  int NLayerOut = Soft[s].NChain - NLayerIn;
470  int HalfLim = (int)(Block[0].Asym*.5);
471  double inc = M_PI * (3. - sqrt(5.));
472  double NInv = 1. / (double)NLayerIn;
473  double Leaflet = -.5*Thickness;
474  for(int c=0,cc=0,p=Soft[s].InitIdx;c<Soft[s].NChain;){
475  printf("%d %d %d %lf \r",p,c,Tries,p/(double)Soft[s].NPart);
476  if(c == NLayerIn){
477  cc=0;
478  NInv = 1./(double)NLayerOut;
479  Leaflet = .5*Thickness;
480  }
481  double x = cc*2.*NInv - 1. + (NInv);
482  double r = sqrt(1.-x*x);
483  double phi = cc*inc;
484  double y = cos(phi)*r;
485  double z = sin(phi)*r;
486  //first part
487  Pm[p+Block[0].Asym].Pos[CLat1] =
488  (Soft[s].Size[0]+Leaflet)*x + Soft[s].Pos[CLat1];
489  Pm[p+Block[0].Asym].Pos[CLat2] =
490  (Soft[s].Size[0]+Leaflet)*y + Soft[s].Pos[CLat2];
491  Pm[p+Block[0].Asym].Pos[CNorm] =
492  (Soft[s].Size[0]+Leaflet)*z + Soft[s].Pos[CNorm];
493  //others
494  int IfContinue = 1;
495  for(int j=1;j<Gen->NPCh;j++){
496  int i = PutPart(j,p,HalfLim,sigma);
497  // absorbing boundary condition
498  double Dist = SQR(Pm[p+i].Pos[CLat1]-Soft[s].Pos[CLat1]) + SQR(Pm[p+i].Pos[CLat2]-Soft[s].Pos[CLat2]) + SQR(Pm[p+i].Pos[CNorm]-Soft[s].Pos[CNorm]);
499  if(arch[i] == 0){
500  if( Dist < SQR(Soft[s].Size[0]-.5*Thickness) || Dist > SQR(Soft[s].Size[0]+.5*Thickness) ) {
501  Tries++;
502  IfContinue = 0;
503  break;
504  }
505  }
506  else if(arch[i] == 1){
507  if(Dist < SQR(Soft[s].Size[0]+.5*Thickness) && Dist > SQR(Soft[s].Size[0] - .5*Thickness) ){
508  Tries++;
509  IfContinue = 0;
510  break;
511  }
512  }
513  if(CheckNano(Pm[p+i].Pos,s)){IfContinue=0;break;}
514  }
515  if(IfContinue){
516  c++;
517  cc++;
518  p += Gen->NPCh;
519  }
520  }
521  // for(int p=0;p<pNPart();p++){
522  // pPos(p);
523  // }
524 }
525 void VarData::CreateCoating(int *arch,double Thickness,int s){
526  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
527  int HalfLim = (int)(Block[0].Asym*.5);
528  double NInv = 1./(double)Soft[s].NChain;
529  double inc = 3.141592654 * (3. - sqrt(5.));
530  double Dist = 0.;
531  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;){
532  printf("%d %d %d %lf \r",p,c,Tries,p/(double)Soft[s].NPart);
533  //first part
534  double Cas1 = Mat->Casuale()*DUE_PI;
535  if(VAR_IF_TYPE(Nano->Shape,SHAPE_HEI)){
536  Pm[p+Block[0].Asym].Pos[CLat1] =
537  (Nano->Rad+.5*Thickness)*cos(Cas1) + Nano->Pos[CLat1];
538  Pm[p+Block[0].Asym].Pos[CLat2] =
539  (Nano->Rad+.5*Thickness)*sin(Cas1) + Nano->Pos[CLat2];
540  Pm[p+Block[0].Asym].Pos[CNorm] =
541  (.5 - Mat->Casuale())*Nano->Height + Nano->Pos[CNorm];
542  }
543  else if(VAR_IF_TYPE(Nano->Shape,SHAPE_SPH)){
544  double x = c*2.*NInv - 1. + (NInv);
545  double r = sqrt(1.-x*x);
546  double y = cos(c*inc)*r;
547  double z = sin(c*inc)*r;
548  Pm[p+Block[0].Asym].Pos[CLat1] =
549  (Nano->Rad+.5*Thickness)*x + Nano->Pos[CLat1];
550  Pm[p+Block[0].Asym].Pos[CLat2] =
551  (Nano->Rad+.5*Thickness)*y + Nano->Pos[CLat2];
552  Pm[p+Block[0].Asym].Pos[CNorm] =
553  (Nano->Rad+.5*Thickness)*z + Nano->Pos[CNorm];
554  }
555  int IfContinue = 1;
556  //others
557  for(int j=1;j<Gen->NPCh;j++){
558  int i = PutPart(j,p,HalfLim,sigma);
559  // absorbing boundary condition
560  if(VAR_IF_TYPE(Nano->Shape,SHAPE_HEI))
561  Dist = SQR(Pm[p+i].Pos[CLat1]-Nano->Pos[CLat1]) + SQR(Pm[p+i].Pos[CLat2]-Nano->Pos[CLat2]);
562  else if(VAR_IF_TYPE(Nano->Shape,SHAPE_SPH))
563  Dist = SQR(Pm[p+i].Pos[CLat1]-Nano->Pos[CLat1]) + SQR(Pm[p+i].Pos[CLat2]-Nano->Pos[CLat2]) + SQR(Pm[p+i].Pos[CNorm]-Nano->Pos[CNorm]);
564  if(arch[i] == 0 && Dist > SQR(Nano->Rad+.5*Thickness) ){
565  Tries++;
566  IfContinue = 0;
567  break;
568  }
569  else if(arch[i] == 1 && Dist < SQR(Nano->Rad+.5*Thickness) ){
570  Tries++;
571  IfContinue = 0;
572  break;
573  }
574  if(CheckNano(Pm[p+i].Pos,s)){IfContinue=0;break;}
575  }
576  if(IfContinue){
577  c++;
578  p += Gen->NPCh;
579  }
580  }
581 }
582 void VarData::CreateTube(int *arch,double Thickness,int s){
583  double sigma = sqrt(QUAD(Gen->ReOverCutOff)/(Gen->NPCh-1)/3.);
584  int HalfLim = (int)(Block[0].Asym*.5);
585  int NLayerIn = (int)(Soft[s].NChain/2.*SQR(Soft[s].Size[0])/(SQR(Soft[s].Size[0])+SQR(Soft[s].Size[0]+Thickness)) );
586  for(int c=0,p=Soft[s].InitIdx;c<Soft[s].NChain;){
587  printf("%d %d %d %lf \r",p,c,Tries,p/(double)Soft[s].NPart);
588  double Leaflet = -.5*Thickness;
589  if(c >= NLayerIn) Leaflet = .5*Thickness;
590  //first part
591  double Cas1 = DUE_PI*Mat->Casuale();
592  //FIXME: this distribution is somehow not uniform
593  Pm[p+Block[0].Asym].Pos[CLat1] =
594  (Soft[s].Size[0]+Leaflet)*cos(Cas1) + Soft[s].Pos[CLat1];
595  Pm[p+Block[0].Asym].Pos[CLat2] =
596  (Soft[s].Size[0]+Leaflet)*sin(Cas1) + Soft[s].Pos[CLat2];
597  Pm[p+Block[0].Asym].Pos[CNorm] = (Mat->Casuale() - .5)*Soft[s].Size[1] + Soft[s].Pos[CNorm];
598 
599  //others
600  int IfContinue = 1;
601  for(int j=1;j<Gen->NPCh;j++){
602  int i = PutPart(j,p,HalfLim,sigma);
603  // absorbing boundary condition
604  double Dist = SQR(Pm[p+i].Pos[CLat1]-Soft[s].Pos[CLat1]) + SQR(Pm[p+i].Pos[CLat2]-Soft[s].Pos[CLat2]);
605  if(arch[i] == 0)
606  if( Dist < SQR(Soft[s].Size[0]-.5*Thickness) || Dist > SQR(Soft[s].Size[0]+.5*Thickness) ) {
607  Tries++;
608  IfContinue = 0;
609  break;
610  }
611  else if(arch[i] == 1)
612  if(Dist > SQR(Soft[s].Size[0]-.5*Thickness) || Dist < SQR(Soft[s].Size[0] + .5*Thickness) ){
613  Tries++;
614  IfContinue = 0;
615  break;
616  }
617  if(CheckNano(Pm[p+i].Pos,s)){IfContinue=0;break;}
618  }
619  if(IfContinue){
620  c++;
621  p += Gen->NPCh;
622  }
623  }
624 }
625 int VarData::CheckNano(double *Pos,int s){
626  for(int n=0;n<pNNano();n++){
627  Point2Shape(Nano[n].Shape);
628  double Add = .0;
629  double Radius2 = 0.;
630  if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_BOUND)) Add = .3;
631  if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_NONE)) continue;
632  else{
633  Radius2 = NanoDist2(Pos,n);
634  }
635  // else{
636  // Vettore PosRel(3);
637  // Vettore NanoAxis(3);
638  // Vettore Distance(3);
639  // for(int d=0;d<3;d++){
640  // PosRel.Set(Nano[n].Pos[d] - Pos[d],d);
641  // NanoAxis.Set(Nano[n].Axis[d],d);
642  // }
643  // Distance.VetV(&NanoAxis,&PosRel);
644  // Radius2 = SQR(Distance.Norm());
645  // double HeiOnAxis = PosRel.ProjOnAxis(&NanoAxis);
646  // if( fabs(HeiOnAxis) > Nano[n].Height*.5)
647  // Radius2 = QUAD(Nano[n].Rad+Add);
648  // }
649  if(Radius2 < QUAD(Nano[n].Rad+Add)){
650  Tries++;
651  return 1;
652  }
653  }
654  return 0;
655 }
656 //Obsolete
657 void VarData::AddProtein(int NCircle,int NHeight,int nNano,char *filename){
658  int NPart = NCircle*NHeight;
659  PART *Pn = (PART *)calloc(NPart,sizeof(PART));
660  double CirInv = 1./(double)NCircle;
661  double HeiInv = Nano[nNano].Height/(double)NHeight;
662  for(int c=0;c<NCircle;c++){
663  double Sin = sin(c*CirInv*DUE_PI);
664  double Cos = cos(c*CirInv*DUE_PI);
665  double x = Nano[nNano].Rad * Cos + Nano[nNano].Pos[0];
666  double y = Nano[nNano].Rad * Sin + Nano[nNano].Pos[1];
667  for(int h=0;h<NHeight;h++){
668  int p = c*NHeight + h;
669  double z = h*HeiInv + Nano[nNano].Pos[2] - .5*Nano[nNano].Height;
670  Pn[p].Pos[CLat1] = x;
671  Pn[p].Pos[CLat2] = y;
672  Pn[p].Pos[CNorm] = z;
673  }
674  }
675  FILE *ReSave = fopen(filename,"a");
676  fprintf(ReSave,"# n=1 N=%d name=PEP1\n",NPart);
677  for(int p=0;p<NPart;p++)
678  fprintf(ReSave,"%lf %lf %lf %lf %lf %lf %d\n",
679  Pn[p].Pos[0],Pn[p].Pos[1],Pn[p].Pos[2],
680  0.0,0.0,0.0,0);
681  Nano[nNano].OffSet = (double)NCircle;
682 }
683 void VarData::CreateProtein(int nNano,int np){
684  int NCircle = Nano[nNano].NCircle;
685  int NHeight = Nano[nNano].NHeight;
686  int NCyl = NCircle*NHeight;
687  int NSph = (NCircle)*(NCircle/2-1) + 2;
688  int NTot = NCyl;// + NSph;
689  int IfDoubleSided = 0;
690  int IfKink = 0;
691  double AsymPhil = -.1*Nano[nNano].Rad;
692  double CirInv = 1./(double)NCircle;
693  double HeiInv = Nano[nNano].Height/(double)(NHeight-1);
694  double HeiSph = Nano[nNano].Height;
695  if(Nano[nNano].NHeight == 0 ) HeiSph + HeiInv;
696  int pWrote = 0;
697  double Shift = .5*sin(1*CirInv*DUE_PI);
698  double SegSide = Nano[nNano].Height/(double)Nano[nNano].NHeight;
699  double SegCirc = DUE_PI*Nano[nNano].Rad/(double)Nano[nNano].NCircle;
700  // Part positions
701  // Cylinder
702  for(int c=0;c<NCircle;c++){
703  for(int h=0;h<NHeight;h++){
704  int NLink = 0;
705  int p = c*NHeight + h;
706  //pos current particle
707  double Sin = sin(c*CirInv*DUE_PI);
708  double Cos = cos(c*CirInv*DUE_PI);
709  double Rad = Nano[nNano].Rad;
710  double Weight = 1.;
711  double Axes[3] = {pNanoPos(nNano,0),pNanoPos(nNano,1),pNanoPos(nNano,2)};
712  if(IfKink){
713  Axes[0] += .25*Nano->Height*(fabs((double)(h-NHeight/2))/(double)(NHeight));
714  // Weight -= .1*Cos;
715  // Weight -= .8*(1.-fabs((double)(h-NHeight/2))/(double)(NHeight));
716  }
717  if( ((h)%2)==0 ){
718  Sin = sin((c+.5)*CirInv*DUE_PI);
719  Cos = cos((c+.5)*CirInv*DUE_PI);
720  Rad = Nano[nNano].Rad - .1;
721  }
722  Ln[p+np].NLink = 0;
723  Pm[p+np].Pos[0] = Rad * Cos + Axes[0];
724  Pm[p+np].Pos[1] = Rad * Sin + Axes[1];
725  Pm[p+np].Pos[2] = h*HeiInv*Weight - .5*Nano[nNano].Height*Weight + Axes[2];
726  Pm[p+np].Typ = 0;
727  if(IfDoubleSided){
728  if(Pm[p+np].Pos[0] < (Nano[nNano].Pos[0] + AsymPhil))
729  Pm[p+np].Typ = 1;
730  }
731  if(h < 2 || h > NHeight - 3) Pm[p+np].Typ = 1;
732  pWrote++;
733  //-------------Connections------------------
734  // Right
735  int pp = (c+1)*NHeight + h;
736  if( pp >= NCyl ) pp -= NCyl;
737  Ln[p+np].Link[NLink++] = pp + np;
738  pp = p + 1;
739  if( (pp%NHeight)!=0 ) Ln[p+np].Link[NLink++] = pp + np;
740  //Left up
741  pp = (c-1)*NHeight + h + 1;
742  if(c==0) pp = (NCircle-1)*NHeight + h + 1;
743  if( ((h)%2)==0 ){
744  pp = (c+1)*NHeight + h + 1;
745  if(c==NCircle-1) pp = 0 + h + 1;
746  }
747  if( (pp%Nano[nNano].NHeight) ) Ln[p+np].Link[NLink++] = pp + np;
748  // Side
749  if( h == 0 ) {
750  pp = (c)*NHeight + NHeight - 2;
751  //if(pp < NCyl) Ln[p+np].Link[NLink++] = pp + np;
752  // Diagonal
753  pp = (c+NCircle/2)*NHeight + NHeight-2;//p + NCyl/2 + NCircle - 1;
754  if(pp > NCyl) pp = (c-NCircle/2)*NHeight + NHeight-2;
755  //if(pp < NCyl) Ln[p+np].Link[NLink++] = pp + np;
756  }
757  // Disks
758  pp = p + NCircle/2*NHeight;
759  //if(pp < NCyl-1) Ln[p+np].Link[NLink++] = pp + np;
760  Ln[p+np].NLink = NLink;
761  }
762  }
763  Vettore Axis(0.,0.,1.);
764  Vettore Axis1(Nano[nNano].Axis[0],Nano[nNano].Axis[1],Nano[nNano].Axis[2]);
765  Vettore Axis2(3);
766  Axis2 = Axis1 + Axis;
767  for(int d=0;d<3;d++){
768  Axis2.Set(Axis1.Val(d)+Axis.Val(d),d);
769  }
770  Vettore Origin(Nano[nNano].Pos[0],Nano[nNano].Pos[1],Nano[nNano].Pos[2]);
771  int b = nNano + NSoft;
772  RotateBlock(&Axis2,&Origin,b);
773  char Filename[60];
774  sprintf(Filename,"Architecture%d.dat",nNano);
775  FILE *CSave = fopen(Filename,"w");
776  fprintf(CSave,"# Cylinder\n");
777  //write the initial mutual distances
778  for(int p=np;p<NTot+np;p++){
779  for(int l=0;l<Ln[p].NLink;l++){
780  int l2 = Ln[p].Link[l];
781  double Dist = sqrt( SQR(Pm[p].Pos[0]-Pm[l2].Pos[0]) + SQR(Pm[p].Pos[1]-Pm[l2].Pos[1]) + SQR(Pm[p].Pos[2]-Pm[l2].Pos[2]) );
782  double kSpr = 10000.;
783  if(Dist > 2.) kSpr = 10000.;
784  if(Dist > 4.) kSpr = 10000.;
785  fprintf(CSave,"%d %d %lf %.0f\n",p-np,l2-np,Dist,kSpr);
786  }
787  }
788  fclose(CSave);return;
789  // FIXME: the links are not closing at the boundaries
790  // Cupola
791  fprintf(CSave,"# Cupola\n");
792  int PType = 2;
793  int NCircHalf = NCircle/2;
794  // Vertical
795  for(int cc=1;cc<NCircHalf-1;cc++){
796  double Sin2 = sin(cc*CirInv*DUE_PI);
797  double Cos2 = cos(cc*CirInv*DUE_PI);
798  // Horizontal
799  for(int c=0;c<NCircle;c++){
800  double Sin = sin(c*CirInv*DUE_PI);
801  double Cos = cos(c*CirInv*DUE_PI);
802  if( (cc%2)==0 ){
803  Sin = sin((c+.5)*CirInv*DUE_PI);
804  Cos = cos((c+.5)*CirInv*DUE_PI);
805  }
806  double x = Nano[nNano].Rad * Cos * Sin2 + Nano[nNano].Pos[0];
807  double y = Nano[nNano].Rad * Sin * Sin2 + Nano[nNano].Pos[1];
808  double Quota = Nano[nNano].Height*.5;
809  if(cc > NCircle/4) Quota = - Nano[nNano].Height*.5;
810  double z = Nano[nNano].Rad * Cos2 + Nano[nNano].Pos[2] + Quota;
811  pWrote++;
812  //fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",x,y,z,0.,0.,0.,PType);
813  // Connections
814  double Dx = x - Nano[nNano].Rad*cos((c+1)*CirInv*DUE_PI)*Sin2 - Nano[nNano].Pos[0];
815  double Dy = y - Nano[nNano].Rad*sin((c+1)*CirInv*DUE_PI)*Sin2 - Nano[nNano].Pos[1];
816  double Elong = sqrt(SQR(Dx)+SQR(Dy));
817  int p = NCyl + c + cc*NCircle;
818  int pp = NCyl + (c+1) + cc*NCircle;
819  if( c+1==NCircle ) pp = NCyl + (0) + cc*NCircle;
820  if(c != 0 && c != NCircle-1) fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
821  if( cc != NCircle/2-1 && cc != NCircle/4){
822  pp = NCyl + (c+1) + (cc+1)*NCircle;
823  Elong = sqrt( SQR(Nano[nNano].Rad*1.*CirInv*DUE_PI) + SQR(.5*Elong));
824  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
825  pp = NCyl + (c) + (cc-1)*NCircle;
826  if( (cc%2)==0 )
827  pp = NCyl + (c+1) + (cc)*NCircle;
828  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
829  }
830  }
831  }
832  // Closing the extrema
833  fprintf(CSave,"# Extrema\n");
834  {
835  double x = Nano[nNano].Pos[0];
836  double y = Nano[nNano].Pos[1];
837  double z = Nano[nNano].Pos[2] - HeiSph*.5 - Nano[nNano].Rad;
838  pWrote++;
839  //fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",x,y,z,0.,0.,0.,PType);
840  z = Nano[nNano].Pos[2] + HeiSph*.5 + Nano[nNano].Rad;
841  pWrote++;
842  //fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",x,y,z,0.,0.,0.,PType);
843  {
844  int p = pWrote - 1;
845  int pp = pWrote - 2;
846  double Elong = Nano[nNano].Height + 2.*Nano[nNano].Rad;
847  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
848  }
849  for(int c=0;c<NCircle;c++){
850  int p = pWrote-2;
851  int pp = c + NCyl + NCircle*(NCircle/2-2);
852  double Elong = Nano[nNano].Rad * 1.*CirInv*DUE_PI;
853  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
854  p = pWrote-1;
855  pp = c + NCyl;
856  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
857  // Conntecting to the cylinder
858  p = c*NHeight;
859  pp = c + NCyl + NCircle*(NCircle/4);
860  Elong = HeiInv;//Actually it is a bit more
861  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
862  p = c*(NHeight) + NHeight - 1;
863  pp = c + NCyl + NCircle*(NCircle/4-1);
864  fprintf(CSave,"%d %d %lf\n",p,pp,Elong);
865  }
866  }
867  fclose(CSave);
868  //fclose(PSave);
869 }
870 void VarData::AddStuffing(char *filename,int NWater,int nNano){
871  if(NWater == 0) return;
872  FILE *PSave = fopen(filename,"a");
873  fprintf(PSave,"# n=%d N=10 name=STUFFING\n",NWater/10);
874  for(int p=0;p<NWater;p++){
875  double Angle = Mat->Casuale()*DUE_PI;
876  double Rad = 2.*(Mat->Casuale()-.5)*Nano[nNano].Rad;
877  double x = cos(Angle)*Rad + Nano[nNano].Pos[0];
878  double y = sin(Angle)*Rad + Nano[nNano].Pos[1];
879  double z = (Mat->Casuale()-.5)*Nano[nNano].Height + Nano[nNano].Pos[2];
880  fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",x,y,z,0.,0.,0.,1);
881  }
882  fclose(PSave);
883 }
884 void VarData::AddSolvent(char *filename,int NWater){
885  if(NSolvent <= 0) return;
886  FILE *PSave = fopen(filename,"a");
887  double sigma = 1./sqrt(pkSpr());
888  fprintf(PSave,"# n=%d N=1 name=SOLVENT\n",NWater);
889  for(int p=0;p<NWater;p++){
890  double x = Mat->Casuale()*Gen->Edge[CLat1];
891  double y = Mat->Casuale()*Gen->Edge[CLat2];
892  double z = Mat->Casuale()*Gen->Edge[CNorm]*.3;
893  fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",x,y,z,Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),2);
894  }
895  fclose(PSave);
896 }
897 void VarData::AddChains(char *filename,double Thickness){
898  if(NAddChain <= 0) return;
899  double sigma = 1./sqrt(pkSpr());
900  PART *Pn = (PART *)calloc(Gen->NPCh,sizeof(PART));
901  FILE *PSave = fopen(filename,"a");
902  int NPCh = (int)(.5*pNPCh());
903  fprintf(PSave,"# n=%d N=%d name=ADDED\n",NAddChain,NPCh);
904  int s = 0;
905  for(int c=0;c<NAddChain;){
906  Pn[0].Pos[CLat1] = Mat->Casuale()*Gen->Edge[CLat1];
907  Pn[0].Pos[CLat2] = Mat->Casuale()*Gen->Edge[CLat2];
908  Pn[0].Pos[CNorm] = .2*(Mat->Casuale()-.5)*Thickness + Soft[s].Pos[CNorm];
909  int IfContinue = !CheckNano(Pn[0].Pos,s);
910  for(int i=1;i<NPCh;i++){
911  for(int d=0;d<3;d++)
912  Pn[i].Pos[d] = Pn[i-1].Pos[d]+Mat->Gaussiano(0.,sigma);
913  if(Pn[i].Pos[CNorm] > Soft[s].Pos[CNorm]+.5*Thickness ||
914  Pn[i].Pos[CNorm] < Soft[s].Pos[CNorm]-.5*Thickness ){
915  Tries++;
916  IfContinue = 0;
917  break;
918  }
919  //if(CheckNano(Pn[i].Pos,0)){IfContinue=0;break;}
920  }
921  if(IfContinue){
922  for(int i=0;i<NPCh;i++){
923  fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",Pn[i].Pos[0],Pn[i].Pos[1],Pn[i].Pos[2],Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),0);
924  }
925  c++;
926  s++;if(s==NSoft)s=0;
927  }
928  }
929  fclose(PSave);
930  free(Pn);
931 }
932 void VarData::AddCholesterol(char *filename,double Thickness,int s){
933  if(NAddChol <= 0) return;
934  double sigma = 1./sqrt(pkSpr());
935  int HalfLim = (int)(Block[0].Asym*.5);
936  int NPCh = 5;
937  int DLim = NPCh-1;
938  PART *Pn = (PART *)calloc(Gen->NPCh,sizeof(PART));
939  FILE *PSave = fopen(filename,"a");
940  fprintf(PSave,"# n=%d N=%d name=CHOL%d\n",NAddChol,NPCh,s);
941  for(int c=0;c<NAddChol;){
942  printf("%d %d %lf \r",c,Tries,c/(double)(NAddChol));
943  double Leaflet = -.5*Thickness;
944  if(c >= NAddChol/2) Leaflet = +.5*Thickness;
945  //first part
946  double Cas1 = Mat->Casuale(),Cas2 = Mat->Casuale();
947  Pn[DLim].Pos[CLat1] = Cas1*Gen->Edge[CLat1];
948  Pn[DLim].Pos[CLat2] = Cas2*Gen->Edge[CLat2];
949  Pn[DLim].Pos[CNorm] = Soft[s].Pos[CNorm]+Leaflet;
950  Pn[DLim].Typ = 1;
951  int IfContinue = !CheckNano(Pn[DLim].Pos,s);
952  //others
953  for(int j=DLim-1;j>=0;j--){
954  Pn[j].Typ = 0;
955  for(int d=0;d<3;d++)
956  Pn[j].Pos[d] = Pn[j+1].Pos[d]+Mat->Gaussiano(0.,sigma);
957  if(Pn[j].Pos[CNorm] > Soft[s].Pos[CNorm]+.5*Thickness ||
958  Pn[j].Pos[CNorm] < Soft[s].Pos[CNorm]-.5*Thickness ){
959  Tries++;
960  IfContinue = 0;
961  break;
962  }
963  if(CheckNano(Pn[j].Pos,s)){IfContinue=0;break;}
964  }
965  if(IfContinue){
966  for(int i=0;i<NPCh;i++){
967  fprintf(PSave,"%lf %lf %lf %lf %lf %lf %d\n",Pn[i].Pos[0],Pn[i].Pos[1],Pn[i].Pos[2],Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),Mat->Gaussiano(0.,sigma),Pn[i].Typ);
968  }
969  c++;
970  }
971  }
972  fclose(PSave);
973 }
974 #include <Cubo.h>
975 typedef DdLinkedList Cubo;
976 void VarData::FindNeighbours(char *FileName){
977  double CutOff = 1.0;
978  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
979  int NeiList[27];
980  int NPair = (int)(pNChain()/10.);
981  int NChOffSet = 0;
982  int bInner = 0;
983  for(int b=0,NCh=0;b<pNBlock();NCh+=Block[b++].NChain){
984  if(strcmp(Block[b].Name,"INNER0")) continue;
985  NChOffSet = NCh;
986  bInner = b;
987  }
988  int NChain = 0;
989  for(int b=0;b<pNBlock();b++){
990  NChain += Block[b].NChain;
991  }
992  SetNChain(NChain);
993  for(int c=0;c<pNChain();c++){
994  for(int d=0;d<3;d++){
995  Ch[c].Pos[d] = Pm[c*pNPCh()+pNPCh()-1].Pos[d];
996  }
997  // for(int p=c*pNPCh();p<(c+1)*pNPCh();p++){
998  // for(int d=0;d<3;d++){
999  // Ch[c].Pos[d] += Pm[p].Pos[d];
1000  // }
1001  // }
1002  // for(int d=0;d<3;d++){
1003  // Ch[c].Pos[d] /= (double)pNPCh();
1004  // }
1005  }
1006  double *cPair = (double *)calloc(3*pNChain(),sizeof(double));
1007  for(int c=0;c<pNChain();c++){
1008  cPair[c*3+0] = c;
1009  cPair[c*3+1] = c;
1010  cPair[c*3+2] = 1000.;
1011  }
1012  Cubo *Pc = new Cubo(Edge,pNChain(),CutOff);
1013  for(int c=0;c<pNChain();c++){
1014  int p = c*pNPCh()+pNPCh()-1;
1015  //Pc->AddPart(c,Pm[p].Pos);
1016  Pc->AddPart(c,Ch[c].Pos);
1017  }
1018  FILE *FWrite = fopen(FileName,"w");
1019  for(int c=NChOffSet;c<NChOffSet+Block[bInner].NChain-1;c++){
1020  cPair[c*3+0] = (double)c;
1021  int p1 = c*pNPCh(bInner)+pNPCh(bInner)-1;
1022  double MinDist = 1000.;
1023  int NNei = Pc->GetNei(Pm[p1].Pos,NeiList);
1024  for(int i=0;i<NNei;i++){
1025  int c1 = NeiList[i];
1026  for(Pc->SetCounters(c1);Pc->IfItCell(c1);Pc->IncrCurr(c1)){
1027  int c2 = Pc->ItCell(c1);
1028  if(c2 <= c) continue;
1029  int p2 = c2*pNPCh(bInner)+pNPCh(bInner)-1;
1030  double Dist2 = 0.;
1031  for(int d=0;d<3;d++){
1032  Dist2 += SQR(Ch[c].Pos[d] - Ch[c2].Pos[d]);
1033  //Dist2 += SQR(Pm[p1].Pos[d] - Pm[p2].Pos[2]);
1034  }
1035  if(MinDist > Dist2){
1036  MinDist = Dist2;
1037  cPair[c*3+1] = (double)c2;
1038  cPair[c*3+2] = MinDist;
1039  }
1040  }
1041  }
1042  }
1043  double Temp[3];
1044  for(int c=1;c<pNChain();c++){
1045  for(int c1=c;c1>=0;c1--){
1046  if(cPair[c1*3+2] >= cPair[(c1-1)*3+2]) break;
1047  for(int d=0;d<3;d++){
1048  Temp[d] = cPair[c1*3+d];
1049  cPair[c1*3+d] = cPair[(c1-1)*3+d];
1050  cPair[(c1-1)*3+d] = Temp[d];
1051  }
1052  // Mat->Swap(cPair,c1*3+2,cPair,(c1-1)*3+2,3);
1053  //printf("%d) %d %lf %lf\n",c,c1,cPair[(c1)*3+2],cPair[(c1-1)*3+2]);
1054  }
1055  }
1056  int pRef = 0;//NChOffSet*pNPCh(bInner);
1057  double KEl = 1.;
1058  for(int c=0;c<NPair;c++){
1059  int p1 = (int)cPair[c*3+0]*pNPCh(bInner)+pNPCh(bInner)-1;
1060  int p2 = (int)cPair[c*3+1]*pNPCh(bInner)+pNPCh(bInner)-1;
1061  double Dist = 0.;
1062  for(int d=0;d<3;d++){
1063  Dist += SQR(Ch[(int)cPair[c*3+0]].Pos[d] - Ch[(int)cPair[c*3+1]].Pos[d]);
1064  //Dist += SQR(Pm[p1].Pos[d] - Pm[p2].Pos[d]);
1065  }
1066  Dist = sqrt(Dist);
1067  //Dist = 0.;
1068  fprintf(FWrite,"%d %d %lf %lf\n",p1-pRef,p2-pRef,Dist,KEl);
1069  }
1070  fclose(FWrite);
1071 }
int CheckNano(double *Pos, int s)
No particle inside the nano.
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int CId
Chain Identifier.
Definition: VarData.h:224
void CreatePlanar(int *arch, double Thickness, int s)
planar membrane
void AddStuffing(char *filename, int nStuffing, int nNano)
Fill the protein with water.
void SetCounters(int c)
Set the counters to the initial position.
Definition: Cubo.cpp:413
void SysDef(char *cSystem)
Print a string with the system definitions.
Definition: VarData.cpp:167
double NanoDist2(double *Pos, int n)
Pointer to a generic function.
Definition: VarData.h:764
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void SysInfo(char *cSystem)
Print a string with the system information.
Definition: VarData.cpp:164
int NAddChol
Additional cholesterol chains into the membrane.
Definition: VarData.h:1068
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
NANO * Nano
Extra particle.
Definition: VarData.h:1044
int SetNChain(int NewNCh)
Set and reallocate the number of chains.
Definition: VarDataComm.cpp:78
void CreateVesicle(int *arch, double Thickness, int s)
vesicle
void AddSolvent(char *filename, int nWater)
Add phantom solvent at the bottom.
void CreateObstacle(int *arch, double Thickness, int s)
Creates obstacles.
void AddCholesterol(char *filename, double Thickness, int s)
Add cholesterol chains in the bilayer.
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
int Asym
Diblock limit of the chain.
Definition: VarData.h:273
bool InizializzaGaussiano(double Scarto, int N)
Initialize the Gaussian number generator.
double Rad
Size.
Definition: VarData.h:445
double Height
Height of the cylinder.
Definition: VarData.h:449
int CLat2
lateral coordinate
Definition: VarData.h:1078
Information for every block.
Definition: VarData.h:255
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
void AddProtein(int NCircle, int NHeight, int nNano, char *filename)
Defines the nanoparticle as a net of monomers.
void DefRest(int *arch, int s)
set the remaining information
int InitIdx
Initial Position.
Definition: VarData.h:294
int NLink
Maximum number of bonds.
Definition: VarData.h:357
double pkSpr()
Spring coupling.
Definition: VarData.h:934
double Val(int N)
Value of the N column.
int IfItCell(const int c)
Stop the loop and set the counter to zero.
Definition: Cubo.cpp:427
int NAddChain
Additional homopolymer chains into the membrane.
Definition: VarData.h:1066
Information of every chain.
Definition: VarData.h:236
int TrialSys()
Creates a trial system.
void CreateCoating(int *arch, double Thickness, int s)
coating around a cylindrical nanoparticle
void FindNeighbours(char *FileName)
Find the couples of most neighbouring chains.
int InitIdx
Initial particle position.
Definition: VarData.h:265
double Pos[3]
Position.
Definition: VarData.h:427
int Shape
0 none, 1 spherical, 2 cylindrical 3 wall
Definition: VarData.h:473
double pNanoPos(int n, int d)
Return back folded nano position.
int nBlock
In which block is the peptide written.
Definition: VarData.h:481
int NChain
chains
Definition: VarData.h:290
int NStuffing
Stuffing for the cylinder.
Definition: VarData.h:1072
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int pNPCh()
Number of particle per chain.
int Idx
Particle identifier.
Definition: VarData.h:222
int pNBlock()
Number of blocks.
void RotateBlock(Vettore *Axis, Vettore *Origin, int b)
Rotate a block wrt to the Axis from the Origin.
Definition: VarDataEl.cpp:444
MatInt * MInt
Matrix of the prefactor of the interactions.
Definition: VarData.h:529
int NChain
Number of chain.
Definition: VarData.h:347
void IncrCurr(const int c)
Increment the current part in the cell.
Definition: Cubo.cpp:435
double Vel[3]
bias velocity
Definition: VarData.h:282
double Casuale()
Random uniform number.
int NPart
particles
Definition: VarData.h:261
int NCircle
Number of monomers per circle.
Definition: VarData.h:479
void Point2Shape(int iShape)
Point to the shape function.
Definition: VarDataEl.cpp:490
int EndIdx
End position.
Definition: VarData.h:296
int NPCh
Number of particle per chain.
Definition: VarData.h:349
int SysType
Contains the definition of the system.
Definition: VarData.h:1086
double Pos[3]
initial position
Definition: VarData.h:280
int pNNano()
Number of nanoparticles.
int NPart
particles
Definition: VarData.h:288
int NPCh
part per chain
Definition: VarData.h:292
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int NPCh
particles per chain
Definition: VarData.h:259
double Gaussiano(double Media, double Scarto)
Gaussian random number.
void CreateTube(int *arch, double Thickness, int s)
Soft in a tube shape.
int GetNei(double *Pos, int *NeiList)
Choose among the different neighbouring lists.
Definition: Cubo.h:208
int ItCell(const int c)
Iterate in the cell.
Definition: Cubo.cpp:421
bool CreateSoft(int *arch, double Thickness, int s)
Creates an initial system.
int NSolvent
Solvent molecules.
Definition: VarData.h:1070
double OffSet
Reference potential.
Definition: VarData.h:457
double pReOverCutOff()
Re/CutOff.
Definition: VarData.h:946
SOFT * Soft
Soft bodies.
Definition: VarData.h:1052
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int NChain
chains
Definition: VarData.h:263
int EndIdx
End particle position.
Definition: VarData.h:267
unsigned long Step
Courrent step.
Definition: VarData.h:343
int HeaderSoft(char *Line)
Header soft.
int pNChain()
Number of chain.
int CLat1
lateral coordinate
Definition: VarData.h:1076
Information of every particle.
Definition: VarData.h:214
void AddChains(char *filename, double Thickness)
Add homopolymer chains in the bilayer.
int DefSoft(char *nome2, char *ConfF)
Define and write the system as described in the conf file.
bool ReadConf(char *InFile)
Reads a "configuration file".
Definition: VarDataRead.cpp:28
void CreateProtein(int nNano, int nStart)
Defines the nanoparticle as a net of monomers.
int NSoft
Number of soft bodies.
Definition: VarData.h:1062
int PutPart(int j, int p, int HalfLim, double sigma)
return the number in the chain of the next particle put
void Rescale(double SFactor, int Order)
Rescale entries.
Definition: VarData.cpp:312
int NPart
Number of particle.
Definition: VarData.h:345
Domain decomposition as pointer to linked particles.
Definition: Cubo.h:162
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
double ReOverCutOff
Convertion unit R_e over CutOff.
Definition: VarData.h:335
int NBlock
Number of blocks.
Definition: VarData.h:359
double Size[3]
dimension xyz/rad height
Definition: VarData.h:284
void Set(double Val, int Col)
Set the N column.
void SetCoeff()
Set the virial coefficients from the known values of density coex...
Definition: VarData.cpp:171
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int NHeight
Number of monomers per side.
Definition: VarData.h:477