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

Comparing trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp (file contents):
Revision 700 by chuckv, Mon Aug 18 20:59:47 2003 UTC vs.
Revision 817 by gezelter, Fri Oct 24 17:36:18 2003 UTC

-+
+#include <iostream>
-+
+#include <vector>
-+
+#include <algorithm>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
-+
+#include "simError.h"
+#include "SimInfo.hpp"
+#include "ReadWrite.hpp"
+#include "latticeBuilder.hpp"
-+
+class SortCond{
-+
-+
+public:
-+
+  bool operator()(const pair<int, double>& p1, const pair<int, double>& p2){
-+
+    return p1.second < p2.second;
-+
+  }
-+
-+
-+
+};
-+
-+
+void buildMap( double &x, double &y, double &z,
+               double boxX, double boxY, double boxZ );
+int buildRandomBilayer( void );
-+
+int buildLatticeBilayer( int isHexLattice,
-+
+                         double hexSpacing,
-+
+                         double aLat,
-+
+                         double bLat,
-+
+                         int targetNlipid,
-+
+                         double targetWaterLipidRatio,
-+
+                         double leafSpacing);
+void getRandomRot( double rot[3][3] );
-+
+void getEulerRot( double theta, double phi, double psi, double rot[3][3] );
-+
+void getUnitRot( double unit[3], double rot[3][3] );
+int buildBilayer( int isRandom ){
+    return buildRandomBilayer();
+  else{
-<
+    sprintf( painCave.errMsg,
-<
+             "Cannot currently create a non-random bilayer.\n" );
-<
+    painCave.isFatal = 1;
-<
+    simError();
-<
+    return 0;
->
->
->
->
+    return buildLatticeBilayer();
+  // find the best box size for the sim
-+
+  int nCellsX, nCellsY, nCellsZ;
-+
-+
+  const double boxTargetX = 66.22752;
-+
+  const double boxTargetY = 60.53088;
-+
-+
+  nCellsX = (int)ceil(boxTargetX / waterCell);
-+
+  nCellsY = (int)ceil(boxTargetY / waterCell);
-+
+  int testTot;
+  int done = 0;
-<
+  ndx = 0;
->
+  nCellsZ = 0;
+  while( !done ){
-<
+    ndx++;
-<
+    testTot = 4 * ndx * ndx * ndx;
->
+    nCellsZ++;
->
+    testTot = 4 * nCellsX * nCellsY * nCellsZ;
+    if( testTot >= targetWaters ) done = 1;
-–
+  nCells = ndx;
-–
-–
+  // create the new water box to the new specifications
-<
+  int newWaters = nCells * nCells * nCells * 4;
->
+  int newWaters = nCellsX * nCellsY * nCellsZ * 4;
+  delete[] waterX;
+  delete[] waterY;
+  coord* waterSites = new coord[newWaters];
-<
+  double box_x = waterCell * nCells;
-<
+  double box_y = waterCell * nCells;
-<
+  double box_z = waterCell * nCells;
->
+  double box_x = waterCell * nCellsX;
->
+  double box_y = waterCell * nCellsY;
->
+  double box_z = waterCell * nCellsZ;
+  // 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++ ){
->
+  for( i=0; i < nCellsX; i++ ){
->
+    for( j=0; j < nCellsY; j++ ){
->
+      for( k=0; k < nCellsZ; k++ ){
+        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
+        for(l=0; l<4; l++){
-+
-+
+  // zSort of the lipid positions
-+
-+
-+
+  vector< pair<int,double> >zSortArray;
-+
+  for(i=0;i<nLipids;i++)
-+
+    zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );
-+
-+
+  sort(zSortArray.begin(),zSortArray.end(),SortCond());
-+
-+
+  ofstream outFile( "./zipper.bass", ios::app);
-+
-+
+  for(i=0; i<nLipids; i++){
-+
+    outFile << "zConstraint[" << i << "]{\n"
-+
+            << "  molIndex = " << zSortArray[i].first <<  ";\n"
-+
+            << "  zPos = ";
-+
-+
+    if(i<32) outFile << "60.0;\n";
-+
+    else outFile << "100.0;\n";
-+
-+
+    outFile << "  kRatio = 0.5;\n"
-+
+            << "}\n";
-+
+  }
-+
-+
+  outFile.close();
-+
-+
+  // cut out the waters that overlap with the lipids.
+  return 1;
-+
+int buildLatticeBilayer(int isHexLattice,
-+
+                        double hexSpacing,
-+
+                        double aLat,
-+
+                        double bLat,
-+
+                        int targetNlipid,
-+
+                        double targetWaterLipidRatio,
-+
+                        double leafSpacing){
-+
-+
+  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
-+
-+
+  double waterCell[3];
-+
-+
+  double *posX, *posY, *posZ;
-+
+  double pos[3], posA[3], posB[3];
-+
-+
+  const double waterFudge = 5.0;
-+
-+
+  int i,j,k,l;
-+
+  int nAtoms, atomIndex, molIndex, molID;
-+
+  int* molSeq;
-+
+  int* molMap;
-+
+  int* molStart;
-+
+  int testTot, done;
-+
+  int nCells, nCellsX, nCellsY, nCellsZ;
-+
+  int nx, ny;
-+
+  double boxX, boxY, boxZ;
-+
+  double unitVector[3];
-+
+  int which;
-+
+  int targetWaters;
-+
-+
-+
-+
+  coord testSite;
-+
-+
+  Atom** atoms;
-+
+  SimInfo* simnfo;
-+
+  SimState* theConfig;
-+
+  DumpWriter* writer;
-+
-+
+  MoleculeStamp* lipidStamp;
-+
+  MoleculeStamp* waterStamp;
-+
+  MoLocator *lipidLocate;
-+
+  MoLocator *waterLocate;
-+
+  int foundLipid, foundWater;
-+
+  int nLipids, lipidNatoms, nWaters, waterNatoms;
-+
-+
+  srand48( RAND_SEED );
-+
-+
+  // create the simInfo objects
-+
-+
+  simnfo = new SimInfo;
-+
-+
+  // 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];
-+
+      lipidNatoms = lipidStamp->getNAtoms();
-+
+    }
-+
+    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
-+
-+
+      foundWater = 1;
-+
-+
+      waterStamp = bsInfo.compStamps[i];
-+
+      nWaters = bsInfo.componentsNmol[i];
-+
+      waterNatoms = waterStamp->getNAtoms();
-+
+    }
-+
+  }
-+
+  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 Molocator arrays
-+
-+
+  lipidLocate = new MoLocator( lipidStamp );
-+
+  waterLocate = new MoLocator( waterStamp );
-+
-+
-+
+  // set up the bilayer leaves
-+
-+
+  if (isHexLattice) {
-+
+    aLat = sqrt(3.0)*hexSpacing;
-+
+    bLat = hexSpacing;
-+
+  }
-+
-+
+  nCells = (int) sqrt( (double)targetNlipid * bLat / (4.0 * aLat) );
-+
-+
+  nx = nCells;
-+
+  ny = (int) ((double)nCells  * aLat / bLat);
-+
-+
+  boxX = nx * aLat;
-+
+  boxY = ny * bLat;
-+
-+
+  nLipids = 4 * nx * ny;
-+
+  coord* lipidSites = new coord[nLipids];
-+
-+
+  unitVector[0] = 0.0;
-+
+  unitVector[1] = 0.0;
-+
-+
+  which = 0;
-+
-+
+  for (i = 0; i < nx; i++) {
-+
+    for (j = 0; j < ny; j++ ) {
-+
+      for (k = 0; k < 2; k++) {
-+
-+
+        lipidSites[which].pos[0] = (double)i * aLat;
-+
+        lipidSites[which].pos[1] = (double)j * bLat;
-+
+        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
-+
-+
+        unitVector[2] = 2.0 * (double)k  - 1.0;
-+
-+
+        getUnitRot( unitVector, lipidSites[which].rot );
-+
-+
+        which++;
-+
-+
+        lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
-+
+        lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
-+
+        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
-+
-+
+        unitVector[2] = 2.0 * (double)k  - 1.0;
-+
-+
+        getUnitRot( unitVector, lipidSites[which].rot );
-+
-+
+        which++;
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  targetWaters = targetWaterLipidRatio * nLipids;
-+
-+
+  // guess the size of the water box
-+
-+
-+
-+
+  nCellsX = (int)ceil(boxX / pow(waterVol, ( 1.0 / 3.0 )) );
-+
+  nCellsY = (int)ceil(boxY / pow(waterVol, ( 1.0 / 3.0 )) );
-+
-+
+  done = 0;
-+
+  nCellsZ = 0;
-+
+  while( !done ){
-+
-+
+    nCellsZ++;
-+
+    testTot = 4 * nCellsX * nCellsY * nCellsZ;
-+
-+
+    if( testTot >= targetWaters ) done = 1;
-+
+  }
-+
-+
+  nWaters = nCellsX * nCellsY * nCellsZ * 4;
-+
-+
+  coord* waterSites = new coord[nWaters];
-+
-+
+  waterCell[0] = boxX / nCellsX;
-+
+  waterCell[1] = boxY / nCellsY;
-+
+  waterCell[2] = 4.0 / (waterRho * waterCell[0] * waterCell[1]);
-+
-+
+  Lattice *myORTHO;
-+
+  myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
-+
+  myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge);
-+
-+
+  boxZ = waterCell[2] * nCellsZ;
-+
-+
+  // create an fcc lattice in the water box.
-+
-+
+  which = 0;
-+
+  for( i=0; i < nCellsX; i++ ){
-+
+    for( j=0; j < nCellsY; j++ ){
-+
+      for( k=0; k < nCellsZ; k++ ){
-+
-+
+        myORTHO->getLatticePoints(&posX, &posY, &posZ, i, j, k);
-+
+        for(l=0; l<4; l++){
-+
+          waterSites[which].pos[0] = posX[l];
-+
+          waterSites[which].pos[1] = posY[l];
-+
+          waterSites[which].pos[2] = posZ[l];
-+
+          which++;
-+
+        }
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  // 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;
-+
+  }
-+
-+
+  theConfig = simnfo->getConfiguration();
-+
+  theConfig->createArrays( nAtoms );
-+
+  simnfo->atoms = new Atom*[nAtoms];
-+
+  atoms = simnfo->atoms;
-+
-+
+  // initialize lipid positions
-+
-+
+  molIndex = 0;
-+
+  for(i=0; i<nLipids; i++ ){
-+
+    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
-+
+                           molStart[molIndex], theConfig );
-+
+    molIndex++;
-+
+  }
-+
-+
+  // initialize the water positions
-+
-+
+  for(i=0; i<nWaters; i++){
-+
-+
+    getRandomRot( waterSites[i].rot );
-+
+    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
-+
+                           molStart[molIndex], theConfig );
-+
+    molIndex++;
-+
+  }
-+
-+
+  // set up the SimInfo object
-+
-+
+  double Hmat[3][3];
-+
-+
+  Hmat[0][0] = boxX;
-+
+  Hmat[0][1] = 0.0;
-+
+  Hmat[0][2] = 0.0;
-+
-+
+  Hmat[1][0] = 0.0;
-+
+  Hmat[1][1] = boxY;
-+
+  Hmat[1][2] = 0.0;
-+
-+
+  Hmat[2][0] = 0.0;
-+
+  Hmat[2][1] = 0.0;
-+
+  Hmat[2][2] = boxZ;
-+
-+
-+
+  bsInfo.boxX = boxX;
-+
+  bsInfo.boxY = boxY;
-+
+  bsInfo.boxZ = boxZ;
-+
-+
+  simnfo->setBoxM( Hmat );
-+
-+
+  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
-+
+  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
-+
-+
+  // set up the writer and write out
-+
-+
+  writer = new DumpWriter( simnfo );
-+
+  writer->writeFinal( 0.0 );
-+
-+
+  return 1;
-+
+}
-+
-+
+void getRandomRot( double rot[3][3] ){
+  double theta, phi, psi;
+  theta = acos( cosTheta );
-+
+  getEulerRot( theta, phi, psi, rot );
-+
+}
-+
-+
-+
+void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
-+
+  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
+  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
+  rot[0][2] = sin(theta) * sin(psi);
+  rot[2][1] = -cos(phi) * sin(theta);
+  rot[2][2] = cos(theta);
-+
-+
-+
+void getUnitRot( double u[3], double rot[3][3] ){
-+
-+
+  double theta, phi, psi;
-+
-+
+  theta = acos(u[2]);
-+
+  phi = atan(u[1] / u[0]);
-+
+  psi = 0.0;
-+
-+
+  getEulerRot( theta, phi, psi, rot );
-+
+}

Diff Legend

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

Comparing trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp (file contents): Revision 700 by chuckv, Mon Aug 18 20:59:47 2003 UTC vs. Revision 817 by gezelter, Fri Oct 24 17:36:18 2003 UTC

Diff Legend

Comparing trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp (file contents):
Revision 700 by chuckv, Mon Aug 18 20:59:47 2003 UTC vs.
Revision 817 by gezelter, Fri Oct 24 17:36:18 2003 UTC