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 / waterRho; // 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, l;
  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;
  DumpWriter* writer;

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


  // 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 = 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
  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++ ){
        
        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 ){
      
      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 );
      
      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 = waterSites[i].pos[0] - atoms[j]->getX();
      dy = waterSites[i].pos[1] - atoms[j]->getY();
      dz = waterSites[i].pos[2] - atoms[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;
  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;
  }
  
  
  Atom::createArrays( nAtoms );
  atoms = new Atom*[nAtoms];

  
  // initialize lipid positions

  molIndex = 0;
  for(i=0; i<nLipids; i++ ){
    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
                           molStart[molIndex] );
    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] );
      molIndex++;
    }
  }  

  // 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:	538
Committed:	Sat May 17 16:57:08 2003 UTC (21 years, 4 months ago) by mmeineke
File size:	17040 byte(s)
Log Message:	all seems to be working
#	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	mmeineke	538	const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
49	mmeineke	536	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	538	int i,j,k, l;
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	mmeineke	538	Atom** atoms;
69			SimInfo* simnfo;
70			DumpWriter* writer;
71
72	mmeineke	536	MoleculeStamp* lipidStamp;
73			MoleculeStamp* waterStamp;
74			MoLocator *lipidLocate;
75	mmeineke	538	MoLocator *waterLocate;
76	mmeineke	536	int foundLipid, foundWater;
77	mmeineke	538	int nLipids, lipidNatoms, nWaters, waterNatoms;
78	mmeineke	536	double testBox, maxLength;
79
80			srand48( RAND_SEED );
81
82
83			// set the the lipidStamp
84
85			foundLipid = 0;
86	mmeineke	538	foundWater = 0;
87	mmeineke	536	for(i=0; i<bsInfo.nComponents; i++){
88			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
89
90	mmeineke	538	foundLipid = 1;
91	mmeineke	536	lipidStamp = bsInfo.compStamps[i];
92			nLipids = bsInfo.componentsNmol[i];
93			}
94			if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
95
96			foundWater = 1;
97	mmeineke	538
98	mmeineke	536	waterStamp = bsInfo.compStamps[i];
99			nWaters = bsInfo.componentsNmol[i];
100			}
101			}
102			if( !foundLipid ){
103			sprintf(painCave.errMsg,
104			"Could not find lipid \"%s\" in the bass file.\n",
105			bsInfo.lipidName );
106			painCave.isFatal = 1;
107			simError();
108			}
109			if( !foundWater ){
110			sprintf(painCave.errMsg,
111	mmeineke	537	"Could not find solvent \"%s\" in the bass file.\n",
112	mmeineke	536	bsInfo.waterName );
113			painCave.isFatal = 1;
114			simError();
115			}
116
117			//create the temp Molocator and atom Arrays
118
119			lipidLocate = new MoLocator( lipidStamp );
120			lipidNatoms = lipidStamp->getNAtoms();
121			maxLength = lipidLocate->getMaxLength();
122
123			waterLocate = new MoLocator( waterStamp );
124	mmeineke	538	waterNatoms = waterStamp->getNAtoms();
125	mmeineke	536
126			nAtoms = nLipids * lipidNatoms;
127
128			Atom::createArrays( nAtoms );
129			atoms = new Atom*[nAtoms];
130
131			// create the test box for initial water displacement
132
133			testBox = maxLength + waterCell * 4.0; // pad with 4 cells
134	mmeineke	538	nCells = (int)( testBox / waterCell + 1.0 );
135	mmeineke	536	int testWaters = 4 * nCells * nCells * nCells;
136
137			double* waterX = new double[testWaters];
138	mmeineke	538	double* waterY = new double[testWaters];
139			double* waterZ = new double[testWaters];
140	mmeineke	536
141			double x0 = 0.0 - ( testBox * 0.5 );
142			double y0 = 0.0 - ( testBox * 0.5 );
143			double z0 = 0.0 - ( testBox * 0.5 );
144
145
146			// create an fcc lattice in the water box.
147
148			int ndx = 0;
149			for( i=0; i < nCells; i++ ){
150			for( j=0; j < nCells; j++ ){
151			for( k=0; k < nCells; k++ ){
152
153			waterX[ndx] = i * waterCell + x0;
154			waterY[ndx] = j * waterCell + y0;
155			waterZ[ndx] = k * waterCell + z0;
156			ndx++;
157
158			waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
159			waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
160			waterZ[ndx] = k * waterCell + z0;
161			ndx++;
162
163			waterX[ndx] = i * waterCell + x0;
164			waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
165			waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
166			ndx++;
167
168			waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
169			waterY[ndx] = j * waterCell + y0;
170			waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
171			ndx++;
172			}
173			}
174			}
175
176			// calculate the number of water's displaced by our lipid.
177
178			testSite.rot[0][0] = 1.0;
179			testSite.rot[0][1] = 0.0;
180			testSite.rot[0][2] = 0.0;
181
182			testSite.rot[1][0] = 0.0;
183			testSite.rot[1][1] = 1.0;
184			testSite.rot[1][2] = 0.0;
185
186			testSite.rot[2][0] = 0.0;
187			testSite.rot[2][1] = 0.0;
188			testSite.rot[2][2] = 1.0;
189
190			testSite.pos[0] = 0.0;
191			testSite.pos[1] = 0.0;
192			testSite.pos[2] = 0.0;
193
194			lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
195
196			int *isActive = new int[testWaters];
197			for(i=0; i<testWaters; i++) isActive[i] = 1;
198
199			int n_deleted = 0;
200			double dx, dy, dz;
201			double dx2, dy2, dz2, dSqr;
202			double rCutSqr = water_padding * water_padding;
203
204			for(i=0; ( (i<testWaters) && isActive[i] ); i++){
205			for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
206
207			dx = waterX[i] - atoms[j]->getX();
208			dy = waterY[i] - atoms[j]->getY();
209			dz = waterZ[i] - atoms[j]->getZ();
210
211			map( dx, dy, dz, testBox, testBox, testBox );
212
213			dx2 = dx * dx;
214			dy2 = dy * dy;
215			dz2 = dz * dz;
216
217			dSqr = dx2 + dy2 + dz2;
218			if( dSqr < rCutSqr ){
219			isActive[i] = 0;
220			n_deleted++;
221			}
222			}
223			}
224
225			int targetWaters = nWaters + n_deleted * nLipids;
226
227			// find the best box size for the sim
228
229			int testTot;
230			int done = 0;
231			ndx = 0;
232			while( !done ){
233
234			ndx++;
235			testTot = 4 * ndx * ndx * ndx;
236
237			if( testTot >= targetWaters ) done = 1;
238			}
239
240			nCells = ndx;
241
242
243			// create the new water box to the new specifications
244
245			int newWaters = nCells * nCells * nCells * 4;
246
247			delete[] waterX;
248			delete[] waterY;
249			delete[] waterZ;
250
251	mmeineke	537	coord* waterSites = new coord[newWaters];
252	mmeineke	536
253			double box_x = waterCell * nCells;
254			double box_y = waterCell * nCells;
255			double box_z = waterCell * nCells;
256
257			// create an fcc lattice in the water box.
258
259			ndx = 0;
260			for( i=0; i < nCells; i++ ){
261			for( j=0; j < nCells; j++ ){
262			for( k=0; k < nCells; k++ ){
263
264	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell;
265			waterSites[ndx].pos[1] = j * waterCell;
266			waterSites[ndx].pos[2] = k * waterCell;
267	mmeineke	536	ndx++;
268
269	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
270			waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
271			waterSites[ndx].pos[2] = k * waterCell;
272	mmeineke	536	ndx++;
273
274	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell;
275			waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
276			waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
277	mmeineke	536	ndx++;
278
279	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
280			waterSites[ndx].pos[1] = j * waterCell;
281			waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
282	mmeineke	536	ndx++;
283			}
284			}
285			}
286
287
288	mmeineke	537	// clear up memory from the test box
289	mmeineke	536
290	mmeineke	537	for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
291	mmeineke	536
292	mmeineke	537	coord* lipidSites = new coord[nLipids];
293	mmeineke	536
294	mmeineke	537	// start a 3D RSA for the for the lipid placements
295
296
297			int reject;
298			int testDX, acceptedDX;
299	mmeineke	536
300	mmeineke	537	rCutSqr = lipid_spaceing * lipid_spaceing;
301	mmeineke	536
302	mmeineke	537	for(i=0; i<nLipids; i++ ){
303			done = 0;
304			while( !done ){
305
306	mmeineke	538	lipidSites[i].pos[0] = drand48() * box_x;
307			lipidSites[i].pos[1] = drand48() * box_y;
308			lipidSites[i].pos[2] = drand48() * box_z;
309	mmeineke	537
310	mmeineke	538	getRandomRot( lipidSites[i].rot );
311	mmeineke	537
312			ndx = i * lipidNatoms;
313	mmeineke	536
314	mmeineke	538	lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
315			ndx );
316	mmeineke	537
317			reject = 0;
318			for( j=0; !reject && j<i; j++){
319			for(k=0; !reject && k<lipidNatoms; k++){
320
321			acceptedDX = j*lipidNatoms + k;
322			for(l=0; !reject && l<lipidNatoms; l++){
323
324			testDX = ndx + l;
325	mmeineke	536
326	mmeineke	537	dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
327			dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
328			dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
329
330			map( dx, dy, dz, box_x, box_y, box_z );
331
332			dx2 = dx * dx;
333			dy2 = dy * dy;
334			dz2 = dz * dz;
335
336			dSqr = dx2 + dy2 + dz2;
337			if( dSqr < rCutSqr ) reject = 1;
338			}
339			}
340			}
341	mmeineke	536
342	mmeineke	537	if( reject ){
343	mmeineke	536
344	mmeineke	537	for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
345			}
346			else{
347			done = 1;
348			std::cout << i << " has been accepted\n";
349			}
350			}
351			}
352
353			// cut out the waters that overlap with the lipids.
354
355			delete[] isActive;
356			isActive = new int[newWaters];
357			for(i=0; i<newWaters; i++) isActive[i] = 1;
358			int n_active = newWaters;
359			rCutSqr = water_padding * water_padding;
360
361			for(i=0; ( (i<newWaters) && isActive[i] ); i++){
362			for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
363	mmeineke	536
364	mmeineke	538	dx = waterSites[i].pos[0] - atoms[j]->getX();
365			dy = waterSites[i].pos[1] - atoms[j]->getY();
366			dz = waterSites[i].pos[2] - atoms[j]->getZ();
367	mmeineke	536
368	mmeineke	537	map( dx, dy, dz, box_x, box_y, box_z );
369	mmeineke	536
370	mmeineke	537	dx2 = dx * dx;
371			dy2 = dy * dy;
372			dz2 = dz * dz;
373
374			dSqr = dx2 + dy2 + dz2;
375			if( dSqr < rCutSqr ){
376			isActive[i] = 0;
377			n_active--;
378			}
379			}
380			}
381
382			if( n_active < nWaters ){
383
384			sprintf( painCave.errMsg,
385			"Too many waters were removed, edit code and try again.\n" );
386
387			painCave.isFatal = 1;
388			simError();
389			}
390
391			int quickKill;
392			while( n_active > nWaters ){
393
394			quickKill = (int)(drand48()*newWaters);
395
396			if( isActive[quickKill] ){
397			isActive[quickKill] = 0;
398			n_active--;
399			}
400			}
401
402			if( n_active != nWaters ){
403
404			sprintf( painCave.errMsg,
405			"QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
406			n_active, nWaters );
407			painCave.isFatal = 1;
408			simError();
409			}
410
411			// clean up our messes before building the final system.
412
413			for(i=0; i<nAtoms; i++){
414
415			delete atoms[i];
416			}
417			Atom::destroyArrays();
418	mmeineke	536
419	mmeineke	537
420			// create the real Atom arrays
421
422	mmeineke	536	nAtoms = 0;
423			molIndex = 0;
424	mmeineke	537	molStart = new int[nLipids + nWaters];
425
426			for(j=0; j<nLipids; j++){
427			molStart[molIndex] = nAtoms;
428			molIndex++;
429			nAtoms += lipidNatoms;
430	mmeineke	536	}
431	mmeineke	537
432	mmeineke	538	for(j=0; j<nWaters; j++){
433	mmeineke	537	molStart[molIndex] = nAtoms;
434			molIndex++;
435			nAtoms += waterNatoms;
436			}
437	mmeineke	536
438	mmeineke	537
439	mmeineke	536	Atom::createArrays( nAtoms );
440			atoms = new Atom*[nAtoms];
441
442	mmeineke	537
443			// initialize lipid positions
444	mmeineke	536
445	mmeineke	538	molIndex = 0;
446			for(i=0; i<nLipids; i++ ){
447			lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
448			molStart[molIndex] );
449			molIndex++;
450			}
451	mmeineke	536
452	mmeineke	538	// initialize the water positions
453
454			for(i=0; i<newWaters; i++){
455
456			if( isActive[i] ){
457
458			getRandomRot( waterSites[i].rot );
459			waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
460			molStart[molIndex] );
461			molIndex++;
462			}
463			}
464
465	mmeineke	537	// set up the SimInfo object
466	mmeineke	536
467	mmeineke	537	bsInfo.boxX = box_x;
468			bsInfo.boxY = box_y;
469			bsInfo.boxZ = box_z;
470
471			simnfo = new SimInfo();
472			simnfo->n_atoms = nAtoms;
473			simnfo->box_x = bsInfo.boxX;
474			simnfo->box_y = bsInfo.boxY;
475			simnfo->box_z = bsInfo.boxZ;
476
477			sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
478			sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
479	mmeineke	536
480	mmeineke	537	simnfo->atoms = atoms;
481
482			// set up the writer and write out
483
484			writer = new DumpWriter( simnfo );
485			writer->writeFinal();
486
487			// clean up the memory
488
489	mmeineke	538	// if( molMap != NULL ) delete[] molMap;
490			// if( cardDeck != NULL ) delete[] cardDeck;
491			// if( locate != NULL ){
492			// for(i=0; i<bsInfo.nComponents; i++){
493			// delete locate[i];
494			// }
495			// delete[] locate;
496			// }
497			// if( atoms != NULL ){
498			// for(i=0; i<nAtoms; i++){
499			// delete atoms[i];
500			// }
501			// Atom::destroyArrays();
502			// delete[] atoms;
503			// }
504			// if( molSeq != NULL ) delete[] molSeq;
505			// if( simnfo != NULL ) delete simnfo;
506			// if( writer != NULL ) delete writer;
507	mmeineke	536
508	mmeineke	537	return 1;
509	mmeineke	536	}
510
511	mmeineke	537
512	mmeineke	536
513			int Old_buildRandomBilayer( void ){
514
515			int i,j,k;
516			int nAtoms, atomIndex, molIndex, molID;
517			int* molSeq;
518			int* molMap;
519			int* molStart;
520			int* cardDeck;
521			int deckSize;
522			int rSite, rCard;
523			double cell;
524			int nCells, nSites, siteIndex;
525	mmeineke	501	double rot[3][3];
526			double pos[3];
527
528			Atom** atoms;
529	mmeineke	498	SimInfo* simnfo;
530	mmeineke	501	DumpWriter* writer;
531			MoLocator** locate;
532
533			// initialize functions and variables
534	mmeineke	498
535	mmeineke	501	srand48( RAND_SEED );
536			molSeq = NULL;
537	mmeineke	504	molStart = NULL;
538	mmeineke	501	molMap = NULL;
539			cardDeck = NULL;
540			atoms = NULL;
541			locate = NULL;
542			simnfo = NULL;
543			writer = NULL;
544
545			// calculate the number of cells in the fcc box
546
547			nCells = 0;
548			nSites = 0;
549			while( nSites < bsInfo.totNmol ){
550			nCells++;
551	mmeineke	502	nSites = 4.0 * pow( (double)nCells, 3.0 );
552	mmeineke	501	}
553
554
555			// create the molMap and cardDeck arrays
556	mmeineke	498
557	mmeineke	501	molMap = new int[nSites];
558			cardDeck = new int[nSites];
559	mmeineke	498
560	mmeineke	501	for(i=0; i<nSites; i++){
561			molMap[i] = -1;
562			cardDeck[i] = i;
563			}
564
565			// randomly place the molecules on the sites
566	mmeineke	498
567	mmeineke	501	deckSize = nSites;
568			for(i=0; i<bsInfo.totNmol; i++){
569			rCard = (int)( deckSize * drand48() );
570			rSite = cardDeck[rCard];
571			molMap[rSite] = i;
572
573			// book keep the card deck;
574
575			deckSize--;
576			cardDeck[rCard] = cardDeck[deckSize];
577			}
578	mmeineke	498
579
580	mmeineke	501	// create the MoLocator and Atom arrays
581
582			nAtoms = 0;
583			molIndex = 0;
584			locate = new MoLocator*[bsInfo.nComponents];
585			molSeq = new int[bsInfo.totNmol];
586	mmeineke	504	molStart = new int[bsInfo.totNmol];
587	mmeineke	501	for(i=0; i<bsInfo.nComponents; i++){
588			locate[i] = new MoLocator( bsInfo.compStamps[i] );
589			for(j=0; j<bsInfo.componentsNmol[i]; j++){
590			molSeq[molIndex] = i;
591	mmeineke	504	molStart[molIndex] = nAtoms;
592	mmeineke	501	molIndex++;
593	mmeineke	502	nAtoms += bsInfo.compStamps[i]->getNAtoms();
594	mmeineke	501	}
595			}
596
597			Atom::createArrays( nAtoms );
598			atoms = new Atom*[nAtoms];
599
600
601			// place the molecules at each FCC site
602
603			cell = 5.0;
604			for(i=0; i<bsInfo.nComponents; i++){
605			if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
606			}
607	mmeineke	506	cell *= 1.2; // add a little buffer
608
609	mmeineke	502	cell *= M_SQRT2;
610	mmeineke	498
611	mmeineke	501	siteIndex = 0;
612			for(i=0; i<nCells; i++){
613			for(j=0; j<nCells; j++){
614			for(k=0; k<nCells; k++){
615
616			if( molMap[siteIndex] >= 0 ){
617			pos[0] = i * cell;
618			pos[1] = j * cell;
619			pos[2] = k * cell;
620
621			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;
631			pos[2] = k * cell + (0.5 * cell);
632	mmeineke	504
633	mmeineke	501	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 + (0.5 * cell);
642			pos[1] = j * cell + (0.5 * cell);
643			pos[2] = k * cell;
644
645			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	mmeineke	498
652	mmeineke	501	if( molMap[siteIndex] >= 0 ){
653			pos[0] = i * cell;
654			pos[1] = j * cell + (0.5 * cell);
655			pos[2] = k * cell + (0.5 * cell);
656	mmeineke	504
657	mmeineke	501	getRandomRot( rot );
658			molID = molSeq[molMap[siteIndex]];
659	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
660			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
661	mmeineke	501	}
662			siteIndex++;
663			}
664			}
665			}
666	mmeineke	498
667	mmeineke	501	// set up the SimInfo object
668
669			bsInfo.boxX = nCells * cell;
670			bsInfo.boxY = nCells * cell;
671			bsInfo.boxZ = nCells * cell;
672
673			simnfo = new SimInfo();
674	mmeineke	502	simnfo->n_atoms = nAtoms;
675			simnfo->box_x = bsInfo.boxX;
676			simnfo->box_y = bsInfo.boxY;
677			simnfo->box_z = bsInfo.boxZ;
678	mmeineke	501
679	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
680	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
681	mmeineke	504
682	mmeineke	502	simnfo->atoms = atoms;
683	mmeineke	501
684			// set up the writer and write out
685
686	mmeineke	502	writer = new DumpWriter( simnfo );
687	mmeineke	501	writer->writeFinal();
688	mmeineke	498
689	mmeineke	501	// clean up the memory
690
691			if( molMap != NULL ) delete[] molMap;
692			if( cardDeck != NULL ) delete[] cardDeck;
693			if( locate != NULL ){
694			for(i=0; i<bsInfo.nComponents; i++){
695			delete locate[i];
696	mmeineke	498	}
697	mmeineke	501	delete[] locate;
698			}
699			if( atoms != NULL ){
700			for(i=0; i<nAtoms; i++){
701			delete atoms[i];
702	mmeineke	498	}
703	mmeineke	501	Atom::destroyArrays();
704			delete[] atoms;
705			}
706			if( molSeq != NULL ) delete[] molSeq;
707			if( simnfo != NULL ) delete simnfo;
708			if( writer != NULL ) delete writer;
709	mmeineke	498
710	mmeineke	501	return 1;
711			}
712	mmeineke	498
713
714	mmeineke	501	void getRandomRot( double rot[3][3] ){
715	mmeineke	498
716	mmeineke	501	double theta, phi, psi;
717			double cosTheta;
718	mmeineke	498
719	mmeineke	501	// select random phi, psi, and cosTheta
720	mmeineke	498
721	mmeineke	501	phi = 2.0 * M_PI * drand48();
722			psi = 2.0 * M_PI * drand48();
723			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
724	mmeineke	498
725	mmeineke	501	theta = acos( cosTheta );
726	mmeineke	498
727	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
728			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
729	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
730	mmeineke	501
731			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
732			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
733			rot[1][2] = sin(theta) * cos(psi);
734
735			rot[2][0] = sin(phi) * sin(theta);
736			rot[2][1] = -cos(phi) * sin(theta);
737			rot[2][2] = cos(theta);
738	mmeineke	498	}
739	mmeineke	501
740	mmeineke	536
741
742			void map( double &x, double &y, double &z,
743			double boxX, double boxY, double boxZ ){
744
745			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
746			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
747
748			if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
749			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
750
751			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
752			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
753			}