OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SquareMatrix3.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45/**
46 * @file SquareMatrix3.hpp
47 * @author Teng Lin
48 * @date 10/11/2004
49 * @version 1.0
50 */
51
52#ifndef MATH_SQUAREMATRIX3_HPP
53#define MATH_SQUAREMATRIX3_HPP
54
55#include <config.h>
56
57#include <cmath>
58#include <limits>
59#include <vector>
60
61#include "Quaternion.hpp"
62#include "SquareMatrix.hpp"
63#include "Vector3.hpp"
64
65namespace OpenMD {
66
67 template<typename Real>
68 class SquareMatrix3 : public SquareMatrix<Real, 3> {
69 public:
70 using ElemType = Real;
71 using ElemPoinerType = Real*;
72
73 /** default constructor */
74 SquareMatrix3() : SquareMatrix<Real, 3>() {}
75
76 /** Constructs and initializes every element of this matrix to a scalar */
77 SquareMatrix3(Real s) : SquareMatrix<Real, 3>(s) {}
78
79 /** Constructs and initializes from an array */
80 SquareMatrix3(Real* array) : SquareMatrix<Real, 3>(array) {}
81
82 /** copy constructor */
84
85 SquareMatrix3(const Vector3<Real>& eulerAngles) {
86 setupRotMat(eulerAngles);
87 }
88
89 SquareMatrix3(Real phi, Real theta, Real psi) {
90 setupRotMat(phi, theta, psi);
91 }
92
93 SquareMatrix3(const Quaternion<Real>& q) { setupRotMat(q); }
94
95 SquareMatrix3(Real w, Real x, Real y, Real z) { setupRotMat(w, x, y, z); }
96
97 /** copy assignment operator */
99 if (this == &m) return *this;
101 return *this;
102 }
103
105 this->setupRotMat(q);
106 return *this;
107 }
108
109 /**
110 * Sets this matrix to a rotation matrix by three euler angles
111 * @ param euler
112 */
113 void setupRotMat(const Vector3<Real>& eulerAngles) {
114 setupRotMat(eulerAngles[0], eulerAngles[1], eulerAngles[2]);
115 }
116
117 /**
118 * Sets this matrix to a rotation matrix by three euler angles
119 * @param phi
120 * @param theta
121 * @param psi
122 */
123 void setupRotMat(Real phi, Real theta, Real psi) {
124 Real sphi, stheta, spsi;
125 Real cphi, ctheta, cpsi;
126
127 sphi = sin(phi);
128 stheta = sin(theta);
129 spsi = sin(psi);
130 cphi = cos(phi);
131 ctheta = cos(theta);
132 cpsi = cos(psi);
133
134 this->data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
135 this->data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
136 this->data_[0][2] = spsi * stheta;
137
138 this->data_[1][0] = -spsi * cphi - ctheta * sphi * cpsi;
139 this->data_[1][1] = -spsi * sphi + ctheta * cphi * cpsi;
140 this->data_[1][2] = cpsi * stheta;
141
142 this->data_[2][0] = stheta * sphi;
143 this->data_[2][1] = -stheta * cphi;
144 this->data_[2][2] = ctheta;
145 }
146
147 /**
148 * Sets this matrix to a rotation matrix by quaternion
149 * @param quat
150 */
151 void setupRotMat(const Quaternion<Real>& quat) {
152 setupRotMat(quat.w(), quat.x(), quat.y(), quat.z());
153 }
154
155 /**
156 * Sets this matrix to a rotation matrix by quaternion
157 * @param w the first element
158 * @param x the second element
159 * @param y the third element
160 * @param z the fourth element
161 */
162 void setupRotMat(Real w, Real x, Real y, Real z) {
163 Quaternion<Real> q(w, x, y, z);
164 *this = q.toRotationMatrix3();
165 }
166
167 void setupSkewMat(Vector3<Real> v) { setupSkewMat(v[0], v[1], v[2]); }
168
169 void setupSkewMat(Real v1, Real v2, Real v3) {
170 this->data_[0][0] = 0;
171 this->data_[0][1] = -v3;
172 this->data_[0][2] = v2;
173 this->data_[1][0] = v3;
174 this->data_[1][1] = 0;
175 this->data_[1][2] = -v1;
176 this->data_[2][0] = -v2;
177 this->data_[2][1] = v1;
178 this->data_[2][2] = 0;
179 }
180
181 /**
182 * Sets this matrix to a symmetric tensor using Voigt Notation
183 * @param vt
184 */
186 setupVoigtTensor(vt[0], vt[1], vt[2], vt[3], vt[4], vt[5]);
187 }
188
189 void setupVoigtTensor(Real v1, Real v2, Real v3, Real v4, Real v5,
190 Real v6) {
191 this->data_[0][0] = v1;
192 this->data_[1][1] = v2;
193 this->data_[2][2] = v3;
194 this->data_[1][2] = v4;
195 this->data_[2][1] = v4;
196 this->data_[0][2] = v5;
197 this->data_[2][0] = v5;
198 this->data_[0][1] = v6;
199 this->data_[1][0] = v6;
200 }
201
202 /**
203 * Sets this matrix to an upper-triangular (asymmetric) tensor
204 * using Voigt Notation
205 * @param vt
206 */
208 setupUpperTriangularVoigtTensor(vt[0], vt[1], vt[2], vt[3], vt[4], vt[5]);
209 }
210
211 void setupUpperTriangularVoigtTensor(Real v1, Real v2, Real v3, Real v4,
212 Real v5, Real v6) {
213 this->data_[0][0] = v1;
214 this->data_[1][1] = v2;
215 this->data_[2][2] = v3;
216 this->data_[1][2] = v4;
217 this->data_[0][2] = v5;
218 this->data_[0][1] = v6;
219 }
220
221 /**
222 * Uses Rodrigues' rotation formula for a rotation matrix.
223 * @param axis the axis to rotate around
224 * @param angle the angle to rotate (in radians)
225 */
226 void axisAngle(Vector3d axis, RealType angle) {
227 axis.normalize();
228 RealType ct = cos(angle);
229 RealType st = sin(angle);
230
231 *this = ct * SquareMatrix3<Real>::identity();
232 *this += st * SquareMatrix3<Real>::setupSkewMat(axis);
233 *this += (1 - ct) * SquareMatrix3<Real>::outProduct(axis, axis);
234 }
235
236 /**
237 * Returns the quaternion from this rotation matrix
238 * @return the quaternion from this rotation matrix
239 * @exception invalid rotation matrix
240 */
243 Real t, s;
244 Real ad1, ad2, ad3;
245 t = this->data_[0][0] + this->data_[1][1] + this->data_[2][2] + 1.0;
246
247 if (t > std::numeric_limits<RealType>::epsilon()) {
248 s = 0.5 / sqrt(t);
249 q[0] = 0.25 / s;
250 q[1] = (this->data_[1][2] - this->data_[2][1]) * s;
251 q[2] = (this->data_[2][0] - this->data_[0][2]) * s;
252 q[3] = (this->data_[0][1] - this->data_[1][0]) * s;
253 } else {
254 ad1 = this->data_[0][0];
255 ad2 = this->data_[1][1];
256 ad3 = this->data_[2][2];
257
258 if (ad1 >= ad2 && ad1 >= ad3) {
259 s = 0.5 / sqrt(1.0 + this->data_[0][0] - this->data_[1][1] -
260 this->data_[2][2]);
261 q[0] = (this->data_[1][2] - this->data_[2][1]) * s;
262 q[1] = 0.25 / s;
263 q[2] = (this->data_[0][1] + this->data_[1][0]) * s;
264 q[3] = (this->data_[0][2] + this->data_[2][0]) * s;
265 } else if (ad2 >= ad1 && ad2 >= ad3) {
266 s = 0.5 / sqrt(1.0 + this->data_[1][1] - this->data_[0][0] -
267 this->data_[2][2]);
268 q[0] = (this->data_[2][0] - this->data_[0][2]) * s;
269 q[1] = (this->data_[0][1] + this->data_[1][0]) * s;
270 q[2] = 0.25 / s;
271 q[3] = (this->data_[1][2] + this->data_[2][1]) * s;
272 } else {
273 s = 0.5 / sqrt(1.0 + this->data_[2][2] - this->data_[0][0] -
274 this->data_[1][1]);
275 q[0] = (this->data_[0][1] - this->data_[1][0]) * s;
276 q[1] = (this->data_[0][2] + this->data_[2][0]) * s;
277 q[2] = (this->data_[1][2] + this->data_[2][1]) * s;
278 q[3] = 0.25 / s;
279 }
280 }
281
282 return q;
283 }
284
285 /**
286 * Returns the euler angles from this rotation matrix
287 * @return the euler angles in a vector
288 * @exception invalid rotation matrix
289 * We use so-called "x-convention", which is the most common definition.
290 * In this convention, the rotation given by Euler angles (phi, theta,
291 * psi), where the first rotation is by an angle phi about the z-axis,
292 * the second is by an angle theta (0 <= theta <= 180) about the x-axis,
293 * and the third is by an angle psi about the z-axis (again).
294 */
296 Vector3<Real> myEuler;
297 Real phi;
298 Real theta;
299 Real psi;
300 Real ctheta;
301 Real stheta;
302
303 // set the tolerance for Euler angles and rotation elements
304
305 theta = acos(
306 std::min((RealType)1.0, std::max((RealType)-1.0, this->data_[2][2])));
307 ctheta = this->data_[2][2];
308 stheta = sqrt(1.0 - ctheta * ctheta);
309
310 // when sin(theta) is close to 0, we need to consider
311 // singularity In this case, we can assign an arbitary value to
312 // phi (or psi), and then determine the psi (or phi) or
313 // vice-versa. We'll assume that phi always gets the rotation,
314 // and psi is 0 in cases of singularity.
315 // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
316 // Since 0 <= theta <= 180, sin(theta) will be always
317 // non-negative. Therefore, it will never change the sign of both of
318 // the parameters passed to atan2.
319
320 if (fabs(stheta) < 1e-6) {
321 psi = 0.0;
322 phi = atan2(-this->data_[1][0], this->data_[0][0]);
323 }
324 // we only have one unique solution
325 else {
326 phi = atan2(this->data_[2][0], -this->data_[2][1]);
327 psi = atan2(this->data_[0][2], this->data_[1][2]);
328 }
329
330 // wrap phi and psi, make sure they are in the range from 0 to 2*Pi
331 if (phi < 0) phi += 2.0 * Constants::PI;
332
333 if (psi < 0) psi += 2.0 * Constants::PI;
334
335 myEuler[0] = phi;
336 myEuler[1] = theta;
337 myEuler[2] = psi;
338
339 return myEuler;
340 }
341
342 bool closeEnough(
343 const Real& a, const Real& b,
344 const Real& epsilon = std::numeric_limits<Real>::epsilon()) {
345 return (epsilon > std::abs(a - b));
346 }
347
348 Vector3<Real> toRPY() {
349 const Real PI = 3.14159265358979323846264;
350 // check for gimbal lock
351 if (closeEnough(this->data_[0][2], -1.0)) {
352 Real x = 0; // gimbal lock, value of x doesn't matter
353 Real y = PI / 2;
354 Real z = x + atan2(this->data_[1][0], this->data_[2][0]);
355 return Vector3<Real>(x, y, z);
356 } else if (closeEnough(this->data_[0][2], 1.0)) {
357 Real x = 0;
358 Real y = -PI / 2;
359 Real z = -x + atan2(-this->data_[1][0], -this->data_[2][0]);
360 return Vector3<Real>(x, y, z);
361 } else {
362 // two solutions exist
363 Real x1 = -asin(this->data_[0][2]);
364 Real x2 = PI - x1;
365
366 Real y1 =
367 atan2(this->data_[1][2] / cos(x1), this->data_[2][2] / cos(x1));
368 Real y2 =
369 atan2(this->data_[1][2] / cos(x2), this->data_[2][2] / cos(x2));
370
371 Real z1 =
372 atan2(this->data_[0][1] / cos(x1), this->data_[0][0] / cos(x1));
373 Real z2 =
374 atan2(this->data_[0][1] / cos(x2), this->data_[0][0] / cos(x2));
375
376 // choose one solution to return
377 // for example the "shortest" rotation
378 if ((std::abs(x1) + std::abs(y1) + std::abs(z1)) <=
379 (std::abs(x2) + std::abs(y2) + std::abs(z2))) {
380 return Vector3<Real>(x1, y1, z1);
381 } else {
382 return Vector3<Real>(x2, y2, z2);
383 }
384 }
385 }
386
387 Vector<Real, 6> toVoigtTensor() {
388 Vector<Real, 6> voigt;
389 voigt[0] = this->data_[0][0];
390 voigt[1] = this->data_[1][1];
391 voigt[2] = this->data_[2][2];
392 voigt[3] = 0.5 * (this->data_[1][2] + this->data_[2][1]);
393 voigt[4] = 0.5 * (this->data_[0][2] + this->data_[2][0]);
394 voigt[5] = 0.5 * (this->data_[0][1] + this->data_[1][0]);
395 return voigt;
396 }
397
398 /** Returns the determinant of this matrix. */
399 Real determinant() const {
400 Real x, y, z;
401
402 x = this->data_[0][0] * (this->data_[1][1] * this->data_[2][2] -
403 this->data_[1][2] * this->data_[2][1]);
404 y = this->data_[0][1] * (this->data_[1][2] * this->data_[2][0] -
405 this->data_[1][0] * this->data_[2][2]);
406 z = this->data_[0][2] * (this->data_[1][0] * this->data_[2][1] -
407 this->data_[1][1] * this->data_[2][0]);
408 return (x + y + z);
409 }
410
411 /** Returns the trace of this matrix. */
412 Real trace() const {
413 return this->data_[0][0] + this->data_[1][1] + this->data_[2][2];
414 }
415
416 /**
417 * Sets the value of this matrix to the inverse of itself.
418 * @note since this simple algorithm can be applied to invert a 3 by 3
419 * matrix, we hide the implementation of inverse in SquareMatrix
420 * class
421 */
424 RealType det = determinant();
425 m(0, 0) = this->data_[1][1] * this->data_[2][2] -
426 this->data_[1][2] * this->data_[2][1];
427 m(1, 0) = this->data_[1][2] * this->data_[2][0] -
428 this->data_[1][0] * this->data_[2][2];
429 m(2, 0) = this->data_[1][0] * this->data_[2][1] -
430 this->data_[1][1] * this->data_[2][0];
431 m(0, 1) = this->data_[2][1] * this->data_[0][2] -
432 this->data_[2][2] * this->data_[0][1];
433 m(1, 1) = this->data_[2][2] * this->data_[0][0] -
434 this->data_[2][0] * this->data_[0][2];
435 m(2, 1) = this->data_[2][0] * this->data_[0][1] -
436 this->data_[2][1] * this->data_[0][0];
437 m(0, 2) = this->data_[0][1] * this->data_[1][2] -
438 this->data_[0][2] * this->data_[1][1];
439 m(1, 2) = this->data_[0][2] * this->data_[1][0] -
440 this->data_[0][0] * this->data_[1][2];
441 m(2, 2) = this->data_[0][0] * this->data_[1][1] -
442 this->data_[0][1] * this->data_[1][0];
443 m /= det;
444 return m;
445 }
446
447 SquareMatrix3<Real> transpose() const {
448 SquareMatrix3<Real> result;
449
450 for (unsigned int i = 0; i < 3; i++)
451 for (unsigned int j = 0; j < 3; j++)
452 result(j, i) = this->data_[i][j];
453
454 return result;
455 }
456 /**
457 * Extract the eigenvalues and eigenvectors from a 3x3 matrix.
458 * The eigenvectors (the columns of V) will be normalized.
459 * The eigenvectors are aligned optimally with the x, y, and z
460 * axes respectively.
461 * @param a symmetric matrix whose eigenvectors are to be computed. On
462 * return, the matrix is overwritten
463 * @param w will contain the eigenvalues of the matrix On return of this
464 * function
465 * @param v the columns of this matrix will contain the eigenvectors. The
466 * eigenvectors are normalized and mutually orthogonal.
467 * @warning a will be overwritten
468 */
471 };
472 /*=========================================================================
473
474 Program: Visualization Toolkit
475 Module: $RCSfile: SquareMatrix3.hpp,v $
476
477 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
478 All rights reserved.
479 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
480
481 This software is distributed WITHOUT ANY WARRANTY; without even
482 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
483 PURPOSE. See the above copyright notice for more information.
484
485 =========================================================================*/
486 template<typename Real>
488 Vector3<Real>& w,
490 int i, j, k, maxI;
491 Real tmp, maxVal;
492 Vector3<Real> v_maxI, v_k, v_j;
493
494 // diagonalize using Jacobi
496 // if all the eigenvalues are the same, return identity matrix
497 if (w[0] == w[1] && w[0] == w[2]) {
499 return;
500 }
501
502 // transpose temporarily, it makes it easier to sort the eigenvectors
503 v = v.transpose();
504
505 // if two eigenvalues are the same, re-orthogonalize to optimally line
506 // up the eigenvectors with the x, y, and z axes
507 for (i = 0; i < 3; i++) {
508 if (w((i + 1) % 3) == w((i + 2) % 3)) { // two eigenvalues are the same
509 // find maximum element of the independant eigenvector
510 maxVal = fabs(v(i, 0));
511 maxI = 0;
512 for (j = 1; j < 3; j++) {
513 if (maxVal < (tmp = fabs(v(i, j)))) {
514 maxVal = tmp;
515 maxI = j;
516 }
517 }
518
519 // swap the eigenvector into its proper position
520 if (maxI != i) {
521 tmp = w(maxI);
522 w(maxI) = w(i);
523 w(i) = tmp;
524
525 v.swapRow(i, maxI);
526 }
527 // maximum element of eigenvector should be positive
528 if (v(maxI, maxI) < 0) {
529 v(maxI, 0) = -v(maxI, 0);
530 v(maxI, 1) = -v(maxI, 1);
531 v(maxI, 2) = -v(maxI, 2);
532 }
533
534 // re-orthogonalize the other two eigenvectors
535 j = (maxI + 1) % 3;
536 k = (maxI + 2) % 3;
537
538 v(j, 0) = 0.0;
539 v(j, 1) = 0.0;
540 v(j, 2) = 0.0;
541 v(j, j) = 1.0;
542
543 /** @todo */
544 v_maxI = v.getRow(maxI);
545 v_j = v.getRow(j);
546 v_k = cross(v_maxI, v_j);
547 v_k.normalize();
548 v_j = cross(v_k, v_maxI);
549 v.setRow(j, v_j);
550 v.setRow(k, v_k);
551
552 // transpose vectors back to columns
553 v = v.transpose();
554 return;
555 }
556 }
557
558 // the three eigenvalues are different, just sort the eigenvectors
559 // to align them with the x, y, and z axes
560
561 // find the vector with the largest x element, make that vector
562 // the first vector
563 maxVal = fabs(v(0, 0));
564 maxI = 0;
565 for (i = 1; i < 3; i++) {
566 if (maxVal < (tmp = fabs(v(i, 0)))) {
567 maxVal = tmp;
568 maxI = i;
569 }
570 }
571
572 // swap eigenvalue and eigenvector
573 if (maxI != 0) {
574 tmp = w(maxI);
575 w(maxI) = w(0);
576 w(0) = tmp;
577 v.swapRow(maxI, 0);
578 }
579 // do the same for the y element
580 if (fabs(v(1, 1)) < fabs(v(2, 1))) {
581 tmp = w(2);
582 w(2) = w(1);
583 w(1) = tmp;
584 v.swapRow(2, 1);
585 }
586
587 // ensure that the sign of the eigenvectors is correct
588 for (i = 0; i < 2; i++) {
589 if (v(i, i) < 0) {
590 v(i, 0) = -v(i, 0);
591 v(i, 1) = -v(i, 1);
592 v(i, 2) = -v(i, 2);
593 }
594 }
595
596 // set sign of final eigenvector to ensure that determinant is positive
597 if (v.determinant() < 0) {
598 v(2, 0) = -v(2, 0);
599 v(2, 1) = -v(2, 1);
600 v(2, 2) = -v(2, 2);
601 }
602
603 // transpose the eigenvectors back again
604 v = v.transpose();
605 return;
606 }
607
608 /**
609 * Return the multiplication of two matrixes (m1 * m2).
610 * @return the multiplication of two matrixes
611 * @param m1 the first matrix
612 * @param m2 the second matrix
613 */
614 template<typename Real>
616 const SquareMatrix3<Real>& m2) {
617 SquareMatrix3<Real> result;
618
619 for (unsigned int i = 0; i < 3; i++)
620 for (unsigned int j = 0; j < 3; j++)
621 for (unsigned int k = 0; k < 3; k++)
622 result(i, j) += m1(i, k) * m2(k, j);
623
624 return result;
625 }
626
627 template<typename Real>
628 inline SquareMatrix3<Real> outProduct(const Vector3<Real>& v1,
629 const Vector3<Real>& v2) {
630 SquareMatrix3<Real> result;
631
632 for (unsigned int i = 0; i < 3; i++) {
633 for (unsigned int j = 0; j < 3; j++) {
634 result(i, j) = v1[i] * v2[j];
635 }
636 }
637
638 return result;
639 }
640
641 using Mat3x3d = SquareMatrix3<RealType>;
642 using RotMat3x3d = SquareMatrix3<RealType>;
643
644 const Mat3x3d M3Zero(0.0);
645} // namespace OpenMD
646
647#endif // MATH_SQUAREMATRIX3_HPP
Quaternion is a sort of a higher-level complex number.
Real x() const
Returns the value of the first element of this quaternion.
Real z() const
Returns the value of the fourth element of this quaternion.
SquareMatrix< Real, 3 > toRotationMatrix3()
Returns the corresponding rotation matrix (3x3)
Real w() const
Returns the value of the first element of this quaternion.
Real y() const
Returns the value of the thirf element of this quaternion.
Vector< Real, Row > getRow(unsigned int row)
Returns a row of this matrix as a vector.
void setRow(unsigned int row, const Vector< Real, Row > &v)
Sets a row of this matrix.
void swapRow(unsigned int i, unsigned int j)
swap two rows of this matrix
static void diagonalize(SquareMatrix3< Real > &a, Vector3< Real > &w, SquareMatrix3< Real > &v)
Extract the eigenvalues and eigenvectors from a 3x3 matrix.
SquareMatrix3< Real > inverse() const
Sets the value of this matrix to the inverse of itself.
SquareMatrix3(Real s)
Constructs and initializes every element of this matrix to a scalar.
Real determinant() const
Returns the determinant of this matrix.
SquareMatrix3()
default constructor
void axisAngle(Vector3d axis, RealType angle)
Uses Rodrigues' rotation formula for a rotation matrix.
SquareMatrix3< Real > & operator=(const SquareMatrix< Real, 3 > &m)
copy assignment operator
void setupRotMat(const Quaternion< Real > &quat)
Sets this matrix to a rotation matrix by quaternion.
SquareMatrix3(Real *array)
Constructs and initializes from an array.
Real trace() const
Returns the trace of this matrix.
Quaternion< Real > toQuaternion()
Returns the quaternion from this rotation matrix.
void setupUpperTriangularVoigtTensor(Vector< Real, 6 > vt)
Sets this matrix to an upper-triangular (asymmetric) tensor using Voigt Notation.
void setupRotMat(const Vector3< Real > &eulerAngles)
Sets this matrix to a rotation matrix by three euler angles @ param euler.
SquareMatrix3(const SquareMatrix< Real, 3 > &m)
copy constructor
Vector3< Real > toEulerAngles()
Returns the euler angles from this rotation matrix.
void setupRotMat(Real w, Real x, Real y, Real z)
Sets this matrix to a rotation matrix by quaternion.
void setupRotMat(Real phi, Real theta, Real psi)
Sets this matrix to a rotation matrix by three euler angles.
void setupVoigtTensor(Vector< Real, 6 > vt)
Sets this matrix to a symmetric tensor using Voigt Notation.
A square matrix class.
SquareMatrix< Real, Dim > & operator=(const RectMatrix< Real, Dim, Dim > &m)
copy assignment operator
static SquareMatrix< Real, Dim > identity()
static int jacobi(SquareMatrix< Real, Dim > &a, Vector< Real, Dim > &d, SquareMatrix< Real, Dim > &v)
Fix length vector class.
Definition Vector.hpp:78
void normalize()
Normalizes this vector in place.
Definition Vector.hpp:402
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136
DynamicRectMatrix< Real > operator*(const DynamicRectMatrix< Real > &m, Real s)
Return the multiplication of scalar and matrix (m * s).