trunk/makeLipid/makeRandomLipid.cpp

#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>

#include "SimSetup.hpp"
#include "SimInfo.hpp"
#include "Atom.hpp"
#include "Integrator.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"


void map( double &x, double &y, double &z,
          double boxX, double boxY, double boxZ );

void rotate( double &x, double &y, double &z,
             double theta, double phi, double psi );

char* program_name;
using namespace std;

int main(int argc,char* argv[]){
  
  int i, j, k, l;
  unsigned int n_atoms, eo, xo;
  char* in_name;
  SimSetup* startMe;
  SimInfo* entry_plug;
  Thermo* tStats;
  
  int lipidNAtoms;
  Atom** lipidAtoms;
  
  int tot_Natoms;
  Atom** totAtoms;
  
  const double water_rho = 0.0334; // number density per cubic angstrom
  const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
  const double water_cell = 4.929; // fcc unit cell length

  int n_lipids = 50;
  double water_ratio = 25.0; // water to lipid ratio
  int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 );
  
  std::cerr << "n_lipids = " << n_lipids << "\n";

  double water_shell = 10.0;
  double water_padding = 2.5;
  double lipid_spaceing = 2.5;
  
  srand48( 1337 ); // initialize the random number generator.
 
  program_name = argv[0]; /*save the program name in case we need it*/
  if( argc < 3 ){
    cerr<< "Error, input and output bass files are needed to run.\n"
        << program_name << " <input.bass> <output.bass>\n"; 
    
    exit(8);
  }

  in_name = argv[1];
  char* out_name = argv[2];
  entry_plug = new SimInfo;
  
  startMe = new SimSetup;
  startMe->setSimInfo( entry_plug );
  startMe->parseFile( in_name );
  startMe->createSim();

  delete startMe;

  lipidAtoms = entry_plug->atoms;
  lipidNAtoms = entry_plug->n_atoms;


  // find the width, height, and length of the molecule

  double min_x, min_y, min_z;
  double max_x, max_y, max_z;
  double test_x, test_y, test_z;

  max_x = min_x = lipidAtoms[0]->getX();
  max_y = min_y = lipidAtoms[0]->getY();
  max_z = min_z = lipidAtoms[0]->getZ();
  
  for(i=0; i<lipidNAtoms; i++){

    test_x = lipidAtoms[i]->getX();
    test_y = lipidAtoms[i]->getY();
    test_z = lipidAtoms[i]->getZ();

    if( test_x < min_x ) min_x = test_x;
    if( test_y < min_y ) min_y = test_y;
    if( test_z < min_z ) min_z = test_z;

    if( test_x > max_x ) max_x = test_x;
    if( test_y > max_y ) max_y = test_y;
    if( test_z > max_z ) max_z = test_z;
  }

  double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 )
    + pow((max_x - min_x), 2 );
  double max_length = sqrt( ml2 );

  
  // from this information, create the test box


  double box_x;
  double box_y;
  double box_z;

  box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells

  int n_cellX = (int)(box_x / water_cell + 1.0 );
  int n_cellY = (int)(box_y / water_cell + 1.0 );
  int n_cellZ = (int)(box_z / water_cell + 1.0 );
  
  box_x = water_cell * n_cellX;
  box_y = water_cell * n_cellY;
  box_z = water_cell * n_cellZ;

  int n_water = n_cellX * n_cellY * n_cellZ * 4;
  
  double *waterX = new double[n_water];
  double *waterY = new double[n_water];
  double *waterZ = new double[n_water];
  
  
  // find the center of the test lipid, and make it the center of our
  // soon to be created water box.


  double cx, cy, cz;

  cx = 0.0;
  cy = 0.0;
  cz = 0.0;
  for(i=0; i<lipidNAtoms; i++){
    cx += lipidAtoms[i]->getX();
    cy += lipidAtoms[i]->getY();
    cz += lipidAtoms[i]->getZ();
  }
  cx /= lipidNAtoms;
  cy /= lipidNAtoms;
  cz /= lipidNAtoms;
   
  double x0 = cx - ( box_x * 0.5 );
  double y0 = cy - ( box_y * 0.5 );
  double z0 = cz - ( box_z * 0.5 );


  // create an fcc lattice in the water box.

  
  int ndx = 0;
  for( i=0; i < n_cellX; i++ ){
    for( j=0; j < n_cellY; j++ ){
      for( k=0; k < n_cellZ; k++ ){
        
        waterX[ndx] = i * water_cell + x0; 
        waterY[ndx] = j * water_cell + y0;
        waterZ[ndx] = k * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
        waterZ[ndx] =   k * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + x0;
        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
        waterY[ndx] = j * water_cell + y0;
        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
        ndx++;
      }
    }
  }


  // calculate the number of water's displaced by our molecule.

  int *isActive = new int[n_water];
  for(i=0; i<n_water; 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<n_water) && isActive[i] ); i++){
    for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){

      dx = waterX[i] - lipidAtoms[j]->getX();
      dy = waterY[i] - lipidAtoms[j]->getY();
      dz = waterZ[i] - lipidAtoms[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_deleted++;
      }
    }
  }
  
  n_h2o_target += n_deleted * n_lipids;

  // find a box size that best suits the number of waters we need.

  int done = 0;
  
  if( n_water < n_h2o_target ){
    
    int n_generated = n_cellX;
    int n_test, nx, ny, nz;
    nx = n_cellX;
    ny = n_cellY;
    nz = n_cellZ;
    
    n_test = 4 * nx * ny * nz;

    while( n_test < n_h2o_target ){
      
      nz++;
      n_test = 4 * nx * ny * nz;
    }
    
    int n_diff, goodX, goodY, goodZ;
    
    n_diff = n_test - n_h2o_target;
    goodX = nx;
    goodY = ny;
    goodZ = nz;
    
    int test_diff;
    int n_limit = nz;
    nz = n_cellZ;
    
    for( i=n_generated; i<=n_limit; i++ ){
      for( j=i; j<=n_limit; j++ ){
        for( k=j; k<=n_limit; k++ ){
          
          n_test = 4 * i * j * k;
          
          if( n_test > n_h2o_target ){
            
            test_diff = n_test - n_h2o_target;
            
            if( test_diff < n_diff ){
              
              n_diff = test_diff;
              goodX = nx;
              goodY = ny;
              goodZ = nz;
            }
          }
        }
      }
    }

    n_cellX = goodX;
    n_cellY = goodY;
    n_cellZ = goodZ;
  }

  // we now have the best box size for the simulation. Next we
  // recreate the water box to the new specifications.

  n_water = n_cellX * n_cellY * n_cellZ * 4;
  
  delete[] waterX;
  delete[] waterY;
  delete[] waterZ;
  
  waterX = new double[n_water];
  waterY = new double[n_water];
  waterZ = new double[n_water];

  box_x = water_cell * n_cellX;
  box_y = water_cell * n_cellY;
  box_z = water_cell * n_cellZ;
     
  x0 = 0.0;
  y0 = 0.0;
  z0 = 0.0;

  cx = ( box_x * 0.5 );
  cy = ( box_y * 0.5 );
  cz = ( box_z * 0.5 );

  // create an fcc lattice in the water box.
  
  ndx = 0;
  for( i=0; i < n_cellX; i++ ){
    for( j=0; j < n_cellY; j++ ){
      for( k=0; k < n_cellZ; k++ ){
        
        waterX[ndx] = i * water_cell + x0; 
        waterY[ndx] = j * water_cell + y0;
        waterZ[ndx] = k * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
        waterZ[ndx] =   k * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + x0;
        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
        ndx++;
        
        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
        waterY[ndx] = j * water_cell + y0;
        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
        ndx++;
      }
    }
  }
  
  // **************************************************************  
  

  // start a 3D RSA for the for the lipid placements
  
  srand48( 1337 );

  int rsaNAtoms = n_lipids * lipidNAtoms;
  Atom** rsaAtoms = new Atom*[rsaNAtoms];
          
  DirectionalAtom* dAtom;
  DirectionalAtom* dAtomNew;

  double rotMat[3][3];
  double unitRotMat[3][3];
  
  unitRotMat[0][0] = 1.0;
  unitRotMat[0][1] = 0.0;
  unitRotMat[0][2] = 0.0;
  
  unitRotMat[1][0] = 0.0;
  unitRotMat[1][1] = 1.0;
  unitRotMat[1][2] = 0.0;
  
  unitRotMat[2][0] = 0.0;
  unitRotMat[2][1] = 0.0;
  unitRotMat[2][2] = 1.0;

  ndx = 0;
  for(i=0; i<n_lipids; i++ ){
    for(j=0; j<lipidNAtoms; j++){
      
      if( lipidAtoms[j]->isDirectional() ){
        dAtom = (DirectionalAtom *)lipidAtoms[j];
        
        dAtomNew = new DirectionalAtom();
        dAtomNew->setSUx( dAtom->getSUx() );
        dAtomNew->setSUx( dAtom->getSUx() );
        dAtomNew->setSUx( dAtom->getSUx() );
        
        dAtom->getA( rotMat );
        dAtomNew->setA( rotMat );
        
        rsaAtoms[ndx] = dAtomNew;
      }
      else{
        
        rsaAtoms[ndx] = new GeneralAtom();
      }
      
      rsaAtoms[ndx]->setType( lipidAtoms[j]->getType() );
      
      ndx++;
    }
  }
  
  double testX, testY, testZ;
  double theta, phi, psi;
  double tempX, tempY, tempZ;
  int reject;
  int testDX, acceptedDX;
  
  rCutSqr = lipid_spaceing * lipid_spaceing;
  
  for(i=0; i<n_lipids; i++ ){
    done = 0;
    while( !done ){
      
      testX = drand48() * box_x;
      testY = drand48() * box_y;
      testZ = drand48() * box_z;
      
      theta = drand48() * 2.0 * M_PI;
      phi   = drand48() * 2.0 * M_PI;
      psi   = drand48() * 2.0 * M_PI;
      
      ndx = i * lipidNAtoms;
      for(j=0; j<lipidNAtoms; j++){
       
        tempX = lipidAtoms[j]->getX();
        tempY = lipidAtoms[j]->getY();
        tempZ = lipidAtoms[j]->getZ();

        rotate( tempX, tempY, tempZ, theta, phi, psi );
        
        rsaAtoms[ndx + j]->setX( tempX + testX );
        rsaAtoms[ndx + j]->setY( tempY + testY );
        rsaAtoms[ndx + j]->setZ( tempZ + testZ );
      }
      
      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 = rsaAtoms[testDX]->getX() - rsaAtoms[acceptedDX]->getX();
            dy = rsaAtoms[testDX]->getY() - rsaAtoms[acceptedDX]->getY();
            dz = rsaAtoms[testDX]->getZ() - rsaAtoms[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 ){
        done = 1;
        std::cerr << i << " has been accepted\n";
      }
    }
  }
        
  // cut out the waters that overlap with the lipids.
  
  delete[] isActive;
  isActive = new int[n_water];
  for(i=0; i<n_water; i++) isActive[i] = 1;
  int n_active = n_water;
  rCutSqr = water_padding * water_padding;
  
  for(i=0; ( (i<n_water) && isActive[i] ); i++){
    for(j=0; ( (j<rsaNAtoms) && isActive[i] ); j++){

      dx = waterX[i] - rsaAtoms[j]->getX();
      dy = waterY[i] - rsaAtoms[j]->getY();
      dz = waterZ[i] - rsaAtoms[j]->getZ();

      map( dx, dy, dz, box_x, box_y, box_z );

      dx2 = dx * dx;
      dy2 = dy * dy;
      dz2 = dz * dz;
      
      dSqr = dx2 + dy2 + dz2;
      if( dSqr < rCutSqr ){
        isActive[i] = 0;
        n_active--;
      }
    }
  }

  std::cerr << "final n_waters = " << n_active << "\n";

  // place all of the waters and lipids into one new array

  int new_nAtoms = rsaNAtoms + n_active;
  Atom** new_atoms = new Atom*[new_nAtoms];
  
  ndx = 0;
  for(i=0; i<rsaNAtoms; i++ ){
    
    if( rsaAtoms[i]->isDirectional() ){
      dAtom = (DirectionalAtom *)rsaAtoms[i];
      
      dAtomNew = new DirectionalAtom();
      dAtomNew->setSUx( dAtom->getSUx() );
      dAtomNew->setSUx( dAtom->getSUx() );
      dAtomNew->setSUx( dAtom->getSUx() );
      
      dAtom->getA( rotMat );
      dAtomNew->setA( rotMat );
      
      new_atoms[ndx] = dAtomNew;
    }
    else{
      
      new_atoms[ndx] = new GeneralAtom();
    }

    new_atoms[ndx]->setType( rsaAtoms[i]->getType() );    
    
    new_atoms[ndx]->setX( rsaAtoms[i]->getX() );    
    new_atoms[ndx]->setY( rsaAtoms[i]->getY() );
    new_atoms[ndx]->setZ( rsaAtoms[i]->getZ() );

    new_atoms[ndx]->set_vx( 0.0 );
    new_atoms[ndx]->set_vy( 0.0 );
    new_atoms[ndx]->set_vz( 0.0 );
    
    ndx++;
  }

  for(i=0; i<n_water; i++){
    if(isActive[i]){
      
      new_atoms[ndx] = new DirectionalAtom();
      new_atoms[ndx]->setType( "SSD" );

      new_atoms[ndx]->setX( waterX[i] );    
      new_atoms[ndx]->setY( waterY[i] );
      new_atoms[ndx]->setZ( waterZ[i] );
      
      new_atoms[ndx]->set_vx( 0.0 );
      new_atoms[ndx]->set_vy( 0.0 );
      new_atoms[ndx]->set_vz( 0.0 );
      
      dAtom = (DirectionalAtom *) new_atoms[ndx];

      dAtom->setSUx( 0.0 );
      dAtom->setSUy( 0.0 );
      dAtom->setSUz( 1.0 );
      
      dAtom->setA( unitRotMat );
      
      ndx++;
    }
  }

  entry_plug->n_atoms = new_nAtoms;
  entry_plug->atoms = new_atoms;
  
  entry_plug->box_x = box_x;
  entry_plug->box_y = box_y;
  entry_plug->box_z = box_z;

  DumpWriter* xyz_out = new DumpWriter( entry_plug );
  xyz_out->writeFinal();
  delete xyz_out;

  FILE* out_file;
  
  out_file = fopen( out_name, "w" );
  
  fprintf(out_file,
          "#include \"water.mdl\"\n"
          "#include \"lipid.mdl\"\n"
          "\n"
          "nComponents = 2;\n"
          "component{\n"
          "  type = \"theLipid\";\n"
          "  nMol = %d;\n"
          "}\n"
          "\n"
          "component{\n"
          "  type = \"SSD_water\";\n"
          "  nMol = %d;\n"
          "}\n"
          "\n"
          "initialConfig = \"%s\";\n"
          "\n"
          "boxX = %lf;\n"
          "boxY = %lf;\n"
          "boxZ = %lf;\n",
          n_lipids, n_active, entry_plug->finalName,
          box_x, box_y, box_z );
 
  fclose( out_file );
  
  return 0;
}


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


