Moving cif::Structure back in as model

CMakeLists.txt

@@ -190,6 +190,8 @@ set(project_sources
 	${PROJECT_SOURCE_DIR}/src/compound.cpp
 	${PROJECT_SOURCE_DIR}/src/point.cpp
 	${PROJECT_SOURCE_DIR}/src/symmetry.cpp
+
+	${PROJECT_SOURCE_DIR}/src/model.cpp
 )
 
 set(project_headers
@@ -212,6 +214,8 @@ set(project_headers
 	${PROJECT_SOURCE_DIR}/include/cif++/compound.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/point.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/symmetry.hpp
+
+	${PROJECT_SOURCE_DIR}/include/cif++/model.hpp
 )
 
 add_library(cifpp ${project_sources} ${project_headers} ${PROJECT_SOURCE_DIR}/src/symop_table_data.hpp)
@@ -366,7 +370,7 @@ if(ENABLE_TESTING)
 
 	find_package(Boost REQUIRED headers)
 
-	list(APPEND CIFPP_tests unit-v2 unit-3d)
+	list(APPEND CIFPP_tests unit-v2 unit-3d format model rename-compound sugar)
 
 	foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 		set(CIFPP_TEST "${CIFPP_TEST}-test")

include/cif++.hpp

@@ -29,6 +29,10 @@
 #include <cif++/utilities.hpp>
 #include <cif++/file.hpp>
 #include <cif++/parser.hpp>
+#include <cif++/format.hpp>
 
 #include <cif++/compound.hpp>
-#include <cif++/point.hpp>
\ No newline at end of file
+#include <cif++/point.hpp>
+#include <cif++/symmetry.hpp>
+
+#include <cif++/model.hpp>
\ No newline at end of file

include/cif++/category.hpp

@@ -595,6 +595,10 @@ class category
 
 	// --------------------------------------------------------------------
 
+	void swap_item(size_t column_ix, row_handle &a, row_handle &b);
+
+	// --------------------------------------------------------------------
+	
 	std::string m_name;
 	std::vector<item_column> m_columns;
 	const validator *m_validator = nullptr;

include/cif++/compound.hpp

@@ -111,7 +111,7 @@ class compound
 	const std::vector<compound_atom> &atoms() const { return m_atoms; }
 	const std::vector<compound_bond> &bonds() const { return m_bonds; }
 
-	compound_atom get_atom_by_id(const std::string &atom_id) const;
+	compound_atom get_atom_by_atom_id(const std::string &atom_id) const;
 
 	bool atoms_bonded(const std::string &atomId_1, const std::string &atomId_2) const;
 	// float atomBondValue(const std::string &atomId_1, const std::string &atomId_2) const;

include/cif++/datablock.hpp

@@ -88,6 +88,10 @@ class datablock : public std::list<category>
 		return os;
 	}
 
+	// --------------------------------------------------------------------
+	
+	bool operator==(const datablock &rhs) const;
+
   private:
 	std::string m_name;
 	const validator *m_validator = nullptr;

include/cif++/file.hpp

@@ -84,6 +84,8 @@ class file : public std::list<datablock>
 	void load_dictionary();
 	void load_dictionary(std::string_view name);
 
+	bool contains(std::string_view name) const;
+
 	datablock &operator[](std::string_view name);
 	const datablock &operator[](std::string_view name) const;
 

include/cif++/format.hpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2022 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#pragma once
+
+#include <string>
+
+/// \file format.hpp
+/// File containing a basic reimplementation of boost::format
+/// but then a bit more simplistic. Still this allowed me to move my code
+/// from using boost::format to something without external dependency easily.
+
+namespace cif
+{
+
+namespace detail
+{
+	template <typename T>
+	struct to_varg
+	{
+		using type = T;
+
+		to_varg(const T &v)
+			: m_value(v)
+		{
+		}
+
+		type operator*() { return m_value; }
+
+		T m_value;
+	};
+
+	// template <>
+	// struct to_varg<char>
+	// {
+	// 	using type = const char *;
+
+	// 	to_varg(const char &v)
+	// 		: m_value({ v })
+	// 	{
+	// 	}
+
+	// 	type operator*() { return m_value.c_str(); }
+
+	// 	std::string m_value;
+	// };
+
+	template <>
+	struct to_varg<const char *>
+	{
+		using type = const char *;
+
+		to_varg(const char *v)
+			: m_value(v)
+		{
+		}
+
+		type operator*() { return m_value.c_str(); }
+
+		std::string m_value;
+	};
+
+	template <>
+	struct to_varg<std::string>
+	{
+		using type = const char *;
+
+		to_varg(const std::string &v)
+			: m_value(v)
+		{
+		}
+
+		type operator*() { return m_value.c_str(); }
+
+		std::string m_value;
+	};
+
+} // namespace
+
+template <typename... Args>
+class format_plus_arg
+{
+  public:
+	using args_vector_type = std::tuple<detail::to_varg<Args>...>;
+	using vargs_vector_type = std::tuple<typename detail::to_varg<Args>::type...>;
+
+	format_plus_arg(const format_plus_arg &) = delete;
+	format_plus_arg &operator=(const format_plus_arg &) = delete;
+
+
+	format_plus_arg(std::string_view fmt, Args... args)
+		: m_fmt(fmt)
+		, m_args(std::forward<Args>(args)...)
+	{
+		auto ix = std::make_index_sequence<sizeof...(Args)>();
+		copy_vargs(ix);
+	}
+
+	std::string str()
+	{
+		char buffer[1024];
+		std::string::size_type r = std::apply(snprintf, std::tuple_cat(std::make_tuple(buffer, sizeof(buffer), m_fmt.c_str()), m_vargs));
+		return { buffer, r };
+	}
+
+	friend std::ostream &operator<<(std::ostream &os, const format_plus_arg &f)
+	{
+		char buffer[1024];
+		std::string::size_type r = std::apply(snprintf, std::tuple_cat(std::make_tuple(buffer, sizeof(buffer), f.m_fmt.c_str()), f.m_vargs));
+		os.write(buffer, r);
+		return os;
+	}
+
+  private:
+
+	template <size_t... I>
+	void copy_vargs(std::index_sequence<I...>)
+	{
+		((std::get<I>(m_vargs) = *std::get<I>(m_args)), ...);
+	}
+
+	std::string m_fmt;
+	args_vector_type m_args;
+	vargs_vector_type m_vargs;
+};
+
+template <typename... Args>
+constexpr auto format(std::string_view fmt, Args... args)
+{
+	return format_plus_arg(fmt, std::forward<Args>(args)...);
+}
+
+// --------------------------------------------------------------------
+/// A streambuf that fills out lines with spaces up until a specified width
+
+class fill_out_streambuf : public std::streambuf
+{
+  public:
+	using base_type = std::streambuf;
+	using int_type = base_type::int_type;
+	using char_type = base_type::char_type;
+	using traits_type = base_type::traits_type;
+
+	fill_out_streambuf(std::ostream &os, int width = 80)
+		: m_os(os)
+		, m_upstream(os.rdbuf())
+		, m_width(width)
+	{
+	}
+
+	~fill_out_streambuf()
+	{
+		m_os.rdbuf(m_upstream);
+	}
+
+	virtual int_type
+	overflow(int_type ic = traits_type::eof())
+	{
+		char ch = traits_type::to_char_type(ic);
+
+		int_type result = ic;
+
+		if (ch == '\n')
+		{
+			for (int i = m_column_count; result != traits_type::eof() and i < m_width; ++i)
+				result = m_upstream->sputc(' ');
+		}
+
+		if (result != traits_type::eof())
+			result = m_upstream->sputc(ch);
+
+		if (result != traits_type::eof())
+		{
+			if (ch == '\n')
+			{
+				m_column_count = 0;
+				++m_line_count;
+			}
+			else
+				++m_column_count;
+		}
+
+		return result;
+	}
+
+	std::streambuf *get_upstream() const { return m_upstream; }
+
+	int get_line_count() const { return m_line_count; }
+
+  private:
+	std::ostream &m_os;
+	std::streambuf *m_upstream;
+	int m_width;
+	int m_line_count = 0;
+	int m_column_count = 0;
+};
+
+} // namespace pdbx

include/cif++/item.hpp

@@ -300,6 +300,11 @@ struct item_handle
 
 	static const item_handle s_null_item;
 
+	friend void swap(item_handle a, item_handle b)
+	{
+		a.swap(b);
+	}
+
   private:
 	item_handle();
 
