OOPSE/utils/bilayerSys.cpp

#include <iostream>

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

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

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


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 ){

  int i,j,k;
  int nAtoms, atomIndex, molIndex, molID;
  int* molSeq;
  int* molMap;
  int* molStart;
  int* cardDeck;
  int deckSize;
  int rSite, rCard;
  double cell;
  int nCells, nSites, siteIndex;
  double rot[3][3];
  double pos[3];

  Atom** atoms;
  SimInfo* simnfo;
  DumpWriter* writer;
  MoLocator** locate;

  // initialize functions and variables
  
  srand48( RAND_SEED );
  molSeq = NULL;
  molStart = NULL;
  molMap = NULL;
  cardDeck = NULL;
  atoms = NULL;
  locate = NULL;
  simnfo = NULL;
  writer = NULL;

  // calculate the number of cells in the fcc box

  nCells = 0;
  nSites = 0;
  while( nSites < bsInfo.totNmol ){
    nCells++;
    nSites = 4.0 * pow( (double)nCells, 3.0 );
  }

 
  // create the molMap and cardDeck arrays
  
  molMap = new int[nSites];
  cardDeck = new int[nSites];
  
  for(i=0; i<nSites; i++){
    molMap[i] = -1;
    cardDeck[i] = i;
  }

  // randomly place the molecules on the sites
  
  deckSize = nSites;
  for(i=0; i<bsInfo.totNmol; i++){
    rCard = (int)( deckSize * drand48() );
    rSite = cardDeck[rCard];
    molMap[rSite] = i;

    // book keep the card deck;
    
    deckSize--;
    cardDeck[rCard] = cardDeck[deckSize];
  }
  
  
  // create the MoLocator and Atom arrays
  
  nAtoms = 0;
  molIndex = 0;
  locate = new MoLocator*[bsInfo.nComponents];
  molSeq = new int[bsInfo.totNmol];
  molStart = new int[bsInfo.totNmol];  
  for(i=0; i<bsInfo.nComponents; i++){
    locate[i] = new MoLocator( bsInfo.compStamps[i] );
    for(j=0; j<bsInfo.componentsNmol[i]; j++){
      molSeq[molIndex] = i;
      molStart[molIndex] = nAtoms;
      molIndex++;
      nAtoms += bsInfo.compStamps[i]->getNAtoms();
    }
  }
  
  Atom::createArrays( nAtoms );
  atoms = new Atom*[nAtoms];
  
  
  // place the molecules at each FCC site
  
  cell = 5.0;
  for(i=0; i<bsInfo.nComponents; i++){
    if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
  }
  cell *= M_SQRT2;

  siteIndex = 0;
  for(i=0; i<nCells; i++){
    for(j=0; j<nCells; j++){
      for(k=0; k<nCells; k++){
        
        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell;
          pos[1] = j * cell;
          pos[2] = k * cell;
          
          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell + (0.5 * cell);
          pos[1] = j * cell;
          pos[2] = k * cell + (0.5 * cell);

          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell + (0.5 * cell);
          pos[1] = j * cell + (0.5 * cell);
          pos[2] = k * cell;
          
          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;

        if( molMap[siteIndex] >= 0 ){
          pos[0] = i * cell;
          pos[1] = j * cell + (0.5 * cell);
          pos[2] = k * cell + (0.5 * cell);

          getRandomRot( rot );
          molID = molSeq[molMap[siteIndex]];
          atomIndex = molStart[ molMap[siteIndex] ];
          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
        }
        siteIndex++;
      }
    }
  }

  // set up the SimInfo object
  
  bsInfo.boxX = nCells * cell;
  bsInfo.boxY = nCells * cell;
  bsInfo.boxZ = nCells * cell;
  
  simnfo = new SimInfo();
  simnfo->n_atoms = nAtoms;
  simnfo->box_x = bsInfo.boxX;
  simnfo->box_y = bsInfo.boxY;
  simnfo->box_z = bsInfo.boxZ;
  
  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );

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

  return 1;
}


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

  double theta, phi, psi;
  double cosTheta;

  // select random phi, psi, and cosTheta

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

  theta = acos( cosTheta );

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

  rot[2][0] = sin(phi) * sin(theta);
  rot[2][1] = -cos(phi) * sin(theta);
  rot[2][2] = cos(theta);  
}
                        
