ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/bilayerSys.cpp
(Generate patch)

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 502 by mmeineke, Tue Apr 15 21:47:54 2003 UTC vs.
Revision 536 by mmeineke, Fri May 16 14:28:27 2003 UTC

# Line 1 | Line 1
1 + #include <iostream>
2 +
3   #include <cstdlib>
4   #include <cstring>
5   #include <cmath>
# Line 11 | Line 13
13   #include "bilayerSys.hpp"
14  
15  
16 +
17 + void map( double &x, double &y, double &z,
18 +          double boxX, double boxY, double boxZ );
19 +
20   int buildRandomBilayer( void );
21  
22   void getRandomRot( double rot[3][3] );
# Line 30 | Line 36 | int buildBilayer( int isRandom ){
36   }
37  
38  
33 –
39   int buildRandomBilayer( void ){
40  
41 +  typedef struct{
42 +    double rot[3][3];
43 +    double pos[3];
44 +  } coord;
45 +
46 +
47 +  const double waterRho = 0.0334; // number density per cubic angstrom
48 +  const double waterVol = 4.0 / water_rho; // volume occupied by 4 waters
49 +  const double waterCell = 4.929; // fcc unit cell length
50 +
51 +  const double water_padding = 2.5;
52 +  const double lipid_spaceing = 2.5;
53 +
54 +
55    int i,j,k;
56    int nAtoms, atomIndex, molIndex, molID;
57    int* molSeq;
58    int* molMap;
59 +  int* molStart;
60    int* cardDeck;
61    int deckSize;
62    int rSite, rCard;
63    double cell;
64    int nCells, nSites, siteIndex;
65 +  
66 +  coord *siteArray;
67 +  coord testSite;
68 +
69 +  MoleculeStamp* lipidStamp;
70 +  MoleculeStamp* waterStamp;  
71 +  MoLocator *lipidLocate;
72 +  MoLocator *waterLocate
73 +  int foundLipid, foundWater;
74 +  int nLipids, lipiNatoms, nWaters;
75 +  double testBox, maxLength;
76 +  
77 +  srand48( RAND_SEED );
78 +
79 +
80 +  // set the the lipidStamp
81 +
82 +  foundLipid = 0;
83 +  for(i=0; i<bsInfo.nComponents; i++){
84 +    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
85 +      
86 +      foundlipid = 1;
87 +      lipidStamp = bsInfo.compStamps[i];
88 +      nLipids = bsInfo.componentsNmol[i];
89 +    }
90 +    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
91 +      
92 +      foundWater = 1;
93 +      waterStamp = bsInfo.compStamps[i];
94 +      nWaters = bsInfo.componentsNmol[i];
95 +    }
96 +  }
97 +  if( !foundLipid ){
98 +    sprintf(painCave.errMsg,
99 +            "Could not find lipid \"%s\" in the bass file.\n",
100 +            bsInfo.lipidName );
101 +    painCave.isFatal = 1;
102 +    simError();
103 +  }
104 +  if( !foundWater ){
105 +    sprintf(painCave.errMsg,
106 +            "Could not find water \"%s\" in the bass file.\n",
107 +            bsInfo.waterName );
108 +    painCave.isFatal = 1;
109 +    simError();
110 +  }
111 +  
112 +  //create the temp Molocator and atom Arrays
113 +  
114 +  lipidLocate = new MoLocator( lipidStamp );
115 +  lipidNatoms = lipidStamp->getNAtoms();
116 +  maxLength = lipidLocate->getMaxLength();
117 +
118 +  waterLocate = new MoLocator( waterStamp );
119 +  
120 +  nAtoms = nLipids * lipidNatoms;
121 +
122 +  Atom::createArrays( nAtoms );
123 +  atoms = new Atom*[nAtoms];
124 +
125 +  // create the test box for initial water displacement
126 +
127 +  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
128 +  int nCells = (int)( testBox / waterCell + 1.0 );
129 +  int testWaters = 4 * nCells * nCells * nCells;
130 +  
131 +  double* waterX = new double[testWaters];
132 +  double* waterX = new double[testWaters];
133 +  double* waterX = new double[testWaters];
134 +  
135 +  double x0 = 0.0 - ( testBox * 0.5 );
136 +  double y0 = 0.0 - ( testBox * 0.5 );
137 +  double z0 = 0.0 - ( testBox * 0.5 );
138 +
139 +
140 +  // create an fcc lattice in the water box.
141 +
142 +  int ndx = 0;
143 +  for( i=0; i < nCells; i++ ){
144 +    for( j=0; j < nCells; j++ ){
145 +      for( k=0; k < nCells; k++ ){
146 +        
147 +        waterX[ndx] = i * waterCell + x0;
148 +        waterY[ndx] = j * waterCell + y0;
149 +        waterZ[ndx] = k * waterCell + z0;
150 +        ndx++;
151 +        
152 +        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
153 +        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
154 +        waterZ[ndx] = k * waterCell + z0;
155 +        ndx++;
156 +        
157 +        waterX[ndx] = i * waterCell + x0;
158 +        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
159 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
160 +        ndx++;
161 +        
162 +        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
163 +        waterY[ndx] = j * waterCell + y0;
164 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
165 +        ndx++;
166 +      }
167 +    }
168 +  }
169 +
170 +  // calculate the number of water's displaced by our lipid.
171 +
172 +  testSite.rot[0][0] = 1.0;
173 +  testSite.rot[0][1] = 0.0;
174 +  testSite.rot[0][2] = 0.0;
175 +
176 +  testSite.rot[1][0] = 0.0;
177 +  testSite.rot[1][1] = 1.0;
178 +  testSite.rot[1][2] = 0.0;
179 +
180 +  testSite.rot[2][0] = 0.0;
181 +  testSite.rot[2][1] = 0.0;
182 +  testSite.rot[2][2] = 1.0;
183 +  
184 +  testSite.pos[0] = 0.0;
185 +  testSite.pos[1] = 0.0;
186 +  testSite.pos[2] = 0.0;
187 +
188 +  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
189 +
190 +  int *isActive = new int[testWaters];
191 +  for(i=0; i<testWaters; i++) isActive[i] = 1;
192 +  
193 +  int n_deleted = 0;
194 +  double dx, dy, dz;
195 +  double dx2, dy2, dz2, dSqr;
196 +  double rCutSqr = water_padding * water_padding;
197 +  
198 +  for(i=0; ( (i<testWaters) && isActive[i] ); i++){
199 +    for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
200 +
201 +      dx = waterX[i] - atoms[j]->getX();
202 +      dy = waterY[i] - atoms[j]->getY();
203 +      dz = waterZ[i] - atoms[j]->getZ();
204 +
205 +      map( dx, dy, dz, testBox, testBox, testBox );
206 +
207 +      dx2 = dx * dx;
208 +      dy2 = dy * dy;
209 +      dz2 = dz * dz;
210 +      
211 +      dSqr = dx2 + dy2 + dz2;
212 +      if( dSqr < rCutSqr ){
213 +        isActive[i] = 0;
214 +        n_deleted++;
215 +      }
216 +    }
217 +  }
218 +  
219 +  int targetWaters = nWaters + n_deleted * nLipids;
220 +
221 +  // find the best box size for the sim
222 +
223 +  int testTot;
224 +  int done = 0;
225 +  ndx = 0;
226 +  while( !done ){
227 +    
228 +    ndx++;
229 +    testTot = 4 * ndx * ndx * ndx;
230 +    
231 +    if( testTot >= targetWaters ) done = 1;
232 +  }
233 +
234 +  nCells = ndx;
235 +    
236 +    
237 +  // create the new water box to the new specifications
238 +  
239 +  int newWaters = nCells * nCells * nCells * 4;
240 +  
241 +  delete[] waterX;
242 +  delete[] waterY;
243 +  delete[] waterZ;
244 +  
245 +  waterX = new double[newWater];
246 +  waterY = new double[newWater];
247 +  waterZ = new double[newWater];
248 +  
249 +  double box_x = waterCell * nCells;
250 +  double box_y = waterCell * nCells;
251 +  double box_z = waterCell * nCells;
252 +    
253 +  // create an fcc lattice in the water box.
254 +  
255 +  ndx = 0;
256 +  for( i=0; i < nCells; i++ ){
257 +    for( j=0; j < nCells; j++ ){
258 +      for( k=0; k < nCells; k++ ){
259 +        
260 +        waterX[ndx] = i * waterCell;
261 +        waterY[ndx] = j * waterCell;
262 +        waterZ[ndx] = k * waterCell;
263 +        ndx++;
264 +        
265 +        waterX[ndx] = i * waterCell + 0.5 * waterCell;
266 +        waterY[ndx] = j * waterCell + 0.5 * waterCell;
267 +        waterZ[ndx] = k * waterCell;
268 +        ndx++;
269 +        
270 +        waterX[ndx] = i * waterCell;
271 +        waterY[ndx] = j * waterCell + 0.5 * waterCell;
272 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
273 +        ndx++;
274 +        
275 +        waterX[ndx] = i * waterCell + 0.5 * waterCell;
276 +        waterY[ndx] = j * waterCell;
277 +        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
278 +        ndx++;
279 +      }
280 +    }
281 +  }  
282 +
283 +
284 +
285 +
286 +
287 +
288 +
289 +
290 +
291 +
292 +
293 +
294 +
295 +
296 +  // create the real MoLocator and Atom arrays
297 +  
298 +  nAtoms = 0;
299 +  molIndex = 0;
300 +  locate = new MoLocator*[bsInfo.nComponents];
301 +  molSeq = new int[bsInfo.totNmol];
302 +  molStart = new int[bsInfo.totNmol];  
303 +  for(i=0; i<bsInfo.nComponents; i++){
304 +    locate[i] = new MoLocator( bsInfo.compStamps[i] );
305 +    for(j=0; j<bsInfo.componentsNmol[i]; j++){
306 +      molSeq[molIndex] = i;
307 +      molStart[molIndex] = nAtoms;
308 +      molIndex++;
309 +      nAtoms += bsInfo.compStamps[i]->getNAtoms();
310 +    }
311 +  }
312 +  
313 +  Atom::createArrays( nAtoms );
314 +  atoms = new Atom*[nAtoms];
315 +
316 +
317 +  // find the width, height, and length of the molecule
318 +
319 +  
320 +
321 +
322 + }
323 +
324 +
325 +
326 + int Old_buildRandomBilayer( void ){
327 +
328 +  int i,j,k;
329 +  int nAtoms, atomIndex, molIndex, molID;
330 +  int* molSeq;
331 +  int* molMap;
332 +  int* molStart;
333 +  int* cardDeck;
334 +  int deckSize;
335 +  int rSite, rCard;
336 +  double cell;
337 +  int nCells, nSites, siteIndex;
338    double rot[3][3];
339    double pos[3];
340  
# Line 54 | Line 347 | int buildRandomBilayer( void ){
347    
348    srand48( RAND_SEED );
349    molSeq = NULL;
350 +  molStart = NULL;
351    molMap = NULL;
352    cardDeck = NULL;
353    atoms = NULL;
# Line 102 | Line 396 | int buildRandomBilayer( void ){
396    molIndex = 0;
397    locate = new MoLocator*[bsInfo.nComponents];
398    molSeq = new int[bsInfo.totNmol];
399 +  molStart = new int[bsInfo.totNmol];  
400    for(i=0; i<bsInfo.nComponents; i++){
401      locate[i] = new MoLocator( bsInfo.compStamps[i] );
402      for(j=0; j<bsInfo.componentsNmol[i]; j++){
403        molSeq[molIndex] = i;
404 +      molStart[molIndex] = nAtoms;
405        molIndex++;
406        nAtoms += bsInfo.compStamps[i]->getNAtoms();
407      }
# Line 121 | Line 417 | int buildRandomBilayer( void ){
417    for(i=0; i<bsInfo.nComponents; i++){
418      if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
419    }
420 +  cell *= 1.2; // add a little buffer
421 +
422    cell *= M_SQRT2;
423  
424    siteIndex = 0;
127 –  atomIndex = 0;
425    for(i=0; i<nCells; i++){
426      for(j=0; j<nCells; j++){
427        for(k=0; k<nCells; k++){
# Line 136 | Line 433 | int buildRandomBilayer( void ){
433            
434            getRandomRot( rot );
435            molID = molSeq[molMap[siteIndex]];
436 <          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
437 <          
141 <          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
436 >          atomIndex = molStart[ molMap[siteIndex] ];
437 >          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
438          }
439          siteIndex++;
440  
# Line 146 | Line 442 | int buildRandomBilayer( void ){
442            pos[0] = i * cell + (0.5 * cell);
443            pos[1] = j * cell;
444            pos[2] = k * cell + (0.5 * cell);
445 <          
445 >
446            getRandomRot( rot );
447            molID = molSeq[molMap[siteIndex]];
448 <          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
449 <          
154 <          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
448 >          atomIndex = molStart[ molMap[siteIndex] ];
449 >          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
450          }
451          siteIndex++;
452  
# Line 162 | Line 457 | int buildRandomBilayer( void ){
457            
458            getRandomRot( rot );
459            molID = molSeq[molMap[siteIndex]];
460 <          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
461 <          
167 <          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
460 >          atomIndex = molStart[ molMap[siteIndex] ];
461 >          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
462          }
463          siteIndex++;
464  
# Line 172 | Line 466 | int buildRandomBilayer( void ){
466            pos[0] = i * cell;
467            pos[1] = j * cell + (0.5 * cell);
468            pos[2] = k * cell + (0.5 * cell);
469 <          
469 >
470            getRandomRot( rot );
471            molID = molSeq[molMap[siteIndex]];
472 <          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
473 <          
180 <          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
472 >          atomIndex = molStart[ molMap[siteIndex] ];
473 >          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
474          }
475          siteIndex++;
476        }
# Line 196 | Line 489 | int buildRandomBilayer( void ){
489    simnfo->box_y = bsInfo.boxY;
490    simnfo->box_z = bsInfo.boxZ;
491    
492 <  sprintf( simnfo->statusName, "%s.dump", bsInfo.outPrefix );
492 >  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
493    sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
494 <  
494 >
495    simnfo->atoms = atoms;
496    
497    // set up the writer and write out
# Line 257 | Line 550 | void getRandomRot( double rot[3][3] ){
550    rot[2][2] = cos(theta);  
551   }
552                          
553 +
554 +
555 + void map( double &x, double &y, double &z,
556 +          double boxX, double boxY, double boxZ ){
557 +  
558 +  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
559 +  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
560 +
561 +  if(y < 0) y -= boxY * (double)( (int)( (y / boxY)  - 0.5 ) );
562 +  else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
563 +  
564 +  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
565 +  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
566 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines