[ViewVC] Diff of: group/trunk/OOPSE/utils/bilayerSys.cpp

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 502 by mmeineke, Tue Apr 15 21:47:54 2003 UTC vs.
Revision 573 by mmeineke, Thu Jul 3 19:41:26 2003 UTC

-+
+#include <iostream>
-+
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include "bilayerSys.hpp"
-+
-+
+void map( double &x, double &y, double &z,
-+
+          double boxX, double boxY, double boxZ );
-+
+int buildRandomBilayer( void );
+void getRandomRot( double rot[3][3] );
-–
+int buildRandomBilayer( void ){
-+
-+
+  typedef struct{
-+
+    double rot[3][3];
-+
+    double pos[3];
-+
+  } coord;
-+
-+
-+
+  const double waterRho = 0.0334; // number density per cubic angstrom
-+
+  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
-+
+  const double waterCell = 4.929; // fcc unit cell length
-+
-+
+  const double water_padding = 6.0;
-+
+  const double lipid_spaceing = 8.0;
-+
-+
-+
+  int i,j,k, l;
-+
+  int nAtoms, atomIndex, molIndex, molID;
-+
+  int* molSeq;
-+
+  int* molMap;
-+
+  int* molStart;
-+
+  int* cardDeck;
-+
+  int deckSize;
-+
+  int rSite, rCard;
-+
+  double cell;
-+
+  int nCells, nSites, siteIndex;
-+
-+
+  coord testSite;
-+
-+
+  Atom** atoms;
-+
+  SimInfo* simnfo;
-+
+  DumpWriter* writer;
-+
-+
+  MoleculeStamp* lipidStamp;
-+
+  MoleculeStamp* waterStamp;
-+
+  MoLocator *lipidLocate;
-+
+  MoLocator *waterLocate;
-+
+  int foundLipid, foundWater;
-+
+  int nLipids, lipidNatoms, nWaters, waterNatoms;
-+
+  double testBox, maxLength;
-+
-+
+  srand48( RAND_SEED );
-+
-+
-+
+  // set the the lipidStamp
-+
-+
+  foundLipid = 0;
-+
+  foundWater = 0;
-+
+  for(i=0; i<bsInfo.nComponents; i++){
-+
+    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
-+
-+
+      foundLipid = 1;
-+
+      lipidStamp = bsInfo.compStamps[i];
-+
+      nLipids = bsInfo.componentsNmol[i];
-+
+    }
-+
+    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
-+
-+
+      foundWater = 1;
-+
-+
+      waterStamp = bsInfo.compStamps[i];
-+
+      nWaters = bsInfo.componentsNmol[i];
-+
+    }
-+
+  }
-+
+  if( !foundLipid ){
-+
+    sprintf(painCave.errMsg,
-+
+            "Could not find lipid \"%s\" in the bass file.\n",
-+
+            bsInfo.lipidName );
-+
+    painCave.isFatal = 1;
-+
+    simError();
-+
+  }
-+
+  if( !foundWater ){
-+
+    sprintf(painCave.errMsg,
-+
+            "Could not find solvent \"%s\" in the bass file.\n",
-+
+            bsInfo.waterName );
-+
+    painCave.isFatal = 1;
-+
+    simError();
-+
+  }
-+
-+
+  //create the temp Molocator and atom Arrays
-+
-+
+  lipidLocate = new MoLocator( lipidStamp );
-+
+  lipidNatoms = lipidStamp->getNAtoms();
-+
+  maxLength = lipidLocate->getMaxLength();
-+
-+
+  waterLocate = new MoLocator( waterStamp );
-+
+  waterNatoms = waterStamp->getNAtoms();
-+
-+
+  nAtoms = nLipids * lipidNatoms;
-+
-+
+  Atom::createArrays( nAtoms );
-+
+  atoms = new Atom*[nAtoms];
-+
-+
+  // create the test box for initial water displacement
-+
-+
+  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
-+
+  nCells = (int)( testBox / waterCell + 1.0 );
-+
+  int testWaters = 4 * nCells * nCells * nCells;
-+
-+
+  double* waterX = new double[testWaters];
-+
+  double* waterY = new double[testWaters];
-+
+  double* waterZ = new double[testWaters];
-+
-+
+  double x0 = 0.0 - ( testBox * 0.5 );
-+
+  double y0 = 0.0 - ( testBox * 0.5 );
-+
+  double z0 = 0.0 - ( testBox * 0.5 );
-+
-+
-+
+  // create an fcc lattice in the water box.
-+
-+
+  int ndx = 0;
-+
+  for( i=0; i < nCells; i++ ){
-+
+    for( j=0; j < nCells; j++ ){
-+
+      for( k=0; k < nCells; k++ ){
-+
-+
+        waterX[ndx] = i * waterCell + x0;
-+
+        waterY[ndx] = j * waterCell + y0;
-+
+        waterZ[ndx] = k * waterCell + z0;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
-+
+        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
-+
+        waterZ[ndx] = k * waterCell + z0;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell + x0;
-+
+        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
-+
+        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
-+
+        waterY[ndx] = j * waterCell + y0;
-+
+        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
-+
+        ndx++;
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  // calculate the number of water's displaced by our lipid.
-+
-+
+  testSite.rot[0][0] = 1.0;
-+
+  testSite.rot[0][1] = 0.0;
-+
+  testSite.rot[0][2] = 0.0;
-+
-+
+  testSite.rot[1][0] = 0.0;
-+
+  testSite.rot[1][1] = 1.0;
-+
+  testSite.rot[1][2] = 0.0;
-+
-+
+  testSite.rot[2][0] = 0.0;
-+
+  testSite.rot[2][1] = 0.0;
-+
+  testSite.rot[2][2] = 1.0;
-+
-+
+  testSite.pos[0] = 0.0;
-+
+  testSite.pos[1] = 0.0;
-+
+  testSite.pos[2] = 0.0;
-+
-+
+  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
-+
-+
+  int *isActive = new int[testWaters];
-+
+  for(i=0; i<testWaters; i++) isActive[i] = 1;
-+
-+
+  int n_deleted = 0;
-+
+  double dx, dy, dz;
-+
+  double dx2, dy2, dz2, dSqr;
-+
+  double rCutSqr = water_padding * water_padding;
-+
-+
+  for(i=0; ( (i<testWaters) && isActive[i] ); i++){
-+
+    for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
-+
-+
+      dx = waterX[i] - atoms[j]->getX();
-+
+      dy = waterY[i] - atoms[j]->getY();
-+
+      dz = waterZ[i] - atoms[j]->getZ();
-+
-+
+      map( dx, dy, dz, testBox, testBox, testBox );
-+
-+
+      dx2 = dx * dx;
-+
+      dy2 = dy * dy;
-+
+      dz2 = dz * dz;
-+
-+
+      dSqr = dx2 + dy2 + dz2;
-+
+      if( dSqr < rCutSqr ){
-+
+        isActive[i] = 0;
-+
+        n_deleted++;
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  int targetWaters = nWaters + n_deleted * nLipids;
-+
-+
+  targetWaters = (int) ( targetWaters * 1.2 );
-+
-+
+  // find the best box size for the sim
-+
-+
+  int testTot;
-+
+  int done = 0;
-+
+  ndx = 0;
-+
+  while( !done ){
-+
-+
+    ndx++;
-+
+    testTot = 4 * ndx * ndx * ndx;
-+
-+
+    if( testTot >= targetWaters ) done = 1;
-+
+  }
-+
-+
+  nCells = ndx;
-+
-+
-+
+  // create the new water box to the new specifications
-+
-+
+  int newWaters = nCells * nCells * nCells * 4;
-+
-+
+  delete[] waterX;
-+
+  delete[] waterY;
-+
+  delete[] waterZ;
-+
-+
+  coord* waterSites = new coord[newWaters];
-+
-+
+  double box_x = waterCell * nCells;
-+
+  double box_y = waterCell * nCells;
-+
+  double box_z = waterCell * nCells;
-+
-+
+  // create an fcc lattice in the water box.
-+
-+
+  ndx = 0;
-+
+  for( i=0; i < nCells; i++ ){
-+
+    for( j=0; j < nCells; j++ ){
-+
+      for( k=0; k < nCells; k++ ){
-+
-+
+        waterSites[ndx].pos[0] = i * waterCell;
-+
+        waterSites[ndx].pos[1] = j * waterCell;
-+
+        waterSites[ndx].pos[2] = k * waterCell;
-+
+        ndx++;
-+
-+
+        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
-+
+        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
-+
+        waterSites[ndx].pos[2] = k * waterCell;
-+
+        ndx++;
-+
-+
+        waterSites[ndx].pos[0] = i * waterCell;
-+
+        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
-+
+        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
-+
+        ndx++;
-+
-+
+        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
-+
+        waterSites[ndx].pos[1] = j * waterCell;
-+
+        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
-+
+        ndx++;
-+
+      }
-+
+    }
-+
+  }
-+
-+
-+
+  // clear up memory from the test box
-+
-+
+  for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
-+
-+
+  coord* lipidSites = new coord[nLipids];
-+
-+
+  // start a 3D RSA for the for the lipid placements
-+
-+
-+
+  int reject;
-+
+  int testDX, acceptedDX;
-+
-+
+  rCutSqr = lipid_spaceing * lipid_spaceing;
-+
-+
+  for(i=0; i<nLipids; i++ ){
-+
+    done = 0;
-+
+    while( !done ){
-+
-+
+      lipidSites[i].pos[0] = drand48() * box_x;
-+
+      lipidSites[i].pos[1] = drand48() * box_y;
-+
+      lipidSites[i].pos[2] = drand48() * box_z;
-+
-+
+      getRandomRot( lipidSites[i].rot );
-+
-+
+      ndx = i * lipidNatoms;
-+
+      lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
-+
+                             ndx );
-+
-+
+      reject = 0;
-+
+      for( j=0; !reject && j<i; j++){
-+
+        for(k=0; !reject && k<lipidNatoms; k++){
-+
-+
+          acceptedDX = j*lipidNatoms + k;
-+
+          for(l=0; !reject && l<lipidNatoms; l++){
-+
-+
+            testDX = ndx + l;
-+
-+
+            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
-+
+            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
-+
+            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
-+
-+
+            map( dx, dy, dz, box_x, box_y, box_z );
-+
-+
+            dx2 = dx * dx;
-+
+            dy2 = dy * dy;
-+
+            dz2 = dz * dz;
-+
-+
+            dSqr = dx2 + dy2 + dz2;
-+
+            if( dSqr < rCutSqr ) reject = 1;
-+
+          }
-+
+        }
-+
+      }
-+
-+
+      if( reject ){
-+
-+
+        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
-+
+      }
-+
+      else{
-+
+        done = 1;
-+
+        std::cout << i << " has been accepted\n";
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  // cut out the waters that overlap with the lipids.
-+
-+
+  delete[] isActive;
-+
+  isActive = new int[newWaters];
-+
+  for(i=0; i<newWaters; i++) isActive[i] = 1;
-+
+  int n_active = newWaters;
-+
+  rCutSqr = water_padding * water_padding;
-+
-+
+  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
-+
+    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
-+
-+
+      dx = waterSites[i].pos[0] - atoms[j]->getX();
-+
+      dy = waterSites[i].pos[1] - atoms[j]->getY();
-+
+      dz = waterSites[i].pos[2] - atoms[j]->getZ();
-+
-+
+      map( dx, dy, dz, box_x, box_y, box_z );
-+
-+
+      dx2 = dx * dx;
-+
+      dy2 = dy * dy;
-+
+      dz2 = dz * dz;
-+
-+
+      dSqr = dx2 + dy2 + dz2;
-+
+      if( dSqr < rCutSqr ){
-+
+        isActive[i] = 0;
-+
+        n_active--;
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  if( n_active < nWaters ){
-+
-+
+    sprintf( painCave.errMsg,
-+
+             "Too many waters were removed, edit code and try again.\n" );
-+
-+
+    painCave.isFatal = 1;
-+
+    simError();
-+
+  }
-+
-+
+  int quickKill;
-+
+  while( n_active > nWaters ){
-+
-+
+    quickKill = (int)(drand48()*newWaters);
-+
-+
+    if( isActive[quickKill] ){
-+
+      isActive[quickKill] = 0;
-+
+      n_active--;
-+
+    }
-+
+  }
-+
-+
+  if( n_active != nWaters ){
-+
-+
+    sprintf( painCave.errMsg,
-+
+             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
-+
+             n_active, nWaters );
-+
+    painCave.isFatal = 1;
-+
+    simError();
-+
+  }
-+
-+
+  // clean up our messes before building the final system.
-+
-+
+  for(i=0; i<nAtoms; i++){
-+
-+
+    delete atoms[i];
-+
+  }
-+
+  Atom::destroyArrays();
-+
-+
-+
+  // create the real Atom arrays
-+
-+
+  nAtoms = 0;
-+
+  molIndex = 0;
-+
+  molStart = new int[nLipids + nWaters];
-+
-+
+  for(j=0; j<nLipids; j++){
-+
+    molStart[molIndex] = nAtoms;
-+
+    molIndex++;
-+
+    nAtoms += lipidNatoms;
-+
+  }
-+
-+
+  for(j=0; j<nWaters; j++){
-+
+    molStart[molIndex] = nAtoms;
-+
+    molIndex++;
-+
+    nAtoms += waterNatoms;
-+
+  }
-+
-+
-+
+  Atom::createArrays( nAtoms );
-+
+  atoms = new Atom*[nAtoms];
-+
-+
-+
+  // initialize lipid positions
-+
-+
+  molIndex = 0;
-+
+  for(i=0; i<nLipids; i++ ){
-+
+    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
-+
+                           molStart[molIndex] );
-+
+    molIndex++;
-+
+  }
-+
-+
+  // initialize the water positions
-+
-+
+  for(i=0; i<newWaters; i++){
-+
-+
+    if( isActive[i] ){
-+
-+
+      getRandomRot( waterSites[i].rot );
-+
+      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
-+
+                             molStart[molIndex] );
-+
+      molIndex++;
-+
+    }
-+
+  }
-+
-+
+  // set up the SimInfo object
-+
-+
+  bsInfo.boxX = box_x;
-+
+  bsInfo.boxY = box_y;
-+
+  bsInfo.boxZ = box_z;
-+
-+
+  double boxVector[3];
-+
+  simnfo = new SimInfo();
-+
+  simnfo->n_atoms = nAtoms;
-+
+  boxVector[0] = bsInfo.boxX;
-+
+  boxVector[1] = bsInfo.boxY;
-+
+  boxVector[2] = bsInfo.boxZ;
-+
+  simnfo->setBox( boxVector );
-+
-+
+  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
-+
+  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
-+
-+
+  simnfo->atoms = atoms;
-+
-+
+  // set up the writer and write out
-+
-+
+  writer = new DumpWriter( simnfo );
-+
+  writer->writeFinal( 0.0 );
-+
-+
+  // clean up the memory
-+
-+
+//   if( molMap != NULL )   delete[] molMap;
-+
+//   if( cardDeck != NULL ) delete[] cardDeck;
-+
+//   if( locate != NULL ){
-+
+//     for(i=0; i<bsInfo.nComponents; i++){
-+
+//       delete locate[i];
-+
+//     }
-+
+//     delete[] locate;
-+
+//   }
-+
+//   if( atoms != NULL ){
-+
+//     for(i=0; i<nAtoms; i++){
-+
+//       delete atoms[i];
-+
+//     }
-+
+//     Atom::destroyArrays();
-+
+//     delete[] atoms;
-+
+//   }
-+
+//   if( molSeq != NULL ) delete[] molSeq;
-+
+//   if( simnfo != NULL ) delete simnfo;
-+
+//   if( writer != NULL ) delete writer;
-+
-+
+  return 1;
-+
+}
-+
-+
-+
-+
+int Old_buildRandomBilayer( void ){
-+
+  int i,j,k;
+  int nAtoms, atomIndex, molIndex, molID;
+  int* molSeq;
+  int* molMap;
-+
+  int* molStart;
+  int* cardDeck;
+  int deckSize;
+  int rSite, rCard;
+  srand48( RAND_SEED );
+  molSeq = NULL;
-+
+  molStart = NULL;
+  molMap = NULL;
+  cardDeck = NULL;
+  atoms = NULL;
+  molIndex = 0;
+  locate = new MoLocator*[bsInfo.nComponents];
+  molSeq = new int[bsInfo.totNmol];
-+
+  molStart = new int[bsInfo.totNmol];
+  for(i=0; i<bsInfo.nComponents; i++){
+    locate[i] = new MoLocator( bsInfo.compStamps[i] );
+    for(j=0; j<bsInfo.componentsNmol[i]; j++){
+      molSeq[molIndex] = i;
-+
+      molStart[molIndex] = nAtoms;
+      molIndex++;
+      nAtoms += bsInfo.compStamps[i]->getNAtoms();
+  for(i=0; i<bsInfo.nComponents; i++){
+    if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
-+
+  cell *= 1.2; // add a little buffer
-+
+  cell *= M_SQRT2;
+  siteIndex = 0;
-–
+  atomIndex = 0;
+  for(i=0; i<nCells; i++){
+    for(j=0; j<nCells; j++){
+      for(k=0; k<nCells; k++){
+          getRandomRot( rot );
+          molID = molSeq[molMap[siteIndex]];
-<
+          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
-<
-<
+          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
->
+          atomIndex = molStart[ molMap[siteIndex] ];
->
+          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
+        siteIndex++;
+          pos[0] = i * cell + (0.5 * cell);
+          pos[1] = j * cell;
+          pos[2] = k * cell + (0.5 * cell);
-<
->
+          getRandomRot( rot );
+          molID = molSeq[molMap[siteIndex]];
-<
+          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
-<
-<
+          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
->
+          atomIndex = molStart[ molMap[siteIndex] ];
->
+          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
+        siteIndex++;
+          getRandomRot( rot );
+          molID = molSeq[molMap[siteIndex]];
-<
+          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
-<
-<
+          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
->
+          atomIndex = molStart[ molMap[siteIndex] ];
->
+          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
+        siteIndex++;
+          pos[0] = i * cell;
+          pos[1] = j * cell + (0.5 * cell);
+          pos[2] = k * cell + (0.5 * cell);
-<
->
+          getRandomRot( rot );
+          molID = molSeq[molMap[siteIndex]];
-<
+          locate[molID]->placeMol( pos, rot, &atoms[atomIndex] );
-<
-<
+          atomIndex += bsInfo.compStamps[molID]->getNAtoms();
->
+          atomIndex = molStart[ molMap[siteIndex] ];
->
+          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
+        siteIndex++;
+  bsInfo.boxY = nCells * cell;
+  bsInfo.boxZ = nCells * cell;
-+
+  double boxVector[3];
+  simnfo = new SimInfo();
+  simnfo->n_atoms = nAtoms;
-<
+  simnfo->box_x = bsInfo.boxX;
-<
+  simnfo->box_y = bsInfo.boxY;
-<
+  simnfo->box_z = bsInfo.boxZ;
-<
-<
+  sprintf( simnfo->statusName, "%s.dump", bsInfo.outPrefix );
->
+  boxVector[0] = bsInfo.boxX;
->
+  boxVector[1] = bsInfo.boxY;
->
+  boxVector[2] = bsInfo.boxZ;
->
+  simnfo->setBox( boxVector );
->
->
+  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
+  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
-<
->
+  simnfo->atoms = atoms;
+  // set up the writer and write out
+  writer = new DumpWriter( simnfo );
-<
+  writer->writeFinal();
->
+  writer->writeFinal(0.0);
+  // clean up the memory
+  rot[2][2] = cos(theta);
-+
-+
-+
+void map( double &x, double &y, double &z,
-+
+          double boxX, double boxY, double boxZ ){
-+
-+
+  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
-+
+  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
-+
-+
+  if(y < 0) y -= boxY * (double)( (int)( (y / boxY)  - 0.5 ) );
-+
+  else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
-+
-+
+  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
-+
+  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
-+
+}

Diff Legend

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

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents): Revision 502 by mmeineke, Tue Apr 15 21:47:54 2003 UTC vs. Revision 573 by mmeineke, Thu Jul 3 19:41:26 2003 UTC

Diff Legend

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 502 by mmeineke, Tue Apr 15 21:47:54 2003 UTC vs.
Revision 573 by mmeineke, Thu Jul 3 19:41:26 2003 UTC