update changelog

* Use `find_package` for FastFloat prior to `FetchContent_Declare`

* Convert space to tab

CMakeLists.txt

@@ -32,7 +32,7 @@ endif()
 # set the project name
 project(
 	libcifpp
-	VERSION 9.0.4
+	VERSION 9.0.6
 	LANGUAGES CXX C)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
@@ -43,6 +43,7 @@ include(GenerateExportHeader)
 include(CTest)
 include(ExternalProject)
 include(FetchContent)
+include(VersionString)
 
 # When building with ninja-multiconfig, build both debug and release by default
 if(CMAKE_GENERATOR STREQUAL "Ninja Multi-Config")
@@ -129,13 +130,6 @@ if(MSVC)
 	# make msvc standards compliant...
 	add_compile_options(/permissive- /bigobj)
 	add_link_options(/NODEFAULTLIB:library)
-
-	# This is dubious...
-	if(BUILD_SHARED_LIBS)
-		set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>DLL")
-	else()
-		set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>")
-	endif()
 endif()
 
 # Libraries
@@ -168,26 +162,23 @@ if(MSVC)
 	endforeach()
 endif()
 
-# First check if <format> is available
-find_file(FMT NAME format)
-
-if(FMT EQUAL "FMT-NOTFOUND")
-	if(NOT (fmt_FOUND OR TARGET fmt))
-		find_package(fmt REQUIRED)
-		message(FATAL_ERROR "cifpp: <format> not found and neither was libfmt, compiler too old, you're out of luck")
-	endif()
-endif()
-
 # Using fast_float for float parsing, but only if needed
-# try_compile(STD_CHARCONV_COMPILING
-# 	SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/cmake/test-charconv.cpp)
+try_compile(STD_CHARCONV_COMPILING
+	SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/cmake/test-charconv.cpp
+	CXX_STANDARD 20
+	CXX_STANDARD_REQUIRED ON)
 
 if(NOT STD_CHARCONV_COMPILING)
 	message(NOTICE "libcifpp: Using fast_float for std::from_chars")
-	FetchContent_Declare(fastfloat
-		GIT_REPOSITORY "https://github.com/fastfloat/fast_float"
-		GIT_TAG v8.0.2)
-	FetchContent_MakeAvailable(fastfloat)
+	find_package(FastFloat 8.0 QUIET CONFIG)
+	if(NOT FastFloat_FOUND)
+		message(STATUS "FastFloat not found in system, fetching from GitHub")
+		FetchContent_Declare(fastfloat
+			GIT_REPOSITORY "https://github.com/fastfloat/fast_float"
+			GIT_TAG v8.0.2
+			EXCLUDE_FROM_ALL)
+		FetchContent_MakeAvailable(fastfloat)
+	endif()
 endif()
 
 find_package(Threads)
@@ -228,10 +219,6 @@ else()
 	set(EIGEN_INCLUDE_DIR ${SOURCE_DIR})
 endif()
 
-# Create a revision file, containing the current git version info
-include(VersionString)
-write_version_header(${CMAKE_CURRENT_SOURCE_DIR}/src/ LIB_NAME "LibCIFPP")
-
 # SymOp data table
 if(CIFPP_RECREATE_SYMOP_DATA)
 	# The tool to create the table
@@ -253,6 +240,9 @@ if(CIFPP_RECREATE_SYMOP_DATA)
 		"$ENV{CLIBD}/symop.lib")
 endif()
 
