right aligning values in cif

CMakeLists.txt

@@ -138,13 +138,6 @@ if(MSVC)
     # make msvc standards compliant...
     add_compile_options(/permissive- /bigobj)
     add_link_options(/NODEFAULTLIB:library)
-
-    # This is dubious...
-    if(BUILD_SHARED_LIBS)
-        set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>DLL")
-    else()
-        set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>")
-    endif()
 endif()
 
 # Libraries
@@ -205,8 +198,6 @@ if(NOT PCRE2_FOUND)
     add_subdirectory(pcre2-simple)
 endif()
 
-find_package(glm REQUIRED)
-
 # Using Eigen3 is a bit of a thing. We don't want to build it completely since
 # we only need a couple of header files. Nothing special. But often, eigen3 is
 # already installed and then we prefer that.
@@ -369,8 +360,9 @@ else()
 endif()
 
 if(NOT STD_CHARCONV_COMPILING)
-     target_include_directories(cifpp PRIVATE ${FastFloat_INCLUDE_DIRS})
-	 target_compile_definitions(cifpp PRIVATE USE_FAST_FLOAT)
+    get_target_property(FF_INC_DIR FastFloat::fast_float INTERFACE_INCLUDE_DIRECTORIES)
+    target_include_directories(cifpp PRIVATE ${FF_INC_DIR})
+	target_compile_definitions(cifpp PRIVATE USE_FAST_FLOAT)
 endif()
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")

changelog

@@ -1,5 +1,8 @@
 Version 10.0.2
 - Fixed regression in reconstruction introduced in 10.0.1
+- Fixed symmetry operations
+- Added validation for pdbx_item_enumeration as well as
+  case-sensitive checks for enumerations when needed.
 
 Version 10.0.1
 - Fixed some regressions, like assigning to items.

docs/symmetry.rst

@@ -8,7 +8,7 @@ Points
 
 The most basic type in use is :cpp:type:`cif::point`. It can be thought of as a point in space with three coordinates, but it is also often used as a vector in 3d space. To keep the interface simple there's no separate vector type.
 
-Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`glm::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
+Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`cif::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
 
 Quaternions
 -----------
@@ -91,7 +91,7 @@ To give an idea how this works, here's a piece of code copied from one of the un
         auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
 
         // The distance between these symmetry atoms should be equal to the distance in the struct_conn record
-        assert(glm::distance(sa1, sa2) == dist);
+        assert(cif::distance(sa1, sa2) == dist);
 
         // And to show how you can obtain the closest symmetry copy of an atom near another one:
         // here we request the symmetry copy of p2 that lies closest to p1

include/cif++/cif++.hpp

@@ -40,6 +40,7 @@
 #include "cif++/gzio.hpp"
 #include "cif++/item.hpp"
 #include "cif++/iterator.hpp"
+#include "cif++/matrix.hpp"
 #include "cif++/model.hpp"
 #include "cif++/parser.hpp"
 #include "cif++/pdb.hpp"

include/cif++/matrix.hpp

@@ -26,4 +26,751 @@
 
 #pragma once
 
-#warning "Using this file is deprecated"
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <ostream>
+#include <type_traits>
+#include <vector>
+
+/**
+ * @file matrix.hpp
+ *
+ * Some basic matrix operations and classes to hold matrices.
+ *
+ * We're using expression templates for optimal performance.
+ *
+ */
+
+namespace cif
+{
+// --------------------------------------------------------------------
+// We're using expression templates here
+
+/**
+ * @brief Base for the matrix expression templates
+ * This all uses the Curiously recurring template pattern
+ *
+ * @tparam M The type of the derived class
+ */
+template <typename M>
+class matrix_expression // NOLINT(bugprone-crtp-constructor-accessibility)
+{
+  public:
+	[[nodiscard]] constexpr std::size_t dim_m() const { return static_cast<const M &>(*this).dim_m(); } ///< Return the size (dimension) in direction m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return static_cast<const M &>(*this).dim_n(); } ///< Return the size (dimension) in direction n
+
+	[[nodiscard]] constexpr bool empty() const { return dim_m() == 0 or dim_n() == 0; } ///< Convenient way to test for empty matrices
+
+	/** Return a reference to element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto &operator()(std::size_t i, std::size_t j)
+	{
+		return static_cast<M &>(*this).operator()(i, j);
+	}
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return static_cast<const M &>(*this).operator()(i, j);
+	}
+
+	/** Swap the contents of rows @a r1 and @a r2 */
+	void swap_row(std::size_t r1, std::size_t r2)
+	{
+		for (std::size_t c = 0; c < dim_m(); ++c)
+		{
+			auto v = operator()(r1, c);
+			operator()(r1, c) = operator()(r2, c);
+			operator()(r2, c) = v;
+		}
+	}
+
+	/** Swap the contents of columns @a c1 and @a c2 */
+	void swap_col(std::size_t c1, std::size_t c2)
+	{
+		for (std::size_t r = 0; r < dim_n(); ++r)
+		{
+			auto &a = operator()(r, c1);
+			auto &b = operator()(r, c2);
+			std::swap(a, b);
+		}
+	}
+
+	/** write the matrix @a m to std::ostream @a os */
+	friend std::ostream &operator<<(std::ostream &os, const matrix_expression &m)
+	{
+		os << '[';
+
+		for (std::size_t i = 0; i < m.dim_m(); ++i)
+		{
+			os << '[';
+
+			for (std::size_t j = 0; j < m.dim_n(); ++j)
+			{
+				os << m(i, j);
+				if (j + 1 < m.dim_n())
+					os << ", ";
+			}
+
+			if (i + 1 < m.dim_m())
+				os << ", ";
+
+			os << ']';
+		}
+
+		os << ']';
+
+		return os;
+	}
+
+	/// compare two matrices
+	template <typename M2>
+	constexpr bool operator==(const matrix_expression<M2> &m) const
+	{
+		bool same = false;
+		if (dim_m() == m.dim_m() and dim_n() == m.dim_n())
+		{
+			same = true;
+			for (std::size_t i = 0; same and i < m.dim_m(); ++i)
+			{
+				for (std::size_t j = 0; same and j < m.dim_n(); ++j)
+					same = operator()(i, j) == m(i, j);
+			}
+		}
+
+		return same;
+	}
+};
+
+// --------------------------------------------------------------------
+
+/**
+ * @brief Storage class implementation of matrix_expression.
+ *
+ * @tparam F The type of the stored values
+ *
+ * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+ * element m i,j is mapped to [i * n + j] and thus storage is row major
+ */
+
+template <typename F = float>
+class matrix : public matrix_expression<matrix<F>>
+{
+  public:
+	/** The value type */
+	using value_type = F;
+
+	/**
+	 * @brief Copy construct a new matrix object using @a m
+	 *
+	 * @tparam M2 Type of @a m
+	 * @param m The matrix expression to copy values from
+	 */
+	template <typename M2>
+	matrix(const matrix_expression<M2> &m)
+		: m_m(m.dim_m())
+		, m_n(m.dim_n())
+		, m_data(m_m * m_n)
+	{
+		for (std::size_t i = 0; i < m_m; ++i)
+		{
+			for (std::size_t j = 0; j < m_n; ++j)
+				operator()(i, j) = m(i, j);
+		}
+	}
+
+	/**
+	 * @brief Construct a new matrix object with dimension @a m and @a n
+	 * setting the values to @a v
+	 *
+	 * @param m Requested dimension M
+	 * @param n Requested dimension N
+	 * @param v Value to store in each element
+	 */
+	matrix(std::size_t m, std::size_t n, value_type v = 0)
+		: m_m(m)
+		, m_n(n)
+		, m_data(m_m * m_n)
+	{
+		std::fill(m_data.begin(), m_data.end(), v);
+	}
+
+	/** @cond */
+	matrix() = default;
+	matrix(matrix &&m) = default;
+	matrix(const matrix &m) = default;
+	matrix &operator=(matrix &&m) = default;
+	matrix &operator=(const matrix &m) = default;
+	/** @endcond */
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+	/** Return a reference to element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+  private:
+	std::size_t m_m = 0, m_n = 0;
+	std::vector<value_type> m_data;
+};
+
+// --------------------------------------------------------------------
+// special case, 3x3 matrix
+
+/**
+ * @brief Storage class implementation of matrix_expression
+ * with compile time fixed size.
+ *
+ * @tparam F The type of the stored values
+ *
+ * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+ * element m i,j is mapped to [i * n + j] and thus storage is row major
+ */
+
+template <typename F, std::size_t M, std::size_t N>
+class matrix_fixed : public matrix_expression<matrix_fixed<F, M, N>>
+{
+  public:
+	/** The value type */
+	using value_type = F;
+
+	/** The storage size */
+	static constexpr std::size_t kSize = M * N;
+
+	/** Copy constructor */
+	template <typename M2>
+	matrix_fixed(const M2 &m)
+	{
+		assert(M == m.dim_m() and N == m.dim_n());
+		for (std::size_t i = 0; i < M; ++i)
+		{
+			for (std::size_t j = 0; j < N; ++j)
+				operator()(i, j) = m(i, j);
+		}
+	}
+
+	/** default constructor */
+	matrix_fixed(value_type v = 0)
+	{
+		m_data.fill(v);
+	}
+
+	/** Alternate constructor taking an array of values to store */
+	matrix_fixed(const F (&v)[kSize])
+	{
+		fill(v, std::make_index_sequence<kSize>{});
+	}
+
+	/** @cond */
+	matrix_fixed(matrix_fixed &&m) = default;
+	matrix_fixed(const matrix_fixed &m) = default;
+	matrix_fixed &operator=(matrix_fixed &&m) = default;
+	matrix_fixed &operator=(const matrix_fixed &m) = default;
+	/** @endcond */
+
+	/** Store the values in @a a in the matrix */
+	template <std::size_t... Ixs>
+	matrix_fixed &fill(const F (&a)[kSize], std::index_sequence<Ixs...>)
+	{
+		m_data = { a[Ixs]... };
+		return *this;
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return M; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return N; } ///< Return dimension n
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
+	{
+		assert(i < M);
+		assert(j < N);
+		return m_data[i * N + j];
+	}
+
+	/** Return a reference to element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
+	{
+		assert(i < M);
+		assert(j < N);
+		return m_data[i * N + j];
+	}
+
+  private:
+	std::array<value_type, M * N> m_data;
+};
+
+/** typedef of a fixed matrix of size 3x3 */
+template <typename F>
+using matrix3x3 = matrix_fixed<F, 3, 3>;
+
+/** typedef of a fixed matrix of size 4x4 */
+template <typename F>
+using matrix4x4 = matrix_fixed<F, 4, 4>;
+
+// --------------------------------------------------------------------
+
+/**
+ * @brief Storage class implementation of symmetric matrix_expression
+ *
+ * @tparam F The type of the stored values
+ *
+ * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+ * element m i,j is mapped to [i * n + j] and thus storage is row major
+ */
+template <typename F = float>
+class symmetric_matrix : public matrix_expression<symmetric_matrix<F>>
+{
+  public:
+	/** The value type */
+	using value_type = F;
+
+	/** constructor for a matrix of size @a n x @a n elements with value @a v */
+	symmetric_matrix(std::size_t n, value_type v = 0)
+		: m_n(n)
+		, m_data((m_n * (m_n + 1)) / 2)
+	{
+		std::fill(m_data.begin(), m_data.end(), v);
+	}
+
+	/** @cond */
+	symmetric_matrix() = default;
+	symmetric_matrix(symmetric_matrix &&m) = default;
+	symmetric_matrix(const symmetric_matrix &m) = default;
+	symmetric_matrix &operator=(symmetric_matrix &&m) = default;
+	symmetric_matrix &operator=(const symmetric_matrix &m) = default;
+	/** @endcond */
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_n; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
+	{
+		return i < j
+		           ? m_data[(j * (j + 1)) / 2 + i]
+		           : m_data[(i * (i + 1)) / 2 + j];
+	}
+
+	/** Return a reference to element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
+	{
+		if (i > j)
+			std::swap(i, j);
+		assert(j < m_n);
+		return m_data[(j * (j + 1)) / 2 + i];
+	}
+
+  private:
+	std::size_t m_n;
+	std::vector<value_type> m_data;
+};
+
+// --------------------------------------------------------------------
+
+/**
+ * @brief Storage class implementation of symmetric matrix_expression
+ * with compile time fixed size.
+ *
+ * @tparam F The type of the stored values
+ *
+ * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+ * element m i,j is mapped to [i * n + j] and thus storage is row major
+ */
+template <typename F, std::size_t M>
+class symmetric_matrix_fixed : public matrix_expression<symmetric_matrix_fixed<F, M>>
+{
+  public:
+	/** The value type */
+	using value_type = F;
+
+	/** constructor with all elements set to value @a v */
+	symmetric_matrix_fixed(value_type v = 0)
+	{
+		std::fill(m_data.begin(), m_data.end(), v);
+	}
+
+	/** @cond */
+	symmetric_matrix_fixed(symmetric_matrix_fixed &&m) = default;
+	symmetric_matrix_fixed(const symmetric_matrix_fixed &m) = default;
+	symmetric_matrix_fixed &operator=(symmetric_matrix_fixed &&m) = default;
+	symmetric_matrix_fixed &operator=(const symmetric_matrix_fixed &m) = default;
+	/** @endcond */
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return M; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return M; } ///< Return dimension n
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
+	{
+		return i < j
+		           ? m_data[(j * (j + 1)) / 2 + i]
+		           : m_data[(i * (i + 1)) / 2 + j];
+	}
+
+	/** Return a reference to element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
+	{
+		if (i > j)
+			std::swap(i, j);
+		assert(j < M);
+		return m_data[(j * (j + 1)) / 2 + i];
+	}
+
+  private:
+	std::array<value_type, (M * (M + 1)) / 2> m_data;
+};
+
+/** typedef of a fixed symmetric matrix of size 3x3 */
+template <typename F>
+using symmetric_matrix3x3 = symmetric_matrix_fixed<F, 3>;
+
+/** typedef of a fixed symmetric matrix of size 4x4 */
+template <typename F>
+using symmetric_matrix4x4 = symmetric_matrix_fixed<F, 4>;
+
+// --------------------------------------------------------------------
+
+/// A transposed matrix view
+
+template <typename M>
+class transposed_matrix : public cif::matrix_expression<transposed_matrix<M>>
+{
+  public:
+	transposed_matrix(const M &m)
+		: m_m(m)
+	{
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_n(); } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_m(); } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(j, i);
+	}
+
+  private:
+	const M &m_m;
+};
+
+// --------------------------------------------------------------------
+
+/**
+ * @brief implementation of symmetric matrix_expression with a value
+ * of 1 for the diagonal values and 0 for all the others.
+ *
+ * @tparam F The type of the stored values
+ *
+ * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+ * element m i,j is mapped to [i * n + j] and thus storage is row major
+ */
+template <typename F = float>
+class identity_matrix : public matrix_expression<identity_matrix<F>>
+{
+  public:
+	/** the value type */
+	using value_type = F;
+
+	/** constructor taking a dimension @a n */
+	identity_matrix(std::size_t n)
+		: m_n(n)
+	{
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_n; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
+
+	/** Return the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
+	{
+		return static_cast<value_type>(i == j ? 1 : 0);
+	}
+
+  private:
+	std::size_t m_n;
+};
+
+// --------------------------------------------------------------------
+// matrix functions, implemented as expression templates
+
+/**
+ * @brief Implementation of a substraction operation as a matrix expression
+ *
+ * @tparam M1 Type of matrix 1
+ * @tparam M2 Type of matrix 2
+ */
+template <typename M1, typename M2>
+class matrix_subtraction : public matrix_expression<matrix_subtraction<M1, M2>>
+{
+  public:
+	/** constructor */
+	matrix_subtraction(const M1 &m1, const M2 &m2)
+		: m_m1(m1)
+		, m_m2(m2)
+	{
+		assert(m_m1.dim_m() == m_m2.dim_m());
+		assert(m_m1.dim_n() == m_m2.dim_n());
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m1.dim_m(); } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m1.dim_n(); } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m1(i, j) - m_m2(i, j);
+	}
+
+  private:
+	const M1 &m_m1;
+	const M2 &m_m2;
+};
+
+/** operator to subtract two matrices and return a matrix expression */
+template <typename M1, typename M2>
+auto operator-(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
+{
+	return matrix_subtraction(m1, m2);
+}
+
+/**
+ * @brief Implementation of a multiplication operation as a matrix expression
+ *
+ * @tparam M1 Type of matrix 1
+ * @tparam M2 Type of matrix 2
+ */
+template <typename M1, typename M2>
+class matrix_matrix_multiplication : public matrix_expression<matrix_matrix_multiplication<M1, M2>>
+{
+  public:
+	/** constructor */
+	matrix_matrix_multiplication(const M1 &m1, const M2 &m2)
+		: m_m1(m1)
+		, m_m2(m2)
+	{
+		assert(m1.dim_n() == m2.dim_m());
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m1.dim_m(); } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m1.dim_n(); } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		using value_type = decltype(m_m1(0, 0));
+
+		value_type result = {};
+
+		for (std::size_t k = 0; k < m_m1.dim_n(); ++k)
+			result += m_m1(i, k) * m_m2(k, j);
+
+		return result;
+	}
+
+  private:
+	const M1 &m_m1;
+	const M2 &m_m2;
+};
+
+/**
+ * @brief Implementation of a multiplication operation of a matrix and a scalar value as a matrix expression
+ *
+ * @tparam M1 Type of matrix
+ * @tparam M2 Type of scalar value
+ */
+template <typename M, typename T>
+class matrix_scalar_multiplication : public matrix_expression<matrix_scalar_multiplication<M, T>>
+{
+  public:
+	/** value type */
+	using value_type = T;
+
+	/** constructor */
+	matrix_scalar_multiplication(const M &m, value_type v)
+		: m_m(m)
+		, m_v(v)
+	{
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m(); } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n(); } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(i, j) * m_v;
+	}
+
+  private:
+	const M &m_m;
+	value_type m_v;
+};
+
+/** First implementation of operator*, enabled if the second parameter is a scalar */
+template <typename M1, typename T>
+auto operator*(const matrix_expression<M1> &m, T v)
+	requires(std::is_floating_point_v<T>)
+{
+	return matrix_scalar_multiplication(m, v);
+}
+
+/** First implementation of operator*, enabled if the second parameter is not a scalar and thus must be a matrix, right? */
+template <typename M1, typename M2>
+auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
+	requires(not std::is_floating_point_v<M2>)
+{
+	return matrix_matrix_multiplication(m1, m2);
+}
+
+// --------------------------------------------------------------------
+
+/// A sub-view on a matrix
+template <typename M2>
+class sub_matrix : public matrix_expression<sub_matrix<M2>>
+{
+  public:
+	/// Constructor
+	sub_matrix(const M2 &m, int i, int j)
+		: m_m(m)
+		, m_i(i)
+		, m_j(j)
+	{
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(
+			i >= m_i ? i + 1 : i,
+			j >= m_j ? j + 1 : j);
+	}
+
+  private:
+	const M2 &m_m;
+	std::size_t m_i, m_j;
+};
+
+// --------------------------------------------------------------------
+
+/** Generic routine to calculate the determinant of a matrix
+ *
+ * @note This is currently only implemented for fixed matrices of size 3x3
+ */
+template <typename M>
+auto determinant(const M &m);
+
+/** Implementation of the determinant function for fixed size matrices of size 3x3 */
+template <typename F = float>
+auto determinant(const matrix3x3<F> &m)
+{
+	return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) +
+			m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) +
+			m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))));
+}
+
+/** Implementation of the determinant function for fixed size matrices of size 4x4 */
+template <typename F = float>
+F determinant(const matrix4x4<F> &m)
+{
+	return m(0, 0) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 0))) -
+	       m(0, 1) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 1))) +
+	       m(0, 2) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 2))) -
+	       m(0, 3) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 3)));
+}
+
+// --------------------------------------------------------------------
+
+/** Generic routine to calculate the inverse of a matrix
+ *
+ * @note This is currently only implemented for fixed matrices of size 3x3
+ */
+template <typename M>
+M inverse(const M &m);
+
+/** Implementation of the inverse function for fixed size matrices of size 3x3 */
+template <typename F = float>
+matrix3x3<F> inverse(const matrix3x3<F> &m)
+{
+	F det = determinant(m);
+
+	matrix3x3<F> result;
+
+	result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) / det;
+	result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) / det;
+	result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) / det;
+	result(0, 1) = (m(2, 1) * m(0, 2) - m(2, 2) * m(0, 1)) / det;
+	result(1, 1) = (m(2, 2) * m(0, 0) - m(2, 0) * m(0, 2)) / det;
+	result(2, 1) = (m(2, 0) * m(0, 1) - m(2, 1) * m(0, 0)) / det;
+	result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) / det;
+	result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) / det;
+	result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) / det;
+
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+/**
+ * @brief Implementation of a cofactor calculation as a matrix expression
+ *
+ * @tparam M Type of matrix
+ */
+template <typename M>
+class matrix_cofactors : public matrix_expression<matrix_cofactors<M>>
+{
+  public:
+	/** constructor */
+	matrix_cofactors(const M &m)
+		: m_m(m)
+	{
+	}
+
+	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m(); } ///< Return dimension m
+	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n(); } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		const std::size_t ixs[4][3] = {
+			{ 1, 2, 3 },
+			{ 0, 2, 3 },
+			{ 0, 1, 3 },
+			{ 0, 1, 2 }
+		};
+
+		const std::size_t *ix = ixs[i];
+		const std::size_t *iy = ixs[j];
+
+		auto result =
+			m_m(ix[0], iy[0]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[2]) +
+			m_m(ix[0], iy[1]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[0]) +
+			m_m(ix[0], iy[2]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[1]) -
+			m_m(ix[0], iy[2]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[0]) -
+			m_m(ix[0], iy[1]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[2]) -
+			m_m(ix[0], iy[0]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[1]);
+
+		return (i + j) % 2 == 1 ? -result : result;
+	}
+
+  private:
+	const M &m_m;
+};
+
+} // namespace cif

