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 543 by mmeineke, Fri May 30 21:32:33 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;
51 >  const double water_padding = 6.0;
52 >  const double lipid_spaceing = 8.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 218 | Line 224 | int buildRandomBilayer( void ){
224    
225    int targetWaters = nWaters + n_deleted * nLipids;
226  
227 +  targetWaters = (int) ( targetWaters * 1.2 );
228 +
229    // find the best box size for the sim
230  
231    int testTot;
# Line 242 | Line 250 | int buildRandomBilayer( void ){
250    delete[] waterY;
251    delete[] waterZ;
252    
253 <  waterX = new double[newWater];
246 <  waterY = new double[newWater];
247 <  waterZ = new double[newWater];
253 >  coord* waterSites = new coord[newWaters];
254    
255    double box_x = waterCell * nCells;
256    double box_y = waterCell * nCells;
# Line 257 | Line 263 | int buildRandomBilayer( void ){
263      for( j=0; j < nCells; j++ ){
264        for( k=0; k < nCells; k++ ){
265          
266 <        waterX[ndx] = i * waterCell;
267 <        waterY[ndx] = j * waterCell;
268 <        waterZ[ndx] = k * waterCell;
266 >        waterSites[ndx].pos[0] = i * waterCell;
267 >        waterSites[ndx].pos[1] = j * waterCell;
268 >        waterSites[ndx].pos[2] = k * waterCell;
269          ndx++;
270          
271 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
272 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
273 <        waterZ[ndx] = k * waterCell;
271 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
272 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
273 >        waterSites[ndx].pos[2] = k * waterCell;
274          ndx++;
275          
276 <        waterX[ndx] = i * waterCell;
277 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
278 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
276 >        waterSites[ndx].pos[0] = i * waterCell;
277 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
278 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
279          ndx++;
280          
281 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
282 <        waterY[ndx] = j * waterCell;
283 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
281 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
282 >        waterSites[ndx].pos[1] = j * waterCell;
283 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
284          ndx++;
285        }
286      }
287    }  
288  
289  
290 +  // clear up memory from the test box
291  
292 +  for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
293  
294 +  coord* lipidSites = new coord[nLipids];
295  
296 +  // start a 3D RSA for the for the lipid placements
297 +  
298 +  
299 +  int reject;
300 +  int testDX, acceptedDX;
301  
302 +  rCutSqr = lipid_spaceing * lipid_spaceing;
303  
304 +  for(i=0; i<nLipids; i++ ){
305 +    done = 0;
306 +    while( !done ){
307 +      
308 +      lipidSites[i].pos[0] = drand48() * box_x;
309 +      lipidSites[i].pos[1] = drand48() * box_y;
310 +      lipidSites[i].pos[2] = drand48() * box_z;
311 +      
312 +      getRandomRot( lipidSites[i].rot );
313 +      
314 +      ndx = i * lipidNatoms;
315  
316 +      lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
317 +                             ndx );
318 +      
319 +      reject = 0;
320 +      for( j=0; !reject && j<i; j++){
321 +        for(k=0; !reject && k<lipidNatoms; k++){
322 +          
323 +          acceptedDX = j*lipidNatoms + k;
324 +          for(l=0; !reject && l<lipidNatoms; l++){
325 +            
326 +            testDX = ndx + l;
327  
328 +            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
329 +            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
330 +            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
331 +            
332 +            map( dx, dy, dz, box_x, box_y, box_z );
333 +            
334 +            dx2 = dx * dx;
335 +            dy2 = dy * dy;
336 +            dz2 = dz * dz;
337 +            
338 +            dSqr = dx2 + dy2 + dz2;
339 +            if( dSqr < rCutSqr ) reject = 1;
340 +          }
341 +        }
342 +      }
343  
344 +      if( reject ){
345  
346 +        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
347 +      }
348 +      else{
349 +        done = 1;
350 +        std::cout << i << " has been accepted\n";
351 +      }
352 +    }
353 +  }
354 +        
355 +  // cut out the waters that overlap with the lipids.
356 +  
357 +  delete[] isActive;
358 +  isActive = new int[newWaters];
359 +  for(i=0; i<newWaters; i++) isActive[i] = 1;
360 +  int n_active = newWaters;
361 +  rCutSqr = water_padding * water_padding;
362 +  
363 +  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
364 +    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
365  
366 +      dx = waterSites[i].pos[0] - atoms[j]->getX();
367 +      dy = waterSites[i].pos[1] - atoms[j]->getY();
368 +      dz = waterSites[i].pos[2] - atoms[j]->getZ();
369  
370 +      map( dx, dy, dz, box_x, box_y, box_z );
371  
372 <  // create the real MoLocator and Atom arrays
372 >      dx2 = dx * dx;
373 >      dy2 = dy * dy;
374 >      dz2 = dz * dz;
375 >      
376 >      dSqr = dx2 + dy2 + dz2;
377 >      if( dSqr < rCutSqr ){
378 >        isActive[i] = 0;
379 >        n_active--;
380 >      }
381 >    }
382 >  }
383 >
384 >  if( n_active < nWaters ){
385 >    
386 >    sprintf( painCave.errMsg,
387 >             "Too many waters were removed, edit code and try again.\n" );
388 >    
389 >    painCave.isFatal = 1;
390 >    simError();
391 >  }
392 >
393 >  int quickKill;
394 >  while( n_active > nWaters ){
395 >
396 >    quickKill = (int)(drand48()*newWaters);
397 >
398 >    if( isActive[quickKill] ){
399 >      isActive[quickKill] = 0;
400 >      n_active--;
401 >    }
402 >  }
403 >
404 >  if( n_active != nWaters ){
405 >    
406 >    sprintf( painCave.errMsg,
407 >             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
408 >             n_active, nWaters );
409 >    painCave.isFatal = 1;
410 >    simError();
411 >  }
412 >
413 >  // clean up our messes before building the final system.
414 >
415 >  for(i=0; i<nAtoms; i++){
416 >    
417 >    delete atoms[i];
418 >  }
419 >  Atom::destroyArrays();
420    
421 +
422 +  // create the real Atom arrays
423 +  
424    nAtoms = 0;
425    molIndex = 0;
426 <  locate = new MoLocator*[bsInfo.nComponents];
427 <  molSeq = new int[bsInfo.totNmol];
428 <  molStart = new int[bsInfo.totNmol];  
429 <  for(i=0; i<bsInfo.nComponents; i++){
430 <    locate[i] = new MoLocator( bsInfo.compStamps[i] );
431 <    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 <    }
426 >  molStart = new int[nLipids + nWaters];  
427 >  
428 >  for(j=0; j<nLipids; j++){
429 >    molStart[molIndex] = nAtoms;
430 >    molIndex++;
431 >    nAtoms += lipidNatoms;
432    }
433 +
434 +  for(j=0; j<nWaters; j++){
435 +    molStart[molIndex] = nAtoms;
436 +    molIndex++;
437 +    nAtoms += waterNatoms;
438 +  }
439    
440 +  
441    Atom::createArrays( nAtoms );
442    atoms = new Atom*[nAtoms];
443  
444 <
445 <  // find the width, height, and length of the molecule
444 >  
445 >  // initialize lipid positions
446  
447 +  molIndex = 0;
448 +  for(i=0; i<nLipids; i++ ){
449 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
450 +                           molStart[molIndex] );
451 +    molIndex++;
452 +  }
453 +
454 +  // initialize the water positions
455 +
456 +  for(i=0; i<newWaters; i++){
457 +    
458 +    if( isActive[i] ){
459 +      
460 +      getRandomRot( waterSites[i].rot );
461 +      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
462 +                             molStart[molIndex] );
463 +      molIndex++;
464 +    }
465 +  }  
466 +
467 +  // set up the SimInfo object
468 +  
469 +  bsInfo.boxX = box_x;
470 +  bsInfo.boxY = box_y;
471 +  bsInfo.boxZ = box_z;
472 +  
473 +  simnfo = new SimInfo();
474 +  simnfo->n_atoms = nAtoms;
475 +  simnfo->box_x = bsInfo.boxX;
476 +  simnfo->box_y = bsInfo.boxY;
477 +  simnfo->box_z = bsInfo.boxZ;
478    
479 +  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
480 +  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
481  
482 +  simnfo->atoms = atoms;
483 +  
484 +  // set up the writer and write out
485 +  
486 +  writer = new DumpWriter( simnfo );
487 +  writer->writeFinal();
488 +    
489 +  // clean up the memory
490 +  
491 + //   if( molMap != NULL )   delete[] molMap;
492 + //   if( cardDeck != NULL ) delete[] cardDeck;
493 + //   if( locate != NULL ){
494 + //     for(i=0; i<bsInfo.nComponents; i++){
495 + //       delete locate[i];
496 + //     }
497 + //     delete[] locate;
498 + //   }
499 + //   if( atoms != NULL ){
500 + //     for(i=0; i<nAtoms; i++){
501 + //       delete atoms[i];
502 + //     }
503 + //     Atom::destroyArrays();
504 + //     delete[] atoms;
505 + //   }
506 + //   if( molSeq != NULL ) delete[] molSeq;
507 + //   if( simnfo != NULL ) delete simnfo;
508 + //   if( writer != NULL ) delete writer;
509  
510 +  return 1;
511   }
512  
513 +  
514  
325 –
515   int Old_buildRandomBilayer( void ){
516  
517    int i,j,k;

Diff Legend

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