- Added more const members

- change PDB writing interface

.gitignore

@@ -11,4 +11,5 @@ data/components.cif*
 CMakeSettings.json
 msvc/
 Testing/
-
+rsrc/feature-request.txt
+test/1cbs.cif

CMakeLists.txt

@@ -25,7 +25,7 @@
 cmake_minimum_required(VERSION 3.16)
 
 # set the project name
-project(cifpp VERSION 2.0.2 LANGUAGES CXX)
+project(cifpp VERSION 4.0.1 LANGUAGES CXX)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
 
@@ -37,17 +37,13 @@ include(CheckIncludeFiles)
 include(CheckLibraryExists)
 include(CMakePackageConfigHelpers)
 include(Dart)
-include(FindFilesystem)
 include(GenerateExportHeader)
 
 set(CXX_EXTENSIONS OFF)
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
-find_package(Filesystem REQUIRED)
-
 if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
-	# https://stackoverflow.com/questions/63902528/program-crashes-when-filesystempath-is-destroyed
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-unused-parameter -Wno-missing-field-initializers")
 elseif(MSVC)
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4")
@@ -74,6 +70,8 @@ if(BUILD_FOR_CCP4)
 			set(BUILD_SHARED_LIBS ON)
 		endif()
 	endif("$ENV{CCP4}" STREQUAL "" OR NOT EXISTS $ENV{CCP4})
+
+	add_definitions(-DCCP4=1)
 endif()
 
 # Check if CCP4 is available
@@ -122,22 +120,6 @@ if(MSVC)
 		message(STATUS "The library and auxiliary files will be installed in $ENV{LOCALAPPDATA}/${PROJECT_NAME}")
 		set(CMAKE_INSTALL_PREFIX "$ENV{LOCALAPPDATA}/${PROJECT_NAME}" CACHE PATH "..." FORCE)
 	endif()
-
-	# Find out the processor type for the target
-	if(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "AMD64")
-		set(COFF_TYPE "x64")
-	elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "i386")
-		set(COFF_TYPE "x86")
-	elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "ARM64")
-		set(COFF_TYPE "arm64")
-	else()
-		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()
 
 if(UNIX AND NOT APPLE AND NOT BUILD_FOR_CCP4 AND CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
@@ -151,12 +133,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,35 +155,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)
+	list(APPEND CIFPP_REQUIRED_LIBRARIES ZLIB::ZLIB)
 endif()
 
+include(FindFilesystem)
+list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})
+
+include(FindAtomic)
+list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPATOMIC_LIBRARY})
+
 # Create a revision file, containing the current git version info
-
-find_package(Git)
-if(GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git")
-	include(GetGitRevisionDescription)
-    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()
-
-# generate version.h
-string(TIMESTAMP BUILD_DATE_TIME "%Y-%m-%dT%H:%M:%SZ" UTC)
-configure_file("${CMAKE_SOURCE_DIR}/src/revision.hpp.in" "${CMAKE_BINARY_DIR}/revision.hpp" @ONLY)
+include(VersionString)
+write_version_header("LibCIFPP")
 
 # SymOp data table
 if(CIFPP_RECREATE_SYMOP_DATA)
@@ -209,7 +178,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} ${CIFPP_REQUIRED_LIBRARIES})
 	if(Boost_INCLUDE_DIR)
 		target_include_directories(symop-map-generator PUBLIC ${Boost_INCLUDE_DIR})
 	endif()
@@ -256,7 +225,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 +249,8 @@ 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 PUBLIC Threads::Threads Boost::regex Boost::iostreams ${CIFPP_REQUIRED_LIBRARIES})
+# target_link_libraries(cifpp PRIVATE)
 
 if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
     target_link_options(cifpp PRIVATE -undefined dynamic_lookup)
