utils/sysbuilder/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"

#include "latticeBuilder.hpp"

void buildMap( 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 / waterRho; // volume occupied by 4 waters
  const double waterCell = 4.929; // fcc unit cell length

  Lattice myFCC( FCC_LATTICE_TYPE, waterCell );
  double *posX, *posY, *posZ;
  double pos[3], posA[3], posB[3];

  const double water_padding = 6.0;
  const double lipid_spaceing = 8.0;


  int i,j,k, l, m;
  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;
  SimState* theConfig;
  DumpWriter* writer;

  MoleculeStamp* lipidStamp;
  MoleculeStamp* waterStamp;  
  MoLocator *lipidLocate;
  MoLocator *waterLocate;
  int foundLipid, foundWater;
  int nLipids, lipidNatoms, nWaters, waterNatoms;
  double testBox, maxLength;
  
  srand48( RAND_SEED );


  // create the simInfo objects

  simnfo = new SimInfo[3];


  // 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 = lipidNatoms;

  simnfo[0].n_atoms = nAtoms;
  simnfo[0].atoms=new Atom*[nAtoms];
 
  theConfig = simnfo[0].getConfiguration();
  theConfig->createArrays( simnfo[0].n_atoms );

  atoms=simnfo[0].atoms;
        

  // 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++ ){

        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
        for(l=0; l<4; l++){
          waterX[ndx]=posX[l];
          waterY[ndx]=posY[l];
          waterZ[ndx]=posZ[l];
          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, theConfig );

  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++){
      
      atoms[j]->getPos( pos );

      dx = waterX[i] - pos[0];
      dy = waterY[i] - pos[1];
      dz = waterZ[i] - pos[2];

      buildMap( 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++ ){

        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
        for(l=0; l<4; l++){
          waterSites[ndx].pos[0] = posX[l];
          waterSites[ndx].pos[1] = posY[l];
          waterSites[ndx].pos[2] = posZ[l];
          ndx++;
        }
      }
    }
  }  

  coord* lipidSites = new coord[nLipids];

  // start a 3D RSA for the for the lipid placements
  
  
  int reject;
  int testDX, acceptedDX;

  nAtoms = nLipids * lipidNatoms;

  simnfo[1].n_atoms = nAtoms;
  simnfo[1].atoms=new Atom*[nAtoms];
 
  theConfig = simnfo[1].getConfiguration();
  theConfig->createArrays( simnfo[1].n_atoms );

  atoms=simnfo[1].atoms;

  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, theConfig );
          
      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;

            atoms[testDX]->getPos( posA );
            atoms[acceptedDX]->getPos( posB );
            
            dx = posA[0] - posB[0];
            dy = posA[1] - posB[1];
            dz = posA[2] - posB[2];
                
            buildMap( 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+1) << " 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++){

      atoms[j]->getPos( pos );

      dx = waterSites[i].pos[0] - pos[0];
      dy = waterSites[i].pos[1] - pos[1];
      dz = waterSites[i].pos[2] - pos[2];

      buildMap( 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.

  simnfo[0].getConfiguration()->destroyArrays();
  simnfo[1].getConfiguration()->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;
  }
  
  theConfig = simnfo[2].getConfiguration();
  theConfig->createArrays( nAtoms );
  simnfo[2].atoms = new Atom*[nAtoms];
  atoms = simnfo[2].atoms;
  simnfo[2].n_atoms = nAtoms;
  
  // 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<newWaters; i++){
        
    if( isActive[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] = box_x;
  Hmat[0][1] = 0.0;
  Hmat[0][2] = 0.0;

  Hmat[1][0] = 0.0;
  Hmat[1][1] = box_y;
  Hmat[1][2] = 0.0;

  Hmat[2][0] = 0.0;
  Hmat[2][1] = 0.0;
  Hmat[2][2] = box_z;
  

  bsInfo.boxX = box_x;
  bsInfo.boxY = box_y;
  bsInfo.boxZ = box_z;
  
  simnfo[2].setBoxM( Hmat );

  sprintf( simnfo[2].sampleName, "%s.dump", bsInfo.outPrefix );
  sprintf( simnfo[2].finalName, "%s.init", bsInfo.outPrefix );

  // set up the writer and write out
  
  writer = new DumpWriter( &simnfo[2] );
  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;
}

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