OOPSE/utils/bilayerSys.cpp

#include <iostream>

#include <cstdlib>
#include <cstring>
#include <cmath>

#include "simError.h"
#include "SimInfo.hpp"
#include "ReadWrite.hpp"

#include "MoLocator.hpp"
#include "sysBuild.hpp"
#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 buildBilayer( int isRandom ){
  
  if( isRandom ){
    return buildRandomBilayer();
  }
  else{
    sprintf( painCave.errMsg,
             "Cannot currently create a non-random bilayer.\n" );
    painCave.isFatal = 1;
    simError();
    return 0;
  }
}


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;
  int* molMap;
  int* molStart;
  int* cardDeck;
  int deckSize;
  int rSite, rCard;
  double cell;
  int nCells, nSites, siteIndex;
  double rot[3][3];
  double pos[3];

  Atom** atoms;
  SimInfo* simnfo;
  DumpWriter* writer;
  MoLocator** locate;

  // initialize functions and variables
  
  srand48( RAND_SEED );
  molSeq = NULL;
  molStart = NULL;
  molMap = NULL;
  cardDeck = NULL;
  atoms = NULL;
  locate = NULL;
  simnfo = NULL;
  writer = NULL;

  // calculate the number of cells in the fcc box

  nCells = 0;
  nSites = 0;
  while( nSites < bsInfo.totNmol ){
    nCells++;
    nSites = 4.0 * pow( (double)nCells, 3.0 );
  }

 
  // create the molMap and cardDeck arrays
  
  molMap = new int[nSites];
  cardDeck = new int[nSites];
  
  for(i=0; i<nSites; i++){
    molMap[i] = -1;
    cardDeck[i] = i;
  }

  // randomly place the molecules on the sites
  
  deckSize = nSites;
  for(i=0; i<bsInfo.totNmol; i++){
    rCard = (int)( deckSize * drand48() );
    rSite = cardDeck[rCard];
    molMap[rSite] = i;

    // book keep the card deck;
    
    deckSize--;
    cardDeck[rCard] = cardDeck[deckSize];
  }
  
  
  // create the 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];
  
  
  // place the molecules at each FCC site
  
  cell = 5.0;
  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;
  for(i=0; i<nCells; i++){
    for(j=0; j<nCells; j++){
      for(k=0; k<nCells; k++){
        
        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell;
          pos[1] = j * cell;
          pos[2] = k * cell;
          
          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell + (0.5 * cell);
          pos[1] = j * cell;
          pos[2] = k * cell + (0.5 * cell);

          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell + (0.5 * cell);
          pos[1] = j * cell + (0.5 * cell);
          pos[2] = k * cell;
          
          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell;
          pos[1] = j * cell + (0.5 * cell);
          pos[2] = k * cell + (0.5 * cell);

          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;
      }
    }
  }

  // set up the SimInfo object
  
  bsInfo.boxX = nCells * cell;
  bsInfo.boxY = nCells * cell;
  bsInfo.boxZ = nCells * cell;
  
  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;
  simnfo->box_x = bsInfo.boxX;
  simnfo->box_y = bsInfo.boxY;
  simnfo->box_z = bsInfo.boxZ;
  
  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();
    
  // 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;
}


