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;

  startX = 0.0;
  startY = 0.0;
  startZ = 0.0;

  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;

  case ORTHORHOMBIC_LATTICE_TYPE:
    hasError = createFccLattice(latticeParameter);
    break;


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

Lattice::Lattice(int latticeType, double latticeParameters[3]){
  int hasError;

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

  startX = 0.0;
  startY = 0.0;
  startZ = 0.0;

  switch(latticeType){
    
  case ORTHORHOMBIC_LATTICE_TYPE:
    hasError = createOrthorhombicLattice(latticeParameters);
    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;
   
   cellLengthX = latticeParameter;
   cellLengthY = latticeParameter;
   cellLengthZ = 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;
 }


 int Lattice::createOrthorhombicLattice(double latticeParameters[3]){
   double cell2;
   double rroot3;
   
   cellLengthX = latticeParameter[0];
   cellLengthY = latticeParameter[1];
   cellLengthZ = latticeParameter[2];

   cellx2  = 0.5 * latticeParameter[0];
   cellx2  = 0.5 * latticeParameter[1];
   cellx2  = 0.5 * latticeParameter[2];
   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]   = celly2;
   thisUnitCell.sz[1]   = cellz2;
   thisUnitCell.s_ex[1] = -rroot3;
   thisUnitCell.s_ey[1] = rroot3;
   thisUnitCell.s_ez[1] = -rroot3;
     
     // Molecule 3
   thisUnitCell.sx[2]   = cellx2; 
   thisUnitCell.sy[2]   = celly2;
   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]   = cellx2;
   thisUnitCell.sy[3]   = 0.0;
   thisUnitCell.sz[3]   = cellz2;
   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] = startX + thisUnitCell.sx[iref] + 
       cellLengthX * (double( ix ) - 0.5);
     latticePosY[iref] = startY + thisUnitCell.sy[iref] + 
       cellLengthY * (double( iy ) - 0.5);
     latticePosZ[iref] = startZ + thisUnitCell.sz[iref] + 
       cellLengthZ * (double( iz ) - 0.5);
     
   }

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

   return nCellSites;
 }

