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 = 5.0;


  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 testSite;

  MoleculeStamp* lipidStamp;
  MoleculeStamp* waterStamp;  
  MoLocator *lipidLocate;
  MoLocator *waterLocate
  int foundLipid, foundWater;
  int nLipids, lipiNatoms, nWaters, waterNatoms;
  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 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
  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;
  
  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 ){
      
      lipidSite[i].pos[0] = drand48() * box_x;
      lipidSite[i].pos[1] = drand48() * box_y;
      lipidSite[i].pos[2] = drand48() * box_z;
      
      getRandomRot( lipidSite[i].rot );
      
      ndx = i * lipidNatoms;

      lipidLocate->placeMol( lipidSite[i].pos, lipidSite[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 = waterSite[i].pos[0] - rsaAtoms[j]->getX();
      dy = waterSite[i].pos[1] - rsaAtoms[j]->getY();
      dz = waterSite[i].pos[2] - rsaAtoms[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;
  locate = new MoLocator*[2];
  molSeq = new int[nLipids + nWaters];
  molStart = new int[nLipids + nWaters];  
  
  locate[0] = lipidLocate;
  for(j=0; j<nLipids; j++){
    molSeq[molIndex] = 0;
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += lipidNatoms;
  }

  locate[1] = waterLocate;
  for(j=0; j<nLipids; j++){
    molSeq[molIndex] = 1;
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += waterNatoms;
  }
  
  
  Atom::createArrays( nAtoms );
  atoms = new Atom*[nAtoms];

  
  // initialize lipid positions

  
  // set up the SimInfo object
  
  bsInfo.boxX = box_x;
  bsInfo.boxY = box_y;
  bsInfo.boxZ = box_z;
  
  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;
}

  
}


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:	537
Committed:	Fri May 16 21:37:56 2003 UTC (21 years, 4 months ago) by mmeineke
File size:	16641 byte(s)
Log Message:	still working on the bilayer code
#	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	mmeineke	537	const double lipid_spaceing = 5.0;
53	mmeineke	536
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 testSite;
67
68			MoleculeStamp* lipidStamp;
69			MoleculeStamp* waterStamp;
70			MoLocator *lipidLocate;
71			MoLocator *waterLocate
72			int foundLipid, foundWater;
73	mmeineke	537	int nLipids, lipiNatoms, nWaters, waterNatoms;
74	mmeineke	536	double testBox, maxLength;
75
76			srand48( RAND_SEED );
77
78
79			// set the the lipidStamp
80
81			foundLipid = 0;
82			for(i=0; i<bsInfo.nComponents; i++){
83			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
84
85			foundlipid = 1;
86			lipidStamp = bsInfo.compStamps[i];
87			nLipids = bsInfo.componentsNmol[i];
88			}
89			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
90
91			foundWater = 1;
92			waterStamp = bsInfo.compStamps[i];
93			nWaters = bsInfo.componentsNmol[i];
94			}
95			}
96			if( !foundLipid ){
97			sprintf(painCave.errMsg,
98			"Could not find lipid \"%s\" in the bass file.\n",
99			bsInfo.lipidName );
100			painCave.isFatal = 1;
101			simError();
102			}
103			if( !foundWater ){
104			sprintf(painCave.errMsg,
105	mmeineke	537	"Could not find solvent \"%s\" in the bass file.\n",
106	mmeineke	536	bsInfo.waterName );
107			painCave.isFatal = 1;
108			simError();
109			}
110
111			//create the temp Molocator and atom Arrays
112
113			lipidLocate = new MoLocator( lipidStamp );
114			lipidNatoms = lipidStamp->getNAtoms();
115			maxLength = lipidLocate->getMaxLength();
116
117			waterLocate = new MoLocator( waterStamp );
118	mmeineke	537	waterNatoms = waterStamp->getNatoms();
119	mmeineke	536
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	mmeineke	537	coord* waterSites = new coord[newWaters];
246	mmeineke	536
247			double box_x = waterCell * nCells;
248			double box_y = waterCell * nCells;
249			double box_z = waterCell * nCells;
250
251			// create an fcc lattice in the water box.
252
253			ndx = 0;
254			for( i=0; i < nCells; i++ ){
255			for( j=0; j < nCells; j++ ){
256			for( k=0; k < nCells; k++ ){
257
258	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell;
259			waterSites[ndx].pos[1] = j * waterCell;
260			waterSites[ndx].pos[2] = k * waterCell;
261	mmeineke	536	ndx++;
262
263	mmeineke	537	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	mmeineke	536	ndx++;
267
268	mmeineke	537	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	mmeineke	536	ndx++;
272
273	mmeineke	537	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	mmeineke	536	ndx++;
277			}
278			}
279			}
280
281
282	mmeineke	537	// clear up memory from the test box
283	mmeineke	536
284	mmeineke	537	for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
285	mmeineke	536
286	mmeineke	537	coord* lipidSites = new coord[nLipids];
287	mmeineke	536
288	mmeineke	537	// start a 3D RSA for the for the lipid placements
289
290
291			int reject;
292			int testDX, acceptedDX;
293	mmeineke	536
294	mmeineke	537	rCutSqr = lipid_spaceing * lipid_spaceing;
295	mmeineke	536
296	mmeineke	537	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	mmeineke	536
308	mmeineke	537	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	mmeineke	536
319	mmeineke	537	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	mmeineke	536
335	mmeineke	537	if( reject ){
336	mmeineke	536
337	mmeineke	537	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	mmeineke	536
357	mmeineke	537	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	mmeineke	536
361	mmeineke	537	map( dx, dy, dz, box_x, box_y, box_z );
362	mmeineke	536
363	mmeineke	537	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			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	mmeineke	536
412	mmeineke	537
413			// create the real Atom arrays
414
415	mmeineke	536	nAtoms = 0;
416			molIndex = 0;
417	mmeineke	537	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	mmeineke	536	}
428	mmeineke	537
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	mmeineke	536
437	mmeineke	537
438	mmeineke	536	Atom::createArrays( nAtoms );
439			atoms = new Atom*[nAtoms];
440
441	mmeineke	537
442			// initialize lipid positions
443	mmeineke	536
444	mmeineke	537
445
446
447	mmeineke	536
448	mmeineke	537	// set up the SimInfo object
449	mmeineke	536
450	mmeineke	537	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	mmeineke	536
463	mmeineke	537	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	mmeineke	536
491	mmeineke	537	return 1;
492	mmeineke	536	}
493
494	mmeineke	537
495	mmeineke	536
496
497	mmeineke	537	}
498
499
500
501	mmeineke	536	int Old_buildRandomBilayer( void ){
502
503			int i,j,k;
504			int nAtoms, atomIndex, molIndex, molID;
505			int* molSeq;
506			int* molMap;
507			int* molStart;
508			int* cardDeck;
509			int deckSize;
510			int rSite, rCard;
511			double cell;
512			int nCells, nSites, siteIndex;
513	mmeineke	501	double rot[3][3];
514			double pos[3];
515
516			Atom** atoms;
517	mmeineke	498	SimInfo* simnfo;
518	mmeineke	501	DumpWriter* writer;
519			MoLocator** locate;
520
521			// initialize functions and variables
522	mmeineke	498
523	mmeineke	501	srand48( RAND_SEED );
524			molSeq = NULL;
525	mmeineke	504	molStart = NULL;
526	mmeineke	501	molMap = NULL;
527			cardDeck = NULL;
528			atoms = NULL;
529			locate = NULL;
530			simnfo = NULL;
531			writer = NULL;
532
533			// calculate the number of cells in the fcc box
534
535			nCells = 0;
536			nSites = 0;
537			while( nSites < bsInfo.totNmol ){
538			nCells++;
539	mmeineke	502	nSites = 4.0 * pow( (double)nCells, 3.0 );
540	mmeineke	501	}
541
542
543			// create the molMap and cardDeck arrays
544	mmeineke	498
545	mmeineke	501	molMap = new int[nSites];
546			cardDeck = new int[nSites];
547	mmeineke	498
548	mmeineke	501	for(i=0; i<nSites; i++){
549			molMap[i] = -1;
550			cardDeck[i] = i;
551			}
552
553			// randomly place the molecules on the sites
554	mmeineke	498
555	mmeineke	501	deckSize = nSites;
556			for(i=0; i<bsInfo.totNmol; i++){
557			rCard = (int)( deckSize * drand48() );
558			rSite = cardDeck[rCard];
559			molMap[rSite] = i;
560
561			// book keep the card deck;
562
563			deckSize--;
564			cardDeck[rCard] = cardDeck[deckSize];
565			}
566	mmeineke	498
567
568	mmeineke	501	// create the MoLocator and Atom arrays
569
570			nAtoms = 0;
571			molIndex = 0;
572			locate = new MoLocator*[bsInfo.nComponents];
573			molSeq = new int[bsInfo.totNmol];
574	mmeineke	504	molStart = new int[bsInfo.totNmol];
575	mmeineke	501	for(i=0; i<bsInfo.nComponents; i++){
576			locate[i] = new MoLocator( bsInfo.compStamps[i] );
577			for(j=0; j<bsInfo.componentsNmol[i]; j++){
578			molSeq[molIndex] = i;
579	mmeineke	504	molStart[molIndex] = nAtoms;
580	mmeineke	501	molIndex++;
581	mmeineke	502	nAtoms += bsInfo.compStamps[i]->getNAtoms();
582	mmeineke	501	}
583			}
584
585			Atom::createArrays( nAtoms );
586			atoms = new Atom*[nAtoms];
587
588
589			// place the molecules at each FCC site
590
591			cell = 5.0;
592			for(i=0; i<bsInfo.nComponents; i++){
593			if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
594			}
595	mmeineke	506	cell *= 1.2; // add a little buffer
596
597	mmeineke	502	cell *= M_SQRT2;
598	mmeineke	498
599	mmeineke	501	siteIndex = 0;
600			for(i=0; i<nCells; i++){
601			for(j=0; j<nCells; j++){
602			for(k=0; k<nCells; k++){
603
604			if( molMap[siteIndex] >= 0 ){
605			pos[0] = i * cell;
606			pos[1] = j * cell;
607			pos[2] = k * cell;
608
609			getRandomRot( rot );
610			molID = molSeq[molMap[siteIndex]];
611	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
612			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
613	mmeineke	501	}
614			siteIndex++;
615	mmeineke	498
616	mmeineke	501	if( molMap[siteIndex] >= 0 ){
617			pos[0] = i * cell + (0.5 * cell);
618			pos[1] = j * cell;
619			pos[2] = k * cell + (0.5 * cell);
620	mmeineke	504
621	mmeineke	501	getRandomRot( rot );
622			molID = molSeq[molMap[siteIndex]];
623	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
624			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
625	mmeineke	501	}
626			siteIndex++;
627	mmeineke	498
628	mmeineke	501	if( molMap[siteIndex] >= 0 ){
629			pos[0] = i * cell + (0.5 * cell);
630			pos[1] = j * cell + (0.5 * cell);
631			pos[2] = k * cell;
632
633			getRandomRot( rot );
634			molID = molSeq[molMap[siteIndex]];
635	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
636			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
637	mmeineke	501	}
638			siteIndex++;
639	mmeineke	498
640	mmeineke	501	if( molMap[siteIndex] >= 0 ){
641			pos[0] = i * cell;
642			pos[1] = j * cell + (0.5 * cell);
643			pos[2] = k * cell + (0.5 * cell);
644	mmeineke	504
645	mmeineke	501	getRandomRot( rot );
646			molID = molSeq[molMap[siteIndex]];
647	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
648			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
649	mmeineke	501	}
650			siteIndex++;
651			}
652			}
653			}
654	mmeineke	498
655	mmeineke	501	// set up the SimInfo object
656
657			bsInfo.boxX = nCells * cell;
658			bsInfo.boxY = nCells * cell;
659			bsInfo.boxZ = nCells * cell;
660
661			simnfo = new SimInfo();
662	mmeineke	502	simnfo->n_atoms = nAtoms;
663			simnfo->box_x = bsInfo.boxX;
664			simnfo->box_y = bsInfo.boxY;
665			simnfo->box_z = bsInfo.boxZ;
666	mmeineke	501
667	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
668	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
669	mmeineke	504
670	mmeineke	502	simnfo->atoms = atoms;
671	mmeineke	501
672			// set up the writer and write out
673
674	mmeineke	502	writer = new DumpWriter( simnfo );
675	mmeineke	501	writer->writeFinal();
676	mmeineke	498
677	mmeineke	501	// clean up the memory
678
679			if( molMap != NULL ) delete[] molMap;
680			if( cardDeck != NULL ) delete[] cardDeck;
681			if( locate != NULL ){
682			for(i=0; i<bsInfo.nComponents; i++){
683			delete locate[i];
684	mmeineke	498	}
685	mmeineke	501	delete[] locate;
686			}
687			if( atoms != NULL ){
688			for(i=0; i<nAtoms; i++){
689			delete atoms[i];
690	mmeineke	498	}
691	mmeineke	501	Atom::destroyArrays();
692			delete[] atoms;
693			}
694			if( molSeq != NULL ) delete[] molSeq;
695			if( simnfo != NULL ) delete simnfo;
696			if( writer != NULL ) delete writer;
697	mmeineke	498
698	mmeineke	501	return 1;
699			}
700	mmeineke	498
701
702	mmeineke	501	void getRandomRot( double rot[3][3] ){
703	mmeineke	498
704	mmeineke	501	double theta, phi, psi;
705			double cosTheta;
706	mmeineke	498
707	mmeineke	501	// select random phi, psi, and cosTheta
708	mmeineke	498
709	mmeineke	501	phi = 2.0 * M_PI * drand48();
710			psi = 2.0 * M_PI * drand48();
711			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
712	mmeineke	498
713	mmeineke	501	theta = acos( cosTheta );
714	mmeineke	498
715	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
716			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
717	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
718	mmeineke	501
719			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
720			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
721			rot[1][2] = sin(theta) * cos(psi);
722
723			rot[2][0] = sin(phi) * sin(theta);
724			rot[2][1] = -cos(phi) * sin(theta);
725			rot[2][2] = cos(theta);
726	mmeineke	498	}
727	mmeineke	501
728	mmeineke	536
729
730			void map( double &x, double &y, double &z,
731			double boxX, double boxY, double boxZ ){
732
733			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
734			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
735
736			if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
737			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
738
739			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
740			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
741			}