utils/sysbuilder/latticeBuilder.cpp

#include "latticeBuilder.hpp"
#include <cmath>
#include <cstdlib>
#include "simError.h"


Lattice::Lattice(int latticeType, double latticeParameter){
  int hasError;

  latticePosX = NULL;
  latticePosY = NULL;
  latticePosZ = NULL;


  switch(latticeType){
    
  case FCC_LATTICE_TYPE:
    hasError = createFccLattice(latticeParameter);
    break;
    
  case BCC_LATTICE_TYPE:
    hasError = createBccLattice(latticeParameter);
    break;
    
  case HCP_LATTICE_TYPE:
    hasError = createHcpLattice(latticeParameter);
    break;

  case HCPWATER_LATTICE_TYPE:
    hasError = createHcpWaterLattice(latticeParameter);
    break;

  default:
    //simerror here.....
    ;
  }
  
  latticePosX = new double[nCellSites];
  latticePosY = new double[nCellSites];
  latticePosZ = new double[nCellSites];
  
  return;
}


Lattice::~Lattice(void){
  if (latticePosX != NULL) delete[] latticePosX;
  if (latticePosY != NULL) delete[] latticePosY;
  if (latticePosZ != NULL) delete[] latticePosZ;
}


 int Lattice::createFccLattice(double latticeParameter){
   double cell2;
   double rroot3;
   
   cellLength = latticeParameter;

   cell2  = 0.5 * latticeParameter;
   rroot3 = 1.0 / sqrt(3.0);

   nCellSites = 4;
   // create new unit cells
   thisUnitCell.sx   = new double[nCellSites];
   thisUnitCell.sy   = new double[nCellSites];
   thisUnitCell.sz   = new double[nCellSites];
   thisUnitCell.s_ex = new double[nCellSites];
   thisUnitCell.s_ey = new double[nCellSites];
   thisUnitCell.s_ez = new double[nCellSites];

   // add members to each unit cell
   // Molecule 1
   thisUnitCell.sx[0]   = 0.0;
   thisUnitCell.sy[0]   = 0.0;
   thisUnitCell.sz[0]   = 0.0;
   thisUnitCell.s_ex[0] = rroot3;
   thisUnitCell.s_ey[0] = rroot3;
   thisUnitCell.s_ez[0] = rroot3;

     // Molecule 2
   thisUnitCell.sx[1]   = 0.0;
   thisUnitCell.sy[1]   = cell2;
   thisUnitCell.sz[1]   = cell2;
   thisUnitCell.s_ex[1] = -rroot3;
   thisUnitCell.s_ey[1] = rroot3;
   thisUnitCell.s_ez[1] = -rroot3;
     
     // Molecule 3
   thisUnitCell.sx[2]   = cell2; 
   thisUnitCell.sy[2]   = cell2;
   thisUnitCell.sz[2]   = 0.0;
   thisUnitCell.s_ex[2] = rroot3;
   thisUnitCell.s_ey[2] = -rroot3;
   thisUnitCell.s_ez[2] = -rroot3;

     // Molecule 4
   thisUnitCell.sx[3]   = cell2;
   thisUnitCell.sy[3]   = 0.0;
   thisUnitCell.sz[3]   = cell2;
   thisUnitCell.s_ex[3] = -rroot3;
   thisUnitCell.s_ey[3] = -rroot3;
   thisUnitCell.s_ez[3] = rroot3;

     return 0;
 }


 // Body centered cubic lattice
 int Lattice::createBccLattice(double latticeParameter){
   return 0;
 }


 // Standard HCP lattice
 int Lattice::createHcpLattice(double latticeParameter){
   return 0;
 }

 // HCP contains tetrahedral sites for waters 
 int Lattice::createHcpWaterLattice(double latticeParameter){
 
   return 0;
   
   double rroot3;
   double cell;
   double cell2;
   double cell4;
   double cell16;
   double cell23;
   double cell34;
   double cell09;
   double cell59;

   nCellSites = 16;

   // create new unit cells
   thisUnitCell.sx   = new double[nCellSites];
   thisUnitCell.sy   = new double[nCellSites];
   thisUnitCell.sz   = new double[nCellSites];
   thisUnitCell.s_ex = new double[nCellSites];
   thisUnitCell.s_ey = new double[nCellSites];
   thisUnitCell.s_ez = new double[nCellSites];
  
   rroot3 = 1.0 / sqrt(3.0);
   cell   = latticeParameter;
   cell2  = 0.5  * cell;
   cell4  = 0.25 * cell;
   cell16 = cell / 6.0;
   cell23 = (2.0 * cell) / 3.0;
   cell34 = 0.75 * cell;
   cell09 = 0.11238 * cell;
   cell59 = 0.61238 * cell;

   // ** build the unit cell **
   
   // molecule 1

   thisUnitCell.sx[0] =  0.0;
   thisUnitCell.sy[0] =  0.0;
   thisUnitCell.sz[0] =  0.0;
 
   // molecule 2

   thisUnitCell.sx[1] =  cell4;
   thisUnitCell.sy[1] =  cell2;
   thisUnitCell.sz[1] =  0.0;

  // molecule 3

   thisUnitCell.sx[2] =  cell4;
   thisUnitCell.sy[2] =  cell16;
   thisUnitCell.sz[2] =  cell09;

   // molecule 4

   thisUnitCell.sx[3] =  0.0;
   thisUnitCell.sy[3] =  cell23;
   thisUnitCell.sz[3] =  cell09;

  // molecule 5

   thisUnitCell.sx[4] =  cell4;
   thisUnitCell.sy[4] =  cell16;
   thisUnitCell.sz[4] =  cell2;

   // molecule 6

   thisUnitCell.sx[5] =  0.0;
   thisUnitCell.sy[5] =  cell23;
   thisUnitCell.sz[5] =  cell2;

   // molecule 7
   
   thisUnitCell.sx[6] =  0.0;
   thisUnitCell.sy[6] =  0.0;
   thisUnitCell.sz[6] =  cell59;


   // molecule 8

   thisUnitCell.sx[7] =  cell4;
   thisUnitCell.sy[7] =  cell2;
   thisUnitCell.sz[7] =  cell59;

   // molecule 9

   thisUnitCell.sz[8] =  cell2;
   thisUnitCell.sy[8] =  0.0;
   thisUnitCell.sz[8] =  0.0;

   // molecule 10

   thisUnitCell.sz[9] =  cell34;
   thisUnitCell.sy[9] =  cell2;
   thisUnitCell.sz[9] =  0.0;

   // molecule 11

   thisUnitCell.sz[10] =  cell34;
   thisUnitCell.sy[10] =  cell16;
   thisUnitCell.sz[10] =  cell09;

   // molecule 12

   thisUnitCell.sz[11] =  cell2;
   thisUnitCell.sy[11] =  cell23;
   thisUnitCell.sz[11] =  cell09;

   // molecule 13

   thisUnitCell.sz[12] =  cell34;
   thisUnitCell.sy[12] =  cell16;
   thisUnitCell.sz[12] =  cell2;

   // molecule 14

   thisUnitCell.sz[13]=  cell2;
   thisUnitCell.sy[13]=  cell23;
   thisUnitCell.sz[13]=  cell2;

   // molecule 15

   thisUnitCell.sz[14] =  cell2;
   thisUnitCell.sy[14] =  0.0;
   thisUnitCell.sz[14] =  cell59;

   // molecule 16

   thisUnitCell.sz[15] =  cell34;
   thisUnitCell.sy[15] =  cell2;
   thisUnitCell.sz[15] =  cell59;

   return 0;
 }


 //Returns lattice points when called repeatedly
 int Lattice::getLatticePoints(double** thePosX, double** thePosY,
                double** thePosZ,
                int ix, int iy, int iz){

   int iref;
  
   for( iref=0;iref < nCellSites;iref++){

     latticePosX[iref] = thisUnitCell.sx[iref] + 
       cellLength * (double( ix ) - 0.5);
     latticePosY[iref] = thisUnitCell.sy[iref] + 
       cellLength * (double( iy ) - 0.5);
     latticePosZ[iref] = thisUnitCell.sz[iref] + 
       cellLength * (double( iz ) - 0.5);
     
   }

   *thePosX = latticePosX;
   *thePosY = latticePosY;
   *thePosZ = latticePosZ;

   return nCellSites;
 }

