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 *= 1.2; // add a little buffer

  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:	506
Committed:	Fri Apr 25 16:02:26 2003 UTC (21 years, 4 months ago) by mmeineke
File size:	5724 byte(s)
Log Message:	i quick fix to th distance in the random bilayer builder
#	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	506	cell *= 1.2; // add a little buffer
131
132	mmeineke	502	cell *= M_SQRT2;
133	mmeineke	498
134	mmeineke	501	siteIndex = 0;
135			for(i=0; i<nCells; i++){
136			for(j=0; j<nCells; j++){
137			for(k=0; k<nCells; k++){
138
139			if( molMap[siteIndex] >= 0 ){
140			pos[0] = i * cell;
141			pos[1] = j * cell;
142			pos[2] = k * cell;
143
144			getRandomRot( rot );
145			molID = molSeq[molMap[siteIndex]];
146	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
147			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
148	mmeineke	501	}
149			siteIndex++;
150	mmeineke	498
151	mmeineke	501	if( molMap[siteIndex] >= 0 ){
152			pos[0] = i * cell + (0.5 * cell);
153			pos[1] = j * cell;
154			pos[2] = k * cell + (0.5 * cell);
155	mmeineke	504
156	mmeineke	501	getRandomRot( rot );
157			molID = molSeq[molMap[siteIndex]];
158	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
159			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
160	mmeineke	501	}
161			siteIndex++;
162	mmeineke	498
163	mmeineke	501	if( molMap[siteIndex] >= 0 ){
164			pos[0] = i * cell + (0.5 * cell);
165			pos[1] = j * cell + (0.5 * cell);
166			pos[2] = k * cell;
167
168			getRandomRot( rot );
169			molID = molSeq[molMap[siteIndex]];
170	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
171			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
172	mmeineke	501	}
173			siteIndex++;
174	mmeineke	498
175	mmeineke	501	if( molMap[siteIndex] >= 0 ){
176			pos[0] = i * cell;
177			pos[1] = j * cell + (0.5 * cell);
178			pos[2] = k * cell + (0.5 * cell);
179	mmeineke	504
180	mmeineke	501	getRandomRot( rot );
181			molID = molSeq[molMap[siteIndex]];
182	mmeineke	504	atomIndex = molStart[ molMap[siteIndex] ];
183			locate[molID]->placeMol( pos, rot, atoms, atomIndex );
184	mmeineke	501	}
185			siteIndex++;
186			}
187			}
188			}
189	mmeineke	498
190	mmeineke	501	// set up the SimInfo object
191
192			bsInfo.boxX = nCells * cell;
193			bsInfo.boxY = nCells * cell;
194			bsInfo.boxZ = nCells * cell;
195
196			simnfo = new SimInfo();
197	mmeineke	502	simnfo->n_atoms = nAtoms;
198			simnfo->box_x = bsInfo.boxX;
199			simnfo->box_y = bsInfo.boxY;
200			simnfo->box_z = bsInfo.boxZ;
201	mmeineke	501
202	mmeineke	504	sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
203	mmeineke	502	sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
204	mmeineke	504
205	mmeineke	502	simnfo->atoms = atoms;
206	mmeineke	501
207			// set up the writer and write out
208
209	mmeineke	502	writer = new DumpWriter( simnfo );
210	mmeineke	501	writer->writeFinal();
211	mmeineke	498
212	mmeineke	501	// clean up the memory
213
214			if( molMap != NULL ) delete[] molMap;
215			if( cardDeck != NULL ) delete[] cardDeck;
216			if( locate != NULL ){
217			for(i=0; i<bsInfo.nComponents; i++){
218			delete locate[i];
219	mmeineke	498	}
220	mmeineke	501	delete[] locate;
221			}
222			if( atoms != NULL ){
223			for(i=0; i<nAtoms; i++){
224			delete atoms[i];
225	mmeineke	498	}
226	mmeineke	501	Atom::destroyArrays();
227			delete[] atoms;
228			}
229			if( molSeq != NULL ) delete[] molSeq;
230			if( simnfo != NULL ) delete simnfo;
231			if( writer != NULL ) delete writer;
232	mmeineke	498
233	mmeineke	501	return 1;
234			}
235	mmeineke	498
236
237	mmeineke	501	void getRandomRot( double rot[3][3] ){
238	mmeineke	498
239	mmeineke	501	double theta, phi, psi;
240			double cosTheta;
241	mmeineke	498
242	mmeineke	501	// select random phi, psi, and cosTheta
243	mmeineke	498
244	mmeineke	501	phi = 2.0 * M_PI * drand48();
245			psi = 2.0 * M_PI * drand48();
246			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
247	mmeineke	498
248	mmeineke	501	theta = acos( cosTheta );
249	mmeineke	498
250	mmeineke	501	rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
251			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
252	mmeineke	502	rot[0][2] = sin(theta) * sin(psi);
253	mmeineke	501
254			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
255			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
256			rot[1][2] = sin(theta) * cos(psi);
257
258			rot[2][0] = sin(phi) * sin(theta);
259			rot[2][1] = -cos(phi) * sin(theta);
260			rot[2][2] = cos(theta);
261	mmeineke	498	}
262	mmeineke	501