include/cif++/model.hpp

@@ -99,7 +99,7 @@ class atom
 		{
 			auto r = row();
 			if (r)
-				std::tie(m_location.x, m_location.y, m_location.z) = r.get<float, float, float>("Cartn_x", "Cartn_y", "Cartn_z");
+				std::tie(m_location.m_x, m_location.m_y, m_location.m_z) = r.get<float, float, float>("Cartn_x", "Cartn_y", "Cartn_z");
 		}
 
 		// constructor for a symmetry copy of an atom
@@ -298,16 +298,16 @@ class atom
 	/// \brief Rotate the position of this atom by \a q
 	void rotate(quaternion q)
 	{
-		set_location(q * get_location());
+		auto loc = get_location();
+		loc.rotate(q);
+		set_location(loc);
 	}
 
 	/// \brief rotate the coordinates of this atom by \a q around point \a p
 	void rotate(quaternion q, point p)
 	{
 		auto loc = get_location();
-		loc -= p;
-		loc = q * loc;
-		loc += p;
+		loc.rotate(q, p);
 		set_location(loc);
 	}
 
@@ -316,7 +316,7 @@ class atom
 	{
 		auto loc = get_location();
 		loc += t;
-		loc = q * loc;
+		loc.rotate(q);
 		set_location(loc);
 	}
 
@@ -325,7 +325,7 @@ class atom
 	{
 		auto loc = get_location();
 		loc += t1;
-		loc = q * loc;
+		loc.rotate(q);
 		loc += t2;
 		set_location(loc);
 	}
@@ -453,6 +453,16 @@ inline float distance(const atom &a, const atom &b)
 	return distance(a.get_location(), b.get_location());
 }
 
+/** Calculate the square of the distance between atoms @a and @a b in ångström
+ *
+ * @note Use this whenever possible instead of simply using distance since
+ * this function does not have to calculate a square root which is expensive.
+ */
+inline float distance_squared(const atom &a, const atom &b)
+{
+	return distance_squared(a.get_location(), b.get_location());
+}
+
 // --------------------------------------------------------------------
 
 /**

include/cif++/point.hpp

@@ -26,19 +26,21 @@
 
 #pragma once
 
+#include <array>
 #include <cmath>
+#include <complex>
+#include <cstdint>
 #include <cstdlib>
 #include <format>
-#include <glm/ext/quaternion_common.hpp>
-#include <glm/ext/quaternion_geometric.hpp>
-#include <glm/ext/quaternion_trigonometric.hpp>
-#include <glm/glm.hpp>
-#include <glm/gtc/quaternion.hpp>
-#include <glm/trigonometric.hpp>
+#include <functional>
 #include <limits>
 #include <numbers>
 #include <optional>
+#include <ostream>
 #include <tuple>
+#include <type_traits>
+#include <utility>
+#include <valarray>
 #include <vector>
 
 #if __has_include(<clipper/core/coords.h>)
@@ -52,36 +54,659 @@
 
 /** \file point.hpp
  *
- * This file contains the definition for *point* as well as
+ * This file contains the definition for *cif::point* as well as
  * lots of routines and classes that can manipulate points.
  */
 
 namespace cif
 {
 
-// Using glm now
+// --------------------------------------------------------------------
+/**
+ * @brief A stripped down quaternion implementation, based on boost::math::quaternion
+ *
+ * We use quaternions to do rotations in 3d space. Quaternions are faster than
+ * matrix calculations and they also suffer less from drift caused by rounding
+ * errors.
+ *
+ * Like complex number, quaternions do have a meaningful notion of "real part",
+ * but unlike them there is no meaningful notion of "imaginary part".
+ * Instead there is an "unreal part" which itself is a quaternion, and usually
+ * nothing simpler (as opposed to the complex number case).
+ * However, for practicality, there are accessors for the other components
+ * (these are necessary for the templated copy constructor, for instance).
+ *
+ * @note Quaternion multiplication is *NOT* commutative;
+ * symbolically, "q *= rhs;" means "q = q * rhs;"
+ * and "q /= rhs;" means "q = q * inverse_of(rhs);"
+ */
 
 template <typename T>
-using quaternion_type = glm::qua<T>;
+class quaternion_type
+{
+  public:
+	/// \brief the value type of the elements, usually this is float
+	using value_type = T;
 
+	/// \brief constructor with the four members
+	constexpr explicit quaternion_type(value_type const &value_a = {}, value_type const &value_b = {}, value_type const &value_c = {}, value_type const &value_d = {})
+		: a(value_a)
+		, b(value_b)
+		, c(value_c)
+		, d(value_d)
+	{
+	}
+
+	/// \brief constructor taking two complex values as input
+	constexpr explicit quaternion_type(std::complex<value_type> const &z0, std::complex<value_type> const &z1 = std::complex<value_type>())
+		: a(z0.real())
+		, b(z0.imag())
+		, c(z1.real())
+		, d(z1.imag())
+	{
+	}
+
+	constexpr quaternion_type(quaternion_type const &) = default; ///< Copy constructor
+	constexpr quaternion_type(quaternion_type &&) = default;      ///< Copy constructor
+
+	/// \brief Copy constructor accepting a quaternion with a different value_type
+	template <typename X>
+	constexpr explicit quaternion_type(quaternion_type<X> const &rhs)
+		: a(static_cast<value_type>(rhs.a))
+		, b(static_cast<value_type>(rhs.b))
+		, c(static_cast<value_type>(rhs.c))
+		, d(static_cast<value_type>(rhs.d))
+	{
+	}
+
+	// accessors
+
+	/// \brief See class description, return the *real* part of the quaternion
+	[[nodiscard]] constexpr value_type real() const
+	{
+		return a;
+	}
+
+	/// \brief See class description, return the *unreal* part of the quaternion
+	[[nodiscard]] constexpr quaternion_type unreal() const
+	{
+		return { 0, b, c, d };
+	}
+
+	/// \brief swap
+	constexpr void swap(quaternion_type &o)
+	{
+		std::swap(a, o.a);
+		std::swap(b, o.b);
+		std::swap(c, o.c);
+		std::swap(d, o.d);
+	}
+
+	// assignment operators
+
+	/// \brief Assignment operator accepting a quaternion with optionally another value_type
+	template <typename X>
+	constexpr quaternion_type &operator=(quaternion_type<X> const &rhs)
+	{
+		a = static_cast<value_type>(rhs.a);
+		b = static_cast<value_type>(rhs.b);
+		c = static_cast<value_type>(rhs.c);
+		d = static_cast<value_type>(rhs.d);
+
+		return *this;
+	}
+
+	/// \brief Assignment operator
+	constexpr quaternion_type &operator=(quaternion_type const &rhs) = default;
+
+	/// \brief Assignment operator that sets the *real* part to @a rhs and the *unreal* parts to zero
+	constexpr quaternion_type &operator=(value_type const &rhs)
+	{
+		a = rhs;
+
+		b = c = d = static_cast<value_type>(0);
+
+		return *this;
+	}
+
+	/// \brief Assignment operator that sets the *real* part to the real part of @a rhs
+	/// and the first *unreal* part to the imaginary part of of @a rhs. The other *unreal*
+	// parts are set to zero.
+	constexpr quaternion_type &operator=(std::complex<value_type> const &rhs)
+	{
+		a = rhs.real();
+		b = rhs.imag();
+
+		c = d = static_cast<value_type>(0);
+
+		return *this;
+	}
+
+	// other assignment-related operators
+
+	/// \brief operator += adding value @a rhs to the *real* part
+	constexpr quaternion_type &operator+=(value_type const &rhs)
+	{
+		a += rhs;
+		return *this;
+	}
+
+	/// \brief operator += adding the real part of @a rhs to the *real* part
+	/// and the imaginary part of @a rhs to the first *unreal* part
+	constexpr quaternion_type &operator+=(std::complex<value_type> const &rhs)
+	{
+		a += std::real(rhs);
+		b += std::imag(rhs);
+		return *this;
+	}
+
+	/// \brief operator += adding the parts of @a rhs to the equivalent part of this
+	template <class X>
+	constexpr quaternion_type &operator+=(quaternion_type<X> const &rhs)
+	{
+		a += rhs.a;
+		b += rhs.b;
+		c += rhs.c;
+		d += rhs.d;
+		return *this;
+	}
+
+	/// \brief operator -= subtracting value @a rhs from the *real* part
+	constexpr quaternion_type &operator-=(value_type const &rhs)
+	{
+		a -= rhs;
+		return *this;
+	}
+
+	/// \brief operator -= subtracting the real part of @a rhs from the *real* part
+	/// and the imaginary part of @a rhs from the first *unreal* part
+	constexpr quaternion_type &operator-=(std::complex<value_type> const &rhs)
+	{
+		a -= std::real(rhs);
+		b -= std::imag(rhs);
+		return *this;
+	}
+
+	/// \brief operator -= subtracting the parts of @a rhs from the equivalent part of this
+	template <class X>
+	constexpr quaternion_type &operator-=(quaternion_type<X> const &rhs)
+	{
+		a -= rhs.a;
+		b -= rhs.b;
+		c -= rhs.c;
+		d -= rhs.d;
+		return *this;
+	}
+
+	/// \brief multiply all parts with value @a rhs
+	constexpr quaternion_type &operator*=(value_type const &rhs)
+	{
+		a *= rhs;
+		b *= rhs;
+		c *= rhs;
+		d *= rhs;
+		return *this;
+	}
+
+	/// \brief multiply with complex number @a rhs
+	constexpr quaternion_type &operator*=(std::complex<value_type> const &rhs)
+	{
+		value_type ar = rhs.real();
+		value_type br = rhs.imag();
+		quaternion_type result(a * ar - b * br, a * br + b * ar, c * ar + d * br, -c * br + d * ar);
+		swap(result);
+		return *this;
+	}
+
+	/// \brief multiply @a a with @a b and return the result
+	friend constexpr quaternion_type operator*(const quaternion_type &a, const quaternion_type &b)
+	{
+		auto result = a;
+		result *= b;
+		return result;
+	}
+
+	/// \brief multiply with quaternion @a rhs
+	template <typename X>
+	constexpr quaternion_type &operator*=(quaternion_type<X> const &rhs)
+	{
+		auto ar = static_cast<value_type>(rhs.a);
+		auto br = static_cast<value_type>(rhs.b);
+		auto cr = static_cast<value_type>(rhs.c);
+		auto dr = static_cast<value_type>(rhs.d);
+
+		quaternion_type result(a * ar - b * br - c * cr - d * dr, a * br + b * ar + c * dr - d * cr, a * cr - b * dr + c * ar + d * br, a * dr + b * cr - c * br + d * ar);
+		swap(result);
+		return *this;
+	}
+
+	/// \brief divide all parts by @a rhs
+	constexpr quaternion_type &operator/=(value_type const &rhs)
+	{
+		a /= rhs;
+		b /= rhs;
+		c /= rhs;
+		d /= rhs;
+		return *this;
+	}
+
+	/// \brief divide by complex number @a rhs
+	constexpr quaternion_type &operator/=(std::complex<value_type> const &rhs)
+	{
+		value_type ar = rhs.real();
+		value_type br = rhs.imag();
+		value_type denominator = ar * ar + br * br;
+		quaternion_type result((+a * ar + b * br) / denominator, (-a * br + b * ar) / denominator, (+c * ar - d * br) / denominator, (+c * br + d * ar) / denominator);
+		swap(result);
+		return *this;
+	}
+
+	/// \brief divide by quaternion @a rhs
+	template <typename X>
+	constexpr quaternion_type &operator/=(quaternion_type<X> const &rhs)
+	{
+		auto ar = static_cast<value_type>(rhs.a);
+		auto br = static_cast<value_type>(rhs.b);
+		auto cr = static_cast<value_type>(rhs.c);
+		auto dr = static_cast<value_type>(rhs.d);
+
+		value_type denominator = ar * ar + br * br + cr * cr + dr * dr;
+		quaternion_type result((+a * ar + b * br + c * cr + d * dr) / denominator, (-a * br + b * ar - c * dr + d * cr) / denominator, (-a * cr + b * dr + c * ar - d * br) / denominator, (-a * dr - b * cr + c * br + d * ar) / denominator);
+		swap(result);
+		return *this;
+	}
+
+	/// \brief normalise the values so that the length of the result is exactly 1
+	friend constexpr quaternion_type normalize(quaternion_type q)
+	{
+		std::valarray<value_type> t(4);
+
+		t[0] = q.a;
+		t[1] = q.b;
+		t[2] = q.c;
+		t[3] = q.d;
+
+		t *= t;
+
+		value_type length = std::sqrt(t.sum());
+
+		if (length > 0.001)
+			q /= static_cast<value_type>(length);
+		else
+			q = quaternion_type(1, 0, 0, 0);
+
+		return q;
+	}
+
+	/// \brief return the conjugate of this
+	friend constexpr quaternion_type conj(quaternion_type q)
+	{
+		return quaternion_type{ +q.a, -q.b, -q.c, -q.d };
+	}
+
+	[[nodiscard]] constexpr value_type get_a() const { return a; } ///< Return part a
+	[[nodiscard]] constexpr value_type get_b() const { return b; } ///< Return part b
+	[[nodiscard]] constexpr value_type get_c() const { return c; } ///< Return part c
+	[[nodiscard]] constexpr value_type get_d() const { return d; } ///< Return part d
+
+	/// \brief compare with @a rhs
+	constexpr bool operator==(const quaternion_type &rhs) const
+	{
+		return a == rhs.a and b == rhs.b and c == rhs.c and d == rhs.d;
+	}
+
+	/// \brief compare with @a rhs
+	constexpr bool operator!=(const quaternion_type &rhs) const
+	{
+		return a != rhs.a or b != rhs.b or c != rhs.c or d != rhs.d;
+	}
+
+	/// \brief test for all zero values
+	constexpr explicit operator bool() const
+	{
+		return a != 0 or b != 0 or c != 0 or d != 0;
+	}
+
+	/// \brief for debugging e.g.
+	friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs)
+	{
+		os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d);
+		return os;
+	}
+
+  private:
+	value_type a, b, c, d;
+};
+
+/**
+ * @brief This code is similar to the code in boost so I copy the documentation as well:
+ *
+ * > spherical is a simple transposition of polar, it takes as inputs a (positive)
+ * > magnitude and a point on the hypersphere, given by three angles. The first of
+ * > these, theta has a natural range of -pi to +pi, and the other two have natural
+ * > ranges of -pi/2 to +pi/2 (as is the case with the usual spherical coordinates in
+ * > **R**<sup>3</sup>). Due to the many symmetries and periodicities, nothing untoward happens if
+ * > the magnitude is negative or the angles are outside their natural ranges. The
+ * > expected degeneracies (a magnitude of zero ignores the angles settings...) do
+ * > happen however.
+ */
+
+template <typename T>
+inline quaternion_type<T> spherical(T const &rho, T const &theta, T const &phi1, T const &phi2)
+{
+	T cos_phi1 = std::cos(phi1);
+	T cos_phi2 = std::cos(phi2);
+
+	T a = std::cos(theta) * cos_phi1 * cos_phi2;
+	T b = std::sin(theta) * cos_phi1 * cos_phi2;
+	T c = std::sin(phi1) * cos_phi2;
+	T d = std::sin(phi2);
+
+	quaternion_type result(a, b, c, d);
+	result *= rho;
+
+	return result;
+}
+
+/// \brief By default we use the float version of a quaternion
 using quaternion = quaternion_type<float>;
 
