Merge branch 'develop' into trunk

Added another rotate/translate method to mmcif::Structure
Fix Structure::createNonPoly to add atoms...
2026-06-04 22:14:24 +08:00 · 2022-01-12 16:41:18 +01:00 · 2022-01-12 14:06:32 +01:00 · 2022-01-11 11:21:56 +01:00 · 2022-01-05 15:54:23 +01:00 · 2022-01-05 15:45:27 +01:00
34 changed files with 3763 additions and 4138 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,7 +25,7 @@
 cmake_minimum_required(VERSION 3.16)

 # set the project name
-project(cifpp VERSION 2.0.3 LANGUAGES CXX)
+project(cifpp VERSION 3.0.2 LANGUAGES CXX)

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

@@ -134,8 +134,6 @@ if(MSVC)
 		message(FATAL_ERROR "Unsupported or unknown processor type ${CMAKE_SYSTEM_PROCESSOR}")
 	endif()	

-	set(COFF_SPEC "--coff=${COFF_TYPE}")
-
 	# for mrc, just in case
 	list(APPEND CMAKE_PREFIX_PATH "$ENV{LOCALAPPDATA}/mrc")
 endif()
@@ -151,12 +149,12 @@ endif()
 if(WIN32 AND BUILD_SHARED_LIBS)
 	message("Not using resources when building shared libraries for Windows")
 else()
-	find_program(MRC mrc)
+	find_package(Mrc)

-	if(MRC)
-		option(CIFPP_USE_RSRC "Use mrc to create resources" ON)
+	if(MRC_FOUND)
+		option(USE_RSRC "Use mrc to create resources" ON)
 	else()
-		message("Using resources not possible since mrc was not found")
+		message(WARNING "Not using resources since mrc was not found")
 	endif()

 	if(CIFPP_USE_RSRC STREQUAL "ON")
@@ -173,29 +171,22 @@ set(CMAKE_THREAD_PREFER_PTHREAD)
 set(THREADS_PREFER_PTHREAD_FLAG)
 find_package(Threads)

-set(Boost_DETAILED_FAILURE_MSG ON)
-if(NOT BUILD_SHARED_LIBS)
-	set(Boost_USE_STATIC_LIBS ON)
-endif()
 find_package(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_options)

 if(NOT MSVC AND Boost_USE_STATIC_LIBS)
 	find_package(ZLIB REQUIRED)
-	find_package(BZip2 REQUIRED)
 endif()

 # Create a revision file, containing the current git version info

-find_package(Git)
-if(GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git")
-	include(GetGitRevisionDescription)
+include(GetGitRevisionDescription)
+option(GENERATE_CUSTOM_VERSION "Generate a custom version string" OFF)
+if(GIT-NOTFOUND OR HEAD-HASH-NOTFOUND OR NOT GENERATE_CUSTOM_VERSION)
    get_git_head_revision(REFSPEC COMMITHASH)

 	# Generate our own version string
 	git_describe_working_tree(BUILD_VERSION_STRING --match=build --dirty)
 else()
-    message(WARNING "Git not found, cannot set version info")
-
    SET(BUILD_VERSION_STRING ${PROJECT_VERSION})
 endif()

@@ -209,7 +200,7 @@ if(CIFPP_RECREATE_SYMOP_DATA)

 	add_executable(symop-map-generator "${CMAKE_SOURCE_DIR}/tools/symop-map-generator.cpp")

-	target_link_libraries(symop-map-generator Threads::Threads ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES} ${BZip2_LIBRARIES})
+	target_link_libraries(symop-map-generator Threads::Threads ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES})
 	if(Boost_INCLUDE_DIR)
 		target_include_directories(symop-map-generator PUBLIC ${Boost_INCLUDE_DIR})
 	endif()
@@ -256,7 +247,6 @@ set(project_headers
 	${PROJECT_SOURCE_DIR}/include/cif++/CifUtils.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/CifValidator.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/Compound.hpp
-	${PROJECT_SOURCE_DIR}/include/cif++/Matrix.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/PDB2Cif.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/PDB2CifRemark3.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/Point.hpp
@@ -281,7 +271,7 @@ target_include_directories(cifpp
 	${CMAKE_BINARY_DIR}
 )

-target_link_libraries(cifpp Threads::Threads ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES} ${BZip2_LIBRARIES})
+target_link_libraries(cifpp Threads::Threads ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES})

 if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
    target_link_options(cifpp PRIVATE -undefined dynamic_lookup)
@@ -341,6 +331,13 @@ install(TARGETS cifpp
 	RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
 	INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

+if(MSVC)
+	install(
+		FILES $<TARGET_PDB_FILE:${PROJECT_NAME}>
+		DESTINATION ${CMAKE_INSTALL_LIBDIR}
+		OPTIONAL)
+endif()
+
 install(EXPORT cifppTargets
 	FILE "cifppTargets.cmake"
 	NAMESPACE cifpp::
@@ -414,13 +411,6 @@ option(CIFPP_BUILD_TESTS "Build test exectuables" OFF)

 if(CIFPP_BUILD_TESTS)

-	if(CIFPP_USE_RSRC)
-		add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj
-			COMMAND ${MRC} -o ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj ${CMAKE_SOURCE_DIR}/rsrc/mmcif_pdbx_v50.dic ${COFF_SPEC}
-		)
-		set(CIFPP_TEST_RESOURCE ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj)
-	endif()
-
 	list(APPEND CIFPP_tests
 		# pdb2cif
 		rename-compound
@@ -431,15 +421,19 @@ if(CIFPP_BUILD_TESTS)
 		set(CIFPP_TEST "${CIFPP_TEST}-test")
 		set(CIFPP_TEST_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/test/${CIFPP_TEST}.cpp")

-		add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE} ${CIFPP_TEST_RESOURCE})
+		add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE})

 		target_include_directories(${CIFPP_TEST} PRIVATE
 			${CMAKE_CURRENT_SOURCE_DIR}/include
 			${CMAKE_CURRENT_BINARY_DIR}  # for config.h
 		)

-		target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES} ${BZip2_LIBRARIES})
-	
+		target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp ${Boost_LIBRARIES} std::filesystem ${ZLIB_LIBRARIES})
+
+		if(CIFPP_USE_RSRC)
+			mrc_target_resources(${CIFPP_TEST} ${CMAKE_SOURCE_DIR}/rsrc/mmcif_pdbx_v50.dic)
+		endif()
+
 		if(MSVC)
 			# Specify unwind semantics so that MSVC knowns how to handle exceptions
 			target_compile_options(${CIFPP_TEST} PRIVATE /EHsc)
@@ -449,10 +443,10 @@ if(CIFPP_BUILD_TESTS)

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

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

 	endforeach()
 endif()
--- a/Config.cmake.in
+++ b/Config.cmake.in
@@ -4,7 +4,6 @@ include(CMakeFindDependencyMacro)
 find_dependency(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_options)
 if(NOT WIN32)
 find_dependency(ZLIB)
-find_dependency(BZip2)
 endif()

 INCLUDE("${CMAKE_CURRENT_LIST_DIR}/cifppTargets.cmake")
--- a/14
+++ b/14
@@ -1,3 +1,17 @@
+Version 3.0.2
+- refactored mmcif::Atom for performance reasons
+
+Version 3.0.1
+- Fixed processing of proline restraints file from CCP4, proline
+  is a peptide, really.
+- Added code to facilitate DSSP
+
+Version 3.0.0
+- Replaced many strings in the API with string_view for
+  performance reasons.
+- Upgraded mmcif::Structure
+- various other small fixes
+
 Version 2.0.4
 - Reverted a too strict test when reading cif files.

