utils/sysbuilder/latticeBuilder.cpp

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


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

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


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

  case HCPWATER_LATTICE_TYPE:
    hasError = createHcpWaterLattice(latticeSpacing);
    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 latticeSpacing){
   double cell2;
   double rroot3;
   
   cell2  = 0.5 * latticeSpacing;
   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 latticeSpacing){
   return 0;
 }


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

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