Revision:	504
Committed:	Thu Apr 17 21:54:18 2003 UTC (21 years, 5 months ago) by mmeineke
File size:	5685 byte(s)
Log Message:	fixed up sysBuild to where it should now build our systems
#	User	Rev	Content
1	mmeineke	504	#include <iostream>
2
3	mmeineke	498	#include <cstdlib>
4			#include <cstring>
5			#include <cmath>
6
7			#include "simError.h"
8			#include "SimInfo.hpp"
9			#include "ReadWrite.hpp"
10
11	mmeineke	501	#include "MoLocator.hpp"
12	mmeineke	498	#include "sysBuild.hpp"
13			#include "bilayerSys.hpp"
14
15
16	mmeineke	501	int buildRandomBilayer( void );
17	mmeineke	498
18	mmeineke	501	void getRandomRot( double rot[3][3] );
19	mmeineke	498
20	mmeineke	501	int buildBilayer( int isRandom ){
21	mmeineke	498
22	mmeineke	501	if( isRandom ){
23			return buildRandomBilayer();
24	mmeineke	498	}
25			else{
26			sprintf( painCave.errMsg,
27			"Cannot currently create a non-random bilayer.\n" );
28			painCave.isFatal = 1;
29			simError();
30	mmeineke	501	return 0;
31	mmeineke	498	}
32			}
33
34
35
36	mmeineke	501	int buildRandomBilayer( void ){
37	mmeineke	498
38	mmeineke	501	int i,j,k;
39	mmeineke	502	int nAtoms, atomIndex, molIndex, molID;
40	mmeineke	501	int* molSeq;
41			int* molMap;
42	mmeineke	504	int* molStart;
43	mmeineke	501	int* cardDeck;
44			int deckSize;
45			int rSite, rCard;
46			double cell;
47			int nCells, nSites, siteIndex;
48			double rot[3][3];
49			double pos[3];
50
51			Atom** atoms;
52	mmeineke	498	SimInfo* simnfo;
53	mmeineke	501	DumpWriter* writer;
54			MoLocator** locate;
55
56			// initialize functions and variables
57	mmeineke	498
58	mmeineke	501	srand48( RAND_SEED );
59			molSeq = NULL;
60	mmeineke	504	molStart = NULL;
61	mmeineke	501	molMap = NULL;
62			cardDeck = NULL;
63			atoms = NULL;
64			locate = NULL;
65			simnfo = NULL;
66			writer = NULL;
67
68			// calculate the number of cells in the fcc box
69
70			nCells = 0;
71			nSites = 0;
72			while( nSites < bsInfo.totNmol ){
73			nCells++;
74	mmeineke	502	nSites = 4.0 * pow( (double)nCells, 3.0 );
75	mmeineke	501	}
76
77
78			// create the molMap and cardDeck arrays
79	mmeineke	498
80	mmeineke	501	molMap = new int[nSites];
81			cardDeck = new int[nSites];
82	mmeineke	498
83	mmeineke	501	for(i=0; i<nSites; i++){
84			molMap[i] = -1;
85			cardDeck[i] = i;
86			}
87
88			// randomly place the molecules on the sites
89	mmeineke	498
90	mmeineke	501	deckSize = nSites;
91			for(i=0; i<bsInfo.totNmol; i++){
92			rCard = (int)( deckSize * drand48() );
93			rSite = cardDeck[rCard];
94			molMap[rSite] = i;
95
96			// book keep the card deck;
97
98			deckSize--;
99			cardDeck[rCard] = cardDeck[deckSize];
100			}
101	mmeineke	498
102
103	mmeineke	501	// create the MoLocator and Atom arrays
104
105			nAtoms = 0;
106			molIndex = 0;
107			locate = new MoLocator*[bsInfo.nComponents];
108			molSeq = new int[bsInfo.totNmol];
109	mmeineke	504	molStart = new int[bsInfo.totNmol];
110	mmeineke	501	for(i=0; i<bsInfo.nComponents; i++){
111			locate[i] = new MoLocator( bsInfo.compStamps[i] );
112			for(j=0; j<bsInfo.componentsNmol[i]; j++){
113			molSeq[molIndex] = i;
114	mmeineke	504	molStart[molIndex] = nAtoms;
115	mmeineke	501	molIndex++;
116	mmeineke	502	nAtoms += bsInfo.compStamps[i]->getNAtoms();
117	mmeineke	501	}
118			}
119
120			Atom::createArrays( nAtoms );
121			atoms = new Atom*[nAtoms];
122
123
124			// place the molecules at each FCC site
125
126			cell = 5.0;
127			for(i=0; i<bsInfo.nComponents; i++){
128			if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
129			}
130	mmeineke	502	cell *= M_SQRT2;
131	mmeineke	498
132	mmeineke	501	siteIndex = 0;
133			for(i=0; i<nCells; i++){
134			for(j=0; j<nCells; j++){
135			for(k=0; k<nCells; k++){
136
137			if( molMap[siteIndex] >= 0 ){
138			pos[0] = i * cell;
139			pos[1] = j * cell;
140			pos[2] = k * cell;
141
142			getRandomRot( rot );
143			molID = molSeq[molMap[siteIndex]];
144	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
145			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
146	mmeineke	501	}
147			siteIndex++;
148	mmeineke	498
149	mmeineke	501	if( molMap[siteIndex] >= 0 ){
150			pos[0] = i * cell + (0.5 * cell);
151			pos[1] = j * cell;
152			pos[2] = k * cell + (0.5 * cell);
153	mmeineke	504
154	mmeineke	501	getRandomRot( rot );
155			molID = molSeq[molMap[siteIndex]];
156	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
157			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
158	mmeineke	501	}
159			siteIndex++;
160	mmeineke	498
161	mmeineke	501	if( molMap[siteIndex] >= 0 ){
162			pos[0] = i * cell + (0.5 * cell);
163			pos[1] = j * cell + (0.5 * cell);
164			pos[2] = k * cell;
165
166			getRandomRot( rot );
167			molID = molSeq[molMap[siteIndex]];
168	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
169			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
170	mmeineke	501	}
171			siteIndex++;
172	mmeineke	498
173	mmeineke	501	if( molMap[siteIndex] >= 0 ){
174			pos[0] = i * cell;
175			pos[1] = j * cell + (0.5 * cell);
176			pos[2] = k * cell + (0.5 * cell);
177	mmeineke	504
178	mmeineke	501	getRandomRot( rot );
179			molID = molSeq[molMap[siteIndex]];
180	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
181			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
182	mmeineke	501	}
183			siteIndex++;
184			}
185			}
186			}
187	mmeineke	498
188	mmeineke	501	// set up the SimInfo object
189
190			bsInfo.boxX = nCells * cell;
191			bsInfo.boxY = nCells * cell;
192			bsInfo.boxZ = nCells * cell;
193
194			simnfo = new SimInfo();
195	mmeineke	502	simnfo->n_atoms = nAtoms;
196			simnfo->box_x = bsInfo.boxX;
197			simnfo->box_y = bsInfo.boxY;
198			simnfo->box_z = bsInfo.boxZ;
199	mmeineke	501
200	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
201	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
202	mmeineke	504
203	mmeineke	502	simnfo->atoms = atoms;
204	mmeineke	501
205			// set up the writer and write out
206
207	mmeineke	502	writer = new DumpWriter( simnfo );
208	mmeineke	501	writer->writeFinal();
209	mmeineke	498
210	mmeineke	501	// clean up the memory
211
212			if( molMap != NULL ) delete[] molMap;
213			if( cardDeck != NULL ) delete[] cardDeck;
214			if( locate != NULL ){
215			for(i=0; i<bsInfo.nComponents; i++){
216			delete locate[i];
217	mmeineke	498	}
218	mmeineke	501	delete[] locate;
219			}
220			if( atoms != NULL ){
221			for(i=0; i<nAtoms; i++){
222			delete atoms[i];
223	mmeineke	498	}
224	mmeineke	501	Atom::destroyArrays();
225			delete[] atoms;
226			}
227			if( molSeq != NULL ) delete[] molSeq;
228			if( simnfo != NULL ) delete simnfo;
229			if( writer != NULL ) delete writer;
230	mmeineke	498
231	mmeineke	501	return 1;
232			}
233	mmeineke	498
234
235	mmeineke	501	void getRandomRot( double rot[3][3] ){
236	mmeineke	498
237	mmeineke	501	double theta, phi, psi;
238			double cosTheta;
239	mmeineke	498
240	mmeineke	501	// select random phi, psi, and cosTheta
241	mmeineke	498
242	mmeineke	501	phi = 2.0 * M_PI * drand48();
243			psi = 2.0 * M_PI * drand48();
244			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
245	mmeineke	498
246	mmeineke	501	theta = acos( cosTheta );
247	mmeineke	498
248	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
249			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
250	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
251	mmeineke	501
252			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
253			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
254			rot[1][2] = sin(theta) * cos(psi);
255
256			rot[2][0] = sin(phi) * sin(theta);
257			rot[2][1] = -cos(phi) * sin(theta);
258			rot[2][2] = cos(theta);
259	mmeineke	498	}
260	mmeineke	501