@@ -341,6 +310,13 @@ install(TARGETS cifpp
 	RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
 	INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
 
+if(MSVC AND BUILD_SHARED_LIBS)
+	install(
+		FILES $<TARGET_PDB_FILE:${PROJECT_NAME}>
+		DESTINATION ${CMAKE_INSTALL_LIBDIR}
+		OPTIONAL)
+endif()
+
 install(EXPORT cifppTargets
 	FILE "cifppTargets.cmake"
 	NAMESPACE cifpp::
@@ -414,32 +390,30 @@ 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
 		structure
+		sugar
 		unit)
 
 	foreach(CIFPP_TEST IN LISTS CIFPP_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 )
+
+		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 +423,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()

Config.cmake.in

@@ -1,10 +1,10 @@
 @PACKAGE_INIT@
 
 include(CMakeFindDependencyMacro)
-find_dependency(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_options)
+find_dependency(Threads)
+find_dependency(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex)
 if(NOT WIN32)
 find_dependency(ZLIB)
-find_dependency(BZip2)
 endif()
 
 INCLUDE("${CMAKE_CURRENT_LIST_DIR}/cifppTargets.cmake")

changelog

@@ -1,3 +1,48 @@
+Version 4.0.1
+- Added a bunch of const methods to Datablock and Category.
+- Changed PDB writing interface to accept Datablock instead of File.
+
+Version 4.0.0
+- getResidue in mmcif::Structure now requires both a
+  sequence ID and an auth sequence ID. As a result the code was cleaned
+  up considerably.
+
+Version 3.0.5
+- mmcif::Structure redesign. It is now a wrapper around a cif::Datablock.
+
+Version 3.0.4
+- Fix in mmCIF parser, now correctly handles the unquoted
+  string ??
+
+Version 3.0.3
+- Better configuration checks, for atomic e.g.
+- Fixed a problem introduced in refactoring mmcif::Atom
+- Version string creation
+
+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.5
+- Backporting updated CMakeLists.txt file
+
+Version 2.0.4
+- Reverted a too strict test when reading cif files.
+
+Version 2.0.3
+- Fixed reading mmCIF files where model numbers are used and
+  model number 1 is missing.
+
 Version 2.0.2
 - Added configuration flag to disable downloading CCD data during build
   Note that there are now two flags for CCD data:

cmake/FindAtomic.cmake

@@ -0,0 +1,64 @@
+# Simple check to see if we need a library for std::atomic
+
+if(TARGET std::atomic)
+	return()
+endif()
+
+cmake_minimum_required(VERSION 3.10)
+
+include(CMakePushCheckState)
+include(CheckIncludeFileCXX)
+include(CheckCXXSourceRuns)
+
+cmake_push_check_state()
+
+set(CMAKE_CXX_STANDARD 17)
+
+check_include_file_cxx("atomic" _CXX_ATOMIC_HAVE_HEADER)
+mark_as_advanced(_CXX_ATOMIC_HAVE_HEADER)
+
+set(code [[
+#include <atomic>
+int main(int argc, char** argv) {
+  std::atomic<long long> s;
+  ++s;
+  return 0;
+}
+]])
+
+check_cxx_source_runs("${code}" _CXX_ATOMIC_BUILTIN)
+
+if(_CXX_ATOMIC_BUILTIN)
+	set(_found 1)
+else()
+  list(APPEND CMAKE_REQUIRED_LIBRARIES atomic)
+  list(APPEND FOLLY_LINK_LIBRARIES atomic)
+
+  check_cxx_source_runs("${code}" _CXX_ATOMIC_LIB_NEEDED)
+  if (NOT _CXX_ATOMIC_LIB_NEEDED)
+    message(FATAL_ERROR "unable to link C++ std::atomic code: you may need \
+      to install GNU libatomic")
+  else()
+	set(_found 1)
+  endif()
+endif()
+
+if(_found)
+	add_library(std::atomic INTERFACE IMPORTED)
+	set_property(TARGET std::atomic APPEND PROPERTY INTERFACE_COMPILE_FEATURES cxx_std_14)
+
+	if(_CXX_ATOMIC_BUILTIN)
+		# Nothing to add...
+	elseif(_CXX_ATOMIC_LIB_NEEDED)
+		set_target_properties(std::atomic PROPERTIES IMPORTED_LIBNAME atomic)
+		set(STDCPPATOMIC_LIBRARY atomic)
+	endif()
+endif()
+
+cmake_pop_check_state()
+
+set(Atomic_FOUND ${_found} CACHE BOOL "TRUE if we can run a program using std::atomic" FORCE)
+
+if(Atomic_FIND_REQUIRED AND NOT Atomic_FOUND)
+    message(FATAL_ERROR "Cannot run simple program using std::atomic")
+endif()

cmake/FindFilesystem.cmake

@@ -59,8 +59,10 @@ if(_found)
 		# Nothing to add...
 	elseif(CXX_FILESYSTEM_STDCPPFS_NEEDED)
 		set_target_properties(std::filesystem PROPERTIES IMPORTED_LIBNAME stdc++fs)
+		set(STDCPPFS_LIBRARY stdc++fs)
 	elseif(CXX_FILESYSTEM_CPPFS_NEEDED)
 		set_target_properties(std::filesystem PROPERTIES IMPORTED_LIBNAME c++fs)
+		set(STDCPPFS_LIBRARY c++fs)
 	endif()
 endif()
 

cmake/VersionString.cmake

@@ -0,0 +1,77 @@
+# SPDX-License-Identifier: BSD-2-Clause
+
+# Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+cmake_minimum_required(VERSION 3.15)
+
+# Create a revision file, containing the current git version info, if any
+function(write_version_header)
+	include(GetGitRevisionDescription)
+	if(NOT(GIT-NOTFOUND OR HEAD-HASH-NOTFOUND))
+		git_describe_working_tree(BUILD_VERSION_STRING --match=build --dirty)
+
+		if(BUILD_VERSION_STRING MATCHES "build-([0-9]+)-g([0-9a-f]+)(-dirty)?")
+			set(BUILD_GIT_TAGREF "${CMAKE_MATCH_2}")
+			if(CMAKE_MATCH_3)
+				set(BUILD_VERSION_STRING "${CMAKE_MATCH_1}*")
+			else()
+				set(BUILD_VERSION_STRING "${CMAKE_MATCH_1}")
+			endif()
+		endif()
+	else()
+		set(BUILD_VERSION_STRING "no git info available")
+	endif()
+
+	include_directories(${CMAKE_BINARY_DIR} PRIVATE)
+	string(TIMESTAMP BUILD_DATE_TIME "%Y-%m-%dT%H:%M:%SZ" UTC)
+
+	if(ARGC GREATER 0)
+		set(VAR_PREFIX "${ARGV0}")
+	endif()
+
+	file(WRITE "${CMAKE_BINARY_DIR}/revision.hpp.in" [[// Generated revision file
+
+#pragma once
+
+#include <ostream>
+
+const char k@VAR_PREFIX@ProjectName[] = "@PROJECT_NAME@";
+const char k@VAR_PREFIX@VersionNumber[] = "@PROJECT_VERSION@";
+const char k@VAR_PREFIX@VersionGitTag[] = "@BUILD_GIT_TAGREF@";
+const char k@VAR_PREFIX@BuildInfo[] = "@BUILD_VERSION_STRING@";
+const char k@VAR_PREFIX@BuildDate[] = "@BUILD_DATE_TIME@";
+
+inline void write_version_string(std::ostream &os, bool verbose)
+{
+	os << k@VAR_PREFIX@ProjectName << " version " << k@VAR_PREFIX@VersionNumber << std::endl;
+	if (verbose)
+	{
+		os << "build: " << k@VAR_PREFIX@BuildInfo << ' ' << k@VAR_PREFIX@BuildDate << std::endl;
+		if (k@VAR_PREFIX@VersionGitTag[0] != 0)
+			os << "git tag: " << k@VAR_PREFIX@VersionGitTag << std::endl;
+	}
+}
+]])
+	configure_file("${CMAKE_BINARY_DIR}/revision.hpp.in" "${CMAKE_BINARY_DIR}/revision.hpp" @ONLY)
+endfunction()
+

data/ccd-subset.cif

@@ -2394,3 +2394,189 @@ VAL "Create component"  1999-07-08 RCSB
 VAL "Modify descriptor" 2011-06-04 RCSB 
 # 
 
+data_NAG
+#
+
+_chem_comp.id                                   NAG
+_chem_comp.name                                 2-acetamido-2-deoxy-beta-D-glucopyranose
+_chem_comp.type                                 "D-saccharide, beta linking"
+_chem_comp.pdbx_type                            ATOMS
+_chem_comp.formula                              "C8 H15 N O6"
+_chem_comp.mon_nstd_parent_comp_id              ?
+_chem_comp.pdbx_synonyms
+;N-acetyl-beta-D-glucosamine; 2-acetamido-2-deoxy-beta-D-glucose; 2-acetamido-2-deoxy-D-glucose;
+2-acetamido-2-deoxy-glucose; N-ACETYL-D-GLUCOSAMINE
+;
+
+_chem_comp.pdbx_formal_charge                   0
+_chem_comp.pdbx_initial_date                    1999-07-08
+_chem_comp.pdbx_modified_date                   2020-07-17
+_chem_comp.pdbx_ambiguous_flag                  N
+_chem_comp.pdbx_release_status                  REL
+_chem_comp.pdbx_replaced_by                     ?
+_chem_comp.pdbx_replaces                        ?
+_chem_comp.formula_weight                       221.208
+_chem_comp.one_letter_code                      ?
+_chem_comp.three_letter_code                    NAG
+_chem_comp.pdbx_model_coordinates_details       ?
+_chem_comp.pdbx_model_coordinates_missing_flag  N
+_chem_comp.pdbx_ideal_coordinates_details       Corina
+_chem_comp.pdbx_ideal_coordinates_missing_flag  N
+_chem_comp.pdbx_model_coordinates_db_code       8PCH
+_chem_comp.pdbx_subcomponent_list               ?
+_chem_comp.pdbx_processing_site                 RCSB
+#   #
+loop_
+_pdbx_chem_comp_synonyms.ordinal
+_pdbx_chem_comp_synonyms.comp_id
+_pdbx_chem_comp_synonyms.name
+_pdbx_chem_comp_synonyms.provenance
+_pdbx_chem_comp_synonyms.type
+1  NAG  N-acetyl-beta-D-glucosamine         PDB  ?
+2  NAG  2-acetamido-2-deoxy-beta-D-glucose  PDB  ?
+3  NAG  2-acetamido-2-deoxy-D-glucose       PDB  ?
+4  NAG  2-acetamido-2-deoxy-glucose         PDB  ?
+5  NAG  N-ACETYL-D-GLUCOSAMINE              PDB  ?
+#   #
+loop_
+_chem_comp_atom.comp_id
+_chem_comp_atom.atom_id
+_chem_comp_atom.alt_atom_id
+_chem_comp_atom.type_symbol
+_chem_comp_atom.charge
+_chem_comp_atom.pdbx_align
+_chem_comp_atom.pdbx_aromatic_flag
+_chem_comp_atom.pdbx_leaving_atom_flag
+_chem_comp_atom.pdbx_stereo_config
+_chem_comp_atom.model_Cartn_x
+_chem_comp_atom.model_Cartn_y
+_chem_comp_atom.model_Cartn_z
+_chem_comp_atom.pdbx_model_Cartn_x_ideal
+_chem_comp_atom.pdbx_model_Cartn_y_ideal
+_chem_comp_atom.pdbx_model_Cartn_z_ideal
+_chem_comp_atom.pdbx_component_atom_id
+_chem_comp_atom.pdbx_component_comp_id
+_chem_comp_atom.pdbx_ordinal
+NAG  C1   C1   C  0  1  N  N  R  7.396  28.163  26.662   0.185   1.082  -0.421  C1   NAG   1
+NAG  C2   C2   C  0  1  N  N  R  6.973  29.233  27.644   0.790  -0.220   0.112  C2   NAG   2
+NAG  C3   C3   C  0  1  N  N  R  7.667  29.055  29.000  -0.124  -1.390  -0.265  C3   NAG   3
+NAG  C4   C4   C  0  1  N  N  S  7.573  27.588  29.490  -1.526  -1.129   0.294  C4   NAG   4
+NAG  C5   C5   C  0  1  N  N  R  7.902  26.592  28.373  -2.042   0.207  -0.246  C5   NAG   5
+NAG  C6   C6   C  0  1  N  N  N  7.599  25.173  28.797  -3.417   0.504   0.355  C6   NAG   6
+NAG  C7   C7   C  0  1  N  N  N  6.291  31.299  26.595   3.197   0.157   0.076  C7   NAG   7
+NAG  C8   C8   C  0  1  N  N  N  6.684  32.649  26.036   4.559  -0.052  -0.533  C8   NAG   8
+NAG  N2   N2   N  0  1  N  N  N  7.268  30.545  27.089   2.114  -0.422  -0.480  N2   NAG   9
+NAG  O1   O1   O  0  1  N  Y  N  6.676  28.363  25.419   1.003   2.185  -0.024  O1   NAG  10
+NAG  O3   O3   O  0  1  N  N  N  7.038  29.909  29.947   0.395  -2.600   0.291  O3   NAG  11
+NAG  O4   O4   O  0  1  N  N  N  8.494  27.358  30.574  -2.405  -2.180  -0.114  O4   NAG  12
+NAG  O5   O5   O  0  1  N  N  N  7.104  26.875  27.206  -1.130   1.248   0.113  O5   NAG  13
+NAG  O6   O6   O  0  1  N  N  N  6.232  25.040  29.165  -3.949   1.691  -0.236  O6   NAG  14
+NAG  O7   O7   O  0  1  N  N  N  5.114  30.936  26.562   3.074   0.845   1.067  O7   NAG  15
+NAG  H1   H1   H  0  1  N  N  N  8.477  28.257  26.481   0.133   1.040  -1.509  H1   NAG  16
+NAG  H2   H2   H  0  1  N  N  N  5.888  29.146  27.803   0.879  -0.163   1.197  H2   NAG  17
+NAG  H3   H3   H  0  1  N  N  N  8.729  29.321  28.892  -0.174  -1.478  -1.350  H3   NAG  18
+NAG  H4   H4   H  0  1  N  N  N  6.544  27.403  29.831  -1.483  -1.091   1.382  H4   NAG  19
+NAG  H5   H5   H  0  1  N  N  N  8.971  26.674  28.128  -2.123   0.154  -1.332  H5   NAG  20
+NAG  H61  H61  H  0  1  N  N  N  7.816  24.492  27.961  -4.088  -0.333   0.157  H61  NAG  21
+NAG  H62  H62  H  0  1  N  N  N  8.232  24.910  29.657  -3.320   0.645   1.431  H62  NAG  22
+NAG  H81  H81  H  0  1  N  N  N  5.791  33.159  25.646   4.560   0.320  -1.558  H81  NAG  23
+NAG  H82  H82  H  0  1  N  N  N  7.136  33.258  26.833   5.305   0.490   0.050  H82  NAG  24
+NAG  H83  H83  H  0  1  N  N  N  7.411  32.511  25.222   4.799  -1.115  -0.532  H83  NAG  25
+NAG  HN2  HN2  H  0  1  N  N  N  8.210  30.881  27.079   2.212  -0.973  -1.273  HN2  NAG  26
+NAG  HO1  HO1  H  0  1  N  Y  N  6.933  27.696  24.793   0.679   3.044  -0.328  HO1  NAG  27
+NAG  HO3  HO3  H  0  1  N  Y  N  7.459  29.809  30.793  -0.135  -3.384   0.091  HO3  NAG  28
+NAG  HO4  HO4  H  0  1  N  Y  N  8.425  26.456  30.863  -3.312  -2.079   0.206  HO4  NAG  29
+NAG  HO6  HO6  H  0  1  N  Y  N  6.060  24.143  29.428  -4.822   1.940   0.099  HO6  NAG  30
+#   #
+loop_
+_chem_comp_bond.comp_id
+_chem_comp_bond.atom_id_1
+_chem_comp_bond.atom_id_2
+_chem_comp_bond.value_order
+_chem_comp_bond.pdbx_aromatic_flag
+_chem_comp_bond.pdbx_stereo_config
+_chem_comp_bond.pdbx_ordinal
+NAG  C1  C2   SING  N  N   1
+NAG  C1  O1   SING  N  N   2
+NAG  C1  O5   SING  N  N   3
+NAG  C1  H1   SING  N  N   4
+NAG  C2  C3   SING  N  N   5
+NAG  C2  N2   SING  N  N   6
+NAG  C2  H2   SING  N  N   7
+NAG  C3  C4   SING  N  N   8
+NAG  C3  O3   SING  N  N   9
+NAG  C3  H3   SING  N  N  10
+NAG  C4  C5   SING  N  N  11
+NAG  C4  O4   SING  N  N  12
+NAG  C4  H4   SING  N  N  13
+NAG  C5  C6   SING  N  N  14
+NAG  C5  O5   SING  N  N  15
+NAG  C5  H5   SING  N  N  16
+NAG  C6  O6   SING  N  N  17
+NAG  C6  H61  SING  N  N  18
+NAG  C6  H62  SING  N  N  19
+NAG  C7  C8   SING  N  N  20
+NAG  C7  N2   SING  N  N  21
+NAG  C7  O7   DOUB  N  N  22
+NAG  C8  H81  SING  N  N  23
+NAG  C8  H82  SING  N  N  24
+NAG  C8  H83  SING  N  N  25
+NAG  N2  HN2  SING  N  N  26
+NAG  O1  HO1  SING  N  N  27
+NAG  O3  HO3  SING  N  N  28
+NAG  O4  HO4  SING  N  N  29
+NAG  O6  HO6  SING  N  N  30
+#   #
+loop_
+_pdbx_chem_comp_descriptor.comp_id
+_pdbx_chem_comp_descriptor.type
+_pdbx_chem_comp_descriptor.program
+_pdbx_chem_comp_descriptor.program_version
+_pdbx_chem_comp_descriptor.descriptor
+NAG  SMILES            ACDLabs               12.01  "O=C(NC1C(O)C(O)C(OC1O)CO)C"
+NAG  InChI             InChI                 1.03   "InChI=1S/C8H15NO6/c1-3(11)9-5-7(13)6(12)4(2-10)15-8(5)14/h4-8,10,12-14H,2H2,1H3,(H,9,11)/t4-,5-,6-,7-,8-/m1/s1"
+NAG  InChIKey          InChI                 1.03   OVRNDRQMDRJTHS-FMDGEEDCSA-N
+NAG  SMILES_CANONICAL  CACTVS                3.370  "CC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O"
+NAG  SMILES            CACTVS                3.370  "CC(=O)N[CH]1[CH](O)O[CH](CO)[CH](O)[CH]1O"
+NAG  SMILES_CANONICAL  "OpenEye OEToolkits"  1.7.6  "CC(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O"
+NAG  SMILES            "OpenEye OEToolkits"  1.7.6  "CC(=O)NC1C(C(C(OC1O)CO)O)O"
+#   #
+loop_
+_pdbx_chem_comp_identifier.comp_id
+_pdbx_chem_comp_identifier.type
+_pdbx_chem_comp_identifier.program
+_pdbx_chem_comp_identifier.program_version
+_pdbx_chem_comp_identifier.identifier
+NAG  "SYSTEMATIC NAME"                      ACDLabs               12.01  "2-(acetylamino)-2-deoxy-beta-D-glucopyranose"
+NAG  "SYSTEMATIC NAME"                      "OpenEye OEToolkits"  1.7.6  "N-[(2R,3R,4R,5S,6R)-6-(hydroxymethyl)-2,4,5-tris(oxidanyl)oxan-3-yl]ethanamide"
+NAG  "CONDENSED IUPAC CARBOHYDRATE SYMBOL"  GMML                  1.0    DGlcpNAcb
+NAG  "COMMON NAME"                          GMML                  1.0    N-acetyl-b-D-glucopyranosamine
+NAG  "IUPAC CARBOHYDRATE SYMBOL"            PDB-CARE              1.0    b-D-GlcpNAc
+NAG  "SNFG CARBOHYDRATE SYMBOL"             GMML                  1.0    GlcNAc
+#   #
+loop_
+_pdbx_chem_comp_feature.comp_id
+_pdbx_chem_comp_feature.type
+_pdbx_chem_comp_feature.value
+_pdbx_chem_comp_feature.source
+_pdbx_chem_comp_feature.support
+NAG  "CARBOHYDRATE ISOMER"                  D         PDB  ?
+NAG  "CARBOHYDRATE RING"                    pyranose  PDB  ?
+NAG  "CARBOHYDRATE ANOMER"                  beta      PDB  ?
+NAG  "CARBOHYDRATE PRIMARY CARBONYL GROUP"  aldose    PDB  ?
+#   #
+loop_
+_pdbx_chem_comp_audit.comp_id
+_pdbx_chem_comp_audit.action_type
+_pdbx_chem_comp_audit.date
+_pdbx_chem_comp_audit.processing_site
+NAG  "Create component"          1999-07-08  RCSB
+NAG  "Modify descriptor"         2011-06-04  RCSB
+NAG  "Modify leaving atom flag"  2011-07-01  RCSB
+NAG  "Modify leaving atom flag"  2012-11-26  RCSB
+NAG  "Other modification"        2019-08-12  RCSB
+NAG  "Other modification"        2019-12-19  RCSB
+NAG  "Other modification"        2020-07-03  RCSB
+NAG  "Modify name"               2020-07-17  RCSB
+NAG  "Modify synonyms"           2020-07-17  RCSB
+##

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,171 +29,179 @@
 #pragma once
 
 #include <cstdint>
-#include <string>
 #include <stdexcept>
+#include <string>
 
 namespace mmcif
 {
 
 enum AtomType : uint8_t
 {
-	Nn = 0,		// Unknown
-	
-	H = 1,		// Hydro­gen
-	He = 2,		// He­lium
+	Nn = 0, // Unknown
 
-	Li = 3,		// Lith­ium
-	Be = 4,		// Beryl­lium
-	B = 5,		// Boron
-	C = 6,		// Carbon
-	N = 7,		// Nitro­gen
-	O = 8,		// Oxy­gen
-	F = 9,		// Fluor­ine
-	Ne = 10,	// Neon
+	H = 1,  // Hydro­gen
+	He = 2, // He­lium
 
-	Na = 11,	// So­dium
-	Mg = 12,	// Magne­sium
-	Al = 13,	// Alumin­ium
-	Si = 14,	// Sili­con
-	P = 15,		// Phos­phorus
-	S = 16,		// Sulfur
-	Cl = 17,	// Chlor­ine
-	Ar = 18,	// Argon
+	Li = 3,  // Lith­ium
+	Be = 4,  // Beryl­lium
+	B = 5,   // Boron
+	C = 6,   // Carbon
+	N = 7,   // Nitro­gen
+	O = 8,   // Oxy­gen
+	F = 9,   // Fluor­ine
+	Ne = 10, // Neon
 
-	K = 19,		// Potas­sium
-	Ca = 20,	// Cal­cium
-	Sc = 21,	// Scan­dium
-	Ti = 22,	// Tita­nium
-	V = 23,		// Vana­dium
-	Cr = 24,	// Chrom­ium
-	Mn = 25,	// Manga­nese
-	Fe = 26,	// Iron
-	Co = 27,	// Cobalt
-	Ni = 28,	// Nickel
-	Cu = 29,	// Copper
-	Zn = 30,	// Zinc
-	Ga = 31,	// Gallium
-	Ge = 32,	// Germa­nium
-	As = 33,	// Arsenic
-	Se = 34,	// Sele­nium
-	Br = 35,	// Bromine
-	Kr = 36,	// Kryp­ton
+	Na = 11, // So­dium
+	Mg = 12, // Magne­sium
+	Al = 13, // Alumin­ium
+	Si = 14, // Sili­con
+	P = 15,  // Phos­phorus
+	S = 16,  // Sulfur
+	Cl = 17, // Chlor­ine
+	Ar = 18, // Argon
 
-	Rb = 37,	// Rubid­ium
-	Sr = 38,	// Stront­ium
-	Y = 39,		// Yttrium
-	Zr = 40,	// Zirco­nium
-	Nb = 41,	// Nio­bium
-	Mo = 42,	// Molyb­denum
-	Tc = 43,	// Tech­netium
-	Ru = 44,	// Ruthe­nium
-	Rh = 45,	// Rho­dium
-	Pd = 46,	// Pallad­ium
-	Ag = 47,	// Silver
-	Cd = 48,	// Cad­mium
-	In = 49,	// Indium
-	Sn = 50,	// Tin
-	Sb = 51,	// Anti­mony
-	Te = 52,	// Tellurium
-	I = 53,		// Iodine
-	Xe = 54,	// Xenon
-	Cs = 55,	// Cae­sium
-	Ba = 56,	// Ba­rium
-	La = 57,	// Lan­thanum
+	K = 19,  // Potas­sium
+	Ca = 20, // Cal­cium
+	Sc = 21, // Scan­dium
+	Ti = 22, // Tita­nium
+	V = 23,  // Vana­dium
+	Cr = 24, // Chrom­ium
+	Mn = 25, // Manga­nese
+	Fe = 26, // Iron
+	Co = 27, // Cobalt
+	Ni = 28, // Nickel
+	Cu = 29, // Copper
+	Zn = 30, // Zinc
+	Ga = 31, // Gallium
+	Ge = 32, // Germa­nium
+	As = 33, // Arsenic
+	Se = 34, // Sele­nium
+	Br = 35, // Bromine
+	Kr = 36, // Kryp­ton
 
-	Hf = 72,	// Haf­nium
-	Ta = 73,	// Tanta­lum
-	W = 74,		// Tung­sten
-	Re = 75,	// Rhe­nium
-	Os = 76,	// Os­mium
-	Ir = 77,	// Iridium
-	Pt = 78,	// Plat­inum
-	Au = 79,	// Gold
-	Hg = 80,	// Mer­cury
-	Tl = 81,	// Thallium
-	Pb = 82,	// Lead
-	Bi = 83,	// Bis­muth
-	Po = 84,	// Polo­nium
-	At = 85,	// Asta­tine
-	Rn = 86,	// Radon
-	Fr = 87,	// Fran­cium
-	Ra = 88,	// Ra­dium
-	Ac = 89,	// Actin­ium
+	Rb = 37, // Rubid­ium
+	Sr = 38, // Stront­ium
+	Y = 39,  // Yttrium
+	Zr = 40, // Zirco­nium
+	Nb = 41, // Nio­bium
+	Mo = 42, // Molyb­denum
+	Tc = 43, // Tech­netium
+	Ru = 44, // Ruthe­nium
+	Rh = 45, // Rho­dium
+	Pd = 46, // Pallad­ium
+	Ag = 47, // Silver
+	Cd = 48, // Cad­mium
+	In = 49, // Indium
+	Sn = 50, // Tin
+	Sb = 51, // Anti­mony
+	Te = 52, // Tellurium
+	I = 53,  // Iodine
+	Xe = 54, // Xenon
+	Cs = 55, // Cae­sium
+	Ba = 56, // Ba­rium
+	La = 57, // Lan­thanum
 
-	Rf = 104,	// Ruther­fordium
-	Db = 105,	// Dub­nium
-	Sg = 106,	// Sea­borgium
-	Bh = 107,	// Bohr­ium
-	Hs = 108,	// Has­sium
-	Mt = 109,	// Meit­nerium
-	Ds = 110,	// Darm­stadtium
-	Rg = 111,	// Roent­genium
-	Cn = 112,	// Coper­nicium
-	Nh = 113,	// Nihon­ium
-	Fl = 114,	// Flerov­ium
-	Mc = 115,	// Moscov­ium
-	Lv = 116,	// Liver­morium
-	Ts = 117,	// Tenness­ine
-	Og = 118,	// Oga­nesson
+	Hf = 72, // Haf­nium
+	Ta = 73, // Tanta­lum
+	W = 74,  // Tung­sten
+	Re = 75, // Rhe­nium
+	Os = 76, // Os­mium
+	Ir = 77, // Iridium
+	Pt = 78, // Plat­inum
+	Au = 79, // Gold
+	Hg = 80, // Mer­cury
+	Tl = 81, // Thallium
+	Pb = 82, // Lead
+	Bi = 83, // Bis­muth
+	Po = 84, // Polo­nium
+	At = 85, // Asta­tine
+	Rn = 86, // Radon
+	Fr = 87, // Fran­cium
+	Ra = 88, // Ra­dium
+	Ac = 89, // Actin­ium
 
-	Ce = 58,	// Cerium
-	Pr = 59,	// Praseo­dymium
-	Nd = 60,	// Neo­dymium
-	Pm = 61,	// Prome­thium
-	Sm = 62,	// Sama­rium
-	Eu = 63,	// Europ­ium
-	Gd = 64,	// Gadolin­ium
-	Tb = 65,	// Ter­bium
-	Dy = 66,	// Dyspro­sium
-	Ho = 67,	// Hol­mium
-	Er = 68,	// Erbium
-	Tm = 69,	// Thulium
-	Yb = 70,	// Ytter­bium
-	Lu = 71,	// Lute­tium
+	Rf = 104, // Ruther­fordium
+	Db = 105, // Dub­nium
+	Sg = 106, // Sea­borgium
+	Bh = 107, // Bohr­ium
+	Hs = 108, // Has­sium
+	Mt = 109, // Meit­nerium
+	Ds = 110, // Darm­stadtium
+	Rg = 111, // Roent­genium
+	Cn = 112, // Coper­nicium
+	Nh = 113, // Nihon­ium
+	Fl = 114, // Flerov­ium
+	Mc = 115, // Moscov­ium
+	Lv = 116, // Liver­morium
+	Ts = 117, // Tenness­ine
+	Og = 118, // Oga­nesson
 
-	Th = 90,	// Thor­ium
-	Pa = 91,	// Protac­tinium
-	U = 92,		// Ura­nium
-	Np = 93,	// Neptu­nium
-	Pu = 94,	// Pluto­nium
-	Am = 95,	// Ameri­cium
-	Cm = 96,	// Curium
-	Bk = 97,	// Berkel­ium
-	Cf = 98,	// Califor­nium
-	Es = 99,	// Einstei­nium
-	Fm = 100,	// Fer­mium
-	Md = 101,	// Mende­levium
-	No = 102,	// Nobel­ium
-	Lr = 103,	// Lawren­cium
+	Ce = 58, // Cerium
+	Pr = 59, // Praseo­dymium
+	Nd = 60, // Neo­dymium
+	Pm = 61, // Prome­thium
+	Sm = 62, // Sama­rium
+	Eu = 63, // Europ­ium
+	Gd = 64, // Gadolin­ium
+	Tb = 65, // Ter­bium
+	Dy = 66, // Dyspro­sium
+	Ho = 67, // Hol­mium
+	Er = 68, // Erbium
+	Tm = 69, // Thulium
+	Yb = 70, // Ytter­bium
+	Lu = 71, // Lute­tium
 
-	D = 129,	// Deuterium
+	Th = 90,  // Thor­ium
+	Pa = 91,  // Protac­tinium
+	U = 92,   // Ura­nium
+	Np = 93,  // Neptu­nium
+	Pu = 94,  // Pluto­nium
+	Am = 95,  // Ameri­cium
+	Cm = 96,  // Curium
+	Bk = 97,  // Berkel­ium
+	Cf = 98,  // Califor­nium
+	Es = 99,  // Einstei­nium
+	Fm = 100, // Fer­mium
+	Md = 101, // Mende­levium
+	No = 102, // Nobel­ium
+	Lr = 103, // Lawren­cium
+
+	D = 129, // Deuterium
 };
 
 // --------------------------------------------------------------------
 // AtomTypeInfo
 
-enum RadiusType {
-	eRadiusCalculated,
-	eRadiusEmpirical,
-	eRadiusCovalentEmpirical,
+enum class RadiusType
+{
+	Calculated,
+	Empirical,
+	CovalentEmpirical,
 
-	eRadiusSingleBond,
-	eRadiusDoubleBond,
-	eRadiusTripleBond,
+	SingleBond,
+	DoubleBond,
+	TripleBond,
 
-	eRadiusVanderWaals,
+	VanderWaals,
 
-	eRadiusTypeCount
+	TypeCount
+};
+
+constexpr size_t RadiusTypeCount = static_cast<size_t>(RadiusType::TypeCount);
+
+enum class IonicRadiusType
+{
+	Effective, Crystal
 };
 
 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[RadiusTypeCount];
 };
 
 extern const AtomTypeInfo kKnownAtoms[];
@@ -205,25 +213,46 @@ 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);
-	
-	float radius(RadiusType type = eRadiusSingleBond) const
+	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 = RadiusType::SingleBond) const
 	{
-		if (type >= eRadiusTypeCount)
+		if (type >= RadiusType::TypeCount)
 			throw std::invalid_argument("invalid radius requested");
-		return mInfo->radii[type] / 100.f;
+		return mInfo->radii[static_cast<size_t>(type)] / 100.f;
 	}
