Working on item_value

CMakeLists.txt

@@ -32,7 +32,7 @@ endif()
 # set the project name
 project(
 	libcifpp
-	VERSION 9.0.4
+	VERSION 9.0.5
 	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")
@@ -168,25 +169,18 @@ 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()
+set(CMAKE_CXX_STANDARD 20)
 
 # 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)
 
 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)
+		GIT_TAG v8.0.2
+		EXCLUDE_FROM_ALL)
 	FetchContent_MakeAvailable(fastfloat)
 endif()
 
@@ -228,10 +222,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 +243,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
@@ -266,18 +259,18 @@ set(project_sources
 	${CMAKE_CURRENT_SOURCE_DIR}/src/validate.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/src/text.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/src/utilities.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/atom_type.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/compound.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/point.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/symmetry.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/model.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/cif2pdb.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb_record.hpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.hpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/reconstruct.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/validate-pdbx.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/atom_type.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/compound.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/point.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/symmetry.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/model.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/cif2pdb.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb_record.hpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.hpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/reconstruct.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/validate-pdbx.cpp
 )
 
 set(project_headers
@@ -356,11 +349,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 +488,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,9 @@
+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++/category.hpp

@@ -28,12 +28,14 @@
 
 #include "cif++/forward_decl.hpp"
 
-#include "cif++/condition.hpp"
+// #include "cif++/condition.hpp"
+#include "cif++/item.hpp"
 #include "cif++/iterator.hpp"
 #include "cif++/row.hpp"
 #include "cif++/text.hpp"
 
-#include <array>
+// #include <array>
+#include <functional>
 
 /** \file category.hpp
  * Documentation for the cif::category class
@@ -108,16 +110,6 @@ class multiple_results_error : public std::runtime_error
 	}
 };
 
-// --------------------------------------------------------------------
-// These should be moved elsewhere, one day.
-
-/// \cond
-template <typename _Tp>
-inline constexpr bool is_optional_v = false;
-template <typename _Tp>
-inline constexpr bool is_optional_v<std::optional<_Tp>> = true;
-/// \endcond
-
 // --------------------------------------------------------------------
 
 /// The class category is a sequence container for rows of data values.
@@ -181,12 +173,8 @@ class category
 
 	const std::string &name() const { return m_name; } ///< Returns the name of the category
 
-	[[deprecated("use key_items instead")]] iset key_fields() const; ///< Returns the cif::iset of key item names. Retrieved from the @ref category_validator for this category
-
 	iset key_items() const; ///< Returns the cif::iset of key item names. Retrieved from the @ref category_validator for this category
 
-	[[deprecated("use key_item_indices instead")]] std::set<uint16_t> key_field_indices() const; ///< Returns a set of indices for the key items.
-
 	std::set<uint16_t> key_item_indices() const; ///< Returns a set of indices for the key items.
 
 	/// @brief Set the validator for this category to @a v
@@ -336,7 +324,7 @@ class category
 	struct key_element_type
 	{
 		std::string name;         ///< Name of the item
-		std::string value;        ///< Value to be found
+		item_value value;         ///< Value to be found
 		bool may_be_null = false; ///< If true, value should be same or empty
 	};
 
@@ -961,7 +949,7 @@ class category
 			for (auto i = b; i != e; ++i)
 			{
 				// item_value *new_item = this->create_item(*i);
-				r->append(add_item(i->name()), { i->value() });
+				r->append(add_item(i->name()), i->value());
 			}
 		}
 		catch (...)
@@ -1005,7 +993,7 @@ class category
 
 	// --------------------------------------------------------------------
 
-	using value_provider_type = std::function<std::string_view(std::string_view)>;
+	using value_provider_type = std::function<item_value(const item_value &)>;
 
 	/// \brief Update a single item named @a item_name in the rows that match
 	/// \a cond to values provided by a callback function \a value_provider
@@ -1036,7 +1024,7 @@ class category
 	/// That means, child categories are updated if the links are absolute
 	/// and unique. If they are not, the child category rows are split.
 
-	void update_value(condition &&cond, std::string_view item_name, std::string_view value)
+	void update_value(condition &&cond, std::string_view item_name, item_value value)
 	{
 		auto rs = find(std::move(cond));
 		std::vector<row_handle> rows;
@@ -1049,66 +1037,12 @@ class category
 	/// That means, child categories are updated if the links are absolute
 	/// and unique. If they are not, the child category rows are split.
 
-	void update_value(const std::vector<row_handle> &rows, std::string_view item_name, std::string_view value)
+	void update_value(const std::vector<row_handle> &rows, std::string_view item_name, item_value value)
 	{
-		update_value(rows, item_name, [value](std::string_view)
+		update_value(rows, item_name, [value](const item_value &v)
 			{ return value; });
 	}
 
-	// --------------------------------------------------------------------
-	// Naming used to be very inconsistent. For backward compatibility,
-	// the old function names are here as deprecated variants.
-
-	/// \brief Return the index number for \a column_name
-	[[deprecated("Use get_item_ix instead")]] uint16_t get_column_ix(std::string_view column_name) const
-	{
-		return get_item_ix(column_name);
-	}
-
-	/// @brief Return the name for column with index @a ix
-	/// @param ix The index number
-	/// @return The name of the column
-	[[deprecated("use get_item_name instead")]] std::string_view get_column_name(uint16_t ix) const
-	{
-		return get_item_name(ix);
-	}
-
-	/// @brief Make sure a item with name @a item_name is known and return its index number
-	/// @param item_name The name of the item
-	/// @return The index number of the item
-	[[deprecated("use add_item instead")]] uint16_t add_column(std::string_view item_name)
-	{
-		return add_item(item_name);
-	}
-
-	/** @brief Remove column name @a colum_name
-	 * @param column_name The column to be removed
-	 */
-	[[deprecated("use remove_item instead")]] void remove_column(std::string_view column_name)
-	{
-		remove_item(column_name);
-	}
-
-	/** @brief Rename column @a from_name to @a to_name */
-	[[deprecated("use rename_item instead")]] void rename_column(std::string_view from_name, std::string_view to_name)
-	{
-		rename_item(from_name, to_name);
-	}
-
-	/// @brief Return whether a column with name @a name exists in this category
-	/// @param name The name of the column
-	/// @return True if the column exists
-	[[deprecated("use has_item instead")]] bool has_column(std::string_view name) const
-	{
-		return has_item(name);
-	}
-
-	/// @brief Return the cif::iset of columns in this category
-	[[deprecated("use get_items instead")]] iset get_columns() const
-	{
-		return get_items();
-	}
-
 	// --------------------------------------------------------------------
 	/// \brief Return the index number for \a item_name
 
@@ -1117,7 +1051,7 @@ class category
 	/// @brief Return the name for item with index @a ix
 	/// @param ix The index number
 	/// @return The name of the item
-	std::string_view get_item_name(uint16_t ix) const
+	const std::string &get_item_name(uint16_t ix) const
 	{
 		if (ix >= m_items.size())
 			throw std::out_of_range("item index is out of range");
@@ -1200,7 +1134,7 @@ class category
 	}
 
   private:
-	void update_value(row *row, uint16_t item, std::string_view value, bool updateLinked, bool validate = true);
+	void update_value(row *row, uint16_t item, item_value value, bool updateLinked, bool validate = true);
 
 	void erase_orphans(condition &&cond, category &parent);
 

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++/condition.hpp

@@ -371,7 +371,7 @@ namespace detail
 	{
 		key_equals_condition_impl(item &&i)
 			: m_item_name(i.name())
-			, m_value(std::forward<item>(i).value())
+			, m_value(std::forward<item_value>(i.value()))
 		{
 		}
 
@@ -409,7 +409,7 @@ namespace detail
 		std::string m_item_name;
 		uint16_t m_item_ix = 0;
 		bool m_icase = false;
-		std::string m_value;
+		item_value m_value;
 		std::optional<row_handle> m_single_hit;
 	};
 
@@ -466,111 +466,11 @@ namespace detail
 
 		std::string m_item_name;
 		uint16_t m_item_ix = 0;
-		std::string m_value;
+		item_value &m_value;
 		bool m_icase = false;
 		std::optional<row_handle> m_single_hit;
 	};
 