+# Create a revision file, containing the current git version info
+write_version_header("${CMAKE_CURRENT_SOURCE_DIR}/src/" LIB_NAME "LibCIFPP")
+
 # Sources
 set(project_sources
 	${CMAKE_CURRENT_SOURCE_DIR}/src/category.cpp
@@ -356,11 +346,7 @@ else()
 endif()
 
 if(NOT STD_CHARCONV_COMPILING)
-	target_link_libraries(cifpp PUBLIC FastFloat::fast_float)
-endif()
-
-if(fmt_FOUND)
-	target_link_libraries(cifpp PUBLIC fmt)
+	target_link_libraries(cifpp PRIVATE FastFloat::fast_float)
 endif()
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
@@ -499,7 +485,6 @@ if(CIFPP_DATA_DIR AND CIFPP_DOWNLOAD_CCD)
 endif()
 
 set(CONFIG_TEMPLATE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cmake/cifpp-config.cmake.in)
-set(REQUIRE_FMT ${fmt_FOUND})
 
 configure_package_config_file(
 	${CONFIG_TEMPLATE_FILE} ${CMAKE_CURRENT_BINARY_DIR}/cifpp/cifpp-config.cmake

changelog

@@ -1,3 +1,12 @@
+Version 9.0.6
+- Various small fixes
+
+Version 9.0.5
+- Added exists to compound_factory
+- Added sub_matrix, fix and extend determinant calculation
+- Added yet another structure::create_non_poly
+- Remove revision.hpp file in make clean (new VersionString.cmake)
+
 Version 9.0.4
 - Fix various stopping and reconstruction errors
 

cmake/VersionString.cmake

@@ -238,7 +238,7 @@ function(write_version_header dir)
 		if(res EQUAL 0)
 			set(REVISION_STRING "${out}")
 		else()
-			message(STATUS "Git hash not found, does this project has a 'build' tag?")
+			message(STATUS "Git hash not found, does this project have a 'build' tag?")
 		endif()
 	else()
 		message(STATUS "Git hash not found")

cmake/cifpp-config.cmake.in

@@ -8,8 +8,5 @@ include(CMakeFindDependencyMacro)
 
 find_dependency(Threads)
 find_dependency(ZLIB REQUIRED)
-if(@REQUIRE_FMT@)
-    find_dependency(fmt REQUIRED)
-endif()
 
 check_required_components(cifpp)

include/cif++/compound.hpp

@@ -179,7 +179,7 @@ class compound
 	friend class compound_factory_impl;
 	friend class local_compound_factory_impl;
 
-	compound(cif::datablock &db);
+	compound(datablock &db);
 
 	std::string m_id;
 	std::string m_name;
@@ -270,11 +270,15 @@ class compound_factory
 		return is_std_base(res_name) or is_std_peptide(res_name);
 	}
 
+	/// Return whether @a res_name is water
 	bool is_water(std::string_view res_name) const
 	{
 		return res_name == "HOH" or res_name == "H2O" or res_name == "WAT";
 	}
 
+	/// Return whether @a res_name already exists, without creating it.
+	bool exists(std::string_view res_name) const;
+
 	/// \brief Create the compound object for \a id
 	///
 	/// This will create the compound instance for \a id if it doesn't exist already.
@@ -331,14 +335,14 @@ class compound_factory
 class compound_source
 {
   public:
-	compound_source(const cif::file &file)
+	compound_source(const file &file)
 	{
-		cif::compound_factory::instance().push_dictionary(file);
+		compound_factory::instance().push_dictionary(file);
 	}
 
 	~compound_source()
 	{
-		cif::compound_factory::instance().pop_dictionary();
+		compound_factory::instance().pop_dictionary();
 	}
 };
 

include/cif++/datablock.hpp

@@ -43,7 +43,7 @@ namespace cif
 
 /**
  * @brief A datablock is a list of category objects with some additional features
- * 
+ *
  */
 
 class datablock : public std::list<category>
@@ -53,7 +53,7 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Construct a new datablock object with name @a name
-	 * 
+	 *
 	 * @param name The name for the new datablock
 	 */
 	datablock(std::string_view name)
@@ -80,7 +80,7 @@ class datablock : public std::list<category>
 	{
 		std::swap(a.m_name, b.m_name);
 		std::swap(a.m_validator, b.m_validator);
-		std::swap(static_cast<std::list<category>&>(a), static_cast<std::list<category>&>(b));
+		std::swap(static_cast<std::list<category> &>(a), static_cast<std::list<category> &>(b));
 	}
 
 	// --------------------------------------------------------------------
@@ -92,7 +92,7 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Set the name of this datablock to @a name
-	 * 
+	 *
 	 * @param name The new name
 	 */
 	void set_name(std::string_view name)
@@ -102,27 +102,27 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Attempt to load the dictionary specified in audit_conform category
-	 * 
+	 *
 	 */
 	void load_dictionary();
 
 	/**
 	 * @brief Set the validator object to @a v
-	 * 
+	 *
 	 * @param v The new validator object, may be null
 	 */
 	void set_validator(const validator *v);
 
 	/**
 	 * @brief Get the validator object
-	 * 
+	 *
 	 * @return const validator* The validator or nullptr if there is none
 	 */
 	const validator *get_validator() const;
 
 	/**
 	 * @brief Validates the content of this datablock and all its content
-	 * 
+	 *
 	 * @return true If the content is valid
 	 * @return false If the content is not valid
 	 */
@@ -131,7 +131,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Validates all contained data for valid links between parents and children
 	 * as defined in the validator
-	 * 
+	 *
 	 * @return true If all links are valid
 	 * @return false If all links are not valid
 	 */
@@ -140,7 +140,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Strip removes all categories and items that are invalid according
 	 * to the assigned validator. Will also add a valid audit_conform block.
-	 * 
+	 *
 	 * @return true if the remaining datablock is valid
 	 */
 	bool strip();
@@ -150,7 +150,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return the category named @a name, will create a new and empty
 	 * category named @a name if it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category to return
 	 * @return category& Reference to the named category
 	 */
@@ -159,7 +159,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return the const category named @a name, will return a reference
 	 * to a static empty category if it was not found.
-	 * 
+	 *
 	 * @param name The name of the category to return
 	 * @return category& Reference to the named category
 	 */
@@ -168,7 +168,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return a pointer to the category named @a name or nullptr if
 	 * it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category
 	 * @return category* Pointer to the category found or nullptr
 	 */
@@ -177,13 +177,12 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return a pointer to the category named @a name or nullptr if
 	 * it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category
 	 * @return category* Pointer to the category found or nullptr
 	 */
 	const category *get(std::string_view name) const;
 
-
 	/**
 	 * @brief Return true if this datablock contains a non-empty category
 	 */
@@ -197,7 +196,7 @@ class datablock : public std::list<category>
 	 * new one if it is not found. The result is a tuple of an iterator
 	 * pointing to the category and a boolean indicating whether the category
 	 * was created or not.
-	 * 
+	 *
 	 * @param name The name for the category
 	 * @return std::tuple<iterator, bool> A tuple containing an iterator pointing
 	 * at the category and a boolean indicating whether the category was newly

include/cif++/matrix.hpp

@@ -124,6 +124,23 @@ class matrix_expression
 
 		return os;
 	}
+
+	template <typename M2>
+	constexpr bool operator==(const matrix_expression<M2> &m) const
+	{
+		bool same = false;
+		if (dim_m() == m.dim_m() and dim_n() == m.dim_n())
+		{
+			same = true;
+			for (std::size_t i = 0; same and i < m.dim_m(); ++i)
+			{
+				for (std::size_t j = 0; same and j < m.dim_n(); ++j)
+					same = operator()(i, j) == m(i, j);
+			}
+		}
+
+		return same;
+	}
 };
 
 // --------------------------------------------------------------------
@@ -594,6 +611,35 @@ auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
 
 // --------------------------------------------------------------------
 
+template <typename M2>
+class sub_matrix : public matrix_expression<sub_matrix<M2>>
+{
+  public:
+	sub_matrix(const M2 &m, int i, int j)
+		: m_m(m)
+		, m_i(i)
+		, m_j(j)
+	{
+	}
+
+	constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m
+	constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(
+			i >= m_i ? i + 1 : i,
+			j >= m_j ? j + 1 : j);
+	}
+
+  private:
+	const M2 &m_m;
+	std::size_t m_i, m_j;
+};
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the determinant of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3
@@ -605,11 +651,23 @@ 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)));
+	return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) +
+			m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) +
+			m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))));
 }
 
