# | Line 48 | Line 48 | void SimInfo::setBox(double newBox[3]) { | |
---|---|---|
48 | ||
49 | void SimInfo::setBox(double newBox[3]) { | |
50 | ||
51 | < | int i; |
52 | < | double tempMat[9]; |
51 | > | int i, j; |
52 | > | double tempMat[3][3]; |
53 | ||
54 | < | for(i=0; i<9; i++) tempMat[i] = 0.0;; |
54 | > | for(i=0; i<3; i++) |
55 | > | for (j=0; j<3; j++) tempMat[i][j] = 0.0;; |
56 | ||
57 | < | tempMat[0] = newBox[0]; |
58 | < | tempMat[4] = newBox[1]; |
59 | < | tempMat[8] = newBox[2]; |
57 | > | tempMat[0][0] = newBox[0]; |
58 | > | tempMat[1][1] = newBox[1]; |
59 | > | tempMat[2][2] = newBox[2]; |
60 | ||
61 | setBoxM( tempMat ); | |
62 | ||
63 | } | |
64 | ||
65 | < | void SimInfo::setBoxM( double theBox[9] ){ |
65 | > | void SimInfo::setBoxM( double theBox[3][3] ){ |
66 | ||
67 | < | int i, status; |
67 | > | int i, j, status; |
68 | double smallestBoxL, maxCutoff; | |
69 | + | double FortranHmat[9]; // to preserve compatibility with Fortran the |
70 | + | // ordering in the array is as follows: |
71 | + | // [ 0 3 6 ] |
72 | + | // [ 1 4 7 ] |
73 | + | // [ 2 5 8 ] |
74 | + | double FortranHmatInv[9]; // the inverted Hmat (for Fortran); |
75 | ||
69 | – | for(i=0; i<9; i++) Hmat[i] = theBox[i]; |
76 | ||
77 | + | for(i=0; i < 3; i++) |
78 | + | for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; |
79 | + | |
80 | cerr | |
81 | << "setting Hmat ->\n" | |
82 | < | << "[ " << Hmat[0] << ", " << Hmat[3] << ", " << Hmat[6] << " ]\n" |
83 | < | << "[ " << Hmat[1] << ", " << Hmat[4] << ", " << Hmat[7] << " ]\n" |
84 | < | << "[ " << Hmat[2] << ", " << Hmat[5] << ", " << Hmat[8] << " ]\n"; |
82 | > | << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n" |
83 | > | << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n" |
84 | > | << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n"; |
85 | ||
77 | – | calcHmatI(); |
86 | calcBoxL(); | |
87 | + | calcHmatInv(); |
88 | ||
89 | + | for(i=0; i < 3; i++) { |
90 | + | for (j=0; j < 3; j++) { |
91 | + | FortranHmat[3*j + i] = Hmat[i][j]; |
92 | + | FortranHmatInv[3*j + i] = HmatInv[i][j]; |
93 | + | } |
94 | + | } |
95 | ||
96 | < | |
82 | < | setFortranBoxSize(Hmat, HmatI, &orthoRhombic); |
96 | > | setFortranBoxSize(FortranHmat, FortranHmatI, &orthoRhombic); |
97 | ||
98 | smallestBoxL = boxLx; | |
99 | if (boxLy < smallestBoxL) smallestBoxL = boxLy; | |
# | Line 127 | Line 141 | void SimInfo::setBoxM( double theBox[9] ){ | |
141 | } | |
142 | ||
143 | ||
144 | < | void SimInfo::getBoxM (double theBox[9]) { |
144 | > | void SimInfo::getBoxM (double theBox[3][3]) { |
145 | ||
146 | < | int i; |
147 | < | for(i=0; i<9; i++) theBox[i] = Hmat[i]; |
146 | > | int i, j; |
147 | > | for(i=0; i<3; i++) |
148 | > | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; |
149 | } | |
150 | ||
151 | ||
152 | void SimInfo::scaleBox(double scale) { | |
153 | < | double theBox[9]; |
154 | < | int i; |
153 | > | double theBox[3][3]; |
154 | > | int i, j; |
155 | ||
156 | cerr << "Scaling box by " << scale << "\n"; | |
157 | ||
158 | < | for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale; |
158 | > | for(i=0; i<3; i++) |
159 | > | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; |
160 | ||
161 | setBoxM(theBox); | |
162 | ||
163 | } | |
164 | ||
165 | < | void SimInfo::calcHmatI( void ) { |
165 | > | void SimInfo::calcHmatInv( void ) { |
166 | ||
151 | – | double C[3][3]; |
152 | – | double detHmat; |
153 | – | int i, j, k; |
167 | double smallDiag; | |
168 | double tol; | |
169 | double sanity[3][3]; | |
170 | ||
171 | < | // calculate the adjunct of Hmat; |
171 | > | invertMat3( Hmat, HmatInv ); |
172 | ||
173 | < | C[0][0] = ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]); |
161 | < | C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]); |
162 | < | C[2][0] = ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]); |
173 | > | // Check the inverse to make sure it is sane: |
174 | ||
175 | < | C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]); |
165 | < | C[1][1] = ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]); |
166 | < | C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]); |
175 | > | matMul3( Hmat, HmatInv, sanity ); |
176 | ||
168 | – | C[0][2] = ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]); |
169 | – | C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]); |
170 | – | C[2][2] = ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]); |
171 | – | |
172 | – | // calcutlate the determinant of Hmat |
173 | – | |
174 | – | detHmat = 0.0; |
175 | – | for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0]; |
176 | – | |
177 | – | |
178 | – | // H^-1 = C^T / det(H) |
179 | – | |
180 | – | i=0; |
181 | – | for(j=0; j<3; j++){ |
182 | – | for(k=0; k<3; k++){ |
183 | – | |
184 | – | HmatI[i] = C[j][k] / detHmat; |
185 | – | i++; |
186 | – | } |
187 | – | } |
188 | – | |
189 | – | // sanity check |
190 | – | |
191 | – | for(i=0; i<3; i++){ |
192 | – | for(j=0; j<3; j++){ |
193 | – | |
194 | – | sanity[i][j] = 0.0; |
195 | – | for(k=0; k<3; k++){ |
196 | – | sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k]; |
197 | – | } |
198 | – | } |
199 | – | } |
200 | – | |
177 | cerr << "sanity => \n" | |
178 | << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n" | |
179 | << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n" | |
180 | << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2] | |
181 | << "\n"; | |
182 | ||
207 | – | |
183 | // check to see if Hmat is orthorhombic | |
184 | ||
185 | < | smallDiag = Hmat[0]; |
186 | < | if(smallDiag > Hmat[4]) smallDiag = Hmat[4]; |
187 | < | if(smallDiag > Hmat[8]) smallDiag = Hmat[8]; |
185 | > | smallDiag = Hmat[0][0]; |
186 | > | if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1]; |
187 | > | if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2]; |
188 | tol = smallDiag * 1E-6; | |
189 | ||
190 | orthoRhombic = 1; | |
191 | < | for(i=0; (i<9) && orthoRhombic; i++){ |
192 | < | |
193 | < | if( (i%4) ){ // ignore the diagonals (0, 4, and 8) |
194 | < | orthoRhombic = (Hmat[i] <= tol); |
191 | > | |
192 | > | for (i = 0; i < 3; i++ ) { |
193 | > | for (j = 0 ; j < 3; j++) { |
194 | > | if (i != j) { |
195 | > | if (orthoRhombic) { |
196 | > | if (Hmat[i][j] >= tol) orthoRhombic = 0; |
197 | > | } |
198 | > | } |
199 | } | |
200 | } | |
222 | – | |
201 | } | |
202 | ||
203 | + | double SimInfo::matDet3(double a[3][3]) { |
204 | + | int i, j, k; |
205 | + | double determinant; |
206 | + | |
207 | + | determinant = 0.0; |
208 | + | |
209 | + | for(i = 0; i < 3; i++) { |
210 | + | j = (i+1)%3; |
211 | + | k = (i+2)%3; |
212 | + | |
213 | + | determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]); |
214 | + | } |
215 | + | |
216 | + | return determinant; |
217 | + | } |
218 | + | |
219 | + | void SimInfo::invertMat3(double a[3][3], double b[3][3]) { |
220 | + | |
221 | + | int i, j, k, l, m, n; |
222 | + | double determinant; |
223 | + | |
224 | + | determinant = matDet3( a ); |
225 | + | |
226 | + | if (determinant == 0.0) { |
227 | + | sprintf( painCave.errMsg, |
228 | + | "Can't invert a matrix with a zero determinant!\n"); |
229 | + | painCave.isFatal = 1; |
230 | + | simError(); |
231 | + | } |
232 | + | |
233 | + | for (i=0; i < 3; i++) { |
234 | + | j = (i+1)%3; |
235 | + | k = (i+2)%3; |
236 | + | for(l = 0; l < 3; l++) { |
237 | + | m = (l+1)%3; |
238 | + | n = (l+2)%3; |
239 | + | |
240 | + | b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant; |
241 | + | } |
242 | + | } |
243 | + | } |
244 | + | |
245 | + | void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) { |
246 | + | double r00, r01, r02, r10, r11, r12, r20, r21, r22; |
247 | + | |
248 | + | r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0]; |
249 | + | r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1]; |
250 | + | r02 = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2]; |
251 | + | |
252 | + | r10 = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0]; |
253 | + | r11 = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1]; |
254 | + | r12 = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2]; |
255 | + | |
256 | + | r20 = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0]; |
257 | + | r21 = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1]; |
258 | + | r22 = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2]; |
259 | + | |
260 | + | c[0][0] = r00; c[0][1] = r01; c[0][2] = r02; |
261 | + | c[1][0] = r10; c[1][1] = r11; c[1][2] = r12; |
262 | + | c[2][0] = r20; c[2][1] = r21; c[2][2] = r22; |
263 | + | } |
264 | + | |
265 | + | void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) { |
266 | + | double a0, a1, a2; |
267 | + | |
268 | + | a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2]; |
269 | + | |
270 | + | outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2; |
271 | + | outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2; |
272 | + | outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2; |
273 | + | } |
274 | + | |
275 | void SimInfo::calcBoxL( void ){ | |
276 | ||
277 | double dx, dy, dz, dsq; | |
278 | int i; | |
279 | ||
280 | < | // boxVol = h1 (dot) h2 (cross) h3 |
280 | > | // boxVol = Determinant of Hmat |
281 | ||
282 | < | boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) ) |
233 | < | + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) ) |
234 | < | + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) ); |
282 | > | boxVol = matDet3( Hmat ); |
283 | ||
236 | – | |
284 | // boxLx | |
285 | ||
286 | < | dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2]; |
286 | > | dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; |
287 | dsq = dx*dx + dy*dy + dz*dz; | |
288 | boxLx = sqrt( dsq ); | |
289 | ||
290 | // boxLy | |
291 | ||
292 | < | dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5]; |
292 | > | dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; |
293 | dsq = dx*dx + dy*dy + dz*dz; | |
294 | boxLy = sqrt( dsq ); | |
295 | ||
296 | // boxLz | |
297 | ||
298 | < | dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8]; |
298 | > | dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; |
299 | dsq = dx*dx + dy*dy + dz*dz; | |
300 | boxLz = sqrt( dsq ); | |
301 | ||
# | Line 262 | Line 309 | void SimInfo::wrapVector( double thePos[3] ){ | |
309 | ||
310 | if( !orthoRhombic ){ | |
311 | // calc the scaled coordinates. | |
312 | + | |
313 | + | |
314 | + | matVecMul3(HmatInv, thePos, scaled); |
315 | ||
316 | for(i=0; i<3; i++) | |
267 | – | scaled[i] = |
268 | – | thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6]; |
269 | – | |
270 | – | // wrap the scaled coordinates |
271 | – | |
272 | – | for(i=0; i<3; i++) |
317 | scaled[i] -= roundMe(scaled[i]); | |
318 | ||
319 | // calc the wrapped real coordinates from the wrapped scaled coordinates | |
320 | ||
321 | < | for(i=0; i<3; i++) |
322 | < | thePos[i] = |
279 | < | scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6]; |
321 | > | matVecMul3(Hmat, scaled, thePos); |
322 | > | |
323 | } | |
324 | else{ | |
325 | // calc the scaled coordinates. | |
326 | ||
327 | for(i=0; i<3; i++) | |
328 | < | scaled[i] = thePos[i]*HmatI[i*4]; |
328 | > | scaled[i] = thePos[i]*HmatInv[i][i]; |
329 | ||
330 | // wrap the scaled coordinates | |
331 | ||
# | Line 292 | Line 335 | void SimInfo::wrapVector( double thePos[3] ){ | |
335 | // calc the wrapped real coordinates from the wrapped scaled coordinates | |
336 | ||
337 | for(i=0; i<3; i++) | |
338 | < | thePos[i] = scaled[i]*Hmat[i*4]; |
338 | > | thePos[i] = scaled[i]*Hmat[i][i]; |
339 | } | |
340 | ||
298 | – | |
341 | } | |
342 | ||
343 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |