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 537 by mmeineke, Fri May 16 21:37:56 2003 UTC

# Line 49 | Line 49 | int buildRandomBilayer( void ){
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;
# Line 63 | Line 63 | int buildRandomBilayer( void ){
63    double cell;
64    int nCells, nSites, siteIndex;
65    
66 –  coord *siteArray;
66    coord testSite;
67  
68    MoleculeStamp* lipidStamp;
# Line 71 | Line 70 | int buildRandomBilayer( void ){
70    MoLocator *lipidLocate;
71    MoLocator *waterLocate
72    int foundLipid, foundWater;
73 <  int nLipids, lipiNatoms, nWaters;
73 >  int nLipids, lipiNatoms, nWaters, waterNatoms;
74    double testBox, maxLength;
75    
76    srand48( RAND_SEED );
# Line 103 | Line 102 | int buildRandomBilayer( void ){
102    }
103    if( !foundWater ){
104      sprintf(painCave.errMsg,
105 <            "Could not find water \"%s\" in the bass file.\n",
105 >            "Could not find solvent \"%s\" in the bass file.\n",
106              bsInfo.waterName );
107      painCave.isFatal = 1;
108      simError();
# Line 116 | Line 115 | int buildRandomBilayer( void ){
115    maxLength = lipidLocate->getMaxLength();
116  
117    waterLocate = new MoLocator( waterStamp );
118 +  waterNatoms = waterStamp->getNatoms();
119    
120    nAtoms = nLipids * lipidNatoms;
121  
# Line 242 | Line 242 | int buildRandomBilayer( void ){
242    delete[] waterY;
243    delete[] waterZ;
244    
245 <  waterX = new double[newWater];
246 <  waterY = new double[newWater];
247 <  waterZ = new double[newWater];
245 >  coord* waterSites = new coord[newWaters];
246    
247    double box_x = waterCell * nCells;
248    double box_y = waterCell * nCells;
# Line 257 | Line 255 | int buildRandomBilayer( void ){
255      for( j=0; j < nCells; j++ ){
256        for( k=0; k < nCells; k++ ){
257          
258 <        waterX[ndx] = i * waterCell;
259 <        waterY[ndx] = j * waterCell;
260 <        waterZ[ndx] = k * waterCell;
258 >        waterSites[ndx].pos[0] = i * waterCell;
259 >        waterSites[ndx].pos[1] = j * waterCell;
260 >        waterSites[ndx].pos[2] = k * waterCell;
261          ndx++;
262          
263 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
264 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
265 <        waterZ[ndx] = k * waterCell;
263 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
264 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
265 >        waterSites[ndx].pos[2] = k * waterCell;
266          ndx++;
267          
268 <        waterX[ndx] = i * waterCell;
269 <        waterY[ndx] = j * waterCell + 0.5 * waterCell;
270 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
268 >        waterSites[ndx].pos[0] = i * waterCell;
269 >        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
270 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
271          ndx++;
272          
273 <        waterX[ndx] = i * waterCell + 0.5 * waterCell;
274 <        waterY[ndx] = j * waterCell;
275 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
273 >        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
274 >        waterSites[ndx].pos[1] = j * waterCell;
275 >        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
276          ndx++;
277        }
278      }
279    }  
280  
281  
282 +  // clear up memory from the test box
283  
284 +  for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
285  
286 +  coord* lipidSites = new coord[nLipids];
287  
288 +  // start a 3D RSA for the for the lipid placements
289 +  
290 +  
291 +  int reject;
292 +  int testDX, acceptedDX;
293  
294 +  rCutSqr = lipid_spaceing * lipid_spaceing;
295  
296 +  for(i=0; i<nLipids; i++ ){
297 +    done = 0;
298 +    while( !done ){
299 +      
300 +      lipidSite[i].pos[0] = drand48() * box_x;
301 +      lipidSite[i].pos[1] = drand48() * box_y;
302 +      lipidSite[i].pos[2] = drand48() * box_z;
303 +      
304 +      getRandomRot( lipidSite[i].rot );
305 +      
306 +      ndx = i * lipidNatoms;
307 +
308 +      lipidLocate->placeMol( lipidSite[i].pos, lipidSite[i].rot, atoms, ndx );
309 +      
310 +      reject = 0;
311 +      for( j=0; !reject && j<i; j++){
312 +        for(k=0; !reject && k<lipidNatoms; k++){
313 +          
314 +          acceptedDX = j*lipidNatoms + k;
315 +          for(l=0; !reject && l<lipidNatoms; l++){
316 +            
317 +            testDX = ndx + l;
318  
319 +            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
320 +            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
321 +            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
322 +            
323 +            map( dx, dy, dz, box_x, box_y, box_z );
324 +            
325 +            dx2 = dx * dx;
326 +            dy2 = dy * dy;
327 +            dz2 = dz * dz;
328 +            
329 +            dSqr = dx2 + dy2 + dz2;
330 +            if( dSqr < rCutSqr ) reject = 1;
331 +          }
332 +        }
333 +      }
334  
335 +      if( reject ){
336  
337 +        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
338 +      }
339 +      else{
340 +        done = 1;
341 +        std::cout << i << " has been accepted\n";
342 +      }
343 +    }
344 +  }
345 +        
346 +  // cut out the waters that overlap with the lipids.
347 +  
348 +  delete[] isActive;
349 +  isActive = new int[newWaters];
350 +  for(i=0; i<newWaters; i++) isActive[i] = 1;
351 +  int n_active = newWaters;
352 +  rCutSqr = water_padding * water_padding;
353 +  
354 +  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
355 +    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
356  
357 +      dx = waterSite[i].pos[0] - rsaAtoms[j]->getX();
358 +      dy = waterSite[i].pos[1] - rsaAtoms[j]->getY();
359 +      dz = waterSite[i].pos[2] - rsaAtoms[j]->getZ();
360  
361 +      map( dx, dy, dz, box_x, box_y, box_z );
362  
363 +      dx2 = dx * dx;
364 +      dy2 = dy * dy;
365 +      dz2 = dz * dz;
366 +      
367 +      dSqr = dx2 + dy2 + dz2;
368 +      if( dSqr < rCutSqr ){
369 +        isActive[i] = 0;
370 +        n_active--;
371 +      }
372 +    }
373 +  }
374  
375 <  // create the real MoLocator and Atom arrays
375 >  if( n_active < nWaters ){
376 >    
377 >    sprintf( painCave.errMsg,
378 >             "Too many waters were removed, edit code and try again.\n" );
379 >    
380 >    painCave.isFatal = 1;
381 >    simError();
382 >  }
383 >
384 >  int quickKill;
385 >  while( n_active > nWaters ){
386 >
387 >    quickKill = (int)(drand48()*newWaters);
388 >
389 >    if( isActive[quickKill] ){
390 >      isActive[quickKill] = 0;
391 >      n_active--;
392 >    }
393 >  }
394 >
395 >  if( n_active != nWaters ){
396 >    
397 >    sprintf( painCave.errMsg,
398 >             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
399 >             n_active, nWaters );
400 >    painCave.isFatal = 1;
401 >    simError();
402 >  }
403 >
404 >  // clean up our messes before building the final system.
405 >
406 >  for(i=0; i<nAtoms; i++){
407 >    
408 >    delete atoms[i];
409 >  }
410 >  Atom::destroyArrays();
411    
412 +
413 +  // create the real Atom arrays
414 +  
415    nAtoms = 0;
416    molIndex = 0;
417 <  locate = new MoLocator*[bsInfo.nComponents];
418 <  molSeq = new int[bsInfo.totNmol];
419 <  molStart = new int[bsInfo.totNmol];  
420 <  for(i=0; i<bsInfo.nComponents; i++){
421 <    locate[i] = new MoLocator( bsInfo.compStamps[i] );
422 <    for(j=0; j<bsInfo.componentsNmol[i]; j++){
423 <      molSeq[molIndex] = i;
424 <      molStart[molIndex] = nAtoms;
425 <      molIndex++;
426 <      nAtoms += bsInfo.compStamps[i]->getNAtoms();
427 <    }
417 >  locate = new MoLocator*[2];
418 >  molSeq = new int[nLipids + nWaters];
419 >  molStart = new int[nLipids + nWaters];  
420 >  
421 >  locate[0] = lipidLocate;
422 >  for(j=0; j<nLipids; j++){
423 >    molSeq[molIndex] = 0;
424 >    molStart[molIndex] = nAtoms;
425 >    molIndex++;
426 >    nAtoms += lipidNatoms;
427 >  }
428 >
429 >  locate[1] = waterLocate;
430 >  for(j=0; j<nLipids; j++){
431 >    molSeq[molIndex] = 1;
432 >    molStart[molIndex] = nAtoms;
433 >    molIndex++;
434 >    nAtoms += waterNatoms;
435    }
436    
437 +  
438    Atom::createArrays( nAtoms );
439    atoms = new Atom*[nAtoms];
440  
441 +  
442 +  // initialize lipid positions
443  
444 <  // find the width, height, and length of the molecule
444 >  
445 >  
446 >  
447 >
448 >  // set up the SimInfo object
449 >  
450 >  bsInfo.boxX = box_x;
451 >  bsInfo.boxY = box_y;
452 >  bsInfo.boxZ = box_z;
453 >  
454 >  simnfo = new SimInfo();
455 >  simnfo->n_atoms = nAtoms;
456 >  simnfo->box_x = bsInfo.boxX;
457 >  simnfo->box_y = bsInfo.boxY;
458 >  simnfo->box_z = bsInfo.boxZ;
459 >  
460 >  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
461 >  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
462 >
463 >  simnfo->atoms = atoms;
464 >  
465 >  // set up the writer and write out
466 >  
467 >  writer = new DumpWriter( simnfo );
468 >  writer->writeFinal();
469 >    
470 >  // clean up the memory
471 >  
472 >  if( molMap != NULL )   delete[] molMap;
473 >  if( cardDeck != NULL ) delete[] cardDeck;
474 >  if( locate != NULL ){
475 >    for(i=0; i<bsInfo.nComponents; i++){
476 >      delete locate[i];
477 >    }
478 >    delete[] locate;
479 >  }
480 >  if( atoms != NULL ){
481 >    for(i=0; i<nAtoms; i++){
482 >      delete atoms[i];
483 >    }
484 >    Atom::destroyArrays();
485 >    delete[] atoms;
486 >  }
487 >  if( molSeq != NULL ) delete[] molSeq;
488 >  if( simnfo != NULL ) delete simnfo;
489 >  if( writer != NULL ) delete writer;
490  
491 +  return 1;
492 + }
493 +
494    
495  
496  

Diff Legend

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