-	
+
+	/// \brief Return the radius for a charged version of this atom in a solid crystal
+	///
+	/// \param charge  The charge of the ion
+	/// \return        The radius of the ion
+	float crystal_ionic_radius(int charge) const;
+
+	/// \brief Return the radius for a charged version of this atom in a non-solid environment
+	///
+	/// \param charge  The charge of the ion
+	/// \return        The radius of the ion
+	float effective_ionic_radius(int charge) const;
+
+	/// \brief Return the radius for a charged version of this atom, returns the effective radius by default
+	///
+	/// \param charge  The charge of the ion
+	/// \return        The radius of the ion
+	float ionic_radius(int charge, IonicRadiusType type = IonicRadiusType::Effective) const
+	{
+		return type == IonicRadiusType::Effective ? effective_ionic_radius(charge) : crystal_ionic_radius(charge);
+	}
+
 	// 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 +260,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

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

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)
 	{
@@ -159,6 +176,7 @@ class Item
 		, mValue(rhs.mValue)
 	{
 	}
+
 	Item(Item &&rhs) noexcept
 		: mName(std::move(rhs.mName))
 		, mValue(std::move(rhs.mValue))
@@ -221,7 +239,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 +248,38 @@ 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);
+	/// \brief Access to the category with name \a name, will create it if it doesn't exist.
+	Category &operator[](std::string_view name);
 
