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.
#	User	Rev	Content
1	chuckv	678	#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