-template <typename T>
-using point_type = glm::vec<3, T>;
-
-using point = point_type<float>;
-
 // --------------------------------------------------------------------
 
-constexpr point rotate(point pt, quaternion q)
-{
-	return q * pt;
-}
+/**
+ * @brief 3D point: a location with x, y and z coordinates as floating point.
+ *
+ * Note that you can simply use structured binding to get access to the
+ * individual parts like so:
+ *
+ * @code{.cpp}
+ * float x, y, z;
+ * tie(x, y, z) = atom.get_location();
+ * @endcode
+ */
 
-constexpr point rotate(point pt, quaternion q, point a)
+template <typename F>
+struct point_type
 {
-	return q * (pt - a) + a;
-}
+	/// \brief the value type of the x, y and z members
+	using value_type = F;
+
+	value_type m_x, ///< The x part of the location
+		m_y,        ///< The y part of the location
+		m_z;        ///< The z part of the location
+
+	/// \brief default constructor, initialises the values to zero
+	constexpr point_type()
+		: m_x(0)
+		, m_y(0)
+		, m_z(0)
+	{
+	}
+
+	/// \brief constructor taking three values
+	constexpr point_type(value_type x, value_type y, value_type z)
+		: m_x(x)
+		, m_y(y)
+		, m_z(z)
+	{
+	}
+
+	/// \brief Copy constructor
+	template <typename PF>
+	constexpr point_type(const point_type<PF> &pt)
+		: m_x(static_cast<F>(pt.m_x))
+		, m_y(static_cast<F>(pt.m_y))
+		, m_z(static_cast<F>(pt.m_z))
+	{
+	}
+
+	/// \brief constructor taking a tuple of three values
+	constexpr point_type(const std::tuple<value_type, value_type, value_type> &pt)
+		: point_type(std::get<0>(pt), std::get<1>(pt), std::get<2>(pt))
+	{
+	}
+
+#if HAVE_LIBCLIPPER
+	/// \brief Construct a point using the values in clipper coordinate @a pt
+	constexpr point_type(const clipper::Coord_orth &pt)
+		: m_x(pt[0])
+		, m_y(pt[1])
+		, m_z(pt[2])
+	{
+	}
+
+	/// \brief Assign a point using the values in clipper coordinate @a rhs
+	constexpr point_type &operator=(const clipper::Coord_orth &rhs)
+	{
+		m_x = rhs[0];
+		m_y = rhs[1];
+		m_z = rhs[2];
+		return *this;
+	}
+#endif
+
+	/// \brief Assignment operator
+	template <typename PF>
+	constexpr point_type &operator=(const point_type<PF> &rhs)
+	{
+		m_x = static_cast<F>(rhs.m_x);
+		m_y = static_cast<F>(rhs.m_y);
+		m_z = static_cast<F>(rhs.m_z);
+		return *this;
+	}
+
+	[[nodiscard]] constexpr value_type &get_x() { return m_x; }      ///< Get a reference to x
+	[[nodiscard]] constexpr value_type get_x() const { return m_x; } ///< Get the value of x
+	constexpr void set_x(value_type x) { m_x = x; }                  ///< Set the value of x to @a x
+
+	[[nodiscard]] constexpr value_type &get_y() { return m_y; }      ///< Get a reference to y
+	[[nodiscard]] constexpr value_type get_y() const { return m_y; } ///< Get the value of y
+	constexpr void set_y(value_type y) { m_y = y; }                  ///< Set the value of y to @a y
+
+	[[nodiscard]] constexpr value_type &get_z() { return m_z; }      ///< Get a reference to z
+	[[nodiscard]] constexpr value_type get_z() const { return m_z; } ///< Get the value of z
+	constexpr void set_z(value_type z) { m_z = z; }                  ///< Set the value of z to @a z
+
+	/// \brief add @a rhs
+	constexpr point_type &operator+=(const point_type &rhs)
+	{
+		m_x += rhs.m_x;
+		m_y += rhs.m_y;
+		m_z += rhs.m_z;
+
+		return *this;
+	}
+
+	/// \brief add @a d to all members
+	constexpr point_type &operator+=(value_type d)
+	{
+		m_x += d;
+		m_y += d;
+		m_z += d;
+
+		return *this;
+	}
+
+	/// \brief Add the points @a lhs and @a rhs and return the result
+	template <typename F2>
+	friend constexpr auto operator+(const point_type &lhs, const point_type<F2> &rhs)
+	{
+		return point_type<std::common_type_t<value_type, F2>>(lhs.m_x + rhs.m_x, lhs.m_y + rhs.m_y, lhs.m_z + rhs.m_z);
+	}
+
+	/// \brief subtract @a rhs
+	constexpr point_type &operator-=(const point_type &rhs)
+	{
+		m_x -= rhs.m_x;
+		m_y -= rhs.m_y;
+		m_z -= rhs.m_z;
+
+		return *this;
+	}
+
+	/// \brief subtract @a d from all members
+	constexpr point_type &operator-=(value_type d)
+	{
+		m_x -= d;
+		m_y -= d;
+		m_z -= d;
+
+		return *this;
+	}
+
+	/// \brief Subtract the points @a lhs and @a rhs and return the result
+	template <typename F2>
+	friend constexpr auto operator-(const point_type &lhs, const point_type<F2> &rhs)
+	{
+		return point_type<std::common_type_t<value_type, F2>>(lhs.m_x - rhs.m_x, lhs.m_y - rhs.m_y, lhs.m_z - rhs.m_z);
+	}
+
+	/// \brief Return the negative copy of @a pt
+	friend constexpr point_type operator-(const point_type &pt)
+	{
+		return point_type(-pt.m_x, -pt.m_y, -pt.m_z);
+	}
+
+	/// \brief multiply all members with @a rhs
+	constexpr point_type &operator*=(value_type rhs)
+	{
+		m_x *= rhs;
+		m_y *= rhs;
+		m_z *= rhs;
+		return *this;
+	}
+
+	/// \brief multiply point @a pt with value @a f and return the result
+	template <typename F2>
+	friend constexpr auto operator*(const point_type &pt, F2 f)
+	{
+		return point_type<std::common_type_t<value_type, F2>>(pt.m_x * f, pt.m_y * f, pt.m_z * f);
+	}
+
+	/// \brief multiply point @a pt with value @a f and return the result
+	template <typename F2>
+	friend constexpr auto operator*(F2 f, const point_type &pt)
+	{
+		return point_type<std::common_type_t<value_type, F2>>(pt.m_x * f, pt.m_y * f, pt.m_z * f);
+	}
+
+	/// \brief divide all members by @a rhs
+	constexpr point_type &operator/=(value_type rhs)
+	{
+		m_x /= rhs;
+		m_y /= rhs;
+		m_z /= rhs;
+		return *this;
+	}
+
+	/// \brief divide point @a pt by value @a f and return the result
+	template <typename F2>
+	friend constexpr auto operator/(const point_type &pt, F2 f)
+	{
+		return point_type<std::common_type_t<value_type, F2>>(pt.m_x / f, pt.m_y / f, pt.m_z / f);
+	}
+
+	/**
+	 * @brief looking at this point as a vector, normalise it which
+	 * means dividing all members by the length making the length
+	 * effectively 1.
+	 *
+	 * @return The previous length of this vector
+	 */
+	constexpr value_type normalize()
+	{
+		auto length = m_x * m_x + m_y * m_y + m_z * m_z;
+		if (length > 0)
+		{
+			length = std::sqrt(length);
+			operator/=(length);
+		}
+		return length;
+	}
+
+	/// \brief Rotate this point using the quaterion @a q
+	constexpr void rotate(const quaternion &q)
+	{
+		quaternion_type<value_type> p(0, m_x, m_y, m_z);
+
+		p = q * p * conj(q);
+
+		m_x = p.get_b();
+		m_y = p.get_c();
+		m_z = p.get_d();
+	}
+
+	/// \brief Rotate this point using the quaterion @a q by first
+	/// moving the point to @a pivot and after rotating moving it
+	/// back
+	constexpr void rotate(const quaternion &q, point_type pivot)
+	{
+		operator-=(pivot);
+		rotate(q);
+		operator+=(pivot);
+	}
+
+#if HAVE_LIBCLIPPER
+	/// \brief Make it possible to pass a point to clipper functions expecting a clipper coordinate
+	operator clipper::Coord_orth() const
+	{
+		return clipper::Coord_orth(m_x, m_y, m_z);
+	}
+#endif
+
+	/// \brief Allow access to this point as if it is a tuple of three const value_type's
+	constexpr operator std::tuple<const value_type &, const value_type &, const value_type &>() const
+	{
+		return std::make_tuple(std::ref(m_x), std::ref(m_y), std::ref(m_z));
+	}
+
+	/// \brief Allow access to this point as if it is a tuple of three value_type's
+	constexpr operator std::tuple<value_type &, value_type &, value_type &>()
+	{
+		return std::make_tuple(std::ref(m_x), std::ref(m_y), std::ref(m_z));
+	}
+
+#if defined(__cpp_impl_three_way_comparison)
+	/// \brief a default spaceship operator
+	constexpr auto operator<=>(const point_type &rhs) const = default;
+#else
+	/// \brief a default equals operator
+	constexpr bool operator==(const point_type &rhs) const
+	{
+		return m_x == rhs.m_x and m_y == rhs.m_y and m_z == rhs.m_z;
+	}
+
+	/// \brief a default not-equals operator
+	constexpr bool operator!=(const point_type &rhs) const
+	{
+		return not operator==(rhs);
+	}
+#endif
+
+	// consider point as a vector... perhaps I should rename point?
+
+	/// \brief looking at the point as if it is a vector, return the squared length
+	[[nodiscard]] constexpr value_type length_sq() const
+	{
+		return m_x * m_x + m_y * m_y + m_z * m_z;
+	}
+
+	/// \brief looking at the point as if it is a vector, return the length
+	[[nodiscard]] constexpr value_type length() const
+	{
+		return std::sqrt(length_sq());
+	}
+
+	/// \brief Print out the point @a pt to @a os
+	friend std::ostream &operator<<(std::ostream &os, const point_type &pt)
+	{
+		os << '(' << pt.m_x << ',' << pt.m_y << ',' << pt.m_z << ')';
+		return os;
+	}
+};
+
+/// \brief By default we use points with float value_type
+using point = point_type<float>;
 
 // --------------------------------------------------------------------
 // several standard 3d operations
@@ -90,16 +715,43 @@ constexpr point rotate(point pt, quaternion q, point a)
 template <typename F1, typename F2>
 constexpr auto distance_squared(const point_type<F1> &a, const point_type<F2> &b)
 {
-	return (a.x - b.x) * (a.x - b.x) +
-	       (a.y - b.y) * (a.y - b.y) +
-	       (a.z - b.z) * (a.z - b.z);
+	return (a.m_x - b.m_x) * (a.m_x - b.m_x) +
+	       (a.m_y - b.m_y) * (a.m_y - b.m_y) +
+	       (a.m_z - b.m_z) * (a.m_z - b.m_z);
+}
+
+/// \brief return the distance between points @a a and @a b
+template <typename F1, typename F2>
+constexpr auto distance(const point_type<F1> &a, const point_type<F2> &b)
+{
+	return std::sqrt(
+		(a.m_x - b.m_x) * (a.m_x - b.m_x) +
+		(a.m_y - b.m_y) * (a.m_y - b.m_y) +
+		(a.m_z - b.m_z) * (a.m_z - b.m_z));
+}
+
+/// \brief return the dot product between the vectors @a a and @a b
+template <typename F1, typename F2>
+inline constexpr auto dot_product(const point_type<F1> &a, const point_type<F2> &b)
+{
+	return a.m_x * b.m_x + a.m_y * b.m_y + a.m_z * b.m_z;
+}
+
+/// \brief return the cross product between the vectors @a a and @a b
+template <typename F1, typename F2>
+inline constexpr auto cross_product(const point_type<F1> &a, const point_type<F2> &b)
+{
+	return point_type<std::common_type_t<F1, F2>>(
+		a.m_y * b.m_z - b.m_y * a.m_z,
+		a.m_z * b.m_x - b.m_z * a.m_x,
+		a.m_x * b.m_y - b.m_x * a.m_y);
 }
 
 /// \brief return the squared norm of point @a p
 template <typename F>
 constexpr F norm_squared(const point_type<F> &p)
 {
-	return p.x * p.x + p.y * p.y + p.z * p.z;
+	return p.m_x * p.m_x + p.m_y * p.m_y + p.m_z * p.m_z;
 }
 
 /// \brief return the norm of point @a p
@@ -111,23 +763,23 @@ constexpr point_type<F> norm(const point_type<F> &p)
 
 /// \brief return the point where two lines intersect, or an empty value if they don't intersect at all
 template <typename F>
-std::optional<point> line_line_intersection(const point_type<F> &p1,
+std::optional<cif::point> line_line_intersection(const point_type<F> &p1,
 	const point_type<F> &p2, const point_type<F> &p3, const point_type<F> &p4)
 {
 	auto p13 = p1 - p3;
 	auto p43 = p4 - p3;
-	if (std::abs(p43.x) < std::numeric_limits<F>::epsilon() and std::abs(p43.y) < std::numeric_limits<F>::epsilon() and std::abs(p43.z) < std::numeric_limits<F>::epsilon())
+	if (std::abs(p43.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_z) < std::numeric_limits<F>::epsilon())
 		return std::nullopt;
 
 	auto p21 = p2 - p1;
-	if (std::abs(p21.x) < std::numeric_limits<F>::epsilon() and std::abs(p21.y) < std::numeric_limits<F>::epsilon() and std::abs(p21.z) < std::numeric_limits<F>::epsilon())
+	if (std::abs(p21.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_z) < std::numeric_limits<F>::epsilon())
 		return std::nullopt;
 
-	auto d1343 = dot(p43, p13);
-	auto d4321 = dot(p43, p21);
-	auto d1321 = dot(p13, p21);
-	auto d4343 = dot(p43, p43);
-	auto d2121 = dot(p21, p21);
+	auto d1343 = cif::dot_product(p43, p13);
+	auto d4321 = cif::dot_product(p43, p21);
+	auto d1321 = cif::dot_product(p13, p21);
+	auto d4343 = cif::dot_product(p43, p43);
+	auto d2121 = cif::dot_product(p21, p21);
 
 	auto denom = d2121 * d4343 - d4321 * d4321;
 	if (std::abs(denom) < std::numeric_limits<F>::epsilon())
@@ -141,7 +793,7 @@ std::optional<point> line_line_intersection(const point_type<F> &p1,
 	auto pa = p1 + mua * p21;
 	auto pb = p3 + mub * p43;
 
-	return { (pa + pb) / 2.0f };
+	return { (pa + pb) / 2 };
 }
 
 /// \brief return the angle in degrees between the vectors from point @a p2 to @a p1 and @a p2 to @a p3
@@ -151,7 +803,7 @@ constexpr auto angle(const point_type<F> &p1, const point_type<F> &p2, const poi
 	point_type<F> v1 = p1 - p2;
 	point_type<F> v2 = p3 - p2;
 
-	return std::acos(glm::dot(v1, v2) / (glm::length(v1) * glm::length(v2))) * 180 / std::numbers::pi_v<F>;
+	return std::acos(dot_product(v1, v2) / (v1.length() * v2.length())) * 180 / std::numbers::pi_v<F>;
 }
 
 /// \brief return the dihedral angle in degrees for the four points @a p1, @a p2, @a p3 and @a p4
@@ -165,18 +817,18 @@ constexpr auto dihedral_angle(const point_type<F> &p1, const point_type<F> &p2,
 
 	point_type<F> z = p2 - p3; // vector from p3 to p2
 
-	point_type<F> p = glm::cross(z, v12);
-	point_type<F> x = glm::cross(z, v43);
-	point_type<F> y = glm::cross(z, x);
+	point_type<F> p = cross_product(z, v12);
+	point_type<F> x = cross_product(z, v43);
+	point_type<F> y = cross_product(z, x);
 
-	auto u = glm::dot(x, x);
-	auto v = glm::dot(y, y);
+	auto u = dot_product(x, x);
+	auto v = dot_product(y, y);
 
 	F result = 360;
 	if (u > 0 and v > 0)
 	{
-		u = glm::dot(p, x) / std::sqrt(u);
-		v = glm::dot(p, y) / std::sqrt(v);
+		u = dot_product(p, x) / std::sqrt(u);
+		v = dot_product(p, y) / std::sqrt(v);
 		if (u != 0 or v != 0)
 			result = std::atan2(v, u) * static_cast<F>(180 / std::numbers::pi_v<F>);
 	}
@@ -191,9 +843,9 @@ constexpr auto cosinus_angle(const point_type<F> &p1, const point_type<F> &p2, c
 	point_type<F> v12 = p1 - p2;
 	point_type<F> v34 = p3 - p4;
 
-	auto x = glm::dot(v12, v12) * glm::dot(v34, v34);
+	auto x = dot_product(v12, v12) * dot_product(v34, v34);
 
-	return x > 0 ? glm::dot(v12, v34) / std::sqrt(x) : 0;
+	return x > 0 ? dot_product(v12, v34) / std::sqrt(x) : 0;
 }
 
 /// \brief return the distance from point @a p to the line from @a l1 to @a l2
@@ -203,26 +855,28 @@ constexpr auto distance_point_to_line(const point_type<F> &l1, const point_type<
 	auto line = l2 - l1;
 	auto p_to_l1 = p - l1;
 	auto p_to_l2 = p - l2;
-	auto cross = glm::cross(p_to_l1, p_to_l2);
-	return glm::length(cross) / glm::length(line);
+	auto cross = cross_product(p_to_l1, p_to_l2);
+	return cross.length() / line.length();
 }
 
 /// \brief return the smallest sphere around the points in @a pts
 std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts);
 
+// --------------------------------------------------------------------
+/**
+ * @brief For e.g. simulated annealing, returns a new point that is moved in
+ * a random direction with a distance randomly chosen from a normal
+ * distribution with a stddev of offset.
+ */
+point nudge(point p, float offset);
+
 // --------------------------------------------------------------------
 
 /// \brief Return a quaternion created from angle @a angle and axis @a axis
-constexpr quaternion construct_from_angle_axis(float angle, point axis)
-{
-	return glm::angleAxis(glm::radians(angle), glm::normalize(axis));
-}
+quaternion construct_from_angle_axis(float angle, point axis);
 
 /// \brief Return a tuple of an angle and an axis for quaternion @a q
-constexpr std::tuple<float, point> quaternion_to_angle_axis(quaternion q)
-{
-	return { glm::degrees(glm::angle(q)), glm::axis(q) };
-}
+std::tuple<float, point> quaternion_to_angle_axis(quaternion q);
 
 /// @brief Given four points and an angle, return the quaternion required to rotate
 /// point p4 along the p2-p3 axis and around point p3 to obtain the required within
@@ -248,19 +902,4 @@ quaternion align_points(const std::vector<point> &a, const std::vector<point> &b
 /// \brief The RMSd for the points in \a a and \a b
 double RMSd(const std::vector<point> &a, const std::vector<point> &b);
 
-
 } // namespace cif
-
-template <>
-struct std::formatter<glm::vec3> : std::formatter<std::string_view> // NOLINT
-{
-	template <class FmtContext>
-	auto format(glm::vec3 pt, FmtContext &ctx) const
-	{
-		std::string temp;
-		std::format_to(std::back_inserter(temp), "( {}, {}, {} )", pt.x, pt.y, pt.z);
-
-		return std::formatter<std::string_view>::format(temp, ctx);
-	}
-};
-

include/cif++/symmetry.hpp

@@ -27,11 +27,11 @@
 #pragma once
 
 #include "cif++/exports.hpp"
+#include "cif++/matrix.hpp"
 #include "cif++/point.hpp"
 
 #include <array>
 #include <cstdint>
-#include <glm/ext/matrix_float3x3.hpp>
 #include <string>
 
 #if defined(__cpp_impl_three_way_comparison)
@@ -49,6 +49,18 @@ namespace cif
 
 // --------------------------------------------------------------------
 
+/// \brief Apply matrix transformation @a m on point @a pt and return the result
+inline point operator*(const matrix3x3<float> &m, const point &pt)
+{
+	return {
+		m(0, 0) * pt.m_x + m(0, 1) * pt.m_y + m(0, 2) * pt.m_z,
+		m(1, 0) * pt.m_x + m(1, 1) * pt.m_y + m(1, 2) * pt.m_z,
+		m(2, 0) * pt.m_x + m(2, 1) * pt.m_y + m(2, 2) * pt.m_z
+	};
+}
+
+// --------------------------------------------------------------------
+
 /// \brief the space groups we know
 enum class space_group_name
 {
@@ -314,7 +326,7 @@ class transformation
 	transformation(const symop_data &data);
 
 	/// \brief constructor taking a rotation matrix @a r and a translation vector @a t
-	transformation(const glm::mat3 &r, const cif::point &t);
+	transformation(const matrix3x3<float> &r, const cif::point &t);
 
 	/** @cond */
 	transformation(const transformation &) = default;
@@ -326,8 +338,8 @@ class transformation
 	/// \brief operator() to perform the transformation on point @a pt and return the result
 	point operator()(point pt) const
 	{
-		if (m_q != quaternion{})
-			pt = m_q * pt;
+		if (m_q)
+			pt.rotate(m_q);
 		else
 			pt = m_rotation * pt;
 
@@ -355,9 +367,9 @@ class transformation
 
 	void try_create_quaternion();
 
-	glm::mat3 m_rotation{};
-	quaternion m_q{};
-	point m_translation{};
+	matrix3x3<float> m_rotation;
+	quaternion m_q;
+	point m_translation;
 };
 
 // --------------------------------------------------------------------
@@ -387,14 +399,14 @@ class cell
 
 	[[nodiscard]] float get_volume() const; ///< return the calculated volume for this cell
 
-	[[nodiscard]] glm::mat3 get_orthogonal_matrix() const { return m_orthogonal; } ///< return the matrix to use to transform coordinates from fractional to orthogonal
-	[[nodiscard]] glm::mat3 get_fractional_matrix() const { return m_fractional; } ///< return the matrix to use to transform coordinates from orthogonal to fractional
+	[[nodiscard]] matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; } ///< return the matrix to use to transform coordinates from fractional to orthogonal
+	[[nodiscard]] matrix3x3<float> get_fractional_matrix() const { return m_fractional; } ///< return the matrix to use to transform coordinates from orthogonal to fractional
 
   private:
 	void init();
 
 	float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
-	glm::mat3 m_orthogonal, m_fractional;
+	matrix3x3<float> m_orthogonal, m_fractional;
 };
 
 // --------------------------------------------------------------------

include/cif++/validate.hpp

@@ -196,7 +196,7 @@ class validation_exception : public std::runtime_error
   public:
 	// Constructors
 	/// @cond
-	
+
 	validation_exception(validation_error err)
 		: validation_exception(make_error_code(err))
 	{
@@ -337,7 +337,7 @@ struct item_validator
 	std::string m_item_name;           ///< The item name
 	bool m_mandatory;                  ///< Flag indicating this item is mandatory
 	const type_validator *m_type;      ///< The type for this item
-	cif::iset m_enums;                 ///< If filled, the set of allowed values
+	std::set<std::string> m_enums;     ///< If filled, the set of allowed values
 	std::string m_default;             ///< If filled, a default value for this item
 	std::string m_category;            ///< The category this item_validator belongs to
 	std::vector<item_alias> m_aliases; ///< The aliases for this item
@@ -502,15 +502,29 @@ class validator
 	/// @brief Bottleneck function to report an error in validation
 	void report_error(std::error_code ec, bool fatal = true) const;
 
+	/// @brief Bottleneck function to report an error in validation
+	void report_error(validation_error err, std::string value, std::string_view category,
+		std::string_view item, bool fatal = true) const
+	{
+		report_error(make_error_code(err), value, category, item, fatal);
+	}
+
 	/// @brief Bottleneck function to report an error in validation
 	void report_error(validation_error err, std::string_view category,
 		std::string_view item, bool fatal = true) const
 	{
-		report_error(make_error_code(err), category, item, fatal);
+		report_error(make_error_code(err), "", category, item, fatal);
 	}
 
 	/// @brief Bottleneck function to report an error in validation
 	void report_error(std::error_code ec, std::string_view category,
+		std::string_view item, bool fatal = true) const
+	{
+		report_error(ec, "", category, item, fatal);
+	}
+
+	/// @brief Bottleneck function to report an error in validation
+	void report_error(std::error_code ec, std::string value, std::string_view category,
 		std::string_view item, bool fatal = true) const;
 
 	/// @brief Write out the audit_conform data for this validator

src/category.cpp

@@ -749,16 +749,38 @@ void category::set_validator(const validator *v, datablock &db)
 
 		bool number = type->m_primitive_type == DDL_PrimitiveType::Numb;
 		if (number)
-			continue;
-
-		for (auto row = m_head; row != nullptr; row = row->m_next)
 		{
-			if (cix >= row->size() or row->operator[](cix).empty())
-				continue;
+			for (auto row = m_head; row != nullptr; row = row->m_next)
+			{
+				if (cix >= row->size() or row->operator[](cix).empty())
+					continue;
 
-			item_value &v = row->operator[](cix);
-			if (v.is_number())
-				v = v.str();
+				item_value &v = row->operator[](cix);
+				if (not v.is_number())
+				{
+					// Try cast the value to a number and throw in case of failure
+					try
+					{
+						v.cast_to_int();
+					}
+					catch (...)
+					{
+						v.cast_to_float();
+					}
+				}
+			}
+		}
+		else
+		{
+			for (auto row = m_head; row != nullptr; row = row->m_next)
+			{
+				if (cix >= row->size() or row->operator[](cix).empty())
+					continue;
+
+				item_value &v = row->operator[](cix);
+				if (v.is_number())
+					v = v.str();
+			}
 		}
 	}
 
@@ -890,11 +912,13 @@ bool category::is_valid() const
 				seen = true;
 				std::error_code ec;
 
-				iv->validate_value(*ri->get(cix), ec);
+				auto &v = *ri->get(cix);
+
+				iv->validate_value(v, ec);
 
 				if (ec != std::errc{})
 				{
-					m_validator->report_error(ec, m_name, m_items[cix].m_name, false);
+					m_validator->report_error(ec, v.str(), m_name, m_items[cix].m_name, false);
 					continue;
 				}
 			}