Revision:	700
Committed:	Mon Aug 18 20:59:47 2003 UTC (20 years, 11 months ago) by chuckv
File size:	6096 byte(s)
Log Message:	Fixed sysBuild -bilayer works. Nanobuilder still broke.
#	User	Rev	Content
1	chuckv	678	#include "latticeBuilder.hpp"
2			#include <cmath>
3			#include <cstdlib>
4			#include "simError.h"
5
6
7	chuckv	700	Lattice::Lattice(int latticeType, double latticeParameter){
8	chuckv	678	int hasError;
9
10			latticePosX = NULL;
11			latticePosY = NULL;
12			latticePosZ = NULL;
13
14
15			switch(latticeType){
16
17			case FCC_LATTICE_TYPE:
18	chuckv	700	hasError = createFccLattice(latticeParameter);
19	chuckv	678	break;
20
21			case BCC_LATTICE_TYPE:
22	chuckv	700	hasError = createBccLattice(latticeParameter);
23	chuckv	678	break;
24
25			case HCP_LATTICE_TYPE:
26	chuckv	700	hasError = createHcpLattice(latticeParameter);
27	chuckv	678	break;
28
29			case HCPWATER_LATTICE_TYPE:
30	chuckv	700	hasError = createHcpWaterLattice(latticeParameter);
31	chuckv	678	break;
32
33			default:
34	chuckv	700	//simerror here.....
35			;
36	chuckv	678	}
37	chuckv	700
38			latticePosX = new double[nCellSites];
39			latticePosY = new double[nCellSites];
40			latticePosZ = new double[nCellSites];
41
42	chuckv	678	return;
43			}
44
45
46			Lattice::~Lattice(void){
47			if (latticePosX != NULL) delete[] latticePosX;
48			if (latticePosY != NULL) delete[] latticePosY;
49			if (latticePosZ != NULL) delete[] latticePosZ;
50			}
51
52
53	chuckv	700	int Lattice::createFccLattice(double latticeParameter){
54	chuckv	678	double cell2;
55			double rroot3;
56
57	chuckv	700	cellLength = latticeParameter;
58
59			cell2 = 0.5 * latticeParameter;
60	chuckv	678	rroot3 = 1.0 / sqrt(3.0);
61
62			nCellSites = 4;
63			// create new unit cells
64			thisUnitCell.sx = new double[nCellSites];
65			thisUnitCell.sy = new double[nCellSites];
66			thisUnitCell.sz = new double[nCellSites];
67			thisUnitCell.s_ex = new double[nCellSites];
68			thisUnitCell.s_ey = new double[nCellSites];
69			thisUnitCell.s_ez = new double[nCellSites];
70
71			// add members to each unit cell
72			// Molecule 1
73			thisUnitCell.sx[0] = 0.0;
74			thisUnitCell.sy[0] = 0.0;
75			thisUnitCell.sz[0] = 0.0;
76			thisUnitCell.s_ex[0] = rroot3;
77			thisUnitCell.s_ey[0] = rroot3;
78			thisUnitCell.s_ez[0] = rroot3;
79
80			// Molecule 2
81			thisUnitCell.sx[1] = 0.0;
82			thisUnitCell.sy[1] = cell2;
83			thisUnitCell.sz[1] = cell2;
84			thisUnitCell.s_ex[1] = -rroot3;
85			thisUnitCell.s_ey[1] = rroot3;
86			thisUnitCell.s_ez[1] = -rroot3;
87
88			// Molecule 3
89			thisUnitCell.sx[2] = cell2;
90			thisUnitCell.sy[2] = cell2;
91			thisUnitCell.sz[2] = 0.0;
92			thisUnitCell.s_ex[2] = rroot3;
93			thisUnitCell.s_ey[2] = -rroot3;
94			thisUnitCell.s_ez[2] = -rroot3;
95
96			// Molecule 4
97			thisUnitCell.sx[3] = cell2;
98			thisUnitCell.sy[3] = 0.0;
99			thisUnitCell.sz[3] = cell2;
100			thisUnitCell.s_ex[3] = -rroot3;
101			thisUnitCell.s_ey[3] = -rroot3;
102			thisUnitCell.s_ez[3] = rroot3;
103
104			return 0;
105			}
106
107
108			// Body centered cubic lattice
109	chuckv	700	int Lattice::createBccLattice(double latticeParameter){
110	chuckv	678	return 0;
111			}
112
113
114			// Standard HCP lattice
115	chuckv	700	int Lattice::createHcpLattice(double latticeParameter){
116	chuckv	678	return 0;
117			}
118
119			// HCP contains tetrahedral sites for waters
120	chuckv	700	int Lattice::createHcpWaterLattice(double latticeParameter){
121
122			return 0;
123
124	chuckv	678	double rroot3;
125			double cell;
126			double cell2;
127			double cell4;
128			double cell16;
129			double cell23;
130			double cell34;
131			double cell09;
132			double cell59;
133
134			nCellSites = 16;
135
136			// create new unit cells
137			thisUnitCell.sx = new double[nCellSites];
138			thisUnitCell.sy = new double[nCellSites];
139			thisUnitCell.sz = new double[nCellSites];
140			thisUnitCell.s_ex = new double[nCellSites];
141			thisUnitCell.s_ey = new double[nCellSites];
142			thisUnitCell.s_ez = new double[nCellSites];
143
144			rroot3 = 1.0 / sqrt(3.0);
145	chuckv	700	cell = latticeParameter;
146	chuckv	678	cell2 = 0.5 * cell;
147			cell4 = 0.25 * cell;
148			cell16 = cell / 6.0;
149			cell23 = (2.0 * cell) / 3.0;
150			cell34 = 0.75 * cell;
151			cell09 = 0.11238 * cell;
152			cell59 = 0.61238 * cell;
153
154			// build the unit cell
155
156			// molecule 1
157
158			thisUnitCell.sx[0] = 0.0;
159			thisUnitCell.sy[0] = 0.0;
160			thisUnitCell.sz[0] = 0.0;
161
162			// molecule 2
163
164			thisUnitCell.sx[1] = cell4;
165			thisUnitCell.sy[1] = cell2;
166			thisUnitCell.sz[1] = 0.0;
167
168			// molecule 3
169
170			thisUnitCell.sx[2] = cell4;
171			thisUnitCell.sy[2] = cell16;
172			thisUnitCell.sz[2] = cell09;
173
174			// molecule 4
175
176			thisUnitCell.sx[3] = 0.0;
177			thisUnitCell.sy[3] = cell23;
178			thisUnitCell.sz[3] = cell09;
179
180			// molecule 5
181
182			thisUnitCell.sx[4] = cell4;
183			thisUnitCell.sy[4] = cell16;
184			thisUnitCell.sz[4] = cell2;
185
186			// molecule 6
187
188			thisUnitCell.sx[5] = 0.0;
189			thisUnitCell.sy[5] = cell23;
190			thisUnitCell.sz[5] = cell2;
191
192			// molecule 7
193
194			thisUnitCell.sx[6] = 0.0;
195			thisUnitCell.sy[6] = 0.0;
196			thisUnitCell.sz[6] = cell59;
197
198
199			// molecule 8
200
201			thisUnitCell.sx[7] = cell4;
202			thisUnitCell.sy[7] = cell2;
203			thisUnitCell.sz[7] = cell59;
204
205			// molecule 9
206
207			thisUnitCell.sz[8] = cell2;
208			thisUnitCell.sy[8] = 0.0;
209			thisUnitCell.sz[8] = 0.0;
210
211			// molecule 10
212
213			thisUnitCell.sz[9] = cell34;
214			thisUnitCell.sy[9] = cell2;
215			thisUnitCell.sz[9] = 0.0;
216
217			// molecule 11
218
219			thisUnitCell.sz[10] = cell34;
220			thisUnitCell.sy[10] = cell16;
221			thisUnitCell.sz[10] = cell09;
222
223			// molecule 12
224
225			thisUnitCell.sz[11] = cell2;
226			thisUnitCell.sy[11] = cell23;
227			thisUnitCell.sz[11] = cell09;
228
229			// molecule 13
230
231			thisUnitCell.sz[12] = cell34;
232			thisUnitCell.sy[12] = cell16;
233			thisUnitCell.sz[12] = cell2;
234
235			// molecule 14
236
237			thisUnitCell.sz[13]= cell2;
238			thisUnitCell.sy[13]= cell23;
239			thisUnitCell.sz[13]= cell2;
240
241			// molecule 15
242
243			thisUnitCell.sz[14] = cell2;
244			thisUnitCell.sy[14] = 0.0;
245			thisUnitCell.sz[14] = cell59;
246
247			// molecule 16
248
249			thisUnitCell.sz[15] = cell34;
250			thisUnitCell.sy[15] = cell2;
251			thisUnitCell.sz[15] = cell59;
252
253			return 0;
254			}
255
256
257			//Returns lattice points when called repeatedly
258			int Lattice::getLatticePoints(double thePosX, double thePosY,
259	chuckv	700	double** thePosZ,
260			int ix, int iy, int iz){
261	chuckv	678
262			int iref;
263
264			for( iref=0;iref < nCellSites;iref++){
265
266			latticePosX[iref] = thisUnitCell.sx[iref] +
267			cellLength * (double( ix ) - 0.5);
268			latticePosY[iref] = thisUnitCell.sy[iref] +
269			cellLength * (double( iy ) - 0.5);
270			latticePosZ[iref] = thisUnitCell.sz[iref] +
271			cellLength * (double( iz ) - 0.5);
272
273			}
274
275			*thePosX = latticePosX;
276			*thePosY = latticePosY;
277			*thePosZ = latticePosZ;
278
279			return nCellSites;
280			}
281