-	std::tuple<iterator, bool> emplace(const std::string &name);
+	// /// \brief Access to the category with name \a name, will throw if it doesn't exist.
+	// const Category &operator[](std::string_view name) const;
+
+	/// \brief Access to the category with name \a name, will return an empty category if is doesn't exist.
+	const Category &operator[](std::string_view name) const;
+
+	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,10 +287,14 @@ 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;
+
+	// for returning empty categories in the const operator[]
+	mutable std::unique_ptr<Category> mNullCategory;
 };
 
 // --------------------------------------------------------------------
@@ -350,14 +377,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 +392,7 @@ namespace detail
 		{
 		}
 
-		const char *mName;
+		std::string_view mName;
 		size_t mColumn;
 		Row &mRow;
 		bool mConst = false;
@@ -746,16 +773,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 +803,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 +825,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;
@@ -1167,13 +1194,16 @@ struct Empty
 {
 };
 
+/// \brief A helper to make it possible to have conditions like ("id"_key == cif::null)
+inline constexpr Empty null = 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 +1846,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>;
@@ -1843,6 +1873,9 @@ class Category
 	Row front() { return Row(mHead); }
 	Row back() { return Row(mTail); }
 
+	const Row front() const { return Row(mHead); }
+	const Row back() const { return Row(mTail); }
+
 	Row operator[](Condition &&cond);
 	const Row operator[](Condition &&cond) const;
 
@@ -2064,7 +2097,7 @@ class Category
 
 	Datablock &db() { return mDb; }
 
-	void setValidator(Validator *v);
+	void setValidator(const Validator *v);
 
 	iset fields() const;
 	iset mandatoryFields() const;
@@ -2077,8 +2110,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 +2152,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;
 };
 
 // --------------------------------------------------------------------
@@ -2142,17 +2186,21 @@ class File
 	File();
 	File(std::istream &is, bool validate = false);
 	File(const std::filesystem::path &path, bool validate = false);
+	File(const char *data, size_t length); // good luck trying to find out what it is...
+
 	File(File &&rhs);
 	File(const File &rhs) = delete;
 	File &operator=(const File &rhs) = delete;
 
-	~File();
+	virtual ~File();
 
-	void load(const std::filesystem::path &p);
-	void save(const std::filesystem::path &p);
+	virtual void load(const std::filesystem::path &p);
+	virtual void save(const std::filesystem::path &p);
 
-	void load(std::istream &is);
-	void save(std::ostream &os);
+	virtual void load(std::istream &is);
+	virtual void save(std::ostream &os);
+
+	virtual void load(const char *data, std::size_t length);
 
 	/// \brief Load only the data block \a datablock from the mmCIF file
 	void load(std::istream &is, const std::string &datablock);
@@ -2162,7 +2210,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 +2232,14 @@ class File
 
 	void append(Datablock *e);
 
-	Datablock *get(const std::string &name) const;
-	Datablock &operator[](const std::string &name);
+	Datablock &emplace(std::string_view name)
+	{
+		append(new Datablock(name));
+		return back();
+	}
+	
+	Datablock *get(std::string_view name) const;
+	Datablock &operator[](std::string_view name);
 
 	struct iterator
 	{
@@ -2220,16 +2275,17 @@ class File
 	iterator begin() const;
 	iterator end() const;
 
+	Datablock &front();
+	Datablock &back();
+
 	bool empty() const { return mHead == nullptr; }
 
 	const Validator &getValidator() const;
 	void getTagOrder(std::vector<std::string> &tags) const;
 
   private:
-	void setValidator(Validator *v);
-
 	Datablock *mHead;
-	Validator *mValidator;
+	const Validator *mValidator;
 };
 
 // --------------------------------------------------------------------

include/cif++/Cif2PDB.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
@@ -28,12 +28,12 @@
 
 #include "cif++/Cif++.hpp"
 
-void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile);
+void WritePDBFile(std::ostream &os, const cif::Datablock &data);
 
 /// \brief Just the HEADER, COMPND, SOURCE and AUTHOR lines
-void WritePDBHeaderLines(std::ostream& os, cif::File& cifFile);
+void WritePDBHeaderLines(std::ostream &os, const cif::Datablock &data);
 
-std::string GetPDBHEADERLine(cif::File& cifFile, std::string::size_type truncate_at = 127);
-std::string GetPDBCOMPNDLine(cif::File& cifFile, std::string::size_type truncate_at = 127);
-std::string GetPDBSOURCELine(cif::File& cifFile, std::string::size_type truncate_at = 127);
-std::string GetPDBAUTHORLine(cif::File& cifFile, std::string::size_type truncate_at = 127);
+std::string GetPDBHEADERLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
+std::string GetPDBCOMPNDLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
+std::string GetPDBSOURCELine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
+std::string GetPDBAUTHORLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);

include/cif++/CifParser.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
@@ -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,30 +76,17 @@ 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)
-{
-	bool result = isOrdinary(*s++);
-	while (result and *s != 0)
-	{
-		result = isNonBlank(*s);
-		++s;
-	}
-	return result;
-}
+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 +94,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 +112,7 @@ class SacParser
 		eCIFTokenValue,
 	};
 
-	static const char* kTokenName[];
+	static const char *kTokenName[];
 
 	enum CIFValueType
 	{
@@ -137,40 +125,39 @@ class SacParser
 		eCIFValueUnknown
 	};
 
-	static const char* kValueName[];
-	
+	static const char *kValueName[];
+
 	int getNextChar();
 
 	void retract();
-	void restart();
-	
+	int restart(int start);
+
 	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 +172,21 @@ class SacParser
 		eStateTextField,
 		eStateFloat = 100,
 		eStateInt = 110,
-//		eStateNumericSuffix = 200,
-		eStateValue = 300
+		eStateValue = 300,
+		eStateDATA,
+		eStateSAVE
 	};
 
-	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;
+	CIFToken mLookahead;
+	std::string mTokenValue;
+	CIFValueType mTokenType;
+	std::stack<int> mBuffer;
 };
 
 // --------------------------------------------------------------------
@@ -207,18 +194,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 +213,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

include/cif++/CifUtils.hpp

@@ -26,11 +26,7 @@
 
 #pragma once
 
-#include <cassert>
 #include <filesystem>
-#include <iostream>
-#include <list>
-#include <memory>
 #include <set>
 #include <vector>
 
@@ -67,8 +63,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 +96,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

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

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;

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

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
 	{
@@ -204,52 +221,52 @@ struct PointF
 
 	FType length() const
 	{
-		return sqrt(mX * mX + mY * mY + mZ * mZ);
+		return std::sqrt(mX * mX + mY * mY + mZ * mZ);
 	}
 };
 
 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,97 +274,96 @@ 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(
+	return std::sqrt(
 		(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 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)
 	{
-		u = DotProduct(p, x) / sqrt(u);
-		v = DotProduct(p, y) / sqrt(v);
+		u = DotProduct(p, x) / std::sqrt(u);
+		v = DotProduct(p, y) / std::sqrt(v);
 		if (u != 0 or v != 0)
-			result = atan2(v, u) * 180 / kPI;
+			result = std::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);
-	
+		result = DotProduct(v12, v34) / std::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,74 +371,86 @@ 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>
-PointF<F> Nudge(PointF<F> p, F offset);
+Point Nudge(Point p, float offset);
 
 // --------------------------------------------------------------------
 // We use quaternions to do rotations in 3d space
 
 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);
+Quaternion ConstructFromAngleAxis(float angle, Point axis);
+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->mX = std::sin(lon) * std::cos(lat);
+			p->mY = std::cos(lon) * std::cos(lat);
+			p->mZ = std::sin(lat);
 
 			++p;
 		}
 	}
 
   private:
-
-	array_type				mPoints;
-	double					mWeight;
+	array_type mPoints;
+	double mWeight;
 };
 
 typedef SphericalDots<50> SphericalDots_50;
 
-}
+} // namespace mmcif

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

include/cif++/Structure.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
@@ -34,16 +34,16 @@
 #include "cif++/Point.hpp"
 
 /*
-	To modify a structure, you will have to use actions.
-	
-	The currently supported actions are:
-	
+    To modify a structure, you will have to use actions.
+
+    The currently supported actions are:
+
 //	- Move atom to new location
-	- Remove atom
+    - Remove atom
 //	- Add new atom that was formerly missing
 //	- Add alternate Residue
-	- 
-	
+    -
+
 */
 
 namespace mmcif
@@ -60,30 +60,112 @@ 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;
+
+		int charge() 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)
+	{
+		if (not mImpl)
+			throw std::logic_error("Error trying to modify an uninitialized atom");
+		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 impl().mID; }
+	AtomType type() const { return impl().mType; }
+
+	Point location() const { return impl().mLocation; }
+	void location(Point p)
+	{
+		if (not mImpl)
+			throw std::logic_error("Error trying to modify an uninitialized atom");
+		mImpl->moveTo(p);
+	}
 
 	/// \brief Translate the position of this atom by \a t
 	void translate(Point t);
@@ -91,47 +173,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 impl().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 impl().mSymmetryCopy; }
+	std::string symmetry() const { return impl().mSymmetryOperator; }
 
-	const Compound &comp() const;
-	bool isWater() const;
+	const Compound &comp() const { return impl().comp(); }
+	bool isWater() const { return impl().mCompID == "HOH" or impl().mCompID == "H2O" or impl().mCompID == "WAT"; }
 	int charge() const;
 
 	float uIso() const;
-	bool getAnisoU(float anisou[6]) const;
+	bool getAnisoU(float anisou[6]) const { return impl().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 impl().mAtomID; }
+	const std::string &labelCompID() const { return impl().mCompID; }
+	const std::string &labelAsymID() const { return impl().mAsymID; }
 	std::string labelEntityID() const;
-	int labelSeqID() const;
-	std::string labelAltID() const;
-	bool isAlternate() const;
+	int labelSeqID() const { return impl().mSeqID; }
+	const std::string &labelAltID() const { return impl().mAltID; }
+	bool isAlternate() const { return not impl().mAltID.empty(); }
 
 	std::string authAtomID() const;
 	std::string authCompID() const;
 	std::string authAsymID() const;
-	std::string authSeqID() const;
+	const std::string &authSeqID() const { return impl().mAuthSeqID; }
 	std::string pdbxAuthInsCode() const;
 	std::string pdbxAuthAltID() const;
 
@@ -139,13 +214,10 @@ class Atom
 	std::string pdbID() const;   // auth_comp_id + '_' + auth_asym_id + '_' + auth_seq_id + pdbx_PDB_ins_code
 
 	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;
+	bool operator!=(const Atom &rhs) const
+	{
+		return not operator==(rhs);
+	}
 
 	// access data in compound for this atom
 
@@ -158,10 +230,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 impl().compare(*b.mImpl); }
 
 	bool operator<(const Atom &rhs) const
 	{
@@ -170,16 +242,45 @@ class Atom
 
 	friend std::ostream &operator<<(std::ostream &os, const Atom &atom);
 
+	/// \brief Synchronize data with underlying cif data
+	void sync()
+	{
+		if (mImpl)
+			mImpl->prefetch();
+	}
+
   private:
 	friend class Structure;
-	void setID(int id);
 
-	AtomImpl *impl();
-	const AtomImpl *impl() const;
+	const AtomImpl &impl() const
+	{
+		if (not mImpl)
+			throw std::runtime_error("Uninitialized atom, not found?");
+		return *mImpl;
+	}
 
-	struct AtomImpl *mImpl_;
+	std::shared_ptr<AtomImpl> mImpl;
 };
 
+template <>
+inline std::string Atom::get_property<std::string>(const std::string_view name) const
+{
+	return impl().get_property(name);
+}
+
+template <>
+inline int Atom::get_property<int>(const std::string_view name) const
+{
+	auto v = impl().get_property(name);
+	return v.empty() ? 0 : stoi(v);
+}
+
+template <>
+inline float Atom::get_property<float>(const std::string_view name) const
+{
+	return stof(impl().get_property(name));
+}
+
 inline void swap(mmcif::Atom &a, mmcif::Atom &b)
 {
 	a.swap(b);
@@ -199,22 +300,26 @@ typedef std::vector<Atom> AtomView;
 
 // --------------------------------------------------------------------
 
+enum class EntityType
+{
+	Polymer, NonPolymer, Macrolide, Water, Branched
+};
+
+// --------------------------------------------------------------------
+
 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, 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;
@@ -225,8 +330,12 @@ class Residue
 	virtual ~Residue();
 
 	const Compound &compound() const;
+
+	AtomView &atoms();
 	const AtomView &atoms() const;
 
+	void addAtom(Atom &atom);
+
 	/// \brief Unique atoms returns only the atoms without alternates and the first of each alternate atom id.
 	AtomView unique_atoms() const;
 
@@ -236,10 +345,14 @@ 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;
 
+	EntityType entityType() const;
+
 	std::string authAsymID() const;
 	std::string authSeqID() const;
 	std::string authInsCode() const;
@@ -275,6 +388,18 @@ class Residue
 
 	friend std::ostream &operator<<(std::ostream &os, const Residue &res);
 
+	friend Structure;
+
+	bool operator==(const mmcif::Residue &rhs) const
+	{
+		return this == &rhs or (
+			mStructure == rhs.mStructure and
+			mSeqID == rhs.mSeqID and
+			mAsymID == rhs.mAsymID and
+			mCompoundID == rhs.mCompoundID and
+			mAuthSeqID == rhs.mAuthSeqID);
+	}
+
   protected:
 	Residue() {}
 
@@ -283,9 +408,6 @@ class Residue
 	const Structure *mStructure = nullptr;
 	std::string mCompoundID, mAsymID;
 	int mSeqID = 0;
-
-	// Watch out, this is used only to label waters... The rest of the code relies on
-	// MapLabelToAuth to get this info. Perhaps we should rename this member field.
 	std::string mAuthSeqID;
 	AtomView mAtoms;
 };
@@ -352,6 +474,11 @@ class Monomer : public Residue
 	// for LEU and VAL
 	float chiralVolume() const;
 
+	bool operator==(const Monomer &rhs) const
+	{
+		return mPolymer == rhs.mPolymer and mIndex == rhs.mIndex;
+	}
+
   private:
 	const Polymer *mPolymer;
 	size_t mIndex;
@@ -389,35 +516,87 @@ class Polymer : public std::vector<Monomer>
 	cif::RowSet mPolySeq;
 };
 
+// --------------------------------------------------------------------
+// Sugar and Branch, to describe glycosylation sites
+
+class Branch;
+
+class Sugar : public Residue
+{
+  public:
+	Sugar(const Branch &branch, const std::string &compoundID,
+		const std::string &asymID, int authSeqID);
+
+	int num() const { return std::stoi(mAuthSeqID); }
+	std::string name() const;
+
+	/// \brief Return the atom the C1 is linked to
+	Atom getLink() const { return mLink; }
+	void setLink(Atom link) { mLink = link; }
+
+  private:
+
+	const Branch &mBranch;
+	Atom mLink;
+};
+
+class Branch : public std::vector<Sugar>
+{
+  public:
+	Branch(Structure &structure, const std::string &asymID);
+
+	void linkAtoms();
+
+	std::string name() const;
+	float weight() const;
+	std::string asymID() const { return mAsymID; }
+
+	Structure &structure() { return *mStructure; }
+	const Structure &structure() const { return *mStructure; }
+
+	Sugar &getSugarByNum(int nr);
+	const Sugar &getSugarByNum(int nr) const;
+
+  private:
+	friend Sugar;
+
+	std::string name(const Sugar &s) const;
+
+	Structure *mStructure;
+	std::string mAsymID;
+};
+
 // --------------------------------------------------------------------
 // file is a reference to the data stored in e.g. the cif file.
 // This object is not copyable.
 
-class File : public std::enable_shared_from_this<File>
+class File : public cif::File
 {
   public:
-	File();
-	File(const std::filesystem::path &path);
-	File(const char *data, size_t length); // good luck trying to find out what it is...
-	~File();
+	File() {}
+
+	File(const std::filesystem::path &path)
+	{
+		load(path);
+	}
+
+	File(const char *data, size_t length)
+	{
+		load(data, length);
+	}
 
 	File(const File &) = delete;
 	File &operator=(const File &) = delete;
 
-	cif::Datablock& createDatablock(const std::string &name);
+	void load(const std::filesystem::path &p) override;
+	void save(const std::filesystem::path &p) override;
 
-	void load(const std::filesystem::path &path);
-	void save(const std::filesystem::path &path);
+	void load(std::istream &is) override;
 
-	Structure *model(size_t nr = 1);
+	using cif::File::load;
+	using cif::File::save;
 
-	struct FileImpl &impl() const { return *mImpl; }
-
-	cif::Datablock &data();
-	cif::File &file();
-
-  private:
-	struct FileImpl *mImpl;
+	cif::Datablock &data() { return front(); }
 };
 
 // --------------------------------------------------------------------
@@ -437,60 +616,115 @@ inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
 class Structure
 {
   public:
-	Structure(File &p, size_t modelNr = 1, StructureOpenOptions options = {});
-	Structure &operator=(const Structure &) = delete;
-	~Structure();
+	Structure(cif::File &p, size_t modelNr = 1, StructureOpenOptions options = {})
+		: Structure(p.front(), modelNr, options)
+	{
+	}
+
+	Structure(cif::Datablock &db, size_t modelNr = 1, StructureOpenOptions options = {});
 
 	// Create a read-only clone of the current structure (for multithreaded calculations that move atoms)
 	Structure(const Structure &);
 
-	File &getFile() const;
+	Structure &operator=(const Structure &) = delete;
+	~Structure();
 
 	const AtomView &atoms() const { return mAtoms; }
+	// AtomView &atoms() { return mAtoms; }
+
+	EntityType getEntityTypeForEntityID(const std::string entityID) const;
+	EntityType getEntityTypeForAsymID(const std::string asymID) const;
+
 	AtomView waters() const;
 
 	const std::list<Polymer> &polymers() const { return mPolymers; }
 	std::list<Polymer> &polymers() { return mPolymers; }
 
-	const std::vector<Residue> &nonPolymers() const { return mNonPolymers; }
-	const std::vector<Residue> &branchResidues() const { return mBranchResidues; }
+	Polymer &getPolymerByAsymID(const std::string &asymID);
 
-	Atom getAtomByID(std::string id) const;
+	const Polymer &getPolymerByAsymID(const std::string &asymID) const
+	{
+		return const_cast<Structure *>(this)->getPolymerByAsymID(asymID);
+	}
+
+	const std::list<Branch> &branches() const { return mBranches; }
+	std::list<Branch> &branches() { return mBranches; }
+
+	Branch &getBranchByAsymID(const std::string &asymID);
+	const Branch &getBranchByAsymID(const std::string &asymID) const;
+
+	const std::vector<Residue> &nonPolymers() const { return mNonPolymers; }
+	
+	Atom getAtomByID(const std::string &id) const;
 	// Atom getAtomByLocation(Point pt, float maxDistance) const;
 
 	Atom getAtomByLabel(const std::string &atomID, const std::string &asymID,
 		const std::string &compID, int seqID, const std::string &altID = "");
 
-	/// \brief 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 Return the atom closest to point \a p
+	Atom getAtomByPosition(Point p) const;
 
-	// map between auth and label locations
+	/// \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;
 
-	std::tuple<std::string, int, std::string> MapAuthToLabel(const std::string &asymID,
-		const std::string &seqID, const std::string &compID, const std::string &insCode = "");
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	Residue &getResidue(const std::string &asymID)
+	{
+		return getResidue(asymID, 0, "");
+	}
 
-	std::tuple<std::string, std::string, std::string, std::string> MapLabelToAuth(
-		const std::string &asymID, int seqID, const std::string &compID);
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	const Residue &getResidue(const std::string &asymID) const
+	{
+		return getResidue(asymID, 0, "");
+	}
 
-	// returns chain, seqnr, icode
-	std::tuple<char, int, char> MapLabelToAuth(
-		const std::string &asymID, int seqID) const;
+	/// \brief Get a residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
+	Residue &getResidue(const std::string &asymID, int seqID, const std::string &authSeqID);
 
-	// returns chain,seqnr,comp,iCode
-	std::tuple<std::string, int, std::string, std::string> MapLabelToPDB(
-		const std::string &asymID, int seqID, const std::string &compID,
-		const std::string &authSeqID) const;
+	/// \brief Get a the single residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
+	const Residue &getResidue(const std::string &asymID, int seqID, const std::string &authSeqID) const
+	{
+		return const_cast<Structure *>(this)->getResidue(asymID, seqID, authSeqID);
+	}
 
-	std::tuple<std::string, int, std::string> MapPDBToLabel(
-		const std::string &asymID, int seqID, const std::string &compID, const std::string &iCode) const;
+	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
+	Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID);
+
+	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
+	const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID) const
+	{
+		return const_cast<Structure *>(this)->getResidue(asymID, compID, seqID, authSeqID);
+	}
+
+	/// \brief Get a the residue for atom \a atom
+	Residue &getResidue(const mmcif::Atom &atom)
+	{
+		return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID(), atom.authSeqID());
+	}
+
+	/// \brief Get a the residue for atom \a atom
+	const Residue &getResidue(const mmcif::Atom &atom) const
+	{
+		return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID(), atom.authSeqID());
+	}
 
 	// Actions
 	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 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(Residue &res, const std::string &newCompound,
 		const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);
 
+	/// \brief Remove a residue, can be monomer or nonpoly
+	///
+	/// \param asym_id     The asym ID
+	/// \param seq_id      The sequence ID
+	void removeResidue(const std::string &asym_id, int seq_id, const std::string &auth_seq_id)
+	{
+		removeResidue(getResidue(asym_id, seq_id, auth_seq_id));
+	}
+
 	/// \brief Create a new non-polymer entity, returns new ID
 	/// \param mon_id	The mon_id for the new nonpoly, must be an existing and known compound from CCD
 	/// \return			The ID of the created entity
@@ -504,9 +738,33 @@ class Structure
 	/// \return				The newly create asym ID
 	std::string createNonpoly(const std::string &entity_id, const std::vector<mmcif::Atom> &atoms);
 
-	/// \brief To sort the atoms in order of model > asym-id > res-id > atom-id
-	/// Will asssign new atom_id's to all atoms. Be carefull
-	void sortAtoms();
+	/// \brief Create a new NonPolymer struct_asym with atoms constructed from info in \a atom_info, returns asym_id.
+	/// This method creates new atom records filled with info from the info.
+	///
+	/// \param entity_id	The entity ID of the new nonpoly
+	/// \param atoms		The array of sets of cif::item data containing the data for the atoms.
+	/// \return				The newly create asym ID
+	std::string createNonpoly(const std::string &entity_id, std::vector<std::vector<cif::Item>> &atom_info);
+
+	/// \brief Create a new (sugar) branch with one first NAG containing atoms constructed from \a nag_atom_info
+	Branch &createBranch(std::vector<std::vector<cif::Item>> &nag_atom_info);
+
+	/// \brief Extend an existing (sugar) branch identified by \a asymID with one sugar containing atoms constructed from \a atom_info
+	///
+	/// \param asym_id      The asym id of the branch to extend
+	/// \param atom_info    Array containing the info for the atoms to construct for the new sugar
+	/// \param link_sugar   The sugar to link to, note: this is the sugar number (1 based)
+	/// \param link_atom    The atom id of the atom linked in the sugar
+	Branch &extendBranch(const std::string &asym_id, std::vector<std::vector<cif::Item>> &atom_info,
+		int link_sugar, const std::string &link_atom);
+
+	/// \brief Remove \a branch
+	void removeBranch(Branch &branch);
+
+	/// \brief Remove residue \a res
+	///
+	/// \param res         The residue to remove
+	void removeResidue(mmcif::Residue &res);
 
 	/// \brief Translate the coordinates of all atoms in the structure by \a t
 	void translate(Point t);
@@ -514,31 +772,56 @@ 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
+	{
+		return mDb[name];
+	}
+
+	cif::Datablock &datablock() const
+	{
+		return mDb;
+	}
+
+	void validateAtoms() 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();
-	void updateAtomIndex();
+	std::string createEntityForBranch(Branch &branch);
 
-	File &mFile;
+	void loadData();
+
+	void loadAtomsForModel(StructureOpenOptions options);
+
+	template<typename... Args>
+	Atom& emplace_atom(Args ...args)
+	{
+		return emplace_atom(Atom{std::forward<Args>(args)...});
+	}
+
+	Atom &emplace_atom(Atom &&atom);
+
+	cif::Datablock &mDb;
 	size_t mModelNr;
 	AtomView mAtoms;
 	std::vector<size_t> mAtomIndex;
 	std::list<Polymer> mPolymers;
-	std::vector<Residue> mNonPolymers, mBranchResidues;
+	std::list<Branch> mBranches;
+	std::vector<Residue> mNonPolymers;
 };
 
 } // namespace mmcif

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
 
 }

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.datablock();
 
 	// 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,84 @@ 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>, std::string> atomMapByAsymSeqAndAtom;
+	for (auto &a : p.atoms())
 	{
-		auto key = make_tuple(a.labelAsymID(), a.labelSeqID(), a.labelAtomID());
+		auto key = make_tuple(a.labelAsymID(), a.labelSeqID(), a.labelAtomID(), a.authSeqID());
 		atomMapByAsymSeqAndAtom[key] = a.id();
 	}