@@ -1450,7 +1474,12 @@ void category::update_value(const std::vector<row_handle> &rows, std::string_vie
 
 			std::error_code ec;
 			col.m_validator->validate_value(value, ec);
-			if (ec)
+			if (ec == cif::make_error_code(cif::validation_error::value_is_not_in_enumeration_list))
+			{
+				if (cif::VERBOSE >= 0)
+					m_validator->report_error(ec, m_name, item_name, false);
+			}
+			else if (ec)
 				throw validation_exception(ec, m_name, item_name);
 		}
 	}
@@ -1577,7 +1606,18 @@ void category::update_value(row *row, uint16_t item, item_value value, bool upda
 
 	// check the value
 	if (col.m_validator and validate)
-		col.m_validator->validate_value(value);
+	{
+		std::error_code ec;
+		col.m_validator->validate_value(value, ec);
+		if (ec == cif::make_error_code(cif::validation_error::value_is_not_in_enumeration_list))
+		{
+			if (cif::VERBOSE >= 0)
+				m_validator->report_error(ec, m_name, m_items[item].m_name, false);
+		}
+		else if (ec)
+			throw validation_exception(ec, m_name, m_items[item].m_name);
+
+	}
 
 	// If the item is part of the Key for this category, remove it from the index
 	// before updating
@@ -1789,6 +1829,12 @@ category::iterator category::insert_impl(const_iterator pos, row *n)
 						}
 						else if (ec == cif::make_error_code(cif::validation_error::value_is_not_a_char_string))
 							v->cast_to_string();
+						else if (ec == cif::make_error_code(cif::validation_error::value_is_not_in_enumeration_list))
+						{
+							if (cif::VERBOSE >= 0)
+								m_validator->report_error(ec, m_name, m_items[ix].m_name, false);
+							continue;
+						}
 						else
 							throw std::system_error(ec, "Attempt to insert invalid value");
 
@@ -2158,7 +2204,7 @@ void category::write_cif(std::ostream &os, const std::vector<uint16_t> &order, b
 					offset = 0;
 				}
 
-				offset = detail::write_value(os, s, offset, w, /* right_aligned[cix] */ iv->is_number());
+				offset = detail::write_value(os, s, offset, w, right_aligned[cix]/*  iv->is_number() */);
 
 				if (offset > 132)
 				{

src/dictionary_parser.cpp

@@ -25,6 +25,8 @@
  */
 
 #include "cif++/cif++.hpp"
+#include "cif++/text.hpp"
+#include "cif++/validate.hpp"
 
 #include <cstddef>
 #include <exception>
@@ -103,7 +105,19 @@ class dictionary_parser : public parser
 				error("Undefined category '" + iv.first);
 
 			for (auto &v : iv.second)
+			{
+				// enums, make lower case if needed
+				auto tv = v.m_type;
+				if (tv and tv->m_primitive_type == DDL_PrimitiveType::UChar)
+				{
+					std::set<std::string> es;
+					for (auto &e : v.m_enums)
+						es.emplace(cif::to_lower_copy(e));
+					std::swap(es, v.m_enums);
+				}
+
 				const_cast<category_validator *>(cv)->add_item_validator(std::move(v));
+			}
 		}
 
 		// check all item validators for having a typeValidator
@@ -275,9 +289,11 @@ class dictionary_parser : public parser
 			if (typeCode.has_value())
 				tv = m_validator.get_validator_for_type(*typeCode);
 
-			iset ess;
+			std::set<std::string> ess;
 			for (auto e : dict["item_enumeration"])
 				ess.insert(e["value"].get<std::string>());
+			for (auto e : dict["pdbx_item_enumeration"])
+				ess.insert(e["value"].get<std::string>());
 
 			std::string defaultValue;
 			if (auto &cat = dict["item_default"]; not cat.empty())

src/model.cpp

@@ -67,9 +67,9 @@ void atom::atom_impl::moveTo(const point &p)
 
 	auto r = row();
 
-	r.assign("Cartn_x", { p.x, 3 }, false, false);
-	r.assign("Cartn_y", { p.y, 3 }, false, false);
-	r.assign("Cartn_z", { p.z, 3 }, false, false);
+	r.assign("Cartn_x", { p.m_x, 3 }, false, false);
+	r.assign("Cartn_y", { p.m_y, 3 }, false, false);
+	r.assign("Cartn_z", { p.m_z, 3 }, false, false);
 
 	m_location = p;
 }
@@ -689,8 +689,8 @@ float monomer::chiral_volume() const
 		auto atom2 = get_atom_by_atom_id("CD1");
 		auto atom3 = get_atom_by_atom_id("CD2");
 
-		result = dot(atom1.get_location() - centre.get_location(),
-			cross(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
+		result = dot_product(atom1.get_location() - centre.get_location(),
+			cross_product(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
 	}
 	else if (m_compound_id == "VAL")
 	{
@@ -699,8 +699,8 @@ float monomer::chiral_volume() const
 		auto atom2 = get_atom_by_atom_id("CG1");
 		auto atom3 = get_atom_by_atom_id("CG2");
 
-		result = dot(atom1.get_location() - centre.get_location(),
-			cross(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
+		result = dot_product(atom1.get_location() - centre.get_location(),
+			cross_product(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
 	}
 
 	return result;
@@ -2178,8 +2178,8 @@ std::string structure::create_non_poly(const std::string &entity_id, const std::
 			{ "label_comp_id", comp_id },
 			{ "label_asym_id", asym_id },
 			{ "label_entity_id", entity_id },
-			{ "label_seq_id", nullptr },
-			{ "pdbx_PDB_ins_code", "" },
+			{ "label_seq_id", cif::item_value_type::INAPPLICABLE },
+			{ "pdbx_PDB_ins_code", cif::item_value_type::MISSING },
 			{ "Cartn_x", atom.get_property_value("Cartn_x") },
 			{ "Cartn_y", atom.get_property_value("Cartn_y") },
 			{ "Cartn_z", atom.get_property_value("Cartn_z") },
@@ -2190,7 +2190,7 @@ std::string structure::create_non_poly(const std::string &entity_id, const std::
 			{ "auth_comp_id", comp_id },
 			{ "auth_asym_id", asym_id },
 			{ "auth_atom_id", atom.get_property_value("label_atom_id") },
-			{ "pdbx_PDB_model_num", 1 } });
+			{ "pdbx_PDB_model_num", m_model_nr } });
 
 		auto &newAtom = emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id));
 		res.add_atom(newAtom);
@@ -2208,7 +2208,7 @@ std::string structure::create_non_poly(const std::string &entity_id, const std::
 		{ "pdb_mon_id", comp_id },
 		{ "auth_mon_id", comp_id },
 		{ "pdb_strand_id", asym_id },
-		{ "pdb_ins_code", nullptr },
+		{ "pdb_ins_code", cif::item_value_type::MISSING },
 	});
 
 	return asym_id;
@@ -2244,12 +2244,14 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 
 		atom.set_value_if_empty({ "group_PDB", "HETATM" });
 		atom.set_value_if_empty({ "label_comp_id", comp_id });
-		atom.set_value_if_empty({ "label_seq_id", nullptr });
+		atom.set_value_if_empty({ "label_seq_id", cif::item_value_type::INAPPLICABLE });
 		atom.set_value_if_empty({ "auth_comp_id", comp_id });
 		atom.set_value_if_empty({ "auth_seq_id", "1" });
-		atom.set_value_if_empty({ "pdbx_PDB_model_num", 1 });
-		atom.set_value_if_empty({ "label_alt_id", "" });
+		atom.set_value_if_empty({ "pdbx_PDB_model_num", m_model_nr });
+		atom.set_value_if_empty({ "label_alt_id", cif::item_value_type::INAPPLICABLE });
 		atom.set_value_if_empty({ "occupancy", { 1.0, 2 } });
+		atom.set_value_if_empty({ "pdbx_formal_charge", cif::item_value_type::MISSING });
+		atom.set_value_if_empty({ "pdbx_tls_group_id", cif::item_value_type::MISSING });
 
 		auto row = atom_site.emplace(atom.begin(), atom.end());
 
@@ -2288,16 +2290,18 @@ std::string structure::create_non_poly(const std::string &compound_id, bool skip
 		if (skip_hydrogen and cif::atom_type_traits(a.type_symbol).symbol() == "H")
 			continue;
 
-		auto aloc = a.get_location();
+		auto ax = a.get_location().get_x();
+		auto ay = a.get_location().get_y();
+		auto az = a.get_location().get_z();
 
 		atoms.emplace_back(cif::row_initializer{
 			{ "type_symbol", cif::atom_type_traits(a.type_symbol).symbol() },
 			{ "label_atom_id", a.id },
 			{ "auth_atom_id", a.id },
-			{ "Cartn_x", aloc.x },
-			{ "Cartn_y", aloc.y },
-			{ "Cartn_z", aloc.z },
-			{ "B_iso_or_equiv", 30.00 } });
+			{ "Cartn_x", ax },
+			{ "Cartn_y", ay },
+			{ "Cartn_z", az },
+			{ "B_iso_or_equiv", { 30.00, 2 } } });
 	}
 
 	return create_non_poly(create_non_poly_entity(compound_id), atoms);

src/pdb/pdb2cif.cpp

@@ -2193,7 +2193,7 @@ void PDBFileParser::ParseRemarks()
 							if (std::regex_match(line, m, rx))
 							{
 								models[0] = std::stoi(m[1].str());
-								models[1] = stoi(m[2].str());
+								models[1] = std::stoi(m[2].str());
 							}
 							else
 								headerSeen = cif::contains(line, "RES C SSSEQI");

src/point.cpp

@@ -34,9 +34,8 @@
 #include <cstdint>
 #include <cstdlib>
 #include <initializer_list>
+#include <numbers>
 #include <optional>
-#include <glm/ext/matrix_double3x3.hpp>
-#include <glm/ext/matrix_double4x4.hpp>
 #include <random>
 #include <stdexcept>
 #include <tuple>
@@ -48,26 +47,84 @@ namespace cif
 
 // --------------------------------------------------------------------
 