+/** Implementation of the determinant function for fixed size matrices of size 4x4 */
+template <typename F = float>
+F determinant(const matrix4x4<F> &m)
+{
+	return m(0, 0) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 0))) -
+	       m(0, 1) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 1))) +
+	       m(0, 2) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 2))) -
+	       m(0, 3) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 3)));
+}
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the inverse of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3

include/cif++/model.hpp

@@ -29,6 +29,7 @@
 #include "cif++/atom_type.hpp"
 #include "cif++/datablock.hpp"
 #include "cif++/point.hpp"
+#include "cif++/row.hpp"
 
 #include <memory>
 #include <numeric>
@@ -134,14 +135,20 @@ class atom
 
 		row_handle row_aniso()
 		{
+			row_handle result{};
 			auto cat = m_db.get("atom_site_anisotrop");
-			return cat ? cat->operator[]({ { "id", m_id } }) : row_handle{};
+			if (cat)
+				result = cat->operator[]({ { "id", m_id } });
+			return result;
 		}
 
 		const row_handle row_aniso() const
 		{
+			row_handle result{};
 			auto cat = m_db.get("atom_site_anisotrop");
-			return cat ? cat->operator[]({ { "id", m_id } }) : row_handle{};
+			if (cat)
+				result = cat->operator[]({ { "id", m_id } });
+			return result;
 		}
 
 		const datablock &m_db;
@@ -1059,6 +1066,14 @@ 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 NonPolymer struct_asym for a compound of type \a compound_id, returns asym_id.
+	/// This method creates new atom records filled with info from the CCD compound info.
+	///
+	/// \param compound_id	 The compound ID of the new nonpoly
+	/// \param skip_hydrogen Do not create hydrogen atoms when true
+	/// \return				 The newly create asym ID
+	std::string create_non_poly(const std::string &compound_id, bool skip_hydrogen);
+
 	/// \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.
 	///

include/cif++/pdb.hpp

@@ -64,7 +64,7 @@ file read(std::istream &is);
  * @brief Read a file in legacy PDB format from std::istream @a is and
  * put the data into @a cifFile
  */
-file read_pdb_file(std::istream &pdbFile);
+void read_pdb_file(std::istream &pdbFile, cif::file &cifFile);
 
 // mmCIF to PDB
 

include/cif++/point.hpp

@@ -30,7 +30,9 @@
 #include <cmath>
 #include <complex>
 #include <cstdint>
+#include <format>
 #include <functional>
+#include <optional>
 #include <valarray>
 
 #if __has_include(<clipper/core/coords.h>)
@@ -365,11 +367,18 @@ class quaternion_type
 	}
 
 	/// \brief test for all zero values
-	constexpr operator bool() const
+	constexpr explicit operator bool() const
 	{
 		return a != 0 or b != 0 or c != 0 or d != 0;
 	}
 
+	/// \brief for debugging e.g.
+	friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs)
+	{
+		os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d);
+		return os;
+	}
+
   private:
 	value_type a, b, c, d;
 };
@@ -743,6 +752,55 @@ inline constexpr auto cross_product(const point_type<F1> &a, const point_type<F2
 		a.m_x * b.m_y - b.m_x * a.m_y);
 }
 