void rotate( double &x, double &y, double &z,
             double theta, double phi, double psi ){

  double newX, newY, newZ;

  double A[3][3];

  A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
  A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
  A[0][2] = sin(theta) * sin(psi);
  
  A[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  A[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  A[1][2] = sin(theta) * cos(psi);

  A[2][0] = sin(phi) * sin(theta);
  A[2][1] = -cos(phi) * sin(theta);
  A[2][2] = cos(theta);

  newX = (x * A[0][0]) + (y * A[0][1]) + (z * A[0][2]);
  newY = (x * A[1][0]) + (y * A[1][1]) + (z * A[1][2]);
  newZ = (x * A[2][0]) + (y * A[2][1]) + (z * A[2][2]);

  x = newX;
  y = newY;
  z = newZ;
}
Revision:	33
Committed:	Wed Jul 17 20:36:25 2002 UTC (21 years, 11 months ago) by mmeineke
File size:	14679 byte(s)
Log Message:	cleaned up the output a little
#	User	Rev	Content
1	mmeineke	28	#include <iostream>
2			#include <cstdlib>
3			#include <cstring>
4			#include <cstdio>
5	mmeineke	29	#include <cmath>
6	mmeineke	28
7			#include "SimSetup.hpp"
8			#include "SimInfo.hpp"
9			#include "Atom.hpp"
10			#include "Integrator.hpp"
11			#include "Thermo.hpp"
12			#include "ReadWrite.hpp"
13
14
15	mmeineke	32	void map( double &x, double &y, double &z,
16			double boxX, double boxY, double boxZ );
17	mmeineke	29
18	mmeineke	32	void rotate( double &x, double &y, double &z,
19			double theta, double phi, double psi );
20
21	mmeineke	28	char* program_name;
22			using namespace std;
23
24			int main(int argc,char* argv[]){
25
26			int i, j, k, l;
27			unsigned int n_atoms, eo, xo;
28			char* in_name;
29			SimSetup* startMe;
30			SimInfo* entry_plug;
31			Thermo* tStats;
32
33			int lipidNAtoms;
34			Atom** lipidAtoms;
35
36			int tot_Natoms;
37			Atom** totAtoms;
38
39			const double water_rho = 0.0334; // number density per cubic angstrom
40			const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
41			const double water_cell = 4.929; // fcc unit cell length
42
43	mmeineke	29	int n_lipids = 50;
44			double water_ratio = 25.0; // water to lipid ratio
45			int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 );
46	mmeineke	28
47			std::cerr << "n_lipids = " << n_lipids << "\n";
48
49			double water_shell = 10.0;
50			double water_padding = 2.5;
51	mmeineke	32	double lipid_spaceing = 2.5;
52	mmeineke	28
53			srand48( 1337 ); // initialize the random number generator.
54
55			program_name = argv[0]; /save the program name in case we need it/
56			if( argc < 3 ){
57			cerr<< "Error, input and output bass files are needed to run.\n"
58			<< program_name << " <input.bass> <output.bass>\n";
59
60			exit(8);
61			}
62
63			in_name = argv[1];
64			char* out_name = argv[2];
65			entry_plug = new SimInfo;
66
67			startMe = new SimSetup;
68			startMe->setSimInfo( entry_plug );
69			startMe->parseFile( in_name );
70			startMe->createSim();
71
72			delete startMe;
73
74			lipidAtoms = entry_plug->atoms;
75			lipidNAtoms = entry_plug->n_atoms;
76
77
78	mmeineke	29	// find the width, height, and length of the molecule
79	mmeineke	28
80			double min_x, min_y, min_z;
81			double max_x, max_y, max_z;
82			double test_x, test_y, test_z;
83
84	mmeineke	29	max_x = min_x = lipidAtoms[0]->getX();
85			max_y = min_y = lipidAtoms[0]->getY();
86			max_z = min_z = lipidAtoms[0]->getZ();
87	mmeineke	28
88	mmeineke	29	for(i=0; i<lipidNAtoms; i++){
89	mmeineke	28
90	mmeineke	29	test_x = lipidAtoms[i]->getX();
91			test_y = lipidAtoms[i]->getY();
92			test_z = lipidAtoms[i]->getZ();
93	mmeineke	28
94			if( test_x < min_x ) min_x = test_x;
95			if( test_y < min_y ) min_y = test_y;
96			if( test_z < min_z ) min_z = test_z;
97
98			if( test_x > max_x ) max_x = test_x;
99			if( test_y > max_y ) max_y = test_y;
100			if( test_z > max_z ) max_z = test_z;
101			}
102
103	mmeineke	29	double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 )
104			+ pow((max_x - min_x), 2 );
105			double max_length = sqrt( ml2 );
106	mmeineke	28
107
108	mmeineke	29	// from this information, create the test box
109
110
111			double box_x;
112			double box_y;
113			double box_z;
114
115			box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells
116
117			int n_cellX = (int)(box_x / water_cell + 1.0 );
118			int n_cellY = (int)(box_y / water_cell + 1.0 );
119			int n_cellZ = (int)(box_z / water_cell + 1.0 );
120
121	mmeineke	28	box_x = water_cell * n_cellX;
122			box_y = water_cell * n_cellY;
123			box_z = water_cell * n_cellZ;
124
125			int n_water = n_cellX * n_cellY * n_cellZ * 4;
126
127			double *waterX = new double[n_water];
128			double *waterY = new double[n_water];
129			double *waterZ = new double[n_water];
130
131	mmeineke	29
132			// find the center of the test lipid, and make it the center of our
133			// soon to be created water box.
134
135
136	mmeineke	28	double cx, cy, cz;
137
138			cx = 0.0;
139			cy = 0.0;
140			cz = 0.0;
141	mmeineke	29	for(i=0; i<lipidNAtoms; i++){
142			cx += lipidAtoms[i]->getX();
143			cy += lipidAtoms[i]->getY();
144			cz += lipidAtoms[i]->getZ();
145	mmeineke	28	}
146	mmeineke	29	cx /= lipidNAtoms;
147			cy /= lipidNAtoms;
148			cz /= lipidNAtoms;
149	mmeineke	28
150			double x0 = cx - ( box_x * 0.5 );
151			double y0 = cy - ( box_y * 0.5 );
152			double z0 = cz - ( box_z * 0.5 );
153	mmeineke	29
154
155			// create an fcc lattice in the water box.
156
157	mmeineke	28
158	mmeineke	32	int ndx = 0;
159	mmeineke	28	for( i=0; i < n_cellX; i++ ){
160			for( j=0; j < n_cellY; j++ ){
161			for( k=0; k < n_cellZ; k++ ){
162
163	mmeineke	32	waterX[ndx] = i * water_cell + x0;
164			waterY[ndx] = j * water_cell + y0;
165			waterZ[ndx] = k * water_cell + z0;
166			ndx++;
167	mmeineke	28
168	mmeineke	32	waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
169			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
170			waterZ[ndx] = k * water_cell + z0;
171			ndx++;
172	mmeineke	28
173	mmeineke	32	waterX[ndx] = i * water_cell + x0;
174			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
175			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
176			ndx++;
177	mmeineke	28
178	mmeineke	32	waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
179			waterY[ndx] = j * water_cell + y0;
180			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
181			ndx++;
182	mmeineke	28	}
183			}
184			}
185
186	mmeineke	29
187			// calculate the number of water's displaced by our molecule.
188
189	mmeineke	28	int *isActive = new int[n_water];
190			for(i=0; i<n_water; i++) isActive[i] = 1;
191
192	mmeineke	29	int n_deleted = 0;
193	mmeineke	28	double dx, dy, dz;
194			double dx2, dy2, dz2, dSqr;
195			double rCutSqr = water_padding * water_padding;
196
197			for(i=0; ( (i<n_water) && isActive[i] ); i++){
198	mmeineke	29	for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){
199	mmeineke	28
200	mmeineke	29	dx = waterX[i] - lipidAtoms[j]->getX();
201			dy = waterY[i] - lipidAtoms[j]->getY();
202			dz = waterZ[i] - lipidAtoms[j]->getZ();
203	mmeineke	28
204	mmeineke	32	map( dx, dy, dz, box_x, box_y, box_z );
205	mmeineke	29
206	mmeineke	28	dx2 = dx * dx;
207			dy2 = dy * dy;
208			dz2 = dz * dz;
209
210			dSqr = dx2 + dy2 + dz2;
211			if( dSqr < rCutSqr ){
212			isActive[i] = 0;
213	mmeineke	29	n_deleted++;
214	mmeineke	28	}
215			}
216			}
217	mmeineke	32
218	mmeineke	29	n_h2o_target += n_deleted * n_lipids;
219	mmeineke	28
220	mmeineke	29	// find a box size that best suits the number of waters we need.
221
222			int done = 0;
223
224	mmeineke	32	if( n_water < n_h2o_target ){
225	mmeineke	29
226			int n_generated = n_cellX;
227			int n_test, nx, ny, nz;
228			nx = n_cellX;
229			ny = n_cellY;
230			nz = n_cellZ;
231
232	mmeineke	32	n_test = 4 * nx * ny * nz;
233
234	mmeineke	29	while( n_test < n_h2o_target ){
235
236			nz++;
237			n_test = 4 * nx * ny * nz;
238			}
239
240			int n_diff, goodX, goodY, goodZ;
241
242	mmeineke	32	n_diff = n_test - n_h2o_target;
243	mmeineke	29	goodX = nx;
244			goodY = ny;
245			goodZ = nz;
246
247			int test_diff;
248	mmeineke	32	int n_limit = nz;
249			nz = n_cellZ;
250	mmeineke	29
251			for( i=n_generated; i<=n_limit; i++ ){
252			for( j=i; j<=n_limit; j++ ){
253			for( k=j; k<=n_limit; k++ ){
254
255			n_test = 4 * i * j * k;
256
257			if( n_test > n_h2o_target ){
258
259			test_diff = n_test - n_h2o_target;
260
261			if( test_diff < n_diff ){
262
263			n_diff = test_diff;
264			goodX = nx;
265			goodY = ny;
266			goodZ = nz;
267			}
268			}
269			}
270			}
271			}
272
273			n_cellX = goodX;
274			n_cellY = goodY;
275			n_cellZ = goodZ;
276			}
277
278	mmeineke	32	// we now have the best box size for the simulation. Next we
279			// recreate the water box to the new specifications.
280
281			n_water = n_cellX * n_cellY * n_cellZ * 4;
282
283			delete[] waterX;
284			delete[] waterY;
285			delete[] waterZ;
286
287			waterX = new double[n_water];
288			waterY = new double[n_water];
289			waterZ = new double[n_water];
290
291			box_x = water_cell * n_cellX;
292			box_y = water_cell * n_cellY;
293			box_z = water_cell * n_cellZ;
294
295			x0 = 0.0;
296			y0 = 0.0;
297			z0 = 0.0;
298
299			cx = ( box_x * 0.5 );
300			cy = ( box_y * 0.5 );
301			cz = ( box_z * 0.5 );
302
303			// create an fcc lattice in the water box.
304
305			ndx = 0;
306			for( i=0; i < n_cellX; i++ ){
307			for( j=0; j < n_cellY; j++ ){
308			for( k=0; k < n_cellZ; k++ ){
309
310			waterX[ndx] = i * water_cell + x0;
311			waterY[ndx] = j * water_cell + y0;
312			waterZ[ndx] = k * water_cell + z0;
313			ndx++;
314
315			waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
316			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
317			waterZ[ndx] = k * water_cell + z0;
318			ndx++;
319
320			waterX[ndx] = i * water_cell + x0;
321			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
322			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
323			ndx++;
324
325			waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
326			waterY[ndx] = j * water_cell + y0;
327			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
328			ndx++;
329			}
330			}
331			}
332
333			// **************************************************************
334
335
336
337			// start a 3D RSA for the for the lipid placements
338
339			srand48( 1337 );
340
341			int rsaNAtoms = n_lipids * lipidNAtoms;
342			Atom** rsaAtoms = new Atom*[rsaNAtoms];
343	mmeineke	29
344	mmeineke	32	DirectionalAtom* dAtom;
345			DirectionalAtom* dAtomNew;
346	mmeineke	29
347	mmeineke	32	double rotMat[3][3];
348			double unitRotMat[3][3];
349
350			unitRotMat[0][0] = 1.0;
351			unitRotMat[0][1] = 0.0;
352			unitRotMat[0][2] = 0.0;
353
354			unitRotMat[1][0] = 0.0;
355			unitRotMat[1][1] = 1.0;
356			unitRotMat[1][2] = 0.0;
357
358			unitRotMat[2][0] = 0.0;
359			unitRotMat[2][1] = 0.0;
360			unitRotMat[2][2] = 1.0;
361	mmeineke	29
362	mmeineke	32	ndx = 0;
363			for(i=0; i<n_lipids; i++ ){
364			for(j=0; j<lipidNAtoms; j++){
365
366			if( lipidAtoms[j]->isDirectional() ){
367			dAtom = (DirectionalAtom *)lipidAtoms[j];
368
369			dAtomNew = new DirectionalAtom();
370			dAtomNew->setSUx( dAtom->getSUx() );
371			dAtomNew->setSUx( dAtom->getSUx() );
372			dAtomNew->setSUx( dAtom->getSUx() );
373
374			dAtom->getA( rotMat );
375			dAtomNew->setA( rotMat );
376
377			rsaAtoms[ndx] = dAtomNew;
378			}
379			else{
380
381			rsaAtoms[ndx] = new GeneralAtom();
382			}
383
384			rsaAtoms[ndx]->setType( lipidAtoms[j]->getType() );
385
386			ndx++;
387			}
388			}
389
390			double testX, testY, testZ;
391			double theta, phi, psi;
392			double tempX, tempY, tempZ;
393			int reject;
394			int testDX, acceptedDX;
395
396			rCutSqr = lipid_spaceing * lipid_spaceing;
397
398			for(i=0; i<n_lipids; i++ ){
399			done = 0;
400			while( !done ){
401
402			testX = drand48() * box_x;
403			testY = drand48() * box_y;
404			testZ = drand48() * box_z;
405
406			theta = drand48() * 2.0 * M_PI;
407			phi = drand48() * 2.0 * M_PI;
408			psi = drand48() * 2.0 * M_PI;
409
410			ndx = i * lipidNAtoms;
411			for(j=0; j<lipidNAtoms; j++){
412
413			tempX = lipidAtoms[j]->getX();
414			tempY = lipidAtoms[j]->getY();
415			tempZ = lipidAtoms[j]->getZ();
416	mmeineke	29
417	mmeineke	32	rotate( tempX, tempY, tempZ, theta, phi, psi );
418
419			rsaAtoms[ndx + j]->setX( tempX + testX );
420			rsaAtoms[ndx + j]->setY( tempY + testY );
421			rsaAtoms[ndx + j]->setZ( tempZ + testZ );
422			}
423
424			reject = 0;
425			for( j=0; !reject && j<i; j++){
426			for(k=0; !reject && k<lipidNAtoms; k++){
427
428			acceptedDX = j*lipidNAtoms + k;
429			for(l=0; !reject && l<lipidNAtoms; l++){
430
431			testDX = ndx + l;
432	mmeineke	29
433	mmeineke	32	dx = rsaAtoms[testDX]->getX() - rsaAtoms[acceptedDX]->getX();
434			dy = rsaAtoms[testDX]->getY() - rsaAtoms[acceptedDX]->getY();
435			dz = rsaAtoms[testDX]->getZ() - rsaAtoms[acceptedDX]->getZ();
436
437			map( dx, dy, dz, box_x, box_y, box_z );
438
439			dx2 = dx * dx;
440			dy2 = dy * dy;
441			dz2 = dz * dz;
442
443			dSqr = dx2 + dy2 + dz2;
444			if( dSqr < rCutSqr ) reject = 1;
445			}
446			}
447			}
448	mmeineke	29
449	mmeineke	32	if( !reject ){
450			done = 1;
451			std::cerr << i << " has been accepted\n";
452			}
453			}
454			}
455
456			// cut out the waters that overlap with the lipids.
457
458			delete[] isActive;
459			isActive = new int[n_water];
460			for(i=0; i<n_water; i++) isActive[i] = 1;
461			int n_active = n_water;
462			rCutSqr = water_padding * water_padding;
463
464			for(i=0; ( (i<n_water) && isActive[i] ); i++){
465			for(j=0; ( (j<rsaNAtoms) && isActive[i] ); j++){
466	mmeineke	29
467	mmeineke	32	dx = waterX[i] - rsaAtoms[j]->getX();
468			dy = waterY[i] - rsaAtoms[j]->getY();
469			dz = waterZ[i] - rsaAtoms[j]->getZ();
470	mmeineke	29
471	mmeineke	32	map( dx, dy, dz, box_x, box_y, box_z );
472
473			dx2 = dx * dx;
474			dy2 = dy * dy;
475			dz2 = dz * dz;
476
477			dSqr = dx2 + dy2 + dz2;
478			if( dSqr < rCutSqr ){
479			isActive[i] = 0;
480			n_active--;
481			}
482			}
483			}
484
485			std::cerr << "final n_waters = " << n_active << "\n";
486
487			// place all of the waters and lipids into one new array
488
489			int new_nAtoms = rsaNAtoms + n_active;
490	mmeineke	28	Atom** new_atoms = new Atom*[new_nAtoms];
491
492	mmeineke	32	ndx = 0;
493			for(i=0; i<rsaNAtoms; i++ ){
494	mmeineke	28
495	mmeineke	32	if( rsaAtoms[i]->isDirectional() ){
496			dAtom = (DirectionalAtom *)rsaAtoms[i];
497	mmeineke	28
498			dAtomNew = new DirectionalAtom();
499			dAtomNew->setSUx( dAtom->getSUx() );
500			dAtomNew->setSUx( dAtom->getSUx() );
501			dAtomNew->setSUx( dAtom->getSUx() );
502
503	mmeineke	32	dAtom->getA( rotMat );
504	mmeineke	28	dAtomNew->setA( rotMat );
505
506	mmeineke	32	new_atoms[ndx] = dAtomNew;
507	mmeineke	28	}
508			else{
509
510	mmeineke	32	new_atoms[ndx] = new GeneralAtom();
511	mmeineke	28	}
512
513	mmeineke	32	new_atoms[ndx]->setType( rsaAtoms[i]->getType() );
514	mmeineke	28
515	mmeineke	32	new_atoms[ndx]->setX( rsaAtoms[i]->getX() );
516			new_atoms[ndx]->setY( rsaAtoms[i]->getY() );
517			new_atoms[ndx]->setZ( rsaAtoms[i]->getZ() );
518	mmeineke	28
519	mmeineke	32	new_atoms[ndx]->set_vx( 0.0 );
520			new_atoms[ndx]->set_vy( 0.0 );
521			new_atoms[ndx]->set_vz( 0.0 );
522	mmeineke	28
523	mmeineke	32	ndx++;
524	mmeineke	28	}
525
526			for(i=0; i<n_water; i++){
527			if(isActive[i]){
528
529	mmeineke	32	new_atoms[ndx] = new DirectionalAtom();
530			new_atoms[ndx]->setType( "SSD" );
531	mmeineke	28
532	mmeineke	32	new_atoms[ndx]->setX( waterX[i] );
533			new_atoms[ndx]->setY( waterY[i] );
534			new_atoms[ndx]->setZ( waterZ[i] );
535	mmeineke	28
536	mmeineke	32	new_atoms[ndx]->set_vx( 0.0 );
537			new_atoms[ndx]->set_vy( 0.0 );
538			new_atoms[ndx]->set_vz( 0.0 );
539	mmeineke	28
540	mmeineke	32	dAtom = (DirectionalAtom *) new_atoms[ndx];
541	mmeineke	28
542			dAtom->setSUx( 0.0 );
543			dAtom->setSUy( 0.0 );
544			dAtom->setSUz( 1.0 );
545
546	mmeineke	32	dAtom->setA( unitRotMat );
547	mmeineke	28
548	mmeineke	32	ndx++;
549	mmeineke	28	}
550			}
551
552			entry_plug->n_atoms = new_nAtoms;
553			entry_plug->atoms = new_atoms;
554
555			entry_plug->box_x = box_x;
556			entry_plug->box_y = box_y;
557			entry_plug->box_z = box_z;
558
559			DumpWriter* xyz_out = new DumpWriter( entry_plug );
560			xyz_out->writeFinal();
561			delete xyz_out;
562
563			FILE* out_file;
564
565			out_file = fopen( out_name, "w" );
566
567			fprintf(out_file,
568			"#include \"water.mdl\"\n"
569			"#include \"lipid.mdl\"\n"
570			"\n"
571			"nComponents = 2;\n"
572			"component{\n"
573			" type = \"theLipid\";\n"
574			" nMol = %d;\n"
575			"}\n"
576			"\n"
577			"component{\n"
578			" type = \"SSD_water\";\n"
579			" nMol = %d;\n"
580			"}\n"
581			"\n"
582			"initialConfig = \"%s\";\n"
583			"\n"
584			"boxX = %lf;\n"
585			"boxY = %lf;\n"
586			"boxZ = %lf;\n",
587			n_lipids, n_active, entry_plug->finalName,
588			box_x, box_y, box_z );
589
590			fclose( out_file );
591
592			return 0;
593			}
594	mmeineke	29
595
596	mmeineke	32	void map( double &x, double &y, double &z,
597			double boxX, double boxY, double boxZ ){
598
599			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
600			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
601	mmeineke	29
602	mmeineke	32	if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
603			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
604
605			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
606			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
607			}
608	mmeineke	29
609
610	mmeineke	32	void rotate( double &x, double &y, double &z,
611			double theta, double phi, double psi ){
612
613			double newX, newY, newZ;
614
615			double A[3][3];
616
617			A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
618			A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
619			A[0][2] = sin(theta) * sin(psi);
620	mmeineke	29
621	mmeineke	32	A[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
622			A[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
623			A[1][2] = sin(theta) * cos(psi);
624
625			A[2][0] = sin(phi) * sin(theta);
626			A[2][1] = -cos(phi) * sin(theta);
627			A[2][2] = cos(theta);
628
629			newX = (x * A[0][0]) + (y * A[0][1]) + (z * A[0][2]);
630			newY = (x * A[1][0]) + (y * A[1][1]) + (z * A[1][2]);
631			newZ = (x * A[2][0]) + (y * A[2][1]) + (z * A[2][2]);
632
633			x = newX;
634			y = newY;
635			z = newZ;
636	mmeineke	29	}