+template <typename T>
+quaternion_type<T> normalize(quaternion_type<T> q)
+{
+	std::valarray<double> t(4);
+
+	t[0] = q.get_a();
+	t[1] = q.get_b();
+	t[2] = q.get_c();
+	t[3] = q.get_d();
+
+	t *= t;
+
+	double length = std::sqrt(t.sum());
+
+	if (length > 0.001)
+		q /= static_cast<quaternion::value_type>(length);
+	else
+		q = quaternion(1, 0, 0, 0);
+
+	return q;
+}
+
+// --------------------------------------------------------------------
+
+quaternion construct_from_angle_axis(float angle, point axis)
+{
+	angle = (angle * std::numbers::pi_v<float> / 180) / 2;
+	auto s = std::sin(angle);
+	auto c = std::cos(angle);
+
+	axis.normalize();
+
+	return normalize(quaternion{
+		static_cast<float>(c),
+		static_cast<float>(s * axis.m_x),
+		static_cast<float>(s * axis.m_y),
+		static_cast<float>(s * axis.m_z) });
+}
+
+std::tuple<float, point> quaternion_to_angle_axis(quaternion q)
+{
+	if (q.get_a() > 1)
+		q = normalize(q);
+
+	// angle:
+	float angle = 2 * std::acos(q.get_a());
+	angle = angle * 180 / std::numbers::pi_v<float>;
+
+	// axis:
+	float s = std::sqrt(1 - q.get_a() * q.get_a());
+	if (s < 0.001)
+		s = 1;
+
+	point axis(q.get_b() / s, q.get_c() / s, q.get_d() / s);
+
+	return { angle, axis };
+}
+
 point center_points(std::vector<point> &Points)
 {
-	point t{};
+	point t;
 
 	for (point &pt : Points)
 	{
-		t.x += pt.x;
-		t.y += pt.y;
-		t.z += pt.z;
+		t.m_x += pt.m_x;
+		t.m_y += pt.m_y;
+		t.m_z += pt.m_z;
 	}
 
-	t.x /= static_cast<float>(Points.size());
-	t.y /= static_cast<float>(Points.size());
-	t.z /= static_cast<float>(Points.size());
+	t.m_x /= static_cast<float>(Points.size());
+	t.m_y /= static_cast<float>(Points.size());
+	t.m_z /= static_cast<float>(Points.size());
 
 	for (point &pt : Points)
 	{
-		pt.x -= t.x;
-		pt.y -= t.y;
-		pt.z -= t.z;
+		pt.m_x -= t.m_x;
+		pt.m_y -= t.m_y;
+		pt.m_z -= t.m_z;
 	}
 
 	return t;
@@ -106,9 +163,9 @@ double RMSd(const std::vector<point> &a, const std::vector<point> &b)
 	{
 		std::valarray<double> d(3);
 
-		d[0] = b[i].x - a[i].x;
-		d[1] = b[i].y - a[i].y;
-		d[2] = b[i].z - a[i].z;
+		d[0] = b[i].m_x - a[i].m_x;
+		d[1] = b[i].m_y - a[i].m_y;
+		d[2] = b[i].m_z - a[i].m_z;
 
 		d *= d;
 
@@ -157,82 +214,43 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
 	return t.max();
 }
 
-/**
- * @brief Implementation of a cofactor calculation as a matrix expression
- *
- * @tparam M Type of matrix
- */
-
-auto cofactors(glm::dmat4 &m)
-{
-	glm::dmat4 result{};
-
-	const std::size_t ixs[4][3] = {
-		{ 1, 2, 3 },
-		{ 0, 2, 3 },
-		{ 0, 1, 3 },
-		{ 0, 1, 2 }
-	};
-
-	for (size_t i = 0; i < 4; ++i)
-	{
-		for (size_t j = 0; j < 4; ++j)
-		{
-			const std::size_t *ix = ixs[i];
-			const std::size_t *iy = ixs[j];
-
-			auto e =
-				m[ix[0]][iy[0]] * m[ix[1]][iy[1]] * m[ix[2]][iy[2]] +
-				m[ix[0]][iy[1]] * m[ix[1]][iy[2]] * m[ix[2]][iy[0]] +
-				m[ix[0]][iy[2]] * m[ix[1]][iy[0]] * m[ix[2]][iy[1]] -
-				m[ix[0]][iy[2]] * m[ix[1]][iy[1]] * m[ix[2]][iy[0]] -
-				m[ix[0]][iy[1]] * m[ix[1]][iy[0]] * m[ix[2]][iy[2]] -
-				m[ix[0]][iy[0]] * m[ix[1]][iy[2]] * m[ix[2]][iy[1]];
-
-			result[i][j] = (i + j) % 2 == 1 ? -e : e;
-		}
-	}
-
-	return result;
-}
-
 quaternion align_points(const std::vector<point> &pa, const std::vector<point> &pb)
 {
 	// First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
-	glm::dmat3x3 M{};
+	matrix3x3<double> M;
 
 	for (uint32_t i = 0; i < pa.size(); ++i)
 	{
 		const point &a = pa[i];
 		const point &b = pb[i];
 
-		M[0][0] += a.x * b.x;
-		M[0][1] += a.x * b.y;
-		M[0][2] += a.x * b.z;
-		M[1][0] += a.y * b.x;
-		M[1][1] += a.y * b.y;
-		M[1][2] += a.y * b.z;
-		M[2][0] += a.z * b.x;
-		M[2][1] += a.z * b.y;
-		M[2][2] += a.z * b.z;
+		M(0, 0) += a.m_x * b.m_x;
+		M(0, 1) += a.m_x * b.m_y;
+		M(0, 2) += a.m_x * b.m_z;
+		M(1, 0) += a.m_y * b.m_x;
+		M(1, 1) += a.m_y * b.m_y;
+		M(1, 2) += a.m_y * b.m_z;
+		M(2, 0) += a.m_z * b.m_x;
+		M(2, 1) += a.m_z * b.m_y;
+		M(2, 2) += a.m_z * b.m_z;
 	}
 
 	// Now calculate N, a symmetric 4x4 matrix
-	glm::dmat4x4 N;
+	symmetric_matrix4x4<double> N(4);
 
-	N[0][0] = M[0][0] + M[1][1] + M[2][2];
-	N[0][1] = N[1][0] = M[1][2] - M[2][1];
-	N[0][2] = N[2][0] = M[2][0] - M[0][2];
-	N[0][3] = N[3][0] = M[0][1] - M[1][0];
+	N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
+	N(0, 1) = M(1, 2) - M(2, 1);
+	N(0, 2) = M(2, 0) - M(0, 2);
+	N(0, 3) = M(0, 1) - M(1, 0);
 
-	N[1][1] = M[0][0] - M[1][1] - M[2][2];
-	N[1][2] = N[2][1] = M[0][1] + M[1][0];
-	N[1][3] = N[3][1] = M[0][2] + M[2][0];
+	N(1, 1) = M(0, 0) - M(1, 1) - M(2, 2);
+	N(1, 2) = M(0, 1) + M(1, 0);
+	N(1, 3) = M(0, 2) + M(2, 0);
 
-	N[2][2] = -M[0][0] + M[1][1] - M[2][2];
-	N[2][3] = N[3][2] = M[1][2] + M[2][1];
+	N(2, 2) = -M(0, 0) + M(1, 1) - M(2, 2);
+	N(2, 3) = M(1, 2) + M(2, 1);
 
-	N[3][3] = -M[0][0] - M[1][1] + M[2][2];
+	N(3, 3) = -M(0, 0) - M(1, 1) + M(2, 2);
 
 	// det(N - λI) = 0
 	// find the largest λ (λm)
@@ -243,41 +261,41 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	// and so this is a so-called depressed quartic
 	// solve it using Ferrari's algorithm
 
-	double C = -2 * (M[0][0] * M[0][0] + M[0][1] * M[0][1] + M[0][2] * M[0][2] +
-						M[1][0] * M[1][0] + M[1][1] * M[1][1] + M[1][2] * M[1][2] +
-						M[2][0] * M[2][0] + M[2][1] * M[2][1] + M[2][2] * M[2][2]);
+	double C = -2 * (M(0, 0) * M(0, 0) + M(0, 1) * M(0, 1) + M(0, 2) * M(0, 2) +
+						M(1, 0) * M(1, 0) + M(1, 1) * M(1, 1) + M(1, 2) * M(1, 2) +
+						M(2, 0) * M(2, 0) + M(2, 1) * M(2, 1) + M(2, 2) * M(2, 2));
 
-	double D = 8 * (M[0][0] * M[1][2] * M[2][1] +
-					   M[1][1] * M[2][0] * M[0][2] +
-					   M[2][2] * M[0][1] * M[1][0]) -
-	           8 * (M[0][0] * M[1][1] * M[2][2] +
-					   M[1][2] * M[2][0] * M[0][1] +
-					   M[2][1] * M[1][0] * M[0][2]);
+	double D = 8 * (M(0, 0) * M(1, 2) * M(2, 1) +
+					   M(1, 1) * M(2, 0) * M(0, 2) +
+					   M(2, 2) * M(0, 1) * M(1, 0)) -
+	           8 * (M(0, 0) * M(1, 1) * M(2, 2) +
+					   M(1, 2) * M(2, 0) * M(0, 1) +
+					   M(2, 1) * M(1, 0) * M(0, 2));
 
 	// E is the determinant of N:
 	double E =
-		(N[0][0] * N[1][1] - N[0][1] * N[0][1]) * (N[2][2] * N[3][3] - N[2][3] * N[2][3]) +
-		(N[0][1] * N[0][2] - N[0][0] * N[2][1]) * (N[2][1] * N[3][3] - N[2][3] * N[1][3]) +
-		(N[0][0] * N[1][3] - N[0][1] * N[0][3]) * (N[2][1] * N[2][3] - N[2][2] * N[1][3]) +
-		(N[0][1] * N[2][1] - N[1][1] * N[0][2]) * (N[0][2] * N[3][3] - N[2][3] * N[0][3]) +
-		(N[1][1] * N[0][3] - N[0][1] * N[1][3]) * (N[0][2] * N[2][3] - N[2][2] * N[0][3]) +
-		(N[0][2] * N[1][3] - N[2][1] * N[0][3]) * (N[0][2] * N[1][3] - N[2][1] * N[0][3]);
+		(N(0, 0) * N(1, 1) - N(0, 1) * N(0, 1)) * (N(2, 2) * N(3, 3) - N(2, 3) * N(2, 3)) +
+		(N(0, 1) * N(0, 2) - N(0, 0) * N(2, 1)) * (N(2, 1) * N(3, 3) - N(2, 3) * N(1, 3)) +
+		(N(0, 0) * N(1, 3) - N(0, 1) * N(0, 3)) * (N(2, 1) * N(2, 3) - N(2, 2) * N(1, 3)) +
+		(N(0, 1) * N(2, 1) - N(1, 1) * N(0, 2)) * (N(0, 2) * N(3, 3) - N(2, 3) * N(0, 3)) +
+		(N(1, 1) * N(0, 3) - N(0, 1) * N(1, 3)) * (N(0, 2) * N(2, 3) - N(2, 2) * N(0, 3)) +
+		(N(0, 2) * N(1, 3) - N(2, 1) * N(0, 3)) * (N(0, 2) * N(1, 3) - N(2, 1) * N(0, 3));
 
 	// solve quartic
 	double lambda = LargestDepressedQuarticSolution(C, D, E);
 
 	// calculate t = (N - λI)
-	auto t = N - glm::dmat4x4(1.0) * lambda;
+	matrix<double> t(N - identity_matrix(4) * lambda);
 
 	// calculate a matrix of cofactors for t
-	auto cf = cofactors(t);
+	auto cf = matrix_cofactors(t);
 
 	int maxR = 0;
-	double maxCF = std::abs(cf[0][0]);
+	double maxCF = std::abs(cf(0, 0));
 
 	for (int r = 1; r < 4; ++r)
 	{
-		auto cfr = std::abs(cf[r][0]);
+		auto cfr = std::abs(cf(r, 0));
 		if (maxCF < cfr)
 		{
 			maxCF = cfr;
@@ -286,10 +304,10 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	}
 
 	quaternion q(
-		static_cast<float>(cf[maxR][0]),
-		static_cast<float>(cf[maxR][1]),
-		static_cast<float>(cf[maxR][2]),
-		static_cast<float>(cf[maxR][3]));
+		static_cast<float>(cf(maxR, 0)),
+		static_cast<float>(cf(maxR, 1)),
+		static_cast<float>(cf(maxR, 2)),
+		static_cast<float>(cf(maxR, 3)));
 	q = normalize(q);
 
 	return q;
@@ -297,29 +315,52 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 
 // --------------------------------------------------------------------
 
-std::tuple<point, float> smallest_sphere_around_2_points(std::array<point, 2> pts)
+point nudge(point p, float offset)
 {
-	return { (pts[0] + pts[1]) / 2.0f, distance(pts[0], pts[1]) / 2.0f };
+	static std::random_device rd;
+	static std::mt19937_64 rng(rd());
+
+	std::uniform_real_distribution<float> randomAngle(0, 2 * std::numbers::pi);
+	std::normal_distribution<float> randomOffset(0, offset);
+
+	float theta = randomAngle(rng);
+	float phi1 = randomAngle(rng) - static_cast<float>(std::numbers::pi);
+	float phi2 = randomAngle(rng) - static_cast<float>(std::numbers::pi);
+
+	quaternion q = spherical(1.0f, theta, phi1, phi2);
+
+	point r{ 0, 0, 1 };
+	r.rotate(q);
+	r *= randomOffset(rng);
+
+	return p + r;
 }
 
-std::tuple<point, float> smallest_sphere_around_3_points(std::array<point, 3> pts)
+// --------------------------------------------------------------------
+
+std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
+{
+	return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
+}
+
+std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
 {
 	// Find two bisectors
-	auto vz = cross(pts[1] - pts[0], pts[2] - pts[0]);
+	auto vz = cross_product(pts[1] - pts[0], pts[2] - pts[0]);
 
-	auto bs1 = cross(vz, pts[1] - pts[0]);
-	bs1 = glm::normalize(bs1);
+	auto bs1 = cross_product(vz, pts[1] - pts[0]);
+	bs1.normalize();
 
 	auto v1 = (pts[1] - pts[0]);
-	v1 = glm::normalize(v1);
+	v1.normalize();
 
 	auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
 
-	auto bs2 = cross(vz, pts[2] - pts[0]);
-	bs2 = glm::normalize(bs2);
+	auto bs2 = cross_product(vz, pts[2] - pts[0]);
+	bs2.normalize();
 
 	auto v2 = (pts[2] - pts[0]);
-	v2 = glm::normalize(v2);
+	v2.normalize();
 
 	auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
 
@@ -340,7 +381,7 @@ std::tuple<point, float> smallest_sphere_around_3_points(std::array<point, 3> pt
 		return smallest_sphere_around_2_points({ pts[1], pts[2] });
 }
 
-std::tuple<point, float> smallest_sphere_around_4_points(std::array<point, 4> pts)
+std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
 {
 	auto t0 = -norm_squared(pts[0]);
 	auto t1 = -norm_squared(pts[1]);
@@ -348,46 +389,46 @@ std::tuple<point, float> smallest_sphere_around_4_points(std::array<point, 4> pt
 	auto t3 = -norm_squared(pts[3]);
 
 	// clang-format off
-	glm::mat4x4 Tm({
-		pts[0].x, pts[0].y, pts[0].z, 1,
-		pts[1].x, pts[1].y, pts[1].z, 1,
-		pts[2].x, pts[2].y, pts[2].z, 1,
-		pts[3].x, pts[3].y, pts[3].z, 1
+	matrix4x4<float> Tm({
+		pts[0].m_x, pts[0].m_y, pts[0].m_z, 1,
+		pts[1].m_x, pts[1].m_y, pts[1].m_z, 1,
+		pts[2].m_x, pts[2].m_y, pts[2].m_z, 1,
+		pts[3].m_x, pts[3].m_y, pts[3].m_z, 1
 	});
 	auto T = determinant(Tm);
 
 	if (T != 0)
 	{
-		glm::mat4x4 Dm({
-			t0, pts[0].y, pts[0].z, 1,
-			t1, pts[1].y, pts[1].z, 1,
-			t2, pts[2].y, pts[2].z, 1,
-			t3, pts[3].y, pts[3].z, 1
+		matrix4x4<float> Dm({
+			t0, pts[0].m_y, pts[0].m_z, 1,
+			t1, pts[1].m_y, pts[1].m_z, 1,
+			t2, pts[2].m_y, pts[2].m_z, 1,
+			t3, pts[3].m_y, pts[3].m_z, 1
 		});
 		auto D = determinant(Dm) / T;
 		
-		glm::mat4x4 Em({
-			pts[0].x, t0, pts[0].z, 1,
-			pts[1].x, t1, pts[1].z, 1,
-			pts[2].x, t2, pts[2].z, 1,
-			pts[3].x, t3, pts[3].z, 1
+		matrix4x4<float> Em({
+			pts[0].m_x, t0, pts[0].m_z, 1,
+			pts[1].m_x, t1, pts[1].m_z, 1,
+			pts[2].m_x, t2, pts[2].m_z, 1,
+			pts[3].m_x, t3, pts[3].m_z, 1
 		});
 		auto E = determinant(Em) / T;
 		
-		glm::mat4x4 Fm({
-			pts[0].x, pts[0].y, t0, 1,
-			pts[1].x, pts[1].y, t1, 1,
-			pts[2].x, pts[2].y, t2, 1,
-			pts[3].x, pts[3].y, t3, 1
+		matrix4x4<float> Fm({
+			pts[0].m_x, pts[0].m_y, t0, 1,
+			pts[1].m_x, pts[1].m_y, t1, 1,
+			pts[2].m_x, pts[2].m_y, t2, 1,
+			pts[3].m_x, pts[3].m_y, t3, 1
 		});
 		
 		auto F = determinant(Fm) / T;
 		
-		glm::mat4x4 Gm({
-			pts[0].x, pts[0].y, pts[0].z, t0,
-			pts[1].x, pts[1].y, pts[1].z, t1,
-			pts[2].x, pts[2].y, pts[2].z, t2,
-			pts[3].x, pts[3].y, pts[3].z, t3
+		matrix4x4<float> Gm({
+			pts[0].m_x, pts[0].m_y, pts[0].m_z, t0,
+			pts[1].m_x, pts[1].m_y, pts[1].m_z, t1,
+			pts[2].m_x, pts[2].m_y, pts[2].m_z, t2,
+			pts[3].m_x, pts[3].m_y, pts[3].m_z, t3
 		});
 		auto G = determinant(Gm) / T;
 		
@@ -467,19 +508,19 @@ bool point_in_circle(point p, std::vector<point> c)
 		case 2:
 		{
 			auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
-			return distance_squared(p, center) <= radius * radius;
+			return cif::distance_squared(p, center) <= radius * radius;
 		}
 
 		case 3:
 		{
 			auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
-			return distance_squared(p, center) <= radius * radius;
+			return cif::distance_squared(p, center) <= radius * radius;
 		}
 
 		case 4:
 		{
 			auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
-			return distance_squared(p, center) <= radius * radius;
+			return cif::distance_squared(p, center) <= radius * radius;
 		}
 
 		default:

src/symmetry.cpp

@@ -24,6 +24,7 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/symmetry.hpp"
 #include "cif++/cif++.hpp"
 #include "symop_table_data.hpp"
 
@@ -81,15 +82,14 @@ void cell::init()
 
 	auto alpha_star = std::acos((std::cos(gamma) * std::cos(beta) - std::cos(alpha)) / (std::sin(beta) * std::sin(gamma)));
 
-	m_orthogonal = glm::mat3(1.0f);
+	m_orthogonal = identity_matrix(3);
 
-	// WARNING: glm matrices are column major, by default
-	m_orthogonal[0][0] = m_a;
-	m_orthogonal[1][0] = m_b * std::cos(gamma);
-	m_orthogonal[2][0] = m_c * std::cos(beta);
-	m_orthogonal[1][1] = m_b * std::sin(gamma);
-	m_orthogonal[2][1] = m_c * std::sin(beta) * std::cos(alpha_star);
-	m_orthogonal[2][2] = m_c * std::sin(beta) * std::sin(alpha_star);
+	m_orthogonal(0, 0) = m_a;
+	m_orthogonal(0, 1) = m_b * std::cos(gamma);
+	m_orthogonal(0, 2) = m_c * std::cos(beta);
+	m_orthogonal(1, 1) = m_b * std::sin(gamma);
+	m_orthogonal(1, 2) = m_c * std::sin(beta) * std::cos(alpha_star);
+	m_orthogonal(2, 2) = m_c * std::sin(beta) * std::sin(alpha_star);
 
 	m_fractional = inverse(m_orthogonal);
 }
@@ -151,25 +151,24 @@ transformation::transformation(const symop_data &data)
 {
 	const auto &d = data.data();
 
-	// WARNING: glm::mat is column major
-	m_rotation[0][0] = static_cast<float>(d[0]);
-	m_rotation[1][0] = static_cast<float>(d[1]);
-	m_rotation[2][0] = static_cast<float>(d[2]);
-	m_rotation[0][1] = static_cast<float>(d[3]);
-	m_rotation[1][1] = static_cast<float>(d[4]);
-	m_rotation[2][1] = static_cast<float>(d[5]);
-	m_rotation[0][2] = static_cast<float>(d[6]);
-	m_rotation[1][2] = static_cast<float>(d[7]);
-	m_rotation[2][2] = static_cast<float>(d[8]);
+	m_rotation(0, 0) = static_cast<float>(d[0]);
+	m_rotation(0, 1) = static_cast<float>(d[1]);
+	m_rotation(0, 2) = static_cast<float>(d[2]);
+	m_rotation(1, 0) = static_cast<float>(d[3]);
+	m_rotation(1, 1) = static_cast<float>(d[4]);
+	m_rotation(1, 2) = static_cast<float>(d[5]);
+	m_rotation(2, 0) = static_cast<float>(d[6]);
+	m_rotation(2, 1) = static_cast<float>(d[7]);
+	m_rotation(2, 2) = static_cast<float>(d[8]);
 
 	try_create_quaternion();
 
-	m_translation.x = static_cast<float>(d[9] == 0 ? 0 : 1.0 * d[9] / d[10]);
-	m_translation.y = static_cast<float>(d[11] == 0 ? 0 : 1.0 * d[11] / d[12]);
-	m_translation.z = static_cast<float>(d[13] == 0 ? 0 : 1.0 * d[13] / d[14]);
+	m_translation.m_x = static_cast<float>(d[9] == 0 ? 0 : 1.0 * d[9] / d[10]);
+	m_translation.m_y = static_cast<float>(d[11] == 0 ? 0 : 1.0 * d[11] / d[12]);
+	m_translation.m_z = static_cast<float>(d[13] == 0 ? 0 : 1.0 * d[13] / d[14]);
 }
 
-transformation::transformation(const glm::mat3 &r, const cif::point &t)
+transformation::transformation(const matrix3x3<float> &r, const cif::point &t)
 	: m_rotation(r)
 	, m_translation(t)
 {
@@ -180,9 +179,9 @@ void transformation::try_create_quaternion()
 {
 	Eigen::Matrix3f rot;
 
-	rot << m_rotation[0][0], m_rotation[0][1], m_rotation[0][2],
-		m_rotation[1][0], m_rotation[1][1], m_rotation[1][2],
-		m_rotation[2][0], m_rotation[2][1], m_rotation[2][2];
+	rot << m_rotation(0, 0), m_rotation(0, 1), m_rotation(0, 2),
+		m_rotation(1, 0), m_rotation(1, 1), m_rotation(1, 2),
+		m_rotation(2, 0), m_rotation(2, 1), m_rotation(2, 2);
 
 	if (rot * rot.transpose() == Eigen::Matrix3f::Identity() and rot.determinant() == 1)
 	{
@@ -243,20 +242,20 @@ point offsetToOrigin(const cell &c, const point &p)
 {
 	point d{};
 
-	while (p.x + d.x < -(c.get_a()))
-		d.x += c.get_a();
-	while (p.x + d.x > (c.get_a()))
-		d.x -= c.get_a();
+	while (p.m_x + d.m_x < -(c.get_a()))
+		d.m_x += c.get_a();
+	while (p.m_x + d.m_x > (c.get_a()))
+		d.m_x -= c.get_a();
 
-	while (p.y + d.y < -(c.get_b()))
-		d.y += c.get_b();
-	while (p.y + d.y > (c.get_b()))
-		d.y -= c.get_b();
+	while (p.m_y + d.m_y < -(c.get_b()))
+		d.m_y += c.get_b();
+	while (p.m_y + d.m_y > (c.get_b()))
+		d.m_y -= c.get_b();
 
-	while (p.z + d.z < -(c.get_c()))
-		d.z += c.get_c();
-	while (p.z + d.z > (c.get_c()))
-		d.z -= c.get_c();
+	while (p.m_z + d.m_z < -(c.get_c()))
+		d.m_z += c.get_c();
+	while (p.m_z + d.m_z > (c.get_c()))
+		d.m_z -= c.get_c();
 
 	return d;
 };
@@ -265,20 +264,20 @@ point offsetToOriginFractional(const point &p)
 {
 	point d{};
 
-	while (p.x + d.x < -0.5f)
-		d.x += 1;
-	while (p.x + d.x > 0.5f)
-		d.x -= 1;
+	while (p.m_x + d.m_x < -0.5f)
+		d.m_x += 1;
+	while (p.m_x + d.m_x > 0.5f)
+		d.m_x -= 1;
 
-	while (p.y + d.y < -0.5f)
-		d.y += 1;
-	while (p.y + d.y > 0.5f)
-		d.y -= 1;
+	while (p.m_y + d.m_y < -0.5f)
+		d.m_y += 1;
+	while (p.m_y + d.m_y > 0.5f)
+		d.m_y -= 1;
 
-	while (p.z + d.z < -0.5f)
-		d.z += 1;
-	while (p.z + d.z > 0.5f)
-		d.z -= 1;
+	while (p.m_z + d.m_z < -0.5f)
+		d.m_z += 1;
+	while (p.m_z + d.m_z > 0.5f)
+		d.m_z -= 1;
 
 	return d;
 };
@@ -290,9 +289,9 @@ point spacegroup::operator()(const point &pt, const cell &c, sym_op symop) const
 
 	transformation t = at(symop.m_nr - 1);
 
-	t.m_translation.x += symop.m_ta - 5;
-	t.m_translation.y += symop.m_tb - 5;
-	t.m_translation.z += symop.m_tc - 5;
+	t.m_translation.m_x += symop.m_ta - 5;
+	t.m_translation.m_y += symop.m_tb - 5;
+	t.m_translation.m_z += symop.m_tc - 5;
 
 	return orthogonal(t(fractional(pt, c)), c);
 }
@@ -304,9 +303,9 @@ point spacegroup::inverse(const point &pt, const cell &c, sym_op symop) const
 
 	transformation t = at(symop.m_nr - 1);
 
-	t.m_translation.x += symop.m_ta - 5;
-	t.m_translation.y += symop.m_tb - 5;
-	t.m_translation.z += symop.m_tc - 5;
+	t.m_translation.m_x += symop.m_ta - 5;
+	t.m_translation.m_y += symop.m_tb - 5;
+	t.m_translation.m_z += symop.m_tc - 5;
 
 	auto it = cif::inverse(t);
 	return orthogonal(it(fractional(pt, c)), c);
@@ -440,7 +439,7 @@ std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b
 	if (m_cell.get_a() == 0 or m_cell.get_b() == 0 or m_cell.get_c() == 0)
 		throw std::runtime_error("Invalid cell, contains a dimension that is zero");
 
-	point result_fsb{};
+	point result_fsb;
 	float result_d = std::numeric_limits<float>::max();
 	sym_op result_s;
 
@@ -454,46 +453,46 @@ std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b
 
 	a = orthogonal(fa, m_cell);
 
-	for (std::size_t i = 0; i < m_spacegroup.size(); ++i)
+	for (uint8_t i = 0; std::cmp_less(i, m_spacegroup.size()); ++i)
 	{
-		sym_op s(static_cast<uint8_t>(i + 1));
+		sym_op s(i + 1);
 		auto &t = m_spacegroup[i];
 
 		auto fsb = t(fb);
 
-		while (fsb.x - 0.5f > fa.x)
+		while (fsb.m_x - 0.5f > fa.m_x)
 		{
-			fsb.x -= 1;
+			fsb.m_x -= 1;
 			s.m_ta -= 1;
 		}
 
-		while (fsb.x + 0.5f < fa.x)
+		while (fsb.m_x + 0.5f < fa.m_x)
 		{
-			fsb.x += 1;
+			fsb.m_x += 1;
 			s.m_ta += 1;
 		}
 
-		while (fsb.y - 0.5f > fa.y)
+		while (fsb.m_y - 0.5f > fa.m_y)
 		{
-			fsb.y -= 1;
+			fsb.m_y -= 1;
 			s.m_tb -= 1;
 		}
 
-		while (fsb.y + 0.5f < fa.y)
+		while (fsb.m_y + 0.5f < fa.m_y)
 		{
-			fsb.y += 1;
+			fsb.m_y += 1;
 			s.m_tb += 1;
 		}
 
-		while (fsb.z - 0.5f > fa.z)
+		while (fsb.m_z - 0.5f > fa.m_z)
 		{
-			fsb.z -= 1;
+			fsb.m_z -= 1;
 			s.m_tc -= 1;
 		}
 
-		while (fsb.z + 0.5f < fa.z)
+		while (fsb.m_z + 0.5f < fa.m_z)
 		{
-			fsb.z += 1;
+			fsb.m_z += 1;
 			s.m_tc += 1;
 		}
 

src/validate.cpp

@@ -27,6 +27,7 @@
 #include "cif++/validate.hpp"
 
 #include "cif++/cif++.hpp"
+#include "cif++/text.hpp"
 
 #include <cassert>
 #include <charconv>
@@ -245,8 +246,17 @@ bool item_validator::validate_value(const item_value &value, std::error_code &ec
 						ec = make_error_code(validation_error::value_does_not_match_rx);
 					}
 
-					if (ec == std::errc{} and not m_enums.empty() and m_enums.count(value.str()) == 0)
-						ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
+					if (ec == std::errc{} and not m_enums.empty())
+					{
+						bool valid =
+							m_type->m_primitive_type == DDL_PrimitiveType::UChar ? //
+								m_enums.contains(cif::to_lower_copy(value.sv()))
+																				 : //
+								m_enums.contains(std::string{ value.sv() });
+
+						if (not valid)
+							ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
+					}
 				}
 			}
 		}
@@ -436,8 +446,8 @@ void validator::report_error(std::error_code ec, bool fatal) const
 		std::cerr << ec.message() << '\n';
 }
 
-void validator::report_error(std::error_code ec, std::string_view category,
-	std::string_view item, bool fatal) const
+void validator::report_error(std::error_code ec, std::string value,
+	std::string_view category, std::string_view item, bool fatal) const
 {
 	if (m_strict or fatal)
 	{
@@ -448,10 +458,19 @@ void validator::report_error(std::error_code ec, std::string_view category,
 	}
 
 	if (VERBOSE > 0)
-		std::cerr << ec.message()
-				  << "; category: " << std::quoted(category)
-				  << " item: " << std::quoted(item)
-				  << '\n';
+	{
+		if (value.empty())
+			std::cerr << ec.message()
+					  << "; category: " << std::quoted(category)
+					  << " item: " << std::quoted(item)
+					  << '\n';
+		else
+			std::cerr << ec.message()
+					  << "; value: " << std::quoted(value)
+					  << "; category: " << std::quoted(category)
+					  << " item: " << std::quoted(item)
+					  << '\n';
+	}
 }
 
 void validator::fill_audit_conform(category &audit_conform) const

test/CMakeLists.txt

@@ -50,6 +50,7 @@ list(
 	reconstruction
 	validate-pdbx
 	cql
+	matrix
 )
 
 add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp")

test/matrix-test.cpp

@@ -0,0 +1,115 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2025 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "cif++/matrix.hpp"
+#include "test-main.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include <cif++/cif++.hpp>
+
+TEST_CASE("m1")
+{
+	cif::matrix3x3<int> m = cif::identity_matrix<int>(3);
+
+	CHECK(cif::determinant(m) == 1);
+}
+
+TEST_CASE("m2")
+{
+	cif::matrix4x4<int> m = cif::identity_matrix<int>(4);
+
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+	CHECK(ms == cif::identity_matrix<int>(3));
+}
+
+TEST_CASE("m3")
+{
+	cif::matrix4x4<int> m{
+		{ 1, 2, 3, 4,      //
+			5, 6, 7, 8,    //
+			9, 10, 11, 12, //
+			13, 14, 15, 16 }
+	};
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+
+	cif::matrix3x3<int> t{
+		{ 1, 3, 4, 9, 11, 12, 13, 15, 16 }
+	};
+
+	CHECK(ms == t);
+}
+
+TEST_CASE("m4")
+{
+	cif::matrix4x4<int> m{
+		{
+			-2,
+			3,
+			1,
+			0,
+			4,
+			1,
+			-3,
+			2,
+			0,
+			-1,
+			2,
+			5,
+			3,
+			2,
+			0,
+			-4,
+		}
+	};
+
+	std::cout << m << "\n\n";
+
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3)) << "\n\n";
+
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3))) << "\n\n";
+
+	CHECK(cif::determinant(m) == 332);
+}
+
+// --------------------------------------------------------------------
+
+TEST_CASE("m5")
+{
+    cif::matrix4x4<float> m = cif::identity_matrix<float>(4);
+    cif::matrix_fixed<float, 1, 4> v({ 0, 0.5f, 0, 1.0f });
+
+    auto mv = v * m;
+
+    CHECK(mv == v);
+
+}
+

test/model-test.cpp

@@ -443,8 +443,6 @@ TEST_CASE("test_alternates_1")
 	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
 	cif::file file(example.string());
 
-	// auto &db = file.front();
-
 	cif::mm::structure s(file);
 
 	for (auto atom : s.atoms())

test/unit-3d-test.cpp

@@ -25,6 +25,7 @@
  */
 
 #include "cif++/cif++.hpp"
+#include "cif++/symmetry.hpp"
 #include "test-main.hpp"
 
 #include <catch2/catch_test_macros.hpp>
@@ -35,7 +36,7 @@
 # pragma warning(disable : 4127) // conditional expression is constant
 #endif
 
-// #include <Eigen/Eigenvalues>
+#include <Eigen/Eigen>
 
 // --------------------------------------------------------------------
 
@@ -87,7 +88,7 @@ TEST_CASE("t1")
 	cif::center_points(p1);
 
 	for (auto &p : p2)
-		p = cif::rotate(p, q);
+		p.rotate(q);
 
 	cif::center_points(p2);
 
@@ -98,7 +99,7 @@ TEST_CASE("t1")
 	CHECK_THAT(std::fmod(360 + angle, 360), Catch::Matchers::WithinRel(std::fmod(360 - angle0, 360), 0.01));
 
 	for (auto &p : p1)
-		p = q2 * p;
+		p.rotate(q2);
 
 	auto rmsd = cif::RMSd(p1, p2);
 
@@ -115,7 +116,7 @@ TEST_CASE("t2")
 		{ 1, 2, 0 }
 	};
 
-	auto xp = glm::cross(p[1] - p[0], p[2] - p[0]);
+	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
 
 	auto q = cif::construct_from_angle_axis(45, xp);
 
@@ -132,16 +133,16 @@ TEST_CASE("t3")
 		{ 1, 2, 0 }
 	};
 
-	cif::point xp = glm::cross(p[1] - p[0], p[2] - p[0]);
+	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
 
 	auto q = cif::construct_from_angle_axis(45, xp);
 
 	auto v = p[1];
 	v -= p[0];
-	v = q * v;
+	v.rotate(q);
 	v += p[0];
 
-	std::println(std::cout, "{}", v);
+	std::cout << v << '\n';
 
 	double a = cif::angle(v, p[0], p[1]);
 
@@ -165,22 +166,22 @@ TEST_CASE("dh_q_0")
 
 	auto q = cif::construct_from_angle_axis(90, axis);
 
-	p = q * p;
+	p.rotate(q);
 
-	REQUIRE(std::abs(p.x - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.y - 0.f) < 0.01f);
-	REQUIRE(std::abs(p.z - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.m_x - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.m_y - 0.f) < 0.01f);
+	REQUIRE(std::abs(p.m_z - 1.f) < 0.01f);
 
 	a = cif::dihedral_angle(t[0], t[1], t[2], p);
 	REQUIRE(std::abs(a - 90.f) < 0.01f);
 
 	q = cif::construct_from_angle_axis(-90, axis);
 
-	p = q * p;
+	p.rotate(q);
 
-	REQUIRE(std::abs(p.x - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.y - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.z - 0.f) < 0.01f);
+	REQUIRE(std::abs(p.m_x - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.m_y - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.m_z - 0.f) < 0.01f);
 
 	a = cif::dihedral_angle(t[0], t[1], t[2], p);
 	REQUIRE(std::abs(a - 0.f) < 0.01f);
@@ -222,109 +223,38 @@ TEST_CASE("dh_q_1")
 	{
 		auto q = cif::construct_for_dihedral_angle(pts[0], pts[1], pts[2], pts[3], angle, 1);
 
-		pts[3] = cif::rotate(pts[3], q, pts[2]);
+		pts[3].rotate(q, pts[2]);
 
 		auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
 		CHECK_THAT(dh, Catch::Matchers::WithinRel(angle, 0.1f));
 	}
 }
 
-/* TEST_CASE("m2q_0a")
-{
-    for (std::size_t i = 0; i < cif::kSymopNrTableSize; ++i)
-    {
-        auto d = cif::kSymopNrTable[i].symop().data();
-
-        Eigen::Matrix3f rot;
-        rot << static_cast<float>(d[0]), static_cast<float>(d[1]), static_cast<float>(d[2]), static_cast<float>(d[3]), static_cast<float>(d[4]), static_cast<float>(d[5]), static_cast<float>(d[6]), static_cast<float>(d[7]), static_cast<float>(d[8]);
-
-        // check to see if this matrix contains a true rotation
-        if (rot * rot.transpose() != Eigen::Matrix3f::Identity() or rot.determinant() != 1)
-            continue;
-
-        Eigen::Quaternionf qe(rot);
-
-        auto q = normalize(cif::quaternion{ qe.w(), qe.x(), qe.y(), qe.z() });
-
-        cif::point p1{ 1, 1, 1 };
-        cif::point p2 = p1;
-        p2 = q * p2;
-
-        cif::glm::mat3 rot_c({ static_cast<float>(d[0]),
-            static_cast<float>(d[1]),
-            static_cast<float>(d[2]),
-            static_cast<float>(d[3]),
-            static_cast<float>(d[4]),
-            static_cast<float>(d[5]),
-            static_cast<float>(d[6]),
-            static_cast<float>(d[7]),
-            static_cast<float>(d[8]) });
-
-        cif::point p3 = rot_c * p1;
-
-        CHECK_THAT(p2.x, Catch::Matchers::WithinRel(p3.x, 0.01f));
-        CHECK_THAT(p2.y, Catch::Matchers::WithinRel(p3.y, 0.01f));
-        CHECK_THAT(p2.z, Catch::Matchers::WithinRel(p3.z, 0.01f));
-    }
-}
- */
-
-// "TEST_CASE(m2q_1")
+// TEST_CASE("m2q_1")
 // {
 // 	for (std::size_t i = 0; i < cif::kSymopNrTableSize; ++i)
 // 	{
 // 		auto d = cif::kSymopNrTable[i].symop().data();
 
-// 		cif::glm::mat3 rot;
-// 		float Qxx = rot(0, 0) = d[0];
-// 		float Qxy = rot(0, 1) = d[1];
-// 		float Qxz = rot(0, 2) = d[2];
-// 		float Qyx = rot(1, 0) = d[3];
-// 		float Qyy = rot(1, 1) = d[4];
-// 		float Qyz = rot(1, 2) = d[5];
-// 		float Qzx = rot(2, 0) = d[6];
-// 		float Qzy = rot(2, 1) = d[7];
-// 		float Qzz = rot(2, 2) = d[8];
+// 		Eigen::Matrix3f rot;
 
-// 		cif::matrix4x4<float> m({
-// 			Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
-// 			Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
-// 			Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
-// 			Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz
-// 		});
+// 		rot << d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8];
 
-// 		auto &&[ev, em] = cif::eigen(m * (1/3.0f), false);
-
-// 		std::size_t bestJ = 0;
-// 		float bestEV = -1;
-
-// 		for (std::size_t j = 0; j < 4; ++j)
+// 		if (rot * rot.transpose() == Eigen::Matrix3f::Identity() and rot.determinant() == 1)
 // 		{
-// 			if (bestEV < ev[j])
-// 			{
-// 				bestEV = ev[j];
-// 				bestJ = j;
-// 			}
+// 			Eigen::Quaternionf qe(rot);
+// 			auto q = normalize(cif::quaternion{ qe.w(), qe.x(), qe.y(), qe.z() });
+
+// 			cif::point p1{ 1, 1, 1 };
+// 			cif::point p2 = p1;
+// 			p2.rotate(q);
+
+// 			auto p3 = rot * Eigen::Vector3f{ p1.m_x, p1.m_y, p1.m_z };
+
+// 			CHECK_THAT(p2.m_x, Catch::Matchers::WithinRel(p3[0], 0.01f));
+// 			CHECK_THAT(p2.m_y, Catch::Matchers::WithinRel(p3[1], 0.01f));
+// 			CHECK_THAT(p2.m_z, Catch::Matchers::WithinRel(p3[2], 0.01f));
 // 		}
-
-// 		if (std::abs(bestEV - 1) > 0.01)
-// 			continue; // not a rotation matrix
-
-// 		auto q = normalize(cif::quaternion{
-// 			static_cast<float>(em(bestJ, 3)),
-// 			static_cast<float>(em(bestJ, 0)),
-// 			static_cast<float>(em(bestJ, 1)),
-// 			static_cast<float>(em(bestJ, 2)) });
-
-// 		cif::point p1{ 1, 1, 1 };
-// 		cif::point p2 = p1;
-// 		p2 = q * p2;
-
-// 		cif::point p3 = rot * p1;
-
-// 		REQUIRE(p2.x == p3.x);
-// 		REQUIRE(p2.y == p3.y);
-// 		REQUIRE(p2.z == p3.z);
 // 	}
 // }
 
@@ -338,15 +268,15 @@ TEST_CASE("symm_1")
 
 	cif::point f = fractional(p, c);
 
-	CHECK_THAT(f.x, Catch::Matchers::WithinRel(0.1f, 0.01f));
-	CHECK_THAT(f.y, Catch::Matchers::WithinRel(0.1f, 0.01f));
-	CHECK_THAT(f.z, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.m_x, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.m_y, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.m_z, Catch::Matchers::WithinRel(0.1f, 0.01f));
 
 	cif::point o = orthogonal(f, c);
 
-	CHECK_THAT(o.x, Catch::Matchers::WithinRel(1.f, 0.01f));
-	CHECK_THAT(o.y, Catch::Matchers::WithinRel(1.f, 0.01f));
-	CHECK_THAT(o.z, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.m_x, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.m_y, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.m_z, Catch::Matchers::WithinRel(1.f, 0.01f));
 }
 
 TEST_CASE("symm_2")
@@ -385,9 +315,9 @@ TEST_CASE("symm_4")
 	CHECK_THAT(distance(a, sg(a, c, "1_554"_symop)), Catch::Matchers::WithinRel(static_cast<float>(c.get_c()), 0.01f));
 
 	auto sb2 = sg(b, c, "4_565"_symop);
-	CHECK_THAT(sb.x, Catch::Matchers::WithinRel(sb2.x, 0.01f));
-	CHECK_THAT(sb.y, Catch::Matchers::WithinRel(sb2.y, 0.01f));
-	CHECK_THAT(sb.z, Catch::Matchers::WithinRel(sb2.z, 0.01f));
+	CHECK_THAT(sb.m_x, Catch::Matchers::WithinRel(sb2.m_x, 0.01f));
+	CHECK_THAT(sb.m_y, Catch::Matchers::WithinRel(sb2.m_y, 0.01f));
+	CHECK_THAT(sb.m_z, Catch::Matchers::WithinRel(sb2.m_z, 0.01f));
 
 	CHECK_THAT(distance(a, sb2), Catch::Matchers::WithinRel(7.42f, 0.01f));
 }
@@ -424,17 +354,17 @@ TEST_CASE("symm_2bi3_1")
 		auto sa1 = c.symmetry_copy(a1.get_location(), cif::sym_op(symm1));
 		auto sa2 = c.symmetry_copy(a2.get_location(), cif::sym_op(symm2));
 
-		CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.01f));
+		CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
 
 		auto pa1 = a1.get_location();
 
 		const auto &[d, p, so] = c.closest_symmetry_copy(pa1, a2.get_location());
 
-		CHECK_THAT(p.x, Catch::Matchers::WithinAbs(sa2.x, 0.01f));
-		CHECK_THAT(p.y, Catch::Matchers::WithinAbs(sa2.y, 0.01f));
-		CHECK_THAT(p.z, Catch::Matchers::WithinAbs(sa2.z, 0.01f));
+		CHECK_THAT(p.m_x, Catch::Matchers::WithinAbs(sa2.m_x, 0.5f));
+		CHECK_THAT(p.m_y, Catch::Matchers::WithinAbs(sa2.m_y, 0.5f));
+		CHECK_THAT(p.m_z, Catch::Matchers::WithinAbs(sa2.m_z, 0.5f));
 
-		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.01f));
+		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.5f));
 		REQUIRE(so.string() == symm2);
 	}
 }