@@ -318,21 +323,24 @@ struct item_handle::item_value_as<T, std::enable_if_t<std::is_arithmetic_v<T> an
 
 	static value_type convert(const item_handle &ref)
 	{
-		auto txt = ref.text();
-
 		value_type result = {};
 
-		std::from_chars_result r = selected_charconv<value_type>::from_chars(txt.data(), txt.data() + txt.size(), result);
-
-		if (r.ec != std::errc())
+		if (not ref.empty())
 		{
-			result = {};
-			if (cif::VERBOSE)
+			auto txt = ref.text();
+
+			std::from_chars_result r = selected_charconv<value_type>::from_chars(txt.data(), txt.data() + txt.size(), result);
+
+			if (r.ec != std::errc())
 			{
-				if (r.ec == std::errc::invalid_argument)
-					std::cerr << "Attempt to convert " << std::quoted(txt) << " into a number" << std::endl;
-				else if (r.ec == std::errc::result_out_of_range)
-					std::cerr << "Conversion of " << std::quoted(txt) << " into a type that is too small" << std::endl;
+				result = {};
+				if (cif::VERBOSE)
+				{
+					if (r.ec == std::errc::invalid_argument)
+						std::cerr << "Attempt to convert " << std::quoted(txt) << " into a number" << std::endl;
+					else if (r.ec == std::errc::result_out_of_range)
+						std::cerr << "Conversion of " << std::quoted(txt) << " into a type that is too small" << std::endl;
+				}
 			}
 		}
 

include/cif++/model.hpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#pragma once
+
+#include <numeric>
+
+#if __cpp_lib_format
+#include <format>
+#endif
+
+/*
+    To modify a structure, you will have to use actions.
+
+    The currently supported actions are:
+
+//	- Move atom to new location
+    - Remove atom
+//	- Add new atom that was formerly missing
+//	- Add alternate residue
+    -
+
+*/
+
+namespace cif::mm
+{
+
+class atom;
+class residue;
+class monomer;
+class polymer;
+class structure;
+
+// --------------------------------------------------------------------
+
+class atom
+{
+  private:
+	struct atom_impl : public std::enable_shared_from_this<atom_impl>
+	{
+		atom_impl(datablock &db, std::string_view id, row_handle row)
+			: m_db(db)
+			, m_id(id)
+			, m_row(row)
+		{
+			prefetch();
+		}
+
+		// constructor for a symmetry copy of an atom
+		atom_impl(const atom_impl &impl, const point &loc, const std::string &sym_op);
+
+		atom_impl(const atom_impl &i) = default;
+
+		void prefetch();
+
+		int compare(const atom_impl &b) const;
+
+		// bool getAnisoU(float anisou[6]) const;
+
+		// int charge() const;
+
+		void moveTo(const point &p)
+		{
+			if (m_symop != "1_555")
+				throw std::runtime_error("Moving symmetry copy");
+
+#if __cpp_lib_format
+			m_row.assign("Cartn_x", std::format("{:.3f}", p.getX()), false, false);
+			m_row.assign("Cartn_y", std::format("{:.3f}", p.getY()), false, false);
+			m_row.assign("Cartn_z", std::format("{:.3f}", p.getZ()), false, false);
+#else
+			m_row.assign("Cartn_x", format("%.3f", p.m_x).str(), false, false);
+			m_row.assign("Cartn_y", format("%.3f", p.m_y).str(), false, false);
+			m_row.assign("Cartn_z", format("%.3f", p.m_z).str(), false, false);
+#endif
+			m_location = p;
+		}
+
+		// const compound *compound() const;
+
+		std::string get_property(std::string_view name) const
+		{
+			return m_row[name].as<std::string>();
+			// 	for (auto &&[tag, ref] : mCachedRefs)
+			// 	{
+			// 		if (tag == name)
+			// 			return ref.as<std::string>();
+			// 	}
+
+			// 	mCachedRefs.emplace_back(name, const_cast<Row &>(mRow)[name]);
+			// return std::get<1>(mCachedRefs.back()).as<std::string>();
+		}
+
+		int get_property_int(std::string_view name) const
+		{
+			int result = 0;
+			if (not m_row[name].empty())
+			{
+				auto s = get_property(name);
+
+				std::from_chars_result r = std::from_chars(s.data(), s.data() + s.length(), result);
+				if (r.ec != std::errc() and VERBOSE > 0)
+					std::cerr << "Error converting " << s << " to number for property " << name << std::endl;
+			}
+			return result;
+		}
+
+		void set_property(const std::string_view name, const std::string &value)
+		{
+			m_row.assign(name, value, true, true);
+		}
+
+		// const datablock &m_db;
+		// std::string mID;
+		// atom_type mType;
+
+		// std::string mAtomID;
+		// std::string mCompID;
+		// std::string m_asym_id;
+		// int m_seq_id;
+		// std::string mAltID;
+		// std::string m_auth_seq_id;
+
+		// point mLocation;
+		// row_handle mRow;
+
+		// // mutable std::vector<std::tuple<std::string, detail::ItemReference>> mCachedRefs;
+
+		// mutable const compound *mcompound = nullptr;
+
+		// bool mSymmetryCopy = false;
+		// bool mClone = false;
+
+		const datablock &m_db;
+		std::string m_id;
+		row_handle m_row;
+		point m_location;
+		std::string m_symop = "1_555";
+	};
+
+  public:
+	atom() {}
+
+	atom(std::shared_ptr<atom_impl> impl)
+		: m_impl(impl)
+	{
+	}
+
+	atom(const atom &rhs)
+		: m_impl(rhs.m_impl)
+	{
+	}
+
+	atom(datablock &db, row_handle &row)
+		: atom(std::make_shared<atom_impl>(db, row["id"].as<std::string>(), row))
+	{
+	}
+
+	// a special constructor to create symmetry copies
+	atom(const atom &rhs, const point &symmmetry_location, const std::string &symmetry_operation)
+		: atom(std::make_shared<atom_impl>(*rhs.m_impl, symmmetry_location, symmetry_operation))
+	{
+	}
+
+	explicit operator bool() const { return (bool)m_impl; }
+
+	// // return a copy of this atom, with data copied instead of referenced
+	// atom clone() const
+	// {
+	// 	auto copy = std::make_shared<atom_impl>(*m_impl);
+	// 	copy->mClone = true;
+	// 	return atom(copy);
+	// }
+
+	atom &operator=(const atom &rhs) = default;
+
+	// template <typename T>
+	// T get_property(const std::string_view name) const;
+
+	std::string get_property(std::string_view name) const
+	{
+		if (not m_impl)
+			throw std::logic_error("Error trying to fetch a property from an uninitialized atom");
+		return m_impl->get_property(name);
+	}
+
+	int get_property_int(std::string_view name) const
+	{
+		if (not m_impl)
+			throw std::logic_error("Error trying to fetch a property from an uninitialized atom");
+		return m_impl->get_property_int(name);
+	}
+
+	void set_property(const std::string_view name, const std::string &value)
+	{
+		if (not m_impl)
+			throw std::logic_error("Error trying to modify an uninitialized atom");
+		m_impl->set_property(name, value);
+	}
+
+	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
+	void set_property(const std::string_view name, const T &value)
+	{
+		set_property(name, std::to_string(value));
+	}
+
+	const std::string &id() const { return impl().m_id; }
+	// AtomType type() const { return impl().mType; }
+
+	point get_location() const { return impl().m_location; }
+	void set_location(point p)
+	{
+		if (not m_impl)
+			throw std::logic_error("Error trying to modify an uninitialized atom");
+		m_impl->moveTo(p);
+	}
+
+	/// \brief Translate the position of this atom by \a t
+	void translate(point t)
+	{
+		set_location(get_location() + t);
+	}
+
+	/// \brief Rotate the position of this atom by \a q
+	void rotate(quaternion q)
+	{
+		auto loc = get_location();
+		loc.rotate(q);
+		set_location(loc);
+	}
+
+	/// \brief Translate and rotate the position of this atom by \a t and \a q
+	void translate_and_rotate(point t, quaternion q)
+	{
+		auto loc = get_location();
+		loc += t;
+		loc.rotate(q);
+		set_location(loc);
+	}
+
+	/// \brief Translate, rotate and translate again the coordinates this atom by \a t1 , \a q and \a t2
+	void translate_rotate_and_translate(point t1, quaternion q, point t2)
+	{
+		auto loc = get_location();
+		loc += t1;
+		loc.rotate(q);
+		loc += t2;
+		set_location(loc);
+	}
+
+	// // for direct access to underlying data, be careful!
+	// const row_handle getRow() const { return impl().mRow; }
+	// const row_handle getRowAniso() const;
+
+	// bool isSymmetryCopy() const { return impl().mSymmetryCopy; }
+	// std::string symmetry() const { return impl().mSymmetryOperator; }
+
+	// const compound &compound() const;
+	// bool isWater() const { return impl().mCompID == "HOH" or impl().mCompID == "H2O" or impl().mCompID == "WAT"; }
+	// int charge() const;
+
+	// float uIso() const;
+	// bool getAnisoU(float anisou[6]) const { return impl().getAnisoU(anisou); }
+	// float occupancy() const;
+
+	// specifications
+
+	std::string get_label_asym_id() const { return get_property("label_asym_id"); }
+	int get_label_seq_id() const { return get_property_int("label_seq_id"); }
+	std::string get_label_atom_id() const { return get_property("label_atom_id"); }
+	std::string get_label_alt_id() const { return get_property("label_alt_id"); }
+	std::string get_label_comp_id() const { return get_property("label_comp_id"); }
+	std::string get_label_entity_id() const { return get_property("label_entity_id"); }
+
+	std::string get_auth_asym_id() const { return get_property("auth_asym_id"); }
+	std::string get_auth_seq_id() const { return get_property("auth_seq_id"); }
+	std::string get_pdb_ins_code() const { return get_property("pdbx_PDB_ins_code"); }
+
+	// const std::string &labelAtomID() const { return impl().mAtomID; }
+	// const std::string &get_label_comp_id() const { return impl().mCompID; }
+	// const std::string &get_label_asym_id() const { return impl().m_asym_id; }
+	// std::string labelEntityID() const;
+	// int get_label_seq_id() const { return impl().m_seq_id; }
+	// const std::string &labelAltID() const { return impl().mAltID; }
+
+	bool is_alternate() const { return not get_label_alt_id().empty(); }
+
+	// std::string authAtomID() const;
+	// std::string authCompID() const;
+	// std::string authAsymID() const;
+	// const std::string &authSeqID() const { return impl().m_auth_seq_id; }
+	// std::string pdbxAuthInsCode() const;
+	// std::string pdbxAuthAltID() const;
+
+	// std::string labelID() const; // label_comp_id + '_' + label_asym_id + '_' + label_seq_id
+	// std::string pdbID() const;   // auth_comp_id + '_' + auth_asym_id + '_' + auth_seq_id + pdbx_PDB_ins_code
+
+	bool operator==(const atom &rhs) const
+	{
+		if (m_impl == rhs.m_impl)
+			return true;
+
+		if (not(m_impl and rhs.m_impl))
+			return false;
+
+		return &m_impl->m_db == &rhs.m_impl->m_db and m_impl->m_id == rhs.m_impl->m_id;
+	}
+
+	bool operator!=(const atom &rhs) const
+	{
+		return not operator==(rhs);
+	}
+
+	// // access data in compound for this atom
+
+	// // convenience routine
+	// bool isBackBone() const
+	// {
+	// 	auto atomID = labelAtomID();
+	// 	return atomID == "N" or atomID == "O" or atomID == "C" or atomID == "CA";
+	// }
+
+	void swap(atom &b)
+	{
+		std::swap(m_impl, b.m_impl);
+	}
+
+	int compare(const atom &b) const { return impl().compare(*b.m_impl); }
+
+	bool operator<(const atom &rhs) const
+	{
+		return compare(rhs) < 0;
+	}
+
+	friend std::ostream &operator<<(std::ostream &os, const atom &atom);
+
+	// /// \brief Synchronize data with underlying cif data
+	// void sync()
+	// {
+	// 	if (m_impl)
+	// 		m_impl->prefetch();
+	// }
+
+  private:
+	friend class structure;
+
+	const atom_impl &impl() const
+	{
+		if (not m_impl)
+			throw std::runtime_error("Uninitialized atom, not found?");
+		return *m_impl;
+	}
+
+	std::shared_ptr<atom_impl> m_impl;
+};
+
+// template <>
+// inline std::string atom::get_property<std::string>(const std::string_view name) const
+// {
+// 	return get_property(name);
+// }
+
+// template <>
+// inline int atom::get_property<int>(const std::string_view name) const
+// {
+// 	auto v = impl().get_property(name);
+// 	return v.empty() ? 0 : stoi(v);
+// }
+
+// template <>
+// inline float atom::get_property<float>(const std::string_view name) const
+// {
+// 	return stof(impl().get_property(name));
+// }
+
+inline void swap(atom &a, atom &b)
+{
+	a.swap(b);
+}
+
+inline float distance(const atom &a, const atom &b)
+{
+	return distance(a.get_location(), b.get_location());
+}
+
+inline float distance_squared(const atom &a, const atom &b)
+{
+	return distance_squared(a.get_location(), b.get_location());
+}
+
+// --------------------------------------------------------------------
+
+enum class EntityType
+{
+	polymer,
+	NonPolymer,
+	Macrolide,
+	Water,
+	Branched
+};
+
+// --------------------------------------------------------------------
+
+class residue
+{
+  public:
+	friend class structure;
+
+	// constructor
+	residue(const structure &structure, const std::string &compoundID,
+		const std::string &asymID, int seqID, const std::string &authSeqID)
+		: m_structure(&structure)
+		, m_compound_id(compoundID)
+		, m_asym_id(asymID)
+		, m_seq_id(seqID)
+		, m_auth_seq_id(authSeqID)
+	{
+	}
+
+	residue(const residue &rhs) = delete;
+	residue &operator=(const residue &rhs) = delete;
+
+	residue(residue &&rhs) = default;
+	residue &operator=(residue &&rhs) = default;
+
+	virtual ~residue() = default;
+
+	const std::string &get_asym_id() const { return m_asym_id; }
+	int get_seq_id() const { return m_seq_id; }
+
+	const std::string get_auth_asym_id() const
+	{
+		return m_atoms.empty() ? m_asym_id : m_atoms.front().get_auth_asym_id();
+	}
+	const std::string get_auth_seq_id() const { return m_auth_seq_id; }
+
+	const std::string &get_compound_id() const { return m_compound_id; }
+	void set_compound_id(const std::string &id) { m_compound_id = id; }
+
+	const structure *get_structure() const { return m_structure; }
+
+	// const compound &compound() const;
+
+	std::vector<atom> &atoms()
+	{
+		return m_atoms;
+	}
+
+	const std::vector<atom> &atoms() const
+	{
+		return m_atoms;
+	}
+
+	void add_atom(atom &atom);
+
+	// /// \brief Unique atoms returns only the atoms without alternates and the first of each alternate atom id.
+	// std::vector<atom> unique_atoms() const;
+
+	// /// \brief The alt ID used for the unique atoms
+	// std::string unique_alt_id() const;
+
+	atom get_atom_by_atom_id(const std::string &atomID) const;
+
+	// const std::string &asymID() const { return m_asym_id; }
+	// int seqID() const { return m_seq_id; }
+	std::string get_entity_id() const;
+
+	EntityType entity_type() const;
+
+	// std::string authAsymID() const;
+	// std::string authSeqID() const;
+	// std::string authInsCode() const;
+
+	// // return a human readable PDB-like auth id (chain+seqnr+iCode)
+	// std::string authID() const;
+
+	// // similar for mmCIF space
+	// std::string labelID() const;
+
+	// Is this residue a single entity?
+	bool is_entity() const;
+
+	// bool isWater() const { return m_compound_id == "HOH"; }
+
+	// bool empty() const { return m_structure == nullptr; }
+
+	// bool hasAlternateAtoms() const;
+
+	// /// \brief Return the list of unique alt ID's present in this residue
+	// std::set<std::string> getAlternateIDs() const;
+
+	// /// \brief Return the list of unique atom ID's
+	// std::set<std::string> getAtomIDs() const;
+
+	// /// \brief Return the list of atoms having ID \a atomID
+	// std::vector<atom> getAtomsByID(const std::string &atomID) const;
+
+	// // some routines for 3d work
+	// std::tuple<point, float> centerAndRadius() const;
+
+	friend std::ostream &operator<<(std::ostream &os, const residue &res);
+
+	bool operator==(const residue &rhs) const
+	{
+		return this == &rhs or (m_structure == rhs.m_structure and
+								   m_seq_id == rhs.m_seq_id and
+								   m_asym_id == rhs.m_asym_id and
+								   m_compound_id == rhs.m_compound_id and
+								   m_auth_seq_id == rhs.m_auth_seq_id);
+	}
+
+  protected:
+	residue() {}
+
+	const structure *m_structure = nullptr;
+	std::string m_compound_id, m_asym_id;
+	int m_seq_id = 0;
+	std::string m_auth_seq_id;
+	std::vector<atom> m_atoms;
+};
+
+// --------------------------------------------------------------------
+// a monomer models a single residue in a protein chain
+
+class monomer : public residue
+{
+  public:
+	//	monomer();
+	monomer(const monomer &rhs) = delete;
+	monomer &operator=(const monomer &rhs) = delete;
+
+	monomer(monomer &&rhs);
+	monomer &operator=(monomer &&rhs);
+
+	monomer(const polymer &polymer, size_t index, int seqID, const std::string &authSeqID,
+		const std::string &compoundID);
+
+	// bool is_first_in_chain() const;
+	// bool is_last_in_chain() const;
+
+	// // convenience
+	// bool has_alpha() const;
+	// bool has_kappa() const;
+
+	// // Assuming this is really an amino acid...
+
+	// float phi() const;
+	// float psi() const;
+	// float alpha() const;
+	// float kappa() const;
+	// float tco() const;
+	// float omega() const;
+
+	// // torsion angles
+	// size_t nrOfChis() const;
+	// float chi(size_t i) const;
+
+	// bool isCis() const;
+
+	// /// \brief Returns true if the four atoms C, CA, N and O are present
+	// bool isComplete() const;
+
+	// /// \brief Returns true if any of the backbone atoms has an alternate
+	// bool hasAlternateBackboneAtoms() const;
+
+	// atom CAlpha() const { return get_atom_by_atom_id("CA"); }
+	// atom C() const { return get_atom_by_atom_id("C"); }
+	// atom N() const { return get_atom_by_atom_id("N"); }
+	// atom O() const { return get_atom_by_atom_id("O"); }
+	// atom H() const { return get_atom_by_atom_id("H"); }
+
+	// bool isBondedTo(const monomer &rhs) const
+	// {
+	// 	return this != &rhs and areBonded(*this, rhs);
+	// }
+
+	// static bool areBonded(const monomer &a, const monomer &b, float errorMargin = 0.5f);
+	// static bool isCis(const monomer &a, const monomer &b);
+	// static float omega(const monomer &a, const monomer &b);
+
+	// // for LEU and VAL
+	// float chiralVolume() const;
+
+	bool operator==(const monomer &rhs) const
+	{
+		return m_polymer == rhs.m_polymer and m_index == rhs.m_index;
+	}
+
+  private:
+	const polymer *m_polymer;
+	size_t m_index;
+};
+
+// --------------------------------------------------------------------
+
+class polymer : public std::vector<monomer>
+{
+  public:
+	polymer(const structure &s, const std::string &entityID, const std::string &asymID);
+
+	polymer(const polymer &) = delete;
+	polymer &operator=(const polymer &) = delete;
+
+	// monomer &getBySeqID(int seqID);
+	// const monomer &getBySeqID(int seqID) const;
+
+	const structure *get_structure() const { return m_structure; }
+
+	std::string get_asym_id() const { return m_asym_id; }
+	std::string get_entity_id() const { return m_entity_id; }
+
+	// std::string chainID() const;
+
+	// int Distance(const monomer &a, const monomer &b) const;
+
+  private:
+	const structure *m_structure;
+	std::string m_entity_id;
+	std::string m_asym_id;
+};
+
+// --------------------------------------------------------------------
+// sugar and branch, to describe glycosylation sites
+
+class branch;
+
+class sugar : public residue
+{
+  public:
+	sugar(const branch &branch, const std::string &compoundID,
+		const std::string &asymID, int authSeqID);
+
+	sugar(sugar &&rhs);
+	sugar &operator=(sugar &&rhs);
+
+	int num() const { return std::stoi(m_auth_seq_id); }
+	std::string name() const;
+
+	/// \brief Return the atom the C1 is linked to
+	atom get_link() const { return m_link; }
+	void set_link(atom link) { m_link = link; }
+
+	size_t get_link_nr() const
+	{
+		return m_link ? std::stoi(m_link.get_auth_seq_id()) : 0;
+	}
+
+  private:
+	const branch *m_branch;
+	atom m_link;
+};
+
+class branch : public std::vector<sugar>
+{
+  public:
+	branch(structure &structure, const std::string &asymID);
+
+	void link_atoms();
+
+	std::string name() const;
+	float weight() const;
+	std::string get_asym_id() const { return m_asym_id; }
+
+	structure &get_structure() { return *m_structure; }
+	const structure &get_structure() const { return *m_structure; }
+
+	sugar &getSugarByNum(int nr);
+	const sugar &getSugarByNum(int nr) const;
+
+  private:
+	friend sugar;
+
+	std::string name(const sugar &s) const;
+
+	structure *m_structure;
+	std::string m_asym_id;
+};
+
+// // --------------------------------------------------------------------
+// // file is a reference to the data stored in e.g. the cif file.
+// // This object is not copyable.
+
+// class File : public file
+// {
+//   public:
+// 	File() {}
+
+// 	// File(const std::filesystem::path &path)
+// 	// {
+// 	// 	load(path);
+// 	// }
+
+// 	// File(const char *data, size_t length)
+// 	// {
+// 	// 	load(data, length);
+// 	// }
+
+// 	File(const File &) = delete;
+// 	File &operator=(const File &) = delete;
+
+// 	// void load(const std::filesystem::path &p) override;
+// 	// void save(const std::filesystem::path &p) override;
+
+// 	// using file::load;
+// 	// using file::save;
+
+// 	datablock &data() { return front(); }
+// };
+
+// --------------------------------------------------------------------
+
+enum class StructureOpenOptions
+{
+	SkipHydrogen = 1 << 0
+};
+
+inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
+{
+	return static_cast<int>(a) bitand static_cast<int>(b);
+}
+
+// --------------------------------------------------------------------
+
+class structure
+{
+  public:
+	structure(file &p, size_t modelNr = 1, StructureOpenOptions options = {});
+
+	structure(datablock &db, size_t modelNr = 1, StructureOpenOptions options = {});
+
+	structure(structure &&s) = default;
+
+	// Create a read-only clone of the current structure (for multithreaded calculations that move atoms)
+	// NOTE: removed, simply create a new structure for each thread
+	structure(const structure &) = delete;
+
+	structure &operator=(const structure &) = delete;
+	// Structure &operator=(Structure &&s) = default;
+
+	~structure() = default;
+
+	const std::vector<atom> &atoms() const { return m_atoms; }
+	// std::vector<atom> &atoms() { return m_atoms; }
+
+	EntityType get_entity_type_for_entity_id(const std::string entityID) const;
+	EntityType get_entity_type_for_asym_id(const std::string asymID) const;
+
+	// std::vector<atom> waters() const;
+
+	const std::list<polymer> &polymers() const { return m_polymers; }
+	std::list<polymer> &polymers() { return m_polymers; }
+
+	// polymer &getPolymerByAsymID(const std::string &asymID);
+
+	// const polymer &getPolymerByAsymID(const std::string &asymID) const
+	// {
+	// 	return const_cast<structure *>(this)->getPolymerByAsymID(asymID);
+	// }
+
+	const std::list<branch> &branches() const { return m_branches; }
+	std::list<branch> &branches() { return m_branches; }
+
+	branch &get_branch_by_asym_id(const std::string &asymID);
+	const branch &get_branch_by_asym_id(const std::string &asymID) const;
+
+	const std::vector<residue> &non_polymers() const { return m_non_polymers; }
+
+	atom get_atom_by_id(const std::string &id) const;
+	// atom getAtomByLocation(point pt, float maxDistance) const;
+
+	// atom getAtomByLabel(const std::string &atomID, const std::string &asymID,
+	// 	const std::string &compID, int seqID, const std::string &altID = "");
+
+	// /// \brief Return the atom closest to point \a p
+	// atom getAtomByPosition(point p) const;
+
+	// /// \brief Return the atom closest to point \a p with atom type \a type in a residue of type \a res_type
+	// atom getAtomByPositionAndType(point p, std::string_view type, std::string_view res_type) const;
+
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	residue &get_residue(const std::string &asymID)
+	{
+		return get_residue(asymID, 0, "");
+	}
+
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	const residue &get_residue(const std::string &asymID) const
+	{
+		return get_residue(asymID, 0, "");
+	}
+
+	/// \brief Get a residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
+	residue &get_residue(const std::string &asymID, int seqID, const std::string &authSeqID);
+
+	/// \brief Get a the single residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
+	const residue &get_residue(const std::string &asymID, int seqID, const std::string &authSeqID) const
+	{
+		return const_cast<structure *>(this)->get_residue(asymID, seqID, authSeqID);
+	}
+
+	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
+	residue &get_residue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID);
+
+	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
+	const residue &get_residue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID) const
+	{
+		return const_cast<structure *>(this)->get_residue(asymID, compID, seqID, authSeqID);
+	}
+
+	/// \brief Get a the residue for atom \a atom
+	residue &get_residue(const atom &atom)
+	{
+		return get_residue(atom.get_label_asym_id(), atom.get_label_comp_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
+	}
+
+	/// \brief Get a the residue for atom \a atom
+	const residue &get_residue(const atom &atom) const
+	{
+		return get_residue(atom.get_label_asym_id(), atom.get_label_comp_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
+	}
+
+	// Actions
+	void remove_atom(atom &a)
+	{
+		remove_atom(a, true);
+	}
+
+	void swap_atoms(atom a1, atom a2); // swap the labels for these atoms
+	void move_atom(atom a, point p);   // move atom to a new location
+	void change_residue(residue &res, const std::string &newcompound,
+		const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);
+
+	/// \brief Remove a residue, can be monomer or nonpoly
+	///
+	/// \param asym_id     The asym ID
+	/// \param seq_id      The sequence ID
+	void remove_residue(const std::string &asym_id, int seq_id, const std::string &auth_seq_id)
+	{
+		remove_residue(get_residue(asym_id, seq_id, auth_seq_id));
+	}
+
+	/// \brief Create a new non-polymer entity, returns new ID
+	/// \param mon_id	The mon_id for the new nonpoly, must be an existing and known compound from CCD
+	/// \return			The ID of the created entity
+	std::string create_non_poly_entity(const std::string &mon_id);
+
+	/// \brief Create a new NonPolymer struct_asym with atoms constructed from \a atoms, returns asym_id.
+	/// This method assumes you are copying data from one cif file to another.
+	///
+	/// \param entity_id	The entity ID of the new nonpoly
+	/// \param atoms		The array of atom_site rows containing the data.
+	/// \return				The newly create asym ID
+	std::string create_non_poly(const std::string &entity_id, const std::vector<atom> &atoms);
+
+	/// \brief Create a new NonPolymer struct_asym with atoms constructed from info in \a atom_info, returns asym_id.
+	/// This method creates new atom records filled with info from the info.
+	///
+	/// \param entity_id	The entity ID of the new nonpoly
+	/// \param atoms		The array of sets of item data containing the data for the atoms.
+	/// \return				The newly create asym ID
+	std::string create_non_poly(const std::string &entity_id, std::vector<std::vector<item>> &atom_info);
+
+	/// \brief Create a new (sugar) branch with one first NAG containing atoms constructed from \a nag_atom_info
+	branch &create_branch(std::vector<std::vector<item>> &nag_atom_info);
+
+	/// \brief Extend an existing (sugar) branch identified by \a asymID with one sugar containing atoms constructed from \a atom_info
+	///
+	/// \param asym_id      The asym id of the branch to extend
+	/// \param atom_info    Array containing the info for the atoms to construct for the new sugar
+	/// \param link_sugar   The sugar to link to, note: this is the sugar number (1 based)
+	/// \param link_atom    The atom id of the atom linked in the sugar
+	branch &extend_branch(const std::string &asym_id, std::vector<std::vector<item>> &atom_info,
+		int link_sugar, const std::string &link_atom);
+
+	/// \brief Remove \a branch
+	void remove_branch(branch &branch);
+
+	/// \brief Remove residue \a res
+	///
+	/// \param res         The residue to remove
+	void remove_residue(residue &res);
+
+	/// \brief Translate the coordinates of all atoms in the structure by \a t
+	void translate(point t);
+
+	/// \brief Rotate the coordinates of all atoms in the structure by \a q
+	void rotate(quaternion t);
+
+	/// \brief Translate and rotate the coordinates of all atoms in the structure by \a t and \a q
+	void translate_and_rotate(point t, quaternion q);
+
+	/// \brief Translate, rotate and translate again the coordinates of all atoms in the structure by \a t1 , \a q and \a t2
+	void translate_rotate_and_translate(point t1, quaternion q, point t2);
+
+	void cleanup_empty_categories();
+
+	// /// \brief Direct access to underlying data
+	// category &category(std::string_view name) const
+	// {
+	// 	return m_db[name];
+	// }
+
+	datablock &get_datablock() const
+	{
+		return m_db;
+	}
+
+	void validate_atoms() const;
+
+  private:
+	friend polymer;
+	friend residue;
+
+	std::string insert_compound(const std::string &compoundID, bool is_entity);
+
+	std::string create_entity_for_branch(branch &branch);
+
+	void loadData();
+
+	void load_atoms_for_model(StructureOpenOptions options);
+
+	template <typename... Args>
+	atom &emplace_atom(Args... args)
+	{
+		return emplace_atom(atom{ std::forward<Args>(args)... });
+	}
+
+	atom &emplace_atom(atom &&atom);
+
+	void remove_atom(atom &a, bool removeFromResidue);
+	void remove_sugar(sugar &sugar);
+
+	datablock &m_db;
+	size_t m_model_nr;
+	std::vector<atom> m_atoms;
+	std::vector<size_t> m_atom_index;
+	std::list<polymer> m_polymers;
+	std::list<branch> m_branches;
+	std::vector<residue> m_non_polymers;
+};
+
+} // namespace cif::mm

include/cif++/row.hpp

@@ -126,6 +126,35 @@ class row
 			m_tail = m_tail->m_next = iv;
 	}
 