-	struct key_equals_number_condition_impl : public condition_impl
-	{
-		key_equals_number_condition_impl(const std::string &name, double v)
-			: m_item_name(name)
-			, m_value(v)
-		{
-		}
-
-		condition_impl *prepare(const category &c) override;
-
-		bool test(row_handle r) const override
-		{
-			return m_single_hit.has_value() ? *m_single_hit == r : r[m_item_ix].compare(m_value) == 0;
-		}
-
-		void str(std::ostream &os) const override
-		{
-			os << m_item_name << " == " << m_value;
-		}
-
-		virtual std::optional<row_handle> single() const override
-		{
-			return m_single_hit;
-		}
-
-		virtual bool equals(const condition_impl *rhs) const override
-		{
-			if (typeid(*rhs) == typeid(key_equals_number_condition_impl))
-			{
-				auto ri = static_cast<const key_equals_number_condition_impl *>(rhs);
-				if (m_single_hit.has_value() or ri->m_single_hit.has_value())
-					return m_single_hit == ri->m_single_hit;
-				else
-					// watch out, both m_item_ix might be the same while item_names might be diffent (in case they both do not exist in the category)
-					return m_item_ix == ri->m_item_ix and m_value == ri->m_value and m_item_name == ri->m_item_name;
-			}
-			return this == rhs;
-		}
-
-		std::string m_item_name;
-		uint16_t m_item_ix = 0;
-		double m_value;
-		std::optional<row_handle> m_single_hit;
-	};
-
-	struct key_equals_number_or_empty_condition_impl : public condition_impl
-	{
-		key_equals_number_or_empty_condition_impl(key_equals_number_condition_impl *equals)
-			: m_item_name(equals->m_item_name)
-			, m_value(equals->m_value)
-			, m_single_hit(equals->m_single_hit)
-		{
-		}
-
-		condition_impl *prepare(const category &c) override
-		{
-			m_item_ix = get_item_ix(c, m_item_name);
-			return this;
-		}
-
-		bool test(row_handle r) const override
-		{
-			bool result = false;
-			if (m_single_hit.has_value())
-				result = *m_single_hit == r;
-			else
-				result = r[m_item_ix].empty() or r[m_item_ix].compare(m_value) == 0;
-			return result;
-		}
-
-		void str(std::ostream &os) const override
-		{
-			os << '(' << m_item_name << " == " << m_value << " OR " << m_item_name << " IS NULL)";
-		}
-
-		virtual std::optional<row_handle> single() const override
-		{
-			return m_single_hit;
-		}
-
-		virtual bool equals(const condition_impl *rhs) const override
-		{
-			if (typeid(*rhs) == typeid(key_equals_number_or_empty_condition_impl))
-			{
-				auto ri = static_cast<const key_equals_number_or_empty_condition_impl *>(rhs);
-				if (m_single_hit.has_value() or ri->m_single_hit.has_value())
-					return m_single_hit == ri->m_single_hit;
-				else
-					// watch out, both m_item_ix might be the same while item_names might be diffent (in case they both do not exist in the category)
-					return m_item_ix == ri->m_item_ix and m_value == ri->m_value and m_item_name == ri->m_item_name;
-			}
-			return this == rhs;
-		}
-
-		std::string m_item_name;
-		uint16_t m_item_ix = 0;
-		double m_value;
-		std::optional<row_handle> m_single_hit;
-	};
-
 	struct key_compare_condition_impl : public condition_impl
 	{
 		template <typename COMP>
@@ -622,7 +522,7 @@ namespace detail
 
 		bool test(row_handle r) const override
 		{
-			std::string_view txt = r[m_item_ix].text();
+			auto txt = r[m_item_ix].get<std::string>();
 			return std::regex_match(txt.begin(), txt.end(), mRx);
 		}
 
@@ -693,7 +593,7 @@ namespace detail
 			{
 				try
 				{
-					std::string_view txt = r[f].text();
+					auto txt = r[f].get<std::string>();
 					if (std::regex_match(txt.begin(), txt.end(), mRx))
 					{
 						result = true;
@@ -970,26 +870,6 @@ inline condition operator or(condition &&a, condition &&b)
 				return condition(new detail::key_equals_or_empty_condition_impl(ci));
 		}
 
-		if (typeid(*a.m_impl) == typeid(detail::key_equals_number_condition_impl) and
-			typeid(*b.m_impl) == typeid(detail::key_is_empty_condition_impl))
-		{
-			auto ci = static_cast<detail::key_equals_number_condition_impl *>(a.m_impl);
-			auto ce = static_cast<detail::key_is_empty_condition_impl *>(b.m_impl);
-
-			if (ci->m_item_name == ce->m_item_name)
-				return condition(new detail::key_equals_number_or_empty_condition_impl(ci));
-		}
-
-		if (typeid(*b.m_impl) == typeid(detail::key_equals_number_condition_impl) and
-			typeid(*a.m_impl) == typeid(detail::key_is_empty_condition_impl))
-		{
-			auto ci = static_cast<detail::key_equals_number_condition_impl *>(b.m_impl);
-			auto ce = static_cast<detail::key_is_empty_condition_impl *>(a.m_impl);
-
-			if (ci->m_item_name == ce->m_item_name)
-				return condition(new detail::key_equals_number_or_empty_condition_impl(ci));
-		}
-
 		return condition(new detail::or_condition_impl(std::move(a), std::move(b)));
 	}
 
@@ -1066,20 +946,10 @@ struct key
 template <typename T>
 concept Numeric = ((std::is_floating_point_v<T> or std::is_integral_v<T>) and not std::is_same_v<T, bool>);
 
-/**
- * @brief Operator to create an equals condition based on a key @a key and a numeric value @a v
- */
-template <Numeric T>
-condition operator==(const key &key, const T &v)
-{
-	// TODO: change key_equals_etc... to use std::variant<double,int64_t> or something
-	return condition(new detail::key_equals_number_condition_impl(key.m_item_name, static_cast<double>(v)));
-}
-
 /**
  * @brief Operator to create an equals condition based on a key @a key and a value @a value
  */
-inline condition operator==(const key &key, std::string_view value)
+inline condition operator==(const key &key, item_value value)
 {
 	if (not value.empty())
 		return condition(new detail::key_equals_condition_impl({ key.m_item_name, value }));
@@ -1087,16 +957,6 @@ inline condition operator==(const key &key, std::string_view value)
 		return condition(new detail::key_is_empty_condition_impl(key.m_item_name));
 }
 
-/**
- * @brief Operator to create an equals condition based on a key @a key and a value @a value
- */
-template <typename T>
-	requires std::is_same_v<T, bool>
-inline condition operator==(const key &key, T value)
-{
-	return condition(new detail::key_equals_condition_impl({ key.m_item_name, value ? "y" : "n" }));
-}
-
 /**
  * @brief Operator to create a not equals condition based on a key @a key and a value @a v
  */

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++/iterator.hpp

@@ -26,6 +26,7 @@
 
 #pragma once
 
+#include "cif++/condition.hpp"
 #include "cif++/row.hpp"
 
 #include <array>
@@ -179,7 +180,7 @@ class iterator_impl
 	template <std::size_t... Is>
 	tuple_type get(std::index_sequence<Is...>) const
 	{
-		return m_current ? tuple_type{ m_current[m_item_ix[Is]].template as<Ts>()... } : tuple_type{};
+		return m_current ? tuple_type{ m_current[m_item_ix[Is]].template get<Ts>()... } : tuple_type{};
 	}
 
 	row_handle m_current;
@@ -422,7 +423,7 @@ class iterator_impl<Category, T>
   private:
 	value_type get() const
 	{
-		return m_current ?  m_current[m_item_ix].template as<value_type>() : value_type{};
+		return m_current ?  m_current[m_item_ix].template get<value_type>() : value_type{};
 	}
 
 	row_handle m_current;

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;
@@ -183,7 +190,7 @@ class atom
 	 * @param row The row containing the data for this atom
 	 */
 	atom(const datablock &db, const row_handle &row)
-		: atom(std::make_shared<atom_impl>(db, row["id"].as<std::string>()))
+		: atom(std::make_shared<atom_impl>(db, row["id"].get<std::string>()))
 	{
 	}
 
@@ -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++/parser.hpp

@@ -164,7 +164,13 @@ class sac_parser
 		SAVE_NAME,
 		STOP,
 		ITEM_NAME,
-		VALUE
+
+		VALUE_INAPPLICABLE,
+		VALUE_UNKNOWN,
+		VALUE_NUMERIC_INTEGER,
+		VALUE_NUMERIC_FLOAT,
+		VALUE_CHARSTRING,
+		VALUE_TEXTFIELD
 	};
 
 	static constexpr const char *get_token_name(CIFToken token)
@@ -180,7 +186,15 @@ class sac_parser
 			case CIFToken::SAVE_NAME: return "SAVE+name";
 			case CIFToken::STOP: return "STOP";
 			case CIFToken::ITEM_NAME: return "Tag";
-			case CIFToken::VALUE: return "Value";
+			// case CIFToken::VALUE: return "Value";
+
+			case CIFToken::VALUE_INAPPLICABLE: return "Inapplicable value";
+			case CIFToken::VALUE_UNKNOWN: return "'Unknown' value (=null)";
+			case CIFToken::VALUE_NUMERIC_INTEGER: return "Integer value";
+			case CIFToken::VALUE_NUMERIC_FLOAT: return "Float value";
+			case CIFToken::VALUE_CHARSTRING: return "Charstring value";
+			case CIFToken::VALUE_TEXTFIELD: return "Textfield value";
+
 			default: return "Invalid token parameter";
 		}
 	}