-	
+
 	// first link all residues in a polyseq
-	
-	std::string lastAsymID;
+
+	std::string lastAsymID, lastAuthSeqID;
 	int lastSeqID = 0;
-	for (auto r: db["pdbx_poly_seq_scheme"])
+	for (const auto &[asymID, seqID, authSeqID] : db["pdbx_poly_seq_scheme"].rows<std::string,int,std::string>("asym_id", "seq_id", "pdb_seq_num"))
 	{
-		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;
+			lastAuthSeqID = authSeqID;
 			continue;
 		}
-		
-		auto c = atomMapByAsymSeqAndAtom[make_tuple(asymID, lastSeqID, "C")];
-		auto n = atomMapByAsymSeqAndAtom[make_tuple(asymID, seqID, "N")];
 
-		if (not (c.empty() or n.empty()))
+		auto kc = make_tuple(asymID, lastSeqID, "C", lastAuthSeqID);
+		auto kn = make_tuple(asymID, seqID, "N", authSeqID);
+
+		if (atomMapByAsymSeqAndAtom.count(kc) and atomMapByAsymSeqAndAtom.count(kn))
+		{
+			auto c = atomMapByAsymSeqAndAtom.at(kc);
+			auto n = atomMapByAsymSeqAndAtom.at(kn);
+
 			bindAtoms(c, n);
-		
+		}
+		// if (not(c.empty() or n.empty()))
+
 		lastSeqID = seqID;
+		lastAuthSeqID = authSeqID;
 	}
 
-	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");
+		std::string authSeqId1, authSeqId2;
+
+		cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2, authSeqId1, authSeqId2) =
+			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_auth_seq_id", "ptnr2_auth_seq_id");
+
+		auto ka = make_tuple(asym1, seqId1, atomId1, authSeqId1);
+		auto kb = make_tuple(asym2, seqId2, atomId2, authSeqId2);
+
+		if (atomMapByAsymSeqAndAtom.count(ka) and atomMapByAsymSeqAndAtom.count(kb))
+		{
+			auto a = atomMapByAsymSeqAndAtom.at(ka);
+			auto b = atomMapByAsymSeqAndAtom.at(kb);
 
-		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()))
 			linkAtoms(a, b);
+		}
+
+		// std::string a = atomMapByAsymSeqAndAtom.at(make_tuple(asym1, seqId1, atomId1, authSeqId1));
+		// std::string b = atomMapByAsymSeqAndAtom.at(make_tuple(asym2, seqId2, atomId2, authSeqId2));
+		// 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 +341,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 +363,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 +381,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 +389,13 @@ BondMap::BondMap(const Structure& p)
 		}
 
 		// loop over pdbx_branch_scheme
-		for (auto r: db["pdbx_branch_scheme"].find(cif::Key("mon_id") == c))
+		for (const auto &[asym_id, pdb_seq_num] : db["pdbx_branch_scheme"].find<std::string,std::string>(cif::Key("mon_id") == c, "asym_id", "pdb_seq_num"))
 		{
-			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; });
-			
+				[&](const Atom &a)
+				{ return a.labelAsymID() == asym_id and a.authSeqID() == pdb_seq_num; });
+
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -531,31 +404,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 +439,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 +453,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

src/Cif2PDB.cpp

@@ -146,7 +146,7 @@ std::string cif2pdbSymmetry(std::string s)
 	return s;
 }
 