+/// \brief return the squared norm of point @a p
+template <typename F>
+constexpr F norm_squared(const point_type<F> &p)
+{
+	return p.m_x * p.m_x + p.m_y * p.m_y + p.m_z * p.m_z;
+}
+
+/// \brief return the norm of point @a p
+template <typename F>
+constexpr point_type<F> norm(const point_type<F> &p)
+{
+	return std::sqrt(norm_squared(p));
+}
+
+/// \brief return the point where two lines intersect, or an empty value if they don't intersect at all
+template <typename F>
+std::optional<cif::point> line_line_intersection(const point_type<F> &p1,
+	const point_type<F> &p2, const point_type<F> &p3, const point_type<F> &p4)
+{
+	auto p13 = p1 - p3;
+	auto p43 = p4 - p3;
+	if (std::abs(p43.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_z) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto p21 = p2 - p1;
+	if (std::abs(p21.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_z) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto d1343 = cif::dot_product(p43, p13);
+	auto d4321 = cif::dot_product(p43, p21);
+	auto d1321 = cif::dot_product(p13, p21);
+	auto d4343 = cif::dot_product(p43, p43);
+	auto d2121 = cif::dot_product(p21, p21);
+
+	auto denom = d2121 * d4343 - d4321 * d4321;
+	if (std::abs(denom) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto numer = d1343 * d4321 - d1321 * d4343;
+
+	auto mua = numer / denom;
+	auto mub = (d1343 + d4321 * mua) / d4343;
+
+	auto pa = p1 + mua * p21;
+	auto pb = p3 + mub * p43;
+
+	return { (pa + pb) / 2 };
+}
+
 /// \brief return the angle in degrees between the vectors from point @a p2 to @a p1 and @a p2 to @a p3
 template <typename F>
 constexpr auto angle(const point_type<F> &p1, const point_type<F> &p2, const point_type<F> &p3)
@@ -806,6 +864,9 @@ constexpr auto distance_point_to_line(const point_type<F> &l1, const point_type<
 	return cross.length() / line.length();
 }
 
+/// \brief return the smallest sphere around the points in @a pts
+std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts);
+
 // --------------------------------------------------------------------
 /**
  * @brief For e.g. simulated annealing, returns a new point that is moved in

src/compound.cpp

@@ -24,14 +24,38 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "cif++.hpp"
+#include "cif++/compound.hpp" // for compound_atom, compound_bond, compoun...
 
-#include <filesystem>
-#include <fstream>
-#include <map>
-#include <mutex>
-#include <numeric>
-#include <shared_mutex>
+#include "cif++/atom_type.hpp" // for atom_type_traits
+#include "cif++/category.hpp"  // for category
+#include "cif++/datablock.hpp" // for datablock
+#include "cif++/file.hpp"      // for file
+#include "cif++/item.hpp"      // for item
+#include "cif++/iterator.hpp"  // for iterator_proxy
+#include "cif++/parser.hpp"    // for parser
+#include "cif++/point.hpp"     // for distance, point
+#include "cif++/row.hpp"       // for tie, row_initializer, tie_wrap
+#include "cif++/text.hpp"      // for iequals, replace_all, iset
+#include "cif++/utilities.hpp" // for load_resource, VERBOSE, colour_type
+
+#include <algorithm>    // for find_if
+#include <cstddef>      // for size_t
+#include <exception>    // for exception, throw_with_nested
+#include <filesystem>   // for path, exists
+#include <fstream>      // for char_traits, basic_ostream, operator<<
+#include <iomanip>      // for operator<<, quoted
+#include <iostream>     // for clog, cout, cerr
+#include <limits>       // for numeric_limits
+#include <list>         // for _List_iterator
+#include <map>          // for allocator, map, _Rb_tree_iterator
+#include <memory>       // for shared_ptr, unique_ptr, __shared_ptr_...
+#include <optional>     // for optional
+#include <shared_mutex> // for shared_lock, shared_timed_mutex
+#include <stdexcept>    // for runtime_error, invalid_argument, out_...
+#include <string>       // for basic_string, string, operator==, ope...
+#include <string_view>  // for string_view, basic_string_view
+#include <utility>      // for pair, exchange, move
+#include <vector>       // for vector
 
 namespace fs = std::filesystem;
 
@@ -140,7 +164,7 @@ compound::compound(cif::datablock &db)
 
 	cif::tie(m_id, m_name, m_type, m_formula, m_formula_weight, m_formal_charge, one_letter_code, m_parent_id) =
 		chemComp.front().get("id", "name", "type", "formula", "formula_weight", "pdbx_formal_charge", "one_letter_code", "mon_nstd_parent_comp_id");
-	
+
 	if (one_letter_code.length() == 1)
 		m_one_letter_code = one_letter_code.front();
 
@@ -159,7 +183,7 @@ compound::compound(cif::datablock &db)
 		if (stereo_config.empty())
 			atom.stereo_config = stereo_config_type::N;
 		else
-		atom.stereo_config = parse_stereo_config_from_string(stereo_config);
+			atom.stereo_config = parse_stereo_config_from_string(stereo_config);
 		m_atoms.push_back(std::move(atom));
 	}
 
@@ -172,7 +196,7 @@ compound::compound(cif::datablock &db)
 		if (valueOrder.empty())
 			bond.type = bond_type::sing;
 		else
-		bond.type = parse_bond_type_from_string(valueOrder);
+			bond.type = parse_bond_type_from_string(valueOrder);
 		m_bonds.push_back(std::move(bond));
 	}
 }
@@ -231,12 +255,12 @@ float compound::bond_length(const std::string &atomId_1, const std::string &atom
 
 bool compound::is_peptide() const
 {
-	return iequals(m_type, "l-peptide linking")	or iequals(m_type, "peptide linking");
+	return iequals(m_type, "l-peptide linking") or iequals(m_type, "peptide linking");
 }
 
 bool compound::is_base() const
 {
-	return iequals(m_type, "dna linking")	or iequals(m_type, "rna linking");
+	return iequals(m_type, "dna linking") or iequals(m_type, "rna linking");
 }
 
 // --------------------------------------------------------------------
@@ -294,12 +318,31 @@ class compound_factory_impl : public std::enable_shared_from_this<compound_facto
 			delete c;
 	}
 
+	virtual bool exists_self(const std::string &id) const
+	{
+		if (m_missing.contains(id))
+			return false;
+
+		if (std::find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+				{ return c->id() == id; }) != m_compounds.end())
+			return true;
+
+		return m_next and m_next->exists_self(id);
+	}
+
+	bool exists(std::string_view id)
+	{
+		std::shared_lock lock(mMutex);
+
+		return exists_self(std::string{ id });
+	}
+
 	compound *get(std::string id)
 	{
 		std::shared_lock lock(mMutex);
 
 		compound *result = nullptr;
-		
+
 		for (auto impl = shared_from_this(); impl; impl = impl->m_next)
 		{
 			result = impl->create(id);
@@ -363,7 +406,9 @@ compound *compound_factory_impl::create(const std::string &id)
 	if (m_missing.contains(id))
 		return nullptr;
 
-	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end())
+	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+			{ return c->id() == id; });
+		i != m_compounds.end())
 		return *i;
 
 	compound *result = nullptr;
@@ -454,7 +499,6 @@ class local_compound_factory_impl : public compound_factory_impl
 	compound *create(const std::string &id) override;
 
   private:
-
 	compound *construct_compound(const datablock &db, const std::string &id, const std::string &name, const std::string &three_letter_code, const std::string &group);
 
 	cif::file m_local_file;
@@ -465,7 +509,9 @@ compound *local_compound_factory_impl::create(const std::string &id)
 	if (m_missing.contains(id))
 		return nullptr;
 
-	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end())
+	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+			{ return c->id() == id; });
+		i != m_compounds.end())
 		return *i;
 
 	compound *result = nullptr;
@@ -510,12 +556,10 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 
 	float formula_weight = 0;
 	int formal_charge = 0;
-	std::map<std::string,std::size_t> formula_data;
+	std::map<std::string, std::size_t> formula_data;
 
 	for (std::size_t ord = 1; const auto &[atom_id, type_symbol, type, charge, x, y, z, xi, yi, zi] :
-		rdb["chem_comp_atom"].rows<std::string, std::string, std::string, int,
-			std::optional<float>, std::optional<float>, std::optional<float>,
-			std::optional<float>, std::optional<float>, std::optional<float>>(
+		rdb["chem_comp_atom"].rows<std::string, std::string, std::string, int, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>>(
 			"atom_id", "type_symbol", "type", "charge",
 			"model_Cartn_x", "model_Cartn_y", "model_Cartn_z",
 			"pdbx_model_Cartn_x_ideal", "pdbx_model_Cartn_y_ideal", "pdbx_model_Cartn_z_ideal"))
@@ -525,16 +569,14 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 
 		formula_data[type_symbol] += 1;
 
-		db["chem_comp_atom"].emplace({
-			{ "comp_id", id },
+		db["chem_comp_atom"].emplace({ { "comp_id", id },
 			{ "atom_id", atom_id },
 			{ "type_symbol", type_symbol },
 			{ "charge", charge },
-			{ "model_Cartn_x",  x.has_value() ? x : xi, 3 },
-			{ "model_Cartn_y",  y.has_value() ? y : yi, 3 },
-			{ "model_Cartn_z",  z.has_value() ? z : zi, 3 },
-			{ "pdbx_ordinal", ord++ }
-		});
+			{ "model_Cartn_x", x.has_value() ? x : xi, 3 },
+			{ "model_Cartn_y", y.has_value() ? y : yi, 3 },
+			{ "model_Cartn_z", z.has_value() ? z : zi, 3 },
+			{ "pdbx_ordinal", ord++ } });
 
 		formal_charge += charge;
 	}
@@ -551,21 +593,19 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 		else if (cif::iequals(type, "triple") or cif::iequals(type, "trip"))
 			value_order = "TRIP";
 
-		db["chem_comp_bond"].emplace({
-			{ "comp_id", id },
+		db["chem_comp_bond"].emplace({ { "comp_id", id },
 			{ "atom_id_1", atom_id_1 },
 			{ "atom_id_2", atom_id_2 },
 			{ "value_order", value_order },
 			{ "pdbx_aromatic_flag", aromatic },
 			// TODO: fetch stereo_config info from chem_comp_chir
-			{ "pdbx_ordinal", ord++ }
-		});
+			{ "pdbx_ordinal", ord++ } });
 	}
 
 	db.emplace_back(rdb["pdbx_chem_comp_descriptor"]);
 
 	std::string formula;
-	for (bool first = true; const auto &[symbol, count]: formula_data)
+	for (bool first = true; const auto &[symbol, count] : formula_data)
 	{
 		if (std::exchange(first, false))
 			formula += ' ';
@@ -584,15 +624,13 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 	else
 		type = "NON-POLYMER";
 
-	db["chem_comp"].emplace({
-		{ "id", id },
+	db["chem_comp"].emplace({ { "id", id },
 		{ "name", name },
 		{ "type", type },
 		{ "formula", formula },
 		{ "pdbx_formal_charge", formal_charge },
 		{ "formula_weight", formula_weight },
-		{ "three_letter_code", three_letter_code }
-	});
+		{ "three_letter_code", three_letter_code } });
 
 	std::shared_lock lock(mMutex);
 
@@ -698,6 +736,11 @@ void compound_factory::pop_dictionary()
 		m_impl = m_impl->next();
 }
 
+bool compound_factory::exists(std::string_view id) const
+{
+	return m_impl and m_impl->exists(id);
+}
+
 const compound *compound_factory::create(std::string_view id)
 {
 	auto result = m_impl ? m_impl->get(std::string{ id }) : nullptr;
@@ -722,7 +765,7 @@ bool compound_factory::is_peptide(std::string_view res_name) const
 	bool result = is_std_peptide(res_name);
 	if (not result and m_impl)
 	{
-		auto compound = const_cast<compound_factory&>(*this).create(res_name);
+		auto compound = const_cast<compound_factory &>(*this).create(res_name);
 		result = compound != nullptr and compound->is_peptide();
 	}
 	return result;
@@ -734,7 +777,7 @@ bool compound_factory::is_base(std::string_view res_name) const
 	bool result = is_std_base(res_name);
 	if (not result and m_impl)
 	{
-		auto compound = const_cast<compound_factory&>(*this).create(res_name);
+		auto compound = const_cast<compound_factory &>(*this).create(res_name);
 		result = compound != nullptr and compound->is_base();
 	}
 	return result;

src/model.cpp

@@ -24,7 +24,9 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/model.hpp"
 #include "cif++.hpp"
+#include "cif++/point.hpp"
 
 #include <filesystem>
 #include <fstream>
@@ -32,6 +34,7 @@
 #include <iomanip>
 #include <numeric>
 #include <stack>
+#include <stdexcept>
 
 namespace fs = std::filesystem;
 
@@ -349,17 +352,7 @@ std::tuple<point, float> residue::center_and_radius() const
 	for (auto &a : m_atoms)
 		pts.push_back(a.get_location());
 
-	auto center = centroid(pts);
-	float radius = 0;
-
-	for (auto &pt : pts)
-	{
-		float d = static_cast<float>(distance(pt, center));
-		if (radius < d)
-			radius = d;
-	}
-
-	return std::make_tuple(center, radius);
+	return smallest_sphere_around_points(pts);
 }
 
 bool residue::has_alternate_atoms() const
@@ -2346,6 +2339,36 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 	return asym_id;
 }
 
+std::string structure::create_non_poly(const std::string &compound_id, bool skip_hydrogen)
+{
+	auto compound = cif::compound_factory::instance().create(compound_id);
+	if (compound == nullptr)
+		throw std::runtime_error(std::format("{} is not a known compound", compound_id));
+	
+	std::vector<cif::row_initializer> atoms;
+	for (auto a : compound->atoms())
+	{
+		// We skip H-atoms, as fitting without H-atoms works better and we avoid conflicts in protonation states between CCD and MONLIB
+		if (skip_hydrogen and cif::atom_type_traits(a.type_symbol).symbol() == "H")
+			continue;
+
+		auto ax = a.get_location().get_x();
+		auto ay = a.get_location().get_y();
+		auto az = a.get_location().get_z();
+
+		atoms.emplace_back(cif::row_initializer{
+			{ "type_symbol", cif::atom_type_traits(a.type_symbol).symbol() },
+			{ "label_atom_id", a.id },
+			{ "auth_atom_id", a.id },
+			{ "Cartn_x", ax },
+			{ "Cartn_y", ay },
+			{ "Cartn_z", az },
+			{ "B_iso_or_equiv", 30.00 } });
+	}
+
+	return create_non_poly(create_non_poly_entity(compound_id), atoms);
+}
+
 void structure::create_water(row_initializer atom)
 {
 	using namespace literals;

src/pdb/pdb2cif.cpp

@@ -3142,7 +3142,6 @@ void PDBFileParser::ParseRemark350()
 	std::map<std::string, std::string> values;
 	std::vector<std::string> asymIdList;
 	std::smatch m;
-	cif::row_handle genR;
 
 	std::vector<double> mat, vec;
 

src/point.cpp

@@ -25,17 +25,19 @@
  */
 
 #include "cif++/point.hpp"
-#include "cif++/matrix.hpp"
 
-#include <cassert>
-#include <random>
+#include "cif++/matrix.hpp" // for matrix_subtraction, matrix_cofactors
+
+#include <initializer_list>
+#include <random> // for uniform_real_distribution, normal_distri...
+#include <stdexcept>
 
 namespace cif
 {
 
 // --------------------------------------------------------------------
 
-template<typename T>
+template <typename T>
 quaternion_type<T> normalize(quaternion_type<T> q)
 {
 	std::valarray<double> t(4);
@@ -126,10 +128,9 @@ quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4,
 	p4 -= p3;
 	p3 -= p3;
 
-	quaternion q;
 	auto axis = -p2;
-
 	float dh = dihedral_angle(p1, p2, p3, p4);
+
 	return construct_from_angle_axis(angle - dh, axis);
 }
 
@@ -293,9 +294,9 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	}
 
 	quaternion q(
-		static_cast<float>(cf(maxR, 0)), 
-		static_cast<float>(cf(maxR, 1)), 
-		static_cast<float>(cf(maxR, 2)), 
+		static_cast<float>(cf(maxR, 0)),
+		static_cast<float>(cf(maxR, 1)),
+		static_cast<float>(cf(maxR, 2)),
 		static_cast<float>(cf(maxR, 3)));
 	q = normalize(q);
 
@@ -327,4 +328,251 @@ point nudge(point p, float offset)
 	return p + r;
 }
 
+// --------------------------------------------------------------------
+
+std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
+{
+	return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
+}
+
+std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
+{
+	// Find two bisectors
+	auto vz = cross_product(pts[1] - pts[0], pts[2] - pts[0]);
+
+	auto bs1 = cross_product(vz, pts[1] - pts[0]);
+	bs1.normalize();
+
+	auto v1 = (pts[1] - pts[0]);
+	v1.normalize();
+
+	auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
+
+	auto bs2 = cross_product(vz, pts[2] - pts[0]);
+	bs2.normalize();
+
+	auto v2 = (pts[2] - pts[0]);
+	v2.normalize();
+
+	auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
+
+	auto c = line_line_intersection(s1, s1 + bs1, s2, s2 + bs2);
+	if (c)
+		return { *c, distance(*c, pts[0]) };
+
+	// Colinear points I guess, try something else
+	auto l1 = distance_squared(pts[0], pts[1]);
+	auto l2 = distance_squared(pts[0], pts[2]);
+	auto l3 = distance_squared(pts[1], pts[2]);
+
+	if (l1 > l2 and l1 > l3)
+		return smallest_sphere_around_2_points({ pts[0], pts[1] });
+	else if (l2 > l1 and l2 > l3)
+		return smallest_sphere_around_2_points({ pts[0], pts[2] });
+	else
+		return smallest_sphere_around_2_points({ pts[1], pts[2] });
+}
+
+std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
+{
+	auto t0 = -norm_squared(pts[0]);
+	auto t1 = -norm_squared(pts[1]);
+	auto t2 = -norm_squared(pts[2]);
+	auto t3 = -norm_squared(pts[3]);
+
+	// clang-format off
+	matrix4x4<float> Tm({
+		pts[0].m_x, pts[0].m_y, pts[0].m_z, 1,
+		pts[1].m_x, pts[1].m_y, pts[1].m_z, 1,
+		pts[2].m_x, pts[2].m_y, pts[2].m_z, 1,
+		pts[3].m_x, pts[3].m_y, pts[3].m_z, 1
+	});
+	auto T = determinant(Tm);
+
+	if (T != 0)
+	{
+		matrix4x4<float> Dm({
+			t0, pts[0].m_y, pts[0].m_z, 1,
+			t1, pts[1].m_y, pts[1].m_z, 1,
+			t2, pts[2].m_y, pts[2].m_z, 1,
+			t3, pts[3].m_y, pts[3].m_z, 1
+		});
+		auto D = determinant(Dm) / T;
+		
+		matrix4x4<float> Em({
+			pts[0].m_x, t0, pts[0].m_z, 1,
+			pts[1].m_x, t1, pts[1].m_z, 1,
+			pts[2].m_x, t2, pts[2].m_z, 1,
+			pts[3].m_x, t3, pts[3].m_z, 1
+		});
+		auto E = determinant(Em) / T;
+		
+		matrix4x4<float> Fm({
+			pts[0].m_x, pts[0].m_y, t0, 1,
+			pts[1].m_x, pts[1].m_y, t1, 1,
+			pts[2].m_x, pts[2].m_y, t2, 1,
+			pts[3].m_x, pts[3].m_y, t3, 1
+		});
+		
+		auto F = determinant(Fm) / T;
+		
+		matrix4x4<float> Gm({
+			pts[0].m_x, pts[0].m_y, pts[0].m_z, t0,
+			pts[1].m_x, pts[1].m_y, pts[1].m_z, t1,
+			pts[2].m_x, pts[2].m_y, pts[2].m_z, t2,
+			pts[3].m_x, pts[3].m_y, pts[3].m_z, t3
+		});
+		auto G = determinant(Gm) / T;
+		
+		point center{ -D / 2, -E / 2, -F / 2 };
+		float radius = std::sqrt(D * D + E * E + F * F - 4 * G) / 2;
+
+		// clang-format on
+
+		return { center, radius };
+	}
+
+	// Perhaps some colinear points, try something else:
+
+	for (auto ix : std::initializer_list<std::array<size_t, 4>>{
+			 { 1, 2, 3, 0 },
+			 { 0, 2, 3, 1 },
+			 { 0, 1, 3, 2 },
+			 { 0, 1, 2, 3 },
+		 })
+	{
+		auto [center, radius] =
+			smallest_sphere_around_3_points({ pts[ix[0]], pts[ix[1]], pts[ix[2]] });
+
+		if (distance(pts[ix[3]], center) <= radius)
+			return { center, radius };
+	}
+
+	assert(false);
+	exit(1);
+}
+
+std::tuple<point, float> smallest_sphere_around_all_points(std::vector<point> P, std::vector<point> R)
+{
+	if (P.empty() or R.size() == 4)
+	{
+		switch (R.size())
+		{
+			case 1:
+				return { R[0], 0 };
+
+			case 2:
+				return smallest_sphere_around_2_points({ R[0], R[1] });
+
+			case 3:
+				return smallest_sphere_around_3_points({ R[0], R[1], R[2] });
+
+			case 4:
+				return smallest_sphere_around_4_points({ R[0], R[1], R[2], R[3] });
+
+			default:
+				assert(false);
+		}
+	}
+
+	auto p = P.back();
+	P.pop_back();
+
+	auto [c, r] = smallest_sphere_around_all_points(P, R);
+	assert(not std::isnan(r));
+	if (distance(c, p) <= r)
+		return { c, r };
+
+	R.emplace_back(p);
+	return smallest_sphere_around_all_points(P, R);
+}
+
+bool point_in_circle(point p, std::vector<point> c)
+{
+	switch (c.size())
+	{
+		case 0:
+			return false;
+
+		case 1:
+			return p == c.front();
+
+		case 2:
+		{
+			auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		case 3:
+		{
+			auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		case 4:
+		{
+			auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		default:
+			assert(false);
+			throw std::runtime_error("Error finding smallest sphere");
+	}
+}
+
+std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts)
+{
+	std::random_device rd;
+	std::mt19937 g(rd());
+
+	std::shuffle(pts.begin(), pts.end(), g);
+
+	std::vector<size_t> cix;
+
+	auto cirle_points = [&]()
+	{
+		std::vector<point> result;
+		for (auto ix : cix)
+			result.emplace_back(pts[ix]);
+		return result;
+	};
+
+	size_t i = 0;
+	while (i < pts.size())
+	{
+		if (std::find(cix.begin(), cix.end(), i) != cix.end() or
+			point_in_circle(pts[i], cirle_points()))
+		{
+			++i;
+		}
+		else
+		{
+			cix.erase(std::remove_if(cix.begin(), cix.end(), [i](size_t j)
+						  { return j < i; }),
+				cix.end());
+			cix.push_back(i);
+			if (cix.size() < 4)
+				i = 0;
+			else
+				++i;
+		}
+	}
+
+	switch (cix.size())
+	{
+		case 1:
+			return { pts[cix[0]], 0 };
+		case 2:
+			return smallest_sphere_around_2_points({ pts[cix[0]], pts[cix[1]] });
+		case 3:
+			return smallest_sphere_around_3_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]] });
+		case 4:
+			return smallest_sphere_around_4_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]], pts[cix[3]] });
+		default:
+			assert(false);
+			throw std::runtime_error("Error finding smallest sphere");
+	}
+}
+
 } // namespace cif