--- a/include/cif++/AtomType.hpp
+++ b/include/cif++/AtomType.hpp
@@ -1,17 +1,17 @@
 /*-
 * SPDX-License-Identifier: BSD-2-Clause
- * 
+ *
 * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- * 
+ *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
- * 
+ *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
- * 
+ *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -29,150 +29,151 @@
 #pragma once

 #include <cstdint>
-#include <string>
 #include <stdexcept>
+#include <string>

 namespace mmcif
 {

 enum AtomType : uint8_t
 {
-	Nn = 0,		// Unknown
-	
-	H = 1,		// Hydrogen
-	He = 2,		// Helium
+	Nn = 0, // Unknown

-	Li = 3,		// Lithium
-	Be = 4,		// Beryllium
-	B = 5,		// Boron
-	C = 6,		// Carbon
-	N = 7,		// Nitrogen
-	O = 8,		// Oxygen
-	F = 9,		// Fluorine
-	Ne = 10,	// Neon
+	H = 1,  // Hydrogen
+	He = 2, // Helium

-	Na = 11,	// Sodium
-	Mg = 12,	// Magnesium
-	Al = 13,	// Aluminium
-	Si = 14,	// Silicon
-	P = 15,		// Phosphorus
-	S = 16,		// Sulfur
-	Cl = 17,	// Chlorine
-	Ar = 18,	// Argon
+	Li = 3,  // Lithium
+	Be = 4,  // Beryllium
+	B = 5,   // Boron
+	C = 6,   // Carbon
+	N = 7,   // Nitrogen
+	O = 8,   // Oxygen
+	F = 9,   // Fluorine
+	Ne = 10, // Neon

-	K = 19,		// Potassium
-	Ca = 20,	// Calcium
-	Sc = 21,	// Scandium
-	Ti = 22,	// Titanium
-	V = 23,		// Vanadium
-	Cr = 24,	// Chromium
-	Mn = 25,	// Manganese
-	Fe = 26,	// Iron
-	Co = 27,	// Cobalt
-	Ni = 28,	// Nickel
-	Cu = 29,	// Copper
-	Zn = 30,	// Zinc
-	Ga = 31,	// Gallium
-	Ge = 32,	// Germanium
-	As = 33,	// Arsenic
-	Se = 34,	// Selenium
-	Br = 35,	// Bromine
-	Kr = 36,	// Krypton
+	Na = 11, // Sodium
+	Mg = 12, // Magnesium
+	Al = 13, // Aluminium
+	Si = 14, // Silicon
+	P = 15,  // Phosphorus
+	S = 16,  // Sulfur
+	Cl = 17, // Chlorine
+	Ar = 18, // Argon

-	Rb = 37,	// Rubidium
-	Sr = 38,	// Strontium
-	Y = 39,		// Yttrium
-	Zr = 40,	// Zirconium
-	Nb = 41,	// Niobium
-	Mo = 42,	// Molybdenum
-	Tc = 43,	// Technetium
-	Ru = 44,	// Ruthenium
-	Rh = 45,	// Rhodium
-	Pd = 46,	// Palladium
-	Ag = 47,	// Silver
-	Cd = 48,	// Cadmium
-	In = 49,	// Indium
-	Sn = 50,	// Tin
-	Sb = 51,	// Antimony
-	Te = 52,	// Tellurium
-	I = 53,		// Iodine
-	Xe = 54,	// Xenon
-	Cs = 55,	// Caesium
-	Ba = 56,	// Barium
-	La = 57,	// Lanthanum
+	K = 19,  // Potassium
+	Ca = 20, // Calcium
+	Sc = 21, // Scandium
+	Ti = 22, // Titanium
+	V = 23,  // Vanadium
+	Cr = 24, // Chromium
+	Mn = 25, // Manganese
+	Fe = 26, // Iron
+	Co = 27, // Cobalt
+	Ni = 28, // Nickel
+	Cu = 29, // Copper
+	Zn = 30, // Zinc
+	Ga = 31, // Gallium
+	Ge = 32, // Germanium
+	As = 33, // Arsenic
+	Se = 34, // Selenium
+	Br = 35, // Bromine
+	Kr = 36, // Krypton

-	Hf = 72,	// Hafnium
-	Ta = 73,	// Tantalum
-	W = 74,		// Tungsten
-	Re = 75,	// Rhenium
-	Os = 76,	// Osmium
-	Ir = 77,	// Iridium
-	Pt = 78,	// Platinum
-	Au = 79,	// Gold
-	Hg = 80,	// Mercury
-	Tl = 81,	// Thallium
-	Pb = 82,	// Lead
-	Bi = 83,	// Bismuth
-	Po = 84,	// Polonium
-	At = 85,	// Astatine
-	Rn = 86,	// Radon
-	Fr = 87,	// Francium
-	Ra = 88,	// Radium
-	Ac = 89,	// Actinium
+	Rb = 37, // Rubidium
+	Sr = 38, // Strontium
+	Y = 39,  // Yttrium
+	Zr = 40, // Zirconium
+	Nb = 41, // Niobium
+	Mo = 42, // Molybdenum
+	Tc = 43, // Technetium
+	Ru = 44, // Ruthenium
+	Rh = 45, // Rhodium
+	Pd = 46, // Palladium
+	Ag = 47, // Silver
+	Cd = 48, // Cadmium
+	In = 49, // Indium
+	Sn = 50, // Tin
+	Sb = 51, // Antimony
+	Te = 52, // Tellurium
+	I = 53,  // Iodine
+	Xe = 54, // Xenon
+	Cs = 55, // Caesium
+	Ba = 56, // Barium
+	La = 57, // Lanthanum

-	Rf = 104,	// Rutherfordium
-	Db = 105,	// Dubnium
-	Sg = 106,	// Seaborgium
-	Bh = 107,	// Bohrium
-	Hs = 108,	// Hassium
-	Mt = 109,	// Meitnerium
-	Ds = 110,	// Darmstadtium
-	Rg = 111,	// Roentgenium
-	Cn = 112,	// Copernicium
-	Nh = 113,	// Nihonium
-	Fl = 114,	// Flerovium
-	Mc = 115,	// Moscovium
-	Lv = 116,	// Livermorium
-	Ts = 117,	// Tennessine
-	Og = 118,	// Oganesson
+	Hf = 72, // Hafnium
+	Ta = 73, // Tantalum
+	W = 74,  // Tungsten
+	Re = 75, // Rhenium
+	Os = 76, // Osmium
+	Ir = 77, // Iridium
+	Pt = 78, // Platinum
+	Au = 79, // Gold
+	Hg = 80, // Mercury
+	Tl = 81, // Thallium
+	Pb = 82, // Lead
+	Bi = 83, // Bismuth
+	Po = 84, // Polonium
+	At = 85, // Astatine
+	Rn = 86, // Radon
+	Fr = 87, // Francium
+	Ra = 88, // Radium
+	Ac = 89, // Actinium

-	Ce = 58,	// Cerium
-	Pr = 59,	// Praseodymium
-	Nd = 60,	// Neodymium
-	Pm = 61,	// Promethium
-	Sm = 62,	// Samarium
-	Eu = 63,	// Europium
-	Gd = 64,	// Gadolinium
-	Tb = 65,	// Terbium
-	Dy = 66,	// Dysprosium
-	Ho = 67,	// Holmium
-	Er = 68,	// Erbium
-	Tm = 69,	// Thulium
-	Yb = 70,	// Ytterbium
-	Lu = 71,	// Lutetium
+	Rf = 104, // Rutherfordium
+	Db = 105, // Dubnium
+	Sg = 106, // Seaborgium
+	Bh = 107, // Bohrium
+	Hs = 108, // Hassium
+	Mt = 109, // Meitnerium
+	Ds = 110, // Darmstadtium
+	Rg = 111, // Roentgenium
+	Cn = 112, // Copernicium
+	Nh = 113, // Nihonium
+	Fl = 114, // Flerovium
+	Mc = 115, // Moscovium
+	Lv = 116, // Livermorium
+	Ts = 117, // Tennessine
+	Og = 118, // Oganesson

-	Th = 90,	// Thorium
-	Pa = 91,	// Protactinium
-	U = 92,		// Uranium
-	Np = 93,	// Neptunium
-	Pu = 94,	// Plutonium
-	Am = 95,	// Americium
-	Cm = 96,	// Curium
-	Bk = 97,	// Berkelium
-	Cf = 98,	// Californium
-	Es = 99,	// Einsteinium
-	Fm = 100,	// Fermium
-	Md = 101,	// Mendelevium
-	No = 102,	// Nobelium
-	Lr = 103,	// Lawrencium
+	Ce = 58, // Cerium
+	Pr = 59, // Praseodymium
+	Nd = 60, // Neodymium
+	Pm = 61, // Promethium
+	Sm = 62, // Samarium
+	Eu = 63, // Europium
+	Gd = 64, // Gadolinium
+	Tb = 65, // Terbium
+	Dy = 66, // Dysprosium
+	Ho = 67, // Holmium
+	Er = 68, // Erbium
+	Tm = 69, // Thulium
+	Yb = 70, // Ytterbium
+	Lu = 71, // Lutetium

-	D = 129,	// Deuterium
+	Th = 90,  // Thorium
+	Pa = 91,  // Protactinium
+	U = 92,   // Uranium
+	Np = 93,  // Neptunium
+	Pu = 94,  // Plutonium
+	Am = 95,  // Americium
+	Cm = 96,  // Curium
+	Bk = 97,  // Berkelium
+	Cf = 98,  // Californium
+	Es = 99,  // Einsteinium
+	Fm = 100, // Fermium
+	Md = 101, // Mendelevium
+	No = 102, // Nobelium
+	Lr = 103, // Lawrencium
+
+	D = 129, // Deuterium
 };

 // --------------------------------------------------------------------
 // AtomTypeInfo

-enum RadiusType {
+enum RadiusType
+{
 	eRadiusCalculated,
 	eRadiusEmpirical,
 	eRadiusCovalentEmpirical,
@@ -188,12 +189,12 @@ enum RadiusType {

 struct AtomTypeInfo
 {
-	AtomType		type;
-	std::string		name;
-	std::string		symbol;
-	float			weight;
-	bool			metal;
-	float			radii[eRadiusTypeCount];
+	AtomType type;
+	std::string name;
+	std::string symbol;
+	float weight;
+	bool metal;
+	float radii[eRadiusTypeCount];
 };

 extern const AtomTypeInfo kKnownAtoms[];
@@ -205,25 +206,25 @@ class AtomTypeTraits
 {
  public:
 	AtomTypeTraits(AtomType a);
-	AtomTypeTraits(const std::string& symbol);
-	
-	AtomType type() const			{ return mInfo->type; }
-	std::string	name() const		{ return mInfo->name; }
-	std::string	symbol() const		{ return mInfo->symbol; }
-	float weight() const			{ return mInfo->weight; }
-	
-	bool isMetal() const			{ return mInfo->metal; }
-	
-	static bool isElement(const std::string& symbol);
-	static bool isMetal(const std::string& symbol);
-	
+	AtomTypeTraits(const std::string &symbol);
+
+	AtomType type() const { return mInfo->type; }
+	std::string name() const { return mInfo->name; }
+	std::string symbol() const { return mInfo->symbol; }
+	float weight() const { return mInfo->weight; }
+
+	bool isMetal() const { return mInfo->metal; }
+
+	static bool isElement(const std::string &symbol);
+	static bool isMetal(const std::string &symbol);
+
 	float radius(RadiusType type = eRadiusSingleBond) const
 	{
 		if (type >= eRadiusTypeCount)
 			throw std::invalid_argument("invalid radius requested");
 		return mInfo->radii[type] / 100.f;
 	}
-	
+
 	// data type encapsulating the Waasmaier & Kirfel scattering factors
 	// in a simplified form (only a and b).
 	// Added the electrion scattering factors as well
@@ -231,15 +232,18 @@ class AtomTypeTraits
 	{
 		double a[6], b[6];
 	};
-	
+
 	// to get the Cval and Siva values, use this constant as charge:
-	enum { kWKSFVal = -99 };
-	
-	const SFData& wksf(int charge = 0) const;
-	const SFData& elsf() const;
+	enum
+	{
+		kWKSFVal = -99
+	};
+
+	const SFData &wksf(int charge = 0) const;
+	const SFData &elsf() const;

  private:
-	const struct AtomTypeInfo*	mInfo;
+	const struct AtomTypeInfo *mInfo;
 };

-}
+} // namespace mmcif
--- a/include/cif++/BondMap.hpp
+++ b/include/cif++/BondMap.hpp
@@ -1,17 +1,17 @@
 /*-
 * SPDX-License-Identifier: BSD-2-Clause
- * 
+ *
 * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- * 
+ *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
- * 
+ *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
- * 
+ *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -26,9 +26,9 @@

 #pragma once

-#include <unordered_map>
 #include <filesystem>
 #include <stdexcept>
+#include <unordered_map>

 #include "cif++/Structure.hpp"

@@ -38,39 +38,40 @@ namespace mmcif
 class BondMapException : public std::runtime_error
 {
  public:
-	BondMapException(const std::string& msg)
-		: runtime_error(msg) {}
+	BondMapException(const std::string &msg)
+		: runtime_error(msg)
+	{
+	}
 };

 class BondMap
 {
  public:
-	BondMap(const Structure& p);
-	
-	BondMap(const BondMap&) = delete;
-	BondMap& operator=(const BondMap&) = delete;
+	BondMap(const Structure &p);

-	bool operator()(const Atom& a, const Atom& b) const
+	BondMap(const BondMap &) = delete;
+	BondMap &operator=(const BondMap &) = delete;
+
+	bool operator()(const Atom &a, const Atom &b) const
 	{
 		return isBonded(index.at(a.id()), index.at(b.id()));
 	}

-	bool is1_4(const Atom& a, const Atom& b) const
+	bool is1_4(const Atom &a, const Atom &b) const
 	{
 		uint32_t ixa = index.at(a.id());
 		uint32_t ixb = index.at(b.id());
-	
+
 		return bond_1_4.count(key(ixa, ixb));
 	}
-	
+
 	// links coming from the struct_conn records:
-	std::vector<std::string> linked(const Atom& a) const;
+	std::vector<std::string> linked(const Atom &a) const;

 	// This list of atomID's is comming from either CCD or the CCP4 dictionaries loaded
-	static std::vector<std::string> atomIDsForCompound(const std::string& compoundID);
-	
-  private:
+	static std::vector<std::string> atomIDsForCompound(const std::string &compoundID);

+  private:
 	bool isBonded(uint32_t ai, uint32_t bi) const
 	{
 		return bond.count(key(ai, bi)) != 0;
@@ -82,20 +83,19 @@ class BondMap
 			std::swap(a, b);
 		return static_cast<uint64_t>(a) | (static_cast<uint64_t>(b) << 32);
 	}
-	
-	std::tuple<uint32_t,uint32_t> dekey(uint64_t k) const
+
+	std::tuple<uint32_t, uint32_t> dekey(uint64_t k) const
 	{
 		return std::make_tuple(
 			static_cast<uint32_t>(k >> 32),
-			static_cast<uint32_t>(k)
-		);
+			static_cast<uint32_t>(k));
 	}
-	
+
 	uint32_t dim;
-	std::unordered_map<std::string,uint32_t> index;
+	std::unordered_map<std::string, uint32_t> index;
 	std::set<uint64_t> bond, bond_1_4;

-	std::map<std::string,std::set<std::string>> link;
+	std::map<std::string, std::set<std::string>> link;
 };

-}
+} // namespace mmcif
--- a/include/cif++/Cif++.hpp
+++ b/include/cif++/Cif++.hpp
@@ -36,6 +36,10 @@
 #include <regex>
 #include <set>
 #include <sstream>
+#include <iomanip>
+#include <shared_mutex>
+
+#include <boost/format.hpp>

 #include "cif++/CifUtils.hpp"

@@ -141,14 +145,27 @@ class Item
  public:
 	Item() {}

+	Item(std::string_view name, char value)
+		: mName(name)
+		, mValue({ value })
+	{
+	}
+
+	template<typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
+	Item(std::string_view name, const T& value, const char *fmt)
+		: mName(name)
+		, mValue((boost::format(fmt) % value).str())
+	{
+	}
+
 	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
-	Item(const std::string &name, const T &value)
+	Item(const std::string_view name, const T &value)
 		: mName(name)
 		, mValue(std::to_string(value))
 	{
 	}

-	Item(const std::string &name, const std::string &value)
+	Item(const std::string_view name, const std::string_view value)
 		: mName(name)
 		, mValue(value)
 	{
@@ -221,7 +238,7 @@ class Datablock
 	using iterator = CategoryList::iterator;
 	using const_iterator = CategoryList::const_iterator;

-	Datablock(const std::string &name);
+	Datablock(const std::string_view name);
 	~Datablock();

 	Datablock(const Datablock &) = delete;
@@ -230,33 +247,31 @@ class Datablock
 	std::string getName() const { return mName; }
 	void setName(const std::string &n) { mName = n; }

-	std::string firstItem(const std::string &tag) const;
-
 	iterator begin() { return mCategories.begin(); }
 	iterator end() { return mCategories.end(); }

 	const_iterator begin() const { return mCategories.begin(); }
 	const_iterator end() const { return mCategories.end(); }

-	Category &operator[](const std::string &name);
+	Category &operator[](std::string_view name);

-	std::tuple<iterator, bool> emplace(const std::string &name);
+	std::tuple<iterator, bool> emplace(std::string_view name);

 	bool isValid();
 	void validateLinks() const;

-	void setValidator(Validator *v);
+	void setValidator(const Validator *v);

 	// this one only looks up a Category, returns nullptr if it does not exist
-	const Category *get(const std::string &name) const;
-	Category *get(const std::string &name);
+	const Category *get(std::string_view name) const;
+	Category *get(std::string_view name);

 	void getTagOrder(std::vector<std::string> &tags) const;
 	void write(std::ostream &os, const std::vector<std::string> &order);
 	void write(std::ostream &os);

 	// convenience function, add a line to the software category
-	void add_software(const std::string &name, const std::string &classification,
+	void add_software(const std::string_view name, const std::string &classification,
 		const std::string &versionNr, const std::string &versionDate);

 	friend bool operator==(const Datablock &lhs, const Datablock &rhs);
@@ -264,9 +279,10 @@ class Datablock
 	friend std::ostream& operator<<(std::ostream &os, const Datablock &data);

  private:
-	std::list<Category> mCategories;
+	CategoryList mCategories;		// LRU
+	mutable std::shared_mutex mLock;
 	std::string mName;
-	Validator *mValidator;
+	const Validator *mValidator;
 	Datablock *mNext;
 };

@@ -350,14 +366,14 @@ namespace detail
 	  private:
 		friend class ::cif::Row;

-		ItemReference(const char *name, size_t column, Row &row)
+		ItemReference(std::string_view name, size_t column, Row &row)
 			: mName(name)
 			, mColumn(column)
 			, mRow(row)
 		{
 		}

-		ItemReference(const char *name, size_t column, const Row &row)
+		ItemReference(std::string_view name, size_t column, const Row &row)
 			: mName(name)
 			, mColumn(column)
 			, mRow(const_cast<Row &>(row))
@@ -365,7 +381,7 @@ namespace detail
 		{
 		}

-		const char *mName;
+		std::string_view mName;
 		size_t mColumn;
 		Row &mRow;
 		bool mConst = false;
@@ -746,16 +762,16 @@ class Row
 		return detail::ItemReference(itemTag, column, *this);
 	}

-	const detail::ItemReference operator[](const std::string &itemTag) const
+	const detail::ItemReference operator[](std::string_view itemTag) const
 	{
-		size_t column = ColumnForItemTag(itemTag.c_str());
-		return detail::ItemReference(itemTag.c_str(), column, *this);
+		size_t column = ColumnForItemTag(itemTag);
+		return detail::ItemReference(itemTag, column, *this);
 	}

-	detail::ItemReference operator[](const std::string &itemTag)
+	detail::ItemReference operator[](std::string_view itemTag)
 	{
-		size_t column = ColumnForItemTag(itemTag.c_str());
-		return detail::ItemReference(itemTag.c_str(), column, *this);
+		size_t column = ColumnForItemTag(itemTag);
+		return detail::ItemReference(itemTag, column, *this);
 	}

 	template <typename... Ts, size_t N>
@@ -776,7 +792,7 @@ class Row
 	}

 	void assign(const std::vector<Item> &values);
-	void assign(const std::string &name, const std::string &value, bool updateLinked);
+	void assign(std::string_view name, const std::string &value, bool updateLinked, bool validate = true);

 	bool operator==(const Row &rhs) const
 	{
@@ -798,12 +814,12 @@ class Row
 	friend std::ostream &operator<<(std::ostream &os, const Row &row);

  private:
-	void assign(size_t column, const std::string &value, bool updateLinked);
+	void assign(size_t column, const std::string &value, bool updateLinked, bool validate = true);
 	void assign(const Item &i, bool updateLinked);

 	static void swap(size_t column, ItemRow *a, ItemRow *b);

-	size_t ColumnForItemTag(const char *itemTag) const;
+	size_t ColumnForItemTag(std::string_view itemTag) const;

 	mutable ItemRow *mData;
 	uint32_t mLineNr = 0;
@@ -1169,11 +1185,11 @@ struct Empty

 struct Key
 {
-	Key(const std::string &itemTag)
+	explicit Key(const std::string &itemTag)
 		: mItemTag(itemTag)
 	{
 	}
-	Key(const char *itemTag)
+	explicit Key(const char *itemTag)
 		: mItemTag(itemTag)
 	{
 	}
@@ -1816,12 +1832,12 @@ class Category
 	friend class Row;
 	friend class detail::ItemReference;

-	Category(Datablock &db, const std::string &name, Validator *Validator);
+	Category(Datablock &db, const std::string_view name, const Validator *Validator);
 	Category(const Category &) = delete;
 	Category &operator=(const Category &) = delete;
 	~Category();

-	const std::string name() const { return mName; }
+	const std::string &name() const { return mName; }

 	using iterator = iterator_impl<Row>;
 	using const_iterator = iterator_impl<const Row>;
@@ -2064,7 +2080,7 @@ class Category

 	Datablock &db() { return mDb; }

-	void setValidator(Validator *v);
+	void setValidator(const Validator *v);

 	iset fields() const;
 	iset mandatoryFields() const;
@@ -2077,8 +2093,8 @@ class Category
 	void getTagOrder(std::vector<std::string> &tags) const;

 	// return index for known column, or the next available column index
-	size_t getColumnIndex(const std::string &name) const;
-	bool hasColumn(const std::string &name) const;
+	size_t getColumnIndex(std::string_view name) const;
+	bool hasColumn(std::string_view name) const;
 	const std::string &getColumnName(size_t columnIndex) const;
 	std::vector<std::string> getColumnNames() const;

@@ -2119,16 +2135,27 @@ class Category
 	void write(std::ostream &os, const std::vector<std::string> &order);
 	void write(std::ostream &os, const std::vector<size_t> &order, bool includeEmptyColumns);

-	size_t addColumn(const std::string &name);
+	size_t addColumn(std::string_view name);
+
+	struct Linked
+	{
+		Category *linked;
+		const ValidateLink *v;
+	};
+
+	void updateLinks();

 	Datablock &mDb;
 	std::string mName;
-	Validator *mValidator;
+	const Validator *mValidator;
 	const ValidateCategory *mCatValidator = nullptr;
 	std::vector<ItemColumn> mColumns;
 	ItemRow *mHead;
 	ItemRow *mTail;
+	size_t mLastUniqueNr = 0;
 	class CatIndex *mIndex;
+
+	std::vector<Linked> mParentLinks, mChildLinks;
 };

 // --------------------------------------------------------------------
@@ -2162,7 +2189,8 @@ class File

 	void loadDictionary();                 // load the default dictionary, that is mmcifDdl in this case
 	void loadDictionary(const char *dict); // load one of the compiled in dictionaries
-	void loadDictionary(std::istream &is); // load dictionary from input stream
+
+	void setValidator(const Validator *v);

 	bool isValid();
 	void validateLinks() const;
@@ -2183,8 +2211,8 @@ class File

 	void append(Datablock *e);

-	Datablock *get(const std::string &name) const;
-	Datablock &operator[](const std::string &name);
+	Datablock *get(std::string_view name) const;
+	Datablock &operator[](std::string_view name);

 	struct iterator
 	{
@@ -2226,10 +2254,8 @@ class File
 	void getTagOrder(std::vector<std::string> &tags) const;

  private:
-	void setValidator(Validator *v);
-
 	Datablock *mHead;
-	Validator *mValidator;
+	const Validator *mValidator;
 };

 // --------------------------------------------------------------------
--- a/include/cif++/CifParser.hpp
+++ b/include/cif++/CifParser.hpp
@@ -28,8 +28,8 @@

 #include "cif++/Cif++.hpp"

-#include <stack>
 #include <map>
+#include <stack>

 namespace cif
 {
@@ -39,7 +39,7 @@ namespace cif
 class CifParserError : public std::runtime_error
 {
  public:
-	CifParserError(uint32_t lineNr, const std::string& message);
+	CifParserError(uint32_t lineNr, const std::string &message);
 };

 // --------------------------------------------------------------------
@@ -48,7 +48,8 @@ extern const uint32_t kMaxLineLength;

 extern const uint8_t kCharTraitsTable[128];

-enum CharTraitsMask: uint8_t {
+enum CharTraitsMask : uint8_t
+{
 	kOrdinaryMask = 1 << 0,
 	kNonBlankMask = 1 << 1,
 	kTextLeadMask = 1 << 2,
@@ -75,13 +76,13 @@ inline bool isTextLead(int ch)
 	return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0;
 }

-inline bool isAnyPrint(int ch)	
+inline bool isAnyPrint(int ch)
 {
-	return ch == '\t' or 
-		(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
+	return ch == '\t' or
+	       (ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
 }

-inline bool isUnquotedString(const char* s)
+inline bool isUnquotedString(const char *s)
 {
 	bool result = isOrdinary(*s++);
 	while (result and *s != 0)
@@ -94,11 +95,7 @@ inline bool isUnquotedString(const char* s)

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

-std::tuple<std::string,std::string> splitTagName(const std::string& tag);
-
-// --------------------------------------------------------------------
-
-using DatablockIndex = std::map<std::string,std::size_t>;
+using DatablockIndex = std::map<std::string, std::size_t>;

 // --------------------------------------------------------------------
 // sac Parser, analogous to SAX Parser (simple api for xml)
@@ -106,15 +103,15 @@ using DatablockIndex = std::map<std::string,std::size_t>;
 class SacParser
 {
  public:
-	SacParser(std::istream& is, bool init = true);
+	SacParser(std::istream &is, bool init = true);
 	virtual ~SacParser() {}

 	enum CIFToken
 	{
 		eCIFTokenUnknown,
-		
+
 		eCIFTokenEOF,
-	
+
 		eCIFTokenDATA,
 		eCIFTokenLOOP,
 		eCIFTokenGLOBAL,
@@ -124,7 +121,7 @@ class SacParser
 		eCIFTokenValue,
 	};

-	static const char* kTokenName[];
+	static const char *kTokenName[];

 	enum CIFValueType
 	{
@@ -137,40 +134,39 @@ class SacParser
 		eCIFValueUnknown
 	};

-	static const char* kValueName[];
-	
+	static const char *kValueName[];
+
 	int getNextChar();

 	void retract();
 	void restart();
-	
+
 	CIFToken getNextToken();
 	void match(CIFToken token);

-	bool parseSingleDatablock(const std::string& datablock);
+	bool parseSingleDatablock(const std::string &datablock);

 	DatablockIndex indexDatablocks();
-	bool parseSingleDatablock(const std::string& datablock, const DatablockIndex &index);
+	bool parseSingleDatablock(const std::string &datablock, const DatablockIndex &index);

 	void parseFile();
 	void parseGlobal();
 	void parseDataBlock();

 	virtual void parseSaveFrame();
-	
+
 	void parseDictionary();
-	
-	void error(const std::string& msg);
-	
+
+	void error(const std::string &msg);
+
 	// production methods, these are pure virtual here
-	
-	virtual void produceDatablock(const std::string& name) = 0;
-	virtual void produceCategory(const std::string& name) = 0;
+
+	virtual void produceDatablock(const std::string &name) = 0;
+	virtual void produceCategory(const std::string &name) = 0;
 	virtual void produceRow() = 0;
-	virtual void produceItem(const std::string& category, const std::string& item, const std::string& value) = 0;
+	virtual void produceItem(const std::string &category, const std::string &item, const std::string &value) = 0;

  protected:
-
 	enum State
 	{
 		eStateStart,
@@ -185,21 +181,21 @@ class SacParser
 		eStateTextField,
 		eStateFloat = 100,
 		eStateInt = 110,
-//		eStateNumericSuffix = 200,
+		//		eStateNumericSuffix = 200,
 		eStateValue = 300
 	};

-	std::istream&			mData;
+	std::istream &mData;

 	// Parser state
-	bool					mValidate;
-	uint32_t					mLineNr;
-	bool					mBol;
-	int						mState, mStart;
-	CIFToken				mLookahead;
-	std::string				mTokenValue;
-	CIFValueType			mTokenType;
-	std::stack<int>			mBuffer;
+	bool mValidate;
+	uint32_t mLineNr;
+	bool mBol;
+	int mState, mStart;
+	CIFToken mLookahead;
+	std::string mTokenValue;
+	CIFValueType mTokenType;
+	std::stack<int> mBuffer;
 };

 // --------------------------------------------------------------------
@@ -207,18 +203,18 @@ class SacParser
 class Parser : public SacParser
 {
  public:
-	Parser(std::istream& is, File& f, bool init = true);
+	Parser(std::istream &is, File &f, bool init = true);

-	virtual void produceDatablock(const std::string& name);
-	virtual void produceCategory(const std::string& name);
+	virtual void produceDatablock(const std::string &name);
+	virtual void produceCategory(const std::string &name);
 	virtual void produceRow();
-	virtual void produceItem(const std::string& category, const std::string& item, const std::string& value);
+	virtual void produceItem(const std::string &category, const std::string &item, const std::string &value);

  protected:
-	File&					mFile;
-	Datablock*				mDataBlock;
-	Datablock::iterator		mCat;
-	Row						mRow;
+	File &mFile;
+	Datablock *mDataBlock;
+	Datablock::iterator mCat;
+	Row mRow;
 };

 // --------------------------------------------------------------------
@@ -226,23 +222,21 @@ class Parser : public SacParser
 class DictParser : public Parser
 {
  public:
-
-	DictParser(Validator& validator, std::istream& is);
+	DictParser(Validator &validator, std::istream &is);
 	~DictParser();
-	
-	void loadDictionary();
-	
-  private:

+	void loadDictionary();
+
+  private:
 	virtual void parseSaveFrame();
-	
+
 	bool collectItemTypes();
 	void linkItems();

-	Validator&						mValidator;
-	File							mFile;
-	struct DictParserDataImpl*		mImpl;
-	bool							mCollectedItemTypes = false;
+	Validator &mValidator;
+	File mFile;
+	struct DictParserDataImpl *mImpl;
+	bool mCollectedItemTypes = false;
 };

-}
+} // namespace cif
--- a/include/cif++/CifUtils.hpp
+++ b/include/cif++/CifUtils.hpp
@@ -67,8 +67,8 @@ std::string get_version_nr();
 // some basic utilities: Since we're using ASCII input only, we define for optimisation
 // our own case conversion routines.

-bool iequals(const std::string &a, const std::string &b);
-int icompare(const std::string &a, const std::string &b);
+bool iequals(std::string_view a, std::string_view b);
+int icompare(std::string_view a, std::string_view b);

 bool iequals(const char *a, const char *b);
 int icompare(const char *a, const char *b);
@@ -100,7 +100,7 @@ inline char tolower(int ch)

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

-std::tuple<std::string, std::string> splitTagName(const std::string &tag);
+std::tuple<std::string, std::string> splitTagName(std::string_view tag);

 // --------------------------------------------------------------------
 // generate a cif name, mainly used to generate asym_id's
--- a/include/cif++/CifValidator.hpp
+++ b/include/cif++/CifValidator.hpp
@@ -36,18 +36,19 @@

 namespace cif
 {
-	
+
 struct ValidateCategory;
+class ValidatorFactory;

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

 class ValidationError : public std::exception
 {
  public:
-	ValidationError(const std::string& msg);
-	ValidationError(const std::string& cat, const std::string& item,
-		const std::string& msg);
-	const char* what() const noexcept		{ return mMsg.c_str(); }
+	ValidationError(const std::string &msg);
+	ValidationError(const std::string &cat, const std::string &item,
+		const std::string &msg);
+	const char *what() const noexcept { return mMsg.c_str(); }
 	std::string mMsg;
 };

@@ -55,58 +56,60 @@ class ValidationError : public std::exception

 enum class DDL_PrimitiveType
 {
-	Char, UChar, Numb
+	Char,
+	UChar,
+	Numb
 };

-DDL_PrimitiveType mapToPrimitiveType(const std::string& s);
+DDL_PrimitiveType mapToPrimitiveType(std::string_view s);

 struct ValidateType
 {
-	std::string				mName;
-	DDL_PrimitiveType		mPrimitiveType;
+	std::string mName;
+	DDL_PrimitiveType mPrimitiveType;
 	// std::regex				mRx;
-	boost::regex			mRx;
+	boost::regex mRx;

-	bool operator<(const ValidateType& rhs) const
+	bool operator<(const ValidateType &rhs) const
 	{
 		return icompare(mName, rhs.mName) < 0;
 	}

-	// compare values based on type	
-//	int compare(const std::string& a, const std::string& b) const
-//	{
-//		return compare(a.c_str(), b.c_str());
-//	}
-	
-	int compare(const char* a, const char* b) const;
+	// compare values based on type
+	//	int compare(const std::string& a, const std::string& b) const
+	//	{
+	//		return compare(a.c_str(), b.c_str());
+	//	}
+
+	int compare(const char *a, const char *b) const;
 };

 struct ValidateItem
 {
-	std::string				mTag;
-	bool					mMandatory;
-	const ValidateType*		mType;
-	cif::iset				mEnums;
-	std::string				mDefault;
-	bool					mDefaultIsNull;
-	ValidateCategory*		mCategory = nullptr;
+	std::string mTag;
+	bool mMandatory;
+	const ValidateType *mType;
+	cif::iset mEnums;
+	std::string mDefault;
+	bool mDefaultIsNull;
+	ValidateCategory *mCategory = nullptr;

 	// ItemLinked is used for non-key links
 	struct ItemLinked
 	{
-		ValidateItem*		mParent;
-		std::string			mParentItem;
-		std::string			mChildItem;
+		ValidateItem *mParent;
+		std::string mParentItem;
+		std::string mChildItem;
 	};

-	std::vector<ItemLinked>	mLinked;
-	
-	bool operator<(const ValidateItem& rhs) const
+	std::vector<ItemLinked> mLinked;
+
+	bool operator<(const ValidateItem &rhs) const
 	{
 		return icompare(mTag, rhs.mTag) < 0;
 	}

-	bool operator==(const ValidateItem& rhs) const
+	bool operator==(const ValidateItem &rhs) const
 	{
 		return iequals(mTag, rhs.mTag);
 	}
@@ -116,22 +119,22 @@ struct ValidateItem

 struct ValidateCategory
 {
-	std::string					mName;
-	std::vector<std::string>	mKeys;
-	cif::iset					mGroups;
-	cif::iset					mMandatoryFields;
-	std::set<ValidateItem>		mItemValidators;
+	std::string mName;
+	std::vector<std::string> mKeys;
+	cif::iset mGroups;
+	cif::iset mMandatoryFields;
+	std::set<ValidateItem> mItemValidators;

-	bool operator<(const ValidateCategory& rhs) const
+	bool operator<(const ValidateCategory &rhs) const
 	{
 		return icompare(mName, rhs.mName) < 0;
 	}

-	void addItemValidator(ValidateItem&& v);
-	
-	const ValidateItem* getValidatorForItem(std::string tag) const;
-	
-	const std::set<ValidateItem>& itemValidators() const
+	void addItemValidator(ValidateItem &&v);
+
+	const ValidateItem *getValidatorForItem(std::string_view tag) const;
+
+	const std::set<ValidateItem> &itemValidators() const
 	{
 		return mItemValidators;
 	}
@@ -139,12 +142,12 @@ struct ValidateCategory

 struct ValidateLink
 {
-	int							mLinkGroupID;
-	std::string					mParentCategory;
-	std::vector<std::string>	mParentKeys;
-	std::string					mChildCategory;
-	std::vector<std::string>	mChildKeys;
-	std::string					mLinkGroupLabel;
+	int mLinkGroupID;
+	std::string mParentCategory;
+	std::vector<std::string> mParentKeys;
+	std::string mChildCategory;
+	std::vector<std::string> mChildKeys;
+	std::string mLinkGroupLabel;
 };

 // --------------------------------------------------------------------
@@ -152,47 +155,72 @@ struct ValidateLink
 class Validator
 {
  public:
-	friend class DictParser;

-	Validator();
+	Validator(std::string_view name, std::istream &is);
 	~Validator();

-	Validator(const Validator& rhs) = delete;
-	Validator& operator=(const Validator& rhs) = delete;
-	
-	Validator(Validator&& rhs);
-	Validator& operator=(Validator&& rhs);
-	
-	void addTypeValidator(ValidateType&& v);
-	const ValidateType* getValidatorForType(std::string typeCode) const;
+	Validator(const Validator &rhs) = delete;
+	Validator &operator=(const Validator &rhs) = delete;

-	void addCategoryValidator(ValidateCategory&& v);
-	const ValidateCategory* getValidatorForCategory(std::string category) const;
+	Validator(Validator &&rhs);
+	Validator &operator=(Validator &&rhs);

-	void addLinkValidator(ValidateLink&& v);
-	std::vector<const ValidateLink*> getLinksForParent(const std::string& category) const;
-	std::vector<const ValidateLink*> getLinksForChild(const std::string& category) const;
+	friend class DictParser;
+	friend class ValidatorFactory;

-	void reportError(const std::string& msg, bool fatal);
-	
-	std::string dictName() const					{ return mName; }
-	void dictName(const std::string& name)			{ mName = name; }
+	void addTypeValidator(ValidateType &&v);
+	const ValidateType *getValidatorForType(std::string_view typeCode) const;

-	std::string dictVersion() const				{ return mVersion; }
-	void dictVersion(const std::string& version)	{ mVersion = version; }
+	void addCategoryValidator(ValidateCategory &&v);
+	const ValidateCategory *getValidatorForCategory(std::string_view category) const;
+
+	void addLinkValidator(ValidateLink &&v);
+	std::vector<const ValidateLink *> getLinksForParent(std::string_view category) const;
+	std::vector<const ValidateLink *> getLinksForChild(std::string_view category) const;
+
+	void reportError(const std::string &msg, bool fatal) const;
+
+	std::string dictName() const { return mName; }
+	void dictName(const std::string &name) { mName = name; }
+
+	std::string dictVersion() const { return mVersion; }
+	void dictVersion(const std::string &version) { mVersion = version; }

  private:

 	// name is fully qualified here:
-	ValidateItem* getValidatorForItem(std::string name) const;
+	ValidateItem *getValidatorForItem(std::string_view name) const;

-	std::string					mName;
-	std::string					mVersion;
-	bool						mStrict = false;
-//	std::set<uint32_t>			mSubCategories;
-	std::set<ValidateType>		mTypeValidators;
-	std::set<ValidateCategory>	mCategoryValidators;
-	std::vector<ValidateLink>	mLinkValidators;
+	std::string mName;
+	std::string mVersion;
+	bool mStrict = false;
+	//	std::set<uint32_t>			mSubCategories;
+	std::set<ValidateType> mTypeValidators;
+	std::set<ValidateCategory> mCategoryValidators;
+	std::vector<ValidateLink> mLinkValidators;
 };

-}
+// --------------------------------------------------------------------
+
+class ValidatorFactory
+{
+  public:
+
+	static ValidatorFactory &instance()
+	{
+		return sInstance;
+	}
+
+	const Validator &operator[](std::string_view dictionary);
+
+  private:
+
+	static ValidatorFactory sInstance;
+
+	ValidatorFactory();
+
+	std::mutex mMutex;
+	std::list<Validator> mValidators;
+};
+
+} // namespace cif
--- a/include/cif++/Compound.hpp
+++ b/include/cif++/Compound.hpp
@@ -104,6 +104,7 @@ class Compound
 	std::string id() const { return mID; }
 	std::string name() const { return mName; }
 	std::string type() const { return mType; }
+	std::string group() const { return mGroup; }
 	std::string formula() const { return mFormula; }
 	float formulaWeight() const { return mFormulaWeight; }
 	int formalCharge() const { return mFormalCharge; }
@@ -130,11 +131,12 @@ class Compound
 	friend class CCP4CompoundFactoryImpl;

 	Compound(cif::Datablock &db);
-	Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type);
+	Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group);

 	std::string mID;
 	std::string mName;
 	std::string mType;
+	std::string mGroup;
 	std::string mFormula;
 	float mFormulaWeight = 0;
 	int mFormalCharge = 0;
--- a/include/cif++/Matrix.hpp
+++ b/include/cif++/Matrix.hpp
@@ -1,391 +0,0 @@
-/*-
- * SPDX-License-Identifier: BSD-2-Clause
- * 
- * Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
- * 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.
- */
-
-// --------------------------------------------------------------------
-// uBlas compatible matrix types
-
-#pragma once
-
-#include <iostream>
-#include <vector>
-
-// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
-// element m i,j is mapped to [i * n + j] and thus storage is row major
-
-template <typename T>
-class MatrixBase
-{
-  public:
-	using value_type = T;
-
-	virtual ~MatrixBase() {}
-
-	virtual uint32_t dim_m() const = 0;
-	virtual uint32_t dim_n() const = 0;
-
-	virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
-	virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
-
-	MatrixBase &operator*=(const value_type &rhs);
-
-	MatrixBase &operator-=(const value_type &rhs);
-};
-
-template <typename T>
-MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
-{
-	for (uint32_t i = 0; i < dim_m(); ++i)
-	{
-		for (uint32_t j = 0; j < dim_n(); ++j)
-		{
-			operator()(i, j) *= rhs;
-		}
-	}
-
-	return *this;
-}
-
-template <typename T>
-MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
-{
-	for (uint32_t i = 0; i < dim_m(); ++i)
-	{
-		for (uint32_t j = 0; j < dim_n(); ++j)
-		{
-			operator()(i, j) -= rhs;
-		}
-	}
-
-	return *this;
-}
-
-template <typename T>
-std::ostream &operator<<(std::ostream &lhs, const MatrixBase<T> &rhs)
-{
-	lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
-	for (uint32_t i = 0; i < rhs.dim_m(); ++i)
-	{
-		lhs << '(';
-		for (uint32_t j = 0; j < rhs.dim_n(); ++j)
-		{
-			if (j > 0)
-				lhs << ',';
-			lhs << rhs(i, j);
-		}
-		lhs << ')';
-	}
-	lhs << ')';
-
-	return lhs;
-}
-
-template <typename T>
-class Matrix : public MatrixBase<T>
-{
-  public:
-	using value_type = T;
-
-	template <typename T2>
-	Matrix(const MatrixBase<T2> &m)
-		: m_m(m.dim_m())
-		, m_n(m.dim_n())
-	{
-		m_data = new value_type[m_m * m_n];
-		for (uint32_t i = 0; i < m_m; ++i)
-		{
-			for (uint32_t j = 0; j < m_n; ++j)
-				operator()(i, j) = m(i, j);
-		}
-	}
-
-	Matrix()
-		: m_data(nullptr)
-		, m_m(0)
-		, m_n(0)
-	{
-	}
-
-	Matrix(const Matrix &m)
-		: m_m(m.m_m)
-		, m_n(m.m_n)
-	{
-		m_data = new value_type[m_m * m_n];
-		std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
-	}
-
-	Matrix &operator=(const Matrix &m)
-	{
-		value_type *t = new value_type[m.m_m * m.m_n];
-		std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
-
-		delete[] m_data;
-		m_data = t;
-		m_m = m.m_m;
-		m_n = m.m_n;
-
-		return *this;
-	}
-
-	Matrix(uint32_t m, uint32_t n, T v = T())
-		: m_m(m)
-		, m_n(n)
-	{
-		m_data = new value_type[m_m * m_n];
-		std::fill(m_data, m_data + (m_m * m_n), v);
-	}
-
-	virtual ~Matrix()
-	{
-		delete[] m_data;
-	}
-
-	virtual uint32_t dim_m() const { return m_m; }
-	virtual uint32_t dim_n() const { return m_n; }
-
-	virtual value_type operator()(uint32_t i, uint32_t j) const
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-	virtual value_type &operator()(uint32_t i, uint32_t j)
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-	template <typename Func>
-	void each(Func f)
-	{
-		for (uint32_t i = 0; i < m_m * m_n; ++i)
-			f(m_data[i]);
-	}
-
-	template <typename U>
-	Matrix &operator/=(U v)
-	{
-		for (uint32_t i = 0; i < m_m * m_n; ++i)
-			m_data[i] /= v;
-
-		return *this;
-	}
-
-  private:
-	value_type *m_data;
-	uint32_t m_m, m_n;
-};
-
-// --------------------------------------------------------------------
-
-template <typename T>
-class SymmetricMatrix : public MatrixBase<T>
-{
-  public:
-	typedef typename MatrixBase<T>::value_type value_type;
-
-	SymmetricMatrix(uint32_t n, T v = T())
-		: m_owner(true)
-		, m_n(n)
-	{
-		uint32_t N = (m_n * (m_n + 1)) / 2;
-		m_data = new value_type[N];
-		std::fill(m_data, m_data + N, v);
-	}
-
-	SymmetricMatrix(const T *data, uint32_t n)
-		: m_owner(false)
-		, m_data(const_cast<T *>(data))
-		, m_n(n)
-	{
-	}
-
-	virtual ~SymmetricMatrix()
-	{
-		if (m_owner)
-			delete[] m_data;
-	}
-
-	virtual uint32_t dim_m() const { return m_n; }
-	virtual uint32_t dim_n() const { return m_n; }
-
-	T operator()(uint32_t i, uint32_t j) const;
-	virtual T &operator()(uint32_t i, uint32_t j);
-
-	// erase two rows, add one at the end (for neighbour joining)
-	void erase_2(uint32_t i, uint32_t j);
-
-	template <typename Func>
-	void each(Func f)
-	{
-		uint32_t N = (m_n * (m_n + 1)) / 2;
-
-		for (uint32_t i = 0; i < N; ++i)
-			f(m_data[i]);
-	}
-
-	template <typename U>
-	SymmetricMatrix &operator/=(U v)
-	{
-		uint32_t N = (m_n * (m_n + 1)) / 2;
-
-		for (uint32_t i = 0; i < N; ++i)
-			m_data[i] /= v;
-
-		return *this;
-	}
-
-  private:
-	bool m_owner;
-	value_type *m_data;
-	uint32_t m_n;
-};
-
-template <typename T>
-inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
-{
-	return i < j
-	           ? m_data[(j * (j + 1)) / 2 + i]
-	           : m_data[(i * (i + 1)) / 2 + j];
-}
-
-template <typename T>
-inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
-{
-	if (i > j)
-		std::swap(i, j);
-	assert(j < m_n);
-	return m_data[(j * (j + 1)) / 2 + i];
-}
-
-template <typename T>
-void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
-{
-	uint32_t s = 0, d = 0;
-	for (uint32_t i = 0; i < m_n; ++i)
-	{
-		for (uint32_t j = 0; j < i; ++j)
-		{
-			if (i != di and j != dj and i != dj and j != di)
-			{
-				if (s != d)
-					m_data[d] = m_data[s];
-				++d;
-			}
-
-			++s;
-		}
-	}
-
-	--m_n;
-}
-
-template <typename T>
-class IdentityMatrix : public MatrixBase<T>
-{
-  public:
-	typedef typename MatrixBase<T>::value_type value_type;
-
-	IdentityMatrix(uint32_t n)
-		: m_n(n)
-	{
-	}
-
-	virtual uint32_t dim_m() const { return m_n; }
-	virtual uint32_t dim_n() const { return m_n; }
-
-	virtual value_type operator()(uint32_t i, uint32_t j) const
-	{
-		value_type result = 0;
-		if (i == j)
-			result = 1;
-		return result;
-	}
-
-  private:
-	uint32_t m_n;
-};
-
-// --------------------------------------------------------------------
-// matrix functions
-
-template <typename T>
-Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
-{
-	Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
-
-	for (uint32_t i = 0; i < result.dim_m(); ++i)
-	{
-		for (uint32_t j = 0; j < result.dim_n(); ++j)
-		{
-			for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
-				result(i, j) += lhs(li, j) * rhs(i, rj);
-		}
-	}
-
-	return result;
-}
-
-template <typename T>
-Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
-{
-	Matrix<T> result(lhs);
-	result *= rhs;
-
-	return result;
-}
-
-template <typename T>
-Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
-{
-	Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
-
-	for (uint32_t i = 0; i < result.dim_m(); ++i)
-	{
-		for (uint32_t j = 0; j < result.dim_n(); ++j)
-		{
-			result(i, j) = lhs(i, j) - rhs(i, j);
-		}
-	}
-
-	return result;
-}
-
-template <typename T>
-Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
-{
-	Matrix<T> result(lhs.dim_m(), lhs.dim_n());
-	result -= rhs;
-	return result;
-}
-
-// template <typename T>
-// symmetric_matrix<T> hammingDistance(const MatrixBase<T> &lhs, T rhs);
-
-// template <typename T>
-// std::vector<T> sum(const MatrixBase<T> &m);
--- a/include/cif++/Point.hpp
+++ b/include/cif++/Point.hpp
@@ -1,17 +1,17 @@
 /*-
 * SPDX-License-Identifier: BSD-2-Clause
- * 
+ *
 * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- * 
+ *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
- * 
+ *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
- * 
+ *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -37,7 +37,7 @@
 namespace mmcif
 {

-typedef boost::math::quaternion<float>	Quaternion;
+typedef boost::math::quaternion<float> Quaternion;

 const double
 	kPI = 3.141592653589793238462643383279502884;
@@ -51,26 +51,43 @@ const double
 //	float x, y, z;
 //	tie(x, y, z) = atom.loc();

-template<typename F>
+template <typename F>
 struct PointF
 {
 	typedef F FType;

 	FType mX, mY, mZ;
-	
-	PointF()							: mX(0), mY(0), mZ(0) {}
-	PointF(FType x, FType y, FType z)	: mX(x), mY(y), mZ(z) {}

-	template<typename PF>
-	PointF(const PointF<PF>& pt)
+	PointF()
+		: mX(0)
+		, mY(0)
+		, mZ(0)
+	{
+	}
+	PointF(FType x, FType y, FType z)
+		: mX(x)
+		, mY(y)
+		, mZ(z)
+	{
+	}
+
+	template <typename PF>
+	PointF(const PointF<PF> &pt)
 		: mX(static_cast<F>(pt.mX))
 		, mY(static_cast<F>(pt.mY))
-		, mZ(static_cast<F>(pt.mZ)) {}
+		, mZ(static_cast<F>(pt.mZ))
+	{
+	}

-#if HAVE_LIBCLIPPER	
-	PointF(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {}
+#if HAVE_LIBCLIPPER
+	PointF(const clipper::Coord_orth &pt)
+		: mX(pt[0])
+		, mY(pt[1])
+		, mZ(pt[2])
+	{
+	}

-	PointF& operator=(const clipper::Coord_orth& rhs)
+	PointF &operator=(const clipper::Coord_orth &rhs)
 	{
 		mX = rhs[0];
 		mY = rhs[1];
@@ -79,72 +96,72 @@ struct PointF
 	}
 #endif

-	template<typename PF>
-	PointF& operator=(const PointF<PF>& rhs)
+	template <typename PF>
+	PointF &operator=(const PointF<PF> &rhs)
 	{
 		mX = static_cast<F>(rhs.mX);
 		mY = static_cast<F>(rhs.mY);
 		mZ = static_cast<F>(rhs.mZ);
 		return *this;
 	}
-	
-	FType& getX()			{ return mX; }
-	FType getX() const		{ return mX; }
-	void setX(FType x)		{ mX = x; }

-	FType& getY()			{ return mY; }
-	FType getY() const		{ return mY; }
-	void setY(FType y)		{ mY = y; }
+	FType &getX() { return mX; }
+	FType getX() const { return mX; }
+	void setX(FType x) { mX = x; }

-	FType& getZ()			{ return mZ; }
-	FType getZ() const		{ return mZ; }
-	void setZ(FType z)		{ mZ = z; }
-	
-	PointF& operator+=(const PointF& rhs)
+	FType &getY() { return mY; }
+	FType getY() const { return mY; }
+	void setY(FType y) { mY = y; }
+
+	FType &getZ() { return mZ; }
+	FType getZ() const { return mZ; }
+	void setZ(FType z) { mZ = z; }
+
+	PointF &operator+=(const PointF &rhs)
 	{
 		mX += rhs.mX;
 		mY += rhs.mY;
 		mZ += rhs.mZ;
-		
+
 		return *this;
 	}
-	
-	PointF& operator+=(FType d)
+
+	PointF &operator+=(FType d)
 	{
 		mX += d;
 		mY += d;
 		mZ += d;
-		
+
 		return *this;
 	}

-	PointF& operator-=(const PointF& rhs)
+	PointF &operator-=(const PointF &rhs)
 	{
 		mX -= rhs.mX;
 		mY -= rhs.mY;
 		mZ -= rhs.mZ;
-		
+
 		return *this;
 	}

-	PointF& operator-=(FType d)
+	PointF &operator-=(FType d)
 	{
 		mX -= d;
 		mY -= d;
 		mZ -= d;
-		
+
 		return *this;
 	}

-	PointF& operator*=(FType rhs)
+	PointF &operator*=(FType rhs)
 	{
 		mX *= rhs;
 		mY *= rhs;
 		mZ *= rhs;
 		return *this;
 	}
-	
-	PointF& operator/=(FType rhs)
+
+	PointF &operator/=(FType rhs)
 	{
 		mX /= rhs;
 		mY /= rhs;
@@ -162,18 +179,18 @@ struct PointF
 		}
 		return length;
 	}
-	
-	void rotate(const boost::math::quaternion<FType>& q)
+
+	void rotate(const boost::math::quaternion<FType> &q)
 	{
 		boost::math::quaternion<FType> p(0, mX, mY, mZ);
-		
+
 		p = q * p * boost::math::conj(q);
-	
+
 		mX = p.R_component_2();
 		mY = p.R_component_3();
 		mZ = p.R_component_4();
 	}
-	
+
 #if HAVE_LIBCLIPPER
 	operator clipper::Coord_orth() const
 	{
@@ -181,21 +198,21 @@ struct PointF
 	}
 #endif

-	operator std::tuple<const FType&, const FType&, const FType&>() const
+	operator std::tuple<const FType &, const FType &, const FType &>() const
 	{
 		return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
 	}

-	operator std::tuple<FType&,FType&,FType&>()
+	operator std::tuple<FType &, FType &, FType &>()
 	{
 		return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
 	}
-	
-	bool operator==(const PointF& rhs) const
+
+	bool operator==(const PointF &rhs) const
 	{
 		return mX == rhs.mX and mY == rhs.mY and mZ == rhs.mZ;
 	}
-	
+
 	// consider point as a vector... perhaps I should rename Point?
 	FType lengthsq() const
 	{
@@ -211,45 +228,45 @@ struct PointF
 typedef PointF<float> Point;
 typedef PointF<double> DPoint;

-template<typename F>
-inline std::ostream& operator<<(std::ostream& os, const PointF<F>& pt)
+template <typename F>
+inline std::ostream &operator<<(std::ostream &os, const PointF<F> &pt)
 {
 	os << '(' << pt.mX << ',' << pt.mY << ',' << pt.mZ << ')';
-	return os; 
+	return os;
 }

-template<typename F>
-inline PointF<F> operator+(const PointF<F>& lhs, const PointF<F>& rhs)
+template <typename F>
+inline PointF<F> operator+(const PointF<F> &lhs, const PointF<F> &rhs)
 {
 	return PointF<F>(lhs.mX + rhs.mX, lhs.mY + rhs.mY, lhs.mZ + rhs.mZ);
 }

-template<typename F>
-inline PointF<F> operator-(const PointF<F>& lhs, const PointF<F>& rhs)
+template <typename F>
+inline PointF<F> operator-(const PointF<F> &lhs, const PointF<F> &rhs)
 {
 	return PointF<F>(lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ);
 }

-template<typename F>
-inline PointF<F> operator-(const PointF<F>& pt)
+template <typename F>
+inline PointF<F> operator-(const PointF<F> &pt)
 {
 	return PointF<F>(-pt.mX, -pt.mY, -pt.mZ);
 }

-template<typename F>
-inline PointF<F> operator*(const PointF<F>& pt, F f)
+template <typename F>
+inline PointF<F> operator*(const PointF<F> &pt, F f)
 {
 	return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
 }

-template<typename F>
-inline PointF<F> operator*(F f, const PointF<F>& pt)
+template <typename F>
+inline PointF<F> operator*(F f, const PointF<F> &pt)
 {
 	return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
 }

-template<typename F>
-inline PointF<F> operator/(const PointF<F>& pt, F f)
+template <typename F>
+inline PointF<F> operator/(const PointF<F> &pt, F f)
 {
 	return PointF<F>(pt.mX / f, pt.mY / f, pt.mZ / f);
 }
@@ -257,17 +274,16 @@ inline PointF<F> operator/(const PointF<F>& pt, F f)
 // --------------------------------------------------------------------
 // several standard 3d operations

-template<typename F>
-inline double DistanceSquared(const PointF<F>& a, const PointF<F>& b)
+template <typename F>
+inline double DistanceSquared(const PointF<F> &a, const PointF<F> &b)
 {
-	return
-		(a.mX - b.mX) * (a.mX - b.mX) +
-		(a.mY - b.mY) * (a.mY - b.mY) +
-		(a.mZ - b.mZ) * (a.mZ - b.mZ);
+	return (a.mX - b.mX) * (a.mX - b.mX) +
+	       (a.mY - b.mY) * (a.mY - b.mY) +
+	       (a.mZ - b.mZ) * (a.mZ - b.mZ);
 }

-template<typename F>
-inline double Distance(const PointF<F>& a, const PointF<F>& b)
+template <typename F>
+inline double Distance(const PointF<F> &a, const PointF<F> &b)
 {
 	return sqrt(
 		(a.mX - b.mX) * (a.mX - b.mX) +
@@ -275,44 +291,44 @@ inline double Distance(const PointF<F>& a, const PointF<F>& b)
 		(a.mZ - b.mZ) * (a.mZ - b.mZ));
 }

-template<typename F>
-inline F DotProduct(const PointF<F>& a, const PointF<F>& b)
+template <typename F>
+inline F DotProduct(const PointF<F> &a, const PointF<F> &b)
 {
 	return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ;
 }

-template<typename F>
-inline PointF<F> CrossProduct(const PointF<F>& a, const PointF<F>& b)
+template <typename F>
+inline PointF<F> CrossProduct(const PointF<F> &a, const PointF<F> &b)
 {
 	return PointF<F>(a.mY * b.mZ - b.mY * a.mZ,
-				  a.mZ * b.mX - b.mZ * a.mX,
-				  a.mX * b.mY - b.mX * a.mY);
+		a.mZ * b.mX - b.mZ * a.mX,
+		a.mX * b.mY - b.mX * a.mY);
 }

-template<typename F>
-double Angle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3)
+template <typename F>
+double Angle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3)
 {
 	PointF<F> v1 = p1 - p2;
 	PointF<F> v2 = p3 - p2;
-	
+
 	return std::acos(DotProduct(v1, v2) / (v1.length() * v2.length())) * 180 / kPI;
 }

-template<typename F>
-double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4)
+template <typename F>
+double DihedralAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
 {
-	PointF<F> v12 = p1 - p2;	// vector from p2 to p1
-	PointF<F> v43 = p4 - p3;	// vector from p3 to p4
-	
-	PointF<F> z = p2 - p3;		// vector from p3 to p2
-	
+	PointF<F> v12 = p1 - p2; // vector from p2 to p1
+	PointF<F> v43 = p4 - p3; // vector from p3 to p4
+
+	PointF<F> z = p2 - p3; // vector from p3 to p2
+
 	PointF<F> p = CrossProduct(z, v12);
 	PointF<F> x = CrossProduct(z, v43);
 	PointF<F> y = CrossProduct(z, x);
-	
+
 	double u = DotProduct(x, x);
 	double v = DotProduct(y, y);
-	
+
 	double result = 360;
 	if (u > 0 and v > 0)
 	{
@@ -321,33 +337,33 @@ double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>&
 		if (u != 0 or v != 0)
 			result = atan2(v, u) * 180 / kPI;
 	}
-	
+
 	return result;
 }

-template<typename F>
-double CosinusAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4)
+template <typename F>
+double CosinusAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
 {
 	PointF<F> v12 = p1 - p2;
 	PointF<F> v34 = p3 - p4;
-	
+
 	double result = 0;
-	
+
 	double x = DotProduct(v12, v12) * DotProduct(v34, v34);
 	if (x > 0)
 		result = DotProduct(v12, v34) / sqrt(x);
-	
+
 	return result;
 }

-template<typename F>
+template <typename F>
 auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<F> &p)
 {
-	auto line       = l2 - l1;
-    auto p_to_l1    = p - l1;
-    auto p_to_l2    = p - l2;
-    auto cross      = CrossProduct(p_to_l1, p_to_l2);
-    return cross.length() / line.length();
+	auto line = l2 - l1;
+	auto p_to_l1 = p - l1;
+	auto p_to_l2 = p - l2;
+	auto cross = CrossProduct(p_to_l1, p_to_l2);
+	return cross.length() / line.length();
 }

 // --------------------------------------------------------------------
@@ -355,7 +371,7 @@ auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<
 // a random direction with a distance randomly chosen from a normal
 // distribution with a stddev of offset.

-template<typename F>
+template <typename F>
 PointF<F> Nudge(PointF<F> p, F offset);

 // --------------------------------------------------------------------
@@ -363,66 +379,77 @@ PointF<F> Nudge(PointF<F> p, F offset);

 Quaternion Normalize(Quaternion q);

-std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q);
-Point Centroid(std::vector<Point>& Points);
-Point CenterPoints(std::vector<Point>& Points);
-Quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
-double RMSd(const std::vector<Point>& a, const std::vector<Point>& b);
+std::tuple<double, Point> QuaternionToAngleAxis(Quaternion q);
+Point Centroid(const std::vector<Point> &Points);
+Point CenterPoints(std::vector<Point> &Points);
+
+/// \brief Returns how the two sets of points \a a and \b b can be aligned
+///        
+/// \param a	The first set of points
+/// \param b    The second set of points
+/// \result     The quaternion which should be applied to the points in \a a to
+///             obtain the best superposition.
+Quaternion AlignPoints(const std::vector<Point> &a, const std::vector<Point> &b);
+
+/// \brief The RMSd for the points in \a a and \a b
+double RMSd(const std::vector<Point> &a, const std::vector<Point> &b);

 // --------------------------------------------------------------------
 // Helper class to generate evenly divided Points on a sphere
 // we use a fibonacci sphere to calculate even distribution of the dots

-template<int N>
+template <int N>
 class SphericalDots
 {
  public:
-	enum { P = 2 * N + 1 };
-	typedef typename std::array<Point,P>	array_type;
-	typedef typename array_type::const_iterator	iterator;
+	enum
+	{
+		P = 2 * N + 1
+	};
+	typedef typename std::array<Point, P> array_type;
+	typedef typename array_type::const_iterator iterator;

-	static SphericalDots& instance()
+	static SphericalDots &instance()
 	{
 		static SphericalDots sInstance;
 		return sInstance;
 	}
-	
-	size_t size() const							{ return mPoints.size(); }
-	const Point operator[](uint32_t inIx) const	{ return mPoints[inIx]; }
-	iterator begin() const						{ return mPoints.begin(); }
-	iterator end() const						{ return mPoints.end(); }

-	double weight() const						{ return mWeight; }
+	size_t size() const { return mPoints.size(); }
+	const Point operator[](uint32_t inIx) const { return mPoints[inIx]; }
+	iterator begin() const { return mPoints.begin(); }
+	iterator end() const { return mPoints.end(); }
+
+	double weight() const { return mWeight; }

 	SphericalDots()
 	{
-				
+
 		const double
 			kGoldenRatio = (1 + std::sqrt(5.0)) / 2;
-		
+
 		mWeight = (4 * kPI) / P;
-		
+
 		auto p = mPoints.begin();
-		
+
 		for (int32_t i = -N; i <= N; ++i)
 		{
 			double lat = std::asin((2.0 * i) / P);
 			double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio;
-			
+
 			p->mX = sin(lon) * cos(lat);
 			p->mY = cos(lon) * cos(lat);
-			p->mZ =            sin(lat);
+			p->mZ = sin(lat);

 			++p;
 		}
 	}

  private:
-
-	array_type				mPoints;
-	double					mWeight;
+	array_type mPoints;
+	double mWeight;
 };

 typedef SphericalDots<50> SphericalDots_50;

-}
+} // namespace mmcif
--- a/include/cif++/Secondary.hpp
+++ b/include/cif++/Secondary.hpp
@@ -137,6 +137,22 @@ class DSSP
 	  public:
 		friend class iterator;

+		ResidueInfo()
+			: mImpl(nullptr)
+		{
+		}
+
+		ResidueInfo(const ResidueInfo &rhs)
+			: mImpl(rhs.mImpl)
+		{
+		}
+
+		ResidueInfo& operator=(const ResidueInfo &rhs)
+		{
+			mImpl = rhs.mImpl;
+			return *this;
+		}
+
 		explicit operator bool() const		{ return not empty(); }
 		bool empty() const					{ return mImpl == nullptr; }

@@ -168,6 +184,15 @@ class DSSP
 		std::tuple<ResidueInfo,double> acceptor(int i) const;
 		std::tuple<ResidueInfo,double> donor(int i) const;

+		/// \brief Simple compare equals
+		bool operator==(const ResidueInfo &rhs) const
+		{
+			return mImpl == rhs.mImpl;
+		}
+
+		/// \brief Returns \result true if there is a bond between two residues
+		friend bool TestBond(ResidueInfo const &a, ResidueInfo const &b);
+
 	  private:
 		ResidueInfo(Res* res) : mImpl(res) {}

@@ -177,7 +202,7 @@ class DSSP
 	class iterator
 	{
 	  public:
-		using iterator_category = std::input_iterator_tag;
+		using iterator_category = std::bidirectional_iterator_tag;
 		using value_type = ResidueInfo;
 		using difference_type = std::ptrdiff_t;
 		using pointer = value_type*;
@@ -198,6 +223,14 @@ class DSSP
 			return tmp;
 		}

+		iterator& operator--();
+		iterator operator--(int)
+		{
+			auto tmp(*this);
+			this->operator--();
+			return tmp;
+		}
+
 		bool operator==(const iterator& rhs) const		{ return mCurrent.mImpl == rhs.mCurrent.mImpl; }
 		bool operator!=(const iterator& rhs) const		{ return mCurrent.mImpl != rhs.mCurrent.mImpl; }

--- a/include/cif++/Structure.hpp
+++ b/include/cif++/Structure.hpp
@@ -60,30 +60,101 @@ class File;

 class Atom
 {
+  private:
+
+	struct AtomImpl : public std::enable_shared_from_this<AtomImpl>
+	{
+		AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row);
+
+		// constructor for a symmetry copy of an atom
+		AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op);
+
+		AtomImpl(const AtomImpl &i) = default;
+
+		void prefetch();
+
+		int compare(const AtomImpl &b) const;
+
+		bool getAnisoU(float anisou[6]) const;
+
+		void moveTo(const Point &p);
+
+		const Compound &comp() const;
+
+		const std::string get_property(const std::string_view name) const;
+		void set_property(const std::string_view name, const std::string &value);
+
+		const cif::Datablock &mDb;
+		std::string mID;
+		AtomType mType;
+
+		std::string mAtomID;
+		std::string mCompID;
+		std::string mAsymID;
+		int mSeqID;
+		std::string mAltID;
+		std::string mAuthSeqID;
+
+		Point mLocation;
+		int mRefcount;
+		cif::Row mRow;
+
+		mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
+
+		mutable const Compound *mCompound = nullptr;
+
+		bool mSymmetryCopy = false;
+		bool mClone = false;
+
+		std::string mSymmetryOperator = "1_555";
+	};
+
  public:
-	Atom();
-	Atom(struct AtomImpl *impl);
-	Atom(const Atom &rhs);
+
+	Atom() {}
+
+	Atom(std::shared_ptr<AtomImpl> impl)
+		: mImpl(impl) {}
+
+	Atom(const Atom &rhs)
+		: mImpl(rhs.mImpl) {}

 	Atom(cif::Datablock &db, cif::Row &row);

 	// a special constructor to create symmetry copies
 	Atom(const Atom &rhs, const Point &symmmetry_location, const std::string &symmetry_operation);

-	~Atom();
-
-	explicit operator bool() const { return mImpl_ != nullptr; }
+	explicit operator bool() const { return (bool)mImpl; }

 	// return a copy of this atom, with data copied instead of referenced
-	Atom clone() const;
+	Atom clone() const
+	{
+		auto copy = std::make_shared<AtomImpl>(*mImpl);
+		copy->mClone = true;
+		return Atom(copy);
+	}

-	Atom &operator=(const Atom &rhs);
+	Atom &operator=(const Atom &rhs) = default;

-	const std::string &id() const;
-	AtomType type() const;
+	template <typename T>
+	T get_property(const std::string_view name) const;

-	Point location() const;
-	void location(Point p);
+	void set_property(const std::string_view name, const std::string &value)
+	{
+		mImpl->set_property(name, value);
+	}
+
+	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
+	void property(const std::string_view name, const T &value)
+	{
+		set_property(name, std::to_string(value));
+	}
+
+	const std::string &id() const			{ return mImpl->mID; }
+	AtomType type() const					{ return mImpl->mType; }
+
+	Point location() const					{ return mImpl->mLocation; }
+	void location(Point p)					{ mImpl->moveTo(p); }

 	/// \brief Translate the position of this atom by \a t
 	void translate(Point t);
@@ -91,47 +162,40 @@ class Atom
 	/// \brief Rotate the position of this atom by \a q
 	void rotate(Quaternion q);

+	/// \brief Translate and rotate the position of this atom by \a t and \a q
+	void translateAndRotate(Point t, Quaternion q);
+
+	/// \brief Translate, rotate and translate again the coordinates this atom by \a t1 , \a q and \a t2
+	void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
+
 	// for direct access to underlying data, be careful!
-	const cif::Row getRow() const;
+	const cif::Row getRow() const			{ return mImpl->mRow; }
 	const cif::Row getRowAniso() const;

-	// Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
-	bool isSymmetryCopy() const;
-	std::string symmetry() const;
-	// const clipper::RTop_orth& symop() const;
+	bool isSymmetryCopy() const				{ return mImpl->mSymmetryCopy; }
+	std::string symmetry() const			{ return mImpl->mSymmetryOperator; }

-	const Compound &comp() const;
-	bool isWater() const;
+	const Compound &comp() const			{ return mImpl->comp(); }
+	bool isWater() const					{ return mImpl->mCompID == "HOH" or mImpl->mCompID == "H2O" or mImpl->mCompID == "WAT"; }
 	int charge() const;

 	float uIso() const;
-	bool getAnisoU(float anisou[6]) const;
+	bool getAnisoU(float anisou[6]) const	{ return mImpl->getAnisoU(anisou); }
 	float occupancy() const;

-	template <typename T>
-	T property(const std::string &name) const;
-
-	void property(const std::string &name, const std::string &value);
-
-	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
-	void property(const std::string &name, const T &value)
-	{
-		property(name, std::to_string(value));
-	}
-
 	// specifications
-	std::string labelAtomID() const;
-	std::string labelCompID() const;
-	std::string labelAsymID() const;
+	const std::string& labelAtomID() const		{ return mImpl->mAtomID; }
+	const std::string& labelCompID() const		{ return mImpl->mCompID; }
+	const std::string& labelAsymID() const		{ return mImpl->mAsymID; }
 	std::string labelEntityID() const;
-	int labelSeqID() const;
-	std::string labelAltID() const;
-	bool isAlternate() const;
+	int labelSeqID() const						{ return mImpl->mSeqID; }
+	const std::string& labelAltID() const		{ return mImpl->mAltID; }
+	bool isAlternate() const					{ return not mImpl->mAltID.empty(); }

 	std::string authAtomID() const;
 	std::string authCompID() const;
 	std::string authAsymID() const;
-	std::string authSeqID() const;
+	const std::string& authSeqID() const		{ return mImpl->mAuthSeqID; }
 	std::string pdbxAuthInsCode() const;
 	std::string pdbxAuthAltID() const;

@@ -140,13 +204,6 @@ class Atom

 	bool operator==(const Atom &rhs) const;

-	// // get clipper format Atom
-	// clipper::Atom toClipper() const;
-
-	// Radius calculation based on integrating the density until perc of electrons is found
-	void calculateRadius(float resHigh, float resLow, float perc);
-	float radius() const;
-
 	// access data in compound for this atom

 	// convenience routine
@@ -158,10 +215,10 @@ class Atom

 	void swap(Atom &b)
 	{
-		std::swap(mImpl_, b.mImpl_);
+		std::swap(mImpl, b.mImpl);
 	}

-	int compare(const Atom &b) const;
+	int compare(const Atom &b) const		{ return mImpl->compare(*b.mImpl); }

 	bool operator<(const Atom &rhs) const
 	{
@@ -172,14 +229,31 @@ class Atom

  private:
 	friend class Structure;
+
 	void setID(int id);

-	AtomImpl *impl();
-	const AtomImpl *impl() const;
-
-	struct AtomImpl *mImpl_;
+	std::shared_ptr<AtomImpl> mImpl;
 };

+template <>
+inline std::string Atom::get_property<std::string>(const std::string_view name) const
+{
+	return mImpl->get_property(name);
+}
+
+template <>
+inline int Atom::get_property<int>(const std::string_view name) const
+{
+	auto v = mImpl->get_property(name);
+	return v.empty() ? 0 : stoi(v);
+}
+
+template <>
+inline float Atom::get_property<float>(const std::string_view name) const
+{
+	return stof(mImpl->get_property(name));
+}
+
 inline void swap(mmcif::Atom &a, mmcif::Atom &b)
 {
 	a.swap(b);
@@ -202,19 +276,16 @@ typedef std::vector<Atom> AtomView;
 class Residue
 {
  public:
-	// constructors should be private, but that's not possible for now (needed in emplace)
-
-	// constructor for waters
+	// constructor
 	Residue(const Structure &structure, const std::string &compoundID,
-		const std::string &asymID, const std::string &authSeqID);
-
-	// constructor for a residue without a sequence number
-	Residue(const Structure &structure, const std::string &compoundID,
-		const std::string &asymID);
-
-	// constructor for a residue with a sequence number
-	Residue(const Structure &structure, const std::string &compoundID,
-		const std::string &asymID, int seqID, const std::string &authSeqID);
+		const std::string &asymID, int seqID = 0, const std::string &authSeqID = {})
+		: mStructure(&structure)
+		, mCompoundID(compoundID)
+		, mAsymID(asymID)
+		, mSeqID(seqID)
+		, mAuthSeqID(authSeqID)
+	{
+	}

 	Residue(const Residue &rhs) = delete;
 	Residue &operator=(const Residue &rhs) = delete;
@@ -227,6 +298,11 @@ class Residue
 	const Compound &compound() const;
 	const AtomView &atoms() const;

+	void addAtom(const Atom &atom)
+	{
+		mAtoms.push_back(atom);
+	}
+
 	/// \brief Unique atoms returns only the atoms without alternates and the first of each alternate atom id.
 	AtomView unique_atoms() const;

@@ -236,6 +312,8 @@ class Residue
 	Atom atomByID(const std::string &atomID) const;

 	const std::string &compoundID() const { return mCompoundID; }
+	void setCompoundID(const std::string &id) { mCompoundID = id; }
+
 	const std::string &asymID() const { return mAsymID; }
 	int seqID() const { return mSeqID; }
 	std::string entityID() const;
@@ -275,6 +353,8 @@ class Residue

 	friend std::ostream &operator<<(std::ostream &os, const Residue &res);

+	friend Structure;
+
  protected:
 	Residue() {}

@@ -404,7 +484,7 @@ class File : public std::enable_shared_from_this<File>
 	File(const File &) = delete;
 	File &operator=(const File &) = delete;

-	cif::Datablock& createDatablock(const std::string &name);
+	cif::Datablock& createDatablock(const std::string_view name);

 	void load(const std::filesystem::path &path);
 	void save(const std::filesystem::path &path);
@@ -461,9 +541,24 @@ class Structure
 	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 residue, if \a seqID is zero, the non-polymers are searched
 	const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0) const;

+	/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
+	Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0);
+
+	/// \brief Get a the single residue for an asym with id \a asymID
+	const Residue &getResidue(const std::string &asymID) const;
+
+	/// \brief Get a the single residue for an asym with id \a asymID
+	Residue &getResidue(const std::string &asymID);
+
 	// map between auth and label locations

 	std::tuple<std::string, int, std::string> MapAuthToLabel(const std::string &asymID,
@@ -488,7 +583,7 @@ class Structure
 	void removeAtom(Atom &a);
 	void swapAtoms(Atom &a1, Atom &a2); // swap the labels for these atoms
 	void moveAtom(Atom &a, Point p);    // move atom to a new location
-	void changeResidue(const Residue &res, const std::string &newCompound,
+	void changeResidue(Residue &res, const std::string &newCompound,
 		const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);

 	/// \brief Create a new non-polymer entity, returns new ID
@@ -514,20 +609,27 @@ class Structure
 	/// \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 translateAndRotate(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 translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
+
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
 	const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }

 	void cleanupEmptyCategories();

+	/// \brief Direct access to underlying data
+	cif::Category &category(std::string_view name) const;
+	cif::Datablock &datablock() const;
+
  private:
 	friend Polymer;
 	friend Residue;
 	// friend residue_view;
 	// friend residue_iterator;

-	cif::Category &category(const char *name) const;
-	cif::Datablock &datablock() const;
-
 	std::string insertCompound(const std::string &compoundID, bool isEntity);

 	void loadData();
--- a/include/cif++/Symmetry.hpp
+++ b/include/cif++/Symmetry.hpp
@@ -37,6 +37,11 @@ namespace mmcif

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

+enum class SpacegroupName
+{
+	full, xHM, Hall
+};
+
 struct Spacegroup
 {
 	const char* name;
@@ -133,6 +138,7 @@ CIFPP_EXPORT extern const std::size_t kSymopNrTableSize;

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

-int GetSpacegroupNumber(std::string spacegroup);	// alternative for clipper's parsing code
+int GetSpacegroupNumber(std::string spacegroup);						// alternative for clipper's parsing code, using SpacegroupName::full
+int GetSpacegroupNumber(std::string spacegroup, SpacegroupName type);	// alternative for clipper's parsing code

 }
--- a/src/AtomType.cpp
+++ b/src/AtomType.cpp
--- a/src/BondMap.cpp
+++ b/src/BondMap.cpp
@@ -1,17 +1,17 @@
 /*-
 * SPDX-License-Identifier: BSD-2-Clause
- * 
+ *
 * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- * 
+ *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
- * 
+ *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
- * 
+ *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -24,14 +24,14 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

-#include <fstream>
 #include <algorithm>
+#include <fstream>
 #include <mutex>

-#include "cif++/Cif++.hpp"
-#include "cif++/Compound.hpp"
-#include "cif++/CifUtils.hpp"
 #include "cif++/BondMap.hpp"
+#include "cif++/Cif++.hpp"
+#include "cif++/CifUtils.hpp"
+#include "cif++/Compound.hpp"

 namespace mmcif
 {
@@ -39,223 +39,79 @@ namespace mmcif
 namespace
 {

-union IDType
-{
-	IDType() : id_n(0){}
-	IDType(const IDType& rhs) : id_n(rhs.id_n) {}
-	IDType(const std::string& s)
-		: IDType()
+	union IDType
 	{
-		assert(s.length() <= 4);
-		if (s.length() > 4)
-			throw BondMapException("Atom ID '" + s + "' is too long");
-		std::copy(s.begin(), s.end(), id_s);
-	}
+		IDType()
+			: id_n(0)
+		{
+		}
+		IDType(const IDType &rhs)
+			: id_n(rhs.id_n)
+		{
+		}
+		IDType(const std::string &s)
+			: IDType()
+		{
+			assert(s.length() <= 4);
+			if (s.length() > 4)
+				throw BondMapException("Atom ID '" + s + "' is too long");
+			std::copy(s.begin(), s.end(), id_s);
+		}

-	IDType& operator=(const IDType& rhs)
-	{
-		id_n = rhs.id_n;
-		return *this;
-	}
+		IDType &operator=(const IDType &rhs)
+		{
+			id_n = rhs.id_n;
+			return *this;
+		}

-	IDType& operator=(const std::string& s)
-	{
-		id_n = 0;
-		assert(s.length() <= 4);
-		if (s.length() > 4)
-			throw BondMapException("Atom ID '" + s + "' is too long");
-		std::copy(s.begin(), s.end(), id_s);
-		return *this;
-	}
+		IDType &operator=(const std::string &s)
+		{
+			id_n = 0;
+			assert(s.length() <= 4);
+			if (s.length() > 4)
+				throw BondMapException("Atom ID '" + s + "' is too long");
+			std::copy(s.begin(), s.end(), id_s);
+			return *this;
+		}

-	bool operator<(const IDType& rhs) const
-	{
-		return id_n < rhs.id_n;
-	}
+		bool operator<(const IDType &rhs) const
+		{
+			return id_n < rhs.id_n;
+		}

-	bool operator<=(const IDType& rhs) const
-	{
-		return id_n <= rhs.id_n;
-	}
+		bool operator<=(const IDType &rhs) const
+		{
+			return id_n <= rhs.id_n;
+		}

-	bool operator==(const IDType& rhs) const
-	{
-		return id_n == rhs.id_n;
-	}
+		bool operator==(const IDType &rhs) const
+		{
+			return id_n == rhs.id_n;
+		}

-	bool operator!=(const IDType& rhs) const
-	{
-		return id_n != rhs.id_n;
-	}
+		bool operator!=(const IDType &rhs) const
+		{
+			return id_n != rhs.id_n;
+		}

-	char		id_s[4];
-	uint32_t	id_n;
-};
+		char id_s[4];
+		uint32_t id_n;
+	};

-static_assert(sizeof(IDType) == 4, "atom_id_type should be 4 bytes");	
-}
-
-// // --------------------------------------------------------------------
-
-// void createBondInfoFile(const fs::path& components, const fs::path& infofile)
-// {
-// 	std::ofstream outfile(infofile.string() + ".tmp", std::ios::binary);
-// 	if (not outfile.is_open())
-// 		throw BondMapException("Could not create bond info file " + infofile.string() + ".tmp");
-
-// 	cif::File infile(components);
-
-// 	std::set<atom_id_type> atomIDs;
-// 	std::vector<atom_id_type> compoundIDs;
-
-// 	for (auto& db: infile)
-// 	{
-// 		auto chem_comp_bond = db.get("chem_comp_bond");
-// 		if (not chem_comp_bond)
-// 		{
-// 			if (cif::VERBOSE > 1)
-// 				std::cerr << "Missing chem_comp_bond category in data block " << db.getName() << std::endl;
-// 			continue;
-// 		}
-
-// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
-// 		{
-// 			atomIDs.insert(atom_id_1);
-// 			atomIDs.insert(atom_id_2);
-// 		}
-
-// 		compoundIDs.push_back({ db.getName() });
-// 	}
-
-// 	if (cif::VERBOSE)
-// 		std::cout << "Number of unique atom names is " << atomIDs.size() << std::endl
-// 				  << "Number of unique residue names is " << compoundIDs.size() << std::endl;
-
-// 	CompoundBondInfoFileHeader header = {};
-// 	header.indexEntries = compoundIDs.size();
-// 	header.atomEntries = atomIDs.size();
-
-// 	outfile << header;
-	
-// 	for (auto atomID: atomIDs)
-// 		outfile << atomID;
-
-// 	auto dataOffset = outfile.tellp();
-
-// 	std::vector<CompoundBondInfo> entries;
-// 	entries.reserve(compoundIDs.size());
-
-// 	std::map<atom_id_type, uint16_t> atomIDMap;
-// 	for (auto& atomID: atomIDs)
-// 		atomIDMap[atomID] = atomIDMap.size();
-
-// 	for (auto& db: infile)
-// 	{
-// 		auto chem_comp_bond = db.get("chem_comp_bond");
-// 		if (not chem_comp_bond)
-// 			continue;
-
-// 		std::set<uint16_t> bondedAtoms;
-
-// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
-// 		{
-// 			bondedAtoms.insert(atomIDMap[atom_id_1]);
-// 			bondedAtoms.insert(atomIDMap[atom_id_2]);
-// 		}
-
-// 		std::map<uint16_t, int32_t> bondedAtomMap;
-// 		for (auto id: bondedAtoms)
-// 			bondedAtomMap[id] = static_cast<int32_t>(bondedAtomMap.size());
-		
-// 		CompoundBondInfo info = {
-// 			db.getName(),
-// 			static_cast<uint32_t>(bondedAtomMap.size()),
-// 			outfile.tellp() - dataOffset
-// 		};
-
-// 		entries.push_back(info);
-
-// 		// An now first write the array of atom ID's in this compound
-// 		for (uint16_t id: bondedAtoms)
-// 			write(outfile, id);
-
-// 		// And then the symmetric matrix with bonds
-// 		size_t N = bondedAtoms.size();
-// 		size_t M = (N * (N - 1)) / 2;
-
-// 		size_t K = M / 8;
-// 		if (M % 8)
-// 			K += 1;
-		
-// 		std::vector<uint8_t> m(K);
-
-// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
-// 		{
-// 			auto a = bondedAtomMap[atomIDMap[atom_id_1]];
-// 			auto b = bondedAtomMap[atomIDMap[atom_id_2]];
-
-// 			assert(a != b);
-// 			assert((int)b < (int)N);
-
-// 			if (a > b)
-// 				std::swap(a, b);
-			
-// 			size_t ix = ((b - 1) * b) / 2 + a;
-// 			assert(ix < M);
-
-// 			auto Bix = ix / 8;
-// 			auto bix = ix % 8;
-
-// 			m[Bix] |= 1 << bix;
-// 		}
-
-// 		outfile.write(reinterpret_cast<char*>(m.data()), m.size());
-// 	}
-
-// 	header.dataSize = outfile.tellp() - dataOffset;
-
-// 	std::sort(entries.begin(), entries.end(), [](CompoundBondInfo& a, CompoundBondInfo& b)
-// 	{
-// 		return a.id < b.id;
-// 	});
-
-// 	for (auto& info: entries)
-// 		outfile << info;
-
-// 	outfile.seekp(0);
-// 	outfile << header;
-
-// 	// compress
-// 	outfile.close();
-
-// 	std::ifstream in(infofile.string() + ".tmp", std::ios::binary);
-// 	std::ofstream out(infofile, std::ios::binary);
-
-// 	{
-// 		io::filtering_stream<io::output> os;
-// 		os.push(io::gzip_compressor());
-// 		os.push(out);
-// 		io::copy(in, os);
-// 	}
-
-// 	in.close();
-// 	out.close();
-
-// 	fs::remove(infofile.string() + ".tmp");
-// }
+	static_assert(sizeof(IDType) == 4, "atom_id_type should be 4 bytes");
+} // namespace

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

 struct CompoundBondInfo
 {
-	IDType	mID;
-	std::set<std::tuple<uint32_t,uint32_t>> mBonded;
+	IDType mID;
+	std::set<std::tuple<uint32_t, uint32_t>> mBonded;

 	bool bonded(uint32_t a1, uint32_t a2) const
 	{
-		return mBonded.count({ a1, a2 }) > 0;
+		return mBonded.count({a1, a2}) > 0;
 	}
-
 };

 // --------------------------------------------------------------------
@@ -263,20 +119,18 @@ struct CompoundBondInfo
 class CompoundBondMap
 {
  public:
-
 	static CompoundBondMap &instance()
 	{
 		static std::unique_ptr<CompoundBondMap> s_instance(new CompoundBondMap);
 		return *s_instance;
 	}

-	bool bonded(const std::string& compoundID, const std::string& atomID1, const std::string& atomID2);
+	bool bonded(const std::string &compoundID, const std::string &atomID1, const std::string &atomID2);

  private:
-
 	CompoundBondMap() {}

-	uint32_t getAtomID(const std::string& atomID)
+	uint32_t getAtomID(const std::string &atomID)
 	{
 		IDType id(atomID);

@@ -290,16 +144,16 @@ class CompoundBondMap
 		}
 		else
 			result = i->second;
-		
+
 		return result;
 	}

-	std::map<IDType,uint32_t> mAtomIDIndex;
+	std::map<IDType, uint32_t> mAtomIDIndex;
 	std::vector<CompoundBondInfo> mCompounds;
 	std::mutex mMutex;
 };

-bool CompoundBondMap::bonded(const std::string &compoundID, const std::string& atomID1, const std::string& atomID2)
+bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &atomID1, const std::string &atomID2)
 {
 	std::lock_guard lock(mMutex);

@@ -310,32 +164,35 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string& a
 	uint32_t a2 = getAtomID(atomID2);
 	if (a1 > a2)
 		std::swap(a1, a2);
-	
-	for (auto &bi: mCompounds)
+
+	for (auto &bi : mCompounds)
 	{
 		if (bi.mID != id)
 			continue;
-		
+
 		return bi.bonded(a1, a2);
 	}

 	bool result = false;

 	// not found in our cache, calculate
-	CompoundBondInfo bondInfo{ id };
+	CompoundBondInfo bondInfo{id};

 	auto compound = mmcif::CompoundFactory::instance().create(compoundID);
 	if (not compound)
-		std::cerr << "Missing compound bond info for " << compoundID << std::endl;
+	{
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Missing compound bond info for " << compoundID << std::endl;
+	}
 	else
 	{
-		for (auto &atom: compound->bonds())
+		for (auto &atom : compound->bonds())
 		{
 			uint32_t ca1 = getAtomID(atom.atomID[0]);
 			uint32_t ca2 = getAtomID(atom.atomID[1]);
 			if (ca1 > ca2)
 				std::swap(ca1, ca2);
-			
+
 			bondInfo.mBonded.insert({ca1, ca2});
 			result = result or (a1 == ca1 and a2 == ca2);
 		}
@@ -348,27 +205,27 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string& a

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

-BondMap::BondMap(const Structure& p)
+BondMap::BondMap(const Structure &p)
 {
-	auto& compoundBondInfo = CompoundBondMap::instance();
+	auto &compoundBondInfo = CompoundBondMap::instance();

 	auto atoms = p.atoms();
 	dim = uint32_t(atoms.size());

-//	bond = std::vector<bool>(dim * (dim - 1), false);
+	//	bond = std::vector<bool>(dim * (dim - 1), false);

-	for (auto& atom: atoms)
+	for (auto &atom : atoms)
 		index[atom.id()] = uint32_t(index.size());
-	
-	auto bindAtoms = [this](const std::string& a, const std::string& b)
+
+	auto bindAtoms = [this](const std::string &a, const std::string &b)
 	{
 		uint32_t ixa = index[a];
 		uint32_t ixb = index[b];
-		
+
 		bond.insert(key(ixa, ixb));
 	};

-	auto linkAtoms = [this,&bindAtoms](const std::string& a, const std::string& b)
+	auto linkAtoms = [this, &bindAtoms](const std::string &a, const std::string &b)
 	{
 		bindAtoms(a, b);

@@ -376,20 +233,20 @@ BondMap::BondMap(const Structure& p)
 		link[b].insert(a);
 	};

-	cif::Datablock& db = p.getFile().data();
+	cif::Datablock &db = p.getFile().data();

 	// collect all compounds first
 	std::set<std::string> compounds;
-	for (auto c: db["chem_comp"])
+	for (auto c : db["chem_comp"])
 		compounds.insert(c["id"].as<std::string>());
-	
+
 	// make sure we also have all residues in the polyseq
-	for (auto m: db["entity_poly_seq"])
+	for (auto m : db["entity_poly_seq"])
 	{
 		std::string c = m["mon_id"].as<std::string>();
 		if (compounds.count(c))
 			continue;
-		
+
 		if (cif::VERBOSE > 1)
 			std::cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << std::endl;
 		compounds.insert(c);
@@ -398,68 +255,68 @@ BondMap::BondMap(const Structure& p)
 	cif::Progress progress(compounds.size(), "Creating bond map");

 	// some helper indices to speed things up a bit
-	std::map<std::tuple<std::string,int,std::string>,std::string> atomMapByAsymSeqAndAtom;
-	for (auto& a: p.atoms())
+	std::map<std::tuple<std::string, int, std::string>, std::string> atomMapByAsymSeqAndAtom;
+	for (auto &a : p.atoms())
 	{
 		auto key = make_tuple(a.labelAsymID(), a.labelSeqID(), a.labelAtomID());
 		atomMapByAsymSeqAndAtom[key] = a.id();
 	}
-	
+
 	// first link all residues in a polyseq
-	
+
 	std::string lastAsymID;
 	int lastSeqID = 0;
-	for (auto r: db["pdbx_poly_seq_scheme"])
+	for (auto r : db["pdbx_poly_seq_scheme"])
 	{
 		std::string asymID;
 		int seqID;

 		cif::tie(asymID, seqID) = r.get("asym_id", "seq_id");

-		if (asymID != lastAsymID)		// first in a new sequece
+		if (asymID != lastAsymID) // first in a new sequece
 		{
 			lastAsymID = asymID;
 			lastSeqID = seqID;
 			continue;
 		}
-		
+
 		auto c = atomMapByAsymSeqAndAtom[make_tuple(asymID, lastSeqID, "C")];
 		auto n = atomMapByAsymSeqAndAtom[make_tuple(asymID, seqID, "N")];

-		if (not (c.empty() or n.empty()))
+		if (not(c.empty() or n.empty()))
 			bindAtoms(c, n);
-		
+
 		lastSeqID = seqID;
 	}

-	for (auto l: db["struct_conn"])
+	for (auto l : db["struct_conn"])
 	{
 		std::string asym1, asym2, atomId1, atomId2;
 		int seqId1 = 0, seqId2 = 0;
 		cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2) =
 			l.get("ptnr1_label_asym_id", "ptnr2_label_asym_id",
-				  "ptnr1_label_atom_id", "ptnr2_label_atom_id",
-				  "ptnr1_label_seq_id", "ptnr2_label_seq_id");
+				"ptnr1_label_atom_id", "ptnr2_label_atom_id",
+				"ptnr1_label_seq_id", "ptnr2_label_seq_id");

 		std::string a = atomMapByAsymSeqAndAtom[make_tuple(asym1, seqId1, atomId1)];
 		std::string b = atomMapByAsymSeqAndAtom[make_tuple(asym2, seqId2, atomId2)];
-			
-		if (not (a.empty() or b.empty()))
+
+		if (not(a.empty() or b.empty()))
 			linkAtoms(a, b);
 	}

 	// then link all atoms in the compounds
-	
-	for (auto c: compounds)
+
+	for (auto c : compounds)
 	{
 		if (c == "HOH" or c == "H2O" or c == "WAT")
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "skipping water in bond map calculation" << std::endl;
 			continue;
 		}

-		auto bonded = [c, &compoundBondInfo](const Atom& a, const Atom& b)
+		auto bonded = [c, &compoundBondInfo](const Atom &a, const Atom &b)
 		{
 			auto label_a = a.labelAtomID();
 			auto label_b = b.labelAtomID();
@@ -468,16 +325,17 @@ BondMap::BondMap(const Structure& p)
 		};

 		// loop over poly_seq_scheme
-		for (auto r: db["pdbx_poly_seq_scheme"].find(cif::Key("mon_id") == c))
+		for (auto r : db["pdbx_poly_seq_scheme"].find(cif::Key("mon_id") == c))
 		{
 			std::string asymID;
 			int seqID;
 			cif::tie(asymID, seqID) = r.get("asym_id", "seq_id");
-			
+
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[&](auto& a) { return a.labelAsymID() == asymID and a.labelSeqID() == seqID; });
-			
+				[&](auto &a)
+				{ return a.labelAsymID() == asymID and a.labelSeqID() == seqID; });
+
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -489,15 +347,16 @@ BondMap::BondMap(const Structure& p)
 		}

 		// loop over pdbx_nonpoly_scheme
-		for (auto r: db["pdbx_nonpoly_scheme"].find(cif::Key("mon_id") == c))
+		for (auto r : db["pdbx_nonpoly_scheme"].find(cif::Key("mon_id") == c))
 		{
 			std::string asymID;
 			cif::tie(asymID) = r.get("asym_id");
-			
+
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[&](auto& a) { return a.labelAsymID() == asymID; });
-			
+				[&](auto &a)
+				{ return a.labelAsymID() == asymID; });
+
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -506,7 +365,7 @@ BondMap::BondMap(const Structure& p)
 					{
 						uint32_t ixa = index[rAtoms[i].id()];
 						uint32_t ixb = index[rAtoms[j].id()];
-						
+
 						bond.insert(key(ixa, ixb));
 					}
 				}
@@ -514,15 +373,16 @@ BondMap::BondMap(const Structure& p)
 		}

 		// loop over pdbx_branch_scheme
-		for (auto r: db["pdbx_branch_scheme"].find(cif::Key("mon_id") == c))
+		for (auto r : db["pdbx_branch_scheme"].find(cif::Key("mon_id") == c))
 		{
 			std::string asymID;
 			cif::tie(asymID) = r.get("asym_id");
-			
+
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[&](auto& a) { return a.labelAsymID() == asymID; });
-			
+				[&](auto &a)
+				{ return a.labelAsymID() == asymID; });
+
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -531,31 +391,31 @@ BondMap::BondMap(const Structure& p)
 					{
 						uint32_t ixa = index[rAtoms[i].id()];
 						uint32_t ixb = index[rAtoms[j].id()];
-						
+
 						bond.insert(key(ixa, ixb));
 					}
 				}
 			}
 		}
 	}
-	
+
 	// start by creating an index for single bonds
-	
-	std::multimap<uint32_t,uint32_t> b1_2;
-	for (auto& bk: bond)
+
+	std::multimap<uint32_t, uint32_t> b1_2;
+	for (auto &bk : bond)
 	{
 		uint32_t a, b;
 		std::tie(a, b) = dekey(bk);
-		
-		b1_2.insert({ a, b });
-		b1_2.insert({ b, a });
+
+		b1_2.insert({a, b});
+		b1_2.insert({b, a});
 	}
-	
-	std::multimap<uint32_t,uint32_t> b1_3;
+
+	std::multimap<uint32_t, uint32_t> b1_3;
 	for (uint32_t i = 0; i < dim; ++i)
 	{
 		auto a = b1_2.equal_range(i);
-		
+
 		std::vector<uint32_t> s;
 		for (auto j = a.first; j != a.second; ++j)
 			s.push_back(j->second);
@@ -566,12 +426,12 @@ BondMap::BondMap(const Structure& p)
 			{
 				uint32_t x = s[si1];
 				uint32_t y = s[si2];
-				
+
 				if (isBonded(x, y))
 					continue;
-				
-				b1_3.insert({ x, y });
-				b1_3.insert({ y, x });
+
+				b1_3.insert({x, y});
+				b1_3.insert({y, x});
 			}
 		}
 	}
@@ -580,48 +440,48 @@ BondMap::BondMap(const Structure& p)
 	{
 		auto a1 = b1_2.equal_range(i);
 		auto a2 = b1_3.equal_range(i);
-		
+
 		for (auto ai1 = a1.first; ai1 != a1.second; ++ai1)
 		{
 			for (auto ai2 = a2.first; ai2 != a2.second; ++ai2)
 			{
 				uint32_t b1 = ai1->second;
 				uint32_t b2 = ai2->second;
-				
+
 				if (isBonded(b1, b2))
 					continue;
-				
+
 				bond_1_4.insert(key(b1, b2));
 			}
 		}
 	}
 }

-std::vector<std::string> BondMap::linked(const Atom& a) const
+std::vector<std::string> BondMap::linked(const Atom &a) const
 {
 	auto i = link.find(a.id());
-	
+
 	std::vector<std::string> result;
-	
+
 	if (i != link.end())
 		result = std::vector<std::string>(i->second.begin(), i->second.end());

 	return result;
 }

-std::vector<std::string> BondMap::atomIDsForCompound(const std::string& compoundID)
+std::vector<std::string> BondMap::atomIDsForCompound(const std::string &compoundID)
 {
 	std::vector<std::string> result;

-	auto* compound = mmcif::CompoundFactory::instance().create(compoundID);
+	auto *compound = mmcif::CompoundFactory::instance().create(compoundID);

 	if (compound == nullptr)
 		throw BondMapException("Missing bond information for compound " + compoundID);

-	for (auto& compAtom: compound->atoms())
+	for (auto &compAtom : compound->atoms())
 		result.push_back(compAtom.id);

 	return result;
 }

-}
+} // namespace mmcif
--- a/src/Cif++.cpp
+++ b/src/Cif++.cpp
--- a/src/Cif2PDB.cpp
+++ b/src/Cif2PDB.cpp
@@ -249,7 +249,7 @@ size_t WriteContinuedLine(std::ostream& pdbFile, std::string header, int& count,

 	for (auto& line: lines)
 	{
-		ba::to_upper(line);
+		// ba::to_upper(line);

 		pdbFile << header;
 		
@@ -753,7 +753,8 @@ class Ff : public FBase
            }
            catch (const std::exception& ex)
            {
-                std::cerr << "Failed to write '" << s << "' as a double, this indicates an error in the code for writing PDB files" << std::endl;
+				if (cif::VERBOSE >= 0)
+	                std::cerr << "Failed to write '" << s << "' as a double, this indicates an error in the code for writing PDB files" << std::endl;
                os << s;
            }
        }
@@ -811,7 +812,7 @@ typedef RM<3> RM3;
 template<int N>
 std::ostream& operator<<(std::ostream& os, RM<N>&& rm)
 {
-	os << "REMARK " << std::setw(3) << std::right << N << " " << rm.mDesc << (rm.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(abs(rm.mWidth)) << std::setprecision(rm.mPrecision);
+	os << "REMARK " << std::setw(3) << std::right << N << " " << rm.mDesc << (rm.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(std::abs(rm.mWidth)) << std::setprecision(rm.mPrecision);
 	return os;
 }

@@ -824,7 +825,7 @@ struct SEP

 std::ostream& operator<<(std::ostream& os, SEP&& sep)
 {
-	os << sep.mText << (sep.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(abs(sep.mWidth)) << std::setprecision(sep.mPrecision);
+	os << sep.mText << (sep.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(std::abs(sep.mWidth)) << std::setprecision(sep.mPrecision);
 	return os;
 }

@@ -2329,7 +2330,8 @@ void WriteRemark200(std::ostream& pdbFile, Datablock& db)
 	}
 	catch (const std::exception& ex)
 	{
-		std::cerr << ex.what() << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << ex.what() << std::endl;
 	}
 }

@@ -2390,7 +2392,8 @@ void WriteRemark280(std::ostream& pdbFile, Datablock& db)
 	}
 	catch (const std::exception& ex)
 	{
-		std::cerr << ex.what() << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << ex.what() << std::endl;
 	}
 }

--- a/src/CifParser.cpp
+++ b/src/CifParser.cpp
--- a/src/CifUtils.cpp
+++ b/src/CifUtils.cpp
@@ -123,11 +123,12 @@ const uint8_t kCharToLowerMap[256] =

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

-bool iequals(const std::string &a, const std::string &b)
+bool iequals(std::string_view a, std::string_view b)
 {
 	bool result = a.length() == b.length();
-	for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end() and bi != b.end(); ++ai, ++bi)
-		result = tolower(*ai) == tolower(*bi);
+	for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end(); ++ai, ++bi)
+		result = kCharToLowerMap[uint8_t(*ai)] == kCharToLowerMap[uint8_t(*bi)];
+		// result = tolower(*ai) == tolower(*bi);
 	return result;
 }

@@ -140,7 +141,7 @@ bool iequals(const char *a, const char *b)
 	return result and *a == *b;
 }

-int icompare(const std::string &a, const std::string &b)
+int icompare(std::string_view a, std::string_view b)
 {
 	int d = 0;
 	auto ai = a.begin(), bi = b.begin();
@@ -193,7 +194,7 @@ std::string toLowerCopy(const std::string &s)

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

-std::tuple<std::string, std::string> splitTagName(const std::string &tag)
+std::tuple<std::string, std::string> splitTagName(std::string_view tag)
 {
 	if (tag.empty())
 		throw std::runtime_error("empty tag");
@@ -1236,7 +1237,7 @@ std::filesystem::path gDataDir;

 void addDataDirectory(std::filesystem::path dataDir)
 {
-	if (VERBOSE and not fs::exists(dataDir))
+	if (VERBOSE > 0 and not fs::exists(dataDir))
 		std::cerr << "The specified data directory " << dataDir << " does not exist" << std::endl;
 	gDataDir = dataDir;
 }
--- a/src/CifValidator.cpp
+++ b/src/CifValidator.cpp
@@ -24,32 +24,39 @@
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

+#include <fstream>
+#include <filesystem>
+
 #include <boost/algorithm/string.hpp>
+#include <boost/iostreams/filtering_stream.hpp>
+#include <boost/iostreams/filter/gzip.hpp>

 #include "cif++/Cif++.hpp"
 #include "cif++/CifParser.hpp"
 #include "cif++/CifValidator.hpp"

 namespace ba = boost::algorithm;
+namespace fs = std::filesystem;
+namespace io = boost::iostreams;

 extern int VERBOSE;

 namespace cif
 {

-ValidationError::ValidationError(const std::string& msg)
+ValidationError::ValidationError(const std::string &msg)
 	: mMsg(msg)
 {
 }

-ValidationError::ValidationError(const std::string& cat, const std::string& item, const std::string& msg)
+ValidationError::ValidationError(const std::string &cat, const std::string &item, const std::string &msg)
 	: mMsg("When validating _" + cat + '.' + item + ": " + msg)
 {
 }

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

-DDL_PrimitiveType mapToPrimitiveType(const std::string& s)
+DDL_PrimitiveType mapToPrimitiveType(std::string_view s)
 {
 	DDL_PrimitiveType result;
 	if (iequals(s, "char"))
@@ -65,10 +72,10 @@ DDL_PrimitiveType mapToPrimitiveType(const std::string& s)

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

-int ValidateType::compare(const char* a, const char* b) const
+int ValidateType::compare(const char *a, const char *b) const
 {
 	int result = 0;
-	
+
 	if (*a == 0)
 		result = *b == 0 ? 0 : -1;
 	else if (*b == 0)
@@ -83,7 +90,7 @@ int ValidateType::compare(const char* a, const char* b) const
 				{
 					double da = strtod(a, nullptr);
 					double db = strtod(b, nullptr);
-					
+
 					auto d = da - db;
 					if (std::abs(d) > std::numeric_limits<double>::epsilon())
 					{
@@ -94,13 +101,13 @@ int ValidateType::compare(const char* a, const char* b) const
 					}
 					break;
 				}
-				
+
 				case DDL_PrimitiveType::UChar:
 				case DDL_PrimitiveType::Char:
 				{
 					// CIF is guaranteed to have ascii only, therefore this primitive code will do
 					// also, we're collapsing spaces
-					
+
 					auto ai = a, bi = b;
 					for (;;)
 					{
@@ -115,7 +122,7 @@ int ValidateType::compare(const char* a, const char* b) const
 							result = 1;
 							break;
 						}
-						
+
 						char ca = *ai;
 						char cb = *bi;

@@ -124,12 +131,12 @@ int ValidateType::compare(const char* a, const char* b) const
 							ca = tolower(ca);
 							cb = tolower(cb);
 						}
-						
+
 						result = ca - cb;
-						
+
 						if (result != 0)
 							break;
-						
+
 						if (ca == ' ')
 						{
 							while (ai[1] == ' ')
@@ -137,21 +144,21 @@ int ValidateType::compare(const char* a, const char* b) const
 							while (bi[1] == ' ')
 								++bi;
 						}
-						
+
 						++ai;
 						++bi;
 					}
-					
+
 					break;
 				}
 			}
 		}
-		catch (const std::invalid_argument& ex)
+		catch (const std::invalid_argument &ex)
 		{
 			result = 1;
 		}
 	}
-	
+
 	return result;
 }

@@ -165,13 +172,13 @@ int ValidateType::compare(const char* a, const char* b) const
 //
 //	if (mType == nullptr and parent != nullptr)
 //		mType = parent->mType;
-//		
+//
 //	if (parent != nullptr)
 //	{
 //		mLinked.push_back({parent, parentItem, childItem});
 //
 //		parent->mChildren.insert(this);
-////	
+////
 ////		if (mCategory->mKeys == std::vector<std::string>{mTag})
 ////			parent->mForeignKeys.insert(this);
 //	}
@@ -194,7 +201,7 @@ void ValidateItem::operator()(std::string value) const

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

-void ValidateCategory::addItemValidator(ValidateItem&& v)
+void ValidateCategory::addItemValidator(ValidateItem &&v)
 {
 	if (v.mMandatory)
 		mMandatoryFields.insert(v.mTag);
@@ -206,10 +213,10 @@ void ValidateCategory::addItemValidator(ValidateItem&& v)
 		std::cout << "Could not add validator for item " << v.mTag << " to category " << mName << std::endl;
 }

-const ValidateItem* ValidateCategory::getValidatorForItem(std::string tag) const
+const ValidateItem *ValidateCategory::getValidatorForItem(std::string_view tag) const
 {
-	const ValidateItem* result = nullptr;
-	auto i = mItemValidators.find(ValidateItem{tag});
+	const ValidateItem *result = nullptr;
+	auto i = mItemValidators.find(ValidateItem{std::string(tag)});
 	if (i != mItemValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -219,26 +226,29 @@ const ValidateItem* ValidateCategory::getValidatorForItem(std::string tag) const

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

-Validator::Validator()
+Validator::Validator(std::string_view name, std::istream &is)
+	: mName(name)
 {
+	DictParser p(*this, is);
+	p.loadDictionary();
 }

 Validator::~Validator()
 {
 }

-void Validator::addTypeValidator(ValidateType&& v)
+void Validator::addTypeValidator(ValidateType &&v)
 {
 	auto r = mTypeValidators.insert(std::move(v));
 	if (not r.second and VERBOSE > 4)
 		std::cout << "Could not add validator for type " << v.mName << std::endl;
 }

-const ValidateType* Validator::getValidatorForType(std::string typeCode) const
+const ValidateType *Validator::getValidatorForType(std::string_view typeCode) const
 {
-	const ValidateType* result = nullptr;
-	
-	auto i = mTypeValidators.find(ValidateType{ typeCode, DDL_PrimitiveType::Char, boost::regex() });
+	const ValidateType *result = nullptr;
+
+	auto i = mTypeValidators.find(ValidateType{std::string(typeCode), DDL_PrimitiveType::Char, boost::regex()});
 	if (i != mTypeValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -246,17 +256,17 @@ const ValidateType* Validator::getValidatorForType(std::string typeCode) const
 	return result;
 }

-void Validator::addCategoryValidator(ValidateCategory&& v)
+void Validator::addCategoryValidator(ValidateCategory &&v)
 {
 	auto r = mCategoryValidators.insert(std::move(v));
 	if (not r.second and VERBOSE > 4)
 		std::cout << "Could not add validator for category " << v.mName << std::endl;
 }

-const ValidateCategory* Validator::getValidatorForCategory(std::string category) const
+const ValidateCategory *Validator::getValidatorForCategory(std::string_view category) const
 {
-	const ValidateCategory* result = nullptr;
-	auto i = mCategoryValidators.find(ValidateCategory{category});
+	const ValidateCategory *result = nullptr;
+	auto i = mCategoryValidators.find(ValidateCategory{std::string(category)});
 	if (i != mCategoryValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -264,16 +274,16 @@ const ValidateCategory* Validator::getValidatorForCategory(std::string category)
 	return result;
 }

-ValidateItem* Validator::getValidatorForItem(std::string tag) const
+ValidateItem *Validator::getValidatorForItem(std::string_view tag) const
 {
-	ValidateItem* result = nullptr;
-	
+	ValidateItem *result = nullptr;
+
 	std::string cat, item;
 	std::tie(cat, item) = splitTagName(tag);

-	auto* cv = getValidatorForCategory(cat);
+	auto *cv = getValidatorForCategory(cat);
 	if (cv != nullptr)
-		result = const_cast<ValidateItem*>(cv->getValidatorForItem(item));
+		result = const_cast<ValidateItem *>(cv->getValidatorForItem(item));

 	if (result == nullptr and VERBOSE > 4)
 		std::cout << "No validator for item " << tag << std::endl;
@@ -281,15 +291,15 @@ ValidateItem* Validator::getValidatorForItem(std::string tag) const
 	return result;
 }

-void Validator::addLinkValidator(ValidateLink&& v)
+void Validator::addLinkValidator(ValidateLink &&v)
 {
 	assert(v.mParentKeys.size() == v.mChildKeys.size());
 	if (v.mParentKeys.size() != v.mChildKeys.size())
 		throw std::runtime_error("unequal number of keys for parent and child in link");
-	
+
 	auto pcv = getValidatorForCategory(v.mParentCategory);
 	auto ccv = getValidatorForCategory(v.mChildCategory);
-	
+
 	if (pcv == nullptr)
 		throw std::runtime_error("unknown parent category " + v.mParentCategory);

@@ -299,53 +309,127 @@ void Validator::addLinkValidator(ValidateLink&& v)
 	for (size_t i = 0; i < v.mParentKeys.size(); ++i)
 	{
 		auto piv = pcv->getValidatorForItem(v.mParentKeys[i]);
-	
+
 		if (piv == nullptr)
 			throw std::runtime_error("unknown parent tag _" + v.mParentCategory + '.' + v.mParentKeys[i]);

 		auto civ = ccv->getValidatorForItem(v.mChildKeys[i]);
 		if (civ == nullptr)
 			throw std::runtime_error("unknown child tag _" + v.mChildCategory + '.' + v.mChildKeys[i]);
-		
+
 		if (civ->mType == nullptr and piv->mType != nullptr)
-			const_cast<ValidateItem*>(civ)->mType = piv->mType;
-	}		
-	
+			const_cast<ValidateItem *>(civ)->mType = piv->mType;
+	}
+
 	mLinkValidators.emplace_back(std::move(v));
 }

-std::vector<const ValidateLink*> Validator::getLinksForParent(const std::string& category) const
+std::vector<const ValidateLink *> Validator::getLinksForParent(std::string_view category) const
 {
-	std::vector<const ValidateLink*> result;
+	std::vector<const ValidateLink *> result;

-	for (auto& l: mLinkValidators)
+	for (auto &l : mLinkValidators)
 	{
 		if (l.mParentCategory == category)
 			result.push_back(&l);
 	}
-	
+
 	return result;
 }

-std::vector<const ValidateLink*> Validator::getLinksForChild(const std::string& category) const
+std::vector<const ValidateLink *> Validator::getLinksForChild(std::string_view category) const
 {
-	std::vector<const ValidateLink*> result;
+	std::vector<const ValidateLink *> result;

-	for (auto& l: mLinkValidators)
+	for (auto &l : mLinkValidators)
 	{
 		if (l.mChildCategory == category)
 			result.push_back(&l);
 	}
-	
+
 	return result;
 }

-void Validator::reportError(const std::string& msg, bool fatal)
+void Validator::reportError(const std::string &msg, bool fatal) const
 {
 	if (mStrict or fatal)
 		throw ValidationError(msg);
-	else if (VERBOSE)
+	else if (VERBOSE > 0)
 		std::cerr << msg << std::endl;
 }

+// --------------------------------------------------------------------
+
+ValidatorFactory ValidatorFactory::sInstance;
+
+ValidatorFactory::ValidatorFactory()
+{
 }
+
+const Validator &ValidatorFactory::operator[](std::string_view dictionary)
+{
+	std::lock_guard lock(mMutex);
+
+	for (auto &validator : mValidators)
+	{
+		if (iequals(validator.mName, dictionary))
+			return validator;
+	}
+
+	// not found, add it
+
+	fs::path dict_name(dictionary);
+
+	auto data = loadResource(dictionary);
+
+	if (not data and dict_name.extension().string() != ".dic")
+		data = loadResource(dict_name.parent_path() / (dict_name.filename().string() + ".dic"));
+
+	if (data)
+		mValidators.emplace_back(dictionary, *data);
+	else
+	{
+		// might be a compressed dictionary on disk
+		fs::path p = dictionary;
+		if (p.extension() == ".dic")
+			p = p.parent_path() / (p.filename().string() + ".gz");
+		else
+			p = p.parent_path() / (p.filename().string() + ".dic.gz");
+
+#if defined(CACHE_DIR) and defined(DATA_DIR)
+		if (not fs::exists(p))
+		{
+			for (const char *dir : {CACHE_DIR, DATA_DIR})
+			{
+				auto p2 = fs::path(dir) / p;
+				if (fs::exists(p2))
+				{
+					swap(p, p2);
+					break;
+				}
+			}
+		}
+#endif
+
+		if (fs::exists(p))
+		{
+			std::ifstream file(p, std::ios::binary);
+			if (not file.is_open())
+				throw std::runtime_error("Could not open dictionary (" + p.string() + ")");
+
+			io::filtering_stream<io::input> in;
+			in.push(io::gzip_decompressor());
+			in.push(file);
+
+			mValidators.emplace_back(dictionary, in);
+		}
+		else
+			throw std::runtime_error("Dictionary not found or defined (" + dict_name.string() + ")");
+	}
+
+	assert(iequals(mValidators.back().mName, dictionary));
+
+	return mValidators.back();
+}
+
+} // namespace cif
--- a/src/Compound.cpp
+++ b/src/Compound.cpp
@@ -128,6 +128,8 @@ Compound::Compound(cif::Datablock &db)
 	// The name should not contain newline characters since that triggers validation errors later on
 	ba::replace_all(mName, "\n", "");

+	mGroup = "non-polymer";
+
 	auto &chemCompAtom = db["chem_comp_atom"];
 	for (auto row : chemCompAtom)
 	{
@@ -151,10 +153,11 @@ Compound::Compound(cif::Datablock &db)
 	}
 }

-Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type)
+Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group)
 	: mID(id)
 	, mName(name)
 	, mType(type)
+	, mGroup(group)
 {
 	auto &chemCompAtom = db["chem_comp_atom"];
 	for (auto row : chemCompAtom)
@@ -190,7 +193,7 @@ Compound::Compound(cif::Datablock &db, const std::string &id, const std::string
 			bond.type = BondType::delo;
 		else
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Unimplemented chem_comp_bond.type " << btype << " in " << id << std::endl;
 			bond.type = BondType::sing;
 		}
@@ -400,7 +403,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const std::filesystem::path &file, std:

 			if (cif::iequals(id, "gly"))
 				type = "peptide linking";
-			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
+			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide") or cif::iequals(group, "p-peptide"))
 				type = "L-peptide linking";
 			else if (cif::iequals(group, "DNA"))
 				type = "DNA linking";
@@ -411,7 +414,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const std::filesystem::path &file, std:

 			auto &db = cifFile["comp_" + id];

-			mCompounds.push_back(new Compound(db, id, name, type));
+			mCompounds.push_back(new Compound(db, id, name, type, group));
 		}
 	}
 	else
@@ -517,7 +520,7 @@ Compound *CCDCompoundFactoryImpl::create(const std::string &id)
 		}
 	}

-	if (result == nullptr and cif::VERBOSE)
+	if (result == nullptr and cif::VERBOSE > 0)
 		std::cerr << "Could not locate compound " << id << " in the CCD components file" << std::endl;

 	return result;
@@ -608,7 +611,7 @@ Compound *CCP4CompoundFactoryImpl::create(const std::string &id)

 			if (cif::iequals(id, "gly"))
 				type = "peptide linking";
-			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
+			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide") or cif::iequals(group, "p-peptide"))
 				type = "L-peptide linking";
 			else if (cif::iequals(group, "DNA"))
 				type = "DNA linking";
@@ -617,7 +620,7 @@ Compound *CCP4CompoundFactoryImpl::create(const std::string &id)
 			else
 				type = "non-polymer";

-			mCompounds.push_back(new Compound(db, id, name, type));
+			mCompounds.push_back(new Compound(db, id, name, type, group));
 			result = mCompounds.back();
 		}
 	}
@@ -642,13 +645,13 @@ CompoundFactory::CompoundFactory()
 	auto ccd = cif::loadResource("components.cif");
 	if (ccd)
 		mImpl.reset(new CCDCompoundFactoryImpl(mImpl));
-	else if (cif::VERBOSE)
+	else if (cif::VERBOSE > 0)
 		std::cerr << "CCD components.cif file was not found" << std::endl;

 	const char *clibd_mon = getenv("CLIBD_MON");
 	if (clibd_mon != nullptr and fs::is_directory(clibd_mon))
 		mImpl.reset(new CCP4CompoundFactoryImpl(clibd_mon));
-	else if (cif::VERBOSE)
+	else if (cif::VERBOSE > 0)
 		std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;

 }
@@ -692,7 +695,8 @@ void CompoundFactory::setDefaultDictionary(const std::filesystem::path &inDictFi
 	}
 	catch (const std::exception &)
 	{
-		std::cerr << "Error loading dictionary " << inDictFile << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Error loading dictionary " << inDictFile << std::endl;
 		throw;
 	}
 }
@@ -712,7 +716,8 @@ void CompoundFactory::pushDictionary(const std::filesystem::path &inDictFile)
 	}
 	catch (const std::exception &)
 	{
-		std::cerr << "Error loading dictionary " << inDictFile << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Error loading dictionary " << inDictFile << std::endl;
 		throw;
 	}
 }
--- a/src/PDB2Cif.cpp
+++ b/src/PDB2Cif.cpp
--- a/src/PDB2CifRemark3.cpp
+++ b/src/PDB2CifRemark3.cpp
@@ -1320,7 +1320,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab

 	if (line != "REFINEMENT.")
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "Unexpected data in REMARK 3" << std::endl;
 		return false;
 	}
@@ -1332,7 +1332,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab

 	if (not std::regex_match(line, m, rxp))
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "Expected valid PROGRAM line in REMARK 3" << std::endl;
 		return false;
 	}
@@ -1367,8 +1367,9 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 		}
 		catch(const std::exception& e)
 		{
-			std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
-					  << e.what() << '\n';
+			if (cif::VERBOSE >= 0)
+				std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
+						  << e.what() << '\n';
 			score = 0;
 		}

@@ -1411,7 +1412,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 			tryParser(new TNT_Remark3Parser(program, expMethod, r, db));
 		else if (ba::starts_with(program, "X-PLOR"))
 			tryParser(new XPLOR_Remark3Parser(program, expMethod, r, db));
-		else if (cif::VERBOSE)
+		else if (cif::VERBOSE > 0)
 			std::cerr << "Skipping unknown program (" << program << ") in REMARK 3" << std::endl;
 	}

@@ -1420,7 +1421,8 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 	bool guessProgram = scores.empty() or scores.front().score < 0.9f;;
 	if (guessProgram)
 	{
-		std::cerr << "Unknown or untrusted program in REMARK 3, trying all parsers to see if there is a match" << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Unknown or untrusted program in REMARK 3, trying all parsers to see if there is a match" << std::endl;
 		
 		tryParser(new BUSTER_TNT_Remark3Parser("BUSTER-TNT", expMethod, r, db));
 		tryParser(new CNS_Remark3Parser("CNS", expMethod, r, db));
@@ -1444,7 +1446,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab

 		auto& best = scores.front();

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "Choosing " << best.parser->program() << " version '" << best.parser->version() << "' as refinement program. Score = " << best.score << std::endl;

 		auto& software = db["software"];
--- a/src/Point.cpp
+++ b/src/Point.cpp
@@ -28,11 +28,248 @@
 #include <valarray>

 #include "cif++/Point.hpp"
-#include "cif++/Matrix.hpp"

 namespace mmcif
 {

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

 Quaternion Normalize(Quaternion q)
@@ -102,14 +339,14 @@ Point CenterPoints(std::vector<Point>& Points)
 	return t;
 }

-Point Centroid(std::vector<Point>& Points)
+Point Centroid(const std::vector<Point>& pts)
 {
 	Point result;
 	
-	for (Point& pt : Points)
+	for (auto &pt : pts)
 		result += pt;
 	
-	result /= static_cast<float>(Points.size());
+	result /= static_cast<float>(pts.size());
 	
 	return result;
 }
@@ -175,7 +412,7 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
 Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
 {
 	// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
-	Matrix<double> M(3, 3, 0);
+	Matrix M(3, 3, 0);

 	for (uint32_t i = 0; i < pa.size(); ++i)
 	{
@@ -188,7 +425,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 	}
 	
 	// Now calculate N, a symmetric 4x4 Matrix
-	SymmetricMatrix<double> N(4);
+	SymmetricMatrix N(4);
 	
 	N(0, 0) =  M(0, 0) + M(1, 1) + M(2, 2);
 	N(0, 1) =  M(1, 2) - M(2, 1);
@@ -225,6 +462,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 					M(1, 2) * M(2, 0) * M(0, 1) +
 					M(2, 1) * M(1, 0) * M(0, 2));
 	
+	// E is the determinant of N:
 	double E = 
 		(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
 		(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
@@ -234,47 +472,22 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 		(N(0,2) * N(1,3) - N(2,1) * N(0,3)) * (N(0,2) * N(1,3) - N(2,1) * N(0,3));
 	
 	// solve quartic
-	double lm = LargestDepressedQuarticSolution(C, D, E);
+	double lambda = LargestDepressedQuarticSolution(C, D, E);
 	
 	// calculate t = (N - λI)
-	Matrix<double> li = IdentityMatrix<double>(4) * lm;
-	Matrix<double> t = N - li;
+	Matrix t = N - IdentityMatrix(4) * lambda;
 	
 	// calculate a Matrix of cofactors for t
-	Matrix<double> cf(4, 4);
+	Matrix cf = Cofactors(t);

-	const uint32_t ixs[4][3] =
+	int maxR = 0;
+	for (int r = 1; r < 4; ++r)
 	{
-		{ 1, 2, 3 },
-		{ 0, 2, 3 },
-		{ 0, 1, 3 },
-		{ 0, 1, 2 }
-	};
-
-	uint32_t maxR = 0;
-	for (uint32_t r = 0; r < 4; ++r)
-	{
-		const uint32_t* ir = ixs[r];
-		
-		for (uint32_t c = 0; c < 4; ++c)
-		{
-			const uint32_t* ic = ixs[c];
-
-			cf(r, c) =
-				t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
-				t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
-				t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
-				t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
-				t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
-				t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
-		}
-		
-		if (r > maxR and cf(r, 0) > cf(maxR, 0))
+		if (std::abs(cf(r, 0)) > std::abs(cf(maxR, 0)))
 			maxR = r;
 	}
 	
-	// NOTE the negation of the y here, why? Maybe I swapped r/c above?
-	Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
+	Quaternion q(cf(maxR, 0), cf(maxR, 1), cf(maxR, 2), cf(maxR, 3));
 	q = Normalize(q);
 	
 	return q;
--- a/src/Secondary.cpp
+++ b/src/Secondary.cpp
--- a/src/Structure.cpp
+++ b/src/Structure.cpp
--- a/src/Symmetry.cpp
+++ b/src/Symmetry.cpp
@@ -90,4 +90,66 @@ int GetSpacegroupNumber(std::string spacegroup)
 	return result;
 }

+// --------------------------------------------------------------------
+
+int GetSpacegroupNumber(std::string spacegroup, SpacegroupName type)
+{
+	if (spacegroup == "P 21 21 2 A")
+		spacegroup = "P 21 21 2 (a)";
+	else if (spacegroup.empty())
+		throw std::runtime_error("No spacegroup, cannot continue");
+
+	int result = 0;
+
+	if (type == SpacegroupName::full)
+	{
+		const size_t N = kNrOfSpaceGroups;
+		int32_t L = 0, R = static_cast<int32_t>(N - 1);
+		while (L <= R)
+		{
+			int32_t i = (L + R) / 2;
+
+			int d = spacegroup.compare(kSpaceGroups[i].name);
+
+			if (d > 0)
+				L = i + 1;
+			else if (d < 0)
+				R = i - 1;
+			else
+			{
+				result = kSpaceGroups[i].nr;
+				break;
+			}
+		}
+	}
+	else if (type == SpacegroupName::xHM)
+	{
+		for (auto &sg : kSpaceGroups)
+		{
+			if (sg.xHM == spacegroup)
+			{
+				result = sg.nr;
+				break;
+			}
+		}
+	}
+	else
+	{
+		for (auto &sg : kSpaceGroups)
+		{
+			if (sg.Hall == spacegroup)
+			{
+				result = sg.nr;
+				break;
+			}
+		}
+	}
+
+	// not found, see if we can find a match based on xHM name
+	if (result == 0)
+		throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
+	
+	return result;
+}
+
 }
--- a/src/TlsParser.cpp
+++ b/src/TlsParser.cpp
@@ -248,7 +248,7 @@ struct TLSSelectionNot : public TLSSelection
 		for (auto& r: residues)
 			r.selected = not r.selected;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "NOT" << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -267,7 +267,7 @@ struct TLSSelectionAll : public TLSSelection
 		for (auto& r: residues)
 			r.selected = true;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "ALL" << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -287,7 +287,7 @@ struct TLSSelectionChain : public TLSSelectionAll
 		for (auto& r: residues)
 			r.selected = allChains or r.chainID == m_chain;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "CHAIN " << m_chain << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -307,7 +307,7 @@ struct TLSSelectionResID : public TLSSelectionAll
 		for (auto& r: residues)
 			r.selected = r.seqNr == m_seq_nr and r.iCode == m_icode;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "ResID " << m_seq_nr << (m_icode ? std::string { m_icode} : "") << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -331,7 +331,7 @@ struct TLSSelectionRangeSeq : public TLSSelectionAll
 						  (r.seqNr <= m_last or m_last == kResidueNrWildcard));
 		}

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Range " << m_first << ':' << m_last << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -374,7 +374,7 @@ struct TLSSelectionRangeID : public TLSSelectionAll
 			}
 		}

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Through " << m_first << ':' << m_last << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -407,7 +407,7 @@ struct TLSSelectionUnion : public TLSSelection
 		for (auto ai = a.begin(), bi = b.begin(), ri = residues.begin(); ri != residues.end(); ++ai, ++bi, ++ri)
 			ri->selected = ai->selected or bi->selected;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Union" << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -440,7 +440,7 @@ struct TLSSelectionIntersection : public TLSSelection
 		for (auto ai = a.begin(), bi = b.begin(), ri = residues.begin(); ri != residues.end(); ++ai, ++bi, ++ri)
 			ri->selected = ai->selected and bi->selected;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Intersection" << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -462,7 +462,7 @@ struct TLSSelectionByName : public TLSSelectionAll
 		for (auto& r: residues)
 			r.selected = r.name == m_name;

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Name " << m_name << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -488,7 +488,7 @@ struct TLSSelectionByElement : public TLSSelectionAll
 		for (auto& r: residues)
 			r.selected = iequals(r.name, m_element);

-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 		{
 			std::cout << std::string(indentLevel * 2, ' ') << "Element " << m_element << std::endl;
 			DumpSelection(residues, indentLevel);
@@ -890,7 +890,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::Parse()

 	Match(pt_EOLN);
 	
-	if (extraParenthesis)
+	if (extraParenthesis and cif::VERBOSE > 0)
 		std::cerr << "WARNING: too many closing parenthesis in TLS selection statement" << std::endl;
 	
 	return result;
@@ -931,7 +931,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::ParseFactor()
 		case '(':
 			Match('(');
 			result = ParseAtomSelection();
-			if (m_lookahead == pt_EOLN)
+			if (m_lookahead == pt_EOLN and cif::VERBOSE > 0)
 				std::cerr << "WARNING: missing closing parenthesis in TLS selection statement" << std::endl;
 			else
 				Match(')');
@@ -1033,7 +1033,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::ParseFactor()
 					result.reset(new TLSSelectionRangeID(from, to, icode_from, icode_to));
 				else
 				{
-					if (cif::VERBOSE and (icode_from or icode_to))
+					if (cif::VERBOSE > 0 and (icode_from or icode_to))
 						std::cerr << "Warning, ignoring insertion codes" << std::endl;
 					
 					result.reset(new TLSSelectionRangeSeq(from, to));
@@ -1231,7 +1231,8 @@ TLSSelectionPtr TLSSelectionParserImplBuster::ParseGroup()
 				std::tie(chain2, seqNr2) = ParseAtom();
 				if (chain1 != chain2)
 				{
-					std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
+					if (cif::VERBOSE > 0)
+						std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
 					
 					TLSSelectionPtr sc1(new TLSSelectionChain(chain1));
 					TLSSelectionPtr sr1(new TLSSelectionRangeSeq(seqNr1, kResidueNrWildcard));
@@ -1289,7 +1290,7 @@ std::tuple<std::string,int> TLSSelectionParserImplBuster::ParseAtom()
 			Match(':');
 			std::string atom = m_value_s;
 			
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Warning: ignoring atom ID '" << atom << "' in TLS selection" << std::endl;
 			
 			Match(bt_IDENT);
@@ -1810,7 +1811,8 @@ class TLSSelectionParser
 		}
 		catch (const std::exception& ex)
 		{
-			std::cerr << "ParseError: " << ex.what() << std::endl;
+			if (cif::VERBOSE >= 0)
+				std::cerr << "ParseError: " << ex.what() << std::endl;
 		}
 		
 		return result;
@@ -1834,14 +1836,14 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str

 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}
 		
 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to PHENIX" << std::endl;
 			result = phenix.Parse(selection);
 		}
@@ -1852,35 +1854,35 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str

 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to BUSTER" << std::endl;
 			result = buster.Parse(selection);
 		}

 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}
 	}
 	else
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "No known program specified, trying PHENIX" << std::endl;

 		result = phenix.Parse(selection);

 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to BUSTER" << std::endl;
 			result = buster.Parse(selection);
 		}

 		if (not result)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}
--- a/test/pdb2cif-test.cpp
+++ b/test/pdb2cif-test.cpp
@@ -18,7 +18,6 @@ int main(int argc, char* argv[])
 	desc.add_options()
 		("input,i",		po::value<std::string>(),	"Input file")
 		("help,h",									"Display help message")
-		("version",									"Print version")
 		("verbose,v",								"Verbose output")
 		("debug,d",		po::value<int>(),			"Debug level (for even more verbose output)");

@@ -29,12 +28,6 @@ int main(int argc, char* argv[])
 	po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
 	po::notify(vm);

-	if (vm.count("version"))
-	{
-		std::cout << argv[0] << " version " PACKAGE_VERSION << std::endl;
-		exit(0);
-	}
-
 	if (vm.count("help") or vm.count("input") == 0)
 	{
 		std::cerr << desc << std::endl;
--- a/test/rename-compound-test.cpp
+++ b/test/rename-compound-test.cpp
@@ -20,8 +20,11 @@ int main(int argc, char* argv[])
 		if (argc == 3)
 			testdir = argv[2];

-		if (std::filesystem::exists(testdir / ".."/"data"/"ccd-subset.cif"))
-			cif::addFileResource("components.cif", testdir / ".."/"data"/"ccd-subset.cif");
+		if (std::filesystem::exists(testdir / ".." / "data" / "ccd-subset.cif"))
+			cif::addFileResource("components.cif", testdir / ".." / "data" / "ccd-subset.cif");
+
+		if (std::filesystem::exists(testdir / ".." / "rsrc" / "mmcif_pdbx_v50.dic"))
+			cif::addFileResource("mmcif_pdbx_v50.dic", testdir / ".." / "rsrc" / "mmcif_pdbx_v50.dic");

 		mmcif::CompoundFactory::instance().pushDictionary(testdir / "REA.cif");
 		mmcif::CompoundFactory::instance().pushDictionary(testdir / "RXA.cif");
--- a/test/structure-test.cpp
+++ b/test/structure-test.cpp
@@ -77,9 +77,6 @@ BOOST_AUTO_TEST_CASE(create_nonpoly_1)
 {
    cif::VERBOSE = 1;

-	// do this now, avoids the need for installing
-	cif::addFileResource("mmcif_pdbx_v50.dic", "../rsrc/mmcif_pdbx_v50.dic");
-
 	mmcif::File file;
 	file.file().loadDictionary("mmcif_pdbx_v50.dic");
 	file.createDatablock("TEST");	// create a datablock
@@ -182,3 +179,23 @@ _struct_asym.details                       ?
 				<< structure.getFile().data() << std::endl;
 	}
 }
+
+// // --------------------------------------------------------------------
+
+// BOOST_AUTO_TEST_CASE(test_load_1)
+// {
+// 	mmcif::File cf(gTestDir / "5v3g.cif.gz");
+// 	mmcif::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;
+// 		}
+// 	}
+// }
--- a/test/unit-test.cpp
+++ b/test/unit-test.cpp
@@ -33,6 +33,10 @@
 // #include "cif++/DistanceMap.hpp"
 #include "cif++/Cif++.hpp"
 #include "cif++/BondMap.hpp"
+#include "cif++/CifValidator.hpp"
+
+namespace tt = boost::test_tools;
+

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

@@ -259,8 +263,10 @@ save__cat_2.desc

    std::istream is_dict(&buffer);

+	cif::Validator validator("test", is_dict);
+
    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -387,8 +393,10 @@ save__cat_1.c

    std::istream is_dict(&buffer);

+	cif::Validator validator("test", is_dict);
+
    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -535,8 +543,10 @@ save__cat_2.desc

    std::istream is_dict(&buffer);

+	cif::Validator validator("test", is_dict);
+
    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -741,8 +751,10 @@ save__cat_2.parent_id3

    std::istream is_dict(&buffer);

+	cif::Validator validator("test", is_dict);
+
    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -963,8 +975,10 @@ cat_2 3 cat_2:cat_1:3

    std::istream is_dict(&buffer);

+	cif::Validator validator("test", is_dict);
+
    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -1389,9 +1403,10 @@ cat_2 1 '_cat_2.num'  '_cat_3.num'  cat_3
    } buffer(const_cast<char*>(dict), sizeof(dict) - 1);

    std::istream is_dict(&buffer);
+	cif::Validator validator("test", is_dict);

    cif::File f;
-    f.loadDictionary(is_dict);
+    f.setValidator(&validator);

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

@@ -1687,4 +1702,58 @@ BOOST_AUTO_TEST_CASE(reading_file_1)

 	cif::File file;
 	BOOST_CHECK_THROW(file.load(is), std::runtime_error);
-}
+}
+
+// 3d tests
+
+using namespace mmcif;
+
+BOOST_AUTO_TEST_CASE(t1)
+{
+	// std::random_device rnd;
+	// std::mt19937 gen(rnd());
+	// std::uniform_real_distribution<float> dis(0, 1);
+
+	// Quaternion q{ dis(gen), dis(gen), dis(gen), dis(gen) };
+	// q = Normalize(q);
+
+	// Quaternion q{ 0.1, 0.2, 0.3, 0.4 };
+	Quaternion q{ 0.5, 0.5, 0.5, 0.5 };
+	q = Normalize(q);
+
+	const auto &&[angle0, axis0] = QuaternionToAngleAxis(q);
+
+	std::vector<Point> p1{
+		{ 16.979, 13.301, 44.555 },
+		{ 18.150, 13.525, 43.680 },
+		{ 18.656, 14.966, 43.784 },
+		{ 17.890, 15.889, 44.078 },
+		{ 17.678, 13.270, 42.255 },
+		{ 16.248, 13.734, 42.347 },
+		{ 15.762, 13.216, 43.724 }
+	};
+
+	auto p2 = p1;
+
+	Point c1 = CenterPoints(p1);
+
+	for (auto &p : p2)
+		p.rotate(q);
+	
+	Point c2 = CenterPoints(p2);
+
+	auto q2 = AlignPoints(p1, p2);
+
+	const auto &&[angle, axis] = QuaternionToAngleAxis(q2);
+
+	BOOST_TEST(std::fmod(360 + angle, 360) == std::fmod(360 - angle0, 360), tt::tolerance(0.01));
+
+	for (auto &p : p1)
+		p.rotate(q2);
+
+	float rmsd = RMSd(p1, p2);
+
+	BOOST_TEST(rmsd < 1e-5);
+
+	// std::cout << "rmsd: " << RMSd(p1, p2) << std::endl;
+}
Author	SHA1	Message	Date
Maarten L. Hekkelman	4732004b67	Merge branch 'develop' into trunk	2022-01-12 16:41:18 +01:00
Maarten L. Hekkelman	faa9cd0431	Added another rotate/translate method to mmcif::Structure	2022-01-12 14:06:32 +01:00
Maarten L. Hekkelman	e0c3c2394d	Fix Structure::createNonPoly to add atoms...	2022-01-11 11:21:56 +01:00
Maarten L. Hekkelman	2dec584f54	clean up code	2022-01-05 15:54:23 +01:00
Maarten L. Hekkelman	5ab2ccae40	avoid calling cif::Category::size() too often	2022-01-05 15:45:27 +01:00
Maarten L. Hekkelman	1017d08626	skip updating links when changing atom location	2022-01-05 15:24:22 +01:00
Maarten L. Hekkelman	32b1bbd943	combine translate and rotate in a single call	2022-01-05 14:27:16 +01:00
Maarten L. Hekkelman	1abf31ffa5	no-validate option in cif::Row::assign	2022-01-05 14:04:28 +01:00
Maarten L. Hekkelman	aec60829d2	more quiet code	2022-01-05 11:29:10 +01:00
Maarten L. Hekkelman	888c3c38c2	Add a 'quiet' mode (cif::VERBOSE < 0)	2022-01-05 10:36:39 +01:00
Maarten L. Hekkelman	e2c4648037	clean up	2022-01-05 10:24:37 +01:00
Maarten L. Hekkelman	f7b98c0530	refactored AtomImpl	2022-01-05 10:23:15 +01:00
Maarten L. Hekkelman	d4bd3faa16	Merge branch 'profiling-structure' into trunk	2022-01-04 10:29:23 +01:00
Maarten L. Hekkelman	c4f3b1cd7b	delay loading atoms in residues	2022-01-04 09:48:41 +01:00
Maarten L. Hekkelman	74add69a83	Finish removing bzip2 support	2022-01-03 15:51:01 +01:00
Maarten L. Hekkelman	a490b19d24	version bump	2022-01-03 15:45:48 +01:00
Maarten L. Hekkelman	44cfa2c1a2	further optimisation	2022-01-03 15:19:50 +01:00
Maarten L. Hekkelman	6dd9522b3f	optimized mmcif::Atom	2022-01-03 14:32:42 +01:00
Maarten L. Hekkelman	5e352cb8e4	Removed erronous dependency in config.cmake.in	2021-12-20 13:33:06 +01:00
Maarten L. Hekkelman	2fad7315b8	make DSSP::iterator bidirectional	2021-12-15 15:12:23 +01:00
Maarten L. Hekkelman	520759dfe8	update changelog	2021-12-14 15:17:50 +01:00
Maarten L. Hekkelman	577b44ae11	Fix in processing CCP4 monomers, proline is a peptide	2021-12-14 15:16:23 +01:00
Maarten L. Hekkelman	66f742d6c0	code to facilitate DSSP	2021-12-14 15:14:45 +01:00
Maarten L. Hekkelman	7ba9f688c7	Merge branch 'develop' into trunk	2021-12-10 10:39:34 +01:00
Maarten L. Hekkelman	883f0307a2	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-12-10 10:38:43 +01:00
Maarten L. Hekkelman	c9719f873f	Merge branch 'develop' into trunk	2021-12-10 10:37:04 +01:00
Maarten L. Hekkelman	123d25f853	formatting of floating points in cif files better verbose info for differences	2021-12-10 10:35:02 +01:00
Maarten L. Hekkelman	56da42db84	formatting of floating points in cif files better verbose info for differences	2021-12-10 10:32:21 +01:00
Maarten L. Hekkelman	7f820449ca	formatting	2021-12-08 09:06:09 +01:00
Maarten L. Hekkelman	ecb2cf5f11	Fix for compiling with gcc 11.2	2021-12-08 09:03:21 +01:00
Maarten L. Hekkelman	7f27da9b3b	Fixed rename-compound-test to work when not using resources	2021-11-25 16:27:45 +01:00
Maarten L. Hekkelman	01eb499c69	attempt to fix running tests in different directory	2021-11-25 16:27:42 +01:00
Maarten L. Hekkelman	1ff6f70682	changelog update	2021-11-25 16:25:53 +01:00
Maarten L. Hekkelman	ddde996e10	strip newlines from compound names read from CCD	2021-11-25 16:24:55 +01:00
Maarten L. Hekkelman	1c9212c7e0	Fixed rename-compound-test to work when not using resources	2021-11-25 16:09:15 +01:00
Maarten L. Hekkelman	a568143991	unneeded loading of resource removed from test	2021-11-24 13:52:23 +01:00
Maarten L. Hekkelman	2b6f1bd9ee	attempt to fix running tests in different directory	2021-11-23 10:51:22 +01:00
Maarten L. Hekkelman	2527aa5ea6	correct version in cmakefile, fix structure-test to no longer load resource	2021-11-24 13:57:36 +01:00
Maarten L. Hekkelman	4c28091ecd	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-11-18 09:56:11 +01:00
Maarten L. Hekkelman	d49725423e	re-added group to Compound. Seems to be important	2021-11-18 09:56:01 +01:00
Maarten L. Hekkelman	fcb4dc61b5	fixed writing PDB files (= remove to_upper for all header lines...)	2021-11-17 08:01:24 +01:00
Maarten L. Hekkelman	b7330c074f	Fixed Structure::changeResidue to actually change the residue itself as well.	2021-11-16 08:38:14 +01:00
Maarten L. Hekkelman	e8f4123030	strip newlines from names in Compound	2021-11-15 12:37:45 +01:00
Maarten L. Hekkelman	975057c4c4	Fixed bug in structure::changeresidue when removing atoms	2021-11-15 11:32:50 +01:00
Maarten L. Hekkelman	a0e01668d1	take largest in value for best quaternion	2021-11-15 10:22:56 +01:00
Maarten L. Hekkelman	2c77491416	cleaner implementation of matrices	2021-11-12 10:36:38 +01:00
Maarten L. Hekkelman	be19e4a9cb	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-11-12 09:20:01 +01:00
Maarten L. Hekkelman	61ce91a9d7	using expression templates for matrices	2021-11-12 09:17:21 +01:00
Maarten L. Hekkelman	18f1d07e85	clean up code	2021-11-12 08:28:10 +01:00
Maarten L. Hekkelman	b596976194	correct implementation of alignpoints	2021-11-12 08:15:33 +01:00
Maarten L. Hekkelman	1f6b86d516	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-11-11 21:41:05 +01:00
Maarten L. Hekkelman	31499b977d	Fix the 3d alignment code	2021-11-11 21:39:52 +01:00
Maarten L. Hekkelman	f83850e380	git and revisions	2021-11-11 09:42:22 +01:00
Maarten L. Hekkelman	7f39d401e2	Optimised assigning data	2021-10-20 14:55:30 +02:00
Maarten L. Hekkelman	af412c284d	clean up	2021-10-20 12:27:27 +02:00
Maarten L. Hekkelman	874cd3bae5	Fix symmetry lookup	2021-10-20 11:36:31 +02:00
Maarten L. Hekkelman	ea28ebdd13	optimized caching of items in mmcif::Atom	2021-10-19 16:10:59 +02:00
Maarten L. Hekkelman	3ba468933f	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-10-18 17:17:20 +02:00
Maarten L. Hekkelman	45f33e4bea	Merge branch 'trunk' into develop	2021-10-18 11:12:56 +02:00
Maarten L. Hekkelman	cb3443ffb1	Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop	2021-10-18 09:42:30 +02:00
Maarten L. Hekkelman	7513cc1947	Merge branch 'trunk' into develop	2021-10-14 09:30:19 +02:00
Maarten L. Hekkelman	ab2dd4b75f	Merge branch 'trunk' into develop	2021-10-13 14:03:29 +02:00
Maarten L. Hekkelman	8bbcba76cf	Performance increase by using std::string_view Updated Structure::changeResidue	2021-10-12 16:04:27 +02:00