@@ -462,28 +392,25 @@ TEST_CASE("symm_2bi3_1a")
 			 "ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
 			 "pdbx_dist_value"))
 	{
-		auto t1 = atom_site.find1<float, float, float>(
+		cif::point p1 = atom_site.find1<float, float, float>(
 			"label_asym_id"_key == asym1 and "label_seq_id"_key == seqid1 and "auth_seq_id"_key == authseqid1 and "label_atom_id"_key == atomid1,
 			"cartn_x", "cartn_y", "cartn_z");
-		auto t2 = atom_site.find1<float, float, float>(
+		cif::point p2 = atom_site.find1<float, float, float>(
 			"label_asym_id"_key == asym2 and "label_seq_id"_key == seqid2 and "auth_seq_id"_key == authseqid2 and "label_atom_id"_key == atomid2,
 			"cartn_x", "cartn_y", "cartn_z");
 
-		cif::point p1{ std::get<0>(t1), std::get<1>(t1), std::get<2>(t1) };
-		cif::point p2{ std::get<0>(t2), std::get<1>(t2), std::get<2>(t2) };
-
 		auto sa1 = c.symmetry_copy(p1, cif::sym_op(symm1));
 		auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
 
-		CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.01f));
+		CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
 
 		const auto &[d, p, so] = c.closest_symmetry_copy(p1, p2);
 
-		CHECK_THAT(p.x, Catch::Matchers::WithinAbs(sa2.x, 0.01f));
-		CHECK_THAT(p.y, Catch::Matchers::WithinAbs(sa2.y, 0.01f));
-		CHECK_THAT(p.z, Catch::Matchers::WithinAbs(sa2.z, 0.01f));
+		CHECK_THAT(p.m_x, Catch::Matchers::WithinAbs(sa2.m_x, 0.5f));
+		CHECK_THAT(p.m_y, Catch::Matchers::WithinAbs(sa2.m_y, 0.5f));
+		CHECK_THAT(p.m_z, Catch::Matchers::WithinAbs(sa2.m_z, 0.5f));
 
-		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.01f));
+		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.5f));
 		REQUIRE(so.string() == symm2);
 	}
 }