@@ -262,7 +276,7 @@ class sac_parser
 	virtual void produce_datablock(std::string_view name) = 0;
 	virtual void produce_category(std::string_view name) = 0;
 	virtual void produce_row() = 0;
-	virtual void produce_item(std::string_view category, std::string_view item, std::string_view value) = 0;
+	virtual void produce_item(std::string_view category, std::string_view item, item_value value) = 0;
 
   protected:
 
@@ -281,7 +295,12 @@ class sac_parser
 		TextItem,
 		TextItemNL,
 		Reserved,
-		Value
+		Value,
+
+		Numeric_Integer,
+		Numeric_Float,
+		Numeric_Exponent1,
+		Numeric_Exponent2
 	};
 
 	std::streambuf &m_source;
@@ -294,6 +313,8 @@ class sac_parser
 	// token buffer
 	std::vector<char> m_token_buffer;
 	std::string_view m_token_value;
+	int64_t m_token_value_int;
+	double m_token_value_float;
 
 	/** @endcond */
 };
@@ -331,7 +352,7 @@ class parser : public sac_parser
 
 	void produce_row() override;
 
-	void produce_item(std::string_view category, std::string_view item, std::string_view value) override;
+	void produce_item(std::string_view category, std::string_view item, item_value value) override;
 
   protected:
 	file &m_file;

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

include/cif++/row.hpp

@@ -41,7 +41,7 @@
  * that's not the case.
  * 
  * You can access the values of stored items by name or index.
- * The return value of operator[] is an cif::item_handle object.
+ * The return value of operator[] is a reference to a cif::item_value object.
  * 
  * @code {.cpp}
  * cif::category &atom_site = my_db["atom_site"];
@@ -93,7 +93,7 @@ namespace detail
 		{
 		}
 
-		const item_handle operator[](uint16_t ix) const
+		const item_value &operator[](uint16_t ix) const
 		{
 			return m_row[m_items[ix]];
 		}
@@ -107,7 +107,7 @@ namespace detail
 		template <typename... Ts, std::size_t... Is>
 		std::tuple<Ts...> get(std::index_sequence<Is...>) const
 		{
-			return std::tuple<Ts...>{ m_row[m_items[Is]].template as<Ts>()... };
+			return std::tuple<Ts...>{ m_row[m_items[Is]].template get<Ts>()... };
 		}
 
 		const row_handle &m_row;
@@ -173,14 +173,14 @@ class row : public std::vector<item_value>
 		return ix < size() ? &data()[ix] : nullptr;
 	}
 
-  private:
+//   private:
 	friend class category;
 	friend class category_index;
 
 	template <typename, typename...>
 	friend class iterator_impl;
 
