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

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 535 by mmeineke, Fri Apr 25 16:02:26 2003 UTC vs.
Revision 536 by mmeineke, Fri May 16 14:28:27 2003 UTC

+#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 / water_rho; // volume occupied by 4 waters
-+
+  const double waterCell = 4.929; // fcc unit cell length
-+
-+
+  const double water_padding = 2.5;
-+
+  const double lipid_spaceing = 2.5;
-+
-+
-+
+  int i,j,k;
-+
+  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 *siteArray;
-+
+  coord testSite;
-+
-+
+  MoleculeStamp* lipidStamp;
-+
+  MoleculeStamp* waterStamp;
-+
+  MoLocator *lipidLocate;
-+
+  MoLocator *waterLocate
-+
+  int foundLipid, foundWater;
-+
+  int nLipids, lipiNatoms, nWaters;
-+
+  double testBox, maxLength;
-+
-+
+  srand48( RAND_SEED );
-+
-+
-+
+  // set the the lipidStamp
-+
-+
+  foundLipid = 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 water \"%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 );
-+
-+
+  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
-+
+  int nCells = (int)( testBox / waterCell + 1.0 );
-+
+  int testWaters = 4 * nCells * nCells * nCells;
-+
-+
+  double* waterX = new double[testWaters];
-+
+  double* waterX = new double[testWaters];
-+
+  double* waterX = 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;
-+
-+
+  // 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;
-+
-+
+  waterX = new double[newWater];
-+
+  waterY = new double[newWater];
-+
+  waterZ = new double[newWater];
-+
-+
+  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++ ){
-+
-+
+        waterX[ndx] = i * waterCell;
-+
+        waterY[ndx] = j * waterCell;
-+
+        waterZ[ndx] = k * waterCell;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell + 0.5 * waterCell;
-+
+        waterY[ndx] = j * waterCell + 0.5 * waterCell;
-+
+        waterZ[ndx] = k * waterCell;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell;
-+
+        waterY[ndx] = j * waterCell + 0.5 * waterCell;
-+
+        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
-+
+        ndx++;
-+
-+
+        waterX[ndx] = i * waterCell + 0.5 * waterCell;
-+
+        waterY[ndx] = j * waterCell;
-+
+        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
-+
+        ndx++;
-+
+      }
-+
+    }
-+
+  }
-+
-+
-+
-+
-+
-+
-+
-+
-+
-+
-+
-+
-+
-+
+  // create the real MoLocator and Atom arrays
-+
-+
+  nAtoms = 0;
-+
+  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();
-+
+    }
-+
+  }
-+
-+
+  Atom::createArrays( nAtoms );
-+
+  atoms = new Atom*[nAtoms];
-+
-+
-+
+  // find the width, height, and length of the molecule
-+
-+
-+
-+
-+
+}
-+
-+
-+
-+
+int Old_buildRandomBilayer( void ){
-+
+  int i,j,k;
+  int nAtoms, atomIndex, molIndex, molID;
+  int* molSeq;
+  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 535 by mmeineke, Fri Apr 25 16:02:26 2003 UTC vs. Revision 536 by mmeineke, Fri May 16 14:28:27 2003 UTC

Diff Legend

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 535 by mmeineke, Fri Apr 25 16:02:26 2003 UTC vs.
Revision 536 by mmeineke, Fri May 16 14:28:27 2003 UTC