src/text.cpp

@@ -28,12 +28,12 @@
 
 #include <algorithm>
 #include <cassert>
+#include <charconv>
 
 #if __has_include("fast_float/fast_float.h")
 #include "fast_float/fast_float.h"
 #endif
 
-
 namespace cif
 {
 
@@ -519,19 +519,29 @@ std::vector<std::string> word_wrap(const std::string &text, std::size_t width)
 
 #if __has_include("fast_float/fast_float.h")
 
-template<>
-std::from_chars_result ff_charconv<float>::from_chars(const char *a, const char *b, float &v)
+template <typename T>
+std::from_chars_result ff_charconv<T, typename std::enable_if_t<std::is_floating_point_v<T>>>::from_chars(const char *a, const char *b, T &v)
 {
 	auto r = fast_float::from_chars(a, b, v);
 	return { r.ptr, r.ec };
 }
 
-template<>
-std::from_chars_result ff_charconv<double>::from_chars(const char *a, const char *b, double &v)
-{
-	auto r = fast_float::from_chars(a, b, v);
-	return { r.ptr, r.ec };
-}
+template struct ff_charconv<float>;
+template struct ff_charconv<double>;
+// template struct ff_charconv<long double>;
+
+#ifdef __STDCPP_FLOAT64_T__
+template struct ff_charconv<std::float64_t>;
+#endif
+#ifdef __STDCPP_FLOAT32_T__
+template struct ff_charconv<std::float32_t>;
+#endif
+#ifdef __STDCPP_FLOAT16_T__
+template struct ff_charconv<std::float16_t>;
+#endif
+#ifdef __STDCPP_BFLOAT16_T__
+template struct ff_charconv<std::bfloat16_t>;
+#endif
 
 #endif
 

src/validate.cpp

@@ -27,12 +27,11 @@
 #include "cif++/validate.hpp"
 #include "cif++/category.hpp"
 #include "cif++/dictionary_parser.hpp"
-#include "cif++/gzio.hpp"
 #include "cif++/utilities.hpp"
 
 #include <cassert>
-#include <fstream>
 #include <iostream>
+#include <mutex>
 
 // The validator depends on regular expressions. Unfortunately,
 // the implementation of std::regex in g++ is buggy and crashes
@@ -75,6 +74,7 @@ struct regex_impl
   private:
 	pcre2_code *m_rx = nullptr;
 	pcre2_match_data *m_data = nullptr;
+	mutable std::mutex m_mutex;
 };
 
 regex_impl::regex_impl(std::string_view rx)
