Merge branch 'trunk' of github.com:PDB-REDO/libcifpp into trunk

improved loading resources from absolute path
better error reporting when loading dictionary
2026-06-06 15:28:04 +08:00 · 2023-05-31 15:17:00 +02:00 · 2023-05-31 15:16:10 +02:00 · 2023-05-16 13:55:07 +02:00 · 2023-05-09 11:46:37 +02:00 · 2023-05-02 13:45:02 +02:00
57 changed files with 5778 additions and 7548 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,7 +25,7 @@
 cmake_minimum_required(VERSION 3.16)

 # set the project name
-project(cifpp VERSION 5.0.8 LANGUAGES CXX)
+project(cifpp VERSION 5.0.9 LANGUAGES CXX)

 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")

@@ -35,14 +35,15 @@ include(CheckIncludeFiles)
 include(CheckLibraryExists)
 include(CMakePackageConfigHelpers)
 include(CheckCXXSourceCompiles)
+include(Dart)
 include(GenerateExportHeader)

 set(CXX_EXTENSIONS OFF)
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
 # set(CMAKE_CXX_VISIBILITY_PRESET hidden)
 # set(CMAKE_VISIBILITY_INLINES_HIDDEN 1)
-
 if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-unused-parameter -Wno-missing-field-initializers")
 elseif(MSVC)
@@ -102,11 +103,22 @@ if(MSVC)
 	add_compile_options(/permissive-)

 	macro(get_WIN32_WINNT version)
-		if(WIN32 AND CMAKE_SYSTEM_VERSION)
+		if(CMAKE_SYSTEM_VERSION)
 			set(ver ${CMAKE_SYSTEM_VERSION})
-			string(REPLACE "." "" ver ${ver})
-			string(REGEX REPLACE "([0-9])" "0\\1" ver ${ver})
+			string(REGEX MATCH "^([0-9]+).([0-9])" ver ${ver})
+			string(REGEX MATCH "^([0-9]+)" verMajor ${ver})

+			# Check for Windows 10, b/c we'll need to convert to hex 'A'.
+			if("${verMajor}" MATCHES "10")
+				set(verMajor "A")
+				string(REGEX REPLACE "^([0-9]+)" ${verMajor} ver ${ver})
+			endif()
+
+			# Remove all remaining '.' characters.
+			string(REPLACE "." "" ver ${ver})
+
+			# Prepend each digit with a zero.
+			string(REGEX REPLACE "([0-9A-Z])" "0\\1" ver ${ver})
 			set(${version} "0x${ver}")
 		endif()
 	endmacro()
@@ -157,17 +169,20 @@ if(MSVC)
 	set(_ZLIB_x86 "(x86)")
 	set(_ZLIB_SEARCH_NORMAL
 		PATHS "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Zlib;InstallPath]"
-			"$ENV{ProgramFiles}/zlib"
-			"$ENV{ProgramFiles${_ZLIB_x86}}/zlib")
+		"$ENV{ProgramFiles}/zlib"
+		"$ENV{ProgramFiles${_ZLIB_x86}}/zlib")
 	unset(_ZLIB_x86)
 	list(APPEND _ZLIB_SEARCHES _ZLIB_SEARCH_NORMAL)

 	foreach(search ${_ZLIB_SEARCHES})
-    	find_library(ZLIB_LIBRARY NAMES zlibstatic NAMES_PER_DIR ${${search}} PATH_SUFFIXES lib)
+		find_library(ZLIB_LIBRARY NAMES zlibstatic NAMES_PER_DIR ${${search}} PATH_SUFFIXES lib)
 	endforeach()
 endif()
+
 find_package(ZLIB REQUIRED)

+find_package(Eigen3 REQUIRED)
+
 include(FindFilesystem)
 list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})

@@ -191,12 +206,12 @@ if(CIFPP_RECREATE_SYMOP_DATA)

 	add_custom_command(
 		OUTPUT ${PROJECT_SOURCE_DIR}/src/symop_table_data.hpp
-		COMMAND $<TARGET_FILE:symop-map-generator> $ENV{CLIBD}/syminfo.lib ${PROJECT_SOURCE_DIR}/src/symop_table_data.hpp
+		COMMAND $<TARGET_FILE:symop-map-generator> $ENV{CLIBD}/syminfo.lib $ENV{CLIBD}/symop.lib ${PROJECT_SOURCE_DIR}/src/symop_table_data.hpp
 	)

 	add_custom_target(
 		OUTPUT ${PROJECT_SOURCE_DIR}/src/symop_table_data.hpp
-		DEPENDS symop-map-generator "$ENV{CLIBD}/syminfo.lib"
+		DEPENDS symop-map-generator "$ENV{CLIBD}/syminfo.lib" "$ENV{CLIBD}/symop.lib"
 	)
 endif()

@@ -273,7 +288,7 @@ set_target_properties(cifpp PROPERTIES POSITION_INDEPENDENT_CODE ON)

 target_include_directories(cifpp
 	PUBLIC
-	"$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${PROJECT_BINARY_DIR}>"
+	"$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${PROJECT_BINARY_DIR};${EIGEN3_INCLUDE_DIR}>"
 	"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )

@@ -416,8 +431,7 @@ if(ENABLE_TESTING)

 	find_package(Boost REQUIRED)

-	list(APPEND CIFPP_tests unit-v2 unit-3d format model rename-compound sugar
-	)
+	list(APPEND CIFPP_tests unit-v2 unit-3d format model rename-compound sugar spinner)

 	foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 		set(CIFPP_TEST "${CIFPP_TEST}-test")
@@ -436,10 +450,10 @@ if(ENABLE_TESTING)

 		add_custom_command(
 			OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/Run${CIFPP_TEST}.touch
-			COMMAND $<TARGET_FILE:${CIFPP_TEST}> -- ${PROJECT_SOURCE_DIR}/test)
+			COMMAND $<TARGET_FILE:${CIFPP_TEST}> -- ${CMAKE_CURRENT_SOURCE_DIR}/test)

 		add_test(NAME ${CIFPP_TEST}
-			COMMAND $<TARGET_FILE:${CIFPP_TEST}> -- ${PROJECT_SOURCE_DIR}/test)
+			COMMAND $<TARGET_FILE:${CIFPP_TEST}> -- ${CMAKE_CURRENT_SOURCE_DIR}/test)
 	endforeach()
 endif()

@@ -447,7 +461,7 @@ endif()
 if(CIFPP_INSTALL_UPDATE_SCRIPT)
 	if(${CMAKE_SYSTEM_NAME} STREQUAL "Linux")
 		set(CIFPP_CRON_DIR "${CIFPP_ETC_DIR}/cron.weekly")
-	elseif(UNIX)	# assume all others are like FreeBSD... 
+	elseif(UNIX) # assume all others are like FreeBSD...
 		set(CIFPP_CRON_DIR "${CIFPP_ETC_DIR}/periodic/weekly")
 	else()
 		message(FATAL_ERROR "Don't know where to install the update script")