-	void append(uint16_t ix, item_value &&iv)
+	void append(uint16_t ix, item_value iv)
 	{
 		if (ix >= size())
 			resize(ix + 1);
@@ -204,7 +204,6 @@ class row_handle
 {
   public:
 	/** @cond */
-	friend struct item_handle;
 	friend class category;
 	friend class category_index;
 	friend class row_initializer;
@@ -245,29 +244,20 @@ class row_handle
 		return not empty();
 	}
 
-	/// \brief return a cif::item_handle to the item in item @a item_ix
-	item_handle operator[](uint16_t item_ix)
-	{
-		return empty() ? item_handle::s_null_item : item_handle(item_ix, *this);
-	}
+	/// \brief return the count of the items
+	[[nodiscard]] size_t size() const { return m_row->size(); }
 
-	/// \brief return a const cif::item_handle to the item in item @a item_ix
-	const item_handle operator[](uint16_t item_ix) const
-	{
-		return empty() ? item_handle::s_null_item : item_handle(item_ix, const_cast<row_handle &>(*this));
-	}
+	/// \brief return a reference to a cif::item_value to the item in item @a item_ix
+	item_value &operator[](uint16_t item_ix);
 
-	/// \brief return a cif::item_handle to the item in the item named @a item_name
-	item_handle operator[](std::string_view item_name)
-	{
-		return empty() ? item_handle::s_null_item : item_handle(add_item(item_name), *this);
-	}
+	/// \brief return a const reference to a cif::item_value to the item in item @a item_ix
+	const item_value &operator[](uint16_t item_ix) const;
 
-	/// \brief return a const cif::item_handle to the item in the item named @a item_name
-	const item_handle operator[](std::string_view item_name) const
-	{
-		return empty() ? item_handle::s_null_item : item_handle(get_item_ix(item_name), const_cast<row_handle &>(*this));
-	}
+	/// \brief return a reference to a cif::item_value to the item in the item named @a item_name
+	item_value &operator[](std::string_view item_name);
+
+	/// \brief return a const reference to a cif::item_value to the item in the item named @a item_name
+	const item_value &operator[](std::string_view item_name) const;
 
 	/// \brief Return an object that can be used in combination with cif::tie
 	/// to assign the values for the items @a items
@@ -288,14 +278,14 @@ class row_handle
 	template <typename T>
 	T get(const char *item) const
 	{
-		return operator[](get_item_ix(item)).template as<T>();
+		return operator[](get_item_ix(item)).template get<T>();
 	}
 
 	/// \brief Get the value of item @a item cast to type @a T
 	template <typename T>
 	T get(std::string_view item) const
 	{
-		return operator[](get_item_ix(item)).template as<T>();
+		return operator[](get_item_ix(item)).template get<T>();
 	}
 
 	/// \brief assign each of the items named in @a values to their respective value
@@ -316,9 +306,9 @@ class row_handle
 	 * checked to see if it conforms to the rules defined in the dictionary
 	 */
 
-	void assign(std::string_view name, std::string_view value, bool updateLinked, bool validate = true)
+	void assign(std::string_view name, item_value value, bool updateLinked, bool validate = true)
 	{
-		assign(add_item(name), value, updateLinked, validate);
+		assign(add_item(name), std::move(value), updateLinked, validate);
 	}
 
 	/** \brief assign the value @a value to item at index @a item
@@ -332,7 +322,7 @@ class row_handle
 	 * checked to see if it conforms to the rules defined in the dictionary
 	 */
 
-	void assign(uint16_t item, std::string_view value, bool updateLinked, bool validate = true);
+	void assign(uint16_t item, item_value value, bool updateLinked, bool validate = true);
 
 	/// \brief compare two rows
 	bool operator==(const row_handle &rhs) const { return m_category == rhs.m_category and m_row == rhs.m_row; }
@@ -356,10 +346,10 @@ class row_handle
 		return m_row;
 	}
 
-	void assign(const item &i, bool updateLinked)
+	void assign(const item &i, bool updateLinked);/* 
 	{
 		assign(i.name(), i.value(), updateLinked);
-	}
+	} */
 
 	void swap(uint16_t item, row_handle &r);
 
@@ -408,7 +398,7 @@ class row_initializer : public std::vector<item>
 
 
 	/// \brief set the value for item name @a name to @a value
-	void set_value(std::string_view name, std::string_view value);
+	void set_value(std::string name, item_value value);
 
 	/// \brief set the value for item based on @a i
 	void set_value(const item &i)
@@ -417,7 +407,7 @@ class row_initializer : public std::vector<item>
 	}
 
 	/// \brief set the value for item name @a name to @a value, but only if the item did not have a value already
-	void set_value_if_empty(std::string_view name, std::string_view value);
+	void set_value_if_empty(std::string name, item_value value);
 
 	/// \brief set the value for item @a i, but only if the item did not have a value already
 	void set_value_if_empty(const item &i)

include/cif++/validate.hpp

@@ -329,10 +329,10 @@ struct item_validator
 
 	/// @brief Validate the value in @a value for this item
 	/// Will throw a std::system_error exception if it fails
-	void operator()(std::string_view value) const;
+	void operator()(const item_value &value) const;
 
 	/// @brief A more gentle version of value validation