+	void remove(item_value *iv)
+	{
+		if (iv == m_head)
+		{
+			if (m_tail == m_head)
+			{
+				assert(iv->m_next == nullptr);
+				m_head = m_tail = nullptr;
+			}
+			else
+				m_head = m_head->m_next;
+		}
+		else
+		{
+			for (auto v = m_head; v->m_next != nullptr; v = v->m_next)
+			{
+				if (v->m_next != iv)
+					continue;
+
+				v->m_next = iv->m_next;
+				iv->m_next = nullptr;
+
+				if (m_tail == iv)
+					m_tail = v;
+				break;
+			}
+		}
+	}
+
 	row *m_next = nullptr;
 	item_value *m_head = nullptr, *m_tail = nullptr;
 };
@@ -240,6 +269,8 @@ class row_handle
 		assign(i.name(), i.value(), updateLinked);
 	}
 
+	void swap(size_t column, row_handle &r);
+
 	category *m_category = nullptr;
 	row *m_row = nullptr;
 };

include/cif++/symmetry.hpp

@@ -50,8 +50,8 @@ struct space_group
 	int nr;
 };
 
-// extern const Spacegroup kSpaceGroups[];
-// extern const std::size_t kNrOfSpaceGroups;
+extern const space_group kSpaceGroups[];
+extern const std::size_t kNrOfSpaceGroups;
 
 // --------------------------------------------------------------------
 
@@ -134,8 +134,8 @@ struct symop_datablock
 
 static_assert(sizeof(symop_datablock) == sizeof(uint64_t), "Size of symop_data is wrong");
 
-// extern const symop_data_block kSymopNrTable[];
-// extern const std::size_t kSymopNrTableSize;
+extern const symop_datablock kSymopNrTable[];
+extern const std::size_t kSymopNrTableSize;
 
 // --------------------------------------------------------------------
 

include/cif++/utilities.hpp

@@ -170,7 +170,7 @@ class Progress
 	Progress(const Progress &) = delete;
 	Progress &operator=(const Progress &) = delete;
 
-	struct ProgressImpl *mImpl;
+	struct ProgressImpl *m_impl;
 };
 
 // --------------------------------------------------------------------

src/category.cpp

@@ -1328,46 +1328,18 @@ void category::update_value(row *row, size_t column, std::string_view value, boo
 	}
 
 	// first remove old value with cix
-
-	if (row->m_head == nullptr)
-		; // nothing to do
-	else if (row->m_head->m_column_ix == column)
-	{
-		auto iv = row->m_head;
-		row->m_head = iv->m_next;
-		iv->m_next = nullptr;
-		delete_item(iv);
-	}
-	else
+	for (auto iv = row->m_head; iv != nullptr; iv = iv->m_next)
 	{
-		for (auto iv = row->m_head; iv->m_next != nullptr; iv = iv->m_next)
-		{
-			if (iv->m_next->m_column_ix != column)
-				continue;
-
-			auto nv = iv->m_next;
-			iv->m_next = nv->m_next;
-			nv->m_next = nullptr;
-			delete_item(nv);
+		if (iv->m_column_ix != column)
+			continue;
 
-			break;
-		}
+		row->remove(iv);
+		delete_item(iv);
+		break;
 	}
 
 	if (not value.empty())
-	{
-		auto nv = create_item(column, value);
-
-		if (row->m_head == nullptr)
-			row->m_head = nv;
-		else
-		{
-			auto iv = row->m_head;
-			while (iv->m_next != nullptr)
-				iv = iv->m_next;
-			iv->m_next = nv;
-		}
-	}
+		row->append(create_item(column, value));
 
 	if (reinsert)
 		m_index->insert(row);
@@ -1671,6 +1643,87 @@ category::iterator category::erase_impl(const_iterator pos)
 	// return iterator(*this, cur);
 }
 
+void category::swap_item(size_t column_ix, row_handle &a, row_handle &b)
+{
+	assert(this == a.m_category);
+	assert(this == b.m_category);
+
+	item_value *va = nullptr, *vb = nullptr;
+
+	auto ra = a.m_row;
+	auto rb = b.m_row;
+
+	if (ra->m_head != nullptr and ra->m_head->m_column_ix == column_ix)
+	{
+		va = ra->m_head;
+		ra->m_head = va->m_next;
+		va->m_next = nullptr;
+
+		if (ra->m_tail == va)
+			ra->m_tail = ra->m_head;
+	}
+	else
+	{
+		for (auto v = ra->m_head; v->m_next != nullptr; v = v->m_next)
+		{
+			if (v->m_next->m_column_ix != column_ix)
+				continue;
+			
+			va = v->m_next;
+			v->m_next = va->m_next;
+			va->m_next = nullptr;
+
+			if (ra->m_tail == va)
+				ra->m_tail = v;
+
+			break;
+		}
+	}
+
+	if (rb->m_head != nullptr and rb->m_head->m_column_ix == column_ix)
+	{
+		vb = rb->m_head;
+		rb->m_head = vb->m_next;
+		vb->m_next = nullptr;
+
+		if (rb->m_tail == vb)
+			rb->m_tail = rb->m_head;
+	}
+	else
+	{
+		for (auto v = rb->m_head; v->m_next != nullptr; v = v->m_next)
+		{
+			if (v->m_next->m_column_ix != column_ix)
+				continue;
+			
+			vb = v->m_next;
+			v->m_next = vb->m_next;
+			vb->m_next = nullptr;
+
+			if (rb->m_tail = vb)
+				rb->m_tail = v;
+
+			break;
+		}
+	}
+
+	if (ra->m_head == nullptr)
+		ra->m_head = ra->m_tail = vb;
+	else
+	{
+		ra->m_tail->m_next = vb;
+		ra->m_tail = vb;
+	}
+
+	if (rb->m_head == nullptr)
+		rb->m_head = rb->m_tail = va;
+	else
+	{
+		rb->m_tail->m_next = va;
+		rb->m_tail = va;
+	}
+}
+
 void category::reorder_by_index()
 {
 	if (m_index)
@@ -1933,4 +1986,100 @@ void category::write(std::ostream &os, const std::vector<uint16_t> &order, bool
 	os << "# " << '\n';
 }
 
+bool category::operator==(const category &rhs) const
+{
+	auto &a = *this;
+	auto &b = rhs;
+
+	using namespace std::placeholders; 
+	
+//	set<std::string> tagsA(a.fields()), tagsB(b.fields());
+//	
+//	if (tagsA != tagsB)
+//		std::cout << "Unequal number of fields" << std::endl;
+
+	auto validator = a.get_validator();
+	auto catValidator = validator->get_validator_for_category(a.name());
+	if (catValidator == nullptr)
+		throw std::runtime_error("missing cat validator");
+	
+	typedef std::function<int(std::string_view,std::string_view)> compType;
+	std::vector<std::tuple<std::string,compType>> tags;
+	auto keys = catValidator->m_keys;
+	std::vector<size_t> keyIx;
+	
+	for (auto& tag: a.fields())
+	{
+		auto iv = catValidator->get_validator_for_item(tag);
+		if (iv == nullptr)
+			throw std::runtime_error("missing item validator");
+		auto tv = iv->m_type;
+		if (tv == nullptr)
+			throw std::runtime_error("missing type validator");
+		tags.push_back(std::make_tuple(tag, std::bind(&cif::type_validator::compare, tv, std::placeholders::_1, std::placeholders::_2)));
+		
+		auto pred = [tag](const std::string& s) -> bool { return cif::iequals(tag, s) == 0; };
+		if (find_if(keys.begin(), keys.end(), pred) == keys.end())
+			keyIx.push_back(tags.size() - 1);
+	}
+	
+	// a.reorderByIndex();
+	// b.reorderByIndex();
+	
+	auto rowEqual = [&](const row_handle& a, const row_handle& b)
+	{
+		int d = 0;
+
+		for (auto kix: keyIx)
+		{
+			std::string tag;
+			compType compare;
+			
+			std::tie(tag, compare) = tags[kix];
+
+			d = compare(a[tag].text(), b[tag].text());
+
+			if (d != 0)
+				break;
+		}
+		
+		return d == 0;
+	};
+
+	auto ai = a.begin(), bi = b.begin();
+	while (ai != a.end() or bi != b.end())
+	{
+		if (ai == a.end() or bi == b.end())
+			return false;
+		
+		auto ra = *ai, rb = *bi;
+		
+		if (not rowEqual(ra, rb))
+			return false;
+		
+		std::vector<std::string> missingA, missingB, different;
+		
+		for (auto& tt: tags)
+		{
+			std::string tag;
+			compType compare;
+			
+			std::tie(tag, compare) = tt;
+			
+			// make it an option to compare unapplicable to empty or something
+			
+			auto ta = ra[tag].text();	if (ta == "." or ta == "?") ta = "";
+			auto tb = rb[tag].text();	if (tb == "." or tb == "?") tb = "";
+			
+			if (compare(ta, tb) != 0)
+				return false;
+		}
+		
+		++ai;
+		++bi;
+	}
+
+	return true;
+}
+
 } // namespace cif
\ No newline at end of file

src/compound.cpp

@@ -194,7 +194,7 @@ compound::compound(cif::datablock &db, const std::string &id, const std::string
 	}
 }
 
