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++;
      }
    }
  }
  
  std::cerr << "nTarget before: " << n_h2o_target;
  
  n_h2o_target += n_deleted * n_lipids;

  std::cerr << ", after: " << n_h2o_target << ", n_deleted: " << n_deleted
            << "\n";

  // 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;
  
  std::cerr << "new waters = " << n_water << "\n";

  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:	32
Committed:	Wed Jul 17 20:33:56 2002 UTC (23 years, 3 months ago) by mmeineke
File size:	14875 byte(s)
Log Message:	fixed makeRandom to set up a correct simulation
#	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			std::cerr << "nTarget before: " << n_h2o_target;
219
220	mmeineke	29	n_h2o_target += n_deleted * n_lipids;
221	mmeineke	28
222	mmeineke	32	std::cerr << ", after: " << n_h2o_target << ", n_deleted: " << n_deleted
223			<< "\n";
224
225	mmeineke	29	// find a box size that best suits the number of waters we need.
226
227			int done = 0;
228
229	mmeineke	32	if( n_water < n_h2o_target ){
230	mmeineke	29
231			int n_generated = n_cellX;
232			int n_test, nx, ny, nz;
233			nx = n_cellX;
234			ny = n_cellY;
235			nz = n_cellZ;
236
237	mmeineke	32	n_test = 4 * nx * ny * nz;
238
239	mmeineke	29	while( n_test < n_h2o_target ){
240
241			nz++;
242			n_test = 4 * nx * ny * nz;
243			}
244
245			int n_diff, goodX, goodY, goodZ;
246
247	mmeineke	32	n_diff = n_test - n_h2o_target;
248	mmeineke	29	goodX = nx;
249			goodY = ny;
250			goodZ = nz;
251
252			int test_diff;
253	mmeineke	32	int n_limit = nz;
254			nz = n_cellZ;
255	mmeineke	29
256			for( i=n_generated; i<=n_limit; i++ ){
257			for( j=i; j<=n_limit; j++ ){
258			for( k=j; k<=n_limit; k++ ){
259
260			n_test = 4 * i * j * k;
261
262			if( n_test > n_h2o_target ){
263
264			test_diff = n_test - n_h2o_target;
265
266			if( test_diff < n_diff ){
267
268			n_diff = test_diff;
269			goodX = nx;
270			goodY = ny;
271			goodZ = nz;
272			}
273			}
274			}
275			}
276			}
277
278			n_cellX = goodX;
279			n_cellY = goodY;
280			n_cellZ = goodZ;
281			}
282
283	mmeineke	32	// we now have the best box size for the simulation. Next we
284			// recreate the water box to the new specifications.
285
286			n_water = n_cellX * n_cellY * n_cellZ * 4;
287
288			std::cerr << "new waters = " << n_water << "\n";
289
290			delete[] waterX;
291			delete[] waterY;
292			delete[] waterZ;
293
294			waterX = new double[n_water];
295			waterY = new double[n_water];
296			waterZ = new double[n_water];
297
298			box_x = water_cell * n_cellX;
299			box_y = water_cell * n_cellY;
300			box_z = water_cell * n_cellZ;
301
302			x0 = 0.0;
303			y0 = 0.0;
304			z0 = 0.0;
305
306			cx = ( box_x * 0.5 );
307			cy = ( box_y * 0.5 );
308			cz = ( box_z * 0.5 );
309
310			// create an fcc lattice in the water box.
311
312			ndx = 0;
313			for( i=0; i < n_cellX; i++ ){
314			for( j=0; j < n_cellY; j++ ){
315			for( k=0; k < n_cellZ; k++ ){
316
317			waterX[ndx] = i * water_cell + x0;
318			waterY[ndx] = j * water_cell + y0;
319			waterZ[ndx] = k * water_cell + z0;
320			ndx++;
321
322			waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
323			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
324			waterZ[ndx] = k * water_cell + z0;
325			ndx++;
326
327			waterX[ndx] = i * water_cell + x0;
328			waterY[ndx] = j * water_cell + 0.5 * water_cell + y0;
329			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
330			ndx++;
331
332			waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
333			waterY[ndx] = j * water_cell + y0;
334			waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
335			ndx++;
336			}
337			}
338			}
339
340			// **************************************************************
341
342
343
344			// start a 3D RSA for the for the lipid placements
345
346			srand48( 1337 );
347
348			int rsaNAtoms = n_lipids * lipidNAtoms;
349			Atom** rsaAtoms = new Atom*[rsaNAtoms];
350	mmeineke	29
351	mmeineke	32	DirectionalAtom* dAtom;
352			DirectionalAtom* dAtomNew;
353	mmeineke	29
354	mmeineke	32	double rotMat[3][3];
355			double unitRotMat[3][3];
356
357			unitRotMat[0][0] = 1.0;
358			unitRotMat[0][1] = 0.0;
359			unitRotMat[0][2] = 0.0;
360
361			unitRotMat[1][0] = 0.0;
362			unitRotMat[1][1] = 1.0;
363			unitRotMat[1][2] = 0.0;
364
365			unitRotMat[2][0] = 0.0;
366			unitRotMat[2][1] = 0.0;
367			unitRotMat[2][2] = 1.0;
368	mmeineke	29
369	mmeineke	32	ndx = 0;
370			for(i=0; i<n_lipids; i++ ){
371			for(j=0; j<lipidNAtoms; j++){
372
373			if( lipidAtoms[j]->isDirectional() ){
374			dAtom = (DirectionalAtom *)lipidAtoms[j];
375
376			dAtomNew = new DirectionalAtom();
377			dAtomNew->setSUx( dAtom->getSUx() );
378			dAtomNew->setSUx( dAtom->getSUx() );
379			dAtomNew->setSUx( dAtom->getSUx() );
380
381			dAtom->getA( rotMat );
382			dAtomNew->setA( rotMat );
383
384			rsaAtoms[ndx] = dAtomNew;
385			}
386			else{
387
388			rsaAtoms[ndx] = new GeneralAtom();
389			}
390
391			rsaAtoms[ndx]->setType( lipidAtoms[j]->getType() );
392
393			ndx++;
394			}
395			}
396
397			double testX, testY, testZ;
398			double theta, phi, psi;
399			double tempX, tempY, tempZ;
400			int reject;
401			int testDX, acceptedDX;
402
403			rCutSqr = lipid_spaceing * lipid_spaceing;
404
405			for(i=0; i<n_lipids; i++ ){
406			done = 0;
407			while( !done ){
408
409			testX = drand48() * box_x;
410			testY = drand48() * box_y;
411			testZ = drand48() * box_z;
412
413			theta = drand48() * 2.0 * M_PI;
414			phi = drand48() * 2.0 * M_PI;
415			psi = drand48() * 2.0 * M_PI;
416
417			ndx = i * lipidNAtoms;
418			for(j=0; j<lipidNAtoms; j++){
419
420			tempX = lipidAtoms[j]->getX();
421			tempY = lipidAtoms[j]->getY();
422			tempZ = lipidAtoms[j]->getZ();
423	mmeineke	29
424	mmeineke	32	rotate( tempX, tempY, tempZ, theta, phi, psi );
425
426			rsaAtoms[ndx + j]->setX( tempX + testX );
427			rsaAtoms[ndx + j]->setY( tempY + testY );
428			rsaAtoms[ndx + j]->setZ( tempZ + testZ );
429			}
430
431			reject = 0;
432			for( j=0; !reject && j<i; j++){
433			for(k=0; !reject && k<lipidNAtoms; k++){
434
435			acceptedDX = j*lipidNAtoms + k;
436			for(l=0; !reject && l<lipidNAtoms; l++){
437
438			testDX = ndx + l;
439	mmeineke	29
440	mmeineke	32	dx = rsaAtoms[testDX]->getX() - rsaAtoms[acceptedDX]->getX();
441			dy = rsaAtoms[testDX]->getY() - rsaAtoms[acceptedDX]->getY();
442			dz = rsaAtoms[testDX]->getZ() - rsaAtoms[acceptedDX]->getZ();
443
444			map( dx, dy, dz, box_x, box_y, box_z );
445
446			dx2 = dx * dx;
447			dy2 = dy * dy;
448			dz2 = dz * dz;
449
450			dSqr = dx2 + dy2 + dz2;
451			if( dSqr < rCutSqr ) reject = 1;
452			}
453			}
454			}
455	mmeineke	29
456	mmeineke	32	if( !reject ){
457			done = 1;
458			std::cerr << i << " has been accepted\n";
459			}
460			}
461			}
462
463			// cut out the waters that overlap with the lipids.
464
465			delete[] isActive;
466			isActive = new int[n_water];
467			for(i=0; i<n_water; i++) isActive[i] = 1;
468			int n_active = n_water;
469			rCutSqr = water_padding * water_padding;
470
471			for(i=0; ( (i<n_water) && isActive[i] ); i++){
472			for(j=0; ( (j<rsaNAtoms) && isActive[i] ); j++){
473	mmeineke	29
474	mmeineke	32	dx = waterX[i] - rsaAtoms[j]->getX();
475			dy = waterY[i] - rsaAtoms[j]->getY();
476			dz = waterZ[i] - rsaAtoms[j]->getZ();
477	mmeineke	29
478	mmeineke	32	map( dx, dy, dz, box_x, box_y, box_z );
479
480			dx2 = dx * dx;
481			dy2 = dy * dy;
482			dz2 = dz * dz;
483
484			dSqr = dx2 + dy2 + dz2;
485			if( dSqr < rCutSqr ){
486			isActive[i] = 0;
487			n_active--;
488			}
489			}
490			}
491
492			std::cerr << "final n_waters = " << n_active << "\n";
493
494			// place all of the waters and lipids into one new array
495
496			int new_nAtoms = rsaNAtoms + n_active;
497	mmeineke	28	Atom** new_atoms = new Atom*[new_nAtoms];
498
499	mmeineke	32	ndx = 0;
500			for(i=0; i<rsaNAtoms; i++ ){
501	mmeineke	28
502	mmeineke	32	if( rsaAtoms[i]->isDirectional() ){
503			dAtom = (DirectionalAtom *)rsaAtoms[i];
504	mmeineke	28
505			dAtomNew = new DirectionalAtom();
506			dAtomNew->setSUx( dAtom->getSUx() );
507			dAtomNew->setSUx( dAtom->getSUx() );
508			dAtomNew->setSUx( dAtom->getSUx() );
509
510	mmeineke	32	dAtom->getA( rotMat );
511	mmeineke	28	dAtomNew->setA( rotMat );
512
513	mmeineke	32	new_atoms[ndx] = dAtomNew;
514	mmeineke	28	}
515			else{
516
517	mmeineke	32	new_atoms[ndx] = new GeneralAtom();
518	mmeineke	28	}
519
520	mmeineke	32	new_atoms[ndx]->setType( rsaAtoms[i]->getType() );
521	mmeineke	28
522	mmeineke	32	new_atoms[ndx]->setX( rsaAtoms[i]->getX() );
523			new_atoms[ndx]->setY( rsaAtoms[i]->getY() );
524			new_atoms[ndx]->setZ( rsaAtoms[i]->getZ() );
525	mmeineke	28
526	mmeineke	32	new_atoms[ndx]->set_vx( 0.0 );
527			new_atoms[ndx]->set_vy( 0.0 );
528			new_atoms[ndx]->set_vz( 0.0 );
529	mmeineke	28
530	mmeineke	32	ndx++;
531	mmeineke	28	}
532
533			for(i=0; i<n_water; i++){
534			if(isActive[i]){
535
536	mmeineke	32	new_atoms[ndx] = new DirectionalAtom();
537			new_atoms[ndx]->setType( "SSD" );
538	mmeineke	28
539	mmeineke	32	new_atoms[ndx]->setX( waterX[i] );
540			new_atoms[ndx]->setY( waterY[i] );
541			new_atoms[ndx]->setZ( waterZ[i] );
542	mmeineke	28
543	mmeineke	32	new_atoms[ndx]->set_vx( 0.0 );
544			new_atoms[ndx]->set_vy( 0.0 );
545			new_atoms[ndx]->set_vz( 0.0 );
546	mmeineke	28
547	mmeineke	32	dAtom = (DirectionalAtom *) new_atoms[ndx];
548	mmeineke	28
549			dAtom->setSUx( 0.0 );
550			dAtom->setSUy( 0.0 );
551			dAtom->setSUz( 1.0 );
552
553	mmeineke	32	dAtom->setA( unitRotMat );
554	mmeineke	28
555	mmeineke	32	ndx++;
556	mmeineke	28	}
557			}
558
559			entry_plug->n_atoms = new_nAtoms;
560			entry_plug->atoms = new_atoms;
561
562			entry_plug->box_x = box_x;
563			entry_plug->box_y = box_y;
564			entry_plug->box_z = box_z;
565
566			DumpWriter* xyz_out = new DumpWriter( entry_plug );
567			xyz_out->writeFinal();
568			delete xyz_out;
569
570			FILE* out_file;
571
572			out_file = fopen( out_name, "w" );
573
574			fprintf(out_file,
575			"#include \"water.mdl\"\n"
576			"#include \"lipid.mdl\"\n"
577			"\n"
578			"nComponents = 2;\n"
579			"component{\n"
580			" type = \"theLipid\";\n"
581			" nMol = %d;\n"
582			"}\n"
583			"\n"
584			"component{\n"
585			" type = \"SSD_water\";\n"
586			" nMol = %d;\n"
587			"}\n"
588			"\n"
589			"initialConfig = \"%s\";\n"
590			"\n"
591			"boxX = %lf;\n"
592			"boxY = %lf;\n"
593			"boxZ = %lf;\n",
594			n_lipids, n_active, entry_plug->finalName,
595			box_x, box_y, box_z );
596
597			fclose( out_file );
598
599			return 0;
600			}
601	mmeineke	29
602
603	mmeineke	32	void map( double &x, double &y, double &z,
604			double boxX, double boxY, double boxZ ){
605
606			if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
607			else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
608	mmeineke	29
609	mmeineke	32	if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
610			else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
611
612			if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
613			else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
614			}
615	mmeineke	29
616
617	mmeineke	32	void rotate( double &x, double &y, double &z,
618			double theta, double phi, double psi ){
619
620			double newX, newY, newZ;
621
622			double A[3][3];
623
624			A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
625			A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
626			A[0][2] = sin(theta) * sin(psi);
627	mmeineke	29
628	mmeineke	32	A[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
629			A[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
630			A[1][2] = sin(theta) * cos(psi);
631
632			A[2][0] = sin(phi) * sin(theta);
633			A[2][1] = -cos(phi) * sin(theta);
634			A[2][2] = cos(theta);
635
636			newX = (x * A[0][0]) + (y * A[0][1]) + (z * A[0][2]);
637			newY = (x * A[1][0]) + (y * A[1][1]) + (z * A[1][2]);
638			newZ = (x * A[2][0]) + (y * A[2][1]) + (z * A[2][2]);
639
640			x = newX;
641			y = newY;
642			z = newZ;
643	mmeineke	29	}