--- a/README.md
+++ b/README.md
@@ -22,25 +22,20 @@ Building
 This library uses [cmake](https://cmake.org). The usual way of building
 and installing is to create a `build` directory and run cmake there.

-On linux e.g. you would issue the following commands:
+On linux e.g. you would issue the following commands to build and install
+libcifpp in your `$HOME/.local` folder:

+```bash
+ git clone https://github.com/PDB-REDO/libcifpp.git
+ cd libcifpp
+ cmake -S . -B build -DCMAKE_INSTALL_PREFIX=$HOME/.local -DCMAKE_BUILD_TYPE=Release
+ cmake --build build
+ cmake --install build
 ```
-	git clone https://github.com/PDB-REDO/libcifpp.git
-	cd libcifpp
-	mkdir build
-	cd build
-	cmake ..
-	cmake --build . --config Release
-	ctest -C Release
-	cmake --install .
-```
+
 This checks out the source code from github, creates a new directory
-where cmake stores its files. Run a configure, build the code and run
-tests. And then it installs the library and auxiliary files.
-
-The default is to install everything in `$HOME/.local` on Linux and
-`%LOCALAPPDATA%` on Windows (the AppData/Local folder in your home directory).
-You can change this by specifying the prefix with the
-[CMAKE_INSTALL_PREFIX](https://cmake.org/cmake/help/v3.21/variable/CMAKE_INSTALL_PREFIX.html)
-variable.
+where cmake stores its files. Run a configure, build the code and then
+it installs the library and auxiliary files.

+If you want to run the tests before installing, you should add `-DENABLE_TESTING=ON`
+to the first cmake command.
--- a/9
+++ b/9
@@ -1,9 +1,18 @@
+Version 5.0.9
+- Fix in dihedral angle calculations
+- Added create_water to model
+- Writing twin domain info in PDB files and more PDB fixes
+- remove_atom improved (remove struct_conn records)
+- Added a specialisation for category::find1<std::optional>
+- fix memory leak in category
+
 Version 5.0.8
 - implemented find_first, find_min, find_max and count in category
 - find1 now throws an exception if condition does not not exactly match one row
 - Change in writing out PDB files, now looking up the original auth_seq_num
  via the pdbx_xxx_scheme categories based on the atom_site.auth_seq_num ->
  pdbx_xxx_scheme.pdb_seq_num relationship.
+- fix memory leak in category

 Version 5.0.7.1
 - Use the implementation from zeep for std::experimental::is_detected
--- a/include/cif++.hpp
+++ b/include/cif++.hpp
@@ -26,16 +26,16 @@

 #pragma once

-#include <cif++/utilities.hpp>
-#include <cif++/file.hpp>
-#include <cif++/parser.hpp>
-#include <cif++/format.hpp>
+#include "cif++/utilities.hpp"
+#include "cif++/file.hpp"
+#include "cif++/parser.hpp"
+#include "cif++/format.hpp"

-#include <cif++/compound.hpp>
-#include <cif++/point.hpp>
-#include <cif++/symmetry.hpp>
+#include "cif++/compound.hpp"
+#include "cif++/point.hpp"
+#include "cif++/symmetry.hpp"

-#include <cif++/model.hpp>
+#include "cif++/model.hpp"

-#include <cif++/pdb/io.hpp>
-#include <cif++/gzio.hpp>
+#include "cif++/pdb/io.hpp"
+#include "cif++/gzio.hpp"
--- a/include/cif++/atom_type.hpp
+++ b/include/cif++/atom_type.hpp
@@ -270,6 +270,10 @@ class atom_type_traits
 	const SFData &wksf(int charge = 0) const;
 	const SFData &elsf() const;

+	// Clipper doesn't like atoms with charges that do not have a scattering factor. And
+	// rightly so, but we need to know in advance if this is the case
+	bool has_sf(int charge) const;
+
  private:
 	const struct atom_type_info *m_info;
 };
--- a/include/cif++/category.hpp
+++ b/include/cif++/category.hpp
@@ -26,12 +26,12 @@

 #pragma once

-#include <cif++/forward_decl.hpp>
+#include "cif++/forward_decl.hpp"

-#include <cif++/condition.hpp>
-#include <cif++/iterator.hpp>
-#include <cif++/row.hpp>
-#include <cif++/validate.hpp>
+#include "cif++/condition.hpp"
+#include "cif++/iterator.hpp"
+#include "cif++/row.hpp"
+#include "cif++/validate.hpp"

 #include <array>

@@ -63,6 +63,12 @@ class multiple_results_error : public std::runtime_error
 	}
 };

+// --------------------------------------------------------------------
+// These should be moved elsewhere, one day.
+
+template<typename _Tp> inline constexpr bool is_optional_v = false;
+template<typename _Tp> inline constexpr bool is_optional_v<std::optional<_Tp>> = true;
+
 // --------------------------------------------------------------------

 class category
@@ -299,7 +305,7 @@ class category
 		return find1<T>(cbegin(), std::move(cond), column);
 	}

-	template <typename T>
+	template <typename T, std::enable_if_t<not is_optional_v<T>, int> = 0>
 	T find1(const_iterator pos, condition &&cond, const char *column) const
 	{
 		auto h = find<T>(pos, std::move(cond), column);
@@ -310,6 +316,20 @@ class category
 		return *h.begin();
 	}

+	template <typename T, std::enable_if_t<is_optional_v<T>, int> = 0>
+	T find1(const_iterator pos, condition &&cond, const char *column) const
+	{
+		auto h = find<typename T::value_type>(pos, std::move(cond), column);
+
+		if (h.size() > 1)
+			throw multiple_results_error();
+
+		if (h.empty())
+			return {};
+
+		return *h.begin();
+	}
+
 	template <typename... Ts, typename... Cs, typename U = std::enable_if_t<sizeof...(Ts) != 1>>
 	std::tuple<Ts...> find1(condition &&cond, Cs... columns) const
 	{
--- a/include/cif++/compound.hpp
+++ b/include/cif++/compound.hpp
@@ -29,9 +29,8 @@
 /// \file This file contains the definition for the class compound, encapsulating
 /// the information found for compounds in the CCD.

-#include <cif++.hpp>
-#include <cif++/atom_type.hpp>
-#include <cif++/point.hpp>
+#include "cif++/atom_type.hpp"
+#include "cif++/point.hpp"

 #include <map>
 #include <set>
--- a/include/cif++/condition.hpp
+++ b/include/cif++/condition.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/row.hpp>
+#include "cif++/row.hpp"

 #include <cassert>
 #include <functional>
@@ -58,7 +58,7 @@ namespace detail
 		virtual void str(std::ostream &) const = 0;
 		virtual std::optional<row_handle> single() const { return {}; };

-		virtual bool equals(const condition_impl *rhs) const { return false; }
+		virtual bool equals([[maybe_unused]] const condition_impl *rhs) const { return false; }
 	};

 	struct all_condition_impl : public condition_impl
--- a/include/cif++/datablock.hpp
+++ b/include/cif++/datablock.hpp
@@ -26,8 +26,8 @@

 #pragma once

-#include <cif++/category.hpp>
-#include <cif++/forward_decl.hpp>
+#include "cif++/category.hpp"
+#include "cif++/forward_decl.hpp"

 namespace cif
 {
--- a/include/cif++/dictionary_parser.hpp
+++ b/include/cif++/dictionary_parser.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/validate.hpp>
+#include "cif++/validate.hpp"

 namespace cif
 {
--- a/include/cif++/file.hpp
+++ b/include/cif++/file.hpp
@@ -28,9 +28,9 @@

 #include <list>

-#include <cif++/exports.hpp>
-#include <cif++/datablock.hpp>
-#include <cif++/parser.hpp>
+#include "cif++/exports.hpp"
+#include "cif++/datablock.hpp"
+#include "cif++/parser.hpp"

 namespace cif
 {
--- a/include/cif++/forward_decl.hpp
+++ b/include/cif++/forward_decl.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/exports.hpp>
+#include "cif++/exports.hpp"

 #include <string>
 #include <vector>
--- a/include/cif++/item.hpp
+++ b/include/cif++/item.hpp
@@ -26,10 +26,10 @@

 #pragma once

-#include <cif++/exports.hpp>
-#include <cif++/forward_decl.hpp>
-#include <cif++/text.hpp>
-#include <cif++/utilities.hpp>
+#include "cif++/exports.hpp"
+#include "cif++/forward_decl.hpp"
+#include "cif++/text.hpp"
+#include "cif++/utilities.hpp"

 #include <cassert>
 #include <charconv>
--- a/include/cif++/iterator.hpp
+++ b/include/cif++/iterator.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/row.hpp>
+#include "cif++/row.hpp"

 #include <array>

--- a/include/cif++/matrix.hpp
+++ b/include/cif++/matrix.hpp
@@ -0,0 +1,531 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2023 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.
+ */
+
+#pragma once
+
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <ostream>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+namespace cif
+{
+// --------------------------------------------------------------------
+// We're using expression templates here
+
+template <typename M>
+class matrix_expression
+{
+  public:
+	constexpr uint32_t dim_m() const { return static_cast<const M &>(*this).dim_m(); }
+	constexpr uint32_t dim_n() const { return static_cast<const M &>(*this).dim_n(); }
+
+	constexpr auto &operator()(uint32_t i, uint32_t j)
+	{
+		return static_cast<M &>(*this).operator()(i, j);
+	}
+
+	constexpr auto operator()(uint32_t i, uint32_t j) const
+	{
+		return static_cast<const M &>(*this).operator()(i, j);
+	}
+
+	void swap_row(uint32_t r1, uint32_t r2)
+	{
+		for (uint32_t c = 0; c < dim_m(); ++c)
+		{
+			auto v = operator()(r1, c);
+			operator()(r1, c) = operator()(r2, c);
+			operator()(r2, c) = v;
+		}
+	}
+
+	void swap_col(uint32_t c1, uint32_t c2)
+	{
+		for (uint32_t r = 0; r < dim_n(); ++r)
+		{
+			auto &a = operator()(r, c1);
+			auto &b = operator()(r, c2);
+			std::swap(a, b);
+		}
+	}
+
+	friend std::ostream &operator<<(std::ostream &os, const matrix_expression &m)
+	{
+		os << '[';
+
+		for (size_t i = 0; i < m.dim_m(); ++i)
+		{
+			os << '[';
+
+			for (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;
+	}
+};
+
+// --------------------------------------------------------------------
+// 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:
+	using value_type = F;
+
+	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 (uint32_t i = 0; i < m_m; ++i)
+		{
+			for (uint32_t j = 0; j < m_n; ++j)
+				operator()(i, j) = m(i, j);
+		}
+	}
+
+	matrix(size_t m, 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);
+	}
+
+	matrix() = default;
+	matrix(matrix &&m) = default;
+	matrix(const matrix &m) = default;
+	matrix &operator=(matrix &&m) = default;
+	matrix &operator=(const matrix &m) = default;
+
+	constexpr size_t dim_m() const { return m_m; }
+	constexpr size_t dim_n() const { return m_n; }
+
+	constexpr value_type operator()(size_t i, size_t j) const
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+	constexpr value_type &operator()(size_t i, size_t j)
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+  private:
+	size_t m_m = 0, m_n = 0;
+	std::vector<value_type> m_data;
+};
+
+// --------------------------------------------------------------------
+// special case, 3x3 matrix
+
+template <typename F, size_t M, size_t N>
+class matrix_fixed : public matrix_expression<matrix_fixed<F, M, N>>
+{
+  public:
+	using value_type = F;
+
+	static constexpr size_t kSize = M * N;
+
+	template <typename M2>
+	matrix_fixed(const M2 &m)
+	{
+		assert(M == m.dim_m() and N == m.dim_n());
+		for (uint32_t i = 0; i < M; ++i)
+		{
+			for (uint32_t j = 0; j < N; ++j)
+				operator()(i, j) = m(i, j);
+		}
+	}
+
+	matrix_fixed(value_type v = 0)
+	{
+		m_data.fill(v);
+	}
+
+	matrix_fixed(const F (&v)[kSize])
+	{
+		fill(v, std::make_index_sequence<kSize>{});
+	}
+
+	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;
+
+	template<size_t... Ixs>
+	matrix_fixed& fill(const F (&a)[kSize], std::index_sequence<Ixs...>)
+	{
+		m_data = { a[Ixs]... };
+		return *this;
+	}
+
+	constexpr size_t dim_m() const { return M; }
+	constexpr size_t dim_n() const { return N; }
+
+	constexpr value_type operator()(size_t i, size_t j) const
+	{
+		assert(i < M);
+		assert(j < N);
+		return m_data[i * N + j];
+	}
+
+	constexpr value_type &operator()(size_t i, size_t j)
+	{
+		assert(i < M);
+		assert(j < N);
+		return m_data[i * N + j];
+	}
+
+  private:
+	std::array<value_type, M * N> m_data;
+};
+
+template <typename F>
+using matrix3x3 = matrix_fixed<F, 3, 3>;
+
+template <typename F>
+using matrix4x4 = matrix_fixed<F, 4, 4>;
+
+// --------------------------------------------------------------------
+
+template <typename F = float>
+class symmetric_matrix : public matrix_expression<symmetric_matrix<F>>
+{
+  public:
+	using value_type = F;
+
+	symmetric_matrix(uint32_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);
+	}
+
+	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;
+
+	constexpr uint32_t dim_m() const { return m_n; }
+	constexpr uint32_t dim_n() const { return m_n; }
+
+	constexpr value_type operator()(uint32_t i, uint32_t j) const
+	{
+		return i < j
+		           ? m_data[(j * (j + 1)) / 2 + i]
+		           : m_data[(i * (i + 1)) / 2 + j];
+	}
+
+	constexpr value_type &operator()(uint32_t i, uint32_t j)
+	{
+		if (i > j)
+			std::swap(i, j);
+		assert(j < m_n);
+		return m_data[(j * (j + 1)) / 2 + i];
+	}
+
+  private:
+	uint32_t m_n;
+	std::vector<value_type> m_data;
+};
+
+// --------------------------------------------------------------------
+
+template <typename F, size_t M>
+class symmetric_matrix_fixed : public matrix_expression<symmetric_matrix_fixed<F, M>>
+{
+  public:
+	using value_type = F;
+
+	symmetric_matrix_fixed(value_type v = 0)
+	{
+		std::fill(m_data.begin(), m_data.end(), v);
+	}
+
+	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;
+
+	constexpr uint32_t dim_m() const { return M; }
+	constexpr uint32_t dim_n() const { return M; }
+
+	constexpr value_type operator()(uint32_t i, uint32_t j) const
+	{
+		return i < j
+		           ? m_data[(j * (j + 1)) / 2 + i]
+		           : m_data[(i * (i + 1)) / 2 + j];
+	}
+
+	constexpr value_type &operator()(uint32_t i, uint32_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;
+};
+
+template <typename F>
+using symmetric_matrix3x3 = symmetric_matrix_fixed<F, 3>;
+
+template <typename F>
+using symmetric_matrix4x4 = symmetric_matrix_fixed<F, 4>;
+
+// --------------------------------------------------------------------
+
+template <typename F = float>
+class identity_matrix : public matrix_expression<identity_matrix<F>>
+{
+  public:
+	using value_type = F;
+
+	identity_matrix(uint32_t n)
+		: m_n(n)
+	{
+	}
+
+	constexpr uint32_t dim_m() const { return m_n; }
+	constexpr uint32_t dim_n() const { return m_n; }
+
+	constexpr value_type operator()(uint32_t i, uint32_t j) const
+	{
+		return i == j ? 1 : 0;
+	}
+
+  private:
+	uint32_t m_n;
+};
+
+// --------------------------------------------------------------------
+// matrix functions, implemented as expression templates
+
+template <typename M1, typename M2>
+class matrix_subtraction : public matrix_expression<matrix_subtraction<M1, M2>>
+{
+  public:
+	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());
+	}
+
+	constexpr uint32_t dim_m() const { return m_m1.dim_m(); }
+	constexpr uint32_t dim_n() const { return m_m1.dim_n(); }
+
+	constexpr auto operator()(uint32_t i, uint32_t j) const
+	{
+		return m_m1(i, j) - m_m2(i, j);
+	}
+
+  private:
+	const M1 &m_m1;
+	const M2 &m_m2;
+};
+
+template <typename M1, typename M2>
+auto operator-(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
+{
+	return matrix_subtraction(m1, m2);
+}
+
+template <typename M1, typename M2>
+class matrix_matrix_multiplication : public matrix_expression<matrix_matrix_multiplication<M1, M2>>
+{
+  public:
+	matrix_matrix_multiplication(const M1 &m1, const M2 &m2)
+		: m_m1(m1)
+		, m_m2(m2)
+	{
+		assert(m1.dim_m() == m2.dim_n());
+	}
+
+	constexpr uint32_t dim_m() const { return m_m1.dim_m(); }
+	constexpr uint32_t dim_n() const { return m_m1.dim_n(); }
+
+	constexpr auto operator()(uint32_t i, uint32_t j) const
+	{
+		using value_type = decltype(m_m1(0, 0));
+
+		value_type result = {};
+
+		for (uint32_t k = 0; k < m_m1.dim_m(); ++k)
+			result += m_m1(i, k) * m_m2(k, j);
+
+		return result;
+	}
+
+  private:
+	const M1 &m_m1;
+	const M2 &m_m2;
+};
+
+template <typename M, typename T>
+class matrix_scalar_multiplication : public matrix_expression<matrix_scalar_multiplication<M, T>>
+{
+  public:
+	using value_type = T;
+
+	matrix_scalar_multiplication(const M &m, value_type v)
+		: m_m(m)
+		, m_v(v)
+	{
+	}
+
+	constexpr uint32_t dim_m() const { return m_m.dim_m(); }
+	constexpr uint32_t dim_n() const { return m_m.dim_n(); }
+
+	constexpr auto operator()(uint32_t i, uint32_t j) const
+	{
+		return m_m(i, j) * m_v;
+	}
+
+  private:
+	const M &m_m;
+	value_type m_v;
+};
+
+template <typename M1, typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
+auto operator*(const matrix_expression<M1> &m, T v)
+{
+	return matrix_scalar_multiplication(m, v);
+}
+
+template <typename M1, typename M2, std::enable_if_t<not std::is_floating_point_v<M2>, int> = 0>
+auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
+{
+	return matrix_matrix_multiplication(m1, m2);
+}
+
+// --------------------------------------------------------------------
+
+template <typename M>
+auto determinant(const M &m);
+
+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)));
+}
+
+template <typename M>
+M inverse(const M &m);
+
+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;
+}
+
+// --------------------------------------------------------------------
+
+template <typename M>
+class matrix_cofactors : public matrix_expression<matrix_cofactors<M>>
+{
+  public:
+	matrix_cofactors(const M &m)
+		: m_m(m)
+	{
+	}
+
+	constexpr uint32_t dim_m() const { return m_m.dim_m(); }
+	constexpr uint32_t dim_n() const { return m_m.dim_n(); }
+
+	constexpr auto operator()(uint32_t i, uint32_t j) const
+	{
+		const size_t ixs[4][3] = {
+			{ 1, 2, 3 },
+			{ 0, 2, 3 },
+			{ 0, 1, 3 },
+			{ 0, 1, 2 }
+		};
+
+		const size_t *ix = ixs[i];
+		const 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
--- a/include/cif++/model.hpp
+++ b/include/cif++/model.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/atom_type.hpp>
+#include "cif++/atom_type.hpp"

 #include <numeric>

@@ -216,6 +216,14 @@ class atom
 		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.rotate(q, p);
+		set_location(loc);
+	}
+
 	/// \brief Translate and rotate the position of this atom by \a t and \a q
 	void translate_and_rotate(point t, quaternion q)
 	{
@@ -653,37 +661,6 @@ class branch : public std::vector<sugar>
 	std::string m_asym_id, m_entity_id;
 };

-// // --------------------------------------------------------------------
-// // file is a reference to the data stored in e.g. the cif file.
-// // This object is not copyable.
-
-// class File : public file
-// {
-//   public:
-// 	File() {}
-
-// 	// File(const std::filesystem::path &path)
-// 	// {
-// 	// 	load(path);
-// 	// }
-
-// 	// File(const char *data, size_t length)
-// 	// {
-// 	// 	load(data, length);
-// 	// }
-
-// 	File(const File &) = delete;
-// 	File &operator=(const File &) = delete;
-
-// 	// void load(const std::filesystem::path &p) override;
-// 	// void save(const std::filesystem::path &p) override;
-
-// 	// using file::load;
-// 	// using file::save;
-
-// 	datablock &data() { return front(); }
-// };
-
 // --------------------------------------------------------------------

 enum class StructureOpenOptions
@@ -691,7 +668,7 @@ enum class StructureOpenOptions
 	SkipHydrogen = 1 << 0
 };

-inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
+constexpr inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
 {
 	return static_cast<int>(a) bitand static_cast<int>(b);
 }
@@ -744,6 +721,7 @@ class structure

 	const std::vector<residue> &non_polymers() const { return m_non_polymers; }

+	bool has_atom_id(const std::string &id) const;
 	atom get_atom_by_id(const std::string &id) const;
 	// atom getAtomByLocation(point pt, float maxDistance) const;

@@ -839,6 +817,12 @@ class structure
 	/// \return				The newly create asym ID
 	std::string create_non_poly(const std::string &entity_id, std::vector<row_initializer> atoms);

+	/// \brief Create a new water with atom constructed from info in \a atom_info
+	/// This method creates a new atom record filled with info from the info.
+	///
+	/// \param atom			The set of item data containing the data for the atoms.
+	void create_water(row_initializer atom);
+
 	/// \brief Create a new and empty (sugar) branch
 	branch &create_branch();

--- a/include/cif++/parser.hpp
+++ b/include/cif++/parser.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/row.hpp>
+#include "cif++/row.hpp"

 #include <map>
 #include <regex>
--- a/include/cif++/pdb/cif2pdb.hpp
+++ b/include/cif++/pdb/cif2pdb.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++.hpp>
+#include "cif++/datablock.hpp"

 namespace cif::pdb
 {
--- a/include/cif++/pdb/io.hpp
+++ b/include/cif++/pdb/io.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++.hpp>
+#include "cif++/datablock.hpp"

 namespace cif::pdb
 {
--- a/include/cif++/pdb/pdb2cif.hpp
+++ b/include/cif++/pdb/pdb2cif.hpp
@@ -26,8 +26,7 @@

 #pragma once

-#include <cif++.hpp>
-
+#include "cif++/file.hpp"
 namespace cif::pdb
 {

--- a/include/cif++/pdb/pdb2cif_remark_3.hpp
+++ b/include/cif++/pdb/pdb2cif_remark_3.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/pdb/pdb2cif.hpp>
+#include "pdb2cif.hpp"

 // --------------------------------------------------------------------

--- a/include/cif++/pdb/tls.hpp
+++ b/include/cif++/pdb/tls.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++.hpp>
+#include "cif++/datablock.hpp"

 #include <string>
 #include <tuple>
--- a/include/cif++/point.hpp
+++ b/include/cif++/point.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/exports.hpp>
+#include "cif++/exports.hpp"

 #include <array>
 #include <cmath>
@@ -320,7 +320,7 @@ class quaternion_type

 	constexpr operator bool() const
 	{
-		return operator!=({});
+		return a != 0 or b != 0 or c != 0 or d != 0;
 	}

  private:
@@ -718,7 +718,6 @@ template <int N>
 class spherical_dots
 {
  public:
-
 	constexpr static int P = 2 * N * 1;

 	using array_type = typename std::array<point, P>;
--- a/include/cif++/row.hpp
+++ b/include/cif++/row.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/item.hpp>
+#include "cif++/item.hpp"

 #include <array>

@@ -112,12 +112,12 @@ class row : public std::vector<item_value>

 	item_value* get(uint16_t ix)
 	{
-		return ix < size() ? &at(ix) : nullptr;
+		return ix < size() ? &data()[ix] : nullptr;
 	}

 	const item_value* get(uint16_t ix) const
 	{
-		return ix < size() ? &at(ix) : nullptr;
+		return ix < size() ? &data()[ix] : nullptr;
 	}

  private:
@@ -210,14 +210,14 @@ class row_handle
 		return detail::get_row_result<C...>(*this, { get_column_ix(columns)... });
 	}

-	template <typename... Ts, typename... C, std::enable_if_t<sizeof...(Ts) == sizeof...(C), int> = 0>
+	template <typename... Ts, typename... C, std::enable_if_t<sizeof...(Ts) == sizeof...(C) and sizeof...(C) != 1, int> = 0>
 	std::tuple<Ts...> get(C... columns) const
 	{
 		return detail::get_row_result<Ts...>(*this, { get_column_ix(columns)... });
 	}

 	template <typename T>
-	T get(const char *column)
+	T get(const char *column) const
 	{
 		return operator[](get_column_ix(column)).template as<T>();
 	}
--- a/include/cif++/symmetry.hpp
+++ b/include/cif++/symmetry.hpp
@@ -26,17 +26,34 @@

 #pragma once

-#include <cif++/exports.hpp>
+#include "cif++/exports.hpp"
+#include "cif++/matrix.hpp"
+#include "cif++/point.hpp"

 #include <array>
 #include <cstdint>
 #include <string>

+/// \file cif++/symmetry.hpp
+/// This file contains code to do symmetry operations based on the
+/// operations as specified in the International Tables.
+
 namespace cif
 {

 // --------------------------------------------------------------------

+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
+	};
+}
+
+// --------------------------------------------------------------------
+
 enum class space_group_name
 {
 	full,
@@ -60,21 +77,21 @@ extern CIFPP_EXPORT const std::size_t kNrOfSpaceGroups;
 struct symop_data
 {
 	constexpr symop_data(const std::array<int, 15> &data)
-		: m_packed((data[0] & 0x03ULL) << 34 bitor
-				   (data[1] & 0x03ULL) << 32 bitor
-				   (data[2] & 0x03ULL) << 30 bitor
-				   (data[3] & 0x03ULL) << 28 bitor
-				   (data[4] & 0x03ULL) << 26 bitor
-				   (data[5] & 0x03ULL) << 24 bitor
-				   (data[6] & 0x03ULL) << 22 bitor
-				   (data[7] & 0x03ULL) << 20 bitor
-				   (data[8] & 0x03ULL) << 18 bitor
-				   (data[9] & 0x07ULL) << 15 bitor
-				   (data[10] & 0x07ULL) << 12 bitor
-				   (data[11] & 0x07ULL) << 9 bitor
-				   (data[12] & 0x07ULL) << 6 bitor
-				   (data[13] & 0x07ULL) << 3 bitor
-				   (data[14] & 0x07ULL) << 0)
+		: m_packed((data[0] bitand 0x03ULL) << 34 bitor
+				   (data[1] bitand 0x03ULL) << 32 bitor
+				   (data[2] bitand 0x03ULL) << 30 bitor
+				   (data[3] bitand 0x03ULL) << 28 bitor
+				   (data[4] bitand 0x03ULL) << 26 bitor
+				   (data[5] bitand 0x03ULL) << 24 bitor
+				   (data[6] bitand 0x03ULL) << 22 bitor
+				   (data[7] bitand 0x03ULL) << 20 bitor
+				   (data[8] bitand 0x03ULL) << 18 bitor
+				   (data[9] bitand 0x07ULL) << 15 bitor
+				   (data[10] bitand 0x07ULL) << 12 bitor
+				   (data[11] bitand 0x07ULL) << 9 bitor
+				   (data[12] bitand 0x07ULL) << 6 bitor
+				   (data[13] bitand 0x07ULL) << 3 bitor
+				   (data[14] bitand 0x07ULL) << 0)
 	{
 	}

@@ -88,24 +105,35 @@ struct symop_data
 		return m_packed < rhs.m_packed;
 	}

-	std::array<int, 15> data() const
+	inline constexpr int unpack3(int offset) const
+	{
+		int result = (m_packed >> offset) bitand 0x03;
+		return result == 3 ? -1 : result;
+	}
+
+	inline constexpr int unpack7(int offset) const
+	{
+		return (m_packed >> offset) bitand 0x07;
+	}
+
+	constexpr std::array<int, 15> data() const
 	{
 		return {
-			static_cast<int>(m_packed >> 34) bitand 0x03,
-			static_cast<int>(m_packed >> 32) bitand 0x03,
-			static_cast<int>(m_packed >> 30) bitand 0x03,
-			static_cast<int>(m_packed >> 28) bitand 0x03,
-			static_cast<int>(m_packed >> 26) bitand 0x03,
-			static_cast<int>(m_packed >> 24) bitand 0x03,
-			static_cast<int>(m_packed >> 22) bitand 0x03,
-			static_cast<int>(m_packed >> 20) bitand 0x03,
-			static_cast<int>(m_packed >> 18) bitand 0x03,
-			static_cast<int>(m_packed >> 15) bitand 0x07,
-			static_cast<int>(m_packed >> 12) bitand 0x07,
-			static_cast<int>(m_packed >> 9) bitand 0x07,
-			static_cast<int>(m_packed >> 6) bitand 0x07,
-			static_cast<int>(m_packed >> 3) bitand 0x07,
-			static_cast<int>(m_packed >> 0) bitand 0x07,
+			unpack3(34),
+			unpack3(32),
+			unpack3(30),
+			unpack3(28),
+			unpack3(26),
+			unpack3(24),
+			unpack3(22),
+			unpack3(20),
+			unpack3(18),
+			unpack7(15),
+			unpack7(12),
+			unpack7(9),
+			unpack7(6),
+			unpack7(3),
+			unpack7(0)
 		};
 	}

@@ -125,8 +153,8 @@ struct symop_data
 struct symop_datablock
 {
 	constexpr symop_datablock(int spacegroup, int rotational_number, const std::array<int, 15> &rt_data)
-		: m_v((spacegroup & 0xffffULL) << 48 bitor
-			  (rotational_number & 0xffULL) << 40 bitor
+		: m_v((spacegroup bitand 0xffffULL) << 48 bitor
+			  (rotational_number bitand 0xffULL) << 40 bitor
 			  symop_data(rt_data).m_packed)
 	{
 	}
@@ -145,8 +173,246 @@ extern CIFPP_EXPORT const symop_datablock kSymopNrTable[];
 extern CIFPP_EXPORT const std::size_t kSymopNrTableSize;

 // --------------------------------------------------------------------
+// Some more symmetry related stuff here.

-int get_space_group_number(std::string spacegroup);                        // alternative for clipper's parsing code, using space_group_name::full
-int get_space_group_number(std::string spacegroup, space_group_name type); // alternative for clipper's parsing code
+class datablock;
+
+class cell;
+class spacegroup;
+class rtop;
+class sym_op;
+
+
+/// @brief A class that encapsulates the symmetry operations as used in PDB files, i.e. a rotational number and a translation vector
+/// The syntax in string format follows the syntax as used in mmCIF files, i.e. rotational number followed by underscore and the
+/// three translations where 5 is no movement.
+struct sym_op
+{
+  public:
+	sym_op(uint8_t nr = 1, uint8_t ta = 5, uint8_t tb = 5, uint8_t tc = 5)
+		: m_nr(nr)
+		, m_ta(ta)
+		, m_tb(tb)
+		, m_tc(tc)
+	{
+	}
+
+	explicit sym_op(std::string_view s);
+
+	sym_op(const sym_op &) = default;
+	sym_op(sym_op &&) = default;
+	sym_op &operator=(const sym_op &) = default;
+	sym_op &operator=(sym_op &&) = default;
+
+	constexpr bool is_identity() const
+	{
+		return m_nr == 1 and m_ta == 5 and m_tb == 5 and m_tc == 5;
+	}
+
+	explicit constexpr operator bool() const
+	{
+		return not is_identity();
+	}
+
+	std::string string() const;
+
+#if defined(__cpp_impl_three_way_comparison)
+	constexpr auto operator<=>(const sym_op &rhs) const = default;
+#else
+	constexpr bool operator==(const sym_op &rhs) const
+	{
+		return m_nr == rhs.m_nr and m_ta == rhs.m_ta and m_tb == rhs.m_tb and m_tc == rhs.m_tc;
+	}
+
+	constexpr bool operator!=(const sym_op &rhs) const
+	{
+		return not operator==(rhs);
+	}
+#endif
+
+	uint8_t m_nr;
+	uint8_t m_ta, m_tb, m_tc;
+};
+
+static_assert(sizeof(sym_op) == 4, "Sym_op should be four bytes");
+
+namespace literals
+{
+	inline sym_op operator""_symop(const char *text, size_t length)
+	{
+		return sym_op({ text, length });
+	}
+} // namespace literals
+
+// --------------------------------------------------------------------
+// The transformation class
+
+class transformation
+{
+  public:
+	transformation(const symop_data &data);
+	transformation(const matrix3x3<float> &r, const cif::point &t);
+
+	transformation(const transformation &) = default;
+	transformation(transformation &&) = default;
+	transformation &operator=(const transformation &) = default;
+	transformation &operator=(transformation &&) = default;
+
+	point operator()(point pt) const
+	{
+		if (m_q)
+			pt.rotate(m_q);
+		else
+			pt = m_rotation * pt;
+
+		return pt + m_translation;
+	}
+
+	friend transformation operator*(const transformation &lhs, const transformation &rhs);
+	friend transformation inverse(const transformation &t);
+
+	transformation operator-() const
+	{
+		return inverse(*this);
+	}
+
+	friend class spacegroup;
+
+  private:
+
+	// Most rotation matrices provided by the International Tables
+	// are really rotation matrices, in those cases we can construct
+	// a quaternion. Unfortunately, that doesn't work for all of them
+
+	void try_create_quaternion();
+
+	matrix3x3<float> m_rotation;
+	quaternion m_q;
+	point m_translation;
+};
+
+// --------------------------------------------------------------------
+// class cell
+
+class cell
+{
+  public:
+	cell(float a, float b, float c, float alpha = 90.f, float beta = 90.f, float gamma = 90.f);
+	cell(const datablock &db);
+
+	float get_a() const { return m_a; }
+	float get_b() const { return m_b; }
+	float get_c() const { return m_c; }
+
+	float get_alpha() const { return m_alpha; }
+	float get_beta() const { return m_beta; }
+	float get_gamma() const { return m_gamma; }
+
+	matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; }
+	matrix3x3<float> get_fractional_matrix() const { return m_fractional; }
+
+  private:
+	void init();
+
+	float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
+	matrix3x3<float> m_orthogonal, m_fractional;
+};
+
+// --------------------------------------------------------------------
+
+int get_space_group_number(const datablock &db);
+int get_space_group_number(std::string_view spacegroup);
+int get_space_group_number(std::string_view spacegroup, space_group_name type);
+
+class spacegroup : public std::vector<transformation>
+{
+  public:
+	spacegroup(const datablock &db)
+		: spacegroup(get_space_group_number(db))
+	{
+
+	}
+
+	spacegroup(std::string_view name)
+		: spacegroup(get_space_group_number(name))
+	{
+	}
+
+	spacegroup(std::string_view name, space_group_name type)
+		: spacegroup(get_space_group_number(name, type))
+	{
+	}
+
+	spacegroup(int nr);
+
+	int get_nr() const { return m_nr; }
+	std::string get_name() const;
+
+	point operator()(const point &pt, const cell &c, sym_op symop) const;
+
+	point inverse(const point &pt, const cell &c, sym_op symop) const;
+
+  private:
+	int m_nr;
+	size_t m_index;
+};
+
+// --------------------------------------------------------------------
+// A crystal combines a cell and a spacegroup.
+
+class crystal
+{
+  public:
+	crystal(const datablock &db)
+		: m_cell(db)
+		, m_spacegroup(db)
+	{
+	}
+
+	crystal(const cell &c, const spacegroup &sg)
+		: m_cell(c)
+		, m_spacegroup(sg)
+	{
+	}
+
+	crystal(const crystal &) = default;
+	crystal(crystal &&) = default;
+	crystal &operator=(const crystal &) = default;
+	crystal &operator=(crystal &&) = default;
+
+	const cell &get_cell() const { return m_cell; }
+	const spacegroup &get_spacegroup() const { return m_spacegroup; }
+
+	point symmetry_copy(const point &pt, sym_op symop) const
+	{
+		return m_spacegroup(pt, m_cell, symop);
+	}
+
+	point inverse_symmetry_copy(const point &pt, sym_op symop) const
+	{
+		return m_spacegroup.inverse(pt, m_cell, symop);
+	}
+	
+	std::tuple<float,point,sym_op> closest_symmetry_copy(point a, point b) const;
+
+  private:
+	cell m_cell;
+	spacegroup m_spacegroup;
+};
+
+// --------------------------------------------------------------------
+// Symmetry operations on points
+
+inline point orthogonal(const point &pt, const cell &c)
+{
+	return c.get_orthogonal_matrix() * pt;
+}
+
+inline point fractional(const point &pt, const cell &c)
+{
+	return c.get_fractional_matrix() * pt;
+}
+
+// --------------------------------------------------------------------

 } // namespace cif
--- a/include/cif++/text.hpp
+++ b/include/cif++/text.hpp
@@ -26,10 +26,11 @@

 #pragma once

-#include <cif++/exports.hpp>
+#include "cif++/exports.hpp"

 #include <charconv>
 #include <cmath>
+#include <cstdint>
 #include <limits>
 #include <set>
 #include <sstream>
--- a/include/cif++/utilities.hpp
+++ b/include/cif++/utilities.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/exports.hpp>
+#include "cif++/exports.hpp"

 #include <filesystem>

@@ -157,11 +157,11 @@ inline auto coloured(std::basic_string<CharT, Traits, Alloc> &s, StringColour fo
 // --------------------------------------------------------------------
 //	A progress bar

-class Progress
+class progress_bar
 {
  public:
-	Progress(int64_t inMax, const std::string &inAction);
-	virtual ~Progress();
+	progress_bar(int64_t inMax, const std::string &inAction);
+	~progress_bar();

 	void consumed(int64_t inConsumed); // consumed is relative
 	void progress(int64_t inProgress); // progress is absolute
@@ -169,10 +169,10 @@ class Progress
 	void message(const std::string &inMessage);

  private:
-	Progress(const Progress &) = delete;
-	Progress &operator=(const Progress &) = delete;
+	progress_bar(const progress_bar &) = delete;
+	progress_bar &operator=(const progress_bar &) = delete;

-	struct ProgressImpl *m_impl;
+	struct progress_bar_impl *m_impl;
 };

 // --------------------------------------------------------------------
--- a/include/cif++/validate.hpp
+++ b/include/cif++/validate.hpp
@@ -26,7 +26,7 @@

 #pragma once

-#include <cif++/text.hpp>
+#include "cif++/text.hpp"

 #include <filesystem>
 #include <list>
--- a/libcifpp.pc.in
+++ b/libcifpp.pc.in
@@ -8,6 +8,6 @@ Name: libcifpp
 Description: C++ library for the manipulation of mmCIF files.
 Version: @PACKAGE_VERSION@

-Requires.private: zlib
+Requires: zlib
 Libs: -L${libdir} -lcifpp
 Cflags: -I${includedir} -pthread
--- a/src/atom_type.cpp
+++ b/src/atom_type.cpp
@@ -26,8 +26,7 @@

 #include <cmath>

-#include <cif++.hpp>
-#include <cif++/atom_type.hpp>
+#include "cif++.hpp"

 namespace cif
 {
@@ -1078,6 +1077,26 @@ bool atom_type_traits::is_metal(const std::string& symbol)
 	return result;
 }

+bool atom_type_traits::has_sf(int charge) const
+{
+	auto type = m_info->type;
+	if (type == D)
+		type = H;
+
+	bool result = false;
+
+	for (auto& sf: data::kWKSFData)
+	{
+		if (sf.symbol == type and sf.charge == charge)
+		{
+			result = true;
+			break;
+		}
+	}
+
+	return result;
+}
+
 auto atom_type_traits::wksf(int charge) const -> const SFData&
 {
 	auto type = m_info->type;
--- a/src/category.cpp
+++ b/src/category.cpp
@@ -24,10 +24,10 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/category.hpp>
-#include <cif++/datablock.hpp>
-#include <cif++/parser.hpp>
-#include <cif++/utilities.hpp>
+#include "cif++/category.hpp"
+#include "cif++/datablock.hpp"
+#include "cif++/parser.hpp"
+#include "cif++/utilities.hpp"

 #include <numeric>
 #include <stack>
@@ -51,7 +51,7 @@ class row_comparator
 	{
 		auto cv = cat.get_cat_validator();

-		for (auto k : cv->m_keys)
+		for (auto &k : cv->m_keys)
 		{
 			uint16_t ix = cat.add_column(k);

@@ -78,13 +78,8 @@ class row_comparator
 		row_handle rhb(m_category, *b);

 		int d = 0;
-		for (auto &c : m_comparator)
+		for (const auto &[k, f] : m_comparator)
 		{
-			uint16_t k;
-			compareFunc f;
-
-			std::tie(k, f) = c;
-
 			std::string_view ka = rha[k].text();
 			std::string_view kb = rhb[k].text();

@@ -103,29 +98,30 @@ class row_comparator

 		row_handle rhb(m_category, *b);

-		int d = 0, i = 0;
-		for (auto &c : m_comparator)
+		int d = 0;
+		auto ai = a.begin();
+
+		for (const auto &[k, f] : m_comparator)
 		{
-			uint16_t k;
-			compareFunc f;
+			assert(ai != a.end());

-			std::tie(k, f) = c;
-
-			std::string_view ka = a[i++].value();
+			std::string_view ka = ai->value();
 			std::string_view kb = rhb[k].text();

 			d = f(ka, kb);

 			if (d != 0)
 				break;
+			
+			++ai;
 		}

 		return d;
 	}

  private:
-	typedef std::function<int(std::string_view, std::string_view)> compareFunc;
-	typedef std::tuple<uint16_t, compareFunc> key_comparator;
+	using compareFunc = std::function<int(std::string_view, std::string_view)>;
+	using key_comparator = std::tuple<uint16_t, compareFunc>;

 	std::vector<key_comparator> m_comparator;
 	category &m_category;
@@ -139,13 +135,7 @@ class row_comparator
 class category_index
 {
  public:
-	category_index(category *cat)
-		: m_category(*cat)
-		, m_row_comparator(m_category)
-		, m_root(nullptr)
-	{
-		reconstruct();
-	}
+	category_index(category *cat);

 	~category_index()
 	{
@@ -158,9 +148,6 @@ class category_index
 	void insert(row *r);
 	void erase(row *r);

-	// batch create
-	void reconstruct();
-
 	// reorder the row's and returns new head and tail
 	std::tuple<row *, row *> reorder()
 	{
@@ -241,7 +228,7 @@ class category_index
 			h->m_right->m_red = not h->m_right->m_red;
 	}

-	bool is_red(entry *h) const
+	constexpr bool is_red(entry *h) const
 	{
 		return h != nullptr and h->m_red;
 	}
@@ -342,6 +329,15 @@ class category_index
 	entry *m_root;
 };

+category_index::category_index(category *cat)
+	: m_category(*cat)
+	, m_row_comparator(m_category)
+	, m_root(nullptr)
+{
+	for (auto r : m_category)
+		insert(r.get_row());
+}
+
 row *category_index::find(row *k) const
 {
 	const entry *r = m_root;
@@ -482,83 +478,6 @@ category_index::entry *category_index::erase(entry *h, row *k)
 	return fix_up(h);
 }

-void category_index::reconstruct()
-{
-	delete m_root;
-	m_root = nullptr;
-
-	for (auto r : m_category)
-		insert(r.get_row());
-
-	// maybe reconstruction can be done quicker by using the following commented code.
-	// however, I've not had the time to think of a way to set the red/black flag correctly in that case.
-
-	//	std::vector<row*> rows;
-	//	transform(mCat.begin(), mCat.end(), backInserter(rows),
-	//		[](Row r) -> row* { assert(r.mData); return r.mData; });
-	//
-	//	assert(std::find(rows.begin(), rows.end(), nullptr) == rows.end());
-	//
-	//	// don't use sort here, it will run out of the stack of something.
-	//	// quicksort is notorious for using excessive recursion.
-	//	// Besides, most of the time, the data is ordered already anyway.
-	//
-	//	stable_sort(rows.begin(), rows.end(), [this](row* a, row* b) -> bool { return this->mComp(a, b) < 0; });
-	//
-	//	for (size_t i = 0; i < rows.size() - 1; ++i)
-	//		assert(mComp(rows[i], rows[i + 1]) < 0);
-	//
-	//	deque<entry*> e;
-	//	transform(rows.begin(), rows.end(), back_inserter(e),
-	//		[](row* r) -> entry* { return new entry(r); });
-	//
-	//	while (e.size() > 1)
-	//	{
-	//		deque<entry*> ne;
-	//
-	//		while (not e.empty())
-	//		{
-	//			entry* a = e.front();
-	//			e.pop_front();
-	//
-	//			if (e.empty())
-	//				ne.push_back(a);
-	//			else
-	//			{
-	//				entry* b = e.front();
-	//				b->mLeft = a;
-	//
-	//				assert(mComp(a->mRow, b->mRow) < 0);
-	//
-	//				e.pop_front();
-	//
-	//				if (not e.empty())
-	//				{
-	//					entry* c = e.front();
-	//					e.pop_front();
-	//
-	//					assert(mComp(b->mRow, c->mRow) < 0);
-	//
-	//					b->mRight = c;
-	//				}
-	//
-	//				ne.push_back(b);
-	//
-	//				if (not e.empty())
-	//				{
-	//					ne.push_back(e.front());
-	//					e.pop_front();
-	//				}
-	//			}
-	//		}
-	//
-	//		swap (e, ne);
-	//	}
-	//
-	//	assert(e.size() == 1);
-	//	mRoot = e.front();
-}
-
 size_t category_index::size() const
 {
 	std::stack<entry *> s;
@@ -600,7 +519,7 @@ category::category(const category &rhs)
 	for (auto r = rhs.m_head; r != nullptr; r = r->m_next)
 		insert_impl(end(), clone_row(*r));

-	if (m_cat_validator != nullptr)
+	if (m_cat_validator != nullptr and m_index == nullptr)
 		m_index = new category_index(this);
 }

@@ -644,7 +563,7 @@ category &category::operator=(const category &rhs)
 		m_validator = rhs.m_validator;
 		m_cat_validator = rhs.m_cat_validator;

-		if (m_cat_validator != nullptr)
+		if (m_cat_validator != nullptr and m_index == nullptr)
 			m_index = new category_index(this);
 	}

@@ -655,9 +574,6 @@ category &category::operator=(category &&rhs)
 {
 	if (this != &rhs)
 	{
-		if (not empty())
-			clear();
-
 		m_name = std::move(rhs.m_name);
 		m_columns = std::move(rhs.m_columns);
 		m_cascade = rhs.m_cascade;
@@ -665,12 +581,10 @@ category &category::operator=(category &&rhs)
 		m_cat_validator = rhs.m_cat_validator;
 		m_parent_links = rhs.m_parent_links;
 		m_child_links = rhs.m_child_links;
-		m_index = rhs.m_index;
-		m_head = rhs.m_head;
-		m_tail = rhs.m_tail;

-		rhs.m_head = rhs.m_tail = nullptr;
-		rhs.m_index = nullptr;
+		std::swap(m_index, rhs.m_index);
+		std::swap(m_head, rhs.m_head);
+		std::swap(m_tail, rhs.m_tail);
 	}

 	return *this;
--- a/src/compound.cpp
+++ b/src/compound.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/compound.hpp>
+#include "cif++.hpp"

 #include <filesystem>
 #include <fstream>
--- a/src/condition.cpp
+++ b/src/condition.cpp
@@ -24,8 +24,8 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/category.hpp>
-#include <cif++/condition.hpp>
+#include "cif++/category.hpp"
+#include "cif++/condition.hpp"

 namespace cif
 {
--- a/src/datablock.cpp
+++ b/src/datablock.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/datablock.hpp>
+#include "cif++/datablock.hpp"

 namespace cif
 {
--- a/src/dictionary_parser.cpp
+++ b/src/dictionary_parser.cpp
@@ -24,10 +24,10 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/condition.hpp>
-#include <cif++/dictionary_parser.hpp>
-#include <cif++/file.hpp>
-#include <cif++/parser.hpp>
+#include "cif++/condition.hpp"
+#include "cif++/dictionary_parser.hpp"
+#include "cif++/file.hpp"
+#include "cif++/parser.hpp"

 namespace cif
 {
--- a/src/file.cpp
+++ b/src/file.cpp
@@ -24,8 +24,8 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/file.hpp>
-#include <cif++/gzio.hpp>
+#include "cif++/file.hpp"
+#include "cif++/gzio.hpp"

 namespace cif
 {
--- a/src/item.cpp
+++ b/src/item.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/row.hpp>
+#include "cif++/row.hpp"

 #include <cassert>

--- a/src/model.cpp
+++ b/src/model.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++.hpp>
+#include "cif++.hpp"

 #include <filesystem>
 #include <fstream>
@@ -290,12 +290,17 @@ int atom::atom_impl::get_charge() const

 std::ostream &operator<<(std::ostream &os, const atom &atom)
 {
-	os << atom.get_label_comp_id() << ' ' << atom.get_label_asym_id() << ':' << atom.get_label_seq_id() << ' ' << atom.get_label_atom_id();
+	if (atom.is_water())
+		os << atom.get_label_comp_id() << ' ' << atom.get_label_asym_id() << ':' << atom.get_auth_seq_id() << ' ' << atom.get_label_atom_id();
+	else
+	{
+		os << atom.get_label_comp_id() << ' ' << atom.get_label_asym_id() << ':' << atom.get_label_seq_id() << ' ' << atom.get_label_atom_id();

-	if (atom.is_alternate())
-		os << '(' << atom.get_label_alt_id() << ')';
-	if (atom.get_auth_asym_id() != atom.get_label_asym_id() or atom.get_auth_seq_id() != std::to_string(atom.get_label_seq_id()) or atom.get_pdb_ins_code().empty() == false)
-		os << " [" << atom.get_auth_asym_id() << ':' << atom.get_auth_seq_id() << atom.get_pdb_ins_code() << ']';
+		if (atom.is_alternate())
+			os << '(' << atom.get_label_alt_id() << ')';
+		if (atom.get_auth_asym_id() != atom.get_label_asym_id() or atom.get_auth_seq_id() != std::to_string(atom.get_label_seq_id()) or atom.get_pdb_ins_code().empty() == false)
+			os << " [" << atom.get_auth_asym_id() << ':' << atom.get_auth_seq_id() << atom.get_pdb_ins_code() << ']';
+	}

 	return os;
 }
@@ -807,11 +812,17 @@ float monomer::chi(size_t nr) const
 					atoms.back() = "CG2";
 			}

-			result = static_cast<float>(dihedral_angle(
-				get_atom_by_atom_id(atoms[nr + 0]).get_location(),
-				get_atom_by_atom_id(atoms[nr + 1]).get_location(),
-				get_atom_by_atom_id(atoms[nr + 2]).get_location(),
-				get_atom_by_atom_id(atoms[nr + 3]).get_location()));
+			auto atom_0 = get_atom_by_atom_id(atoms[nr + 0]);
+			auto atom_1 = get_atom_by_atom_id(atoms[nr + 1]);
+			auto atom_2 = get_atom_by_atom_id(atoms[nr + 2]);
+			auto atom_3 = get_atom_by_atom_id(atoms[nr + 3]);
+
+			if (atom_0 and atom_1 and atom_2 and atom_3)
+				result = static_cast<float>(dihedral_angle(
+					atom_0.get_location(),
+					atom_1.get_location(),
+					atom_2.get_location(),
+					atom_3.get_location()));
 		}
 	}
 	catch (const std::exception &e)
@@ -1365,23 +1376,16 @@ structure::structure(datablock &db, size_t modelNr, StructureOpenOptions options

 void structure::load_atoms_for_model(StructureOpenOptions options)
 {
+	using namespace literals;
+
 	auto &atomCat = m_db["atom_site"];

-	for (const auto &a : atomCat)
-	{
-		std::string id, type_symbol;
-		std::optional<size_t> model_nr;
-
-		cif::tie(id, type_symbol, model_nr) = a.get("id", "type_symbol", "pdbx_PDB_model_num");
-
-		if (model_nr and *model_nr != m_model_nr)
-			continue;
-
-		if ((options bitand StructureOpenOptions::SkipHydrogen) and (type_symbol == "H" or type_symbol == "D"))
-			continue;
+	condition c = "pdbx_PDB_model_num"_key == null or "pdbx_PDB_model_num"_key == m_model_nr;
+	if (options bitand StructureOpenOptions::SkipHydrogen)
+		c = std::move(c) and ("type_symbol"_key != "H" and "type_symbol"_key != "D");

+	for (auto id : atomCat.find<std::string>(std::move(c), "id"))
 		emplace_atom(std::make_shared<atom::atom_impl>(m_db, id));
-	}
 }

 // structure::structure(const structure &s)
@@ -1539,6 +1543,36 @@ EntityType structure::get_entity_type_for_asym_id(const std::string asym_id) con
 // 	return result;
 // }

+bool structure::has_atom_id(const std::string &id) const
+{
+	assert(m_atoms.size() == m_atom_index.size());
+
+	bool result = false;
+
+	int L = 0, R = static_cast<int>(m_atoms.size() - 1);
+	while (L <= R)
+	{
+		int i = (L + R) / 2;
+
+		const atom &atom = m_atoms[m_atom_index[i]];
+
+		int d = atom.id().compare(id);
+
+		if (d == 0)
+		{
+			result = true;
+			break;
+		}
+
+		if (d < 0)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	return result;
+}
+
 atom structure::get_atom_by_id(const std::string &id) const
 {
 	assert(m_atoms.size() == m_atom_index.size());
@@ -1845,9 +1879,19 @@ void structure::remove_atom(atom &a, bool removeFromResidue)
 {
 	using namespace literals;

-	auto &atomSites = m_db["atom_site"];
+	auto &atomSite = m_db["atom_site"];

-	if (removeFromResidue)
+	if (a.is_water())
+	{
+		auto ra = atomSite.find1("id"_key == a.id());
+		if (ra)
+		{
+			auto &nps = m_db["pdbx_nonpoly_scheme"];
+			for (auto rnp : atomSite.get_children(ra, nps))
+				nps.erase(rnp);
+		}
+	}
+	else if (removeFromResidue)
 	{
 		try
 		{
@@ -1861,7 +1905,37 @@ void structure::remove_atom(atom &a, bool removeFromResidue)
 		}
 	}

-	atomSites.erase("id"_key == a.id());
+	for (auto ri : atomSite.find("id"_key == a.id()))
+	{
+		// also remove struct_conn records for this atom
+		auto &structConn = m_db["struct_conn"];
+
+		condition cond;
+
+		for (std::string prefix : { "ptnr1_", "ptnr2_", "pdbx_ptnr3_" })
+		{
+			if (a.get_label_seq_id() == 0)
+				cond = std::move(cond) or (
+					cif::key(prefix + "label_asym_id") == a.get_label_asym_id() and
+					cif::key(prefix + "label_seq_id") == null and
+					cif::key(prefix + "auth_seq_id") == a.get_auth_seq_id() and
+					cif::key(prefix + "label_atom_id") == a.get_label_atom_id()
+				);
+			else
+				cond = std::move(cond) or (
+					cif::key(prefix + "label_asym_id") == a.get_label_asym_id() and
+					cif::key(prefix + "label_seq_id") == a.get_label_seq_id() and
+					cif::key(prefix + "auth_seq_id") == a.get_auth_seq_id() and
+					cif::key(prefix + "label_atom_id") == a.get_label_atom_id()
+				);
+		}
+
+		if (cond)
+			structConn.erase(std::move(cond));
+
+		atomSite.erase(ri);
+		break;
+	}

 	assert(m_atom_index.size() == m_atoms.size());

@@ -2224,6 +2298,7 @@ void structure::remove_branch(branch &branch)

 	m_db["pdbx_branch_scheme"].erase("asym_id"_key == branch.get_asym_id());
 	m_db["struct_asym"].erase("id"_key == branch.get_asym_id());
+	m_db["struct_conn"].erase("ptnr1_label_asym_id"_key == branch.get_asym_id() or "ptnr2_label_asym_id"_key == branch.get_asym_id());

 	m_branches.erase(remove(m_branches.begin(), m_branches.end(), branch), m_branches.end());
 }
@@ -2329,18 +2404,19 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 	{
 		auto atom_id = atom_site.get_unique_id("");

-		atom.set_value("name", atom_id);
+		atom.set_value("id", atom_id);
 		atom.set_value("label_asym_id", asym_id);
 		atom.set_value("auth_asym_id", asym_id);
 		atom.set_value("label_entity_id", entity_id);

 		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", ""});
+		atom.set_value_if_empty({"label_seq_id", "."});
 		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({"occupancy", 1.0, 2});

 		auto row = atom_site.emplace(atom.begin(), atom.end());

@@ -2366,6 +2442,72 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 	return asym_id;
 }

+void structure::create_water(row_initializer atom)
+{
+	using namespace literals;
+
+	auto entity_id = insert_compound("HOH", true);
+
+	auto &struct_asym = m_db["struct_asym"];
+	std::string asym_id;
+	try
+	{
+		asym_id = struct_asym.find1<std::string>("entity_id"_key == entity_id, "id");
+	}
+	catch (const std::exception &)
+	{
+		asym_id = struct_asym.get_unique_id();
+
+		struct_asym.emplace({
+			{"id", asym_id},
+			{"pdbx_blank_PDB_chainid_flag", "N"},
+			{"pdbx_modified", "N"},
+			{"entity_id", entity_id},
+			{"details", "?"}
+		});
+	}
+	
+	auto &atom_site = m_db["atom_site"];
+	auto auth_seq_id = atom_site.find_max<int>("auth_seq_id", "label_entity_id"_key == entity_id) + 1;
+	if (auth_seq_id < 0)
+		auth_seq_id = 1;
+
+	auto atom_id = atom_site.get_unique_id("");
+
+	atom.set_value("id", atom_id);
+	atom.set_value("label_asym_id", asym_id);
+	atom.set_value("auth_asym_id", asym_id);
+	atom.set_value("label_entity_id", entity_id);
+	atom.set_value("auth_seq_id", std::to_string(auth_seq_id));
+
+	atom.set_value_if_empty({"group_PDB", "HETATM"});
+	atom.set_value_if_empty({"label_comp_id", "HOH"});
+	atom.set_value_if_empty({"label_seq_id", "."});
+	atom.set_value_if_empty({"auth_comp_id", "HOH"});
+	atom.set_value_if_empty({"pdbx_PDB_model_num", 1});
+	atom.set_value_if_empty({"label_alt_id", ""});
+	atom.set_value_if_empty({"occupancy", 1.0, 2});
+
+	auto row = atom_site.emplace(atom.begin(), atom.end());
+
+	emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id));
+
+	auto &pdbx_nonpoly_scheme = m_db["pdbx_nonpoly_scheme"];
+	int ndb_nr = pdbx_nonpoly_scheme.find_max<int>("ndb_seq_num") + 1;
+	pdbx_nonpoly_scheme.emplace({
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"mon_id", "HOH"},
+		{"ndb_seq_num", ndb_nr},
+		{"pdb_seq_num", auth_seq_id},
+		{"auth_seq_num", auth_seq_id},
+		{"pdb_mon_id", "HOH"},
+		{"auth_mon_id", "HOH"},
+		{"pdb_strand_id", asym_id},
+		{"pdb_ins_code", "."},
+	});
+}
+
 branch &structure::create_branch()
 {
 	auto &entity = m_db["entity"];
@@ -2605,10 +2747,10 @@ std::string structure::create_entity_for_branch(branch &branch)
 		{
 			auto l2 = s1.get_link();

-			if (not l2)
+			if (not l2 or l2.get_auth_seq_id().empty())
 				continue;

-			auto &s2 = branch.at(std::stoi(l2.get_auth_seq_id()) - 1);
+			auto &s2 = branch.at(stoi(l2.get_auth_seq_id()) - 1);
 			auto l1 = s2.get_atom_by_atom_id("C1");

 			pdbx_entity_branch_link.emplace({
--- a/src/parser.cpp
+++ b/src/parser.cpp
@@ -24,10 +24,10 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/utilities.hpp>
-#include <cif++/forward_decl.hpp>
-#include <cif++/parser.hpp>
-#include <cif++/file.hpp>
+#include "cif++/utilities.hpp"
+#include "cif++/forward_decl.hpp"
+#include "cif++/parser.hpp"
+#include "cif++/file.hpp"

 #include <cassert>
 #include <iostream>
--- a/src/pdb/cif2pdb.cpp
+++ b/src/pdb/cif2pdb.cpp
@@ -24,9 +24,9 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++.hpp>
-#include <cif++/pdb/cif2pdb.hpp>
-#include <cif++/gzio.hpp>
+#include "cif++.hpp"
+#include "cif++/pdb/cif2pdb.hpp"
+#include "cif++/gzio.hpp"

 #include <cmath>
 #include <deque>
@@ -1406,10 +1406,24 @@ void WriteRemark3Refmac(std::ostream &pdbFile, const datablock &db)
 		}
 	}

-	// TODO: add twin information
+	pdbFile << RM3("") << std::endl
+			<< RM3(" TWIN DETAILS") << std::endl;

-	//	{ R"(TWIN DETAILS)", "", {} },
-	//	{ R"(NUMBER OF TWIN DOMAINS)", "", {} },
+	auto &twins = db["pdbx_reflns_twin"];
+	if (twins.empty())
+		pdbFile << RM3("  NUMBER OF TWIN DOMAINS  : NULL") << std::endl;
+	else
+	{
+		pdbFile << RM3("  NUMBER OF TWIN DOMAINS  :    ") << twins.size() << std::endl;
+
+		int nr = 1;
+		for (auto twin : twins)
+		{
+			pdbFile << RM3("     TWIN DOMAIN   : ") << nr++ << std::endl
+					<< RM3("     TWIN OPERATOR : ") << Fs(twin, "operator") << std::endl
+					<< RM3("     TWIN FRACTION : ") << SEP("", -6, 3) << Ff(twin, "fraction") << std::endl;
+		}
+	}

 	auto &tls = db["pdbx_refine_tls"];

@@ -3358,8 +3372,8 @@ std::tuple<int, int> WriteCoordinatesForModel(std::ostream &pdbFile, const datab
 	auto &atom_site = db["atom_site"];
 	auto &atom_site_anisotrop = db["atom_site_anisotrop"];
 	auto &entity = db["entity"];
-	auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
-	auto &pdbx_nonpoly_scheme = db["pdbx_nonpoly_scheme"];
+	// auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
+	// auto &pdbx_nonpoly_scheme = db["pdbx_nonpoly_scheme"];
 	auto &pdbx_branch_scheme = db["pdbx_branch_scheme"];

 	int serial = 1;
@@ -3434,16 +3448,26 @@ std::tuple<int, int> WriteCoordinatesForModel(std::ostream &pdbFile, const datab
 			r.get("id", "group_PDB", "label_atom_id", "label_alt_id", "auth_comp_id", "auth_asym_id", "auth_seq_id",
 				"pdbx_PDB_ins_code", "Cartn_x", "Cartn_y", "Cartn_z", "occupancy", "B_iso_or_equiv", "type_symbol", "pdbx_formal_charge");

-		int entity_id = r.get<int>("label_entity_id");
-		auto type = entity.find1<std::string>("id"_key == entity_id, "type");
-
-		if (type == "branched")	// find the real auth_seq_num, since sugars have their auth_seq_num reused as sugar number... sigh.
-			resSeq = pdbx_branch_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
-		// else if (type == "non-polymer")	// same for non-polymers
-		// 	resSeq = pdbx_nonpoly_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
-		else if (type == "polymer")
-			resSeq = pdbx_poly_seq_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
+		if (resName != "HOH")
+		{
+			int entity_id = r.get<int>("label_entity_id");
+			try
+			{
+				auto type = entity.find1<std::string>("id"_key == entity_id, "type");

+				if (type == "branched")	// find the real auth_seq_num, since sugars have their auth_seq_num reused as sugar number... sigh.
+					resSeq = pdbx_branch_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
+				// else if (type == "non-polymer")	// same for non-polymers
+				// 	resSeq = pdbx_nonpoly_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
+				// else if (type == "polymer")
+				// 	resSeq = pdbx_poly_seq_scheme.find1<int>("asym_id"_key == r.get<std::string>("label_asym_id") and "pdb_seq_num"_key == resSeq, "auth_seq_num");
+			}
+			catch (const std::exception &ex)
+			{
+				std::cerr << "Oops, there was not exactly one entity with id " << entity_id << std::endl;
+			}
+		}
+		
 		if (chainID.length() > 1)
 			throw std::runtime_error("Chain ID " + chainID + " won't fit into a PDB file");

--- a/src/pdb/pdb2cif_remark_3.cpp
+++ b/src/pdb/pdb2cif_remark_3.cpp
@@ -24,8 +24,8 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++.hpp>
-#include <cif++/pdb/pdb2cif_remark_3.hpp>
+#include "cif++.hpp"
+#include "cif++/pdb/pdb2cif_remark_3.hpp"

 #include <map>
 #include <set>
@@ -1233,7 +1233,9 @@ void Remark3Parser::storeCapture(const char *category, std::initializer_list<con
 			{
 				cat.emplace({ // #warning("crystal id, diffrn id, what should be put here?")
 					{ "crystal_id", 1 },
-					{ "diffrn_id", 1 } });
+					{ "diffrn_id", 1 },
+					{ "operator", "" },
+					{ "fraction", 0.f } });
 			}
 			else if (iequals(category, "reflns"))
 				cat.emplace({ { "pdbx_ordinal", cat.size() + 1 },
--- a/src/pdb/tls.cpp
+++ b/src/pdb/tls.cpp
@@ -27,8 +27,8 @@
 // #include <sys/ioctl.h>
 // #include <termios.h>

-#include <cif++.hpp>
-#include <cif++/pdb/tls.hpp>
+#include "cif++.hpp"
+#include "cif++/pdb/tls.hpp"

 #include <iomanip>
 #include <iostream>
--- a/src/point.cpp
+++ b/src/point.cpp
@@ -24,7 +24,8 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/point.hpp>
+#include "cif++/point.hpp"
+#include "cif++/matrix.hpp"

 #include <cassert>
 #include <random>
@@ -32,245 +33,6 @@
 namespace cif
 {

-// --------------------------------------------------------------------
-// We're using expression templates here
-
-template <typename M>
-class MatrixExpression
-{
-  public:
-	uint32_t dim_m() const { return static_cast<const M &>(*this).dim_m(); }
-	uint32_t dim_n() const { return static_cast<const M &>(*this).dim_n(); }
-
-	double &operator()(uint32_t i, uint32_t j)
-	{
-		return static_cast<M &>(*this).operator()(i, j);
-	}
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return static_cast<const M &>(*this).operator()(i, j);
-	}
-};
-
-// --------------------------------------------------------------------
-// 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
-
-class Matrix : public MatrixExpression<Matrix>
-{
-  public:
-	template <typename M2>
-	Matrix(const MatrixExpression<M2> &m)
-		: m_m(m.dim_m())
-		, m_n(m.dim_n())
-		, m_data(m_m * m_n)
-	{
-		for (uint32_t i = 0; i < m_m; ++i)
-		{
-			for (uint32_t j = 0; j < m_n; ++j)
-				operator()(i, j) = m(i, j);
-		}
-	}
-
-	Matrix(size_t m, size_t n, double v = 0)
-		: m_m(m)
-		, m_n(n)
-		, m_data(m_m * m_n)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	Matrix() = default;
-	Matrix(Matrix &&m) = default;
-	Matrix(const Matrix &m) = default;
-	Matrix &operator=(Matrix &&m) = default;
-	Matrix &operator=(const Matrix &m) = default;
-
-	size_t dim_m() const { return m_m; }
-	size_t dim_n() const { return m_n; }
-
-	double operator()(size_t i, size_t j) const
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-	double &operator()(size_t i, size_t j)
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-  private:
-	size_t m_m = 0, m_n = 0;
-	std::vector<double> m_data;
-};
-
-// --------------------------------------------------------------------
-
-class SymmetricMatrix : public MatrixExpression<SymmetricMatrix>
-{
-  public:
-	SymmetricMatrix(uint32_t n, double v = 0)
-		: m_n(n)
-		, m_data((m_n * (m_n + 1)) / 2)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	SymmetricMatrix() = default;
-	SymmetricMatrix(SymmetricMatrix &&m) = default;
-	SymmetricMatrix(const SymmetricMatrix &m) = default;
-	SymmetricMatrix &operator=(SymmetricMatrix &&m) = default;
-	SymmetricMatrix &operator=(const SymmetricMatrix &m) = default;
-
-	uint32_t dim_m() const { return m_n; }
-	uint32_t dim_n() const { return m_n; }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return i < j
-		           ? m_data[(j * (j + 1)) / 2 + i]
-		           : m_data[(i * (i + 1)) / 2 + j];
-	}
-
-	double &operator()(uint32_t i, uint32_t j)
-	{
-		if (i > j)
-			std::swap(i, j);
-		assert(j < m_n);
-		return m_data[(j * (j + 1)) / 2 + i];
-	}
-
-  private:
-	uint32_t m_n;
-	std::vector<double> m_data;
-};
-
-class IdentityMatrix : public MatrixExpression<IdentityMatrix>
-{
-  public:
-	IdentityMatrix(uint32_t n)
-		: m_n(n)
-	{
-	}
-
-	uint32_t dim_m() const { return m_n; }
-	uint32_t dim_n() const { return m_n; }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return i == j ? 1 : 0;
-	}
-
-  private:
-	uint32_t m_n;
-};
-
-// --------------------------------------------------------------------
-// matrix functions, implemented as expression templates
-
-template <typename M1, typename M2>
-class MatrixSubtraction : public MatrixExpression<MatrixSubtraction<M1, M2>>
-{
-  public:
-	MatrixSubtraction(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());
-	}
-
-	uint32_t dim_m() const { return m_m1.dim_m(); }
-	uint32_t dim_n() const { return m_m1.dim_n(); }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return m_m1(i, j) - m_m2(i, j);
-	}
-
-  private:
-	const M1 &m_m1;
-	const M2 &m_m2;
-};
-
-template <typename M1, typename M2>
-MatrixSubtraction<M1, M2> operator-(const MatrixExpression<M1> &m1, const MatrixExpression<M2> &m2)
-{
-	return MatrixSubtraction(*static_cast<const M1 *>(&m1), *static_cast<const M2 *>(&m2));
-}
-
-template <typename M>
-class MatrixMultiplication : public MatrixExpression<MatrixMultiplication<M>>
-{
-  public:
-	MatrixMultiplication(const M &m, double v)
-		: m_m(m)
-		, m_v(v)
-	{
-	}
-
-	uint32_t dim_m() const { return m_m.dim_m(); }
-	uint32_t dim_n() const { return m_m.dim_n(); }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return m_m(i, j) * m_v;
-	}
-
-  private:
-	const M &m_m;
-	double m_v;
-};
-
-template <typename M>
-MatrixMultiplication<M> operator*(const MatrixExpression<M> &m, double v)
-{
-	return MatrixMultiplication(*static_cast<const M *>(&m), v);
-}
-
-// --------------------------------------------------------------------
-
-template <class M1>
-Matrix Cofactors(const M1 &m)
-{
-	Matrix cf(m.dim_m(), m.dim_m());
-
-	const size_t ixs[4][3] = {
-		{ 1, 2, 3 },
-		{ 0, 2, 3 },
-		{ 0, 1, 3 },
-		{ 0, 1, 2 }
-	};
-
-	for (size_t x = 0; x < 4; ++x)
-	{
-		const size_t *ix = ixs[x];
-
-		for (size_t y = 0; y < 4; ++y)
-		{
-			const size_t *iy = ixs[y];
-
-			cf(x, y) =
-				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]);
-
-			if ((x + y) % 2 == 1)
-				cf(x, y) *= -1;
-		}
-	}
-
-	return cf;
-}
-
 // --------------------------------------------------------------------

 template<typename T>
@@ -299,13 +61,14 @@ quaternion_type<T> normalize(quaternion_type<T> q)

 quaternion construct_from_angle_axis(float angle, point axis)
 {
-	auto q = std::cos((angle * kPI / 180) / 2);
-	auto s = std::sqrt(1 - q * q);
+	angle = (angle * kPI / 180) / 2;
+	auto s = std::sin(angle);
+	auto c = std::cos(angle);

 	axis.normalize();

 	return normalize(quaternion{
-		static_cast<float>(q),
+		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) });
@@ -364,33 +127,10 @@ quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4,
 	p3 -= p3;

 	quaternion q;
-	auto axis = p2;
+	auto axis = -p2;

 	float dh = dihedral_angle(p1, p2, p3, p4);
-	for (int iteration = 0; iteration < 100; ++iteration)
-	{
-		float delta = std::fmod(angle - dh, 360.0f);
-
-		if (delta < -180)
-			delta += 360;
-		if (delta > 180)
-			delta -= 360;
-
-		if (std::abs(delta) < esd)
-			break;
-
-		// if (iteration > 0)
-		// 	std::cout << cif::coloured(("iteration " + std::to_string(iteration)).c_str(), cif::scBLUE, cif::scBLACK) << " delta: " << delta << std::endl;
-
-		auto q2 = construct_from_angle_axis(delta, axis);
-		q = iteration == 0 ? q2 : q * q2;
-			
-		p4.rotate(q2);
-
-		dh = dihedral_angle(p1, p2, p3, p4);
-	}
-
-	return q;
+	return construct_from_angle_axis(angle - dh, axis);
 }

 point centroid(const std::vector<point> &pts)
@@ -465,8 +205,8 @@ double LargestDepressedQuarticSolution(double a, double b, double c)

 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
-	Matrix M(3, 3, 0);
+	// First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
+	matrix3x3<double> M;

 	for (uint32_t i = 0; i < pa.size(); ++i)
 	{
@@ -484,8 +224,8 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 		M(2, 2) += a.m_z * b.m_z;
 	}

-	// Now calculate N, a symmetric 4x4 Matrix
-	SymmetricMatrix N(4);
+	// Now calculate N, a symmetric 4x4 matrix
+	symmetric_matrix4x4<double> N(4);

 	N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
 	N(0, 1) = M(1, 2) - M(2, 1);
@@ -534,16 +274,22 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	double lambda = LargestDepressedQuarticSolution(C, D, E);

 	// calculate t = (N - λI)
-	Matrix t = N - IdentityMatrix(4) * lambda;
+	matrix<double> t(N - identity_matrix(4) * lambda);

-	// calculate a Matrix of cofactors for t
-	Matrix cf = Cofactors(t);
+	// calculate a matrix of cofactors for t
+	auto cf = matrix_cofactors(t);

 	int maxR = 0;
+	double maxCF = std::abs(cf(0, 0));
+
 	for (int r = 1; r < 4; ++r)
 	{
-		if (std::abs(cf(r, 0)) > std::abs(cf(maxR, 0)))
+		auto cfr = std::abs(cf(r, 0));
+		if (maxCF < cfr)
+		{
+			maxCF = cfr;
 			maxR = r;
+		}
 	}

 	quaternion q(
--- a/src/row.cpp
+++ b/src/row.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/category.hpp>
+#include "cif++/category.hpp"

 namespace cif
 {
--- a/src/symmetry.cpp
+++ b/src/symmetry.cpp
@@ -1,17 +1,17 @@
 /*-
 * SPDX-License-Identifier: BSD-2-Clause
- * 
+ *
 * Copyright (c) 2020 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
@@ -24,23 +24,301 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/symmetry.hpp>
+#include "cif++/symmetry.hpp"
+#include "cif++/datablock.hpp"
+#include "cif++/point.hpp"

 #include <stdexcept>

-#include "./symop_table_data.hpp"
+#include "symop_table_data.hpp"
+
+#include <Eigen/Eigenvalues>

 namespace cif
 {

 // --------------------------------------------------------------------
-// Unfortunately, clipper has a different numbering scheme than PDB
-// for rotation numbers. So we created a table to map those.
-// Perhaps a bit over the top, but hey....
+
+cell::cell(float a, float b, float c, float alpha, float beta, float gamma)
+	: m_a(a)
+	, m_b(b)
+	, m_c(c)
+	, m_alpha(alpha)
+	, m_beta(beta)
+	, m_gamma(gamma)
+{
+	init();
+}
+
+cell::cell(const datablock &db)
+{
+	auto &_cell = db["cell"];
+
+	tie(m_a, m_b, m_c, m_alpha, m_beta, m_gamma) =
+		_cell.front().get("length_a", "length_b", "length_c", "angle_alpha", "angle_beta", "angle_gamma");
+
+	init();
+}
+
+void cell::init()
+{
+	auto alpha = (m_alpha * kPI) / 180;
+	auto beta = (m_beta * kPI) / 180;
+	auto gamma = (m_gamma * kPI) / 180;
+
+	auto alpha_star = std::acos((std::cos(gamma) * std::cos(beta) - std::cos(alpha)) / (std::sin(beta) * std::sin(gamma)));
+
+	m_orthogonal = identity_matrix(3);
+
+	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);
+}

 // --------------------------------------------------------------------

-int get_space_group_number(std::string spacegroup)
+sym_op::sym_op(std::string_view s)
+{
+	auto b = s.data();
+	auto e = b + s.length();
+
+	int rnri = 256;	// default to unexisting number
+	auto r = std::from_chars(b, e, rnri);
+	
+	m_nr = rnri;
+	m_ta = r.ptr[1] - '0';
+	m_tb = r.ptr[2] - '0';
+	m_tc = r.ptr[3] - '0';
+
+	if (r.ec != std::errc() or rnri > 192 or r.ptr[0] != '_' or m_ta > 9 or m_tb > 9 or m_tc > 9)
+		throw std::invalid_argument("Could not convert string into sym_op");
+}
+
+std::string sym_op::string() const
+{
+	char b[9];
+	auto r = std::to_chars(b, b + sizeof(b), m_nr);
+	if (r.ec != std::errc() or r.ptr > b + 4)
+		throw std::runtime_error("Could not write out symmetry operation to string");
+	
+	*r.ptr++ = '_';
+	*r.ptr++ = '0' + m_ta;
+	*r.ptr++ = '0' + m_tb;
+	*r.ptr++ = '0' + m_tc;
+	*r.ptr = 0;
+
+	return { b, static_cast<size_t>(r.ptr - b) };
+}
+
+// --------------------------------------------------------------------
+
+transformation::transformation(const symop_data &data)
+{
+	const auto &d = data.data();
+
+	m_rotation(0, 0) = d[0];
+	m_rotation(0, 1) = d[1];
+	m_rotation(0, 2) = d[2];
+	m_rotation(1, 0) = d[3];
+	m_rotation(1, 1) = d[4];
+	m_rotation(1, 2) = d[5];
+	m_rotation(2, 0) = d[6];
+	m_rotation(2, 1) = d[7];
+	m_rotation(2, 2) = d[8];
+
+	try_create_quaternion();
+
+	m_translation.m_x = d[9] == 0 ? 0 : 1.0 * d[9] / d[10];
+	m_translation.m_y = d[11] == 0 ? 0 : 1.0 * d[11] / d[12];
+	m_translation.m_z = d[13] == 0 ? 0 : 1.0 * d[13] / d[14];
+}
+
+transformation::transformation(const matrix3x3<float> &r, const cif::point &t)
+	: m_rotation(r)
+	, m_translation(t)
+{
+	try_create_quaternion();
+}
+
+void transformation::try_create_quaternion()
+{
+	float Qxx = m_rotation(0, 0);
+	float Qxy = m_rotation(0, 1);
+	float Qxz = m_rotation(0, 2);
+	float Qyx = m_rotation(1, 0);
+	float Qyy = m_rotation(1, 1);
+	float Qyz = m_rotation(1, 2);
+	float Qzx = m_rotation(2, 0);
+	float Qzy = m_rotation(2, 1);
+	float Qzz = m_rotation(2, 2);
+
+	Eigen::Matrix4f em;
+
+	em << 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;
+
+	Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
+
+	auto ev = es.eigenvalues();
+
+	for (size_t j = 0; j < 4; ++j)
+	{
+		if (std::abs(ev[j].real() - 1) > 0.01)
+			continue;
+		
+		auto col = es.eigenvectors().col(j);
+
+		m_q = normalize(cif::quaternion{
+			static_cast<float>(col(3).real()),
+			static_cast<float>(col(0).real()),
+			static_cast<float>(col(1).real()),
+			static_cast<float>(col(2).real()) });
+
+		break;
+	}
+}
+
+transformation operator*(const transformation &lhs, const transformation &rhs)
+{
+	auto r = lhs.m_rotation * rhs.m_rotation;
+	auto t = lhs.m_rotation * rhs.m_translation;
+	t = t + lhs.m_translation;
+
+	return transformation(r, t);
+}
+
+transformation inverse(const transformation &t)
+{
+	auto inv_matrix = inverse(t.m_rotation);
+	return { inv_matrix, -(inv_matrix * t.m_translation) };
+}
+
+// --------------------------------------------------------------------
+
+spacegroup::spacegroup(int nr)
+	: m_nr(nr)
+{
+	const size_t N = kSymopNrTableSize;
+	int32_t L = 0, R = static_cast<int32_t>(N - 1);
+	while (L <= R)
+	{
+		int32_t i = (L + R) / 2;
+		if (kSymopNrTable[i].spacegroup() < m_nr)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	m_index = L;
+
+	for (size_t i = L; i < N and kSymopNrTable[i].spacegroup() == m_nr; ++i)
+		emplace_back(kSymopNrTable[i].symop().data());
+}
+
+std::string spacegroup::get_name() const
+{
+	for (auto &s : kSpaceGroups)
+	{
+		if (s.nr == m_nr)
+			return s.name;
+	}
+
+	throw std::runtime_error("Spacegroup has an invalid number: " + std::to_string(m_nr));
+}
+
+point offsetToOrigin(const cell &c, const point &p)
+{
+	point d{};
+
+	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.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.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;
+};
+
+point offsetToOriginFractional(const point &p)
+{
+	point d{};
+
+	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.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.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;
+};
+
+point spacegroup::operator()(const point &pt, const cell &c, sym_op symop) const
+{
+	if (symop.m_nr < 1 or symop.m_nr > size())
+		throw std::out_of_range("symmetry operator number out of range");
+	
+	transformation t = at(symop.m_nr - 1);
+
+	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 fpt = fractional(pt, c);
+	auto o = offsetToOriginFractional(fpt);
+
+	auto spt = t(fpt + o) - o;
+
+	return orthogonal(spt, c);
+}
+
+point spacegroup::inverse(const point &pt, const cell &c, sym_op symop) const
+{
+	if (symop.m_nr < 1 or symop.m_nr > size())
+		throw std::out_of_range("symmetry operator number out of range");
+	
+	transformation t = at(symop.m_nr - 1);
+
+	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 fpt = fractional(pt, c);
+	auto o = offsetToOriginFractional(fpt);
+
+	auto it = cif::inverse(t);
+	auto spt = it(fpt + o) - o;
+
+	return orthogonal(spt, c);
+}
+
+// --------------------------------------------------------------------
+
+int get_space_group_number(std::string_view spacegroup)
 {
 	if (spacegroup == "P 21 21 2 A")
 		spacegroup = "P 21 21 2 (a)";
@@ -73,7 +351,7 @@ int get_space_group_number(std::string spacegroup)
 	{
 		for (size_t i = 0; i < kNrOfSpaceGroups; ++i)
 		{
-			auto& sp = kSpaceGroups[i];
+			auto &sp = kSpaceGroups[i];
 			if (sp.xHM == spacegroup)
 			{
 				result = sp.nr;
@@ -83,14 +361,14 @@ int get_space_group_number(std::string spacegroup)
 	}

 	if (result == 0)
-		throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
-	
+		throw std::runtime_error("Spacegroup name " + std::string(spacegroup) + " was not found in table");
+
 	return result;
 }

 // --------------------------------------------------------------------

-int get_space_group_number(std::string spacegroup, space_group_name type)
+int get_space_group_number(std::string_view spacegroup, space_group_name type)
 {
 	if (spacegroup == "P 21 21 2 A")
 		spacegroup = "P 21 21 2 (a)";
@@ -145,9 +423,99 @@ int get_space_group_number(std::string spacegroup, space_group_name type)

 	// not found, see if we can find a match based on xHM name
 	if (result == 0)
-		throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
-	
+		throw std::runtime_error("Spacegroup name " + std::string(spacegroup) + " was not found in table");
+
 	return result;
 }

+int get_space_group_number(const datablock &db)
+{
+	auto &_symmetry = db["symmetry"];
+
+	if (_symmetry.size() != 1)
+		throw std::runtime_error("Could not find a unique symmetry in this mmCIF file");
+	
+	return _symmetry.front().get<int>("Int_Tables_number");
 }
+
+// --------------------------------------------------------------------
+
+std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b) const
+{
+	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;
+	float result_d = std::numeric_limits<float>::max();
+	sym_op result_s;
+
+	auto fa = fractional(a, m_cell);
+	auto fb = fractional(b, m_cell);
+
+	auto o = offsetToOriginFractional(fa);
+
+	fa = fa + o;
+	fb = fb + o;
+
+	a = orthogonal(fa, m_cell);
+
+	for (size_t i = 0; i < m_spacegroup.size(); ++i)
+	{
+		sym_op s(i + 1);
+		auto &t = m_spacegroup[i];
+
+		auto fsb = t(fb);
+
+		while (fsb.m_x - 0.5f > fa.m_x)
+		{
+			fsb.m_x -= 1;
+			s.m_ta -= 1;
+		}
+
+		while (fsb.m_x + 0.5f < fa.m_x)
+		{
+			fsb.m_x += 1;
+			s.m_ta += 1;			
+		}
+
+		while (fsb.m_y - 0.5f > fa.m_y)
+		{
+			fsb.m_y -= 1;
+			s.m_tb -= 1;
+		}
+
+		while (fsb.m_y + 0.5f < fa.m_y)
+		{
+			fsb.m_y += 1;
+			s.m_tb += 1;			
+		}
+
+		while (fsb.m_z - 0.5f > fa.m_z)
+		{
+			fsb.m_z -= 1;
+			s.m_tc -= 1;
+		}
+
+		while (fsb.m_z + 0.5f < fa.m_z)
+		{
+			fsb.m_z += 1;
+			s.m_tc += 1;			
+		}
+
+		auto p = orthogonal(fsb, m_cell);
+		auto dsq = distance_squared(a, p);
+
+		if (result_d > dsq)
+		{
+			result_d = dsq;
+			result_fsb = fsb;
+			result_s = s;
+		}
+	}
+
+	auto p = orthogonal(result_fsb - o, m_cell);
+
+	return { std::sqrt(result_d), p, result_s };
+}
+
+} // namespace cif
--- a/src/symop-map-generator.cpp
+++ b/src/symop-map-generator.cpp
@@ -27,6 +27,7 @@
 #include <cassert>

 #include <array>
+#include <charconv>
 #include <iostream>
 #include <iomanip>
 #include <fstream>
@@ -169,7 +170,7 @@ class SymopParser
 	}

 	Token m_lookahead;
-	int m_nr;
+	int m_nr = -1;

 	std::string m_s;
 	std::string::const_iterator m_p, m_e;
@@ -230,14 +231,15 @@ int main(int argc, char* const argv[])

 	try
 	{
-		if (argc != 3)
+		if (argc != 4)
 		{
-			std::cerr << "Usage symop-map-generator <input-file> <output-file>" << std::endl;
+			std::cerr << "Usage symop-map-generator <syminfo.lib-file> <symop.lib-file> < <output-file>" << std::endl;
 			exit(1);
 		}

-		fs::path input(argv[1]);
-		fs::path output(argv[2]);
+		fs::path syminfolib(argv[1]);
+		fs::path symoplib(argv[2]);
+		fs::path output(argv[3]);
 		
 		tmpFile = output.parent_path() / (output.filename().string() + ".tmp");

@@ -261,22 +263,51 @@ int main(int argc, char* const argv[])
 		};

 		std::map<int,SymInfoBlock> symInfo;
-		int symopnr, mysymnr = 10000;

-		std::ifstream file(input);
+		std::ifstream file(symoplib);
+		if (not file.is_open())
+			throw std::runtime_error("Could not open symop.lib file");
+
+		std::string line;
+		int sgnr = 0;
+		int rnr = 0;
+
+		while (getline(file, line))
+		{
+			if (line.empty())
+				continue;
+			
+			if (std::isdigit(line[0]))	// start of new spacegroup
+			{
+				auto r = std::from_chars(line.data(), line.data() + line.length(), sgnr);
+				if (r.ec != std::errc())
+					throw std::runtime_error("Error parsing symop.lib file");
+				rnr = 1;
+				continue;
+			}
+
+			if (not std::isspace(line[0]) or sgnr == 0)
+				throw std::runtime_error("Error parsing symop.lib file");
+			
+			SymopParser p;
+			data.emplace_back(sgnr, rnr, p.parse(line));
+			++rnr;
+		}
+
+		file.close();
+
+		file.open(syminfolib);
 		if (not file.is_open())
 			throw std::runtime_error("Could not open syminfo.lib file");

 		enum class State { skip, spacegroup } state = State::skip;

-		std::string line;
-
 		const std::regex rx(R"(^symbol +(Hall|xHM|old) +'(.+?)'(?: +'(.+?)')?$)"),
 			rx2(R"(symbol ccp4 (\d+))");;

 		SymInfoBlock cur = {};

-		std::vector<std::array<int,15>> symops, cenops;
+		// std::vector<std::array<int,15>> symops, cenops;

 		while (getline(file, line))
 		{
@@ -286,9 +317,7 @@ int main(int argc, char* const argv[])
 					if (line == "begin_spacegroup")
 					{
 						state = State::spacegroup;
-						symopnr = 1;
-						++mysymnr;
-						cur = { mysymnr };
+						cur = {};
 					}
 					break;
 				
@@ -314,34 +343,34 @@ int main(int argc, char* const argv[])
 						if (nr != 0)
 							cur.nr = nr;
 					}
-					else if (line.compare(0, 6, "symop ") == 0)
-					{
-						SymopParser p;
-						symops.emplace_back(p.parse(line.substr(6)));
-					}
-					else if (line.compare(0, 6, "cenop ") == 0)
-					{
-						SymopParser p;
-						cenops.emplace_back(p.parse(line.substr(6)));
-					}
+					// else if (line.compare(0, 6, "symop ") == 0)
+					// {
+					// 	SymopParser p;
+					// 	symops.emplace_back(p.parse(line.substr(6)));
+					// }
+					// else if (line.compare(0, 6, "cenop ") == 0)
+					// {
+					// 	SymopParser p;
+					// 	cenops.emplace_back(p.parse(line.substr(6)));
+					// }
 					else if (line == "end_spacegroup")
 					{
-						for (auto& cenop: cenops)
-						{
-							for (auto symop: symops)
-							{
-								symop = move_symop(symop, cenop);
+					// 	for (auto& cenop: cenops)
+					// 	{
+					// 		for (auto symop: symops)
+					// 		{
+					// 			symop = move_symop(symop, cenop);

-								data.emplace_back(cur.nr, symopnr, symop);
-								++symopnr;
-							}
-						}
+					// 			data.emplace_back(cur.nr, symopnr, symop);
+					// 			++symopnr;
+					// 		}
+					// 	}

 						symInfo.emplace(cur.nr, cur);
 						state = State::skip;

-						symops.clear();
-						cenops.clear();
+					// 	symops.clear();
+					// 	cenops.clear();
 					}
 					break;
 				}
@@ -358,7 +387,7 @@ int main(int argc, char* const argv[])
 // and $CLIBD/syminfo.lib using symop-map-generator,
 // part of the PDB-REDO suite of programs.

-#include <cif++/symmetry.hpp>
+#include "cif++/symmetry.hpp"

 namespace cif
 {
@@ -383,10 +412,10 @@ const space_group kSpaceGroups[] =
 			old = '"' + old + '"' + std::string(20 - old.length(), ' ');
 			xHM = '"' + xHM + '"' + std::string(30 - xHM.length(), ' ');

-			for (std::string::size_type p = Hall.length(); p > 0; --p)
+			for (auto p = Hall.begin(); p != Hall.end(); ++p)
 			{
-				if (Hall[p - 1] == '"')
-					Hall.insert(p - 1, "\\", 1);
+				if (*p == '"')
+					p = Hall.insert(p, '\\') + 1;
 			}

 			Hall = '"' + Hall + '"' + std::string(40 - Hall.length(), ' ');
--- a/src/symop_table_data.hpp
+++ b/src/symop_table_data.hpp
--- a/src/text.cpp
+++ b/src/text.cpp
@@ -24,7 +24,7 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/text.hpp>
+#include "cif++/text.hpp"

 #include <algorithm>
 #include <cassert>
--- a/src/utilities.cpp
+++ b/src/utilities.cpp
@@ -24,12 +24,14 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/utilities.hpp>
+#include "cif++/utilities.hpp"
+
 #include "revision.hpp"

 #include <atomic>
 #include <cassert>
 #include <cmath>
+#include <condition_variable>
 #include <cstring>
 #include <deque>
 #include <fstream>
@@ -134,41 +136,53 @@ std::string get_executable_path()

 // --------------------------------------------------------------------

-struct ProgressImpl
+struct progress_bar_impl
 {
-	ProgressImpl(int64_t inMax, const std::string &inAction)
-		: mMax(inMax)
-		, mConsumed(0)
-		, mAction(inAction)
-		, mMessage(inAction)
-		, mThread(std::bind(&ProgressImpl::Run, this))
+	progress_bar_impl(int64_t inMax, const std::string &inAction)
+		: m_max_value(inMax)
+		, m_consumed(0)
+		, m_action(inAction)
+		, m_message(inAction)
+		, m_thread(std::bind(&progress_bar_impl::run, this))
 	{
 	}

-	void Run();
-	void Stop()
-	{
-		mStop = true;
-		if (mThread.joinable())
-			mThread.join();
-	}
+	progress_bar_impl(const progress_bar_impl&) = delete;
+	progress_bar_impl &operator=(const progress_bar_impl &) = delete;

-	void PrintProgress();
-	void PrintDone();
+	~progress_bar_impl();

-	int64_t mMax;
-	std::atomic<int64_t> mConsumed;
-	int64_t mLastConsumed = 0;
-	int mSpinnerIndex = 0;
-	std::string mAction, mMessage;
-	std::mutex mMutex;
-	std::thread mThread;
+	void run();
+
+	void consumed(int64_t n);
+	void progress(int64_t p);
+	void message(const std::string &msg);
+
+	void print_progress();
+	void print_done();
+
+	int64_t m_max_value;
+	std::atomic<int64_t> m_consumed;
+	int64_t m_last_consumed = 0;
+	int m_spinner_index = 0;
+	std::string m_action, m_message;
+	std::mutex m_mutex;
+	std::thread m_thread;
 	std::chrono::time_point<std::chrono::system_clock>
-		mStart = std::chrono::system_clock::now();
-	bool mStop = false;
+		m_start = std::chrono::system_clock::now();
+	bool m_stop = false;
 };

-void ProgressImpl::Run()
+progress_bar_impl::~progress_bar_impl()
+{
+	using namespace std::literals;
+	assert(m_thread.joinable());
+
+	m_stop = true;
+	m_thread.join();
+}
+
+void progress_bar_impl::run()
 {
 	using namespace std::literals;

@@ -176,21 +190,21 @@ void ProgressImpl::Run()

 	try
 	{
-		for (;;)
+		while (not m_stop)
 		{
-			std::this_thread::sleep_for(2s);
-
-			std::unique_lock lock(mMutex);
-
-			if (mStop or mConsumed == mMax)
-				break;
-
-			auto elapsed = std::chrono::system_clock::now() - mStart;
-
-			if (elapsed < std::chrono::seconds(5))
+			if (std::chrono::system_clock::now() - m_start < 2s)
+			{
+				std::this_thread::sleep_for(10ms);
 				continue;
+			}
+
+			std::lock_guard lock(m_mutex);
+
+			if (not printedAny and isatty(STDOUT_FILENO))
+				std::cout << "\e[?25l";
+
+			print_progress();

-			PrintProgress();
 			printedAny = true;
 		}
 	}
@@ -199,93 +213,93 @@ void ProgressImpl::Run()
 	}

 	if (printedAny)
-		PrintDone();
+	{
+		print_done();
+		if (isatty(STDOUT_FILENO))
+			std::cout << "\e[?25h";
+	}
 }

-void ProgressImpl::PrintProgress()
+void progress_bar_impl::consumed(int64_t n)
 {
-	//	const char* kBlocks[] = {
-	//		" ",				// 0
-	//		u8"\u258F",			// 1
-	//		u8"\u258E",			// 2
-	//		u8"\u258D",			// 3
-	//		u8"\u258C",			// 4
-	//		u8"\u258B",			// 5
-	//		u8"\u258A",			// 6
-	//		u8"\u2589",			// 7
-	//		u8"\u2588",			// 8
-	//	};
+	m_consumed += n;
+}

+void progress_bar_impl::progress(int64_t p)
+{
+	m_consumed = p;
+}
+
+void progress_bar_impl::message(const std::string &msg)
+{
+	std::unique_lock lock(m_mutex);
+	m_message = msg;
+}
+
+const char* kSpinner[] = {
+	// "▉", "▊", "▋", "▌", "▍", "▎", "▏", "▎", "▍", "▌", "▋", "▊", "▉"
+	".", "o", "O", "0", "O", "o", ".", " "
+};
+const size_t kSpinnerCount = sizeof(kSpinner) / sizeof(char*);
+const uint32_t kMinBarWidth = 40, kMinMsgWidth = 12;
+
+void progress_bar_impl::print_progress()
+{
 	const char *kBlocks[] = {
-		" ", // 0
-		" ", // 1
-		" ", // 2
-		"-", // 3
-		"-", // 4
-		"-", // 5
-		"=", // 6
-		"=", // 7
-		"=", // 8
+		// "▯", // 0
+		// "▮", // 1
+		"=",
+		"-"
 	};

 	uint32_t width = get_terminal_width();

-	std::string msg;
-	msg.reserve(width + 1);
-	if (mMessage.length() <= 20)
-	{
-		msg = mMessage;
-		if (msg.length() < 20)
-			msg.append(20 - msg.length(), ' ');
-	}
+	float progress = static_cast<float>(m_consumed) / m_max_value;
+	
+	if (width < kMinBarWidth)
+		std::cout << (100 * progress) << '%' << std::endl;
 	else
-		msg = mMessage.substr(0, 17) + "...";
-
-	msg += " |";
-
-	int64_t consumed = mConsumed;
-	float progress = static_cast<float>(consumed) / mMax;
-	int pi = static_cast<int>(std::ceil(progress * 33 * 8));
-	//	int tw = width - 28;
-	//	int twd = static_cast<int>(tw * progress + 0.5f);
-	//	msg.append(twd, '=');
-	//	msg.append(tw - twd, ' ');
-
-	for (int i = 0; i < 33; ++i)
 	{
-		if (pi <= 0)
-			msg += kBlocks[0];
-		else if (pi >= 8)
-			msg += kBlocks[8];
+		uint32_t bar_width = 7 * width / 10;
+		uint32_t pct_width = 7;
+		uint32_t msg_width = width - bar_width - pct_width - 1;
+
+		if (msg_width < kMinMsgWidth)
+		{
+			bar_width += kMinMsgWidth - msg_width;
+			msg_width = kMinMsgWidth;
+		}
+
+		std::ostringstream msg;
+
+		if (m_message.length() <= msg_width)
+		{
+			msg << m_message;
+			if (m_message.length() < msg_width)
+				msg << std::string(msg_width - m_message.length(), ' ');
+		}
 		else
-			msg += kBlocks[pi];
-		pi -= 8;
+			msg << m_message.substr(0, msg_width - 3) << "...";
+
+		msg << ' ';
+
+		uint32_t pi = static_cast<uint32_t>(std::ceil(progress * bar_width));
+
+		for (uint32_t i = 0; i < bar_width; ++i)
+			msg << kBlocks[i > pi ? 1 : 0];
+
+		msg << ' ';
+
+		msg << std::setw(3) << static_cast<int>(std::ceil(progress * 100)) << "% ";
+
+		auto now = std::chrono::system_clock::now();
+		m_spinner_index = (std::chrono::duration_cast<std::chrono::milliseconds>(now - m_start).count() / 200) % kSpinnerCount;
+
+		msg << kSpinner[m_spinner_index];
+
+		std::cout << '\r' << msg.str();
+		std::cout.flush();
 	}
-
-	msg.append("| ");
-
-	const char kSpinner[] = {' ', '.', 'o', 'O', '0', 'O', 'o', '.'};
-	const size_t kSpinnerCount = sizeof(kSpinner);
-
-	if (mLastConsumed < consumed)
-	{
-		mLastConsumed = consumed;
-		mSpinnerIndex = (mSpinnerIndex + 1) % kSpinnerCount;
-	}
-
-	const char spinner[2] = {kSpinner[mSpinnerIndex], 0};
-	msg.append(spinner);
-
-	//	int perc = static_cast<int>(100 * progress);
-	//	if (perc < 100)
-	//		msg += ' ';
-	//	if (perc < 10)
-	//		msg += ' ';
-	//	msg += to_string(perc);
-	//	msg += '%';
-
-	std::cout << '\r' << msg;
-	std::cout.flush();
 }

 namespace
@@ -324,12 +338,12 @@ namespace

 } // namespace

-void ProgressImpl::PrintDone()
+void progress_bar_impl::print_done()
 {
-	std::chrono::duration<double> elapsed = std::chrono::system_clock::now() - mStart;
+	std::chrono::duration<double> elapsed = std::chrono::system_clock::now() - m_start;

 	std::ostringstream msgstr;
-	msgstr << mAction << " done in " << elapsed << " seconds";
+	msgstr << m_action << " done in " << elapsed << " seconds";
 	auto msg = msgstr.str();

 	uint32_t width = get_terminal_width();
@@ -340,46 +354,34 @@ void ProgressImpl::PrintDone()
 	std::cout << '\r' << msg << std::endl;
 }

-Progress::Progress(int64_t inMax, const std::string &inAction)
+progress_bar::progress_bar(int64_t inMax, const std::string &inAction)
 	: m_impl(nullptr)
 {
 	if (isatty(STDOUT_FILENO) and VERBOSE >= 0)
-		m_impl = new ProgressImpl(inMax, inAction);
+		m_impl = new progress_bar_impl(inMax, inAction);
 }

-Progress::~Progress()
+progress_bar::~progress_bar()
 {
-	if (m_impl != nullptr)
-		m_impl->Stop();
-
 	delete m_impl;
 }

-void Progress::consumed(int64_t inConsumed)
-{
-	if (m_impl != nullptr and
-		(m_impl->mConsumed += inConsumed) >= m_impl->mMax)
-	{
-		m_impl->Stop();
-	}
-}
-
-void Progress::progress(int64_t inProgress)
-{
-	if (m_impl != nullptr and
-		(m_impl->mConsumed = inProgress) >= m_impl->mMax)
-	{
-		m_impl->Stop();
-	}
-}
-
-void Progress::message(const std::string &inMessage)
+void progress_bar::consumed(int64_t inConsumed)
 {
 	if (m_impl != nullptr)
-	{
-		std::unique_lock lock(m_impl->mMutex);
-		m_impl->mMessage = inMessage;
-	}
+		m_impl->consumed(inConsumed);
+}
+
+void progress_bar::progress(int64_t inProgress)
+{
+	if (m_impl != nullptr)
+		m_impl->progress(inProgress);
+}
+
+void progress_bar::message(const std::string &inMessage)
+{
+	if (m_impl != nullptr)
+		m_impl->message(inMessage);
 }

 } // namespace cif
@@ -821,12 +823,12 @@ namespace cif

 // --------------------------------------------------------------------

-class ResourcePool
+class resource_pool
 {
  public:
-	static ResourcePool &instance()
+	static resource_pool &instance()
 	{
-		static std::unique_ptr<ResourcePool> s_instance(new ResourcePool);
+		static std::unique_ptr<resource_pool> s_instance(new resource_pool);
 		return *s_instance;
 	}

@@ -856,7 +858,7 @@ class ResourcePool
 	std::unique_ptr<std::istream> load(fs::path name);

  private:
-	ResourcePool();
+	resource_pool();

 	std::unique_ptr<std::ifstream> open(fs::path &p)
 	{
@@ -882,7 +884,7 @@ class ResourcePool
 	std::deque<fs::path> mDirs;
 };

-ResourcePool::ResourcePool()
+resource_pool::resource_pool()
 {
 #if defined(DATA_DIR)
 	pushDir(DATA_DIR);
@@ -899,7 +901,7 @@ ResourcePool::ResourcePool()
 #endif
 }

-std::unique_ptr<std::istream> ResourcePool::load(fs::path name)
+std::unique_ptr<std::istream> resource_pool::load(fs::path name)
 {
 	std::unique_ptr<std::istream> result;
 	std::error_code ec;
@@ -909,6 +911,9 @@ std::unique_ptr<std::istream> ResourcePool::load(fs::path name)
 	if (mLocalResources.count(name.string()))
 		result = open(mLocalResources[name.string()]);

+	if (fs::exists(p, ec) and not ec)
+		result = open(p);
+
 	for (auto di = mDirs.begin(); not result and di != mDirs.end(); ++di)
 	{
 		auto p2 = *di / p;
@@ -931,17 +936,17 @@ std::unique_ptr<std::istream> ResourcePool::load(fs::path name)

 void add_data_directory(std::filesystem::path dataDir)
 {
-	ResourcePool::instance().pushDir(dataDir);
+	resource_pool::instance().pushDir(dataDir);
 }

 void add_file_resource(const std::string &name, std::filesystem::path dataFile)
 {
-	ResourcePool::instance().pushAlias(name, dataFile);
+	resource_pool::instance().pushAlias(name, dataFile);
 }

 std::unique_ptr<std::istream> load_resource(std::filesystem::path name)
 {
-	return ResourcePool::instance().load(name);
+	return resource_pool::instance().load(name);
 }

 } // namespace cif
--- a/src/validate.cpp
+++ b/src/validate.cpp
@@ -24,10 +24,10 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <cif++/dictionary_parser.hpp>
-#include <cif++/validate.hpp>
-#include <cif++/utilities.hpp>
-#include <cif++/gzio.hpp>
+#include "cif++/validate.hpp"
+#include "cif++/dictionary_parser.hpp"
+#include "cif++/gzio.hpp"
+#include "cif++/utilities.hpp"

 #include <cassert>
 #include <fstream>
@@ -401,87 +401,94 @@ void validator::report_error(const std::string &msg, bool fatal) const

 const validator &validator_factory::operator[](std::string_view dictionary_name)
 {
-	std::lock_guard lock(m_mutex);
-
-	for (auto &validator : m_validators)
+	try
 	{
-		if (iequals(validator.name(), dictionary_name))
-			return validator;
-	}
-
-	// not found, try to see if it helps if we tweak the name a little
-
-	// too bad clang version 10 did not have a constructor for std::filesystem::path that accepts a std::string_view
-	std::filesystem::path dictionary(dictionary_name.data(), dictionary_name.data() + dictionary_name.length());
-
-	if (dictionary.extension() != ".dic")
-	{
-		auto dict_name = dictionary.filename().string() + ".dic";
+		std::lock_guard lock(m_mutex);

 		for (auto &validator : m_validators)
 		{
-			if (iequals(validator.name(), dict_name))
+			if (iequals(validator.name(), dictionary_name))
 				return validator;
 		}
-	}

-	// not found, add it
+		// not found, try to see if it helps if we tweak the name a little

+		// too bad clang version 10 did not have a constructor for std::filesystem::path that accepts a std::string_view
+		std::filesystem::path dictionary(dictionary_name.data(), dictionary_name.data() + dictionary_name.length());

-	auto data = load_resource(dictionary_name);
-
-	if (not data and dictionary.extension().string() != ".dic")
-		data = load_resource(dictionary.parent_path() / (dictionary.filename().string() + ".dic"));
-
-	if (data)
-		construct_validator(dictionary_name, *data);
-	else
-	{
-		std::error_code ec;
-
-		// might be a compressed dictionary on disk
-		std::filesystem::path p = dictionary;
-		if (p.extension() == ".dic")
-			p = p.parent_path() / (p.filename().string() + ".gz");
-		else
-			p = p.parent_path() / (p.filename().string() + ".dic.gz");
-
-#if defined(CACHE_DIR) or defined(DATA_DIR)
-		if (not std::filesystem::exists(p, ec) or ec)
+		if (dictionary.extension() != ".dic")
 		{
-			for (const char *dir : {
-#if defined(CACHE_DIR)
-					 CACHE_DIR,
-#endif
-#if defined(DATA_DIR)
-						 DATA_DIR
-#endif
-				 })
+			auto dict_name = dictionary.filename().string() + ".dic";
+
+			for (auto &validator : m_validators)
 			{
-				auto p2 = std::filesystem::path(dir) / p;
-				if (std::filesystem::exists(p2, ec) and not ec)
-				{
-					swap(p, p2);
-					break;
-				}
+				if (iequals(validator.name(), dict_name))
+					return validator;
 			}
 		}
+
+		// not found, add it
+		auto data = load_resource(dictionary_name);
+
+		if (not data and dictionary.extension().string() != ".dic")
+			data = load_resource(dictionary.parent_path() / (dictionary.filename().string() + ".dic"));
+
+		if (data)
+			construct_validator(dictionary_name, *data);
+		else
+		{
+			std::error_code ec;
+
+			// might be a compressed dictionary on disk
+			std::filesystem::path p = dictionary;
+			if (p.extension() == ".dic")
+				p = p.parent_path() / (p.filename().string() + ".gz");
+			else
+				p = p.parent_path() / (p.filename().string() + ".dic.gz");
+
+#if defined(CACHE_DIR) or defined(DATA_DIR)
+			if (not std::filesystem::exists(p, ec) or ec)
+			{
+				for (const char *dir : {
+#if defined(CACHE_DIR)
+						 CACHE_DIR,
+#endif
+#if defined(DATA_DIR)
+							 DATA_DIR
+#endif
+					 })
+				{
+					auto p2 = std::filesystem::path(dir) / p;
+					if (std::filesystem::exists(p2, ec) and not ec)
+					{
+						swap(p, p2);
+						break;
+					}
+				}
+			}
 #endif

-		if (std::filesystem::exists(p, ec) and not ec)
-		{
-			gzio::ifstream in(p);
+			if (std::filesystem::exists(p, ec) and not ec)
+			{
+				gzio::ifstream in(p);

-			if (not in.is_open())
-				throw std::runtime_error("Could not open dictionary (" + p.string() + ")");
+				if (not in.is_open())
+					throw std::runtime_error("Could not open dictionary (" + p.string() + ")");

-			construct_validator(dictionary_name, in);
+				construct_validator(dictionary_name, in);
+			}
+			else
+				throw std::runtime_error("Dictionary not found or defined (" + dictionary.string() + ")");
 		}
-		else
-			throw std::runtime_error("Dictionary not found or defined (" + dictionary.string() + ")");
-	}

-	return m_validators.back();
+		return m_validators.back();
+	}
+	catch (const std::exception &ex)
+	{
+		std::string msg = "Error while loading dictionary ";
+		msg += dictionary_name;
+		std::throw_with_nested(std::runtime_error(msg));
+	}
 }

 void validator_factory::construct_validator(std::string_view name, std::istream &is)
--- a/test/2bi3.cif.gz
+++ b/test/2bi3.cif.gz
--- a/test/3bwh.cif.gz
+++ b/test/3bwh.cif.gz
--- a/test/4wvp.cif.gz
+++ b/test/4wvp.cif.gz
--- a/test/model-test.cpp
+++ b/test/model-test.cpp
@@ -200,6 +200,123 @@ _atom_type.symbol   C
 	}
 }

+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_nonpoly_2)
+{
+    cif::VERBOSE = 1;
+
+	cif::file file;
+	file.load_dictionary("mmcif_pdbx.dic");
+	file.emplace("TEST");	// create a datablock
+	
+	cif::mm::structure structure(file);
+
+	cif::file lig(gTestDir / "HEM.cif");
+	auto &chem_comp_atom = lig["HEM"]["chem_comp_atom"];
+
+	std::vector<cif::row_initializer> atoms;
+
+	for (const auto &[type_symbol, label_atom_id, Cartn_x, Cartn_y, Cartn_z] :
+		chem_comp_atom.rows<std::string,std::string,float,float,float>(
+				"type_symbol", "atom_id", "model_Cartn_x", "model_Cartn_y", "model_Cartn_z"))
+	{
+		atoms.emplace_back(cif::row_initializer{
+			{ "type_symbol", type_symbol },
+			{ "label_atom_id", label_atom_id },
+			{ "auth_atom_id", label_atom_id },
+			{ "Cartn_x", Cartn_x },
+			{ "Cartn_y", Cartn_y },
+			{ "Cartn_z", Cartn_z }
+		});
+
+		if (atoms.size() == 4)
+			break;
+	}
+
+	std::string entity_id = structure.create_non_poly_entity("HEM");
+	structure.create_non_poly(entity_id, atoms);
+
+	auto expected = R"(
+data_TEST
+# 
+_pdbx_nonpoly_scheme.asym_id         A 
+_pdbx_nonpoly_scheme.ndb_seq_num     1 
+_pdbx_nonpoly_scheme.entity_id       1 
+_pdbx_nonpoly_scheme.mon_id          HEM 
+_pdbx_nonpoly_scheme.pdb_seq_num     1 
+_pdbx_nonpoly_scheme.auth_seq_num    1 
+_pdbx_nonpoly_scheme.pdb_mon_id      HEM 
+_pdbx_nonpoly_scheme.auth_mon_id     HEM 
+_pdbx_nonpoly_scheme.pdb_strand_id   A 
+_pdbx_nonpoly_scheme.pdb_ins_code    . 
+#
+loop_
+_atom_site.id
+_atom_site.auth_asym_id
+_atom_site.label_alt_id
+_atom_site.label_asym_id
+_atom_site.label_atom_id
+_atom_site.label_comp_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.type_symbol
+_atom_site.group_PDB
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+1 A ? A CHA HEM 1 . C HETATM ? 2.748 -19.531 39.896 1.00 ? 1 HEM CHA 1
+2 A ? A CHB HEM 1 . C HETATM ? 3.258 -17.744 35.477 1.00 ? 1 HEM CHB 1
+3 A ? A CHC HEM 1 . C HETATM ? 1.703 -21.9   33.637 1.00 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? 1.149 -23.677 38.059 1.00 ? 1 HEM CHD 1
+#
+_chem_comp.id               HEM
+_chem_comp.type             NON-POLYMER
+_chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
+_chem_comp.formula          'C34 H32 Fe N4 O4'
+_chem_comp.formula_weight   616.487000
+#
+_pdbx_entity_nonpoly.entity_id   1
+_pdbx_entity_nonpoly.name        'PROTOPORPHYRIN IX CONTAINING FE'
+_pdbx_entity_nonpoly.comp_id     HEM
+#
+_entity.id                 1
+_entity.type               non-polymer
+_entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
+_entity.formula_weight     616.487000
+#
+_struct_asym.id                            A
+_struct_asym.entity_id                     1
+_struct_asym.pdbx_blank_PDB_chainid_flag   N
+_struct_asym.pdbx_modified                 N
+_struct_asym.details                       ?
+#
+_atom_type.symbol   C
+)"_cf;
+
+	expected.load_dictionary("mmcif_pdbx.dic");
+
+	if (not (expected.front() == structure.get_datablock()))
+	{
+		BOOST_TEST(false);
+		std::cout << expected.front() << std::endl
+				<< std::endl
+				<< structure.get_datablock() << std::endl;
+		
+
+		expected.save("/tmp/a");
+		file.save("/tmp/b");
+	}
+}
+
+
 // --------------------------------------------------------------------

 BOOST_AUTO_TEST_CASE(test_atom_id)
--- a/test/spinner-test.cpp
+++ b/test/spinner-test.cpp
@@ -0,0 +1,49 @@
+#include <cif++/utilities.hpp>
+
+#include <random>
+#include <thread>
+
+void test_one()
+{
+	std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> distrib(100, 1000);
+
+	cif::progress_bar pb(10, "test");
+
+	for (int i = 0; i < 10; ++i)
+	{
+		std::this_thread::sleep_for(std::chrono::milliseconds(distrib(gen)));
+
+		pb.message("step " + std::to_string(i));
+		pb.consumed(1);
+	}
+}
+
+void test_two()
+{
+	cif::progress_bar pb(10, "test");
+
+
+	for (int i = 0; i < 5; ++i)
+		pb.consumed(1);
+}
+
+void test_three()
+{
+	using namespace std::literals;
+
+	cif::progress_bar pb(10, "test");
+	pb.consumed(10);
+
+	std::this_thread::sleep_for(100ms);
+}
+
+int main()
+{
+	test_one();
+	test_two();
+	test_three();
+
+	return 0;
+}
--- a/test/unit-3d-test.cpp
+++ b/test/unit-3d-test.cpp
@@ -31,10 +31,10 @@

 #include <cif++.hpp>

-#include <cif++/dictionary_parser.hpp>
-#include <cif++/parser.hpp>
+#include <Eigen/Eigenvalues>

 namespace tt = boost::test_tools;
+namespace utf = boost::unit_test;

 std::filesystem::path gTestDir = std::filesystem::current_path(); // filled in first test

@@ -168,6 +168,46 @@ BOOST_AUTO_TEST_CASE(t3)
 	BOOST_TEST(a == 45, tt::tolerance(0.01));
 }

+BOOST_AUTO_TEST_CASE(dh_q_0)
+{
+	cif::point axis(1, 0, 0);
+
+	cif::point p(1, 1, 0);
+	
+	cif::point t[3] =
+	{
+		{ 0, 1, 0 },
+		{ 0, 0, 0 },
+		{ 1, 0, 0 }
+	};
+
+	auto a = cif::dihedral_angle(t[0], t[1], t[2], p);
+	BOOST_TEST(a == 0, tt::tolerance(0.01f));
+
+	auto q = cif::construct_from_angle_axis(90, axis);
+
+	p.rotate(q);
+
+	BOOST_TEST(p.m_x == 1, tt::tolerance(0.01f));
+	BOOST_TEST(p.m_y == 0, tt::tolerance(0.01f));
+	BOOST_TEST(p.m_z == 1, tt::tolerance(0.01f));
+
+	a = cif::dihedral_angle(t[0], t[1], t[2], p);
+	BOOST_TEST(a == 90, tt::tolerance(0.01f));
+
+	q = cif::construct_from_angle_axis(-90, axis);
+
+	p.rotate(q);
+
+	BOOST_TEST(p.m_x == 1, tt::tolerance(0.01f));
+	BOOST_TEST(p.m_y == 1, tt::tolerance(0.01f));
+	BOOST_TEST(p.m_z == 0, tt::tolerance(0.01f));
+
+	a = cif::dihedral_angle(t[0], t[1], t[2], p);
+	BOOST_TEST(a == 0, tt::tolerance(0.01f));
+
+}
+
 BOOST_AUTO_TEST_CASE(dh_q_1)
 {
 	struct
@@ -204,11 +244,299 @@ BOOST_AUTO_TEST_CASE(dh_q_1)
 	{
 		auto q = cif::construct_for_dihedral_angle(pts[0], pts[1], pts[2], pts[3], angle, 1);

-		pts[3] -= pts[2];
-		pts[3].rotate(q);
-		pts[3] += pts[2];
+		pts[3].rotate(q, pts[2]);

 		auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
 		BOOST_TEST(dh == angle, tt::tolerance(0.1f));
 	}
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(m2q_0, *utf::tolerance(0.001f))
+{
+	for (size_t i = 0; i < cif::kSymopNrTableSize; ++i)
+	{
+		auto d = cif::kSymopNrTable[i].symop().data();
+
+		cif::matrix3x3<float> 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::Matrix4f em;
+
+		em << 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;
+
+		Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
+
+		auto ev = es.eigenvalues();
+
+		size_t bestJ = 0;
+		float bestEV = -1;
+
+		for (size_t j = 0; j < 4; ++j)
+		{
+			if (bestEV < ev[j].real())
+			{
+				bestEV = ev[j].real();
+				bestJ = j;
+			}
+		}
+
+		if (std::abs(bestEV - 1) > 0.01)
+			continue; // not a rotation matrix
+
+		auto col = es.eigenvectors().col(bestJ);
+
+		auto q = normalize(cif::quaternion{
+			static_cast<float>(col(3).real()),
+			static_cast<float>(col(0).real()),
+			static_cast<float>(col(1).real()),
+			static_cast<float>(col(2).real()) });
+		
+		cif::point p1{ 1, 1, 1 };
+		cif::point p2 = p1;
+		p2.rotate(q);
+
+		cif::point p3 = rot * p1;
+
+		BOOST_TEST(p2.m_x == p3.m_x);
+		BOOST_TEST(p2.m_y == p3.m_y);
+		BOOST_TEST(p2.m_z == p3.m_z);
+	}
+}
+
+// BOOST_AUTO_TEST_CASE(m2q_1, *utf::tolerance(0.001f))
+// {
+// 	for (size_t i = 0; i < cif::kSymopNrTableSize; ++i)
+// 	{
+// 		auto d = cif::kSymopNrTable[i].symop().data();
+
+// 		cif::matrix3x3<float> 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];
+
+// 		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
+// 		});
+
+// 		auto &&[ev, em] = cif::eigen(m * (1/3.0f), false);
+
+// 		size_t bestJ = 0;
+// 		float bestEV = -1;
+
+// 		for (size_t j = 0; j < 4; ++j)
+// 		{
+// 			if (bestEV < ev[j])
+// 			{
+// 				bestEV = ev[j];
+// 				bestJ = j;
+// 			}
+// 		}
+
+// 		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.rotate(q);
+
+// 		cif::point p3 = rot * p1;
+
+// 		BOOST_TEST(p2.m_x == p3.m_x);
+// 		BOOST_TEST(p2.m_y == p3.m_y);
+// 		BOOST_TEST(p2.m_z == p3.m_z);
+// 	}
+// }
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(symm_1)
+{
+	cif::cell c(10, 10, 10);
+
+	cif::point p{ 1, 1, 1 };
+
+	cif::point f = fractional(p, c);
+
+	BOOST_TEST(f.m_x == 0.1f, tt::tolerance(0.01));
+	BOOST_TEST(f.m_y == 0.1f, tt::tolerance(0.01));
+	BOOST_TEST(f.m_z == 0.1f, tt::tolerance(0.01));
+
+	cif::point o = orthogonal(f, c);
+
+	BOOST_TEST(o.m_x == 1.f, tt::tolerance(0.01));
+	BOOST_TEST(o.m_y == 1.f, tt::tolerance(0.01));
+	BOOST_TEST(o.m_z == 1.f, tt::tolerance(0.01));
+}
+
+BOOST_AUTO_TEST_CASE(symm_2)
+{
+	using namespace cif::literals;
+
+	auto symop = "1_555"_symop;
+
+	BOOST_TEST(symop.is_identity() == true);
+}
+
+BOOST_AUTO_TEST_CASE(symm_3)
+{
+	using namespace cif::literals;
+
+	cif::spacegroup sg(18);
+
+	BOOST_TEST(sg.size() == 4);
+	BOOST_TEST(sg.get_name() == "P 21 21 2");
+}
+
+BOOST_AUTO_TEST_CASE(symm_4, *utf::tolerance(0.1f))
+{
+	using namespace cif::literals;
+
+	// based on 2b8h
+	auto sg = cif::spacegroup(154); // p 32 2 1
+	auto c = cif::cell(107.516, 107.516, 338.487, 90.00, 90.00, 120.00);
+	
+	cif::point a{   -8.688,  79.351, 10.439 }; // O6 NAG A 500
+	cif::point b{  -35.356,  33.693, -3.236 }; // CG2 THR D 400
+	cif::point sb(  -6.916,   79.34,   3.236); // 4_565 copy of b
+
+	BOOST_TEST(distance(a, sg(a, c, "1_455"_symop)) == static_cast<float>(c.get_a()));
+	BOOST_TEST(distance(a, sg(a, c, "1_545"_symop)) == static_cast<float>(c.get_b()));
+	BOOST_TEST(distance(a, sg(a, c, "1_554"_symop)) == static_cast<float>(c.get_c()));
+
+	auto sb2 = sg(b, c, "4_565"_symop);
+	BOOST_TEST(sb.m_x == sb2.m_x);
+	BOOST_TEST(sb.m_y == sb2.m_y);
+	BOOST_TEST(sb.m_z == sb2.m_z);
+
+	BOOST_TEST(distance(a, sb2) == 7.42f);	
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(symm_4wvp_1, *utf::tolerance(0.1f))
+{
+	using namespace cif::literals;
+
+	cif::file f(gTestDir / "4wvp.cif.gz");
+
+	auto &db = f.front();
+	cif::mm::structure s(db);
+
+	cif::crystal c(db);
+
+	cif::point p{ -78.722, 98.528,  11.994 };
+	auto a = s.get_residue("A", 10, "").get_atom_by_atom_id("O");
+
+	auto sp1 = c.symmetry_copy(a.get_location(), "2_565"_symop);
+	BOOST_TEST(sp1.m_x == p.m_x);
+	BOOST_TEST(sp1.m_y == p.m_y);
+	BOOST_TEST(sp1.m_z == p.m_z);
+
+	const auto &[d, sp2, so] = c.closest_symmetry_copy(p, a.get_location());
+
+	BOOST_TEST(d < 1);
+
+	BOOST_TEST(sp2.m_x == p.m_x);
+	BOOST_TEST(sp2.m_y == p.m_y);
+	BOOST_TEST(sp2.m_z == p.m_z);
+
+}
+
+BOOST_AUTO_TEST_CASE(symm_2bi3_1, *utf::tolerance(0.1f))
+{
+	cif::file f(gTestDir / "2bi3.cif.gz");
+
+	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 sa1 = c.symmetry_copy(a1.get_location(), cif::sym_op(symm1));
+		auto sa2 = c.symmetry_copy(a2.get_location(), cif::sym_op(symm2));
+
+		BOOST_TEST(cif::distance(sa1, sa2) == dist);
+
+		auto pa1 = a1.get_location();
+
+		const auto &[d, p, so] = c.closest_symmetry_copy(pa1, a2.get_location());
+
+		BOOST_TEST(p.m_x == sa2.m_x);
+		BOOST_TEST(p.m_y == sa2.m_y);
+		BOOST_TEST(p.m_z == sa2.m_z);
+
+		BOOST_TEST(d == dist);
+		BOOST_TEST(so.string() == symm2);
+	}
+}
+
+
+BOOST_AUTO_TEST_CASE(symm_3bwh_1, *utf::tolerance(0.1f))
+{
+	cif::file f(gTestDir / "3bwh.cif.gz");
+
+	auto &db = f.front();
+
+	cif::crystal c(db);
+	cif::mm::structure s(db);
+
+	for (auto a1 : s.atoms())
+	{
+		for (auto a2 : s.atoms())
+		{
+			if (a1 == a2)
+				continue;
+			
+			const auto&[ d, p, so ] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
+
+			BOOST_TEST(d == distance(a1.get_location(), p));
+		}
+	}
 }
--- a/test/unit-v2-test.cpp
+++ b/test/unit-v2-test.cpp
@@ -31,8 +31,8 @@

 #include <cif++.hpp>

-#include <cif++/dictionary_parser.hpp>
-#include <cif++/parser.hpp>
+#include "cif++/dictionary_parser.hpp"
+

 namespace tt = boost::test_tools;

@@ -3105,3 +3105,29 @@ _date    today
 	BOOST_TEST(db == db2);
 }

+BOOST_AUTO_TEST_CASE(find1_opt_1)
+{
+	using namespace cif::literals;
+	using namespace std::literals;
+
+	auto f = R"(data_TEST
+#
+loop_
+_test.id
+_test.name
+_test.value
+1 aap   1.0
+2 noot  1.1
+3 mies  1.2
+    )"_cf;
+
+	auto &db = f.front();
+	auto &test = db["test"];
+
+	auto v = test.find1<std::optional<float>>("id"_key == 1, "value");
+	BOOST_CHECK(v.has_value());
+	BOOST_TEST(*v == 1.0f);
+
+	v = test.find1<std::optional<float>>("id"_key == 4, "value");
+	BOOST_CHECK(v.has_value() == false);
+}
Author	SHA1	Message	Date
Maarten L. Hekkelman	bf9bdd2aae	Merge branch 'trunk' of github.com:PDB-REDO/libcifpp into trunk	2023-05-31 15:17:00 +02:00
Maarten L. Hekkelman	ce14593f0b	improved loading resources from absolute path better error reporting when loading dictionary	2023-05-31 15:16:10 +02:00
Maarten L. Hekkelman	1c02a451e1	improved has_atom_id added couple of comparison operators to sym_op class	2023-05-16 13:55:07 +02:00
Maarten L. Hekkelman	448855a2d3	catch error in create entity for branch	2023-05-09 11:46:37 +02:00
Maarten L. Hekkelman	8ac8e89f2b	Fix progress bar by removing conditional variable	2023-05-02 13:45:02 +02:00
Maarten L. Hekkelman	2281f59401	Remove struct_conn records as well in remove_branch	2023-05-02 13:44:36 +02:00
Maarten L. Hekkelman	4cb0673370	small change to matrix	2023-04-25 10:13:30 +02:00
Maarten L. Hekkelman	76c5706f7c	moving to eigen3 eigensolver, fixing include and dependencies	2023-04-22 14:14:48 +02:00
Maarten L. Hekkelman	2bf4284ff4	cleanup	2023-04-21 14:52:12 +02:00
Maarten L. Hekkelman	d9e2fc97f3	Added missing spinner test	2023-04-21 14:50:22 +02:00
Maarten L. Hekkelman	85dfdf4174	Better progress bar	2023-04-21 14:49:54 +02:00
Martin Salinas	1bede3efda	Removed unused argument warning (#36 ) As argument rhs is not being used in that equals (should the equals function always return false), I added that flag so the compiler skips that warning.	2023-04-21 09:18:47 +02:00
Maarten L. Hekkelman	505f0fdd31	oops	2023-04-20 16:33:24 +02:00
Maarten L. Hekkelman	eed7ec3a4a	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2023-04-20 13:38:48 +02:00
Maarten L. Hekkelman	fdb057e0e2	for now, require eigen3 add inverse_symmetry_copy to crystal	2023-04-20 13:14:52 +02:00
Maarten L. Hekkelman	3fddd1a628	Using quaternions, when possible	2023-04-20 11:37:36 +02:00
Maarten L. Hekkelman	2440706b87	backup	2023-04-19 18:51:41 +02:00
Maarten L. Hekkelman	cf628fa95c	backup	2023-04-19 18:36:33 +02:00
Maarten L. Hekkelman	2b0b47d20d	Fix special case	2023-04-19 16:04:59 +02:00
Maarten L. Hekkelman	a8abf2804f	attempt to use quaternions	2023-04-19 16:01:52 +02:00
Maarten L. Hekkelman	22d7757949	Introduced cif::crystal	2023-04-19 10:17:38 +02:00
Maarten L. Hekkelman	0b0d170c96	a bit of documentation	2023-04-19 09:57:49 +02:00
Maarten L. Hekkelman	1e8e9adf62	Merge branch 'trunk' into develop	2023-04-19 09:22:49 +02:00
Maarten L. Hekkelman	0f03fc31e0	added required include	2023-04-19 09:22:32 +02:00
Maarten L. Hekkelman	518432e0fb	test data	2023-04-17 20:54:57 +02:00
Maarten L. Hekkelman	10ef3464ef	Fix symmetry issue	2023-04-17 20:52:10 +02:00
Maarten L. Hekkelman	226abbd577	Merge branch 'develop' of s4.hekkelman.net:git-repo/libcifpp into develop	2023-04-17 18:56:46 +02:00
Maarten L. Hekkelman	8d66f42ab1	more test cases	2023-04-17 18:56:02 +02:00
Maarten L. Hekkelman	0f14d06f9a	Added inverse symmetry operation	2023-04-14 19:38:39 +02:00
Maarten L. Hekkelman	c53be78496	symmetry fixes	2023-04-14 19:04:16 +02:00
Maarten L. Hekkelman	a38f31ce48	fix closest_symmetry_copy	2023-04-14 17:56:59 +02:00
Maarten L. Hekkelman	1258bd5047	eigen, fixed	2023-04-14 14:08:52 +02:00
Maarten L. Hekkelman	d25cbeb14c	matrix eigen value work	2023-04-14 11:47:18 +02:00
Maarten L. Hekkelman	9b60a07fb6	calculating eigen values	2023-04-13 19:55:32 +02:00
Maarten L. Hekkelman	c0dd41ce50	added inverse symmetry operation	2023-04-13 15:49:15 +02:00
Maarten L. Hekkelman	4cff92bbcc	symmetry operations now working correctly	2023-04-13 11:42:59 +02:00
Maarten L. Hekkelman	9aa8a223c7	symmetry work	2023-04-12 17:00:09 +02:00
Maarten L. Hekkelman	fb59adcfdd	Fix symmetry rotational numbers	2023-04-12 10:59:23 +02:00
Maarten L. Hekkelman	4acca8a3e3	Merge branch 'trunk' into develop	2023-04-07 09:31:11 +02:00
Maarten L. Hekkelman	c1030d2b08	Merge branch 'MartinSalinas98-trunk' into trunk	2023-04-07 09:26:00 +02:00
Maarten L. Hekkelman	16a185c6c0	More include changes	2023-04-07 09:16:38 +02:00
Maarten L. Hekkelman	174e818bd0	Merge branch 'trunk' of github.com:MartinSalinas98/libcifpp into MartinSalinas98-trunk	2023-04-07 08:43:44 +02:00
Maarten L. Hekkelman	7f829bf5df	Merge branch 'trunk' of github.com:PDB-REDO/libcifpp into trunk	2023-04-07 08:43:01 +02:00
Maarten L. Hekkelman	71908282bb	merged from trunk	2023-04-07 08:42:46 +02:00
MartinSalinas98	db3ae446af	Imported local files with relative path	2023-04-07 03:39:02 +02:00
Martin Salinas	bc7d291307	Merge branch 'PDB-REDO:trunk' into trunk	2023-04-07 03:33:57 +02:00
Maarten L. Hekkelman	cfd4702279	Fix memory leak	2023-04-05 20:46:18 +02:00
Maarten L. Hekkelman	54eefb546d	Fix memory leak	2023-04-05 20:44:47 +02:00
Maarten L. Hekkelman	6af0d96a4e	Fix memory leak in category	2023-04-05 20:28:47 +02:00
Maarten L. Hekkelman	eb50bee4a3	atom_type_traits changes	2023-04-04 19:19:30 +02:00
Maarten L. Hekkelman	b6143f3652	optimise load atom data	2023-03-30 20:49:03 +02:00
Maarten L. Hekkelman	348aa7afb6	fix test (use gTestDir)	2023-03-30 20:48:41 +02:00
Martin Salinas	66912b68cc	Commented unused variable ``` /home/msalinas/Documents/standalone-installations/cifParsers/libcifpp/src/pdb/cif2pdb.cpp: In function ‘std::tuple<int, int> cif::pdb::WriteCoordinatesForModel(std::ostream&, const cif::datablock&, const std::map<std::__cxx11::basic_string<char>, std::tuple<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, int, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > >&, std::set<std::__cxx11::basic_string<char> >&, int)’: /home/msalinas/Documents/standalone-installations/cifParsers/libcifpp/src/pdb/cif2pdb.cpp:3362:15: warning: unused variable ‘pdbx_nonpoly_scheme’ [-Wunused-variable] 3362 \| auto &pdbx_nonpoly_scheme = db["pdbx_nonpoly_scheme"]; \| ``` This warning appeared while compiling the library. The use of that variable below has been commented, so I think it's appropiate to do the same thing with the unused variable.	2023-03-28 13:54:45 +02:00
Maarten L. Hekkelman	84dd218758	Merge branch 'MartinSalinas98-patch-1' into trunk	2023-03-28 11:33:05 +02:00
Maarten L. Hekkelman	106ae38976	Update readme	2023-03-28 11:32:11 +02:00
Maarten L. Hekkelman	f1a52245ea	Merge branch 'patch-1' of github.com:MartinSalinas98/libcifpp into MartinSalinas98-patch-1	2023-03-28 11:27:48 +02:00
Maarten L. Hekkelman	cea38e5bb2	Merge branch 'trunk' into develop	2023-03-28 10:31:24 +02:00
Maarten L. Hekkelman	ed5aac358c	libcifpp really requires zlib, not only private.	2023-03-28 10:27:58 +02:00
Maarten L. Hekkelman	5eb128251e	Added category::find1<std::optional>	2023-03-27 10:36:47 +02:00
Maarten L. Hekkelman	cfa46ec954	Added model::has_atom_id	2023-03-23 14:32:21 +01:00
Martin Salinas	07cc60e264	Fixed installation commands Installation commands in the readme cause an error when running last command `cmake --install .` because of the lack of sudo privileges. The following commands don't require sudo to run successfully and install the library.	2023-03-23 11:42:10 +01:00
Maarten L. Hekkelman	90973dc547	version bump, update changelog	2023-03-22 12:39:16 +01:00
Maarten L. Hekkelman	12e3d71b00	fix construct_from_angle_axis	2023-03-21 19:49:21 +01:00
Maarten L. Hekkelman	9addc8f873	fix remove_atom add create_water	2023-03-21 19:49:11 +01:00
Maarten L. Hekkelman	343465cef0	Added test for create_non_poly with initializers	2023-03-08 16:00:40 +01:00
Maarten L. Hekkelman	bec5159415	residue numbering in pdb, again...	2023-02-14 08:28:40 +01:00
Maarten L. Hekkelman	f8da8360e6	write twin info in pdb format	2023-02-14 08:27:13 +01:00
Maarten L. Hekkelman	fb2ad7b75d	Fix in REMARK3 parser for more strict mmcif_pdbx dictionary	2023-02-10 16:24:53 +01:00