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 (22 years, 2 months ago) by chuckv
File size:	6096 byte(s)
Log Message:	Fixed sysBuild -bilayer works. Nanobuilder still broke.
#	Content
1	#include "latticeBuilder.hpp"
2	#include <cmath>
3	#include <cstdlib>
4	#include "simError.h"
5
6
7	Lattice::Lattice(int latticeType, double latticeParameter){
8	int hasError;
9
10	latticePosX = NULL;
11	latticePosY = NULL;
12	latticePosZ = NULL;
13
14
15	switch(latticeType){
16
17	case FCC_LATTICE_TYPE:
18	hasError = createFccLattice(latticeParameter);
19	break;
20
21	case BCC_LATTICE_TYPE:
22	hasError = createBccLattice(latticeParameter);
23	break;
24
25	case HCP_LATTICE_TYPE:
26	hasError = createHcpLattice(latticeParameter);
27	break;
28
29	case HCPWATER_LATTICE_TYPE:
30	hasError = createHcpWaterLattice(latticeParameter);
31	break;
32
33	default:
34	//simerror here.....
35	;
36	}
37
38	latticePosX = new double[nCellSites];
39	latticePosY = new double[nCellSites];
40	latticePosZ = new double[nCellSites];
41
42	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	int Lattice::createFccLattice(double latticeParameter){
54	double cell2;
55	double rroot3;
56
57	cellLength = latticeParameter;
58
59	cell2 = 0.5 * latticeParameter;
60	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	int Lattice::createBccLattice(double latticeParameter){
110	return 0;
111	}
112
113
114	// Standard HCP lattice
115	int Lattice::createHcpLattice(double latticeParameter){
116	return 0;
117	}
118
119	// HCP contains tetrahedral sites for waters
120	int Lattice::createHcpWaterLattice(double latticeParameter){
121
122	return 0;
123
124	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	cell = latticeParameter;
146	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	double** thePosZ,
260	int ix, int iy, int iz){
261
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