@@ -95,6 +95,8 @@ regex_impl::regex_impl(std::string_view rx)
 
 regex_impl::~regex_impl()
 {
+	std::unique_lock lock(m_mutex);
+
 	if (m_data)
 		pcre2_match_data_free(m_data);
 
@@ -104,6 +106,8 @@ regex_impl::~regex_impl()
 
 bool regex_impl::match(std::string_view v) const
 {
+	std::unique_lock lock(m_mutex);
+
 	bool result = false;
 
 	if (int rc = pcre2_match(m_rx, (PCRE2_SPTR)v.data(), v.length(), 0, 0, m_data, nullptr); rc >= 0)

test/CMakeLists.txt

@@ -29,7 +29,7 @@ if(NOT (Catch2_FOUND OR TARGET Catch2))
 		FetchContent_Declare(
 			Catch2
 			GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-            GIT_TAG v3.4.0)
+			GIT_TAG v3.4.0)
 
 		FetchContent_MakeAvailable(Catch2)
 
@@ -49,6 +49,7 @@ list(
 	spinner
 	reconstruction
 	validate-pdbx
+	matrix
 )
 
 add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp")
@@ -70,6 +71,8 @@ foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 		target_compile_options(${CIFPP_TEST} PRIVATE /EHsc)
 	endif()
 
-	add_test(NAME ${CIFPP_TEST}
-		COMMAND $<TARGET_FILE:${CIFPP_TEST}> --data-dir ${CMAKE_CURRENT_SOURCE_DIR})
+	if(NOT (CIFPP_TEST STREQUAL "spinner-test"))
+		add_test(NAME ${CIFPP_TEST}
+			COMMAND $<TARGET_FILE:${CIFPP_TEST}> --data-dir ${CMAKE_CURRENT_SOURCE_DIR})
+	endif()
 endforeach()

test/matrix-test.cpp

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

test/query-test.cpp

@@ -49,6 +49,4 @@ TEST_CASE("q-1")
 
 		CHECK(pdbx_poly_seq_scheme.count("asym_id"_key == "A" and "entity_id"_key == 1 and "seq_id"_key == 1 and "mon_id"_key == "PRO" and "hetero"_key == false) == 1);
 	}
