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 536 by mmeineke, Fri May 16 14:28:27 2003 UTC vs.
Revision 538 by mmeineke, Sat May 17 16:57:08 2003 UTC

# Line 45 | Line 45 | int buildRandomBilayer( void ){
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
48 >  const double waterVol = 4.0 / waterRho; // 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;
52 >  const double lipid_spaceing = 5.0;
53  
54  
55 <  int i,j,k;
55 >  int i,j,k, l;
56    int nAtoms, atomIndex, molIndex, molID;
57    int* molSeq;
58    int* molMap;
# Line 63 | Line 63 | int buildRandomBilayer( void ){
63    double cell;
64    int nCells, nSites, siteIndex;
65    
66 –  coord *siteArray;
66    coord testSite;
67  
68 +  Atom** atoms;
69 +  SimInfo* simnfo;
70 +  DumpWriter* writer;
71 +
72    MoleculeStamp* lipidStamp;
73    MoleculeStamp* waterStamp;  
74    MoLocator *lipidLocate;
75 <  MoLocator *waterLocate
75 >  MoLocator *waterLocate;
76    int foundLipid, foundWater;
77 <  int nLipids, lipiNatoms, nWaters;
77 >  int nLipids, lipidNatoms, nWaters, waterNatoms;
78    double testBox, maxLength;
79    
80    srand48( RAND_SEED );
# Line 80 | Line 83 | int buildRandomBilayer( void ){
83    // set the the lipidStamp
84  
85    foundLipid = 0;
86 +  foundWater = 0;
87    for(i=0; i<bsInfo.nComponents; i++){
88      if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
89        
90 <      foundlipid = 1;
90 >      foundLipid = 1;
91        lipidStamp = bsInfo.compStamps[i];
92        nLipids = bsInfo.componentsNmol[i];
93      }
94      if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
95        
96        foundWater = 1;
97 +      
98        waterStamp = bsInfo.compStamps[i];
99        nWaters = bsInfo.componentsNmol[i];
100      }
# Line 103 | Line 108 | int buildRandomBilayer( void ){
108    }
109    if( !foundWater ){
110      sprintf(painCave.errMsg,
111 <            "Could not find water \"%s\" in the bass file.\n",
111 >            "Could not find solvent \"%s\" in the bass file.\n",
112              bsInfo.waterName );
113      painCave.isFatal = 1;
114      simError();
# Line 116 | Line 121 | int buildRandomBilayer( void ){
121    maxLength = lipidLocate->getMaxLength();
122  
123    waterLocate = new MoLocator( waterStamp );
124 +  waterNatoms = waterStamp->getNAtoms();
125    
126    nAtoms = nLipids * lipidNatoms;
127  
# Line 125 | Line 131 | int buildRandomBilayer( void ){
131    // create the test box for initial water displacement
132  
133    testBox = maxLength + waterCell * 4.0; // pad with 4 cells
134 <  int nCells = (int)( testBox / waterCell + 1.0 );
134 >  nCells = (int)( testBox / waterCell + 1.0 );
135    int testWaters = 4 * nCells * nCells * nCells;
136    
137    double* waterX = new double[testWaters];
138 <  double* waterX = new double[testWaters];
139 <  double* waterX = new double[testWaters];
138 >  double* waterY = new double[testWaters];
139 >  double* waterZ = new double[testWaters];
140    
141    double x0 = 0.0 - ( testBox * 0.5 );
142    double y0 = 0.0 - ( testBox * 0.5 );
# Line 242 | Line 248 | int buildRandomBilayer( void ){
248    delete[] waterY;
249    delete[] waterZ;
250    
251 <  waterX = new double[newWater];
246 <  waterY = new double[newWater];
247 <  waterZ = new double[newWater];
251 >  coord* waterSites = new coord[newWaters];
252    
253    double box_x = waterCell * nCells;
254    double box_y = waterCell * nCells;
# Line 257 | Line 261 | int buildRandomBilayer( void ){
261      for( j=0; j < nCells; j++ ){
262        for( k=0; k < nCells; k++ ){
263          
264 <        waterX[ndx] = i * waterCell;
265 <        waterY[ndx] = j * waterCell;
266 <        waterZ[ndx] = k * waterCell;
264 >        waterSites[ndx].pos[0] = i * waterCell;
265 >        waterSites[ndx].pos[1] = j * waterCell;
266 >        waterSites[ndx].pos[2] = k * waterCell;
267          ndx++;
268          
269 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
270 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
271 <        waterZ[ndx] = k * waterCell;
269 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
270 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
271 >        waterSites[ndx].pos[2] = k * waterCell;
272          ndx++;
273          
274 <        waterX[ndx] = i * waterCell;
275 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
276 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
274 >        waterSites[ndx].pos[0] = i * waterCell;
275 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
276 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
277          ndx++;
278          
279 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
280 <        waterY[ndx] = j * waterCell;
281 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
279 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
280 >        waterSites[ndx].pos[1] = j * waterCell;
281 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
282          ndx++;
283        }
284      }
285    }  
286  
287  
288 +  // clear up memory from the test box
289  
290 +  for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
291  
292 +  coord* lipidSites = new coord[nLipids];
293  
294 +  // start a 3D RSA for the for the lipid placements
295 +  
296 +  
297 +  int reject;
298 +  int testDX, acceptedDX;
299  
300 +  rCutSqr = lipid_spaceing * lipid_spaceing;
301  
302 +  for(i=0; i<nLipids; i++ ){
303 +    done = 0;
304 +    while( !done ){
305 +      
306 +      lipidSites[i].pos[0] = drand48() * box_x;
307 +      lipidSites[i].pos[1] = drand48() * box_y;
308 +      lipidSites[i].pos[2] = drand48() * box_z;
309 +      
310 +      getRandomRot( lipidSites[i].rot );
311 +      
312 +      ndx = i * lipidNatoms;
313  
314 +      lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
315 +                             ndx );
316 +      
317 +      reject = 0;
318 +      for( j=0; !reject && j<i; j++){
319 +        for(k=0; !reject && k<lipidNatoms; k++){
320 +          
321 +          acceptedDX = j*lipidNatoms + k;
322 +          for(l=0; !reject && l<lipidNatoms; l++){
323 +            
324 +            testDX = ndx + l;
325  
326 +            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
327 +            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
328 +            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
329 +            
330 +            map( dx, dy, dz, box_x, box_y, box_z );
331 +            
332 +            dx2 = dx * dx;
333 +            dy2 = dy * dy;
334 +            dz2 = dz * dz;
335 +            
336 +            dSqr = dx2 + dy2 + dz2;
337 +            if( dSqr < rCutSqr ) reject = 1;
338 +          }
339 +        }
340 +      }
341  
342 +      if( reject ){
343  
344 +        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
345 +      }
346 +      else{
347 +        done = 1;
348 +        std::cout << i << " has been accepted\n";
349 +      }
350 +    }
351 +  }
352 +        
353 +  // cut out the waters that overlap with the lipids.
354 +  
355 +  delete[] isActive;
356 +  isActive = new int[newWaters];
357 +  for(i=0; i<newWaters; i++) isActive[i] = 1;
358 +  int n_active = newWaters;
359 +  rCutSqr = water_padding * water_padding;
360 +  
361 +  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
362 +    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
363  
364 +      dx = waterSites[i].pos[0] - atoms[j]->getX();
365 +      dy = waterSites[i].pos[1] - atoms[j]->getY();
366 +      dz = waterSites[i].pos[2] - atoms[j]->getZ();
367  
368 +      map( dx, dy, dz, box_x, box_y, box_z );
369  
370 <  // create the real MoLocator and Atom arrays
370 >      dx2 = dx * dx;
371 >      dy2 = dy * dy;
372 >      dz2 = dz * dz;
373 >      
374 >      dSqr = dx2 + dy2 + dz2;
375 >      if( dSqr < rCutSqr ){
376 >        isActive[i] = 0;
377 >        n_active--;
378 >      }
379 >    }
380 >  }
381 >
382 >  if( n_active < nWaters ){
383 >    
384 >    sprintf( painCave.errMsg,
385 >             "Too many waters were removed, edit code and try again.\n" );
386 >    
387 >    painCave.isFatal = 1;
388 >    simError();
389 >  }
390 >
391 >  int quickKill;
392 >  while( n_active > nWaters ){
393 >
394 >    quickKill = (int)(drand48()*newWaters);
395 >
396 >    if( isActive[quickKill] ){
397 >      isActive[quickKill] = 0;
398 >      n_active--;
399 >    }
400 >  }
401 >
402 >  if( n_active != nWaters ){
403 >    
404 >    sprintf( painCave.errMsg,
405 >             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
406 >             n_active, nWaters );
407 >    painCave.isFatal = 1;
408 >    simError();
409 >  }
410 >
411 >  // clean up our messes before building the final system.
412 >
413 >  for(i=0; i<nAtoms; i++){
414 >    
415 >    delete atoms[i];
416 >  }
417 >  Atom::destroyArrays();
418    
419 +
420 +  // create the real Atom arrays
421 +  
422    nAtoms = 0;
423    molIndex = 0;
424 <  locate = new MoLocator*[bsInfo.nComponents];
425 <  molSeq = new int[bsInfo.totNmol];
426 <  molStart = new int[bsInfo.totNmol];  
427 <  for(i=0; i<bsInfo.nComponents; i++){
428 <    locate[i] = new MoLocator( bsInfo.compStamps[i] );
429 <    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 <    }
424 >  molStart = new int[nLipids + nWaters];  
425 >  
426 >  for(j=0; j<nLipids; j++){
427 >    molStart[molIndex] = nAtoms;
428 >    molIndex++;
429 >    nAtoms += lipidNatoms;
430    }
431 +
432 +  for(j=0; j<nWaters; j++){
433 +    molStart[molIndex] = nAtoms;
434 +    molIndex++;
435 +    nAtoms += waterNatoms;
436 +  }
437    
438 +  
439    Atom::createArrays( nAtoms );
440    atoms = new Atom*[nAtoms];
441  
442 <
443 <  // find the width, height, and length of the molecule
442 >  
443 >  // initialize lipid positions
444  
445 +  molIndex = 0;
446 +  for(i=0; i<nLipids; i++ ){
447 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
448 +                           molStart[molIndex] );
449 +    molIndex++;
450 +  }
451 +
452 +  // initialize the water positions
453 +
454 +  for(i=0; i<newWaters; i++){
455 +    
456 +    if( isActive[i] ){
457 +      
458 +      getRandomRot( waterSites[i].rot );
459 +      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
460 +                             molStart[molIndex] );
461 +      molIndex++;
462 +    }
463 +  }  
464 +
465 +  // set up the SimInfo object
466 +  
467 +  bsInfo.boxX = box_x;
468 +  bsInfo.boxY = box_y;
469 +  bsInfo.boxZ = box_z;
470 +  
471 +  simnfo = new SimInfo();
472 +  simnfo->n_atoms = nAtoms;
473 +  simnfo->box_x = bsInfo.boxX;
474 +  simnfo->box_y = bsInfo.boxY;
475 +  simnfo->box_z = bsInfo.boxZ;
476    
477 +  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
478 +  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
479  
480 +  simnfo->atoms = atoms;
481 +  
482 +  // set up the writer and write out
483 +  
484 +  writer = new DumpWriter( simnfo );
485 +  writer->writeFinal();
486 +    
487 +  // clean up the memory
488 +  
489 + //   if( molMap != NULL )   delete[] molMap;
490 + //   if( cardDeck != NULL ) delete[] cardDeck;
491 + //   if( locate != NULL ){
492 + //     for(i=0; i<bsInfo.nComponents; i++){
493 + //       delete locate[i];
494 + //     }
495 + //     delete[] locate;
496 + //   }
497 + //   if( atoms != NULL ){
498 + //     for(i=0; i<nAtoms; i++){
499 + //       delete atoms[i];
500 + //     }
501 + //     Atom::destroyArrays();
502 + //     delete[] atoms;
503 + //   }
504 + //   if( molSeq != NULL ) delete[] molSeq;
505 + //   if( simnfo != NULL ) delete simnfo;
506 + //   if( writer != NULL ) delete writer;
507  
508 +  return 1;
509   }
510  
511 +  
512  
325 –
513   int Old_buildRandomBilayer( void ){
514  
515    int i,j,k;

Diff Legend

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