-	bool validate_value(std::string_view value, std::error_code &ec) const noexcept;
+	bool validate_value(const item_value &value, std::error_code &ec) const noexcept;
 };
 
 /**

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/condition.cpp

@@ -123,21 +123,6 @@ namespace detail
 		return this;
 	}
 
-	condition_impl *key_equals_number_condition_impl::prepare(const category &c)
-	{
-		m_item_ix = c.get_item_ix(m_item_name);
-
-		if (c.get_cat_validator() != nullptr and
-			c.key_item_indices().contains(m_item_ix) and
-			c.key_item_indices().size() == 1)
-		{
-			item v(m_item_name, m_value);
-			m_single_hit = c[{ { m_item_name, std::string{ v.value() }, false } }];
-		}
-
-		return this;
-	}
-
 	bool found_in_range(condition_impl *c, std::vector<and_condition_impl *>::iterator b, std::vector<and_condition_impl *>::iterator e)
 	{
 		bool result = true;
@@ -234,17 +219,6 @@ namespace detail
 					continue;
 				}
 
-				if (auto s = dynamic_cast<const key_equals_number_condition_impl *>(sub); s != nullptr)
-				{
-					if (keys.contains(s->m_item_name))
-					{
-						item v{ s->m_item_name, s->m_value };
-						lookup.emplace_back(s->m_item_name, std::string{ v.value() } );
-						subs.emplace_back(sub);
-					}
-					continue;
-				}
-
 				if (auto s = dynamic_cast<const key_equals_or_empty_condition_impl *>(sub); s != nullptr)
 				{
 					if (keys.contains(s->m_item_name))
@@ -255,17 +229,6 @@ namespace detail
 					}
 					continue;
 				}
-
-				if (auto s = dynamic_cast<const key_equals_number_or_empty_condition_impl *>(sub); s != nullptr)
-				{
-					if (keys.contains(s->m_item_name))
-					{
-						item v{ s->m_item_name, s->m_value };
-						lookup.emplace_back(s->m_item_name, std::string{ v.value() }, true );
-						subs.emplace_back(sub);
-					}
-					continue;
-				}
 			}
 
 			if (lookup.size() == keys.size())

src/dictionary_parser.cpp

@@ -153,7 +153,7 @@ class dictionary_parser : public parser
 					match(CIFToken::ITEM_NAME);
 				}
 
-				while (m_lookahead == CIFToken::VALUE)
+				while (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
 				{
 					cat->emplace({});
 					auto row = cat->back();
@@ -161,7 +161,7 @@ class dictionary_parser : public parser
 					for (auto item_name : item_names)
 					{
 						row[item_name] = m_token_value;
-						match(CIFToken::VALUE);
+						match(m_lookahead);
 					}
 				}
 
@@ -181,7 +181,7 @@ class dictionary_parser : public parser
 					cat->emplace({});
 				cat->back()[item_name] = m_token_value;
 
-				match(CIFToken::VALUE);
+				match(m_lookahead >= CIFToken::VALUE_INAPPLICABLE ? m_lookahead : CIFToken::VALUE_CHARSTRING);
 			}
 		}
 
@@ -193,11 +193,11 @@ class dictionary_parser : public parser
 
 			std::vector<std::string> keys;
 			for (auto k : dict["category_key"])
-				keys.push_back(std::get<1>(split_item_name(k["name"].as<std::string>())));
+				keys.push_back(std::get<1>(split_item_name(k["name"].get<std::string>())));
 
 			iset groups;
 			for (auto g : dict["category_group"])
-				groups.insert(g["id"].as<std::string>());
+				groups.insert(g["id"].get<std::string>());
 
 			mCategoryValidators.push_back(category_validator{ category, keys, groups });
 		}
@@ -212,7 +212,7 @@ class dictionary_parser : public parser
 
 			iset ess;
 			for (auto e : dict["item_enumeration"])
-				ess.insert(e["value"].as<std::string>());
+				ess.insert(e["value"].get<std::string>());
 
 			std::string defaultValue = dict["item_default"].front().get<std::string>("value");
 			// bool defaultIsNull = false;
@@ -405,7 +405,7 @@ class dictionary_parser : public parser
 			// look up the label
 			for (auto r : linkedGroup.find("category_id"_key == link.m_child_category and "link_group_id"_key == link.m_link_group_id))
 			{
-				link.m_link_group_label = r["label"].as<std::string>();
+				link.m_link_group_label = r["label"].get<std::string>();
 				break;
 			}
 

src/item.cpp

@@ -24,45 +24,109 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/item.hpp"
+
 #include "cif++/row.hpp"
 
 #include <cassert>
+#include <compare>
+#include <ios>
 
 namespace cif
 {
 
-const item_handle item_handle::s_null_item;
-row_handle s_null_row_handle;
-
-item_handle::item_handle()
-	: m_item_ix(std::numeric_limits<uint16_t>::max())
-	, m_row_handle(s_null_row_handle)
+int item_value::compare(const item_value &b, bool ignore_case) const noexcept
 {
-}
+	int d = static_cast<int>(m_data.m_type) - static_cast<int>(b.m_data.m_type);
 
-std::string_view item_handle::text() const
-{
-	if (not m_row_handle.empty())
+	if (d == 0)
 	{
-		auto iv = m_row_handle.m_row->get(m_item_ix);
-		if (iv != nullptr)
-			return iv->text();
+		switch (m_data.m_type)
+		{
+			case cif::item_value_type::BOOLEAN:
+				d = static_cast<int>(m_data.m_value.m_boolean) - static_cast<int>(b.m_data.m_value.m_boolean);
+				break;
+			case cif::item_value_type::INT:
+				d = m_data.m_value.m_integer - b.m_data.m_value.m_integer;
+				break;
+			case cif::item_value_type::FLOAT:
+			{
+				auto dp = (m_data.m_value.m_float <=> b.m_data.m_value.m_float);
+				if (dp == std::partial_ordering::less)
+					d = -1;
+				else if (dp == std::partial_ordering::greater)
+					d = 1;
+				break;
+			}
+			case cif::item_value_type::TEXT:
+				d = m_data.sv().compare(b.m_data.sv());
+				break;
+			default:;
+		}
 	}
 
-	return {};
+	return d;
 }
 
-void item_handle::assign_value(std::string_view value)
+// const item_handle item_handle::s_null_item;
+// row_handle s_null_row_handle;
+
+// item_handle::item_handle()
+// 	: m_item_ix(std::numeric_limits<uint16_t>::max())
+// 	, m_row_handle(s_null_row_handle)
+// {
+// }
+
+// std::string_view item_handle::text() const
+// {
+// 	if (not m_row_handle.empty())
+// 	{
+// 		auto iv = m_row_handle.m_row->get(m_item_ix);
+// 		if (iv != nullptr)
+// 			return iv->text();
+// 	}
+
+// 	return {};
+// }
+
+// void item_handle::assign_value(std::string_view value)
+// {
+// 	assert(not m_row_handle.empty());
+// 	m_row_handle.assign(m_item_ix, value, true);
+// }
+
+// void item_handle::swap(item_handle &b)
+// {
+// 	assert(m_item_ix == b.m_item_ix);
+// 	// assert(&m_row_handle.m_category == &b.m_row_handle.m_category);
+// 	m_row_handle.swap(m_item_ix, b.m_row_handle);
+// }
+
+std::ostream &operator<<(std::ostream &os, const item_value &v)
 {
-	assert(not m_row_handle.empty());
-	m_row_handle.assign(m_item_ix, value, true);
+	switch (v.type())
+	{
+		case cif::item_value_type::BOOLEAN:
+			os << std::boolalpha << v.m_data.m_value.m_boolean;
+			break;
+		case cif::item_value_type::INT:
+			os << v.m_data.m_value.m_integer;
+			break;
+		case cif::item_value_type::FLOAT:
+			os << v.m_data.m_value.m_float;
+			break;
+		case cif::item_value_type::TEXT:
+			os << v.m_data.sv();
+			break;
+		case cif::item_value_type::MISSING:
+			os << '?';
+			break;
+		case cif::item_value_type::EMPTY:
+			os << '.';
+			break;
+	}
+
+	return os;
 }
 
-void item_handle::swap(item_handle &b)
-{
-	assert(m_item_ix == b.m_item_ix);
-	// assert(&m_row_handle.m_category == &b.m_row_handle.m_category);
-	m_row_handle.swap(m_item_ix, b.m_row_handle);
-}
-
-}
+} // namespace cif

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;
 
@@ -59,7 +62,7 @@ void atom::atom_impl::moveTo(const point &p)
 
 std::string atom::atom_impl::get_property(std::string_view name) const
 {
-	return row()[name].as<std::string>();
+	return row()[name].get<std::string>();
 }
 
 int atom::atom_impl::get_property_int(std::string_view name) const
@@ -132,7 +135,7 @@ int atom::atom_impl::compare(const atom_impl &b) const
 
 int atom::atom_impl::get_charge() const
 {
-	auto formalCharge = row()["pdbx_formal_charge"].as<std::optional<int>>();
+	auto formalCharge = row()["pdbx_formal_charge"].get<std::optional<int>>();
 
 	if (not formalCharge.has_value())
 	{
@@ -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
@@ -1898,11 +1891,7 @@ void structure::swap_atoms(atom a1, atom a2)
 		auto r2 = atomSites.find1(key("id") == a2.id());
 
 		for (std::string fld : std::initializer_list<std::string>{ "label_atom_id", "auth_atom_id", "type_symbol" })
-		{
-			auto l1 = r1[fld];
-			auto l2 = r2[fld];
-			l1.swap(l2);
-		}
+			swap(r1[fld], r2[fld]);
 	}
 	catch (const std::exception &ex)
 	{
@@ -2346,6 +2335,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;
@@ -2735,7 +2754,7 @@ void structure::cleanup_empty_categories()
 
 	for (auto chemComp : chem_comp)
 	{
-		std::string compID = chemComp["id"].as<std::string>();
+		std::string compID = chemComp["id"].get<std::string>();
 		if (atomSite.contains("label_comp_id"_key == compID or "auth_comp_id"_key == compID) or
 			pdbxPolySeqScheme.contains("mon_id"_key == compID or "auth_mon_id"_key == compID or "pdb_mon_id"_key == compID) or
 			entityPolySeq.contains("mon_id"_key == compID))
@@ -2756,7 +2775,7 @@ void structure::cleanup_empty_categories()
 
 	for (auto entity : entities)
 	{
-		std::string entityID = entity["id"].as<std::string>();
+		std::string entityID = entity["id"].get<std::string>();
 		if (atomSite.contains("label_entity_id"_key == entityID))
 			continue;
 

src/parser.cpp

@@ -24,15 +24,19 @@
  * 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++/forward_decl.hpp"
+#include "cif++/item.hpp"
+#include "cif++/utilities.hpp"
 
 #include <cassert>
+#include <charconv>
 #include <iostream>
 #include <map>
 #include <stack>
+#include <system_error>
 
 namespace cif
 {
@@ -58,12 +62,12 @@ class reserved_words_automaton
 
 	constexpr bool finished() const
 	{
-		return m_state <= 0; 
+		return m_state <= 0;
 	}
 
 	constexpr bool matched() const
 	{
-		return m_state < 0; 
+		return m_state < 0;
 	}
 
 	constexpr move_result move(int ch)
@@ -75,7 +79,7 @@ class reserved_words_automaton
 			case 0:
 				break;
 
-			case -1:		// data_
+			case -1: // data_
 				if (sac_parser::is_non_blank(ch))
 					m_seen_trailing_chars = true;
 				else if (m_seen_trailing_chars)
@@ -84,15 +88,15 @@ class reserved_words_automaton
 					result = no_keyword;
 				break;
 
-			case -2:		// global_
+			case -2: // global_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : global;
 				break;
 
-			case -3:		// loop_
+			case -3: // loop_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : loop;
 				break;
 
-			case -4:		// save_
+			case -4: // save_
 				if (sac_parser::is_non_blank(ch))
 					m_seen_trailing_chars = true;
 				else if (m_seen_trailing_chars)
@@ -101,10 +105,10 @@ class reserved_words_automaton
 					result = save;
 				break;
 
-			case -5:		// stop_
+			case -5: // stop_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : stop;
 				break;
-			
+
 			default:
 				assert(m_state > 0 and m_state < NODE_COUNT);
 
@@ -141,13 +145,13 @@ class reserved_words_automaton
 		int8_t next_nomatch;
 	} s_dag[] = {
 		{ 0 },
-		{ 'D',  5, 2 },
-		{ 'G',  9, 3 },
+		{ 'D', 5, 2 },
+		{ 'G', 9, 3 },
 		{ 'L', 15, 4 },
 		{ 'S', 19, 0 },
-		{ 'A',  6, 0 },
-		{ 'T',  7, 0 },
-		{ 'A',  8, 0 },
+		{ 'A', 6, 0 },
+		{ 'T', 7, 0 },
+		{ 'A', 8, 0 },
 		{ '_', -1, 0 },
 		{ 'L', 10, 0 },
 		{ 'O', 11, 0 },
@@ -155,7 +159,7 @@ class reserved_words_automaton
 		{ 'A', 13, 0 },
 		{ 'L', 14, 0 },
 		{ '_', -2, 0 },
-		{ 'O', 16, 0},
+		{ 'O', 16, 0 },
 		{ 'O', 17, 0 },
 		{ 'P', 18, 0 },
 		{ '_', -3, 0 },
@@ -238,7 +242,7 @@ int sac_parser::get_next_char()
 		}
 		else if (result == '\n')
 			++m_line_nr;
-		
+
 		m_token_buffer.push_back(std::char_traits<char>::to_char_type(result));
 	}
 
@@ -277,6 +281,8 @@ sac_parser::CIFToken sac_parser::get_next_token()
 	m_token_buffer.clear();
 	m_token_value = {};
 
+	bool negative = false;
+
 	reserved_words_automaton dag;
 
 	while (result == CIFToken::UNKNOWN)
@@ -310,6 +316,15 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				}
 				else if (dag.move(ch) == reserved_words_automaton::undefined)
 					state = State::Reserved;
+				else if (ch == '+' or ch == '-')
+				{
+					negative = true;
+					state = State::Numeric_Integer;
+				}
+				else if (ch >= '0' and ch <= '9')
+					state = State::Numeric_Integer;
+				else if (ch == '.')
+					state = State::Numeric_Float;
 				else
 					state = State::Value;
 				break;
@@ -326,7 +341,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else
 					m_bol = (ch == '\n');
 				break;
-			
+
 			case State::Comment:
 				if (ch == '\n')
 				{
@@ -339,12 +354,12 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (not is_any_print(ch))
 					error("invalid character in comment");
 				break;
-			
+
 			case State::QuestionMark:
 				if (not is_non_blank(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_UNKNOWN;
 				}
 				else
 					state = State::Value;
@@ -356,7 +371,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (ch == kEOF)
 					error("unterminated textfield");
 				else if (not is_any_print(ch) and cif::VERBOSE > 2)
-					warning("invalid character in text field '" + std::string({static_cast<char>(ch)}) + "' (" + std::to_string((int)ch) + ")");
+					warning("invalid character in text field '" + std::string({ static_cast<char>(ch) }) + "' (" + std::to_string((int)ch) + ")");
 				break;
 
 			case State::TextItemNL:
@@ -366,7 +381,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				{
 					assert(m_token_buffer.size() >= 2);
 					m_token_value = std::string_view(m_token_buffer.data() + 1, m_token_buffer.size() - 3);
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_TEXTFIELD;
 				}
 				else if (ch == kEOF)
 					error("unterminated textfield");
@@ -380,14 +395,14 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (ch == quoteChar)
 					state = State::QuotedStringQuote;
 				else if (not is_any_print(ch) and cif::VERBOSE > 2)
-					warning("invalid character in quoted string: '" + std::string({static_cast<char>(ch)}) + "' (" + std::to_string((int)ch) + ")");
+					warning("invalid character in quoted string: '" + std::string({ static_cast<char>(ch) }) + "' (" + std::to_string((int)ch) + ")");
 				break;
 
 			case State::QuotedStringQuote:
 				if (is_white(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_CHARSTRING;
 					if (m_token_buffer.size() < 2)
 						error("Invalid quoted string token");
 
@@ -422,7 +437,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 						if (not is_non_blank(ch))
 						{
 							retract();
-							result = CIFToken::VALUE;
+							result = CIFToken::VALUE_CHARSTRING;
 							m_token_value = std::string_view(m_token_buffer.data(), m_token_buffer.size());
 						}
 						else
@@ -463,11 +478,65 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				}
 				break;
 
+			case State::Numeric_Integer:
+				if (ch == '.')
+					state = State::Numeric_Float;
+				else if (ch == 'e' or ch == 'E')
+					state = State::Numeric_Exponent1;
+				else if (not is_non_blank(ch))
+				{
+					retract();
+					result = CIFToken::VALUE_NUMERIC_INTEGER;
+				}
+				else if (ch < '0' or ch > '9')
+					state = State::Value;
+				break;
+
+			case State::Numeric_Float:
+				if (not is_non_blank(ch))
+				{
+					retract();
+					if (m_token_buffer.size() == 1)
+						result = CIFToken::VALUE_INAPPLICABLE;
+					else
+						result = CIFToken::VALUE_NUMERIC_FLOAT;
+				}
+				else if (ch == 'e' or ch == 'E')
+					state = State::Numeric_Exponent1;
+				else if (ch < '0' or ch > '9')
+					state = State::Value;
+				break;
+
+			case State::Numeric_Exponent1:
+				if (ch == '+' or ch == '-' or (ch >= '0' and ch <= '9'))
+					state = State::Numeric_Exponent2;
+				else
+				{
+					if (VERBOSE > 0)
+						std::cerr << "parsing " << std::string_view{ m_token_buffer.data(), m_token_buffer.size() } << " Invalid floating point value, expected digit or sign character\n";
+					state = State::Value;
+				}
+				break;
+
+			case State::Numeric_Exponent2:
+				if (not is_non_blank(ch))
+				{
+					retract();
+					result = CIFToken::VALUE_NUMERIC_FLOAT;
+				}
+				else if (ch < '0' or ch > '9')
+				{
+					if (VERBOSE > 0)
+						std::cerr << "parsing " << std::string_view{ m_token_buffer.data(), m_token_buffer.size() } << " Invalid floating point value, expected exponent digit\n";
+					state = State::Value;
+				}
+				break;
+
 			case State::Value:
 				if (not is_non_blank(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_CHARSTRING;
 					m_token_value = std::string_view(m_token_buffer.data(), m_token_buffer.size());
 					break;
 				}
@@ -480,12 +549,25 @@ sac_parser::CIFToken sac_parser::get_next_token()
 		}
 	}
 
-	if (VERBOSE >= 5)
+	// if (VERBOSE >= 5)
+	// {
+	// 	std::cerr << get_token_name(result);
+	// 	if (result != CIFToken::END_OF_FILE)
+	// 		std::cerr << " " << std::quoted(m_token_value);
+	// 	std::cerr << '\n';
+	// }
+
+	if (result == CIFToken::VALUE_NUMERIC_INTEGER)
 	{
-		std::cerr << get_token_name(result);
-		if (result != CIFToken::END_OF_FILE)
-			std::cerr << " " << std::quoted(m_token_value);
-		std::cerr << '\n';
+		auto [ptr, ec] = std::from_chars(m_token_buffer.data(), m_token_buffer.data() + m_token_buffer.size(), m_token_value_int);
+		if (ec != std::errc{})
+			error("Invalid integer value: " + std::make_error_code(ec).message());
+	}
+	else if (result == CIFToken::VALUE_NUMERIC_FLOAT)
+	{
+		auto [ptr, ec] = std::from_chars(m_token_buffer.data(), m_token_buffer.data() + m_token_buffer.size(), m_token_value_float);
+		if (ec != std::errc{})
+			error("Invalid integer value: " + std::make_error_code(ec).message());
 	}
 
 	return result;
@@ -661,7 +743,7 @@ sac_parser::datablock_index sac_parser::index_datablocks()
 			case data:
 				if (dblk[si] == 0 and is_non_blank(ch))
 				{
-					datablock = {static_cast<char>(ch)};
+					datablock = { static_cast<char>(ch) };
 					state = data_name;
 				}
 				else if (dblk[si++] != ch)
@@ -738,14 +820,17 @@ void sac_parser::parse_global()
 	while (m_lookahead == CIFToken::ITEM_NAME)
 	{
 		match(CIFToken::ITEM_NAME);
-		match(CIFToken::VALUE);
+		if (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
+			match(m_lookahead);
+		else
+			match(CIFToken::VALUE_CHARSTRING);
 	}
 }
 
 void sac_parser::parse_datablock()
 {
 	static const std::string kUnitializedCategory("<invalid>");
-	std::string cat = kUnitializedCategory;	// intial value acts as a guard for empty category names
+	std::string cat = kUnitializedCategory; // intial value acts as a guard for empty category names
 
 	while (m_lookahead == CIFToken::LOOP or m_lookahead == CIFToken::ITEM_NAME or m_lookahead == CIFToken::SAVE_NAME)
 	{
@@ -777,14 +862,38 @@ void sac_parser::parse_datablock()
 					match(CIFToken::ITEM_NAME);
 				}
 
-				while (m_lookahead == CIFToken::VALUE)
+				while (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
 				{
 					produce_row();
 
 					for (auto item_name : item_names)
 					{
-						produce_item(cat, item_name, m_token_value);
-						match(CIFToken::VALUE);
+						switch (m_lookahead)
+						{
+							case CIFToken::VALUE_INAPPLICABLE:
+								produce_item(cat, item_name, nullptr);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_UNKNOWN:
+								produce_item(cat, item_name, std::optional<std::string>{});
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_NUMERIC_INTEGER:
+								produce_item(cat, item_name, m_token_value_int);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_NUMERIC_FLOAT:
+								produce_item(cat, item_name, m_token_value_float);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_CHARSTRING:
+							case CIFToken::VALUE_TEXTFIELD:
+								produce_item(cat, item_name, m_token_value);
+								match(m_lookahead);
+								break;
+							default:;
+								match(CIFToken::VALUE_CHARSTRING);
+						}
 					}
 				}
 
@@ -806,9 +915,33 @@ void sac_parser::parse_datablock()
 
 				match(CIFToken::ITEM_NAME);
 
-				produce_item(cat, itemName, m_token_value);
+				switch (m_lookahead)
+				{
+					case CIFToken::VALUE_INAPPLICABLE:
+						produce_item(cat, itemName, nullptr);
+						match(CIFToken::VALUE_INAPPLICABLE);
+						break;
+					case CIFToken::VALUE_UNKNOWN:
+						produce_item(cat, itemName, item_value{ std::optional<std::string>{} });
+						match(CIFToken::VALUE_UNKNOWN);
+						break;
+					case CIFToken::VALUE_NUMERIC_INTEGER:
+						produce_item(cat, itemName, m_token_value_int);
+						match(CIFToken::VALUE_NUMERIC_INTEGER);
+						break;
+					case CIFToken::VALUE_NUMERIC_FLOAT:
+						produce_item(cat, itemName, m_token_value_float);
+						match(CIFToken::VALUE_NUMERIC_FLOAT);
+						break;
+					case CIFToken::VALUE_CHARSTRING:
+					case CIFToken::VALUE_TEXTFIELD:
+						produce_item(cat, itemName, m_token_value);
+						match(m_lookahead);
+						break;
+					default:
+						match(CIFToken::VALUE_CHARSTRING);
+				}
 
-				match(CIFToken::VALUE);
 				break;
 			}
 
@@ -864,7 +997,7 @@ void parser::produce_row()
 	// m_row.lineNr(m_line_nr);
 }
 
-void parser::produce_item(std::string_view category, std::string_view item, std::string_view value)
+void parser::produce_item(std::string_view category, std::string_view item, item_value value)
 {
 	if (VERBOSE >= 4)
 		std::cerr << "producing _" << category << '.' << item << " -> " << value << '\n';
@@ -872,7 +1005,7 @@ void parser::produce_item(std::string_view category, std::string_view item, std:
 	if (m_category == nullptr or not iequals(category, m_category->name()))
 		error("inconsistent categories in loop_");
 
-	m_row[item] = m_token_value;
+	m_row[item] = std::move(value);
 }
 
 } // namespace cif

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/row.cpp

@@ -24,17 +24,50 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/row.hpp"
+
 #include "cif++/category.hpp"
+#include "cif++/item.hpp"
+#include <stdexcept>
 
 namespace cif
 {
+// --------------------------------------------------------------------
 
-void row_handle::assign(uint16_t item, std::string_view value, bool updateLinked, bool validate)
+item_value s_null_item;
+
+item_value &row_handle::operator[](uint16_t item_ix)
+{
+	return empty() or item_ix >= m_row->size() ? s_null_item : m_row->operator[](item_ix);
+}
+
+const item_value &row_handle::operator[](uint16_t item_ix) const
+{
+	return empty() or item_ix >= m_row->size() ? s_null_item : m_row->operator[](item_ix);
+}
+
+item_value &row_handle::operator[](std::string_view item_name)
+{
+	auto ix = add_item(item_name);
+	if (ix >= size())
+		m_row->resize(ix + 1);
+	
+	return m_row->operator[](ix);
+}
+
+const item_value &row_handle::operator[](std::string_view item_name) const
+{
+	return operator[](get_item_ix(item_name));
+}
+
+// --------------------------------------------------------------------
+
+void row_handle::assign(uint16_t item, item_value value, bool updateLinked, bool validate)
 {
 	if (not m_category)
 		throw std::runtime_error("uninitialized row");
 
-	m_category->update_value(m_row, item, value, updateLinked, validate);
+	m_category->update_value(m_row, item, std::move(value), updateLinked, validate);
 }
 
 uint16_t row_handle::get_item_ix(std::string_view name) const
@@ -86,28 +119,29 @@ row_initializer::row_initializer(row_handle rh)
 		auto &i = r->operator[](ix);
 		if (not i)
 			continue;
-		emplace_back(cat.get_item_name(ix), i.text());
+		emplace_back(cat.get_item_name(ix), i);
 	}
 }
 
-void row_initializer::set_value(std::string_view name, std::string_view value)
+void row_initializer::set_value(std::string name, item_value value)
 {
 	for (auto &i : *this)
 	{
 		if (i.name() == name)
 		{
-			i.value(value);
+			i.value(std::move(value));
 			return;
 		}
 	}
 
-	emplace_back(name, value);
+	emplace_back(std::move(name), std::move(value));
 }
 
-void row_initializer::set_value_if_empty(std::string_view name, std::string_view value)
+void row_initializer::set_value_if_empty(std::string name, item_value value)
 {
-	if (find_if(begin(), end(), [name](auto &i) { return i.name() == name; }) == end())
-		emplace_back(name, value);
+	if (std::ranges::find_if(*this, [name](auto &i)
+			{ return i.name() == name; }) == end())
+		emplace_back(std::move(name), std::move(value));
 }
 
 } // 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)
@@ -259,24 +263,24 @@ int type_validator::compare(std::string_view a, std::string_view b) const
 
 // --------------------------------------------------------------------
 
-void item_validator::operator()(std::string_view value) const
+void item_validator::operator()(const item_value &value) const
 {
-	std::error_code ec;
-	if (not validate_value(value, ec))
-		throw std::system_error(ec, std::string{ value } + " does not match rx for " + m_item_name);
+	// std::error_code ec;
+	// if (not validate_value(value, ec))
+	// 	throw std::system_error(ec, std::string{ value } + " does not match rx for " + m_item_name);
 }
 
-bool item_validator::validate_value(std::string_view value, std::error_code &ec) const noexcept
+bool item_validator::validate_value(const item_value &value, std::error_code &ec) const noexcept
 {
 	ec.clear();
 
-	if (not value.empty() and value != "?" and value != ".")
-	{
-		if (m_type != nullptr and not m_type->m_rx->match(value))
-			ec = make_error_code(validation_error::value_does_not_match_rx);
-		else if (not m_enums.empty() and m_enums.count(std::string{ value }) == 0)
-			ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
-	}
+	// if (not value.empty() and value != "?" and value != ".")
+	// {
+	// 	if (m_type != nullptr and not m_type->m_rx->match(value))
+	// 		ec = make_error_code(validation_error::value_does_not_match_rx);
+	// 	else if (not m_enums.empty() and m_enums.count(std::string{ value }) == 0)
+	// 		ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
+	// }
 
 	return not(bool) ec;
 }

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)
 
@@ -41,14 +41,15 @@ list(
 	APPEND
 	CIFPP_tests
 	unit-v2
-	unit-3d
-	model
-	query
-	rename-compound
-	sugar
-	spinner
-	reconstruction
-	validate-pdbx
+	# unit-3d
+	# model
+	# query
+	# rename-compound
+	# sugar
+	# 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/test-main.cpp

@@ -2,7 +2,7 @@
 
 #include "test-main.hpp"
 
-#include <cif++.hpp>
+#include <cif++/utilities.hpp>
 
 std::filesystem::path gTestDir = std::filesystem::current_path();
 
@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
 	// initialize CCD location
 	cif::add_file_resource("components.cif", gTestDir / ".." / "rsrc" / "ccd-subset.cif");
 
-	cif::compound_factory::instance().push_dictionary(gTestDir / "HEM.cif");
+	// cif::compound_factory::instance().push_dictionary(gTestDir / "HEM.cif");
 
 	return session.run();
 }

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));
+	}
+}