-	
-
 }

test/unit-3d-test.cpp

@@ -24,15 +24,17 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/point.hpp"
 #include "test-main.hpp"
 
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <cif++.hpp>
 #include <stdexcept>
 
-#include <cif++.hpp>
-
 #if defined(_MSC_VER)
-#pragma warning (disable : 5054)	// warning C5054: operator '&': deprecated between enumerations of different types
-#pragma warning (disable : 4127)	// conditional expression is constant
+# pragma warning(disable : 5054) // warning C5054: operator '&': deprecated between enumerations of different types
+# pragma warning(disable : 4127) // conditional expression is constant
 #endif
 
 #include <Eigen/Eigenvalues>
@@ -315,8 +317,7 @@ TEST_CASE("m2q_0a")
 		cif::point p2 = p1;
 		p2.rotate(q);
 
-		cif::matrix3x3<float> rot_c({
-			static_cast<float>(d[0]),
+		cif::matrix3x3<float> rot_c({ static_cast<float>(d[0]),
 			static_cast<float>(d[1]),
 			static_cast<float>(d[2]),
 			static_cast<float>(d[3]),
@@ -324,8 +325,7 @@ TEST_CASE("m2q_0a")
 			static_cast<float>(d[5]),
 			static_cast<float>(d[6]),
 			static_cast<float>(d[7]),
-			static_cast<float>(d[8])
-		});
+			static_cast<float>(d[8]) });
 
 		cif::point p3 = rot_c * p1;
 