-std::string cif2pdbAtomName(std::string name, std::string resName, Datablock& db)
+std::string cif2pdbAtomName(std::string name, std::string resName, const Datablock& db)
 {
 	if (name.length() < 4)
 	{
@@ -166,7 +166,7 @@ std::string cif2pdbAtomName(std::string name, std::string resName, Datablock& db
 
 enum SoftwareType { eRefinement, eDataScaling, eDataExtraction, eDataReduction, ePhasing };
 
-std::string cifSoftware(Datablock& db, SoftwareType sw)
+std::string cifSoftware(const Datablock& db, SoftwareType sw)
 {
 	std::string result = "NULL";
 	
@@ -210,7 +210,7 @@ std::string cifSoftware(Datablock& db, SoftwareType sw)
 }
 
 // Map asym ID's back to PDB Chain ID's
-std::vector<std::string> MapAsymIDs2ChainIDs(const std::vector<std::string>& asymIDs, Datablock& db)
+std::vector<std::string> MapAsymIDs2ChainIDs(const std::vector<std::string>& asymIDs, const Datablock& db)
 {
 	std::set<std::string> result;
 	
@@ -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;
 		
@@ -275,7 +275,7 @@ size_t WriteOneContinuedLine(std::ostream& pdbFile, std::string header, int cLen
 	return WriteContinuedLine(pdbFile, header, count, cLen, line, lStart);
 }
 
-size_t WriteCitation(std::ostream& pdbFile, Datablock& db, Row r, int reference)
+size_t WriteCitation(std::ostream& pdbFile, const Datablock& db, Row r, int reference)
 {
 	size_t result = 0;
 	
@@ -361,7 +361,7 @@ size_t WriteCitation(std::ostream& pdbFile, Datablock& db, Row r, int reference)
 	return result;
 }
 
-void WriteHeaderLines(std::ostream& pdbFile, Datablock& db)
+void WriteHeaderLines(std::ostream& pdbFile, const Datablock& db)
 {
 	//    0         1         2         3         4         5         6         7         8
 	//    HEADER    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxDDDDDDDDD   IIII
@@ -592,7 +592,7 @@ void WriteHeaderLines(std::ostream& pdbFile, Datablock& db)
 		WriteOneContinuedLine(pdbFile, "AUTHOR  ", 2, ba::join(authors, ","));
 }
 
-void WriteTitle(std::ostream& pdbFile, Datablock& db)
+void WriteTitle(std::ostream& pdbFile, const Datablock& db)
 {
 	WriteHeaderLines(pdbFile, db);
 	
@@ -644,7 +644,7 @@ void WriteTitle(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark1(std::ostream& pdbFile, Datablock& db)
+void WriteRemark1(std::ostream& pdbFile, const Datablock& db)
 {
 	int reference = 0;
 	
@@ -662,7 +662,7 @@ void WriteRemark1(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark2(std::ostream& pdbFile, Datablock& db)
+void WriteRemark2(std::ostream& pdbFile, const Datablock& db)
 {
 	auto& refine = db["refine"];
 	if (refine.empty())
@@ -698,7 +698,7 @@ class FBase
   protected:
 	FBase(Row r, const char* f)
 		: mRow(r), mField(f) {}
-	FBase(Category& cat, cif::Condition&& cond, const char* f)
+	FBase(const Category& cat, cif::Condition&& cond, const char* f)
 		: mField(f)
 	{
 		auto r = cat.find(std::move(cond));
@@ -714,7 +714,7 @@ class Fi : public FBase
 {
   public:
 	Fi(Row r, const char* f) : FBase(r, f) {}
-	Fi(Category& cat, cif::Condition&& cond, const char* f) : FBase(cat, std::move(cond), f) {}
+	Fi(const Category& cat, cif::Condition&& cond, const char* f) : FBase(cat, std::move(cond), f) {}
 	
 	virtual void out(std::ostream& os)
 	{
@@ -734,7 +734,7 @@ class Ff : public FBase
 {
   public:
 	Ff(Row r, const char* f) : FBase(r, f) {}
-	Ff(Category& cat, cif::Condition&& cond, const char* f) : FBase(cat, std::move(cond), f) {}
+	Ff(const Category& cat, cif::Condition&& cond, const char* f) : FBase(cat, std::move(cond), f) {}
 	
 	virtual void out(std::ostream& os)
 	{
@@ -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;
             }
         }
@@ -764,7 +765,7 @@ class Fs : public FBase
 {
   public:
 	Fs(Row r, const char* f, int remarkNr = 3) : FBase(r, f), mNr(remarkNr) {}
-	Fs(Category& cat, cif::Condition&& cond, const char* f, int remarkNr = 3) : FBase(cat, std::move(cond), f), mNr(remarkNr) {}
+	Fs(const Category& cat, cif::Condition&& cond, const char* f, int remarkNr = 3) : FBase(cat, std::move(cond), f), mNr(remarkNr) {}
 	
 	virtual void out(std::ostream& os)
 	{
@@ -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,13 +825,13 @@ 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;
 }
 
 // --------------------------------------------------------------------
 
-void WriteRemark3BusterTNT(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3BusterTNT(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto ls_shell = db["refine_ls_shell"].front();
@@ -1008,7 +1009,7 @@ void WriteRemark3BusterTNT(std::ostream& pdbFile, Datablock& db)
 
 // --------------------------------------------------------------------
 
-void WriteRemark3CNS(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3CNS(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto ls_shell = db["refine_ls_shell"].front();
@@ -1148,7 +1149,7 @@ void WriteRemark3CNS(std::ostream& pdbFile, Datablock& db)
 
 // --------------------------------------------------------------------
 
-void WriteRemark3Refmac(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3Refmac(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto ls_shell = db["refine_ls_shell"].front();
@@ -1475,7 +1476,7 @@ void WriteRemark3Refmac(std::ostream& pdbFile, Datablock& db)
 			<< RM3("") << std::endl;
 }
 
-void WriteRemark3Shelxl(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3Shelxl(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 //	auto ls_shell = db["refine_ls_shell"].front();
@@ -1554,7 +1555,7 @@ void WriteRemark3Shelxl(std::ostream& pdbFile, Datablock& db)
 			<< RM3("") << std::endl;
 }
 
-void WriteRemark3Phenix(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3Phenix(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 //	auto ls_shell = db["refine_ls_shell"].front();
@@ -1781,7 +1782,7 @@ void WriteRemark3Phenix(std::ostream& pdbFile, Datablock& db)
 	pdbFile << RM3("") << std::endl;
 }
 
-void WriteRemark3XPlor(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3XPlor(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto ls_shell = db["refine_ls_shell"].front();
@@ -1896,7 +1897,7 @@ void WriteRemark3XPlor(std::ostream& pdbFile, Datablock& db)
 			<< RM3("") << std::endl;
 }
 
-void WriteRemark3NuclSQ(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3NuclSQ(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto pdbx_refine = db["pdbx_refine"].front();
@@ -2000,7 +2001,7 @@ void WriteRemark3NuclSQ(std::ostream& pdbFile, Datablock& db)
 			<< RM3("") << std::endl;
 }
 
-void WriteRemark3ProlSQ(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3ProlSQ(std::ostream& pdbFile, const Datablock& db)
 {
 	auto refine = db["refine"].front();
 	auto pdbx_refine = db["pdbx_refine"].front();
@@ -2119,7 +2120,7 @@ void WriteRemark3ProlSQ(std::ostream& pdbFile, Datablock& db)
 			<< RM3("") << std::endl;
 }
 
-void WriteRemark3(std::ostream& pdbFile, Datablock& db)
+void WriteRemark3(std::ostream& pdbFile, const Datablock& db)
 {
 	std::string program, authors;
 
@@ -2199,7 +2200,7 @@ void WriteRemark3(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark200(std::ostream& pdbFile, Datablock& db)
+void WriteRemark200(std::ostream& pdbFile, const Datablock& db)
 {
 	typedef RM<200> RM;
 	
@@ -2329,11 +2330,12 @@ 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;
 	}
 }
 
-void WriteRemark280(std::ostream& pdbFile, Datablock& db)
+void WriteRemark280(std::ostream& pdbFile, const Datablock& db)
 {
 	typedef RM<280> RM;
 	
@@ -2390,11 +2392,12 @@ 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;
 	}
 }
 
-void WriteRemark350(std::ostream& pdbFile, Datablock& db)
+void WriteRemark350(std::ostream& pdbFile, const Datablock& db)
 {
 	auto& c1 = db["pdbx_struct_assembly"];
 	if (c1.empty())
@@ -2514,7 +2517,7 @@ void WriteRemark350(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark400(std::ostream& pdbFile, Datablock& db)
+void WriteRemark400(std::ostream& pdbFile, const Datablock& db)
 {
 	for (auto& r: db["pdbx_entry_details"])
 	{
@@ -2524,7 +2527,7 @@ void WriteRemark400(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark450(std::ostream& pdbFile, Datablock& db)
+void WriteRemark450(std::ostream& pdbFile, const Datablock& db)
 {
 	for (auto& r: db["pdbx_entry_details"])
 	{
@@ -2535,7 +2538,7 @@ void WriteRemark450(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark465(std::ostream& pdbFile, Datablock& db)
+void WriteRemark465(std::ostream& pdbFile, const Datablock& db)
 {
 	bool first = true;
 	typedef RM<465> RM;
@@ -2584,7 +2587,7 @@ void WriteRemark465(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark470(std::ostream& pdbFile, Datablock& db)
+void WriteRemark470(std::ostream& pdbFile, const Datablock& db)
 {
 	typedef RM<470> RM;
 	boost::format fmt("REMARK 470 %3.3s %3.3s %1.1s%4.4d%1.1s  ");
@@ -2643,12 +2646,12 @@ void WriteRemark470(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark610(std::ostream& pdbFile, Datablock& db)
+void WriteRemark610(std::ostream& pdbFile, const Datablock& db)
 {
 	// #warning("unimplemented!");
 }
 
-void WriteRemark800(std::ostream& pdbFile, Datablock& db)
+void WriteRemark800(std::ostream& pdbFile, const Datablock& db)
 {
 	int nr = 0;
 	for (auto r: db["struct_site"])
@@ -2673,7 +2676,7 @@ void WriteRemark800(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemark999(std::ostream& pdbFile, Datablock& db)
+void WriteRemark999(std::ostream& pdbFile, const Datablock& db)
 {
 	for (auto& r: db["pdbx_entry_details"])
 	{
@@ -2684,7 +2687,7 @@ void WriteRemark999(std::ostream& pdbFile, Datablock& db)
 	}
 }
 
-void WriteRemarks(std::ostream& pdbFile, Datablock& db)
+void WriteRemarks(std::ostream& pdbFile, const Datablock& db)
 {
 	WriteRemark1(pdbFile, db);
 	WriteRemark2(pdbFile, db);
@@ -2706,7 +2709,7 @@ void WriteRemarks(std::ostream& pdbFile, Datablock& db)
 	WriteRemark999(pdbFile, db);
 }
 
-int WritePrimaryStructure(std::ostream& pdbFile, Datablock& db)
+int WritePrimaryStructure(std::ostream& pdbFile, const Datablock& db)
 {
 	int numSeq = 0;
 	
@@ -2860,7 +2863,7 @@ int WritePrimaryStructure(std::ostream& pdbFile, Datablock& db)
 	return numSeq;
 }
 
-int WriteHeterogen(std::ostream& pdbFile, Datablock& db)
+int WriteHeterogen(std::ostream& pdbFile, const Datablock& db)
 {
 	int numHet = 0;
 	
@@ -3129,7 +3132,7 @@ int WriteHeterogen(std::ostream& pdbFile, Datablock& db)
 	return numHet;
 }
 
-std::tuple<int,int> WriteSecondaryStructure(std::ostream& pdbFile, Datablock& db)
+std::tuple<int,int> WriteSecondaryStructure(std::ostream& pdbFile, const Datablock& db)
 {
 	int numHelix = 0, numSheet = 0;
 	
@@ -3287,7 +3290,7 @@ std::tuple<int,int> WriteSecondaryStructure(std::ostream& pdbFile, Datablock& db
 	return std::make_tuple(numHelix, numSheet);
 }
 
-void WriteConnectivity(std::ostream& pdbFile, cif::Datablock& db)
+void WriteConnectivity(std::ostream& pdbFile, const cif::Datablock& db)
 {
 	// SSBOND
 	// have to filter out alts
@@ -3403,7 +3406,7 @@ void WriteConnectivity(std::ostream& pdbFile, cif::Datablock& db)
 	}
 }
 
-int WriteMiscellaneousFeatures(std::ostream& pdbFile, Datablock& db)
+int WriteMiscellaneousFeatures(std::ostream& pdbFile, const Datablock& db)
 {
 	int numSite = 0;
 	
@@ -3456,7 +3459,7 @@ int WriteMiscellaneousFeatures(std::ostream& pdbFile, Datablock& db)
 	return numSite;
 }
 
-void WriteCrystallographic(std::ostream& pdbFile, Datablock& db)
+void WriteCrystallographic(std::ostream& pdbFile, const Datablock& db)
 {
 	auto r = db["symmetry"][cif::Key("entry_id") == db.getName()];
 	std::string symmetry = r["space_group_name_H-M"].as<std::string>();
@@ -3476,7 +3479,7 @@ void WriteCrystallographic(std::ostream& pdbFile, Datablock& db)
 		% r["Z_PDB"].as<double>()) << std::endl;
 }
 
-int WriteCoordinateTransformation(std::ostream& pdbFile, Datablock& db)
+int WriteCoordinateTransformation(std::ostream& pdbFile, const Datablock& db)
 {
 	int result = 0;
 	
@@ -3517,7 +3520,7 @@ int WriteCoordinateTransformation(std::ostream& pdbFile, Datablock& db)
 	return result;
 }
 
-std::tuple<int,int> WriteCoordinatesForModel(std::ostream& pdbFile, Datablock& db,
+std::tuple<int,int> WriteCoordinatesForModel(std::ostream& pdbFile, const Datablock& db,
 	const std::map<std::string,std::tuple<std::string,int,std::string>>& last_resseq_for_chain_map,
 	std::set<std::string>& TERminatedChains, int model_nr)
 {
@@ -3655,7 +3658,7 @@ std::tuple<int,int> WriteCoordinatesForModel(std::ostream& pdbFile, Datablock& d
 	return std::make_tuple(numCoord, numTer);
 }
 
-std::tuple<int,int> WriteCoordinate(std::ostream& pdbFile, Datablock& db)
+std::tuple<int,int> WriteCoordinate(std::ostream& pdbFile, const Datablock& db)
 {
 	// residues known from seqres
 //	map<tuple<std::string,int,std::string>,std::string> res2chain_map;
@@ -3715,7 +3718,7 @@ std::tuple<int,int> WriteCoordinate(std::ostream& pdbFile, Datablock& db)
 	return result;
 }
 
-void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile)
+void WritePDBFile(std::ostream& pdbFile, const Datablock &db)
 {
 	io::filtering_ostream out;
 	out.push(FillOutLineFilter());
@@ -3723,8 +3726,6 @@ void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile)
 
 	auto filter = out.component<FillOutLineFilter>(0);
 	assert(filter);
-
-	auto& db = cifFile.firstDatablock();
 	
 	int numRemark = 0, numHet = 0, numHelix = 0, numSheet = 0, numTurn = 0, numSite = 0, numXform = 0, numCoord = 0, numTer = 0, numConect = 0, numSeq = 0;
 	
@@ -3750,7 +3751,7 @@ void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile)
 		<< "END" << std::endl;
 }
 
-void WritePDBHeaderLines(std::ostream& os, cif::File& cifFile)
+void WritePDBHeaderLines(std::ostream& os, const Datablock &db)
 {
 	// io::filtering_ostream out;
 	// out.push(FillOutLineFilter());
@@ -3759,8 +3760,6 @@ void WritePDBHeaderLines(std::ostream& os, cif::File& cifFile)
 	// auto filter = out.component<FillOutLineFilter>(0);
 	// assert(filter);
 
-	auto& db = cifFile.firstDatablock();
-
 	WriteHeaderLines(os, db);
 }
 
@@ -3776,10 +3775,8 @@ std::string FixStringLength(const std::string& s, std::string::size_type l)
 	return result;
 }
 
-std::string GetPDBHEADERLine(cif::File& cifFile, std::string::size_type truncate_at)
+std::string GetPDBHEADERLine(const Datablock &db, std::string::size_type truncate_at)
 {
-	auto& db = cifFile.firstDatablock();
-
 	//    0         1         2         3         4         5         6         7         8
 	//    HEADER    xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxDDDDDDDDD   IIII
 	const char kHeader[] =
@@ -3822,10 +3819,8 @@ std::string GetPDBHEADERLine(cif::File& cifFile, std::string::size_type truncate
 	return FixStringLength((boost::format(kHeader) % keywords % date % db.getName()).str(), truncate_at);
 }
 
-std::string GetPDBCOMPNDLine(cif::File& cifFile, std::string::size_type truncate_at)
+std::string GetPDBCOMPNDLine(const Datablock &db, std::string::size_type truncate_at)
 {
-	auto& db = cifFile.firstDatablock();
-
 	// COMPND
 	using namespace std::placeholders;
 	
@@ -3883,10 +3878,8 @@ std::string GetPDBCOMPNDLine(cif::File& cifFile, std::string::size_type truncate
 	return FixStringLength("COMPND    " + ba::join(cmpnd, "; "), truncate_at);
 }
 
-std::string GetPDBSOURCELine(cif::File& cifFile, std::string::size_type truncate_at)
+std::string GetPDBSOURCELine(const Datablock &db, std::string::size_type truncate_at)
 {
-	auto& db = cifFile.firstDatablock();
-
 	// SOURCE
 	
 	int molID = 0;
@@ -3968,10 +3961,8 @@ std::string GetPDBSOURCELine(cif::File& cifFile, std::string::size_type truncate
 	return FixStringLength("SOURCE    " + ba::join(source, "; "), truncate_at);
 }
 
-std::string GetPDBAUTHORLine(cif::File& cifFile, std::string::size_type truncate_at)
+std::string GetPDBAUTHORLine(const Datablock &db, std::string::size_type truncate_at)
 {
-	auto& db = cifFile.firstDatablock();
-
 	// AUTHOR
 	std::vector<std::string> author;
 	for (auto r: db["audit_author"])

src/CifUtils.cpp

@@ -25,10 +25,7 @@
  */
 
 #include <atomic>
-#include <chrono>
 #include <cmath>
-#include <cstdio>
-#include <filesystem>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -37,7 +34,6 @@
 #include <regex>
 #include <sstream>
 #include <thread>
-#include <tuple>
 
 #if defined(_MSC_VER)
 #define TERM_WIDTH 80
@@ -50,6 +46,8 @@
 
 #include "cif++/CifUtils.hpp"
 
+#include "revision.hpp"
+
 namespace ba = boost::algorithm;
 namespace fs = std::filesystem;
 
@@ -64,39 +62,9 @@ extern int VERBOSE;
 
 std::string get_version_nr()
 {
-	const std::regex
-		rxVersionNr1(R"(build-(\d+)-g[0-9a-f]{7}(-dirty)?)"),
-		rxVersionNr2(R"(libcifpp-version: (\d+\.\d+\.\d+))");
-
-#include "revision.hpp"
-
-	struct membuf : public std::streambuf
-	{
-		membuf(char *data, size_t length) { this->setg(data, data, data + length); }
-	} buffer(const_cast<char *>(kRevision), sizeof(kRevision));
-
-	std::istream is(&buffer);
-
-	std::string line, result;
-
-	while (getline(is, line))
-	{
-		std::smatch m;
-
-		if (std::regex_match(line, m, rxVersionNr1))
-		{
-			result = m[1];
-			if (m[2].matched)
-				result += '*';
-			break;
-		}
-
-		// always the first, replace with more specific if followed by the other info
-		if (std::regex_match(line, m, rxVersionNr2))
-			result = m[1];
-	}
-
-	return result;
+	std::ostringstream s;
+	write_version_string(s, false);
+	return s.str();
 }
 
 // --------------------------------------------------------------------
@@ -123,11 +91,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 +109,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 +162,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");
@@ -674,7 +643,6 @@ void ProgressImpl::PrintProgress()
 
 	msg.append("| ");
 
-	//	const char		kSpinner[] = { '|', '/', '-', '\\' };
 	const char kSpinner[] = {' ', '.', 'o', 'O', '0', 'O', 'o', '.'};
 	const size_t kSpinnerCount = sizeof(kSpinner);
 
@@ -1236,7 +1204,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;
 }
@@ -1283,6 +1251,19 @@ std::unique_ptr<std::istream> loadResource(std::filesystem::path name)
 	}
 #endif
 
+#if defined(CCP4) and CCP4
+	if (not result and not fs::exists(p))
+	{
+		const char* CCP4_DIR = getenv("CCP4");
+		if (CCP4_DIR != nullptr and fs::exists(CCP4_DIR))
+		{
+			auto p2 = fs::path(DATA_DIR) / p;
+			if (fs::exists(p2))
+				swap(p, p2);
+		}
+	}
+#endif
+
 	if (not result and fs::exists(p))
 	{
 		std::unique_ptr<std::ifstream> file(new std::ifstream(p, std::ios::binary));

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

src/Compound.cpp

@@ -125,6 +125,11 @@ Compound::Compound(cif::Datablock &db)
 	cif::tie(mID, mName, mType, mFormula, mFormulaWeight, mFormalCharge) =
 		chemComp.front().get("id", "name", "type", "formula", "formula_weight", "pdbx_formal_charge");
 
+	// 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)
 	{
@@ -148,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)
@@ -187,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;
 		}
@@ -397,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";
@@ -408,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
@@ -514,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;
@@ -605,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";
@@ -614,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();
 		}
 	}
@@ -636,17 +642,18 @@ void CompoundFactory::init(bool useThreadLocalInstanceOnly)
 CompoundFactory::CompoundFactory()
 	: mImpl(nullptr)
 {
-	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)
-		std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;
-
 	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 > 0)
+		std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;
+
 }
 
 CompoundFactory::~CompoundFactory()
@@ -688,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;
 	}
 }
@@ -708,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;
 	}
 }

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"];

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)
@@ -58,13 +295,27 @@ Quaternion Normalize(Quaternion q)
 
 // --------------------------------------------------------------------
 
+Quaternion ConstructFromAngleAxis(float angle, Point axis)
+{
+	auto q = std::cos((angle * mmcif::kPI / 180) / 2);
+	auto s = std::sqrt(1 - q * q);
+
+	axis.normalize();
+
+	return Normalize(Quaternion{
+		static_cast<float>(q),
+		static_cast<float>(s * axis.mX),
+		static_cast<float>(s * axis.mY),
+		static_cast<float>(s * axis.mZ)});
+}
+
 std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q)
 {
 	if (q.R_component_1() > 1)
 		q = Normalize(q);
 
 	// angle:
-	double angle = 2 * acos(q.R_component_1());
+	double angle = 2 * std::acos(q.R_component_1());
 	angle = angle * 180 / kPI;
 
 	// axis:
@@ -102,14 +353,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 +426,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 +439,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 +476,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 +486,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;

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

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

src/revision.hpp.in

@@ -1,5 +0,0 @@
-const char kRevision[] = R"(
-lib@PROJECT_NAME@-version: @PROJECT_VERSION@
-@BUILD_VERSION_STRING@
-Date: @BUILD_DATE_TIME@
-)";

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;

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");
@@ -29,12 +32,12 @@ int main(int argc, char* argv[])
 		mmcif::File f(testdir / ".."/"examples"/"1cbs.cif.gz");
 		mmcif::Structure structure(f);
 
-		auto &res = structure.getResidue("B", "REA");
+		auto &res = structure.getResidue("B");
 		structure.changeResidue(res, "RXA", {});
 
 		structure.cleanupEmptyCategories();
 
-		f.file().save(std::cout);
+		f.save(std::cout);
 	}
 	catch (const std::exception& e)
 	{

test/structure-test.cpp

@@ -77,12 +77,9 @@ 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
+	file.loadDictionary("mmcif_pdbx_v50.dic");
+	file.emplace("TEST");	// create a datablock
 	
 	mmcif::Structure structure(file);
 
@@ -121,6 +118,17 @@ HETATM C  CHD . ? -4.342  36.262 -3.536  1.00 8.00  ?
 
 	auto expected = R"(
 data_TEST
+# 
+_pdbx_nonpoly_scheme.asym_id         A 
+_pdbx_nonpoly_scheme.ndb_seq_num     1 
+_pdbx_nonpoly_scheme.entity_id       1 
+_pdbx_nonpoly_scheme.mon_id          HEM 
+_pdbx_nonpoly_scheme.pdb_seq_num     1 
+_pdbx_nonpoly_scheme.auth_seq_num    1 
+_pdbx_nonpoly_scheme.pdb_mon_id      HEM 
+_pdbx_nonpoly_scheme.auth_mon_id     HEM 
+_pdbx_nonpoly_scheme.pdb_strand_id   A 
+_pdbx_nonpoly_scheme.pdb_ins_code    . 
 #
 loop_
 _atom_site.id
@@ -144,10 +152,10 @@ _atom_site.auth_seq_id
 _atom_site.auth_comp_id
 _atom_site.auth_atom_id
 _atom_site.pdbx_PDB_model_num
-1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? ? HEM CHA 1
-2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? ? HEM CHB 1
-3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? ? HEM CHC 1
-4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? ? HEM CHD 1
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? 1 HEM CHA 1
+2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? 1 HEM CHB 1
+3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? 1 HEM CHD 1
 #
 _chem_comp.id               HEM
 _chem_comp.type             NON-POLYMER
@@ -174,11 +182,178 @@ _struct_asym.details                       ?
 
 	expected.loadDictionary("mmcif_pdbx_v50.dic");
 
-	if (not (expected.firstDatablock() == structure.getFile().data()))
+	if (not (expected.firstDatablock() == structure.datablock()))
 	{
 		BOOST_TEST(false);
 		std::cout << expected.firstDatablock() << std::endl
 				<< std::endl
-				<< structure.getFile().data() << std::endl;
+				<< structure.datablock() << 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;
+// 		}
+// 	}
+// }
+
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(test_atom_id)
+{
+	auto data = R"(
+data_TEST
+# 
+_pdbx_nonpoly_scheme.asym_id         A 
+_pdbx_nonpoly_scheme.ndb_seq_num     1 
+_pdbx_nonpoly_scheme.entity_id       1 
+_pdbx_nonpoly_scheme.mon_id          HEM 
+_pdbx_nonpoly_scheme.pdb_seq_num     1 
+_pdbx_nonpoly_scheme.auth_seq_num    1 
+_pdbx_nonpoly_scheme.pdb_mon_id      HEM 
+_pdbx_nonpoly_scheme.auth_mon_id     HEM 
+_pdbx_nonpoly_scheme.pdb_strand_id   A 
+_pdbx_nonpoly_scheme.pdb_ins_code    . 
+#
+loop_
+_atom_site.id
+_atom_site.auth_asym_id
+_atom_site.label_alt_id
+_atom_site.label_asym_id
+_atom_site.label_atom_id
+_atom_site.label_comp_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.type_symbol
+_atom_site.group_PDB
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? 1 HEM CHA 1
+3 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? 1 HEM CHB 1
+2 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? 1 HEM CHD 1
+#
+_chem_comp.id               HEM
+_chem_comp.type             NON-POLYMER
+_chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
+_chem_comp.formula          'C34 H32 Fe N4 O4'
+_chem_comp.formula_weight   616.487000
+#
+_pdbx_entity_nonpoly.entity_id   1
+_pdbx_entity_nonpoly.name        'PROTOPORPHYRIN IX CONTAINING FE'
+_pdbx_entity_nonpoly.comp_id     HEM
+#
+_entity.id                 1
+_entity.type               non-polymer
+_entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
+_entity.formula_weight     616.487000
+#
+_struct_asym.id                            A
+_struct_asym.entity_id                     1
+_struct_asym.pdbx_blank_PDB_chainid_flag   N
+_struct_asym.pdbx_modified                 N
+_struct_asym.details                       ?
+#
+)"_cf;
+
+	data.loadDictionary("mmcif_pdbx_v50.dic");
+
+	mmcif::Structure s(data);
+
+	BOOST_CHECK_EQUAL(s.getAtomByID("1").authAtomID(), "CHA");
+	BOOST_CHECK_EQUAL(s.getAtomByID("2").authAtomID(), "CHC");
+	BOOST_CHECK_EQUAL(s.getAtomByID("3").authAtomID(), "CHB");
+	BOOST_CHECK_EQUAL(s.getAtomByID("4").authAtomID(), "CHD");
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(atom_numbers_1)
+{
+	const std::filesystem::path test1(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	mmcif::File file(test1.string());
+	mmcif::Structure structure(file);
+
+	auto &db = file.data();
+
+	auto &atoms = structure.atoms();
+	auto ai = atoms.begin();
+
+	for (const auto &[id, label_asym_id, label_seq_id, label_atom_id, auth_seq_id, label_comp_id] :
+		db["atom_site"].rows<std::string,std::string,int,std::string,std::string,std::string>("id", "label_asym_id", "label_seq_id", "label_atom_id", "auth_seq_id", "label_comp_id"))
+	{
+		auto atom = structure.getAtomByID(id);
+
+		BOOST_CHECK_EQUAL(atom.labelAsymID(), label_asym_id);
+		BOOST_CHECK_EQUAL(atom.labelSeqID(), label_seq_id);
+		BOOST_CHECK_EQUAL(atom.labelAtomID(), label_atom_id);
+		BOOST_CHECK_EQUAL(atom.authSeqID(), auth_seq_id);
+		BOOST_CHECK_EQUAL(atom.labelCompID(), label_comp_id);
+
+		BOOST_ASSERT(ai != atoms.end());
+
+		BOOST_CHECK_EQUAL(ai->id(), id);
+		++ai;
+	}
+
+	BOOST_ASSERT(ai == atoms.end());
+}
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(test_load_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	mmcif::File file(example.string());
+
+	auto &db = file.data();
+
+	mmcif::Structure s(file);
+
+	BOOST_CHECK(s.polymers().size() == 1);
+
+	auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
+
+	for (auto &poly : s.polymers())
+	{
+		BOOST_CHECK_EQUAL(poly.size(), pdbx_poly_seq_scheme.find("asym_id"_key == poly.asymID()).size());
+	}
+}
+
+BOOST_AUTO_TEST_CASE(remove_residue_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+	mmcif::File file(example.string());
+
+	mmcif::Structure s(file);
+	s.removeResidue(s.getResidue("B"));
+
+	BOOST_CHECK_NO_THROW(s.validateAtoms());
+}

test/sugar-test.cpp

@@ -0,0 +1,168 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ * 
+ * Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+#include <boost/test/included/unit_test.hpp>
+
+#include <stdexcept>
+
+#include "cif++/Cif++.hpp"
+#include "cif++/Structure.hpp"
+#include "cif++/CifValidator.hpp"
+
+// --------------------------------------------------------------------
+
+cif::File operator""_cf(const char* text, size_t length)
+{
+    struct membuf : public std::streambuf
+    {
+        membuf(char* text, size_t length)
+        {
+            this->setg(text, text, text + length);
+        }
+    } buffer(const_cast<char*>(text), length);
+
+    std::istream is(&buffer);
+    return cif::File(is);
+}
+
+// --------------------------------------------------------------------
+
+std::filesystem::path gTestDir = std::filesystem::current_path();
+
+bool init_unit_test()
+{
+    cif::VERBOSE = 1;
+
+	// not a test, just initialize test dir
+	if (boost::unit_test::framework::master_test_suite().argc == 2)
+		gTestDir = boost::unit_test::framework::master_test_suite().argv[1];
+
+	// do this now, avoids the need for installing
+	cif::addFileResource("mmcif_pdbx_v50.dic", gTestDir / ".." / "rsrc" / "mmcif_pdbx_v50.dic");
+
+	// initialize CCD location
+	cif::addFileResource("components.cif", gTestDir / ".." / "data" / "ccd-subset.cif");
+
+	mmcif::CompoundFactory::instance().pushDictionary(gTestDir / "HEM.cif");
+
+	return true;
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(sugar_name_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::File file(example.string());
+	mmcif::Structure s(file);
+
+	auto &db = s.datablock();
+	auto &entity = db["entity"];
+
+	auto &branches = s.branches();
+
+	BOOST_CHECK_EQUAL(branches.size(), 4);
+
+	for (auto &branch : branches)
+	{
+		auto entityID = branch.front().entityID();
+
+		auto name = entity.find1<std::string>("id"_key == entityID, "pdbx_description");
+		BOOST_CHECK_EQUAL(branch.name(), name);
+	}
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_sugar_1)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::File file(example.string());
+	mmcif::Structure s(file);
+
+	// collect atoms from asym L first
+	auto &NAG = s.getResidue("L");
+	auto nagAtoms = NAG.atoms();
+
+	std::vector<std::vector<cif::Item>> ai;
+
+	auto &db = s.datablock();
+	auto &as = db["atom_site"];
+
+	for (auto r : as.find("label_asym_id"_key == "L"))
+		ai.emplace_back(r.begin(), r.end());
+
+	s.removeResidue(NAG);
+
+	auto &branch = s.createBranch(ai);
+
+	BOOST_CHECK_EQUAL(branch.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose");
+	BOOST_CHECK_EQUAL(branch.size(), 1);
+
+	BOOST_CHECK_EQUAL(branch[0].atoms().size(), nagAtoms.size());
+}
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(create_sugar_2)
+{
+	using namespace cif::literals;
+
+	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
+	mmcif::File file(example.string());
+	mmcif::Structure s(file);
+
+	// Get branch for H
+	auto &bH = s.getBranchByAsymID("H");
+
+	BOOST_CHECK_EQUAL(bH.size(), 2);
+
+	std::vector<std::vector<cif::Item>> ai[2];
+
+	auto &db = s.datablock();
+	auto &as = db["atom_site"];
+
+	for (size_t i = 0; i < 2; ++i)
+	{
+		for (auto r : as.find("label_asym_id"_key == "H" and "auth_seq_id"_key == i + 1))
+			ai[i].emplace_back(r.begin(), r.end());
+	}
+
+	s.removeBranch(bH);
+
+	auto &bN = s.createBranch(ai[0]);
+	s.extendBranch(bN.asymID(), ai[1], 1, "O4");
+
+	BOOST_CHECK_EQUAL(bN.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose-(1-4)-2-acetamido-2-deoxy-beta-D-glucopyranose");
+	BOOST_CHECK_EQUAL(bN.size(), 2);
+
+	file.save(gTestDir / "test-create_sugar_2.cif");
+}