1 /******************************************************************************
2
3 This source file is part of the Avogadro project.
4
5 Copyright 2013 Kitware, Inc.
6
7 This source code is released under the New BSD License, (the "License").
8
9 Unless required by applicable law or agreed to in writing, software
10 distributed under the License is distributed on an "AS IS" BASIS,
11 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 See the License for the specific language governing permissions and
13 limitations under the License.
14
15 ******************************************************************************/
16
17 #ifndef AVOGADRO_CORE_UNITCELL_H
18 #define AVOGADRO_CORE_UNITCELL_H
19
20 #include "avogadrocore.h"
21
22 #include "matrix.h"
23 #include "vector.h"
24
25 namespace Avogadro {
26 namespace Core {
27
28 /**
29 * @class UnitCell unitcell.h <avogadro/core/unitcell.h>
30 * @brief The UnitCell class provides a representation of a crystal's unit cell.
31 */
32 class AVOGADROCORE_EXPORT UnitCell
33 {
34 public:
35 UnitCell();
36 UnitCell(Real a, Real b, Real c, Real alpha, Real beta, Real gamma);
37 UnitCell(const Vector3& a, const Vector3& b, const Vector3& c);
38 explicit UnitCell(const Matrix3& cellMatrix);
39 UnitCell(const UnitCell& other);
40 ~UnitCell();
41 UnitCell& operator=(UnitCell other);
42 friend void swap(UnitCell& lhs, UnitCell& rhs);
43
44 /** The lattice vector in the unit cell. Units: Angstrom @{ */
aVector()45 Vector3 aVector() const { return m_cellMatrix.col(0); }
bVector()46 Vector3 bVector() const { return m_cellMatrix.col(1); }
cVector()47 Vector3 cVector() const { return m_cellMatrix.col(2); }
48 void setAVector(const Vector3& v);
49 void setBVector(const Vector3& v);
50 void setCVector(const Vector3& v);
51 /** @} */
52
53 /** The length of the lattice vector in the unit cell. Units: Angstrom @{ */
a()54 Real a() const { return m_cellMatrix.col(0).norm(); }
b()55 Real b() const { return m_cellMatrix.col(1).norm(); }
c()56 Real c() const { return m_cellMatrix.col(2).norm(); }
57 /** @} */
58
59 /** The angle (radians) between the 'b' and 'c' lattice vectors. */
60 Real alpha() const;
61
62 /** The angle (radians) between the 'c' and 'a' lattice vectors. */
63 Real beta() const;
64
65 /** The angle (radians) between the 'a' and 'b' lattice vectors. */
66 Real gamma() const;
67
68 /**
69 * Set the cell parameters defining the unit cell. @a a, @a b, and @a c are
70 * in Angstrom, @a alpha, @a beta, and @a gamma are in radians.
71 */
72 void setCellParameters(Real a, Real b, Real c, Real alpha, Real beta,
73 Real gamma);
74
75 /**
76 * The volume of the unit cell in cubic Angstroms.
77 */
78 Real volume() const;
79
80 /**
81 * @return A vector pointing to the origin of the translational image that is
82 * @a i images in the a() direction, @a j images in the b() direction, and
83 * @a k images in the c() direction.
84 */
85 Vector3 imageOffset(int i, int j, int k) const;
86
87 /**
88 * The cell matrix with lattice vectors as columns. Units: Angstrom @{
89 */
90 const Matrix3& cellMatrix() const;
91 void setCellMatrix(const Matrix3& m);
92 /** @} */
93
94 /**
95 * The matrix used to convert cartesian to fractional coordinates.
96 */
97 const Matrix3& fractionalMatrix() const;
98 void setFractionalMatrix(const Matrix3& m);
99 /** @} */
100
101 /**
102 * Convert the cartesian coordinate @a cart to fractional (lattice) units. @{
103 */
104 Vector3 toFractional(const Vector3& cart) const;
105 void toFractional(const Vector3& cart, Vector3& frac) const;
106
107 /**
108 * Convert the fractional (lattice) coordinate @a frac to cartesian units. @{
109 */
110 Vector3 toCartesian(const Vector3& frac) const;
111 void toCartesian(const Vector3& frac, Vector3& cart) const;
112 /** @} */
113
114 /**
115 * Adjust the fractional (lattice) coordinate @a frac so that it lies within
116 * the unit cell. @{
117 */
118 Vector3 wrapFractional(const Vector3& frac) const;
119 void wrapFractional(const Vector3& frac, Vector3& wrapped) const;
120 /** @} */
121
122 /**
123 * Adjust the cartesian coordinate @a cart so that it lies within the unit
124 * cell. @{
125 */
126 Vector3 wrapCartesian(const Vector3& cart) const;
127 void wrapCartesian(const Vector3& cart, Vector3& wrapped) const;
128 /** @} */
129
130 /**
131 * Find the minimum fractional image of a fractional vector @a v.
132 * A minimum image has all fractional coordinates between -0.5 and 0.5.
133 */
134 static Vector3 minimumImageFractional(const Vector3& v);
135
136 /**
137 * Find the minimum image of a Cartesian vector @a v.
138 * A minimum image has all fractional coordinates between -0.5 and 0.5
139 */
140 Vector3 minimumImage(const Vector3& v) const;
141
142 /**
143 * Find the shortest distance between vectors @a v1 and @a v2.
144 */
145 Real distance(const Vector3& v1, const Vector3& v2) const;
146
147 private:
148 static Real signedAngleRadians(const Vector3& v1, const Vector3& v2,
149 const Vector3& axis);
computeCellMatrix()150 void computeCellMatrix() { m_cellMatrix = m_fractionalMatrix.inverse(); }
computeFractionalMatrix()151 void computeFractionalMatrix()
152 {
153 m_fractionalMatrix = m_cellMatrix.inverse();
154 }
155
156 Matrix3 m_cellMatrix;
157 Matrix3 m_fractionalMatrix;
158 };
159
UnitCell()160 inline UnitCell::UnitCell()
161 : m_cellMatrix(Matrix3::Identity()), m_fractionalMatrix(Matrix3::Identity())
162 {
163 }
164
UnitCell(Real a_,Real b_,Real c_,Real alpha_,Real beta_,Real gamma_)165 inline UnitCell::UnitCell(Real a_, Real b_, Real c_, Real alpha_, Real beta_,
166 Real gamma_)
167 {
168 setCellParameters(a_, b_, c_, alpha_, beta_, gamma_);
169 }
170
UnitCell(const Vector3 & a_,const Vector3 & b_,const Vector3 & c_)171 inline UnitCell::UnitCell(const Vector3& a_, const Vector3& b_,
172 const Vector3& c_)
173 {
174 m_cellMatrix.col(0) = a_;
175 m_cellMatrix.col(1) = b_;
176 m_cellMatrix.col(2) = c_;
177 computeFractionalMatrix();
178 }
179
UnitCell(const Matrix3 & cellMatrix_)180 inline UnitCell::UnitCell(const Matrix3& cellMatrix_)
181 {
182 m_cellMatrix = cellMatrix_;
183 computeFractionalMatrix();
184 }
185
UnitCell(const UnitCell & other)186 inline UnitCell::UnitCell(const UnitCell& other)
187 : m_cellMatrix(other.m_cellMatrix),
188 m_fractionalMatrix(other.m_fractionalMatrix)
189 {
190 }
191
~UnitCell()192 inline UnitCell::~UnitCell()
193 {
194 }
195
196 inline UnitCell& UnitCell::operator=(UnitCell other)
197 {
198 using std::swap;
199 swap(*this, other);
200 return *this;
201 }
202
swap(UnitCell & lhs,UnitCell & rhs)203 inline void swap(UnitCell& lhs, UnitCell& rhs)
204 {
205 using std::swap;
206 swap(lhs.m_cellMatrix, rhs.m_cellMatrix);
207 swap(lhs.m_fractionalMatrix, rhs.m_fractionalMatrix);
208 }
209
setAVector(const Vector3 & v)210 inline void UnitCell::setAVector(const Vector3& v)
211 {
212 m_cellMatrix.col(0) = v;
213 computeFractionalMatrix();
214 }
215
setBVector(const Vector3 & v)216 inline void UnitCell::setBVector(const Vector3& v)
217 {
218 m_cellMatrix.col(1) = v;
219 computeFractionalMatrix();
220 }
221
setCVector(const Vector3 & v)222 inline void UnitCell::setCVector(const Vector3& v)
223 {
224 m_cellMatrix.col(2) = v;
225 computeFractionalMatrix();
226 }
227
alpha()228 inline Real UnitCell::alpha() const
229 {
230 return signedAngleRadians(bVector(), cVector(), aVector());
231 }
232
beta()233 inline Real UnitCell::beta() const
234 {
235 return signedAngleRadians(cVector(), aVector(), bVector());
236 }
237
gamma()238 inline Real UnitCell::gamma() const
239 {
240 return signedAngleRadians(aVector(), bVector(), cVector());
241 }
242
volume()243 inline Real UnitCell::volume() const
244 {
245 return std::fabs(aVector().cross(bVector()).dot(cVector()));
246 }
247
imageOffset(int i,int j,int k)248 inline Vector3 UnitCell::imageOffset(int i, int j, int k) const
249 {
250 return (static_cast<Real>(i) * m_cellMatrix.col(0) +
251 static_cast<Real>(j) * m_cellMatrix.col(1) +
252 static_cast<Real>(k) * m_cellMatrix.col(2));
253 }
254
cellMatrix()255 inline const Matrix3& UnitCell::cellMatrix() const
256 {
257 return m_cellMatrix;
258 }
259
setCellMatrix(const Matrix3 & m)260 inline void UnitCell::setCellMatrix(const Matrix3& m)
261 {
262 m_cellMatrix = m;
263 computeFractionalMatrix();
264 }
265
fractionalMatrix()266 inline const Matrix3& UnitCell::fractionalMatrix() const
267 {
268 return m_fractionalMatrix;
269 }
270
setFractionalMatrix(const Matrix3 & m)271 inline void UnitCell::setFractionalMatrix(const Matrix3& m)
272 {
273 m_fractionalMatrix = m;
274 computeCellMatrix();
275 }
276
toFractional(const Vector3 & cart)277 inline Vector3 UnitCell::toFractional(const Vector3& cart) const
278 {
279 return m_fractionalMatrix * cart;
280 }
281
toFractional(const Vector3 & cart,Vector3 & frac)282 inline void UnitCell::toFractional(const Vector3& cart, Vector3& frac) const
283 {
284 frac = m_fractionalMatrix * cart;
285 }
286
toCartesian(const Vector3 & f)287 inline Vector3 UnitCell::toCartesian(const Vector3& f) const
288 {
289 return m_cellMatrix * f;
290 }
291
toCartesian(const Vector3 & frac,Vector3 & cart)292 inline void UnitCell::toCartesian(const Vector3& frac, Vector3& cart) const
293 {
294 cart = m_cellMatrix * frac;
295 }
296
wrapFractional(const Vector3 & f)297 inline Vector3 UnitCell::wrapFractional(const Vector3& f) const
298 {
299 const Real one = static_cast<Real>(1.0);
300 Vector3 result(std::fmod(f[0], one), std::fmod(f[1], one),
301 std::fmod(f[2], one));
302 if (result[0] < static_cast<Real>(0.0))
303 ++result[0];
304 if (result[1] < static_cast<Real>(0.0))
305 ++result[1];
306 if (result[2] < static_cast<Real>(0.0))
307 ++result[2];
308 return result;
309 }
310
wrapFractional(const Vector3 & f,Vector3 & wrapped)311 inline void UnitCell::wrapFractional(const Vector3& f, Vector3& wrapped) const
312 {
313 const Real one = static_cast<Real>(1.0);
314 wrapped =
315 Vector3(std::fmod(f[0], one), std::fmod(f[1], one), std::fmod(f[2], one));
316 if (wrapped[0] < static_cast<Real>(0.0))
317 ++wrapped[0];
318 if (wrapped[1] < static_cast<Real>(0.0))
319 ++wrapped[1];
320 if (wrapped[2] < static_cast<Real>(0.0))
321 ++wrapped[2];
322 }
323
wrapCartesian(const Vector3 & cart)324 inline Vector3 UnitCell::wrapCartesian(const Vector3& cart) const
325 {
326 Vector3 result = toFractional(cart);
327 wrapFractional(result, result);
328 toCartesian(result, result);
329 return result;
330 }
331
wrapCartesian(const Vector3 & cart,Vector3 & wrapped)332 inline void UnitCell::wrapCartesian(const Vector3& cart, Vector3& wrapped) const
333 {
334 toFractional(cart, wrapped);
335 wrapFractional(wrapped, wrapped);
336 toCartesian(wrapped, wrapped);
337 }
338
minimumImageFractional(const Vector3 & v)339 inline Vector3 UnitCell::minimumImageFractional(const Vector3& v)
340 {
341 Real x = v[0] - rint(v[0]);
342 Real y = v[1] - rint(v[1]);
343 Real z = v[2] - rint(v[2]);
344 return Vector3(x, y, z);
345 }
346
minimumImage(const Vector3 & v)347 inline Vector3 UnitCell::minimumImage(const Vector3& v) const
348 {
349 return toCartesian(minimumImageFractional(toFractional(v)));
350 }
351
distance(const Vector3 & v1,const Vector3 & v2)352 inline Real UnitCell::distance(const Vector3& v1, const Vector3& v2) const
353 {
354 return std::fabs(minimumImage(v1 - v2).norm());
355 }
356
357 } // namespace Core
358 } // namespace Avogadro
359
360 #endif // AVOGADRO_CORE_UNITCELL_H
361