OOPSE/utils/bilayerSys.cpp

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