Revision:	817
Committed:	Fri Oct 24 17:36:18 2003 UTC (20 years, 8 months ago) by gezelter
File size:	8509 byte(s)
Log Message:	work on bilayer builder
#	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	gezelter	817	startX = 0.0;
15			startY = 0.0;
16			startZ = 0.0;
17	chuckv	678
18			switch(latticeType){
19
20			case FCC_LATTICE_TYPE:
21	chuckv	700	hasError = createFccLattice(latticeParameter);
22	chuckv	678	break;
23
24			case BCC_LATTICE_TYPE:
25	chuckv	700	hasError = createBccLattice(latticeParameter);
26	chuckv	678	break;
27
28			case HCP_LATTICE_TYPE:
29	chuckv	700	hasError = createHcpLattice(latticeParameter);
30	chuckv	678	break;
31
32			case HCPWATER_LATTICE_TYPE:
33	chuckv	700	hasError = createHcpWaterLattice(latticeParameter);
34	chuckv	678	break;
35
36	gezelter	817	case ORTHORHOMBIC_LATTICE_TYPE:
37			hasError = createFccLattice(latticeParameter);
38			break;
39
40
41	chuckv	678	default:
42	chuckv	700	//simerror here.....
43			;
44	chuckv	678	}
45	chuckv	700
46			latticePosX = new double[nCellSites];
47			latticePosY = new double[nCellSites];
48			latticePosZ = new double[nCellSites];
49
50	chuckv	678	return;
51			}
52
53	gezelter	817	Lattice::Lattice(int latticeType, double latticeParameters[3]){
54			int hasError;
55	chuckv	678
56	gezelter	817	latticePosX = NULL;
57			latticePosY = NULL;
58			latticePosZ = NULL;
59
60			startX = 0.0;
61			startY = 0.0;
62			startZ = 0.0;
63
64			switch(latticeType){
65
66			case ORTHORHOMBIC_LATTICE_TYPE:
67			hasError = createOrthorhombicLattice(latticeParameters);
68			break;
69
70
71			default:
72			//simerror here.....
73			;
74			}
75
76			latticePosX = new double[nCellSites];
77			latticePosY = new double[nCellSites];
78			latticePosZ = new double[nCellSites];
79
80			return;
81			}
82
83
84	chuckv	678	Lattice::~Lattice(void){
85			if (latticePosX != NULL) delete[] latticePosX;
86			if (latticePosY != NULL) delete[] latticePosY;
87			if (latticePosZ != NULL) delete[] latticePosZ;
88			}
89
90
91	chuckv	700	int Lattice::createFccLattice(double latticeParameter){
92	chuckv	678	double cell2;
93			double rroot3;
94
95	gezelter	817	cellLengthX = latticeParameter;
96			cellLengthY = latticeParameter;
97			cellLengthZ = latticeParameter;
98	chuckv	700
99			cell2 = 0.5 * latticeParameter;
100	chuckv	678	rroot3 = 1.0 / sqrt(3.0);
101
102			nCellSites = 4;
103			// create new unit cells
104			thisUnitCell.sx = new double[nCellSites];
105			thisUnitCell.sy = new double[nCellSites];
106			thisUnitCell.sz = new double[nCellSites];
107			thisUnitCell.s_ex = new double[nCellSites];
108			thisUnitCell.s_ey = new double[nCellSites];
109			thisUnitCell.s_ez = new double[nCellSites];
110
111			// add members to each unit cell
112			// Molecule 1
113			thisUnitCell.sx[0] = 0.0;
114			thisUnitCell.sy[0] = 0.0;
115			thisUnitCell.sz[0] = 0.0;
116			thisUnitCell.s_ex[0] = rroot3;
117			thisUnitCell.s_ey[0] = rroot3;
118			thisUnitCell.s_ez[0] = rroot3;
119
120			// Molecule 2
121			thisUnitCell.sx[1] = 0.0;
122			thisUnitCell.sy[1] = cell2;
123			thisUnitCell.sz[1] = cell2;
124			thisUnitCell.s_ex[1] = -rroot3;
125			thisUnitCell.s_ey[1] = rroot3;
126			thisUnitCell.s_ez[1] = -rroot3;
127
128			// Molecule 3
129			thisUnitCell.sx[2] = cell2;
130			thisUnitCell.sy[2] = cell2;
131			thisUnitCell.sz[2] = 0.0;
132			thisUnitCell.s_ex[2] = rroot3;
133			thisUnitCell.s_ey[2] = -rroot3;
134			thisUnitCell.s_ez[2] = -rroot3;
135
136			// Molecule 4
137			thisUnitCell.sx[3] = cell2;
138			thisUnitCell.sy[3] = 0.0;
139			thisUnitCell.sz[3] = cell2;
140			thisUnitCell.s_ex[3] = -rroot3;
141			thisUnitCell.s_ey[3] = -rroot3;
142			thisUnitCell.s_ez[3] = rroot3;
143
144			return 0;
145			}
146
147
148	gezelter	817	int Lattice::createOrthorhombicLattice(double latticeParameters[3]){
149			double cell2;
150			double rroot3;
151
152			cellLengthX = latticeParameter[0];
153			cellLengthY = latticeParameter[1];
154			cellLengthZ = latticeParameter[2];
155
156			cellx2 = 0.5 * latticeParameter[0];
157			cellx2 = 0.5 * latticeParameter[1];
158			cellx2 = 0.5 * latticeParameter[2];
159			rroot3 = 1.0 / sqrt(3.0);
160
161			nCellSites = 4;
162			// create new unit cells
163			thisUnitCell.sx = new double[nCellSites];
164			thisUnitCell.sy = new double[nCellSites];
165			thisUnitCell.sz = new double[nCellSites];
166			thisUnitCell.s_ex = new double[nCellSites];
167			thisUnitCell.s_ey = new double[nCellSites];
168			thisUnitCell.s_ez = new double[nCellSites];
169
170			// add members to each unit cell
171			// Molecule 1
172			thisUnitCell.sx[0] = 0.0;
173			thisUnitCell.sy[0] = 0.0;
174			thisUnitCell.sz[0] = 0.0;
175			thisUnitCell.s_ex[0] = rroot3;
176			thisUnitCell.s_ey[0] = rroot3;
177			thisUnitCell.s_ez[0] = rroot3;
178
179			// Molecule 2
180			thisUnitCell.sx[1] = 0.0;
181			thisUnitCell.sy[1] = celly2;
182			thisUnitCell.sz[1] = cellz2;
183			thisUnitCell.s_ex[1] = -rroot3;
184			thisUnitCell.s_ey[1] = rroot3;
185			thisUnitCell.s_ez[1] = -rroot3;
186
187			// Molecule 3
188			thisUnitCell.sx[2] = cellx2;
189			thisUnitCell.sy[2] = celly2;
190			thisUnitCell.sz[2] = 0.0;
191			thisUnitCell.s_ex[2] = rroot3;
192			thisUnitCell.s_ey[2] = -rroot3;
193			thisUnitCell.s_ez[2] = -rroot3;
194
195			// Molecule 4
196			thisUnitCell.sx[3] = cellx2;
197			thisUnitCell.sy[3] = 0.0;
198			thisUnitCell.sz[3] = cellz2;
199			thisUnitCell.s_ex[3] = -rroot3;
200			thisUnitCell.s_ey[3] = -rroot3;
201			thisUnitCell.s_ez[3] = rroot3;
202
203			return 0;
204			}
205
206
207	chuckv	678	// Body centered cubic lattice
208	chuckv	700	int Lattice::createBccLattice(double latticeParameter){
209	chuckv	678	return 0;
210			}
211
212
213			// Standard HCP lattice
214	chuckv	700	int Lattice::createHcpLattice(double latticeParameter){
215	chuckv	678	return 0;
216			}
217
218			// HCP contains tetrahedral sites for waters
219	chuckv	700	int Lattice::createHcpWaterLattice(double latticeParameter){
220
221			return 0;
222
223	chuckv	678	double rroot3;
224			double cell;
225			double cell2;
226			double cell4;
227			double cell16;
228			double cell23;
229			double cell34;
230			double cell09;
231			double cell59;
232
233			nCellSites = 16;
234
235			// create new unit cells
236			thisUnitCell.sx = new double[nCellSites];
237			thisUnitCell.sy = new double[nCellSites];
238			thisUnitCell.sz = new double[nCellSites];
239			thisUnitCell.s_ex = new double[nCellSites];
240			thisUnitCell.s_ey = new double[nCellSites];
241			thisUnitCell.s_ez = new double[nCellSites];
242
243			rroot3 = 1.0 / sqrt(3.0);
244	chuckv	700	cell = latticeParameter;
245	chuckv	678	cell2 = 0.5 * cell;
246			cell4 = 0.25 * cell;
247			cell16 = cell / 6.0;
248			cell23 = (2.0 * cell) / 3.0;
249			cell34 = 0.75 * cell;
250			cell09 = 0.11238 * cell;
251			cell59 = 0.61238 * cell;
252
253			// build the unit cell
254
255			// molecule 1
256
257			thisUnitCell.sx[0] = 0.0;
258			thisUnitCell.sy[0] = 0.0;
259			thisUnitCell.sz[0] = 0.0;
260
261			// molecule 2
262
263			thisUnitCell.sx[1] = cell4;
264			thisUnitCell.sy[1] = cell2;
265			thisUnitCell.sz[1] = 0.0;
266
267			// molecule 3
268
269			thisUnitCell.sx[2] = cell4;
270			thisUnitCell.sy[2] = cell16;
271			thisUnitCell.sz[2] = cell09;
272
273			// molecule 4
274
275			thisUnitCell.sx[3] = 0.0;
276			thisUnitCell.sy[3] = cell23;
277			thisUnitCell.sz[3] = cell09;
278
279			// molecule 5
280
281			thisUnitCell.sx[4] = cell4;
282			thisUnitCell.sy[4] = cell16;
283			thisUnitCell.sz[4] = cell2;
284
285			// molecule 6
286
287			thisUnitCell.sx[5] = 0.0;
288			thisUnitCell.sy[5] = cell23;
289			thisUnitCell.sz[5] = cell2;
290
291			// molecule 7
292
293			thisUnitCell.sx[6] = 0.0;
294			thisUnitCell.sy[6] = 0.0;
295			thisUnitCell.sz[6] = cell59;
296
297
298			// molecule 8
299
300			thisUnitCell.sx[7] = cell4;
301			thisUnitCell.sy[7] = cell2;
302			thisUnitCell.sz[7] = cell59;
303
304			// molecule 9
305
306			thisUnitCell.sz[8] = cell2;
307			thisUnitCell.sy[8] = 0.0;
308			thisUnitCell.sz[8] = 0.0;
309
310			// molecule 10
311
312			thisUnitCell.sz[9] = cell34;
313			thisUnitCell.sy[9] = cell2;
314			thisUnitCell.sz[9] = 0.0;
315
316			// molecule 11
317
318			thisUnitCell.sz[10] = cell34;
319			thisUnitCell.sy[10] = cell16;
320			thisUnitCell.sz[10] = cell09;
321
322			// molecule 12
323
324			thisUnitCell.sz[11] = cell2;
325			thisUnitCell.sy[11] = cell23;
326			thisUnitCell.sz[11] = cell09;
327
328			// molecule 13
329
330			thisUnitCell.sz[12] = cell34;
331			thisUnitCell.sy[12] = cell16;
332			thisUnitCell.sz[12] = cell2;
333
334			// molecule 14
335
336			thisUnitCell.sz[13]= cell2;
337			thisUnitCell.sy[13]= cell23;
338			thisUnitCell.sz[13]= cell2;
339
340			// molecule 15
341
342			thisUnitCell.sz[14] = cell2;
343			thisUnitCell.sy[14] = 0.0;
344			thisUnitCell.sz[14] = cell59;
345
346			// molecule 16
347
348			thisUnitCell.sz[15] = cell34;
349			thisUnitCell.sy[15] = cell2;
350			thisUnitCell.sz[15] = cell59;
351
352			return 0;
353			}
354
355
356			//Returns lattice points when called repeatedly
357			int Lattice::getLatticePoints(double thePosX, double thePosY,
358	chuckv	700	double** thePosZ,
359			int ix, int iy, int iz){
360	chuckv	678
361			int iref;
362
363			for( iref=0;iref < nCellSites;iref++){
364
365	gezelter	817	latticePosX[iref] = startX + thisUnitCell.sx[iref] +
366			cellLengthX * (double( ix ) - 0.5);
367			latticePosY[iref] = startY + thisUnitCell.sy[iref] +
368			cellLengthY * (double( iy ) - 0.5);
369			latticePosZ[iref] = startZ + thisUnitCell.sz[iref] +
370			cellLengthZ * (double( iz ) - 0.5);
371	chuckv	678
372			}
373
374			*thePosX = latticePosX;
375			*thePosY = latticePosY;
376			*thePosZ = latticePosZ;
377
378			return nCellSites;
379			}
380