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 = 6.0;
  const double lipid_spaceing = 8.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;

  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++ ){
        
        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;
  
  double boxVector[3];
  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;
  boxVector[0] = bsInfo.boxX;
  boxVector[1] = bsInfo.boxY;
  boxVector[2] = bsInfo.boxZ;
  simnfo->setBox( boxVector );

  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( 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;
}

  
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;
  
  double boxVector[3];
  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;
  boxVector[0] = bsInfo.boxX;
  boxVector[1] = bsInfo.boxY;
  boxVector[2] = bsInfo.boxZ;
  simnfo->setBox( boxVector );

  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(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 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:	573
Committed:	Thu Jul 3 19:41:26 2003 UTC (21 years, 2 months ago) by mmeineke
File size:	17194 byte(s)
Log Message:	cleaned up the dependecy scripts in the makefiles
#	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	mmeineke	543	const double water_padding = 6.0;
52			const double lipid_spaceing = 8.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	mmeineke	543	targetWaters = (int) ( targetWaters * 1.2 );
228
229	mmeineke	536	// find the best box size for the sim
230
231			int testTot;
232			int done = 0;
233			ndx = 0;
234			while( !done ){
235
236			ndx++;
237			testTot = 4 * ndx * ndx * ndx;
238
239			if( testTot >= targetWaters ) done = 1;
240			}
241
242			nCells = ndx;
243
244
245			// create the new water box to the new specifications
246
247			int newWaters = nCells * nCells * nCells * 4;
248
249			delete[] waterX;
250			delete[] waterY;
251			delete[] waterZ;
252
253	mmeineke	537	coord* waterSites = new coord[newWaters];
254	mmeineke	536
255			double box_x = waterCell * nCells;
256			double box_y = waterCell * nCells;
257			double box_z = waterCell * nCells;
258
259			// create an fcc lattice in the water box.
260
261			ndx = 0;
262			for( i=0; i < nCells; i++ ){
263			for( j=0; j < nCells; j++ ){
264			for( k=0; k < nCells; k++ ){
265
266	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell;
267			waterSites[ndx].pos[1] = j * waterCell;
268			waterSites[ndx].pos[2] = k * waterCell;
269	mmeineke	536	ndx++;
270
271	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
272			waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
273			waterSites[ndx].pos[2] = k * waterCell;
274	mmeineke	536	ndx++;
275
276	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell;
277			waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
278			waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
279	mmeineke	536	ndx++;
280
281	mmeineke	537	waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
282			waterSites[ndx].pos[1] = j * waterCell;
283			waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
284	mmeineke	536	ndx++;
285			}
286			}
287			}
288
289
290	mmeineke	537	// clear up memory from the test box
291	mmeineke	536
292	mmeineke	537	for(i=0; i<lipidNatoms; i++ ) delete atoms[i];
293	mmeineke	536
294	mmeineke	537	coord* lipidSites = new coord[nLipids];
295	mmeineke	536
296	mmeineke	537	// start a 3D RSA for the for the lipid placements
297
298
299			int reject;
300			int testDX, acceptedDX;
301	mmeineke	536
302	mmeineke	537	rCutSqr = lipid_spaceing * lipid_spaceing;
303	mmeineke	536
304	mmeineke	537	for(i=0; i<nLipids; i++ ){
305			done = 0;
306			while( !done ){
307
308	mmeineke	538	lipidSites[i].pos[0] = drand48() * box_x;
309			lipidSites[i].pos[1] = drand48() * box_y;
310			lipidSites[i].pos[2] = drand48() * box_z;
311	mmeineke	537
312	mmeineke	538	getRandomRot( lipidSites[i].rot );
313	mmeineke	537
314			ndx = i * lipidNatoms;
315	mmeineke	536
316	mmeineke	538	lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
317			ndx );
318	mmeineke	537
319			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
326			testDX = ndx + l;
327	mmeineke	536
328	mmeineke	537	dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
329			dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
330			dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
331
332			map( dx, dy, dz, box_x, box_y, box_z );
333
334			dx2 = dx * dx;
335			dy2 = dy * dy;
336			dz2 = dz * dz;
337
338			dSqr = dx2 + dy2 + dz2;
339			if( dSqr < rCutSqr ) reject = 1;
340			}
341			}
342			}
343	mmeineke	536
344	mmeineke	537	if( reject ){
345	mmeineke	536
346	mmeineke	537	for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
347			}
348			else{
349			done = 1;
350			std::cout << i << " has been accepted\n";
351			}
352			}
353			}
354
355			// cut out the waters that overlap with the lipids.
356
357			delete[] isActive;
358			isActive = new int[newWaters];
359			for(i=0; i<newWaters; i++) isActive[i] = 1;
360			int n_active = newWaters;
361			rCutSqr = water_padding * water_padding;
362
363			for(i=0; ( (i<newWaters) && isActive[i] ); i++){
364			for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
365	mmeineke	536
366	mmeineke	538	dx = waterSites[i].pos[0] - atoms[j]->getX();
367			dy = waterSites[i].pos[1] - atoms[j]->getY();
368			dz = waterSites[i].pos[2] - atoms[j]->getZ();
369	mmeineke	536
370	mmeineke	537	map( dx, dy, dz, box_x, box_y, box_z );
371	mmeineke	536
372	mmeineke	537	dx2 = dx * dx;
373			dy2 = dy * dy;
374			dz2 = dz * dz;
375
376			dSqr = dx2 + dy2 + dz2;
377			if( dSqr < rCutSqr ){
378			isActive[i] = 0;
379			n_active--;
380			}
381			}
382			}
383
384			if( n_active < nWaters ){
385
386			sprintf( painCave.errMsg,
387			"Too many waters were removed, edit code and try again.\n" );
388
389			painCave.isFatal = 1;
390			simError();
391			}
392
393			int quickKill;
394			while( n_active > nWaters ){
395
396			quickKill = (int)(drand48()*newWaters);
397
398			if( isActive[quickKill] ){
399			isActive[quickKill] = 0;
400			n_active--;
401			}
402			}
403
404			if( n_active != nWaters ){
405
406			sprintf( painCave.errMsg,
407			"QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
408			n_active, nWaters );
409			painCave.isFatal = 1;
410			simError();
411			}
412
413			// clean up our messes before building the final system.
414
415			for(i=0; i<nAtoms; i++){
416
417			delete atoms[i];
418			}
419			Atom::destroyArrays();
420	mmeineke	536
421	mmeineke	537
422			// create the real Atom arrays
423
424	mmeineke	536	nAtoms = 0;
425			molIndex = 0;
426	mmeineke	537	molStart = new int[nLipids + nWaters];
427
428			for(j=0; j<nLipids; j++){
429			molStart[molIndex] = nAtoms;
430			molIndex++;
431			nAtoms += lipidNatoms;
432	mmeineke	536	}
433	mmeineke	537
434	mmeineke	538	for(j=0; j<nWaters; j++){
435	mmeineke	537	molStart[molIndex] = nAtoms;
436			molIndex++;
437			nAtoms += waterNatoms;
438			}
439	mmeineke	536
440	mmeineke	537
441	mmeineke	536	Atom::createArrays( nAtoms );
442			atoms = new Atom*[nAtoms];
443
444	mmeineke	537
445			// initialize lipid positions
446	mmeineke	536
447	mmeineke	538	molIndex = 0;
448			for(i=0; i<nLipids; i++ ){
449			lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
450			molStart[molIndex] );
451			molIndex++;
452			}
453	mmeineke	536
454	mmeineke	538	// initialize the water positions
455
456			for(i=0; i<newWaters; i++){
457
458			if( isActive[i] ){
459
460			getRandomRot( waterSites[i].rot );
461			waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
462			molStart[molIndex] );
463			molIndex++;
464			}
465			}
466
467	mmeineke	537	// set up the SimInfo object
468	mmeineke	536
469	mmeineke	537	bsInfo.boxX = box_x;
470			bsInfo.boxY = box_y;
471			bsInfo.boxZ = box_z;
472
473	mmeineke	573	double boxVector[3];
474	mmeineke	537	simnfo = new SimInfo();
475			simnfo->n_atoms = nAtoms;
476	mmeineke	573	boxVector[0] = bsInfo.boxX;
477			boxVector[1] = bsInfo.boxY;
478			boxVector[2] = bsInfo.boxZ;
479			simnfo->setBox( boxVector );
480
481	mmeineke	537	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
482			sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
483	mmeineke	536
484	mmeineke	537	simnfo->atoms = atoms;
485
486			// set up the writer and write out
487
488			writer = new DumpWriter( simnfo );
489	mmeineke	573	writer->writeFinal( 0.0 );
490	mmeineke	537
491			// clean up the memory
492
493	mmeineke	538	// if( molMap != NULL ) delete[] molMap;
494			// if( cardDeck != NULL ) delete[] cardDeck;
495			// if( locate != NULL ){
496			// for(i=0; i<bsInfo.nComponents; i++){
497			// delete locate[i];
498			// }
499			// delete[] locate;
500			// }
501			// if( atoms != NULL ){
502			// for(i=0; i<nAtoms; i++){
503			// delete atoms[i];
504			// }
505			// Atom::destroyArrays();
506			// delete[] atoms;
507			// }
508			// if( molSeq != NULL ) delete[] molSeq;
509			// if( simnfo != NULL ) delete simnfo;
510			// if( writer != NULL ) delete writer;
511	mmeineke	536
512	mmeineke	537	return 1;
513	mmeineke	536	}
514
515	mmeineke	537
516	mmeineke	536
517			int Old_buildRandomBilayer( void ){
518
519			int i,j,k;
520			int nAtoms, atomIndex, molIndex, molID;
521			int* molSeq;
522			int* molMap;
523			int* molStart;
524			int* cardDeck;
525			int deckSize;
526			int rSite, rCard;
527			double cell;
528			int nCells, nSites, siteIndex;
529	mmeineke	501	double rot[3][3];
530			double pos[3];
531
532			Atom** atoms;
533	mmeineke	498	SimInfo* simnfo;
534	mmeineke	501	DumpWriter* writer;
535			MoLocator** locate;
536
537			// initialize functions and variables
538	mmeineke	498
539	mmeineke	501	srand48( RAND_SEED );
540			molSeq = NULL;
541	mmeineke	504	molStart = NULL;
542	mmeineke	501	molMap = NULL;
543			cardDeck = NULL;
544			atoms = NULL;
545			locate = NULL;
546			simnfo = NULL;
547			writer = NULL;
548
549			// calculate the number of cells in the fcc box
550
551			nCells = 0;
552			nSites = 0;
553			while( nSites < bsInfo.totNmol ){
554			nCells++;
555	mmeineke	502	nSites = 4.0 * pow( (double)nCells, 3.0 );
556	mmeineke	501	}
557
558
559			// create the molMap and cardDeck arrays
560	mmeineke	498
561	mmeineke	501	molMap = new int[nSites];
562			cardDeck = new int[nSites];
563	mmeineke	498
564	mmeineke	501	for(i=0; i<nSites; i++){
565			molMap[i] = -1;
566			cardDeck[i] = i;
567			}
568
569			// randomly place the molecules on the sites
570	mmeineke	498
571	mmeineke	501	deckSize = nSites;
572			for(i=0; i<bsInfo.totNmol; i++){
573			rCard = (int)( deckSize * drand48() );
574			rSite = cardDeck[rCard];
575			molMap[rSite] = i;
576
577			// book keep the card deck;
578
579			deckSize--;
580			cardDeck[rCard] = cardDeck[deckSize];
581			}
582	mmeineke	498
583
584	mmeineke	501	// create the MoLocator and Atom arrays
585
586			nAtoms = 0;
587			molIndex = 0;
588			locate = new MoLocator*[bsInfo.nComponents];
589			molSeq = new int[bsInfo.totNmol];
590	mmeineke	504	molStart = new int[bsInfo.totNmol];
591	mmeineke	501	for(i=0; i<bsInfo.nComponents; i++){
592			locate[i] = new MoLocator( bsInfo.compStamps[i] );
593			for(j=0; j<bsInfo.componentsNmol[i]; j++){
594			molSeq[molIndex] = i;
595	mmeineke	504	molStart[molIndex] = nAtoms;
596	mmeineke	501	molIndex++;
597	mmeineke	502	nAtoms += bsInfo.compStamps[i]->getNAtoms();
598	mmeineke	501	}
599			}
600
601			Atom::createArrays( nAtoms );
602			atoms = new Atom*[nAtoms];
603
604
605			// place the molecules at each FCC site
606
607			cell = 5.0;
608			for(i=0; i<bsInfo.nComponents; i++){
609			if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
610			}
611	mmeineke	506	cell *= 1.2; // add a little buffer
612
613	mmeineke	502	cell *= M_SQRT2;
614	mmeineke	498
615	mmeineke	501	siteIndex = 0;
616			for(i=0; i<nCells; i++){
617			for(j=0; j<nCells; j++){
618			for(k=0; k<nCells; k++){
619
620			if( molMap[siteIndex] >= 0 ){
621			pos[0] = i * cell;
622			pos[1] = j * cell;
623			pos[2] = k * cell;
624
625			getRandomRot( rot );
626			molID = molSeq[molMap[siteIndex]];
627	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
628			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
629	mmeineke	501	}
630			siteIndex++;
631	mmeineke	498
632	mmeineke	501	if( molMap[siteIndex] >= 0 ){
633			pos[0] = i * cell + (0.5 * cell);
634			pos[1] = j * cell;
635			pos[2] = k * cell + (0.5 * cell);
636	mmeineke	504
637	mmeineke	501	getRandomRot( rot );
638			molID = molSeq[molMap[siteIndex]];
639	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
640			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
641	mmeineke	501	}
642			siteIndex++;
643	mmeineke	498
644	mmeineke	501	if( molMap[siteIndex] >= 0 ){
645			pos[0] = i * cell + (0.5 * cell);
646			pos[1] = j * cell + (0.5 * cell);
647			pos[2] = k * cell;
648
649			getRandomRot( rot );
650			molID = molSeq[molMap[siteIndex]];
651	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
652			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
653	mmeineke	501	}
654			siteIndex++;
655	mmeineke	498
656	mmeineke	501	if( molMap[siteIndex] >= 0 ){
657			pos[0] = i * cell;
658			pos[1] = j * cell + (0.5 * cell);
659			pos[2] = k * cell + (0.5 * cell);
660	mmeineke	504
661	mmeineke	501	getRandomRot( rot );
662			molID = molSeq[molMap[siteIndex]];
663	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
664			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
665	mmeineke	501	}
666			siteIndex++;
667			}
668			}
669			}
670	mmeineke	498
671	mmeineke	501	// set up the SimInfo object
672
673			bsInfo.boxX = nCells * cell;
674			bsInfo.boxY = nCells * cell;
675			bsInfo.boxZ = nCells * cell;
676
677	mmeineke	573	double boxVector[3];
678	mmeineke	501	simnfo = new SimInfo();
679	mmeineke	502	simnfo->n_atoms = nAtoms;
680	mmeineke	573	boxVector[0] = bsInfo.boxX;
681			boxVector[1] = bsInfo.boxY;
682			boxVector[2] = bsInfo.boxZ;
683			simnfo->setBox( boxVector );
684
685	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
686	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
687	mmeineke	504
688	mmeineke	502	simnfo->atoms = atoms;
689	mmeineke	501
690			// set up the writer and write out
691
692	mmeineke	502	writer = new DumpWriter( simnfo );
693	mmeineke	573	writer->writeFinal(0.0);
694	mmeineke	498
695	mmeineke	501	// clean up the memory
696
697			if( molMap != NULL ) delete[] molMap;
698			if( cardDeck != NULL ) delete[] cardDeck;
699			if( locate != NULL ){
700			for(i=0; i<bsInfo.nComponents; i++){
701			delete locate[i];
702	mmeineke	498	}
703	mmeineke	501	delete[] locate;
704			}
705			if( atoms != NULL ){
706			for(i=0; i<nAtoms; i++){
707			delete atoms[i];
708	mmeineke	498	}
709	mmeineke	501	Atom::destroyArrays();
710			delete[] atoms;
711			}
712			if( molSeq != NULL ) delete[] molSeq;
713			if( simnfo != NULL ) delete simnfo;
714			if( writer != NULL ) delete writer;
715	mmeineke	498
716	mmeineke	501	return 1;
717			}
718	mmeineke	498
719
720	mmeineke	501	void getRandomRot( double rot[3][3] ){
721	mmeineke	498
722	mmeineke	501	double theta, phi, psi;
723			double cosTheta;
724	mmeineke	498
725	mmeineke	501	// select random phi, psi, and cosTheta
726	mmeineke	498
727	mmeineke	501	phi = 2.0 * M_PI * drand48();
728			psi = 2.0 * M_PI * drand48();
729			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
730	mmeineke	498
731	mmeineke	501	theta = acos( cosTheta );
732	mmeineke	498
733	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
734			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
735	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
736	mmeineke	501
737			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
738			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
739			rot[1][2] = sin(theta) * cos(psi);
740
741			rot[2][0] = sin(phi) * sin(theta);
742			rot[2][1] = -cos(phi) * sin(theta);
743			rot[2][2] = cos(theta);
744	mmeineke	498	}
745	mmeineke	501
746	mmeineke	536
747
748			void map( double &x, double &y, double &z,
749			double boxX, double boxY, double boxZ ){
750
751			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
752			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
753
754			if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
755			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
756
757			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
758			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
759			}