utils/sysbuilder/bilayerSys.cpp


#include <iostream>
#include <vector>
#include <algorithm>

#include <cstdlib>
#include <cstring>
#include <cmath>


#include "simError.h"
#include "SimInfo.hpp"
#include "ReadWrite.hpp"

#include "MoLocator.hpp"
#include "sysBuild.hpp"
#include "bilayerSys.hpp"

#include "latticeBuilder.hpp"

class SortCond{
  
public:
  bool operator()(const pair<int, double>& p1, const pair<int, double>& p2){
    return p1.second < p2.second;
  }
  
  
};


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

int buildRandomBilayer( void );

void getRandomRot( double rot[3][3] );

int buildBilayer( int isRandom ){
  
  if( isRandom ){
    return buildRandomBilayer();
  }
  else{
    sprintf( painCave.errMsg,
             "Cannot currently create a non-random bilayer.\n" );
    painCave.isFatal = 1;
    simError();
    return 0;
  }
}


int buildRandomBilayer( void ){

  typedef struct{
    double rot[3][3];
    double pos[3];
  } coord;


  const double waterRho = 0.0334; // number density per cubic angstrom
  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
  const double waterCell = 4.929; // fcc unit cell length

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

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


  int i,j,k, l, m;
  int nAtoms, atomIndex, molIndex, molID;
  int* molSeq;
  int* molMap;
  int* molStart;
  int* cardDeck;
  int deckSize;
  int rSite, rCard;
  double cell;
  int nCells, nSites, siteIndex;
  
  coord testSite;

  Atom** atoms;
  SimInfo* simnfo;
  SimState* theConfig;
  DumpWriter* writer;

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


  // create the simInfo objects

  simnfo = new SimInfo[3];


  // set the the lipidStamp

  foundLipid = 0;
  foundWater = 0;
  for(i=0; i<bsInfo.nComponents; i++){
    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
          
      foundLipid = 1;
      lipidStamp = bsInfo.compStamps[i];
      nLipids = bsInfo.componentsNmol[i];
    }
    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
          
      foundWater = 1;
          
      waterStamp = bsInfo.compStamps[i];
      nWaters = bsInfo.componentsNmol[i];
    }
  }
  if( !foundLipid ){
    sprintf(painCave.errMsg,
            "Could not find lipid \"%s\" in the bass file.\n",
            bsInfo.lipidName );
    painCave.isFatal = 1;
    simError();
  }
  if( !foundWater ){
    sprintf(painCave.errMsg,
            "Could not find solvent \"%s\" in the bass file.\n",
            bsInfo.waterName );
    painCave.isFatal = 1;
    simError();
  }
  
  //create the temp Molocator and atom Arrays
  
  lipidLocate = new MoLocator( lipidStamp );
  lipidNatoms = lipidStamp->getNAtoms();
  maxLength = lipidLocate->getMaxLength();

  waterLocate = new MoLocator( waterStamp );
  waterNatoms = waterStamp->getNAtoms();
  
  nAtoms = lipidNatoms;

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

  atoms=simnfo[0].atoms;
        

  // create the test box for initial water displacement

  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
  nCells = (int)( testBox / waterCell + 1.0 );
  int testWaters = 4 * nCells * nCells * nCells;
  
  double* waterX = new double[testWaters];
  double* waterY = new double[testWaters];
  double* waterZ = new double[testWaters];
  
  double x0 = 0.0 - ( testBox * 0.5 );
  double y0 = 0.0 - ( testBox * 0.5 );
  double z0 = 0.0 - ( testBox * 0.5 );


  // create an fcc lattice in the water box.

  int ndx = 0;
  for( i=0; i < nCells; i++ ){
    for( j=0; j < nCells; j++ ){
      for( k=0; k < nCells; k++ ){

        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
        for(l=0; l<4; l++){
          waterX[ndx]=posX[l];
          waterY[ndx]=posY[l];
          waterZ[ndx]=posZ[l];
          ndx++;
        }
      }
    }
  }

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

  testSite.rot[0][0] = 1.0;
  testSite.rot[0][1] = 0.0;
  testSite.rot[0][2] = 0.0;

  testSite.rot[1][0] = 0.0;
  testSite.rot[1][1] = 1.0;
  testSite.rot[1][2] = 0.0;

  testSite.rot[2][0] = 0.0;
  testSite.rot[2][1] = 0.0;
  testSite.rot[2][2] = 1.0;
  
  testSite.pos[0] = 0.0;
  testSite.pos[1] = 0.0;
  testSite.pos[2] = 0.0;

  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );

  int *isActive = new int[testWaters];
  for(i=0; i<testWaters; i++) isActive[i] = 1;
  
  int n_deleted = 0;
  double dx, dy, dz;
  double dx2, dy2, dz2, dSqr;
  double rCutSqr = water_padding * water_padding;
  
  for(i=0; ( (i<testWaters) && isActive[i] ); i++){
    for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
      
      atoms[j]->getPos( pos );

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

      buildMap( dx, dy, dz, testBox, testBox, testBox );

      dx2 = dx * dx;
      dy2 = dy * dy;
      dz2 = dz * dz;
          
      dSqr = dx2 + dy2 + dz2;
      if( dSqr < rCutSqr ){
        isActive[i] = 0;
        n_deleted++;
      }
    }
  }
  
  int targetWaters = nWaters + n_deleted * nLipids;

  targetWaters = (int) ( targetWaters * 1.2 );

  // find the best box size for the sim

  int nCellsX, nCellsY, nCellsZ;

  const double boxTargetX = 20;
  const double boxTargetY = 20;

  nCellsX = (int)ceil(boxTargetX / waterCell);
  nCellsY = (int)ceil(boxTargetY / waterCell);
  
  int testTot;
  int done = 0;
  nCellsZ = 0;
  while( !done ){
        
    nCellsZ++;
    testTot = 4 * nCellsX * nCellsY * nCellsZ;
        
    if( testTot >= targetWaters ) done = 1;
  }

  // create the new water box to the new specifications
  
  int newWaters = nCellsX * nCellsY * nCellsZ * 4;
  
  delete[] waterX;
  delete[] waterY;
  delete[] waterZ;
  
  coord* waterSites = new coord[newWaters];
  
  double box_x = waterCell * nCellsX;
  double box_y = waterCell * nCellsY;
  double box_z = waterCell * nCellsZ;
         
  
  // create an fcc lattice in the water box.
  
  ndx = 0;
  for( i=0; i < nCellsX; i++ ){
    for( j=0; j < nCellsY; j++ ){
      for( k=0; k < nCellsZ; k++ ){

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

  coord* lipidSites = new coord[nLipids];

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

  nAtoms = nLipids * lipidNatoms;

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

  atoms=simnfo[1].atoms;

  rCutSqr = lipid_spaceing * lipid_spaceing;

  for(i=0; i<nLipids; i++ ){
    done = 0;
    while( !done ){
          
      lipidSites[i].pos[0] = drand48() * box_x;
      lipidSites[i].pos[1] = drand48() * box_y;
      lipidSites[i].pos[2] = drand48() * box_z;
          
      getRandomRot( lipidSites[i].rot );
          
      ndx = i * lipidNatoms;

      lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, 
                             ndx, theConfig );
          
      reject = 0;
      for( j=0; !reject && j<i; j++){
        for(k=0; !reject && k<lipidNatoms; k++){
          
          acceptedDX = j*lipidNatoms + k;
          for(l=0; !reject && l<lipidNatoms; l++){
                
            testDX = ndx + l;

            atoms[testDX]->getPos( posA );
            atoms[acceptedDX]->getPos( posB );
            
            dx = posA[0] - posB[0];
            dy = posA[1] - posB[1];
            dz = posA[2] - posB[2];
                
            buildMap( dx, dy, dz, box_x, box_y, box_z );
                
            dx2 = dx * dx;
            dy2 = dy * dy;
            dz2 = dz * dz;
                
            dSqr = dx2 + dy2 + dz2;
            if( dSqr < rCutSqr ) reject = 1;
          }
        }
      }

      if( reject ){

        for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
      }
      else{
        done = 1;
        std::cout << (i+1) << " has been accepted\n";
      }
    }
  }
        

  // zSort of the lipid positions


  vector< pair<int,double> >zSortArray;
  for(i=0;i<nLipids;i++) 
    zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );

  sort(zSortArray.begin(),zSortArray.end(),SortCond());

  ofstream outFile( "./zipper.bass", ios::app);

  for(i=0; i<nLipids; i++){
    outFile << "zConstraint[" << i << "]{\n"
            << "  molIndex = " << zSortArray[i].first <<  ";\n"
            << "  zPos = ";

    if(i<32) outFile << "60.0;\n";
    else outFile << "100.0;\n";

    outFile << "  kRatio = 0.5;\n"
            << "}\n";
  }
  
  outFile.close();


  // cut out the waters that overlap with the lipids.
  

  delete[] isActive;
  isActive = new int[newWaters];
  for(i=0; i<newWaters; i++) isActive[i] = 1;
  int n_active = newWaters;
  rCutSqr = water_padding * water_padding;
  
  for(i=0; ( (i<newWaters) && isActive[i] ); i++){
    for(j=0; ( (j<nAtoms) && isActive[i] ); j++){

      atoms[j]->getPos( pos );

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

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

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


      }
    }
  }


  if( n_active < nWaters ){
        
    sprintf( painCave.errMsg,
             "Too many waters were removed, edit code and try again.\n" );
        
    painCave.isFatal = 1;
    simError();
  }

  int quickKill;
  while( n_active > nWaters ){

    quickKill = (int)(drand48()*newWaters);

    if( isActive[quickKill] ){
      isActive[quickKill] = 0;
      n_active--;

    }
  }

  if( n_active != nWaters ){
        
    sprintf( painCave.errMsg,
             "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
             n_active, nWaters );
    painCave.isFatal = 1;
    simError();
  }

  // clean up our messes before building the final system.

  simnfo[0].getConfiguration()->destroyArrays();
  simnfo[1].getConfiguration()->destroyArrays();

  // create the real Atom arrays
  
  nAtoms = 0;
  molIndex = 0;
  molStart = new int[nLipids + nWaters];  
  
  for(j=0; j<nLipids; j++){
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += lipidNatoms;
  }

  for(j=0; j<nWaters; j++){
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += waterNatoms;
  }
  
  theConfig = simnfo[2].getConfiguration();
  theConfig->createArrays( nAtoms );
  simnfo[2].atoms = new Atom*[nAtoms];
  atoms = simnfo[2].atoms;
  simnfo[2].n_atoms = nAtoms;
  
  // initialize lipid positions

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

  // initialize the water positions

  for(i=0; i<newWaters; i++){
        
    if( isActive[i] ){
          
      getRandomRot( waterSites[i].rot );
      waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
                             molStart[molIndex], theConfig );
      molIndex++;
    }
  }  

  // set up the SimInfo object
  
  double Hmat[3][3];

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

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

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

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

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

  // set up the writer and write out
  
  writer = new DumpWriter( &simnfo[2] );
  writer->writeFinal( 0.0 );
        
  // clean up the memory
  
  //     if( molMap != NULL )   delete[] molMap;
  //     if( cardDeck != NULL ) delete[] cardDeck;
  //     if( locate != NULL ){
  //       for(i=0; i<bsInfo.nComponents; i++){
  //             delete locate[i];
  //       }
  //       delete[] locate;
  //     }
  //     if( atoms != NULL ){
  //       for(i=0; i<nAtoms; i++){
  //             delete atoms[i];
  //       }
  //       Atom::destroyArrays();
  //       delete[] atoms;
  //     }
  //     if( molSeq != NULL ) delete[] molSeq;
  //     if( simnfo != NULL ) delete simnfo;
  //     if( writer != NULL ) delete writer;

  return 1;
}

void getRandomRot( double rot[3][3] ){

  double theta, phi, psi;
  double cosTheta;

  // select random phi, psi, and cosTheta

  phi = 2.0 * M_PI * drand48();
  psi = 2.0 * M_PI * drand48();
  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1

  theta = acos( cosTheta );

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

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

void buildMap( double &x, double &y, double &z, 
               double boxX, double boxY, double boxZ ){ 
  
  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));

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