-compound_atom compound::get_atom_by_id(const std::string &atom_id) const
+compound_atom compound::get_atom_by_atom_id(const std::string &atom_id) const
 {
 	compound_atom result = {};
 	for (auto &a : m_atoms)
@@ -367,9 +367,9 @@ compound_factory_impl::compound_factory_impl(const fs::path &file, std::shared_p
 {
 	cif::file cifFile(file);
 
-	auto &compList = cifFile["comp_list"];
-	if (not compList.empty()) // So this is a CCP4 restraints file, special handling
+	if (cifFile.contains("comp_list")) // So this is a CCP4 restraints file, special handling
 	{
+		auto &compList = cifFile["comp_list"];
 		auto &chemComp = compList["chem_comp"];
 
 		for (const auto &[id, name, group] : chemComp.rows<std::string, std::string, std::string>("id", "name", "group"))
@@ -726,7 +726,7 @@ const compound *compound_factory::create(std::string id)
 {
 	// static bool warned = false;
 
-	// if (mImpl and warned == false)
+	// if (m_impl and warned == false)
 	// {
 	// 	std::cerr << "Warning: no compound information library was found, resulting data may be incorrect or incomplete" << std::endl;
 	// 	warned = true;

src/datablock.cpp

@@ -241,4 +241,71 @@ void datablock::write(std::ostream &os, const std::vector<std::string> &tag_orde
 	}
 }
 
+bool datablock::operator==(const datablock &rhs) const
+{
+	auto &dbA = *this;
+	auto &dbB = rhs;
+
+	std::vector<std::string> catA, catB;
+
+	for (auto &cat : dbA)
+	{
+		if (not cat.empty())
+			catA.push_back(cat.name());
+	}
+	std::sort(catA.begin(), catA.end());
+
+	for (auto &cat : dbB)
+	{
+		if (not cat.empty())
+			catB.push_back(cat.name());
+	}
+	std::sort(catB.begin(), catB.end());
+
+	// loop over categories twice, to group output
+	// First iteration is to list missing categories.
+
+	std::vector<std::string> missingA, missingB;
+
+	auto catA_i = catA.begin(), catB_i = catB.begin();
+
+	while (catA_i != catA.end() and catB_i != catB.end())
+	{
+		if (not iequals(*catA_i, *catB_i))
+			return false;
+
+		++catA_i, ++catB_i;
+	}
+
+	if (catA_i != catA.end() or catB_i != catB.end())
+		return false;
+
+	// Second loop, now compare category values
+	catA_i = catA.begin(), catB_i = catB.begin();
+
+	while (catA_i != catA.end() and catB_i != catB.end())
+	{
+		std::string nA = *catA_i;
+		to_lower(nA);
+
+		std::string nB = *catB_i;
+		to_lower(nB);
+
+		int d = nA.compare(nB);
+		if (d > 0)
+			++catB_i;
+		else if (d < 0)
+			++catA_i;
+		else
+		{
+			if (not (*dbA.get(*catA_i) == *dbB.get(*catB_i)))
+				return false;
+			++catA_i;
+			++catB_i;
+		}
+	}
+
+	return true;
+}
+
 } // namespace cif::cif
\ No newline at end of file

src/file.cpp

@@ -111,6 +111,11 @@ void file::load_dictionary(std::string_view name)
 	set_validator(&validator_factory::instance()[name]);
 }
 
+bool file::contains(std::string_view name) const
+{
+	return std::find_if(begin(), end(), [name](const datablock &db) { return db.name() == name; }) != end();
+}
+
 datablock &file::operator[](std::string_view name)
 {
 	auto i = std::find_if(begin(), end(), [name](const datablock &c)

src/item.cpp

@@ -60,4 +60,11 @@ void item_handle::assign_value(const item &v)
 	m_row_handle.assign(m_column, v.value(), true);
 }
 
+void item_handle::swap(item_handle &b)
+{
+	assert(m_column == b.m_column);
+	// assert(&m_row_handle.m_category == &b.m_row_handle.m_category);
+	m_row_handle.swap(m_column, b.m_row_handle);
+}
+
 }

src/model.cpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * 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++.hpp>
+
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <numeric>
+
+#include <gxrio.hpp>
+
+namespace fs = std::filesystem;
+
+extern int VERBOSE;
+
+namespace cif::mm
+{
+
+// --------------------------------------------------------------------
+// atom
+
+// constructor for a symmetry copy of an atom
+atom::atom_impl::atom_impl(const atom_impl &impl, const point &loc, const std::string &sym_op)
+	: m_db(impl.m_db)
+	, m_id(impl.m_id)
+	, m_row(impl.m_row)
+	, m_location(loc)
+	, m_symop(sym_op)
+{
+}
+
+void atom::atom_impl::prefetch()
+{
+	tie(m_location.m_x, m_location.m_y, m_location.m_z) = m_row.get("Cartn_x", "Cartn_y", "Cartn_z");
+}
+
+// int atom::atom_impl::compare(const atom_impl &b) const
+// {
+// 	int d = m_asym_id.compare(b.m_asym_id);
+// 	if (d == 0)
+// 		d = m_seq_id - b.m_seq_id;
+// 	if (d == 0)
+// 		d = m_auth_seq_id.compare(b.m_auth_seq_id);
+// 	if (d == 0)
+// 		d = mAtom_id.compare(b.mAtom_id);
+
+// 	return d;
+// }
+
+// bool atom::atom_impl::getAnisoU(float anisou[6]) const
+// {
+// 	bool result = false;
+
+// 	auto cat = m_db.get("atom_site_anisotrop");
+// 	if (cat)
+// 	{
+// 		for (auto r : cat->find(key("id") == m_id))
+// 		{
+// 			tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
+// 				r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
+// 			result = true;
+// 			break;
+// 		}
+// 	}
+
+// 	return result;
+// }
+
+// int atom::atom_impl::charge() const
+// {
+// 	auto formalCharge = mRow["pdbx_formal_charge"].as<std::optional<int>>();
+
+// 	if (not formalCharge.has_value())
+// 	{
+// 		auto c = compound();
+
+// 		if (c != nullptr and c->atoms().size() == 1)
+// 			formalCharge = c->atoms().front().charge;
+// 	}
+
+// 	return formalCharge.value_or(0);
+// }
+
+// const Compound *atom::atom_impl::compound() const
+// {
+// 	if (mCompound == nullptr)
+// 	{
+// 		std::string compID = get_property("label_comp_id");
+
+// 		mCompound = compound_factory::instance().create(compID);
+// 	}
+
+// 	return mCompound;
+// }
+
+// const std::string atom::atom_impl::get_property(const std::string_view name) const
+// {
+// 	for (auto &&[tag, ref] : mCachedRefs)
+// 	{
+// 		if (tag == name)
+// 			return ref.as<std::string>();
+// 	}
+
+// 	mCachedRefs.emplace_back(name, const_cast<Row &>(mRow)[name]);
+// 	return std::get<1>(mCachedRefs.back()).as<std::string>();
+// }
+
+// void atom::atom_impl::set_property(const std::string_view name, const std::string &value)
+// {
+// 	for (auto &&[tag, ref] : mCachedRefs)
+// 	{
+// 		if (tag != name)
+// 			continue;
+
+// 		ref = value;
+// 		return;
+// 	}
+
+// 	mCachedRefs.emplace_back(name, mRow[name]);
+// 	std::get<1>(mCachedRefs.back()) = value;
+// }
+
+// const Row atom::getRowAniso() const
+// {
+// 	auto &db = m_impl->m_db;
+// 	auto cat = db.get("atom_site_anisotrop");
+// 	if (not cat)
+// 		return {};
+// 	else
+// 		return cat->find1(key("id") == m_impl->m_id);
+// }
+
+// float atom::uIso() const
+// {
+// 	float result;
+
+// 	if (not get_property<std::string>("U_iso_or_equiv").empty())
+// 		result = get_property<float>("U_iso_or_equiv");
+// 	else if (not get_property<std::string>("B_iso_or_equiv").empty())
+// 		result = get_property<float>("B_iso_or_equiv") / static_cast<float>(8 * kPI * kPI);
+// 	else
+// 		throw std::runtime_error("Missing B_iso or U_iso");
+
+// 	return result;
+// }
+
+// std::string atom::labelID() const
+// {
+// 	return m_impl->mCompID + '_' + m_impl->m_asym_id + '_' + std::to_string(m_impl->m_seq_id) + ':' + m_impl->mAtom_id;
+// }
+
+// std::string atom::pdbID() const
+// {
+// 	return get_property<std::string>("auth_comp_id") + '_' +
+// 	       get_property<std::string>("auth_asym_id") + '_' +
+// 	       get_property<std::string>("auth_seq_id") +
+// 	       get_property<std::string>("pdbx_PDB_ins_code");
+// }
+
+// const Compound &atom::compound() const
+// {
+// 	auto result = impl().compound();
+
+// 	if (result == nullptr)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << "Compound not found: '" << get_property<std::string>("label_comp_id") << '\'' << std::endl;
+
+// 		throw std::runtime_error("no compound");
+// 	}
+
+// 	return *result;
+// }
+
+// int atom::charge() const
+// {
+// 	return impl().charge();
+// }
+
+// float atom::occupancy() const
+// {
+// 	return get_property<float>("occupancy");
+// }
+
+// std::string atom::labelEntityID() const
+// {
+// 	return get_property<std::string>("label_entity_id");
+// }
+
+// std::string atom::authAtom_id() const
+// {
+// 	return get_property<std::string>("auth_atom_id");
+// }
+
+// std::string atom::authCompID() const
+// {
+// 	return get_property<std::string>("auth_comp_id");
+// }
+
+// std::string atom::get_auth_asym_id() const
+// {
+// 	return get_property<std::string>("auth_asym_id");
+// }
+
+// std::string atom::get_pdb_ins_code() const
+// {
+// 	return get_property<std::string>("pdbx_PDB_ins_code");
+// }
+
+// std::string atom::pdbxAuthAltID() const
+// {
+// 	return get_property<std::string>("pdbx_auth_alt_id");
+// }
+
+// void atom::translate(point t)
+// {
+// 	auto loc = location();
+// 	loc += t;
+// 	location(loc);
+// }
+
+// void atom::rotate(quaternion q)
+// {
+// 	auto loc = location();
+// 	loc.rotate(q);
+// 	location(loc);
+// }
+
+// void atom::translate_and_rotate(point t, quaternion q)
+// {
+// 	auto loc = location();
+// 	loc += t;
+// 	loc.rotate(q);
+// 	location(loc);
+// }
+
+// void atom::translate_rotate_and_translate(point t1, quaternion q, point t2)
+// {
+// 	auto loc = location();
+// 	loc += t1;
+// 	loc.rotate(q);
+// 	loc += t2;
+// 	location(loc);
+// }
+
+std::ostream &operator<<(std::ostream &os, const atom &atom)
+{
+	os << atom.get_label_comp_id() << ' ' << atom.get_label_asym_id() << ':' << atom.get_label_seq_id() << ' ' << atom.get_label_atom_id();
+
+	if (atom.is_alternate())
+		os << '(' << atom.get_label_alt_id() << ')';
+	if (atom.get_auth_asym_id() != atom.get_label_asym_id() or atom.get_auth_seq_id() != std::to_string(atom.get_label_seq_id()) or atom.get_pdb_ins_code().empty() == false)
+		os << " [" << atom.get_auth_asym_id() << ':' << atom.get_auth_seq_id() << atom.get_pdb_ins_code() << ']';
+
+	return os;
+}
+
+// --------------------------------------------------------------------
+// residue
+
+// residue::residue(residue &&rhs)
+// 	: m_structure(rhs.m_structure)
+// 	, m_compound_id(std::move(rhs.m_compound_id))
+// 	, m_asym_id(std::move(rhs.m_asym_id))
+// 	, m_seq_id(rhs.m_seq_id)
+// 	, m_auth_seq_id(rhs.m_auth_seq_id)
+// 	, m_atoms(std::move(rhs.m_atoms))
+// {
+// 	// std::cerr << "move constructor residue" << std::endl;
+// 	rhs.m_structure = nullptr;
+// }
+
+// residue &residue::operator=(residue &&rhs)
+// {
+// 	// std::cerr << "move assignment residue" << std::endl;
+// 	m_structure = rhs.m_structure;
+// 	rhs.m_structure = nullptr;
+// 	m_compound_id = std::move(rhs.m_compound_id);
+// 	m_asym_id = std::move(rhs.m_asym_id);
+// 	m_seq_id = rhs.m_seq_id;
+// 	m_auth_seq_id = rhs.m_auth_seq_id;
+// 	m_atoms = std::move(rhs.m_atoms);
+
+// 	return *this;
+// }
+
+// residue::~residue()
+// {
+// 	// std::cerr << "~residue" << std::endl;
+// }
+
+std::string residue::get_entity_id() const
+{
+	std::string result;
+
+	if (not m_atoms.empty())
+		result = m_atoms.front().get_label_entity_id();
+	else if (m_structure != nullptr and not m_asym_id.empty())
+	{
+		using namespace literals;
+
+		auto &db = m_structure->get_datablock();
+		result = db["struct_asym"].find1<std::string>("id"_key == m_asym_id, "entity_id");
+	}
+
+	return result;
+}
+
+EntityType residue::entity_type() const
+{
+	assert(m_structure);
+	return m_structure->get_entity_type_for_entity_id(get_entity_id());
+}
+
+// std::string residue::authInsCode() const
+// {
+// 	assert(m_structure);
+
+// 	std::string result;
+// 	if (not m_atoms.empty())
+// 		result = m_atoms.front().get_property("pdbx_PDB_ins_code");
+
+// 	return result;
+// }
+
+// std::string residue::get_auth_asym_id() const
+// {
+// 	assert(m_structure);
+
+// 	std::string result;
+// 	if (not m_atoms.empty())
+// 		result = m_atoms.front().get_property("auth_asym_id");
+
+// 	return result;
+// }
+
+// std::string residue::authSeqID() const
+// {
+// 	return m_auth_seq_id;
+// }
+
+// const Compound &residue::compound() const
+// {
+// 	auto result = compound_factory::instance().create(m_compound_id);
+// 	if (result == nullptr)
+// 		throw std::runtime_error("Failed to create compound " + m_compound_id);
+// 	return *result;
+// }
+
+// std::string residue::unique_alt_id() const
+// {
+// 	if (m_structure == nullptr)
+// 		throw std::runtime_error("Invalid residue object");
+
+// 	auto firstAlt = std::find_if(m_atoms.begin(), m_atoms.end(), [](auto &a)
+// 		{ return not a.get_label_alt_id().empty(); });
+
+// 	return firstAlt != m_atoms.end() ? firstAlt->get_label_alt_id() : "";
+// }
+
+void residue::add_atom(atom &atom)
+{
+	// atom.set_property("label_comp_id", m_compound_id);
+	// atom.set_property("label_asym_id", m_asym_id);
+	// if (m_seq_id != 0)
+	// 	atom.set_property("label_seq_id", std::to_string(m_seq_id));
+	// atom.set_property("auth_seq_id", m_auth_seq_id);
+
+	m_atoms.push_back(atom);
+}
+
+// std::vector<atom> residue::unique_atoms() const
+// {
+// 	if (m_structure == nullptr)
+// 		throw std::runtime_error("Invalid residue object");
+
+// 	std::vector<atom> result;
+// 	std::string firstAlt;
+
+// 	for (auto &atom : m_atoms)
+// 	{
+// 		auto alt = atom.get_label_alt_id();
+// 		if (alt.empty())
+// 		{
+// 			result.push_back(atom);
+// 			continue;
+// 		}
+
+// 		if (firstAlt.empty())
+// 			firstAlt = alt;
+// 		else if (alt != firstAlt)
+// 		{
+// 			if (VERBOSE > 0)
+// 				std::cerr << "skipping alternate atom " << atom << std::endl;
+// 			continue;
+// 		}
+
+// 		result.push_back(atom);
+// 	}
+
+// 	return result;
+// }
+
+// std::set<std::string> residue::getAlternateIDs() const
+// {
+// 	std::set<std::string> result;
+
+// 	for (auto a : m_atoms)
+// 	{
+// 		auto alt = a.get_label_alt_id();
+// 		if (not alt.empty())
+// 			result.insert(alt);
+// 	}
+
+// 	return result;
+// }
+
+atom residue::get_atom_by_atom_id(const std::string &atom_id) const
+{
+	atom result;
+
+	for (auto &a : m_atoms)
+	{
+		if (a.get_label_atom_id() == atom_id)
+		{
+			result = a;
+			break;
+		}
+	}
+
+	if (not result and VERBOSE > 1)
+		std::cerr << "atom with atom_id " << atom_id << " not found in residue " << m_asym_id << ':' << m_seq_id << std::endl;
+
+	return result;
+}
+
+// residue is a single entity if the atoms for the asym with m_asym_id is equal
+// to the number of atoms in this residue...  hope this is correct....
+bool residue::is_entity() const
+{
+	auto &db = m_structure->get_datablock();
+
+	auto a1 = db["atom_site"].find(key("label_asym_id") == m_asym_id);
+	//	auto a2 = atoms();
+	auto &a2 = m_atoms;
+
+	return a1.size() == a2.size();
+}
+
+// std::string residue::authID() const
+// {
+// 	return get_auth_asym_id() + authSeqID() + authInsCode();
+// }
+
+// std::string residue::labelID() const
+// {
+// 	if (m_compound_id == "HOH")
+// 		return m_asym_id + m_auth_seq_id;
+// 	else
+// 		return m_asym_id + std::to_string(m_seq_id);
+// }
+
+// std::tuple<point, float> residue::centerAndRadius() const
+// {
+// 	std::vector<point> pts;
+// 	for (auto &a : m_atoms)
+// 		pts.push_back(a.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);
+// }
+
+// bool residue::hasAlternateAtoms() const
+// {
+// 	return std::find_if(m_atoms.begin(), m_atoms.end(), [](const atom &atom)
+// 			   { return atom.is_alternate(); }) != m_atoms.end();
+// }
+
+// std::set<std::string> residue::getAtom_ids() const
+// {
+// 	std::set<std::string> ids;
+// 	for (auto a : m_atoms)
+// 		ids.insert(a.get_label_atom_id());
+
+// 	return ids;
+// }
+
+// std::vector<atom> residue::getAtomsByID(const std::string &atom_id) const
+// {
+// 	std::vector<atom> atoms;
+// 	for (auto a : m_atoms)
+// 	{
+// 		if (a.get_label_atom_id() == atom_id)
+// 			atoms.push_back(a);
+// 	}
+// 	return atoms;
+// }
+
+std::ostream &operator<<(std::ostream &os, const residue &res)
+{
+	os << res.get_compound_id() << ' ' << res.get_asym_id() << ':' << res.get_seq_id();
+
+	if (res.get_auth_asym_id() != res.get_asym_id() or res.get_auth_seq_id() != std::to_string(res.get_seq_id()))
+		os << " [" << res.get_auth_asym_id() << ':' << res.get_auth_seq_id() << ']';
+
+	return os;
+}
+
+// --------------------------------------------------------------------
+// monomer
+
+monomer::monomer(const polymer &polymer, size_t index, int seqID, const std::string &authSeqID, const std::string &compoundID)
+	: residue(*polymer.get_structure(), compoundID, polymer.get_asym_id(), seqID, authSeqID)
+	, m_polymer(&polymer)
+	, m_index(index)
+{
+}
+
+monomer::monomer(monomer &&rhs)
+	: residue(std::move(rhs))
+	, m_polymer(rhs.m_polymer)
+	, m_index(rhs.m_index)
+{
+	rhs.m_polymer = nullptr;
+}
+
+monomer &monomer::operator=(monomer &&rhs)
+{
+	residue::operator=(std::move(rhs));
+	m_polymer = rhs.m_polymer;
+	rhs.m_polymer = nullptr;
+	m_index = rhs.m_index;
+
+	return *this;
+}
+
+// bool monomer::is_first_in_chain() const
+// {
+// 	return m_index == 0;
+// }
+
+// bool monomer::is_last_in_chain() const
+// {
+// 	return m_index + 1 == m_polymer->size();
+// }
+
+// bool monomer::has_alpha() const
+// {
+// 	return m_index >= 1 and m_index + 2 < m_polymer->size();
+// }
+
+// bool monomer::has_kappa() const
+// {
+// 	return m_index >= 2 and m_index + 2 < m_polymer->size();
+// }
+
+// float monomer::phi() const
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		if (m_index > 0)
+// 		{
+// 			auto &prev = m_polymer->operator[](m_index - 1);
+// 			if (prev.m_seq_id + 1 == m_seq_id)
+// 				result = static_cast<float>(DihedralAngle(prev.C().location(), N().location(), CAlpha().location(), C().location()));
+// 		}
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::psi() const
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		if (m_index + 1 < m_polymer->size())
+// 		{
+// 			auto &next = m_polymer->operator[](m_index + 1);
+// 			if (m_seq_id + 1 == next.m_seq_id)
+// 				result = static_cast<float>(DihedralAngle(N().location(), CAlpha().location(), C().location(), next.N().location()));
+// 		}
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::alpha() const
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		if (m_index >= 1 and m_index + 2 < m_polymer->size())
+// 		{
+// 			auto &prev = m_polymer->operator[](m_index - 1);
+// 			auto &next = m_polymer->operator[](m_index + 1);
+// 			auto &nextNext = m_polymer->operator[](m_index + 2);
+
+// 			result = static_cast<float>(DihedralAngle(prev.CAlpha().location(), CAlpha().location(), next.CAlpha().location(), nextNext.CAlpha().location()));
+// 		}
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::kappa() const
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		if (m_index >= 2 and m_index + 2 < m_polymer->size())
+// 		{
+// 			auto &prevPrev = m_polymer->operator[](m_index - 2);
+// 			auto &nextNext = m_polymer->operator[](m_index + 2);
+
+// 			if (prevPrev.m_seq_id + 4 == nextNext.m_seq_id)
+// 			{
+// 				double ckap = CosinusAngle(CAlpha().location(), prevPrev.CAlpha().location(), nextNext.CAlpha().location(), CAlpha().location());
+// 				double skap = std::sqrt(1 - ckap * ckap);
+// 				result = static_cast<float>(std::atan2(skap, ckap) * 180 / kPI);
+// 			}
+// 		}
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << "When trying to calculate kappa for " << asym_id() << ':' << seqID() << ": "
+// 					  << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::tco() const
+// {
+// 	float result = 0.0;
+
+// 	try
+// 	{
+// 		if (m_index > 0)
+// 		{
+// 			auto &prev = m_polymer->operator[](m_index - 1);
+// 			if (prev.m_seq_id + 1 == m_seq_id)
+// 				result = static_cast<float>(CosinusAngle(C().location(), O().location(), prev.C().location(), prev.O().location()));
+// 		}
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << "When trying to calculate tco for " << asym_id() << ':' << seqID() << ": "
+// 					  << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::omega() const
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		if (not is_last_in_chain())
+// 			result = omega(*this, m_polymer->operator[](m_index + 1));
+// 	}
+// 	catch (const std::exception &ex)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << "When trying to calculate omega for " << asym_id() << ':' << seqID() << ": "
+// 					  << ex.what() << std::endl;
+// 	}
+
+// 	return result;
+// }
+
+// const std::map<std::string, std::vector<std::string>> kChiAtomsMap = {
+// 	{"ASP", {"CG", "OD1"}},
+// 	{"ASN", {"CG", "OD1"}},
+// 	{"ARG", {"CG", "CD", "NE", "CZ"}},
+// 	{"HIS", {"CG", "ND1"}},
+// 	{"GLN", {"CG", "CD", "OE1"}},
+// 	{"GLU", {"CG", "CD", "OE1"}},
+// 	{"SER", {"OG"}},
+// 	{"THR", {"OG1"}},
+// 	{"LYS", {"CG", "CD", "CE", "NZ"}},
+// 	{"TYR", {"CG", "CD1"}},
+// 	{"PHE", {"CG", "CD1"}},
+// 	{"LEU", {"CG", "CD1"}},
+// 	{"TRP", {"CG", "CD1"}},
+// 	{"CYS", {"SG"}},
+// 	{"ILE", {"CG1", "CD1"}},
+// 	{"MET", {"CG", "SD", "CE"}},
+// 	{"MSE", {"CG", "SE", "CE"}},
+// 	{"PRO", {"CG", "CD"}},
+// 	{"VAL", {"CG1"}}};
+
+// size_t monomer::nrOfChis() const
+// {
+// 	size_t result = 0;
+
+// 	auto i = kChiAtomsMap.find(m_compound_id);
+// 	if (i != kChiAtomsMap.end())
+// 		result = i->second.size();
+
+// 	return result;
+// }
+
+// float monomer::chi(size_t nr) const
+// {
+// 	float result = 0;
+
+// 	try
+// 	{
+// 		auto i = kChiAtomsMap.find(m_compound_id);
+// 		if (i != kChiAtomsMap.end() and nr < i->second.size())
+// 		{
+// 			std::vector<std::string> atoms{"N", "CA", "CB"};
+
+// 			atoms.insert(atoms.end(), i->second.begin(), i->second.end());
+
+// 			// in case we have a positive chiral volume we need to swap atoms
+// 			if (chiralVolume() > 0)
+// 			{
+// 				if (m_compound_id == "LEU")
+// 					atoms.back() = "CD2";
+// 				if (m_compound_id == "VAL")
+// 					atoms.back() = "CG2";
+// 			}
+
+// 			result = static_cast<float>(DihedralAngle(
+// 				get_atom_by_id(atoms[nr + 0]).location(),
+// 				get_atom_by_id(atoms[nr + 1]).location(),
+// 				get_atom_by_id(atoms[nr + 2]).location(),
+// 				get_atom_by_id(atoms[nr + 3]).location()));
+// 		}
+// 	}
+// 	catch (const std::exception &e)
+// 	{
+// 		if (VERBOSE > 0)
+// 			std::cerr << e.what() << std::endl;
+// 		result = 0;
+// 	}
+
+// 	return result;
+// }
+
+// bool monomer::isCis() const
+// {
+// 	bool result = false;
+
+// 	if (m_index + 1 < m_polymer->size())
+// 	{
+// 		auto &next = m_polymer->operator[](m_index + 1);
+
+// 		result = monomer::isCis(*this, next);
+// 	}
+
+// 	return result;
+// }
+
+// bool monomer::isComplete() const
+// {
+// 	int seen = 0;
+// 	for (auto &a : m_atoms)
+// 	{
+// 		if (a.get_label_atom_id() == "CA")
+// 			seen |= 1;
+// 		else if (a.get_label_atom_id() == "C")
+// 			seen |= 2;
+// 		else if (a.get_label_atom_id() == "N")
+// 			seen |= 4;
+// 		else if (a.get_label_atom_id() == "O")
+// 			seen |= 8;
+// 		// else if (a.get_label_atom_id() == "OXT")		seen |= 16;
+// 	}
+// 	return seen == 15;
+// }
+
+// bool monomer::hasAlternateBackboneAtoms() const
+// {
+// 	bool result = false;
+
+// 	for (auto &a : m_atoms)
+// 	{
+// 		if (not a.is_alternate())
+// 			continue;
+
+// 		auto atom_id = a.get_label_atom_id();
+// 		if (atom_id == "CA" or atom_id == "C" or atom_id == "N" or atom_id == "O")
+// 		{
+// 			result = true;
+// 			break;
+// 		}
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::chiralVolume() const
+// {
+// 	float result = 0;
+
+// 	if (m_compound_id == "LEU")
+// 	{
+// 		auto centre = get_atom_by_id("CG");
+// 		auto atom1 = get_atom_by_id("CB");
+// 		auto atom2 = get_atom_by_id("CD1");
+// 		auto atom3 = get_atom_by_id("CD2");
+
+// 		result = DotProduct(atom1.location() - centre.location(),
+// 			CrossProduct(atom2.location() - centre.location(), atom3.location() - centre.location()));
+// 	}
+// 	else if (m_compound_id == "VAL")
+// 	{
+// 		auto centre = get_atom_by_id("CB");
+// 		auto atom1 = get_atom_by_id("CA");
+// 		auto atom2 = get_atom_by_id("CG1");
+// 		auto atom3 = get_atom_by_id("CG2");
+
+// 		result = DotProduct(atom1.location() - centre.location(),
+// 			CrossProduct(atom2.location() - centre.location(), atom3.location() - centre.location()));
+// 	}
+
+// 	return result;
+// }
+
+// bool monomer::areBonded(const monomer &a, const monomer &b, float errorMargin)
+// {
+// 	bool result = false;
+
+// 	try
+// 	{
+// 		point atoms[4] = {
+// 			a.get_atom_by_id("CA").location(),
+// 			a.get_atom_by_id("C").location(),
+// 			b.get_atom_by_id("N").location(),
+// 			b.get_atom_by_id("CA").location()};
+
+// 		auto distanceCACA = Distance(atoms[0], atoms[3]);
+// 		double omega = DihedralAngle(atoms[0], atoms[1], atoms[2], atoms[3]);
+
+// 		bool cis = std::abs(omega) <= 30.0;
+// 		float maxCACADistance = cis ? 3.0f : 3.8f;
+
+// 		result = std::abs(distanceCACA - maxCACADistance) < errorMargin;
+// 	}
+// 	catch (...)
+// 	{
+// 	}
+
+// 	return result;
+// }
+
+// float monomer::omega(const mmmonomer &a, const mmmonomer &b)
+// {
+// 	float result = 360;
+
+// 	try
+// 	{
+// 		result = static_cast<float>(DihedralAngle(
+// 			a.get_atom_by_id("CA").location(),
+// 			a.get_atom_by_id("C").location(),
+// 			b.get_atom_by_id("N").location(),
+// 			b.get_atom_by_id("CA").location()));
+// 	}
+// 	catch (...)
+// 	{
+// 	}
+
+// 	return result;
+// }
+
+// bool monomer::isCis(const mmmonomer &a, const mmmonomer &b)
+// {
+// 	return omega(a, b) < 30.0f;
+// }
+
+// --------------------------------------------------------------------
+// polymer
+
+polymer::polymer(const structure &s, const std::string &entityID, const std::string &asym_id)
+	: m_structure(const_cast<structure *>(&s))
+	, m_entity_id(entityID)
+	, m_asym_id(asym_id)
+	// , mPolySeq(s.m_db["pdbx_poly_seq_scheme"), key("asym_id") == m_asym_id and key("entity_id") == m_entity_i])
+{
+	std::map<size_t, size_t> ix;
+
+	auto &poly_seq_scheme = s.get_datablock()["pdbx_poly_seq_scheme"];
+	reserve(poly_seq_scheme.size());
+
+	for (auto r : poly_seq_scheme)
+	{
+		int seqID;
+		std::string compoundID, authSeqID;
+		cif::tie(seqID, authSeqID, compoundID) = r.get("seq_id", "auth_seq_num", "mon_id");
+
+		size_t index = size();
+
+		// store only the first
+		if (not ix.count(seqID))
+		{
+			ix[seqID] = index;
+			emplace_back(*this, index, seqID, authSeqID, compoundID);
+		}
+		else if (VERBOSE > 0)
+		{
+			monomer m{*this, index, seqID, authSeqID, compoundID};
+			std::cerr << "Dropping alternate residue " << m << std::endl;
+		}
+	}
+}
+
+// std::string polymer::chainID() const
+// {
+// 	return mPolySeq.front()["pdb_strand_id"].as<std::string>();
+// }
+
+// monomer &polymer::getBySeqID(int seqID)
+// {
+// 	for (auto &m : *this)
+// 		if (m.get_seq_id() == seqID)
+// 			return m;
+
+// 	throw std::runtime_error("monomer with seqID " + std::to_string(seqID) + " not found in polymer " + m_asym_id);
+// }
+
+// const monomer &polymer::getBySeqID(int seqID) const
+// {
+// 	for (auto &m : *this)
+// 		if (m.get_seq_id() == seqID)
+// 			return m;
+
+// 	throw std::runtime_error("monomer with seqID " + std::to_string(seqID) + " not found in polymer " + m_asym_id);
+// }
+
+// int polymer::Distance(const monomer &a, const monomer &b) const
+// {
+// 	int result = std::numeric_limits<int>::max();
+
+// 	if (a.get_asym_id() == b.get_asym_id())
+// 	{
+// 		int ixa = std::numeric_limits<int>::max(), ixb = std::numeric_limits<int>::max();
+
+// 		int ix = 0, f = 0;
+// 		for (auto &m : *this)
+// 		{
+// 			if (m.get_seq_id() == a.get_seq_id())
+// 				ixa = ix, ++f;
+// 			if (m.get_seq_id() == b.get_seq_id())
+// 				ixb = ix, ++f;
+// 			if (f == 2)
+// 			{
+// 				result = std::abs(ixa - ixb);
+// 				break;
+// 			}
+// 		}
+// 	}
+
+// 	return result;
+// }
+
+// --------------------------------------------------------------------
+
+sugar::sugar(const branch &branch, const std::string &compoundID,
+	const std::string &asym_id, int authSeqID)
+	: residue(branch.get_structure(), compoundID, asym_id, 0, std::to_string(authSeqID))
+	, m_branch(&branch)
+{
+}
+
+sugar::sugar(sugar &&rhs)
+	: residue(std::forward<residue>(rhs))
+	, m_branch(rhs.m_branch)
+{
+
+}
+
+sugar &sugar::operator=(sugar &&rhs)
+{
+	if (this != &rhs)
+	{
+		residue::operator=(std::forward<residue>(rhs));
+		m_branch = rhs.m_branch;
+	}
+
+	return *this;
+}
+
+// bool sugar::hasLinkedSugarAtLeavingO(int leavingO) const
+// {
+// 	return false;
+// }
+
+// sugar &sugar::operator[](int leavingO)
+// {
+// 	throw std::logic_error("not implemented");
+// }
+
+// const sugar &sugar::operator[](int leavingO) const
+// {
+// 	throw std::logic_error("not implemented");
+// }
+
+std::string sugar::name() const
+{
+	std::string result;
+
+	if (m_compound_id == "MAN")
+		result += "alpha-D-mannopyranose";
+	else if (m_compound_id == "BMA")
+		result += "beta-D-mannopyranose";
+	else if (m_compound_id == "NAG")
+		result += "2-acetamido-2-deoxy-beta-D-glucopyranose";
+	else if (m_compound_id == "NDG")
+		result += "2-acetamido-2-deoxy-alpha-D-glucopyranose";
+	else if (m_compound_id == "FUC")
+		result += "alpha-L-fucopyranose";
+	else if (m_compound_id == "FUL")
+		result += "beta-L-fucopyranose";
+	else
+	{
+		auto compound = compound_factory::instance().create(m_compound_id);
+		if (compound)
+			result += compound->name();
+		else
+			result += m_compound_id;
+	}
+
+	return result;
+}
+
+branch::branch(structure &structure, const std::string &asym_id)
+	: m_structure(&structure)
+	, m_asym_id(asym_id)
+{
+	using namespace literals;
+
+	auto &db = structure.get_datablock();
+	auto &struct_asym = db["struct_asym"];
+	auto &branch_scheme = db["pdbx_branch_scheme"];
+	auto &branch_link = db["pdbx_entity_branch_link"];
+
+	for (const auto &[entity_id] : struct_asym.find<std::string>("id"_key == asym_id, "entity_id"))
+	{
+		for (const auto &[comp_id, num] : branch_scheme.find<std::string, int>(
+				 "asym_id"_key == asym_id, "mon_id", "pdb_seq_num"))
+		{
+			emplace_back(*this, comp_id, asym_id, num);
+		}
+
+		for (const auto &[num1, num2, atom1, atom2] : branch_link.find<size_t, size_t, std::string, std::string>(
+				 "entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"))
+		{
+			if (not iequals(atom1, "c1"))
+				throw std::runtime_error("invalid pdbx_entity_branch_link");
+
+			auto &s1 = at(num1 - 1);
+			auto &s2 = at(num2 - 1);
+
+			s1.set_link(s2.get_atom_by_atom_id(atom2));
+		}
+
+		break;
+	}
+}
+
+void branch::link_atoms()
+{
+	using namespace literals;
+
+	auto &db = m_structure->get_datablock();
+	auto &branch_link = db["pdbx_entity_branch_link"];
+
+	auto entity_id = front().get_entity_id();
+
+	for (const auto &[num1, num2, atom1, atom2] : branch_link.find<size_t, size_t, std::string, std::string>(
+			 "entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"))
+	{
+		if (not iequals(atom1, "c1"))
+			throw std::runtime_error("invalid pdbx_entity_branch_link");
+
+		auto &s1 = at(num1 - 1);
+		auto &s2 = at(num2 - 1);
+
+		s1.set_link(s2.get_atom_by_atom_id(atom2));
+	}
+}
+
+std::string branch::name() const
+{
+	return empty() ? "" : name(front());
+}
+
+std::string branch::name(const sugar &s) const
+{
+	using namespace literals;
+
+	std::string result;
+
+	for (auto &sn : *this)
+	{
+		if (not sn.get_link() or sn.get_link().get_auth_seq_id() != s.get_auth_seq_id())
+			continue;
+
+		auto n = name(sn) + "-(1-" + sn.get_link().get_label_atom_id().substr(1) + ')';
+
+		result = result.empty() ? n : result + "-[" + n + ']';
+	}
+
+	if (not result.empty() and result.back() != ']')
+		result += '-';
+
+	return result + s.name();
+}
+
+float branch::weight() const
+{
+	return std::accumulate(begin(), end(), 0.f, [](float sum, const sugar &s)
+		{
+		auto compound = compound_factory::instance().create(s.get_compound_id());
+		if (compound)
+			sum += compound->formula_weight();
+		return sum; });
+}
+
+// // --------------------------------------------------------------------
+// // File
+
+// void File::load(const std::filesystem::path &path)
+// {
+// 	gxrio::ifstream in(path);
+
+// 	auto ext = path.extension().string();
+// 	if (ext == ".gz" or ext = ".xz")
+// 		ext = path.stem().extension().string();
+
+// 	if (ext == ".pdb" or ext == ".ent")
+// 		ReadPDBFile(in, *this);
+// 	else
+// 		file::load(in);
+
+// 	// validate, otherwise lots of functionality won't work
+// 	loadDictionary("mmcif_pdbx");
+// 	if (not isValid() and VERBOSE >= 0)
+// 		std::cerr << "Invalid mmCIF file" << (VERBOSE > 0 ? "." : " use --verbose option to see errors") << std::endl;
+// }
+
+// void File::save(const std::filesystem::path &path)
+// {
+// 	gxrio::ostream outFile(path);
+
+// 	auto ext = path.extension().string();
+// 	if (ext == ".gz" or ext = ".xz")
+// 		ext = path.stem().extension().string();
+
+// 	if (ext == ".pdb" or ext == ".ent")
+// 		WritePDBFile(outFile, data());
+// 	else
+// 		file::save(outFile);
+// }
+
+// --------------------------------------------------------------------
+//	structure
+
+structure::structure(file &p, size_t modelNr, StructureOpenOptions options)
+	: structure(p.front(), modelNr, options)
+{
+}
+
+structure::structure(datablock &db, size_t modelNr, StructureOpenOptions options)
+	: m_db(db)
+	, m_model_nr(modelNr)
+{
+	auto &atomCat = db["atom_site"];
+
+	load_atoms_for_model(options);
+
+	// Check to see if we should actually load another model?
+	if (m_atoms.empty() and m_model_nr == 1)
+	{
+		std::optional<size_t> model_nr;
+		cif::tie(model_nr) = atomCat.front().get("pdbx_PDB_model_num");
+		if (model_nr and *model_nr != m_model_nr)
+		{
+			if (VERBOSE > 0)
+				std::cerr << "No atoms loaded for model 1, trying model " << *model_nr << std::endl;
+			m_model_nr = *model_nr;
+			load_atoms_for_model(options);
+		}
+	}
+
+	if (m_atoms.empty())
+	{
+		if (VERBOSE >= 0)
+			std::cerr << "Warning: no atoms loaded" << std::endl;
+	}
+	else
+		loadData();
+}
+
+void structure::load_atoms_for_model(StructureOpenOptions options)
+{
+	auto &atomCat = m_db["atom_site"];
+
+	for (const auto &a : atomCat)
+	{
+		std::string id, type_symbol;
+		std::optional<size_t> model_nr;
+
+		cif::tie(id, type_symbol, model_nr) = a.get("id", "type_symbol", "pdbx_PDB_model_num");
+
+		if (model_nr and *model_nr != m_model_nr)
+			continue;
+
+		if ((options bitand StructureOpenOptions::SkipHydrogen) and type_symbol == "H")
+			continue;
+
+		emplace_atom(std::make_shared<atom::atom_impl>(m_db, id, a));
+	}
+}
+
+// structure::structure(const structure &s)
+// 	: m_db(s.m_db)
+// 	, m_model_nr(s.m_model_nr)
+// {
+// 	m_atoms.reserve(s.m_atoms.size());
+// 	for (auto &atom : s.m_atoms)
+// 		emplace_atom(atom.clone());
+
+// 	loadData();
+// }
+
+// structure::~structure()
+// {
+// }
+
+void structure::loadData()
+{
+	auto &polySeqScheme = m_db["pdbx_poly_seq_scheme"];
+
+	for (const auto &[asym_id, entityID] : polySeqScheme.rows<std::string,std::string>("asym_id", "entity_id"))
+	{
+		if (m_polymers.empty() or m_polymers.back().get_asym_id() != asym_id or m_polymers.back().get_entity_id() != entityID)
+			m_polymers.emplace_back(*this, entityID, asym_id);
+	}
+
+	auto &branchScheme = m_db["pdbx_branch_scheme"];
+
+	for (const auto &[asym_id] : branchScheme.rows<std::string>("asym_id"))
+	{
+		if (m_branches.empty() or m_branches.back().get_asym_id() != asym_id)
+			m_branches.emplace_back(*this, asym_id);
+	}
+
+	auto &nonPolyScheme = m_db["pdbx_nonpoly_scheme"];
+
+	for (const auto&[asym_id, monID, pdbSeqNum] : nonPolyScheme.rows<std::string,std::string,std::string>("asym_id", "mon_id", "pdb_seq_num"))
+		m_non_polymers.emplace_back(*this, monID, asym_id, 0, pdbSeqNum);
+
+	// place atoms in residues
+
+	using key_type = std::tuple<std::string, int, std::string>;
+	std::map<key_type, residue *> resMap;
+
+	for (auto &poly : m_polymers)
+	{
+		for (auto &res : poly)
+			resMap[{res.get_asym_id(), res.get_seq_id(), res.get_auth_seq_id()}] = &res;
+	}
+
+	for (auto &res : m_non_polymers)
+		resMap[{res.get_asym_id(), res.get_seq_id(), res.get_auth_seq_id()}] = &res;
+
+	std::set<std::string> sugars;
+	for (auto &branch : m_branches)
+	{
+		for (auto &sugar : branch)
+		{
+			resMap[{sugar.get_asym_id(), sugar.get_seq_id(), sugar.get_auth_seq_id()}] = &sugar;
+			sugars.insert(sugar.get_compound_id());
+		}
+	}
+
+	for (auto &atom : m_atoms)
+	{
+		key_type k(atom.get_label_asym_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
+		auto ri = resMap.find(k);
+
+		if (ri == resMap.end())
+		{
+			if (VERBOSE > 0)
+				std::cerr << "Missing residue for atom " << atom << std::endl;
+
+			// see if it might match a non poly
+			for (auto &res : m_non_polymers)
+			{
+				if (res.get_asym_id() != atom.get_label_asym_id())
+					continue;
+
+				res.add_atom(atom);
+				break;
+			}
+
+			continue;
+		}
+
+		ri->second->add_atom(atom);
+	}
+
+	for (auto &branch : m_branches)
+		branch.link_atoms();
+}
+
+EntityType structure::get_entity_type_for_entity_id(const std::string entityID) const
+{
+	using namespace literals;
+
+	auto &entity = m_db["entity"];
+	auto entity_type = entity.find1<std::string>("id"_key == entityID, "type");
+
+	EntityType result;
+
+	if (iequals(entity_type, "polymer"))
+		result = EntityType::polymer;
+	else if (iequals(entity_type, "non-polymer"))
+		result = EntityType::NonPolymer;
+	else if (iequals(entity_type, "macrolide"))
+		result = EntityType::Macrolide;
+	else if (iequals(entity_type, "water"))
+		result = EntityType::Water;
+	else if (iequals(entity_type, "branched"))
+		result = EntityType::Branched;
+	else
+		throw std::runtime_error("Unknown entity type " + entity_type);
+
+	return result;
+}
+
+EntityType structure::get_entity_type_for_asym_id(const std::string asym_id) const
+{
+	using namespace literals;
+
+	auto &struct_asym = m_db["struct_asym"];
+	auto entityID = struct_asym.find1<std::string>("id"_key == asym_id, "entity_id");
+
+	return get_entity_type_for_entity_id(entityID);
+}
+
+// std::vector<atom> structure::waters() const
+// {
+// 	using namespace literals;
+
+// 	std::vector<atom> result;
+
+// 	auto &db = datablock();
+
+// 	// Get the entity id for water. Watch out, structure may not have water at all
+// 	auto &entityCat = db["entity"];
+// 	for (const auto &[waterEntityID] : entityCat.find<std::string>("type"_key == "water", "id"))
+// 	{
+// 		for (auto &a : m_atoms)
+// 		{
+// 			if (a.get_property("label_entity_id") == waterEntityID)
+// 				result.push_back(a);
+// 		}
+
+// 		break;
+// 	}
+
+// 	return result;
+// }
+
+atom structure::get_atom_by_id(const std::string &id) const
+{
+	assert(m_atoms.size() == m_atom_index.size());
+
+	int L = 0, R = m_atoms.size() - 1;
+	while (L <= R)
+	{
+		int i = (L + R) / 2;
+
+		const atom &atom = m_atoms[m_atom_index[i]];
+
+		int d = atom.id().compare(id);
+
+		if (d == 0)
+			return atom;
+
+		if (d < 0)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	throw std::out_of_range("Could not find atom with id " + id);
+}
+
+// atom structure::getAtomByLabel(const std::string &atom_id, const std::string &asym_id, const std::string &compID, int seqID, const std::string &altID)
+// {
+// 	for (auto &a : m_atoms)
+// 	{
+// 		if (a.get_label_atom_id() == atom_id and
+// 			a.get_label_asym_id() == asym_id and
+// 			a.labelCompID() == compID and
+// 			a.get_label_seq_id() == seqID and
+// 			a.get_label_alt_id() == altID)
+// 		{
+// 			return a;
+// 		}
+// 	}
+
+// 	throw std::out_of_range("Could not find atom with specified label");
+// }
+
+// atom structure::getAtomByPosition(point p) const
+// {
+// 	double distance = std::numeric_limits<double>::max();
+// 	size_t index = std::numeric_limits<size_t>::max();
+
+// 	for (size_t i = 0; i < m_atoms.size(); ++i)
+// 	{
+// 		auto &a = m_atoms.at(i);
+
+// 		auto d = Distance(a.location(), p);
+// 		if (d < distance)
+// 		{
+// 			distance = d;
+// 			index = i;
+// 		}
+// 	}
+
+// 	if (index < m_atoms.size())
+// 		return m_atoms.at(index);
+
+// 	return {};
+// }
+
+// atom structure::getAtomByPositionAndType(point p, std::string_view type, std::string_view res_type) const
+// {
+// 	double distance = std::numeric_limits<double>::max();
+// 	size_t index = std::numeric_limits<size_t>::max();
+
+// 	for (size_t i = 0; i < m_atoms.size(); ++i)
+// 	{
+// 		auto &a = m_atoms.at(i);
+
+// 		if (a.labelCompID() != res_type)
+// 			continue;
+
+// 		if (a.get_label_atom_id() != type)
+// 			continue;
+
+// 		auto d = Distance(a.location(), p);
+// 		if (d < distance)
+// 		{
+// 			distance = d;
+// 			index = i;
+// 		}
+// 	}
+
+// 	if (index < m_atoms.size())
+// 		return m_atoms.at(index);
+
+// 	return {};
+// }
+
+// polymer &structure::getPolymerByAsym_id(const std::string &asym_id)
+// {
+// 	for (auto &poly : m_polymers)
+// 	{
+// 		if (poly.get_asym_id() != asym_id)
+// 			continue;
+
+// 		return poly;
+// 	}
+
+// 	throw std::runtime_error("polymer with asym id " + asym_id + " not found");
+// }
+
+residue &structure::get_residue(const std::string &asym_id, int seqID, const std::string &authSeqID)
+{
+	if (seqID == 0)
+	{
+		for (auto &res : m_non_polymers)
+		{
+			if (res.get_asym_id() == asym_id and (authSeqID.empty() or res.get_auth_seq_id() == authSeqID))
+				return res;
+		}
+	}
+
+	for (auto &poly : m_polymers)
+	{
+		if (poly.get_asym_id() != asym_id)
+			continue;
+
+		for (auto &res : poly)
+		{
+			if (res.get_seq_id() == seqID)
+				return res;
+		}
+	}
+
+	for (auto &branch : m_branches)
+	{
+		if (branch.get_asym_id() != asym_id)
+			continue;
+
+		for (auto &sugar : branch)
+		{
+			if (sugar.get_asym_id() == asym_id and sugar.get_auth_seq_id() == authSeqID)
+				return sugar;
+		}
+	}
+
+	std::string desc = asym_id;
+
+	if (seqID != 0)
+		desc += "/" + std::to_string(seqID);
+	
+	if (not authSeqID.empty())
+		desc += "-" + authSeqID;
+
+	throw std::out_of_range("Could not find residue " + desc);
+}
+
+residue &structure::get_residue(const std::string &asym_id, const std::string &compID, int seqID, const std::string &authSeqID)
+{
+	if (seqID == 0)
+	{
+		for (auto &res : m_non_polymers)
+		{
+			if (res.get_asym_id() == asym_id and res.get_auth_seq_id() == authSeqID and res.get_compound_id() == compID)
+				return res;
+		}
+	}
+
+	for (auto &poly : m_polymers)
+	{
+		if (poly.get_asym_id() != asym_id)
+			continue;
+
+		for (auto &res : poly)
+		{
+			if (res.get_seq_id() == seqID and res.get_compound_id() == compID)
+				return res;
+		}
+	}
+
+	for (auto &branch : m_branches)
+	{
+		if (branch.get_asym_id() != asym_id)
+			continue;
+
+		for (auto &sugar : branch)
+		{
+			if (sugar.get_asym_id() == asym_id and sugar.get_auth_seq_id() == authSeqID and sugar.get_compound_id() == compID)
+				return sugar;
+		}
+	}
+
+	std::string desc = asym_id;
+
+	if (seqID != 0)
+		desc += "/" + std::to_string(seqID);
+	
+	if (not authSeqID.empty())
+		desc += "-" + authSeqID;
+
+	throw std::out_of_range("Could not find residue " + desc + " of type " + compID);
+}
+
+branch &structure::get_branch_by_asym_id(const std::string &asym_id)
+{
+	for (auto &branch : m_branches)
+	{
+		if (branch.get_asym_id() == asym_id)
+			return branch;
+	}
+
+	throw std::runtime_error("branch not found for asym id " + asym_id);
+}
+
+std::string structure::insert_compound(const std::string &compoundID, bool is_entity)
+{
+	using namespace literals;
+
+	auto compound = compound_factory::instance().create(compoundID);
+	if (compound == nullptr)
+		throw std::runtime_error("Trying to insert unknown compound " + compoundID + " (not found in CCD)");
+
+	auto &chemComp = m_db["chem_comp"];
+	auto r = chemComp.find(key("id") == compoundID);
+	if (r.empty())
+	{
+		chemComp.emplace({
+			{"id", compoundID},
+			{"name", compound->name()},
+			{"formula", compound->formula()},
+			{"formula_weight", compound->formula_weight()},
+			{"type", compound->type()}});
+	}
+
+	std::string entity_id;
+
+	if (is_entity)
+	{
+		auto &pdbxEntityNonpoly = m_db["pdbx_entity_nonpoly"];
+
+		entity_id = pdbxEntityNonpoly.find1<std::string>("comp_id"_key == compoundID, "entity_id");
+
+		if (entity_id.empty())
+		{
+			auto &entity = m_db["entity"];
+			entity_id = entity.get_unique_id("");
+
+			entity.emplace({
+				{"id", entity_id},
+				{"type", "non-polymer"},
+				{"pdbx_description", compound->name()},
+				{"formula_weight", compound->formula_weight()}});
+
+			pdbxEntityNonpoly.emplace({
+				{"entity_id", entity_id},
+				{"name", compound->name()},
+				{"comp_id", compoundID}});
+		}
+	}
+
+	return entity_id;
+}
+
+// --------------------------------------------------------------------
+
+atom &structure::emplace_atom(atom &&atom)
+{
+	int L = 0, R = m_atom_index.size() - 1;
+	while (L <= R)
+	{
+		int i = (L + R) / 2;
+
+		auto &ai = m_atoms[m_atom_index[i]];
+
+		int d = ai.id().compare(atom.id());
+
+		if (d == 0)
+			throw std::runtime_error("Duplicate atom ID " + atom.id());
+
+		if (d < 0)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	m_atom_index.insert(m_atom_index.begin() + R + 1, m_atoms.size());
+
+	return m_atoms.emplace_back(std::move(atom));
+}
+
+void structure::remove_atom(atom &a, bool removeFromResidue)
+{
+	using namespace literals;
+
+	auto &atomSites = m_db["atom_site"];
+	atomSites.erase("id"_key == a.id());
+
+	if (removeFromResidue)
+	{
+		try
+		{
+			auto &res = get_residue(a);
+			res.m_atoms.erase(std::remove(res.m_atoms.begin(), res.m_atoms.end(), a), res.m_atoms.end());
+		}
+		catch (const std::exception &ex)
+		{
+			if (VERBOSE > 0)
+				std::cerr << "Error removing atom from residue: " << ex.what() << std::endl;
+		}
+	}
+
+	assert(m_atom_index.size() == m_atoms.size());
+
+#ifndef NDEBUG
+	bool removed = false;
+#endif
+
+	int L = 0, R = m_atom_index.size() - 1;
+	while (L <= R)
+	{
+		int i = (L + R) / 2;
+
+		const atom &atom = m_atoms[m_atom_index[i]];
+
+		int d = atom.id().compare(a.id());
+
+		if (d == 0)
+		{
+			m_atoms.erase(m_atoms.begin() + m_atom_index[i]);
+
+			auto ai = m_atom_index[i];
+			m_atom_index.erase(m_atom_index.begin() + i);
+
+			for (auto &j : m_atom_index)
+			{
+				if (j > ai)
+					--j;
+			}
+#ifndef NDEBUG
+			removed = true;
+#endif
+			break;
+		}
+
+		if (d < 0)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+#ifndef NDEBUG
+	assert(removed);
+#endif
+}
+
+void structure::swap_atoms(atom a1, atom a2)
+{
+	auto &atomSites = m_db["atom_site"];
+
+	try
+	{
+		auto r1 = atomSites.find1(key("id") == a1.id());
+		auto r2 = atomSites.find1(key("id") == a2.id());
+
+		auto l1 = r1["label_atom_id"];
+		auto l2 = r2["label_atom_id"];
+		l1.swap(l2);
+
+		auto l3 = r1["auth_atom_id"];
+		auto l4 = r2["auth_atom_id"];
+		l3.swap(l4);
+	}
+	catch (const std::exception &ex)
+	{
+		std::throw_with_nested(std::runtime_error("Failed to swap atoms"));
+	}
+}
+
+void structure::move_atom(atom a, point p)
+{
+	a.set_location(p);
+}
+
+void structure::change_residue(residue &res, const std::string &newCompound,
+	const std::vector<std::tuple<std::string, std::string>> &remappedAtoms)
+{
+	using namespace literals;
+
+	std::string asym_id = res.get_asym_id();
+
+	const auto compound = compound_factory::instance().create(newCompound);
+	if (not compound)
+		throw std::runtime_error("Unknown compound " + newCompound);
+
+	// First make sure the compound is already known or insert it.
+	// And if the residue is an entity, we must make sure it exists
+
+	std::string entityID;
+	if (res.is_entity())
+	{
+		// create a copy of the entity first
+		auto &entity = m_db["entity"];
+
+		try
+		{
+			entityID = entity.find1<std::string>("type"_key == "non-polymer" and "pdbx_description"_key == compound->name(), "id");
+		}
+		catch (const std::exception &ex)
+		{
+			entityID = entity.get_unique_id("");
+			entity.emplace({{"id", entityID},
+				{"type", "non-polymer"},
+				{"pdbx_description", compound->name()},
+				{"formula_weight", compound->formula_weight()}});
+		}
+
+		auto &pdbxEntityNonpoly = m_db["pdbx_entity_nonpoly"];
+		pdbxEntityNonpoly.emplace({{"entity_id", entityID},
+			{"name", compound->name()},
+			{"comp_id", newCompound}});
+
+		auto &pdbxNonPolyScheme = m_db["pdbx_nonpoly_scheme"];
+		for (auto nps : pdbxNonPolyScheme.find("asym_id"_key == asym_id))
+		{
+			nps.assign("mon_id", newCompound, true);
+			nps.assign("auth_mon_id", newCompound, true);
+			nps.assign("entity_id", entityID, true);
+		}
+
+		// create rest
+		auto &chemComp = m_db["chem_comp"];
+		if (not chemComp.exists(key("id") == newCompound))
+		{
+			chemComp.emplace({{"id", newCompound},
+				{"name", compound->name()},
+				{"formula", compound->formula()},
+				{"formula_weight", compound->formula_weight()},
+				{"type", compound->type()}});
+		}
+
+		// update the struct_asym for the new entity
+		m_db["struct_asym"].update_value("id"_key == asym_id, "entity_id", entityID);
+	}
+	else
+		insert_compound(newCompound, false);
+
+	res.set_compound_id(newCompound);
+
+	auto &atomSites = m_db["atom_site"];
+	auto atoms = res.atoms();
+
+	for (const auto &[a1, a2] : remappedAtoms)
+	{
+		auto i = find_if(atoms.begin(), atoms.end(), [id = a1](const atom &a)
+			{ return a.get_label_atom_id() == id; });
+		if (i == atoms.end())
+		{
+			if (VERBOSE >= 0)
+				std::cerr << "Missing atom for atom ID " << a1 << std::endl;
+			continue;
+		}
+
+		auto r = atomSites.find(key("id") == i->id());
+
+		if (r.size() != 1)
+			continue;
+
+		if (a2.empty() or a2 == ".")
+			remove_atom(*i);
+		else if (a1 != a2)
+		{
+			auto ra = r.front();
+			ra["label_atom_id"] = a2;
+			ra["auth_atom_id"] = a2;
+			ra["type_symbol"] = atom_type_traits(compound->get_atom_by_atom_id(a2).type_symbol).symbol();
+		}
+	}
+
+	for (auto a : atoms)
+	{
+		atomSites.update_value(key("id") == a.id(), "label_comp_id", newCompound);
+		atomSites.update_value(key("id") == a.id(), "auth_comp_id", newCompound);
+	}
+}
+
+void structure::remove_residue(residue &res)
+{
+	using namespace literals;
+
+	auto atoms = res.atoms();
+
+	switch (res.entity_type())
+	{
+		case EntityType::polymer:
+		{
+			auto &m = dynamic_cast<monomer &>(res);
+
+			m_db["pdbx_poly_seq_scheme"].erase(
+				"asym_id"_key == res.get_asym_id() and
+				"seq_id"_key == res.get_seq_id());
+
+			for (auto &poly : m_polymers)
+				poly.erase(std::remove(poly.begin(), poly.end(), m), poly.end());
+			break;
+		}
+
+		case EntityType::NonPolymer:
+			m_db["pdbx_nonpoly_scheme"].erase("asym_id"_key == res.get_asym_id());
+			m_db["struct_asym"].erase("id"_key == res.get_asym_id());
+			m_non_polymers.erase(std::remove(m_non_polymers.begin(), m_non_polymers.end(), res), m_non_polymers.end());
+			break;
+
+		case EntityType::Water:
+			m_db["pdbx_nonpoly_scheme"].erase("asym_id"_key == res.get_asym_id());
+			m_non_polymers.erase(std::remove(m_non_polymers.begin(), m_non_polymers.end(), res), m_non_polymers.end());
+			break;
+
+		case EntityType::Branched:
+		{
+			auto &s = dynamic_cast<sugar&>(res);
+
+			remove_sugar(s);
+
+			atoms.clear();
+			break;
+		}
+
+		case EntityType::Macrolide:
+			// TODO: Fix this?
+			throw std::runtime_error("no support for macrolides yet");
+	}
+
+	for (auto atom : atoms)
+		remove_atom(atom, false);
+}
+
+void structure::remove_sugar(sugar &s)
+{
+	using namespace literals;
+
+	std::string asym_id = s.get_asym_id();
+	branch &branch = get_branch_by_asym_id(asym_id);
+	auto si = std::find(branch.begin(), branch.end(), s);
+	if (si == branch.end())
+		throw std::runtime_error("sugar not part of branch");
+	size_t six = si - branch.begin();
+
+	if (six == 0)	// first sugar, means the death of this branch
+		remove_branch(branch);
+	else
+	{
+		std::set<size_t> dix;
+		std::stack<size_t> test;
+		test.push(s.num());
+
+		while (not test.empty())
+		{
+			auto tix = test.top();
+			test.pop();
+
+			if (dix.count(tix))
+				continue;
+			
+			dix.insert(tix);
+
+			for (auto atom : branch[tix - 1].atoms())
+				remove_atom(atom, false);
+
+			for (auto &s : branch)
+			{
+				if (s.get_link_nr() == tix)
+					test.push(s.num());
+			}
+		}
+
+		branch.erase(remove_if(branch.begin(), branch.end(), [dix](const sugar &s) { return dix.count(s.num()); }), branch.end());
+
+		auto entity_id = create_entity_for_branch(branch);
+
+		// Update the entity id of the asym
+		auto &struct_asym = m_db["struct_asym"];
+		auto r = struct_asym.find1("id"_key == asym_id);
+		r["entity_id"] = entity_id;
+
+		for (auto &sugar : branch)
+		{
+			for (auto atom : sugar.atoms())
+				atom.set_property("label_entity_id", entity_id);
+		}
+
+		auto &pdbx_branch_scheme = m_db["pdbx_branch_scheme"];
+		pdbx_branch_scheme.erase("asym_id"_key == asym_id);
+
+		for (auto &sugar : branch)
+		{
+			pdbx_branch_scheme.emplace({
+				{"asym_id", asym_id},
+				{"entity_id", entity_id},
+				{"num", sugar.num()},
+				{"mon_id", sugar.get_compound_id()},
+
+				{"pdb_asym_id", asym_id},
+				{"pdb_seq_num", sugar.num()},
+				{"pdb_mon_id", sugar.get_compound_id()},
+
+				// TODO: need fix, collect from nag_atoms?
+				{"auth_asym_id", asym_id},
+				{"auth_mon_id", sugar.get_compound_id()},
+				{"auth_seq_num", sugar.get_auth_seq_id()},
+
+				{"hetero", "n"}
+			});
+		}
+	}
+}
+
+void structure::remove_branch(branch &branch)
+{
+	using namespace literals;
+
+	m_db["pdbx_branch_scheme"].erase("asym_id"_key == branch.get_asym_id());
+	m_db["struct_asym"].erase("id"_key == branch.get_asym_id());
+
+	for (auto &sugar : branch)
+	{
+		auto atoms = sugar.atoms();
+		for (auto atom : atoms)
+			remove_atom(atom);
+	}
+
+	m_branches.erase(remove(m_branches.begin(), m_branches.end(), branch), m_branches.end());
+}
+
+std::string structure::create_non_poly_entity(const std::string &comp_id)
+{
+	return insert_compound(comp_id, true);
+}
+
+std::string structure::create_non_poly(const std::string &entity_id, const std::vector<atom> &atoms)
+{
+	using namespace cif::literals;
+
+	auto &struct_asym = m_db["struct_asym"];
+	std::string asym_id = struct_asym.get_unique_id();
+
+	struct_asym.emplace({
+		{"id", asym_id},
+		{"pdbx_blank_PDB_chainid_flag", "N"},
+		{"pdbx_modified", "N"},
+		{"entity_id", entity_id},
+		{"details", "?"}
+	});
+
+	std::string comp_id = m_db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id");
+
+	auto &atom_site = m_db["atom_site"];
+
+	auto &res = m_non_polymers.emplace_back(*this, comp_id, asym_id, 0, "1");
+
+	for (auto &atom : atoms)
+	{
+		auto atom_id = atom_site.get_unique_id("");
+
+		auto row = atom_site.emplace({
+			{"group_PDB", atom.get_property("group_PDB")},
+			{"id", atom_id},
+			{"type_symbol", atom.get_property("type_symbol")},
+			{"label_atom_id", atom.get_property("label_atom_id")},
+			{"label_alt_id", atom.get_property("label_alt_id")},
+			{"label_comp_id", comp_id},
+			{"label_asym_id", asym_id},
+			{"label_entity_id", entity_id},
+			{"label_seq_id", "."},
+			{"pdbx_PDB_ins_code", ""},
+			{"Cartn_x", atom.get_property("Cartn_x")},
+			{"Cartn_y", atom.get_property("Cartn_y")},
+			{"Cartn_z", atom.get_property("Cartn_z")},
+			{"occupancy", atom.get_property("occupancy")},
+			{"B_iso_or_equiv", atom.get_property("B_iso_or_equiv")},
+			{"pdbx_formal_charge", atom.get_property("pdbx_formal_charge")},
+			{"auth_seq_id", 1},
+			{"auth_comp_id", comp_id},
+			{"auth_asym_id", asym_id},
+			{"auth_atom_id", atom.get_property("label_atom_id")},
+			{"pdbx_PDB_model_num", 1}
+		});
+
+		auto &newAtom = emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id, row));
+		res.add_atom(newAtom);
+	}
+
+	auto &pdbx_nonpoly_scheme = m_db["pdbx_nonpoly_scheme"];
+	int ndb_nr = pdbx_nonpoly_scheme.find("asym_id"_key == asym_id and "entity_id"_key == entity_id).size() + 1;
+	pdbx_nonpoly_scheme.emplace({
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"mon_id", comp_id},
+		{"ndb_seq_num", ndb_nr},
+		{"pdb_seq_num", res.get_auth_seq_id()},
+		{"auth_seq_num", res.get_auth_seq_id()},
+		{"pdb_mon_id", comp_id},
+		{"auth_mon_id", comp_id},
+		{"pdb_strand_id", asym_id},
+		{"pdb_ins_code", "."},
+	});
+
+	return asym_id;
+}
+
+std::string structure::create_non_poly(const std::string &entity_id, std::vector<std::vector<item>> &atom_info)
+{
+	using namespace literals;
+
+	auto &struct_asym = m_db["struct_asym"];
+	std::string asym_id = struct_asym.get_unique_id();
+
+	struct_asym.emplace({
+		{"id", asym_id},
+		{"pdbx_blank_PDB_chainid_flag", "N"},
+		{"pdbx_modified", "N"},
+		{"entity_id", entity_id},
+		{"details", "?"}
+	});
+
+	std::string comp_id = m_db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id");
+
+	auto &atom_site = m_db["atom_site"];
+
+	auto &res = m_non_polymers.emplace_back(*this, comp_id, asym_id, 0, "1");
+
+	auto appendUnlessSet = [](std::vector<item> &ai, item &&i)
+	{
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](item &ci)
+				{ return ci.name() == name; }) == ai.end())
+			ai.emplace_back(std::move(i));
+	};
+
+	for (auto &atom : atom_info)
+	{
+		auto atom_id = atom_site.get_unique_id("");
+
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_comp_id", comp_id});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_seq_id", ""});
+		appendUnlessSet(atom, {"label_entity_id", entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", comp_id});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"auth_seq_id", 1});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
+
+		auto row = atom_site.emplace(atom.begin(), atom.end());
+
+		auto &newAtom = emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id, row));
+		res.add_atom(newAtom);
+	}
+
+	auto &pdbx_nonpoly_scheme = m_db["pdbx_nonpoly_scheme"];
+	int ndb_nr = pdbx_nonpoly_scheme.find("asym_id"_key == asym_id and "entity_id"_key == entity_id).size() + 1;
+	pdbx_nonpoly_scheme.emplace({
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"mon_id", comp_id},
+		{"ndb_seq_num", ndb_nr},
+		{"pdb_seq_num", res.get_auth_seq_id()},
+		{"auth_seq_num", res.get_auth_seq_id()},
+		{"pdb_mon_id", comp_id},
+		{"auth_mon_id", comp_id},
+		{"pdb_strand_id", asym_id},
+		{"pdb_ins_code", "."},
+	});
+
+	return asym_id;
+}
+
+branch &structure::create_branch(std::vector<std::vector<item>> &nag_atoms)
+{
+	// sanity check
+	for (auto &nag_atom : nag_atoms)
+	{
+		for (auto info : nag_atom)
+		{
+			if (info.name() == "label_comp_id" and info.value() != "NAG")
+				throw std::logic_error("The first sugar in a branch should be a NAG");
+		}
+	}
+
+	using namespace literals;
+
+	auto &struct_asym = m_db["struct_asym"];
+	std::string asym_id = struct_asym.get_unique_id();
+
+	auto &branch = m_branches.emplace_back(*this, asym_id);
+	auto &sugar = branch.emplace_back(branch, "NAG", asym_id, 1);
+	auto tmp_entity_id = m_db["entity"].get_unique_id("");
+
+	auto &atom_site = m_db["atom_site"];
+
+	auto appendUnlessSet = [](std::vector<item> &ai, item &&i)
+	{
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](item &ci)
+				{ return ci.name() == name; }) == ai.end())
+			ai.emplace_back(std::move(i));
+	};
+
+	for (auto &atom : nag_atoms)
+	{
+		auto atom_id = atom_site.get_unique_id("");
+
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_comp_id", "NAG"});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_seq_id", "."});
+		appendUnlessSet(atom, {"label_entity_id", tmp_entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", "NAG"});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"auth_seq_id", 1});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
+
+		auto row = atom_site.emplace(atom.begin(), atom.end());
+
+		auto &newAtom = emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id, row));
+		sugar.add_atom(newAtom);
+	}
+
+	// now we can create the entity and get the real ID
+	auto entity_id = create_entity_for_branch(branch);
+
+	struct_asym.emplace({
+		{"id", asym_id},
+		{"pdbx_blank_PDB_chainid_flag", "N"},
+		{"pdbx_modified", "N"},
+		{"entity_id", entity_id},
+		{"details", "?"}
+	});
+
+	for (auto &a : sugar.atoms())
+		a.set_property("label_entity_id", entity_id);
+
+	m_db["pdbx_branch_scheme"].emplace({
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"num", 1},
+		{"mon_id", "NAG"},
+
+		{"pdb_asym_id", asym_id},
+		{"pdb_seq_num", 1},
+		{"pdb_mon_id", "NAG"},
+
+		// TODO: need fix, collect from nag_atoms?
+		{"auth_asym_id", asym_id},
+		{"auth_mon_id", "NAG"},
+		{"auth_seq_num", 1},
+
+		{"hetero", "n"}
+	});
+
+	return branch;
+}
+
+branch &structure::extend_branch(const std::string &asym_id, std::vector<std::vector<item>> &atom_info,
+	int link_sugar, const std::string &link_atom)
+{
+	// sanity check
+	std::string compoundID;
+
+	for (auto &atom : atom_info)
+	{
+		for (auto info : atom)
+		{
+			if (info.name() != "label_comp_id")
+				continue;
+
+			if (compoundID.empty())
+				compoundID = info.value();
+			else if (info.value() != compoundID)
+				throw std::logic_error("All atoms should be of the same type");
+		}
+	}
+
+	using namespace literals;
+
+	// auto &branch = m_branches.emplace_back(*this, asym_id);
+	auto tmp_entity_id = m_db["entity"].get_unique_id("");
+
+	auto &atom_site = m_db["atom_site"];
+
+	auto appendUnlessSet = [](std::vector<item> &ai, item &&i)
+	{
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](item &ci)
+				{ return ci.name() == name; }) == ai.end())
+			ai.emplace_back(std::move(i));
+	};
+
+	auto bi = std::find_if(m_branches.begin(), m_branches.end(), [asym_id](branch &b)
+		{ return b.get_asym_id() == asym_id; });
+	if (bi == m_branches.end())
+		throw std::logic_error("Create a branch first!");
+
+	branch &branch = *bi;
+
+	int sugarNum = branch.size() + 1;
+
+	auto &sugar = branch.emplace_back(branch, compoundID, asym_id, sugarNum);
+
+	for (auto &atom : atom_info)
+	{
+		auto atom_id = atom_site.get_unique_id("");
+
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_comp_id", compoundID});
+		appendUnlessSet(atom, {"label_entity_id", tmp_entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", compoundID});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
+
+		auto row = atom_site.emplace(atom.begin(), atom.end());
+
+		auto &newAtom = emplace_atom(std::make_shared<atom::atom_impl>(m_db, atom_id, row));
+		sugar.add_atom(newAtom);
+	}
+
+	sugar.set_link(branch.at(link_sugar - 1).get_atom_by_atom_id(link_atom));
+
+	auto entity_id = create_entity_for_branch(branch);
+
+	// Update the entity id of the asym
+	auto &struct_asym = m_db["struct_asym"];
+	auto r = struct_asym.find1("id"_key == asym_id);
+	r["entity_id"] = entity_id;
+
+	for (auto &sugar : branch)
+	{
+		for (auto atom : sugar.atoms())
+			atom.set_property("label_entity_id", entity_id);
+	}
+
+	auto &pdbx_branch_scheme = m_db["pdbx_branch_scheme"];
+	pdbx_branch_scheme.erase("asym_id"_key == asym_id);
+
+	for (auto &sugar : branch)
+	{
+		pdbx_branch_scheme.emplace({
+			{"asym_id", asym_id},
+			{"entity_id", entity_id},
+			{"num", sugar.num()},
+			{"mon_id", sugar.get_compound_id()},
+
+			{"pdb_asym_id", asym_id},
+			{"pdb_seq_num", sugar.num()},
+			{"pdb_mon_id", sugar.get_compound_id()},
+
+			// TODO: need fix, collect from nag_atoms?
+			{"auth_asym_id", asym_id},
+			{"auth_mon_id", sugar.get_compound_id()},
+			{"auth_seq_num", sugar.get_auth_seq_id()},
+
+			{"hetero", "n"}
+		});
+	}
+
+	return branch;
+}
+
+std::string structure::create_entity_for_branch(branch &branch)
+{
+	using namespace literals;
+
+	std::string entityName = branch.name(), entityID;
+
+	auto &entity = m_db["entity"];
+
+	try
+	{
+		entityID = entity.find1<std::string>("type"_key == "branched" and "pdbx_description"_key == entityName, "id");
+	}
+	catch (const std::exception &e)
+	{
+		entityID = entity.get_unique_id("");
+
+		if (VERBOSE)
+			std::cout << "Creating new entity " << entityID << " for branched sugar " << entityName << std::endl;
+
+		entity.emplace({{"id", entityID},
+			{"type", "branched"},
+			{"src_method", "man"},
+			{"pdbx_description", entityName},
+			{"formula_weight", branch.weight()}});
+		auto &pdbx_entity_branch_list = m_db["pdbx_entity_branch_list"];
+		for (auto &sugar : branch)
+		{
+			pdbx_entity_branch_list.emplace({
+				{"entity_id", entityID},
+				{"comp_id", sugar.get_compound_id()},
+				{"num", sugar.num()},
+				{"hetero", "n"}
+			});
+		}
+
+		auto &pdbx_entity_branch_link = m_db["pdbx_entity_branch_link"];
+		for (auto &s1 : branch)
+		{
+			auto l2 = s1.get_link();
+
+			if (not l2)
+				continue;
+
+			auto &s2 = branch.at(std::stoi(l2.get_auth_seq_id()) - 1);
+			auto l1 = s2.get_atom_by_atom_id("C1");
+
+			pdbx_entity_branch_link.emplace({
+				{"link_id", pdbx_entity_branch_link.get_unique_id("")},
+				{"entity_id", entityID},
+				{"entity_branch_list_num_1", s1.get_auth_seq_id()},
+				{"comp_id_1", s1.get_compound_id()},
+				{"atom_id_1", l1.get_label_atom_id()},
+				{"leaving_atom_id_1", "O1"},
+				{"entity_branch_list_num_2", s2.get_auth_seq_id()},
+				{"comp_id_2", s2.get_compound_id()},
+				{"atom_id_2", l2.get_label_atom_id()},
+				{"leaving_atom_id_2", "H" + l2.get_label_atom_id()},
+				{"value_order", "sing"}
+			});
+		}
+	}
+
+	return entityID;
+}
+
+void structure::cleanup_empty_categories()
+{
+	using namespace literals;
+
+	auto &atomSite = m_db["atom_site"];
+
+	// Remove chem_comp's for which there are no atoms at all
+	auto &chem_comp = m_db["chem_comp"];
+	std::vector<row_handle> obsoleteChemComps;
+
+	for (auto chemComp : chem_comp)
+	{
+		std::string compID = chemComp["id"].as<std::string>();
+		if (atomSite.exists("label_comp_id"_key == compID or "auth_comp_id"_key == compID))
+			continue;
+
+		obsoleteChemComps.push_back(chemComp);
+	}
+
+	for (auto chemComp : obsoleteChemComps)
+		chem_comp.erase(chemComp);
+
+	// similarly, remove entities not referenced by any atom
+
+	auto &entities = m_db["entity"];
+	std::vector<row_handle> obsoleteEntities;
+
+	for (auto entity : entities)
+	{
+		std::string entityID = entity["id"].as<std::string>();
+		if (atomSite.exists("label_entity_id"_key == entityID))
+			continue;
+
+		obsoleteEntities.push_back(entity);
+	}
+
+	for (auto entity : obsoleteEntities)
+		entities.erase(entity);
+
+	// the rest?
+
+	for (const char *cat : {"pdbx_entity_nonpoly"})
+	{
+		auto &category = m_db[cat];
+
+		std::vector<row_handle> empty;
+		for (auto row : category)
+		{
+			if (not category.has_children(row) and not category.has_parents(row))
+				empty.push_back(row);
+		}
+
+		for (auto row : empty)
+			category.erase(row);
+	}
+
+	// count molecules
+	for (auto entity : entities)
+	{
+		std::string type, id;
+		cif::tie(type, id) = entity.get("type", "id");
+
+		std::optional<size_t> count;
+		if (type == "polymer")
+			count = m_db["entity_poly"].find("entity_id"_key == id).size();
+		else if (type == "non-polymer" or type == "water")
+			count = m_db["pdbx_nonpoly_scheme"].find("entity_id"_key == id).size();
+		else if (type == "branched")
+		{
+			// is this correct?
+			std::set<std::string> asym_ids;
+			for (const auto &[asym_id] : m_db["pdbx_branch_scheme"].find<std::string>("entity_id"_key == id, "asym_id"))
+				asym_ids.insert(asym_id);
+			count = asym_ids.size();
+		}
+
+		entity["pdbx_number_of_molecules"] = count.value_or(0);
+	}
+}
+
+void structure::translate(point t)
+{
+	for (auto &a : m_atoms)
+		a.translate(t);
+}
+
+void structure::rotate(quaternion q)
+{
+	for (auto &a : m_atoms)
+		a.rotate(q);
+}
+
+void structure::translate_and_rotate(point t, quaternion q)
+{
+	for (auto &a : m_atoms)
+		a.translate_and_rotate(t, q);
+}
+
+void structure::translate_rotate_and_translate(point t1, quaternion q, point t2)
+{
+	for (auto &a : m_atoms)
+		a.translate_rotate_and_translate(t1, q, t2);
+}
+
+void structure::validate_atoms() const
+{
+	// validate order
+	assert(m_atoms.size() == m_atom_index.size());
+	for (size_t i = 0; i + i < m_atoms.size(); ++i)
+		assert(m_atoms[m_atom_index[i]].id().compare(m_atoms[m_atom_index[i + 1]].id()) < 0);
+
+	std::vector<atom> atoms = m_atoms;
+
+	auto removeAtomFromList = [&atoms](const atom &a)
+	{
+		auto i = std::find(atoms.begin(), atoms.end(), a);
+		assert(i != atoms.end());
+		atoms.erase(i);
+	};
+
+	for (auto &poly : m_polymers)
+	{
+		for (auto &monomer : poly)
+		{
+			for (auto &atom : monomer.atoms())
+				removeAtomFromList(atom);
+		}
+	}
+
+	for (auto &branch : m_branches)
+	{
+		for (auto &sugar : branch)
+		{
+			for (auto &atom : sugar.atoms())
+				removeAtomFromList(atom);
+		}
+	}
+
+	for (auto &res : m_non_polymers)
+	{
+		for (auto &atom : res.atoms())
+			removeAtomFromList(atom);
+	}
+
+	assert(atoms.empty());
+}
+
+} // namespace pdbx

