16constexpr double pi() {
return 3.1415926535897932384626433832795029; }
20constexpr double hc() {
return 12398.4197386209; }
31constexpr double deg(
double angle) {
return 180.0 /
pi() * angle; }
32constexpr double rad(
double angle) {
return pi() / 180.0 * angle; }
34constexpr float sq(
float x) {
return x * x; }
35constexpr double sq(
double x) {
return x * x; }
37inline int iround(
double d) {
return static_cast<int>(std::round(
d)); }
40 double d = std::fabs(a - b);
43 return std::min(
d,
full -
d);
46template <
typename Real>
50 Vec3_() : x(0), y(0), z(0) {}
51 Vec3_(Real x_, Real y_, Real z_) : x(x_), y(y_), z(z_) {}
52 explicit Vec3_(std::array<int, 3> h) : x(h[0]), y(h[1]), z(h[2]) {}
59 default:
throw std::out_of_range(
"Vec3 index must be 0, 1 or 2.");
62 Real at(
int i)
const {
return const_cast<Vec3_*
>(
this)->at(i); }
64 Vec3_ operator-()
const {
return {-x, -y, -z}; }
65 Vec3_ operator-(
const Vec3_& o)
const {
return {x-o.x, y-o.y, z-o.z}; }
66 Vec3_ operator+(
const Vec3_& o)
const {
return {x+o.x, y+o.y, z+o.z}; }
67 Vec3_
operator*(Real d)
const {
return {x*d, y*d, z*d}; }
68 Vec3_ operator/(Real d)
const {
return *
this * (1.0/d); }
69 Vec3_& operator-=(
const Vec3_& o) { *
this = *
this - o;
return *
this; }
70 Vec3_& operator+=(
const Vec3_& o) { *
this = *
this + o;
return *
this; }
71 Vec3_&
operator*=(Real d) { *
this = *
this * d;
return *
this; }
72 Vec3_& operator/=(Real d) {
return operator*=(1.0/d); }
74 Vec3_ negated()
const {
return {-x, -y, -z}; }
75 Real dot(
const Vec3_& o)
const {
return x*o.x + y*o.y + z*o.z; }
76 Vec3_ cross(
const Vec3_& o)
const {
77 return {y*o.z - z*o.y, z*o.x - x*o.z, x*o.y - y*o.x};
79 Real length_sq()
const {
return x * x + y * y + z * z; }
80 Real length()
const {
return std::sqrt(length_sq()); }
81 Vec3_ changed_magnitude(Real m)
const {
return operator*(m / length()); }
82 Vec3_ normalized()
const {
return changed_magnitude(1.0); }
83 Real dist_sq(
const Vec3_& o)
const {
return (*
this - o).length_sq(); }
84 Real dist(
const Vec3_& o)
const {
return std::sqrt(dist_sq(o)); }
85 Real cos_angle(
const Vec3_& o)
const {
86 return dot(o) / std::sqrt(length_sq() * o.length_sq());
88 Real angle(
const Vec3_& o)
const {
89 return std::acos(std::max(-1., std::min(1., cos_angle(o))));
91 bool approx(
const Vec3_& o, Real epsilon)
const {
92 return std::fabs(x - o.x) <= epsilon &&
93 std::fabs(y - o.y) <= epsilon &&
94 std::fabs(z - o.z) <= epsilon;
96 bool has_nan()
const {
97 return std::isnan(x) || std::isnan(y) || std::isnan(z);
116 double a[3][3] = { {1.,0.,0.}, {0.,1.,0.}, {0.,0.,1.} };
126 double c1,
double c2,
double c3)
135 throw std::out_of_range(
"Mat33 row index must be 0, 1 or 2.");
136 return Vec3(a[
i][0], a[
i][1], a[
i][2]);
141 throw std::out_of_range(
"Mat33 column index must be 0, 1 or 2.");
142 return Vec3(a[0][
i], a[1][
i], a[2][
i]);
146 return Mat33(a[0][0] + b[0][0], a[0][1] + b[0][1], a[0][2] + b[0][2],
147 a[1][0] + b[1][0], a[1][1] + b[1][1], a[1][2] + b[1][2],
148 a[2][0] + b[2][0], a[2][1] + b[2][1], a[2][2] + b[2][2]);
151 return Mat33(a[0][0] - b[0][0], a[0][1] - b[0][1], a[0][2] - b[0][2],
152 a[1][0] - b[1][0], a[1][1] - b[1][1], a[1][2] - b[1][2],
153 a[2][0] - b[2][0], a[2][1] - b[2][1], a[2][2] - b[2][2]);
157 return {a[0][0] *
p.x + a[0][1] *
p.y + a[0][2] *
p.z,
158 a[1][0] *
p.x + a[1][1] *
p.y + a[1][2] *
p.z,
159 a[2][0] *
p.x + a[2][1] *
p.y + a[2][2] *
p.z};
162 return {a[0][0] *
p.x + a[1][0] *
p.y + a[2][0] *
p.z,
163 a[0][1] *
p.x + a[1][1] *
p.y + a[2][1] *
p.z,
164 a[0][2] *
p.x + a[1][2] *
p.y + a[2][2] *
p.z};
168 return Mat33(a[0][0] *
p.x, a[0][1] *
p.y, a[0][2] *
p.z,
169 a[1][0] *
p.x, a[1][1] *
p.y, a[1][2] *
p.z,
170 a[2][0] *
p.x, a[2][1] *
p.y, a[2][2] *
p.z);
174 for (
int i = 0;
i != 3; ++
i)
175 for (
int j = 0;
j != 3; ++
j)
176 r[
i][
j] = a[
i][0] * b[0][
j] + a[
i][1] * b[1][
j] + a[
i][2] * b[2][
j];
180 return Mat33(a[0][0], a[1][0], a[2][0],
181 a[0][1], a[1][1], a[2][1],
182 a[0][2], a[1][2], a[2][2]);
184 double trace()
const {
return a[0][0] + a[1][1] + a[2][2]; }
187 for (
int i = 0;
i < 3; ++
i)
188 for (
int j = 0;
j < 3; ++
j)
194 for (
int i = 0;
i < 3; ++
i)
195 for (
int j = 0;
j < 3; ++
j)
196 if (std::isnan(a[
i][
j]))
202 return a[0][0] * (a[1][1]*a[2][2] - a[2][1]*a[1][2]) +
203 a[0][1] * (a[1][2]*a[2][0] - a[2][2]*a[1][0]) +
204 a[0][2] * (a[1][0]*a[2][1] - a[2][0]*a[1][1]);
208 double inv_det = 1.0 / determinant();
209 inv[0][0] =
inv_det * (a[1][1] * a[2][2] - a[2][1] * a[1][2]);
210 inv[0][1] =
inv_det * (a[0][2] * a[2][1] - a[0][1] * a[2][2]);
211 inv[0][2] =
inv_det * (a[0][1] * a[1][2] - a[0][2] * a[1][1]);
212 inv[1][0] =
inv_det * (a[1][2] * a[2][0] - a[1][0] * a[2][2]);
213 inv[1][1] =
inv_det * (a[0][0] * a[2][2] - a[0][2] * a[2][0]);
214 inv[1][2] =
inv_det * (a[1][0] * a[0][2] - a[0][0] * a[1][2]);
215 inv[2][0] =
inv_det * (a[1][0] * a[2][1] - a[2][0] * a[1][1]);
216 inv[2][1] =
inv_det * (a[2][0] * a[0][1] - a[0][0] * a[2][1]);
217 inv[2][2] =
inv_det * (a[0][0] * a[1][1] - a[1][0] * a[0][1]);
221 return a[0][0] == 1 && a[0][1] == 0 && a[0][2] == 0 &&
222 a[1][0] == 0 && a[1][1] == 1 && a[1][2] == 0 &&
223 a[2][0] == 0 && a[2][1] == 0 && a[2][2] == 1;
227 return a[0][
i] * a[0][
j] + a[1][
i] * a[1][
j] + a[2][
i] * a[2][
j];
231 return a[1][0] == 0 && a[2][0] == 0 && a[2][1] == 0;
236 double a11 = 0, a12 = 0, a13 = 0;
237 double a22 = 0, a23 = 0;
241 if (
m.is_upper_triangular()) {
249 a11 = a12 = a13 = a22 = a23 = a33 =
NAN;
253 return {a11 *
p.x + a12 *
p.y + a13 *
p.z,
254 a22 *
p.y + a23 *
p.z,
261 T u11, u22, u33, u12, u13, u23;
265 std::array<T, 6>
elements_pdb()
const {
return {{u11, u22, u33, u12, u13, u23}}; }
266 std::array<T, 6>
elements_voigt()
const {
return {{u11, u22, u33, u23, u13, u12}}; }
269 return Mat33(u11, u12, u13, u12, u22, u23, u13, u23, u33);
272 T trace()
const {
return u11 + u22 + u33; }
276 return u11 == 0 && u22 == 0 && u33 == 0 && u12 == 0 && u13 == 0 && u23 == 0;
280 u11 *= s; u22 *= s; u33 *= s; u12 *= s; u13 *= s; u23 *= s;
283 template<
typename Real>
285 return SMat33<Real>{u11*s, u22*s, u33*s, u12*s, u13*s, u23*s};
290 return {u11+
k, u22+
k, u33+
k, u12, u13, u23};
294 template<
typename VT>
296 return r.x *
r.x * u11 +
r.y *
r.y * u22 +
r.z *
r.z * u33 +
297 2 * (
r.x *
r.y * u12 +
r.x *
r.z * u13 +
r.y *
r.z * u23);
299 double r_u_r(
const std::array<int,3>&
h)
const {
301 return r_u_r(
Vec3(
h));
305 return {u11 *
p.x + u12 *
p.y + u13 *
p.z,
306 u12 *
p.x + u22 *
p.y + u23 *
p.z,
307 u13 *
p.x + u23 *
p.y + u33 *
p.z};
311 return {u11-
o.u11, u22-
o.u22, u33-
o.u33, u12-
o.u12, u13-
o.u13, u23-
o.u23};
314 return {u11+
o.u11, u22+
o.u22, u33+
o.u33, u12+
o.u12, u13+
o.u13, u23+
o.u23};
318 template<
typename Real=
double>
321 auto elem = [&](
int i,
int j) {
322 return static_cast<Real>(
323 m[
i][0] * (
m[
j][0] * u11 +
m[
j][1] * u12 +
m[
j][2] * u13) +
324 m[
i][1] * (
m[
j][0] * u12 +
m[
j][1] * u22 +
m[
j][2] * u23) +
325 m[
i][2] * (
m[
j][0] * u13 +
m[
j][1] * u23 +
m[
j][2] * u33));
328 elem(0, 1), elem(0, 2), elem(1, 2)};
332 return u11 * (u22*u33 - u23*u23) +
333 u12 * (u23*u13 - u33*u12) +
334 u13 * (u12*u23 - u13*u22);
352 double p1 = u12*u12 + u13*u13 + u23*u23;
354 return {{u11, u22, u33}};
355 double q = (1./3.) * trace();
357 double p2 =
sq(b.u11) +
sq(b.u22) +
sq(b.u33) + 2 *
p1;
358 double p = std::sqrt((1./6.) *
p2);
359 double r = b.determinant() / ((1./3.) *
p2 *
p);
362 phi = (1./3.) *
pi();
364 phi = (1./3.) * std::acos(
r);
366 double eig3 =
q + 2 *
p * std::cos(
phi + 2./3.*
pi());
377 return {
minv,
minv.multiply(vec).negated()};
387 return mat.
is_identity() && vec.x == 0. && vec.y == 0. && vec.z == 0.;
392 return mat.
has_nan() || vec.has_nan();
400template<
typename Pos>
405 if (
p.x < minimum.x) minimum.x =
p.x;
406 if (
p.y < minimum.y) minimum.y =
p.y;
407 if (
p.z < minimum.z) minimum.z =
p.z;
408 if (
p.x > maximum.x) maximum.x =
p.x;
409 if (
p.y > maximum.y) maximum.y =
p.y;
410 if (
p.z > maximum.z) maximum.z =
p.z;
421typename std::enable_if<std::is_integral<T>::value,
bool>::type
422is_nan(T) {
return false; }
424typename std::enable_if<std::is_floating_point<T>::value,
bool>::type
425is_nan(T a) {
return std::isnan(a); }
428typename std::enable_if<std::is_integral<T>::value,
bool>::type
429is_same(T a, T b) {
return a == b; }
431typename std::enable_if<std::is_floating_point<T>::value,
bool>::type
432is_same(T a, T b) {
return std::isnan(b) ? std::isnan(a) : a == b; }
constexpr double deg(double angle)
constexpr double bohrradius()
Vec3 rotate_about_axis(const Vec3 &v, const Vec3 &axis, double theta)
Rodrigues' rotation formula: rotate vector v about given axis of rotation (which must be a unit vecto...
constexpr double mott_bethe_const()
Vec3 operator*(double d, const Vec3 &v)
Op & operator*=(Op &a, const Op &b)
double angle_abs_diff(double a, double b, double full=360.0)
constexpr double rad(double angle)
constexpr float sq(float x)
constexpr double u_to_b()
void extend(const Pos &p)
void add_margins(const Pos &p)
void add_margin(double m)
bool is_upper_triangular() const
Vec3 column_copy(int i) const
row_t & operator[](int i)
double determinant() const
Mat33 operator+(const Mat33 &b) const
Mat33 multiply(const Mat33 &b) const
Vec3 left_multiply(const Vec3 &p) const
Mat33(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
bool approx(const Mat33 &other, double epsilon) const
Mat33 multiply_by_diagonal(const Vec3 &p) const
const row_t & operator[](int i) const
Vec3 row_copy(int i) const
static Mat33 from_columns(const Vec3 &c1, const Vec3 &c2, const Vec3 &c3)
Vec3 multiply(const Vec3 &p) const
double column_dot(int i, int j) const
Mat33 operator-(const Mat33 &b) const
SMat33 operator+(const SMat33 &o) const
std::array< T, 6 > elements_voigt() const
SMat33 operator-(const SMat33 &o) const
Vec3 multiply(const Vec3 &p) const
SMat33< Real > scaled(Real s) const
SMat33< Real > transformed_by(const Mat33 &m) const
std::array< double, 3 > calculate_eigenvalues() const
Based on https://en.wikipedia.org/wiki/Eigenvalue_algorithm To calculate both eigenvalues and eigenve...
double r_u_r(const std::array< int, 3 > &h) const
auto r_u_r(const Vec3_< VT > &r) const -> decltype(r.x+u11)
SMat33< T > added_kI(T k) const
std::array< T, 6 > elements_pdb() const
Vec3 multiply(const Vec3 &p) const
UpperTriangularMat33()=default
void operator=(const Mat33 &m)