@@ -507,11 +434,101 @@ TEST_CASE("symm_3bwh_1")
 
 			const auto &[d, p, so] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
 
-			CHECK_THAT(d, Catch::Matchers::WithinAbs(distance(a1.get_location(), p), 0.01f));
+			CHECK_THAT(d, Catch::Matchers::WithinAbs(distance(a1.get_location(), p), 0.5f));
 		}
 	}
 }
 
+TEST_CASE("symm_476d")
+{
+	cif::file f(gTestDir / "476d.cif.gz");
+	f.front().set_validator(cif::validator_factory::instance().get("mmcif_pdbx.dic"));
+
+	auto &db = f.front();
+	cif::mm::structure s(db);
+
+	cif::crystal c(db);
+
+	auto struct_conn = db["struct_conn"];
+	for (const auto &[asym1, seqid1, authseqid1, atomid1, symm1,
+			 asym2, seqid2, authseqid2, atomid2, symm2,
+			 dist] : struct_conn.find<std::string, int, std::string, std::string, std::string,
+			 std::string, int, std::string, std::string, std::string,
+			 float>(
+			 cif::key("ptnr1_symmetry") != "1_555" or cif::key("ptnr2_symmetry") != "1_555",
+			 "ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr1_auth_seq_id", "ptnr1_label_atom_id", "ptnr1_symmetry",
+			 "ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
+			 "pdbx_dist_value"))
+	{
+		auto &r1 = s.get_residue(asym1, seqid1, authseqid1);
+		auto &r2 = s.get_residue(asym2, seqid2, authseqid2);
+
+		auto a1 = r1.get_atom_by_atom_id(atomid1);
+		auto a2 = r2.get_atom_by_atom_id(atomid2);
+
+		auto p1 = a1.get_location();
+		auto p2 = a2.get_location();
+
+		cif::sym_op so1(symm1);
+		cif::sym_op so2(symm2);
+
+		auto sa1 = c.symmetry_copy(p1, so1);
+		auto sa2 = c.symmetry_copy(p2, so2);
+
+		CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.01f));
+	}
+}
+
+TEST_CASE("symm-P_32_2_1_a")
+{
+	cif::cell c{ 80, 80, 120, 90, 90, 120 };
+	cif::spacegroup sg{ 154 };
+
+	cif::crystal crystal{ c, sg };
+
+	cif::point a{ 1, 90, 1 };
+	cif::point p1{ 2, 2, 2 };
+
+	auto d = distance(a, p1);
+
+	auto [d2, p, so] = crystal.closest_symmetry_copy(a, p1);
+
+	std::cout << "d: " << d2 << " p: " << p << " so: " << so.string() << '\n';
+
+	auto p2 = crystal.symmetry_copy(p1, so);
+	auto d3 = distance(p2, a);
+
+	CHECK_THAT(cif::distance(p2, p), Catch::Matchers::WithinAbs(0.f, 0.01f));
+	CHECK_THAT(d3, Catch::Matchers::WithinAbs(d2, 0.01f));
+
+	CHECK(d2 <= d);
+}
+
+TEST_CASE("symm-P_32_2_1")
+{
+	cif::cell c{ 82.162, 82.162, 135.202, 90, 90, 120 };
+	cif::spacegroup sg{ 154 };
+
+	cif::crystal crystal{ c, sg };
+
+	cif::point a{ 1.73727,89.1813,11.1388 };
+	cif::point p1{ -8.98574, 50.3861, -11.6447 };
+
+	auto d = distance(a, p1);
+
+	auto [d2, p, so] = crystal.closest_symmetry_copy(a, p1);
+
+	std::cout << "d: " << d2 << " p: " << p << " so: " << so.string() << '\n';
+
+	auto p2 = crystal.symmetry_copy(p1, so);
+	auto d3 = distance(p2, a);
+
+	CHECK_THAT(cif::distance(p2, p), Catch::Matchers::WithinAbs(0.f, 0.01f));
+	CHECK_THAT(d3, Catch::Matchers::WithinAbs(d2, 0.01f));
+
+	CHECK(d2 <= d);
+}
+
 TEST_CASE("volume_3bwh_1")
 {
 	cif::file f(gTestDir / "1juh.cif.gz");
@@ -556,7 +573,7 @@ TEST_CASE("smallest_sphere-1")
 	for (int i = 0; i < 1000; ++i)
 	{
 		auto [c, r] = cif::smallest_sphere_around_points(pts);
-		CHECK_THAT(glm::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
+		CHECK_THAT(cif::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
 		CHECK_THAT(r, Catch::Matchers::WithinAbs(7.31248331f, 0.01f));
 	}
 }