src/row.cpp

@@ -53,6 +53,11 @@ uint16_t row_handle::add_column(std::string_view name)
 	return m_category->add_column(name);
 }
 
+void row_handle::swap(size_t column, row_handle &b)
+{
+	m_category->swap_item(column, *this, b);
+}
+
 // --------------------------------------------------------------------
 
 row_initializer::row_initializer(row_handle rh)

src/utilities.cpp

@@ -338,44 +338,44 @@ void ProgressImpl::PrintDone()
 }
 
 Progress::Progress(int64_t inMax, const std::string &inAction)
-	: mImpl(nullptr)
+	: m_impl(nullptr)
 {
 	if (isatty(STDOUT_FILENO))
-		mImpl = new ProgressImpl(inMax, inAction);
+		m_impl = new ProgressImpl(inMax, inAction);
 }
 
 Progress::~Progress()
 {
-	if (mImpl != nullptr)
-		mImpl->Stop();
+	if (m_impl != nullptr)
+		m_impl->Stop();
 
-	delete mImpl;
+	delete m_impl;
 }
 
 void Progress::consumed(int64_t inConsumed)
 {
-	if (mImpl != nullptr and
-		(mImpl->mConsumed += inConsumed) >= mImpl->mMax)
+	if (m_impl != nullptr and
+		(m_impl->mConsumed += inConsumed) >= m_impl->mMax)
 	{
-		mImpl->Stop();
+		m_impl->Stop();
 	}
 }
 
 void Progress::progress(int64_t inProgress)
 {
-	if (mImpl != nullptr and
-		(mImpl->mConsumed = inProgress) >= mImpl->mMax)
+	if (m_impl != nullptr and
+		(m_impl->mConsumed = inProgress) >= m_impl->mMax)
 	{
-		mImpl->Stop();
+		m_impl->Stop();
 	}
 }
 
 void Progress::message(const std::string &inMessage)
 {
-	if (mImpl != nullptr)
+	if (m_impl != nullptr)
 	{
-		std::unique_lock lock(mImpl->mMutex);
-		mImpl->mMessage = inMessage;
+		std::unique_lock lock(m_impl->mMutex);
+		m_impl->mMessage = inMessage;
 	}
 }
 