void getRandomRot( double rot[3][3] ){

  double theta, phi, psi;
  double cosTheta;

  // select random phi, psi, and cosTheta

  phi = 2.0 * M_PI * drand48();
  psi = 2.0 * M_PI * drand48();
  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1

  theta = acos( cosTheta );

  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[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  rot[1][2] = sin(theta) * cos(psi);

  rot[2][0] = sin(phi) * sin(theta);
  rot[2][1] = -cos(phi) * sin(theta);
  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));
}
Revision:	536
Committed:	Fri May 16 14:28:27 2003 UTC (21 years, 4 months ago) by mmeineke
File size:	12645 byte(s)
Log Message:	doing some work to overhaul sysbuild.
#	User	Rev	Content
1	mmeineke	504	#include <iostream>
2
3	mmeineke	498	#include <cstdlib>
4			#include <cstring>
5			#include <cmath>
6
7			#include "simError.h"
8			#include "SimInfo.hpp"
9			#include "ReadWrite.hpp"
10
11	mmeineke	501	#include "MoLocator.hpp"
12	mmeineke	498	#include "sysBuild.hpp"
13			#include "bilayerSys.hpp"
14
15
16	mmeineke	536
17			void map( double &x, double &y, double &z,
18			double boxX, double boxY, double boxZ );
19
20	mmeineke	501	int buildRandomBilayer( void );
21	mmeineke	498
22	mmeineke	501	void getRandomRot( double rot[3][3] );
23	mmeineke	498
24	mmeineke	501	int buildBilayer( int isRandom ){
25	mmeineke	498
26	mmeineke	501	if( isRandom ){
27			return buildRandomBilayer();
28	mmeineke	498	}
29			else{
30			sprintf( painCave.errMsg,
31			"Cannot currently create a non-random bilayer.\n" );
32			painCave.isFatal = 1;
33			simError();
34	mmeineke	501	return 0;
35	mmeineke	498	}
36			}
37
38
39	mmeineke	501	int buildRandomBilayer( void ){
40	mmeineke	498
41	mmeineke	536	typedef struct{
42			double rot[3][3];
43			double pos[3];
44			} coord;
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
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;
53
54
55	mmeineke	501	int i,j,k;
56	mmeineke	502	int nAtoms, atomIndex, molIndex, molID;
57	mmeineke	501	int* molSeq;
58			int* molMap;
59	mmeineke	504	int* molStart;
60	mmeineke	501	int* cardDeck;
61			int deckSize;
62			int rSite, rCard;
63			double cell;
64			int nCells, nSites, siteIndex;
65	mmeineke	536
66			coord *siteArray;
67			coord testSite;
68
69			MoleculeStamp* lipidStamp;
70			MoleculeStamp* waterStamp;
71			MoLocator *lipidLocate;
72			MoLocator *waterLocate
73			int foundLipid, foundWater;
74			int nLipids, lipiNatoms, nWaters;
75			double testBox, maxLength;
76
77			srand48( RAND_SEED );
78
79
80			// set the the lipidStamp
81
82			foundLipid = 0;
83			for(i=0; i<bsInfo.nComponents; i++){
84			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
85
86			foundlipid = 1;
87			lipidStamp = bsInfo.compStamps[i];
88			nLipids = bsInfo.componentsNmol[i];
89			}
90			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
91
92			foundWater = 1;
93			waterStamp = bsInfo.compStamps[i];
94			nWaters = bsInfo.componentsNmol[i];
95			}
96			}
97			if( !foundLipid ){
98			sprintf(painCave.errMsg,
99			"Could not find lipid \"%s\" in the bass file.\n",
100			bsInfo.lipidName );
101			painCave.isFatal = 1;
102			simError();
103			}
104			if( !foundWater ){
105			sprintf(painCave.errMsg,
106			"Could not find water \"%s\" in the bass file.\n",
107			bsInfo.waterName );
108			painCave.isFatal = 1;
109			simError();
110			}
111
112			//create the temp Molocator and atom Arrays
113
114			lipidLocate = new MoLocator( lipidStamp );
115			lipidNatoms = lipidStamp->getNAtoms();
116			maxLength = lipidLocate->getMaxLength();
117
118			waterLocate = new MoLocator( waterStamp );
119
120			nAtoms = nLipids * lipidNatoms;
121
122			Atom::createArrays( nAtoms );
123			atoms = new Atom*[nAtoms];
124
125			// create the test box for initial water displacement
126
127			testBox = maxLength + waterCell * 4.0; // pad with 4 cells
128			int nCells = (int)( testBox / waterCell + 1.0 );
129			int testWaters = 4 * nCells * nCells * nCells;
130
131			double* waterX = new double[testWaters];
132			double* waterX = new double[testWaters];
133			double* waterX = new double[testWaters];
134
135			double x0 = 0.0 - ( testBox * 0.5 );
136			double y0 = 0.0 - ( testBox * 0.5 );
137			double z0 = 0.0 - ( testBox * 0.5 );
138
139
140			// create an fcc lattice in the water box.
141
142			int ndx = 0;
143			for( i=0; i < nCells; i++ ){
144			for( j=0; j < nCells; j++ ){
145			for( k=0; k < nCells; k++ ){
146
147			waterX[ndx] = i * waterCell + x0;
148			waterY[ndx] = j * waterCell + y0;
149			waterZ[ndx] = k * waterCell + z0;
150			ndx++;
151
152			waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
153			waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
154			waterZ[ndx] = k * waterCell + z0;
155			ndx++;
156
157			waterX[ndx] = i * waterCell + x0;
158			waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
159			waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
160			ndx++;
161
162			waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
163			waterY[ndx] = j * waterCell + y0;
164			waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
165			ndx++;
166			}
167			}
168			}
169
170			// calculate the number of water's displaced by our lipid.
171
172			testSite.rot[0][0] = 1.0;
173			testSite.rot[0][1] = 0.0;
174			testSite.rot[0][2] = 0.0;
175
176			testSite.rot[1][0] = 0.0;
177			testSite.rot[1][1] = 1.0;
178			testSite.rot[1][2] = 0.0;
179
180			testSite.rot[2][0] = 0.0;
181			testSite.rot[2][1] = 0.0;
182			testSite.rot[2][2] = 1.0;
183
184			testSite.pos[0] = 0.0;
185			testSite.pos[1] = 0.0;
186			testSite.pos[2] = 0.0;
187
188			lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
189
190			int *isActive = new int[testWaters];
191			for(i=0; i<testWaters; i++) isActive[i] = 1;
192
193			int n_deleted = 0;
194			double dx, dy, dz;
195			double dx2, dy2, dz2, dSqr;
196			double rCutSqr = water_padding * water_padding;
197
198			for(i=0; ( (i<testWaters) && isActive[i] ); i++){
199			for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
200
201			dx = waterX[i] - atoms[j]->getX();
202			dy = waterY[i] - atoms[j]->getY();
203			dz = waterZ[i] - atoms[j]->getZ();
204
205			map( dx, dy, dz, testBox, testBox, testBox );
206
207			dx2 = dx * dx;
208			dy2 = dy * dy;
209			dz2 = dz * dz;
210
211			dSqr = dx2 + dy2 + dz2;
212			if( dSqr < rCutSqr ){
213			isActive[i] = 0;
214			n_deleted++;
215			}
216			}
217			}
218
219			int targetWaters = nWaters + n_deleted * nLipids;
220
221			// find the best box size for the sim
222
223			int testTot;
224			int done = 0;
225			ndx = 0;
226			while( !done ){
227
228			ndx++;
229			testTot = 4 * ndx * ndx * ndx;
230
231			if( testTot >= targetWaters ) done = 1;
232			}
233
234			nCells = ndx;
235
236
237			// create the new water box to the new specifications
238
239			int newWaters = nCells * nCells * nCells * 4;
240
241			delete[] waterX;
242			delete[] waterY;
243			delete[] waterZ;
244
245			waterX = new double[newWater];
246			waterY = new double[newWater];
247			waterZ = new double[newWater];
248
249			double box_x = waterCell * nCells;
250			double box_y = waterCell * nCells;
251			double box_z = waterCell * nCells;
252
253			// create an fcc lattice in the water box.
254
255			ndx = 0;
256			for( i=0; i < nCells; i++ ){
257			for( j=0; j < nCells; j++ ){
258			for( k=0; k < nCells; k++ ){
259
260			waterX[ndx] = i * waterCell;
261			waterY[ndx] = j * waterCell;
262			waterZ[ndx] = k * waterCell;
263			ndx++;
264
265			waterX[ndx] = i * waterCell + 0.5 * waterCell;
266			waterY[ndx] = j * waterCell + 0.5 * waterCell;
267			waterZ[ndx] = k * waterCell;
268			ndx++;
269
270			waterX[ndx] = i * waterCell;
271			waterY[ndx] = j * waterCell + 0.5 * waterCell;
272			waterZ[ndx] = k * waterCell + 0.5 * waterCell;
273			ndx++;
274
275			waterX[ndx] = i * waterCell + 0.5 * waterCell;
276			waterY[ndx] = j * waterCell;
277			waterZ[ndx] = k * waterCell + 0.5 * waterCell;
278			ndx++;
279			}
280			}
281			}
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296			// create the real MoLocator and Atom arrays
297
298			nAtoms = 0;
299			molIndex = 0;
300			locate = new MoLocator*[bsInfo.nComponents];
301			molSeq = new int[bsInfo.totNmol];
302			molStart = new int[bsInfo.totNmol];
303			for(i=0; i<bsInfo.nComponents; i++){
304			locate[i] = new MoLocator( bsInfo.compStamps[i] );
305			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			}
311			}
312
313			Atom::createArrays( nAtoms );
314			atoms = new Atom*[nAtoms];
315
316
317			// find the width, height, and length of the molecule
318
319
320
321
322			}
323
324
325
326			int Old_buildRandomBilayer( void ){
327
328			int i,j,k;
329			int nAtoms, atomIndex, molIndex, molID;
330			int* molSeq;
331			int* molMap;
332			int* molStart;
333			int* cardDeck;
334			int deckSize;
335			int rSite, rCard;
336			double cell;
337			int nCells, nSites, siteIndex;
338	mmeineke	501	double rot[3][3];
339			double pos[3];
340
341			Atom** atoms;
342	mmeineke	498	SimInfo* simnfo;
343	mmeineke	501	DumpWriter* writer;
344			MoLocator** locate;
345
346			// initialize functions and variables
347	mmeineke	498
348	mmeineke	501	srand48( RAND_SEED );
349			molSeq = NULL;
350	mmeineke	504	molStart = NULL;
351	mmeineke	501	molMap = NULL;
352			cardDeck = NULL;
353			atoms = NULL;
354			locate = NULL;
355			simnfo = NULL;
356			writer = NULL;
357
358			// calculate the number of cells in the fcc box
359
360			nCells = 0;
361			nSites = 0;
362			while( nSites < bsInfo.totNmol ){
363			nCells++;
364	mmeineke	502	nSites = 4.0 * pow( (double)nCells, 3.0 );
365	mmeineke	501	}
366
367
368			// create the molMap and cardDeck arrays
369	mmeineke	498
370	mmeineke	501	molMap = new int[nSites];
371			cardDeck = new int[nSites];
372	mmeineke	498
373	mmeineke	501	for(i=0; i<nSites; i++){
374			molMap[i] = -1;
375			cardDeck[i] = i;
376			}
377
378			// randomly place the molecules on the sites
379	mmeineke	498
380	mmeineke	501	deckSize = nSites;
381			for(i=0; i<bsInfo.totNmol; i++){
382			rCard = (int)( deckSize * drand48() );
383			rSite = cardDeck[rCard];
384			molMap[rSite] = i;
385
386			// book keep the card deck;
387
388			deckSize--;
389			cardDeck[rCard] = cardDeck[deckSize];
390			}
391	mmeineke	498
392
393	mmeineke	501	// create the MoLocator and Atom arrays
394
395			nAtoms = 0;
396			molIndex = 0;
397			locate = new MoLocator*[bsInfo.nComponents];
398			molSeq = new int[bsInfo.totNmol];
399	mmeineke	504	molStart = new int[bsInfo.totNmol];
400	mmeineke	501	for(i=0; i<bsInfo.nComponents; i++){
401			locate[i] = new MoLocator( bsInfo.compStamps[i] );
402			for(j=0; j<bsInfo.componentsNmol[i]; j++){
403			molSeq[molIndex] = i;
404	mmeineke	504	molStart[molIndex] = nAtoms;
405	mmeineke	501	molIndex++;
406	mmeineke	502	nAtoms += bsInfo.compStamps[i]->getNAtoms();
407	mmeineke	501	}
408			}
409
410			Atom::createArrays( nAtoms );
411			atoms = new Atom*[nAtoms];
412
413
414			// place the molecules at each FCC site
415
416			cell = 5.0;
417			for(i=0; i<bsInfo.nComponents; i++){
418			if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
419			}
420	mmeineke	506	cell *= 1.2; // add a little buffer
421
422	mmeineke	502	cell *= M_SQRT2;
423	mmeineke	498
424	mmeineke	501	siteIndex = 0;
425			for(i=0; i<nCells; i++){
426			for(j=0; j<nCells; j++){
427			for(k=0; k<nCells; k++){
428
429			if( molMap[siteIndex] >= 0 ){
430			pos[0] = i * cell;
431			pos[1] = j * cell;
432			pos[2] = k * cell;
433
434			getRandomRot( rot );
435			molID = molSeq[molMap[siteIndex]];
436	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
437			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
438	mmeineke	501	}
439			siteIndex++;
440	mmeineke	498
441	mmeineke	501	if( molMap[siteIndex] >= 0 ){
442			pos[0] = i * cell + (0.5 * cell);
443			pos[1] = j * cell;
444			pos[2] = k * cell + (0.5 * cell);
445	mmeineke	504
446	mmeineke	501	getRandomRot( rot );
447			molID = molSeq[molMap[siteIndex]];
448	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
449			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
450	mmeineke	501	}
451			siteIndex++;
452	mmeineke	498
453	mmeineke	501	if( molMap[siteIndex] >= 0 ){
454			pos[0] = i * cell + (0.5 * cell);
455			pos[1] = j * cell + (0.5 * cell);
456			pos[2] = k * cell;
457
458			getRandomRot( rot );
459			molID = molSeq[molMap[siteIndex]];
460	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
461			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
462	mmeineke	501	}
463			siteIndex++;
464	mmeineke	498
465	mmeineke	501	if( molMap[siteIndex] >= 0 ){
466			pos[0] = i * cell;
467			pos[1] = j * cell + (0.5 * cell);
468			pos[2] = k * cell + (0.5 * cell);
469	mmeineke	504
470	mmeineke	501	getRandomRot( rot );
471			molID = molSeq[molMap[siteIndex]];
472	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
473			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
474	mmeineke	501	}
475			siteIndex++;
476			}
477			}
478			}
479	mmeineke	498
480	mmeineke	501	// set up the SimInfo object
481
482			bsInfo.boxX = nCells * cell;
483			bsInfo.boxY = nCells * cell;
484			bsInfo.boxZ = nCells * cell;
485
486			simnfo = new SimInfo();
487	mmeineke	502	simnfo->n_atoms = nAtoms;
488			simnfo->box_x = bsInfo.boxX;
489			simnfo->box_y = bsInfo.boxY;
490			simnfo->box_z = bsInfo.boxZ;
491	mmeineke	501
492	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
493	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
494	mmeineke	504
495	mmeineke	502	simnfo->atoms = atoms;
496	mmeineke	501
497			// set up the writer and write out
498
499	mmeineke	502	writer = new DumpWriter( simnfo );
500	mmeineke	501	writer->writeFinal();
501	mmeineke	498
502	mmeineke	501	// clean up the memory
503
504			if( molMap != NULL ) delete[] molMap;
505			if( cardDeck != NULL ) delete[] cardDeck;
506			if( locate != NULL ){
507			for(i=0; i<bsInfo.nComponents; i++){
508			delete locate[i];
509	mmeineke	498	}
510	mmeineke	501	delete[] locate;
511			}
512			if( atoms != NULL ){
513			for(i=0; i<nAtoms; i++){
514			delete atoms[i];
515	mmeineke	498	}
516	mmeineke	501	Atom::destroyArrays();
517			delete[] atoms;
518			}
519			if( molSeq != NULL ) delete[] molSeq;
520			if( simnfo != NULL ) delete simnfo;
521			if( writer != NULL ) delete writer;
522	mmeineke	498
523	mmeineke	501	return 1;
524			}
525	mmeineke	498
526
527	mmeineke	501	void getRandomRot( double rot[3][3] ){
528	mmeineke	498
529	mmeineke	501	double theta, phi, psi;
530			double cosTheta;
531	mmeineke	498
532	mmeineke	501	// select random phi, psi, and cosTheta
533	mmeineke	498
534	mmeineke	501	phi = 2.0 * M_PI * drand48();
535			psi = 2.0 * M_PI * drand48();
536			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
537	mmeineke	498
538	mmeineke	501	theta = acos( cosTheta );
539	mmeineke	498
540	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
541			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
542	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
543	mmeineke	501
544			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
545			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
546			rot[1][2] = sin(theta) * cos(psi);
547
548			rot[2][0] = sin(phi) * sin(theta);
549			rot[2][1] = -cos(phi) * sin(theta);
550			rot[2][2] = cos(theta);
551	mmeineke	498	}
552	mmeineke	501
553	mmeineke	536
554
555			void map( double &x, double &y, double &z,
556			double boxX, double boxY, double boxZ ){
557
558			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
559			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
560
561			if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
562			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
563
564			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
565			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
566			}