@@ -610,3 +610,40 @@ TEST_CASE("volume_3bwh_1")
 
 	CHECK_THAT(c.get_cell().get_volume(), Catch::Matchers::WithinRel(741009.625f, 0.01f));
 }
+
+// --------------------------------------------------------------------
+
+TEST_CASE("smallest_sphere-1")
+{
+	std::vector<cif::point> pts{
+		{ 0.9295, 4.9006, 46.9706 },
+		{ -0.1215, 5.5936, 46.0726 },
+		{ -0.7975, 4.7046, 45.0796 },
+		{ -1.4875, 3.5486, 45.7196 },
+		{ -0.6535, 2.8816, 46.8186 },
+		{ 0.3825, 3.5156, 47.4496 },
+		{ 1.1995, 2.9206, 48.5286 },
+		{ 0.8255, 2.0466, 49.4716 },
+		{ 1.6625, 1.5036, 50.5176 },
+		{ 1.1165, 0.6056, 51.3626 },
+		{ 1.8325, -0.0064, 52.4656 },
+		{ 1.1945, -0.9044, 53.2216 },
+		{ 1.8135, -1.5534, 54.3566 },
+		{ 1.0925, -2.4574, 55.0656 },
+		{ 1.5205, -3.2204, 56.2476 },
+		{ 1.1955, 5.8066, 48.1796 },
+		{ 2.2495, 4.6896, 46.1796 },
+		{ -1.2515, 1.5186, 47.1786 },
+		{ 3.1385, 1.9106, 50.6166 },
+		{ 3.2605, -1.1834, 54.7206 },
+		{ 2.5975, -3.8554, 56.2096 },
+		{ 0.7975, -3.2184, 57.2686 }
+	};
+
+	for (int i = 0; i < 1000; ++i)
+	{
+		auto [c, r] = cif::smallest_sphere_around_points(pts);
+		CHECK_THAT(cif::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
+		CHECK_THAT(r, Catch::Matchers::WithinAbs(7.31248331f, 0.01f));
+	}
+}