test/format-test.cpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * 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.
+ */
+
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+#include <boost/test/included/unit_test.hpp>
+
+#include <stdexcept>
+
+#include <cif++.hpp>
+
+namespace tt = boost::test_tools;
+
+std::filesystem::path gTestDir = std::filesystem::current_path(); // filled in first test
+
+// --------------------------------------------------------------------
+
+cif::file operator""_cf(const char *text, size_t length)
+{
+	struct membuf : public std::streambuf
+	{
+		membuf(char *text, size_t length)
+		{
+			this->setg(text, text, text + length);
+		}
+	} buffer(const_cast<char *>(text), length);
+
+	std::istream is(&buffer);
+	return cif::file(is);
+}
+
+// --------------------------------------------------------------------
+
+bool init_unit_test()
+{
+	cif::VERBOSE = 1;
+
+	// // not a test, just initialize test dir
+	// if (boost::unit_test::framework::master_test_suite().argc == 2)
+	// 	gTestDir = boost::unit_test::framework::master_test_suite().argv[1];
+
+	// // do this now, avoids the need for installing
+	// cif::add_file_resource("mmcif_pdbx.dic", gTestDir / ".." / "rsrc" / "mmcif_pdbx.dic");
+
+	// // initialize CCD location
+	// cif::add_file_resource("components.cif", gTestDir / ".." / "data" / "ccd-subset.cif");
+
+	return true;
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(fmt_1)
+{
+	std::ostringstream os;
+
+	std::string world("world");
+	os << cif::format("Hello, %-10.10s, the magic number is %d and pi is %g", world, 42, M_PI);
+	BOOST_CHECK_EQUAL(os.str(), "Hello, world     , the magic number is 42 and pi is 3.14159");
+
+	BOOST_CHECK_EQUAL(cif::format("Hello, %-10.10s, the magic number is %d and pi is %g", world, 42, M_PI).str(),
+		"Hello, world     , the magic number is 42 and pi is 3.14159");
+}
\ No newline at end of file

test/model-test.cpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ * 
+ * Copyright (c) 2021 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.
+ */
+
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+#include <boost/test/included/unit_test.hpp>
+
+#include <stdexcept>
+
+#include <cif++.hpp>
+
+// --------------------------------------------------------------------
+
+cif::file operator""_cf(const char* text, size_t length)
+{
+    struct membuf : public std::streambuf
+    {
+        membuf(char* text, size_t length)
+        {
+            this->setg(text, text, text + length);
+        }
+    } buffer(const_cast<char*>(text), length);
+
+    std::istream is(&buffer);
+    return cif::file(is);
+}
+
+// --------------------------------------------------------------------
+
+std::filesystem::path gTestDir = std::filesystem::current_path();
+
+bool init_unit_test()
+{
+    cif::VERBOSE = 1;
+
+	// not a test, just initialize test dir
+	if (boost::unit_test::framework::master_test_suite().argc == 2)
+		gTestDir = boost::unit_test::framework::master_test_suite().argv[1];
+
+	// do this now, avoids the need for installing
+	cif::add_file_resource("mmcif_pdbx.dic", gTestDir / ".." / "rsrc" / "mmcif_pdbx.dic");
+
+	// initialize CCD location
+	cif::add_file_resource("components.cif", gTestDir / ".." / "data" / "ccd-subset.cif");
+
+	cif::compound_factory::instance().push_dictionary(gTestDir / "HEM.cif");
+
+	return true;
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_nonpoly_1)
+{
+    cif::VERBOSE = 1;
+
+	cif::file file;
+	file.load_dictionary("mmcif_pdbx.dic");
+	file.emplace("TEST");	// create a datablock
+	
+	cif::mm::structure structure(file);
+
+	std::string entity_id = structure.create_non_poly_entity("HEM");
+
+	auto atoms = R"(
+data_HEM
+loop_
+_atom_site.group_PDB
+_atom_site.type_symbol
+_atom_site.label_atom_id
+_atom_site.label_alt_id
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+HETATM C  CHA . ? -5.248  39.769 -0.250  1.00 7.67  ?
+HETATM C  CHB . ? -3.774  36.790 3.280   1.00 7.05  ?
+HETATM C  CHC . ? -2.879  33.328 0.013   1.00 7.69  ?
+HETATM C  CHD . ? -4.342  36.262 -3.536  1.00 8.00  ?
+# that's enough to test with
+)"_cf;
+
+	auto &hem_data = atoms["HEM"];
+	auto &atom_site = hem_data["atom_site"];
+
+	auto hem_atoms = atom_site.rows();
+	std::vector<cif::mm::atom> atom_data;
+	for (auto hem_atom: hem_atoms)
+		atom_data.emplace_back(hem_data, hem_atom);
+
+	structure.create_non_poly(entity_id, atom_data);
+
+	auto expected = R"(
+data_TEST
+# 
+_pdbx_nonpoly_scheme.asym_id         A 
+_pdbx_nonpoly_scheme.ndb_seq_num     1 
+_pdbx_nonpoly_scheme.entity_id       1 
+_pdbx_nonpoly_scheme.mon_id          HEM 
+_pdbx_nonpoly_scheme.pdb_seq_num     1 
+_pdbx_nonpoly_scheme.auth_seq_num    1 
+_pdbx_nonpoly_scheme.pdb_mon_id      HEM 
+_pdbx_nonpoly_scheme.auth_mon_id     HEM 
+_pdbx_nonpoly_scheme.pdb_strand_id   A 
+_pdbx_nonpoly_scheme.pdb_ins_code    . 
+#
+loop_
+_atom_site.id
+_atom_site.auth_asym_id
+_atom_site.label_alt_id
+_atom_site.label_asym_id
+_atom_site.label_atom_id
+_atom_site.label_comp_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.type_symbol
+_atom_site.group_PDB
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? 1 HEM CHA 1
+2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? 1 HEM CHB 1
+3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? 1 HEM CHD 1
+#
+_chem_comp.id               HEM
+_chem_comp.type             NON-POLYMER
+_chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
+_chem_comp.formula          'C34 H32 Fe N4 O4'
+_chem_comp.formula_weight   616.487000
+#
+_pdbx_entity_nonpoly.entity_id   1
+_pdbx_entity_nonpoly.name        'PROTOPORPHYRIN IX CONTAINING FE'
+_pdbx_entity_nonpoly.comp_id     HEM
+#
+_entity.id                 1
+_entity.type               non-polymer
+_entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
+_entity.formula_weight     616.487000
+#
+_struct_asym.id                            A
+_struct_asym.entity_id                     1
+_struct_asym.pdbx_blank_PDB_chainid_flag   N
+_struct_asym.pdbx_modified                 N
+_struct_asym.details                       ?
+#
+)"_cf;
+
+	expected.load_dictionary("mmcif_pdbx.dic");
+
+	if (not (expected.front() == structure.get_datablock()))
+	{
+		BOOST_TEST(false);
+		std::cout << expected.front() << std::endl
+				<< std::endl
+				<< structure.get_datablock() << std::endl;
+	}
+}
+
+// // --------------------------------------------------------------------
+
+// BOOST_AUTO_TEST_CASE(test_load_1)
+// {
+// 	cif::file cf(gTestDir / "5v3g.cif.gz");
+// 	cif::mm::structure s(cf);
+
+// 	for (auto &poly : s.polymers())
+// 	{
+// 		std::cout << std::string(80, '=') << std::endl;
+// 		for (auto &res : poly)
+// 		{
+// 			std::cout << res << std::endl;
+
+// 			for (auto &atom : res.atoms())
+// 				std::cout << "  " << atom << std::endl;
+// 		}
+// 	}
+// }
+
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(test_atom_id)
+{
+	auto data = R"(
+data_TEST
+# 
+_pdbx_nonpoly_scheme.asym_id         A 
+_pdbx_nonpoly_scheme.ndb_seq_num     1 
+_pdbx_nonpoly_scheme.entity_id       1 
+_pdbx_nonpoly_scheme.mon_id          HEM 
+_pdbx_nonpoly_scheme.pdb_seq_num     1 
+_pdbx_nonpoly_scheme.auth_seq_num    1 
+_pdbx_nonpoly_scheme.pdb_mon_id      HEM 
+_pdbx_nonpoly_scheme.auth_mon_id     HEM 
+_pdbx_nonpoly_scheme.pdb_strand_id   A 
+_pdbx_nonpoly_scheme.pdb_ins_code    . 
+#
+loop_
+_atom_site.id
+_atom_site.auth_asym_id
+_atom_site.label_alt_id
+_atom_site.label_asym_id
+_atom_site.label_atom_id
+_atom_site.label_comp_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.type_symbol
+_atom_site.group_PDB
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? 1 HEM CHA 1
+3 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? 1 HEM CHB 1
+2 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? 1 HEM CHD 1
+#
+_chem_comp.id               HEM
+_chem_comp.type             NON-POLYMER
+_chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
+_chem_comp.formula          'C34 H32 Fe N4 O4'
+_chem_comp.formula_weight   616.487000
+#
+_pdbx_entity_nonpoly.entity_id   1
+_pdbx_entity_nonpoly.name        'PROTOPORPHYRIN IX CONTAINING FE'
+_pdbx_entity_nonpoly.comp_id     HEM
+#
+_entity.id                 1
+_entity.type               non-polymer
+_entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
+_entity.formula_weight     616.487000
+#
+_struct_asym.id                            A
+_struct_asym.entity_id                     1
+_struct_asym.pdbx_blank_PDB_chainid_flag   N
+_struct_asym.pdbx_modified                 N
+_struct_asym.details                       ?
+#
+)"_cf;
+
+	data.load_dictionary("mmcif_pdbx.dic");
+
+	cif::mm::structure s(data);
+
+	BOOST_CHECK_EQUAL(s.get_atom_by_id("1").get_label_atom_id(), "CHA");
+	BOOST_CHECK_EQUAL(s.get_atom_by_id("2").get_label_atom_id(), "CHC");
+	BOOST_CHECK_EQUAL(s.get_atom_by_id("3").get_label_atom_id(), "CHB");
+	BOOST_CHECK_EQUAL(s.get_atom_by_id("4").get_label_atom_id(), "CHD");
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(atom_numbers_1)
+{
+	const std::filesystem::path test1(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	cif::file file(test1.string());
+	cif::mm::structure structure(file);
+
+	auto &db = file.front();
+
+	auto &atoms = structure.atoms();
+	auto ai = atoms.begin();
+
+	for (const auto &[id, label_asym_id, label_seq_id, label_atom_id, auth_seq_id, label_comp_id] :
+		db["atom_site"].rows<std::string,std::string,int,std::string,std::string,std::string>("id", "label_asym_id", "label_seq_id", "label_atom_id", "auth_seq_id", "label_comp_id"))
+	{
+		auto atom = structure.get_atom_by_id(id);
+
+		BOOST_CHECK_EQUAL(atom.get_label_asym_id(), label_asym_id);
+		BOOST_CHECK_EQUAL(atom.get_label_seq_id(), label_seq_id);
+		BOOST_CHECK_EQUAL(atom.get_label_atom_id(), label_atom_id);
+		BOOST_CHECK_EQUAL(atom.get_auth_seq_id(), auth_seq_id);
+		BOOST_CHECK_EQUAL(atom.get_label_comp_id(), label_comp_id);
+
+		BOOST_ASSERT(ai != atoms.end());
+
+		BOOST_CHECK_EQUAL(ai->id(), id);
+		++ai;
+	}
+
+	BOOST_ASSERT(ai == atoms.end());
+}
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(test_load_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	cif::file file(example.string());
+
+	auto &db = file.front();
+
+	cif::mm::structure s(file);
+
+	BOOST_CHECK(s.polymers().size() == 1);
+
+	auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
+
+	for (auto &poly : s.polymers())
+	{
+		BOOST_CHECK_EQUAL(poly.size(), pdbx_poly_seq_scheme.find("asym_id"_key == poly.get_asym_id()).size());
+	}
+}
+
+BOOST_AUTO_TEST_CASE(remove_residue_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	cif::file file(example.string());
+
+	cif::mm::structure s(file);
+	s.remove_residue(s.get_residue("B"));
+
+	BOOST_CHECK_NO_THROW(s.validate_atoms());
+}

test/rename-compound-test.cpp

+#include "../include/cif++/Cif++.hpp"
+#include "../include/cif++/PDB2Cif.hpp"
+#include "../include/cif++/Structure.hpp"
+
+#include <iostream>
+#include <fstream>
+
+int main(int argc, char* argv[])
+{
+	cif::VERBOSE = 3;
+
+	try
+	{
+		std::filesystem::path testdir = std::filesystem::current_path();
+
+		if (argc == 3)
+			testdir = argv[2];
+
+		if (std::filesystem::exists(testdir / ".." / "data" / "ccd-subset.cif"))
+			cif::add_file_resource("components.cif", testdir / ".." / "data" / "ccd-subset.cif");
+
+		if (std::filesystem::exists(testdir / ".." / "rsrc" / "mmcif_pdbx.dic"))
+			cif::add_file_resource("mmcif_pdbx.dic", testdir / ".." / "rsrc" / "mmcif_pdbx.dic");
+
+		mmcif::CompoundFactory::instance().pushDictionary(testdir / "REA.cif");
+		mmcif::CompoundFactory::instance().pushDictionary(testdir / "RXA.cif");
+
+		mmcif::file f(testdir / ".."/"examples"/"1cbs.cif.gz");
+		mmcif::Structure structure(f);
+
+		auto &res = structure.getResidue("B");
+		structure.change_residue(res, "RXA", {});
+
+		structure.cleanupEmptyCategories();
+
+		f.save(std::cout);
+	}
+	catch (const std::exception& e)
+	{
+		std::cerr << e.what() << std::endl;
+		exit(1);
+	}
+	
+	return 0;	
+}

test/sugar-test.cpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ * 
+ * Copyright (c) 2021 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.
+ */
+
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+#include <boost/test/included/unit_test.hpp>
+
+#include <stdexcept>
+
+#include <cif++/cif.hpp>
+#include <cif++/structure/Structure.hpp>
+
+// --------------------------------------------------------------------
+
+cif::file operator""_cf(const char* text, size_t length)
+{
+    struct membuf : public std::streambuf
+    {
+        membuf(char* text, size_t length)
+        {
+            this->setg(text, text, text + length);
+        }
+    } buffer(const_cast<char*>(text), length);
+
+    std::istream is(&buffer);
+    return cif::file(is);
+}
+
+// --------------------------------------------------------------------
+
+std::filesystem::path gTestDir = std::filesystem::current_path();
+
+bool init_unit_test()
+{
+    cif::VERBOSE = 1;
+
+	// not a test, just initialize test dir
+	if (boost::unit_test::framework::master_test_suite().argc == 2)
+		gTestDir = boost::unit_test::framework::master_test_suite().argv[1];
+
+	// do this now, avoids the need for installing
+	cif::add_file_resource("mmcif_pdbx.dic", gTestDir / ".." / "rsrc" / "mmcif_pdbx.dic");
+
+	// initialize CCD location
+	cif::add_file_resource("components.cif", gTestDir / ".." / "data" / "ccd-subset.cif");
+
+	mmcif::CompoundFactory::instance().pushDictionary(gTestDir / "HEM.cif");
+
+	return true;
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(sugar_name_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::file file(example.string());
+	mmcif::Structure s(file);
+
+	auto &db = s.datablock();
+	auto &entity = db["entity"];
+
+	auto &branches = s.branches();
+
+	BOOST_CHECK_EQUAL(branches.size(), 4);
+
+	for (auto &branch : branches)
+	{
+		auto entityID = branch.front().entityID();
+
+		auto name = entity.find1<std::string>("id"_key == entityID, "pdbx_description");
+		BOOST_CHECK_EQUAL(branch.name(), name);
+	}
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_sugar_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::file file(example.string());
+	mmcif::Structure s(file);
+
+	// collect atoms from asym L first
+	auto &NAG = s.getResidue("L");
+	auto nagAtoms = NAG.atoms();
+
+	std::vector<std::vector<cif::Item>> ai;
+
+	auto &db = s.datablock();
+	auto &as = db["atom_site"];
+
+	for (auto r : as.find("label_asym_id"_key == "L"))
+		ai.emplace_back(r.begin(), r.end());
+
+	s.remove_residue(NAG);
+
+	auto &branch = s.create_branch(ai);
+
+	BOOST_CHECK_EQUAL(branch.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose");
+	BOOST_CHECK_EQUAL(branch.size(), 1);
+
+	BOOST_CHECK_EQUAL(branch[0].atoms().size(), nagAtoms.size());
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_sugar_2)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::file file(example.string());
+	mmcif::Structure s(file);
+
+	// Get branch for H
+	auto &bH = s.getBranchByAsymID("H");
+
+	BOOST_CHECK_EQUAL(bH.size(), 2);
+
+	std::vector<std::vector<cif::Item>> ai[2];
+
+	auto &db = s.datablock();
+	auto &as = db["atom_site"];
+
+	for (size_t i = 0; i < 2; ++i)
+	{
+		for (auto r : as.find("label_asym_id"_key == "H" and "auth_seq_id"_key == i + 1))
+			ai[i].emplace_back(r.begin(), r.end());
+	}
+
+	s.remove_branch(bH);
+
+	auto &bN = s.create_branch(ai[0]);
+	s.extend_branch(bN.asymID(), ai[1], 1, "O4");
+
+	BOOST_CHECK_EQUAL(bN.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose-(1-4)-2-acetamido-2-deoxy-beta-D-glucopyranose");
+	BOOST_CHECK_EQUAL(bN.size(), 2);
+
+	file.save(gTestDir / "test-create_sugar_2.cif");
+
+	BOOST_CHECK_NO_THROW(mmcif::Structure s2(file));
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(delete_sugar_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::file file(example.string());
+	mmcif::Structure s(file);
+
+	// Get branch for H
+	auto &bG = s.getBranchByAsymID("G");
+
+	BOOST_CHECK_EQUAL(bG.size(), 4);
+
+	s.remove_residue(bG[1]);
+
+	BOOST_CHECK_EQUAL(bG.size(), 1);
+
+	auto &bN = s.getBranchByAsymID("G");
+
+	BOOST_CHECK_EQUAL(bN.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose");
+	BOOST_CHECK_EQUAL(bN.size(), 1);
+
+	file.save(gTestDir / "test-create_sugar_3.cif");
+
+	BOOST_CHECK_NO_THROW(mmcif::Structure s2(file));
+}

test/unit-v2-test.cpp

@@ -524,6 +524,37 @@ _test.value
 
 // --------------------------------------------------------------------
 
+BOOST_AUTO_TEST_CASE(sw_1)
+{
+	using namespace cif::literals;
+
+	auto f = R"(data_TEST
+#
+loop_
+_test.id
+_test.name
+_test.value
+1 aap   1.0
+2 noot  1.1
+3 mies  1.2
+    )"_cf;
+
+	auto &db = f.front();
+	auto &test = db["test"];
+
+	swap(test.front()["name"], test.back()["name"]);
+
+	BOOST_CHECK_EQUAL(test.find1<std::string>("id"_key == 1, "name"), "mies");
+	BOOST_CHECK_EQUAL(test.find1<std::string>("id"_key == 3, "name"), "aap");
+
+	swap(test.front()["name"], test.back()["name"]);
+
+	BOOST_CHECK_EQUAL(test.find1<std::string>("id"_key == 1, "name"), "aap");
+	BOOST_CHECK_EQUAL(test.find1<std::string>("id"_key == 3, "name"), "mies");
+}
+
+// --------------------------------------------------------------------
+
 BOOST_AUTO_TEST_CASE(d1)
 {
 	const char dict[] = R"(
@@ -2017,7 +2048,7 @@ _cat_3.num
 
 // 	for (const auto &[compound, seqnr] : std::initializer_list<std::tuple<std::string, int>>{{"PRO", 1}, {"ASN", 2}, {"PHE", 3}})
 // 	{
-// 		auto &res = structure.getResidue("A", compound, seqnr, "");
+// 		auto &res = structure.get_residue("A", compound, seqnr, "");
 // 		auto atoms = res.atoms();
 
 // 		auto dc = components.get(compound);