Fixed all tests

- Optimize isUnquotedString

.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 3.0.1 LANGUAGES CXX)
+project(cifpp VERSION 4.0.0 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")
@@ -177,25 +159,18 @@ find_package(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_opt
 
 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
-
-include(GetGitRevisionDescription)
-option(GENERATE_CUSTOM_VERSION "Generate a custom version string" OFF)
-if(GIT-NOTFOUND OR HEAD-HASH-NOTFOUND OR NOT GENERATE_CUSTOM_VERSION)
-    get_git_head_revision(REFSPEC COMMITHASH)
-
-	# Generate our own version string
-	git_describe_working_tree(BUILD_VERSION_STRING --match=build --dirty)
-else()
-    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)
@@ -203,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()
@@ -274,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)
@@ -334,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::
@@ -407,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)

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,9 +1,37 @@
+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.
 

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

@@ -172,19 +172,26 @@ enum AtomType : uint8_t
 // --------------------------------------------------------------------
 // AtomTypeInfo
 
-enum RadiusType
+enum class RadiusType
 {
-	eRadiusCalculated,
-	eRadiusEmpirical,
-	eRadiusCovalentEmpirical,
+	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
@@ -194,7 +201,7 @@ struct AtomTypeInfo
 	std::string symbol;
 	float weight;
 	bool metal;
-	float radii[eRadiusTypeCount];
+	float radii[RadiusTypeCount];
 };
 
 extern const AtomTypeInfo kKnownAtoms[];
@@ -218,11 +225,32 @@ class AtomTypeTraits
 	static bool isElement(const std::string &symbol);
 	static bool isMetal(const std::string &symbol);
 
-	float radius(RadiusType type = eRadiusSingleBond) const
+	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

include/cif++/Cif++.hpp

@@ -176,6 +176,7 @@ class Item
 		, mValue(rhs.mValue)
 	{
 	}
+
 	Item(Item &&rhs) noexcept
 		: mName(std::move(rhs.mName))
 		, mValue(std::move(rhs.mValue))
@@ -792,7 +793,7 @@ class Row
 	}
 
 	void assign(const std::vector<Item> &values);
-	void assign(std::string_view 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
 	{
@@ -814,7 +815,7 @@ 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);
@@ -1183,6 +1184,9 @@ struct Empty
 {
 };
 
+/// \brief A helper to make it possible to have conditions like ("id"_key == cif::null)
+inline constexpr Empty null = Empty();
+
 struct Key
 {
 	explicit Key(const std::string &itemTag)
@@ -2152,6 +2156,7 @@ class Category
 	std::vector<ItemColumn> mColumns;
 	ItemRow *mHead;
 	ItemRow *mTail;
+	size_t mLastUniqueNr = 0;
 	class CatIndex *mIndex;
 
 	std::vector<Linked> mParentLinks, mChildLinks;
@@ -2168,17 +2173,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);
@@ -2210,6 +2219,12 @@ class File
 
 	void append(Datablock *e);
 
+	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);
 
@@ -2247,6 +2262,9 @@ class File
 	iterator begin() const;
 	iterator end() const;
 
+	Datablock &front();
+	Datablock &back();
+
 	bool empty() const { return mHead == nullptr; }
 
 	const Validator &getValidator() const;

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
@@ -82,16 +82,7 @@ inline bool isAnyPrint(int ch)
 	       (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);
 
 // --------------------------------------------------------------------
 
@@ -139,7 +130,7 @@ class SacParser
 	int getNextChar();
 
 	void retract();
-	void restart();
+	int restart(int start);
 
 	CIFToken getNextToken();
 	void match(CIFToken token);
@@ -181,8 +172,9 @@ class SacParser
 		eStateTextField,
 		eStateFloat = 100,
 		eStateInt = 110,
-		//		eStateNumericSuffix = 200,
-		eStateValue = 300
+		eStateValue = 300,
+		eStateDATA,
+		eStateSAVE
 	};
 
 	std::istream &mData;
@@ -191,7 +183,6 @@ class SacParser
 	bool mValidate;
 	uint32_t mLineNr;
 	bool mBol;
-	int mState, mStart;
 	CIFToken mLookahead;
 	std::string mTokenValue;
 	CIFValueType mTokenType;

include/cif++/CifUtils.hpp

@@ -26,11 +26,7 @@
 
 #pragma once
 
-#include <cassert>
 #include <filesystem>
-#include <iostream>
-#include <list>
-#include <memory>
 #include <set>
 #include <vector>
 

include/cif++/Point.hpp

@@ -221,7 +221,7 @@ struct PointF
 
 	FType length() const
 	{
-		return sqrt(mX * mX + mY * mY + mZ * mZ);
+		return std::sqrt(mX * mX + mY * mY + mZ * mZ);
 	}
 };
 
@@ -285,7 +285,7 @@ inline double DistanceSquared(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));
@@ -332,10 +332,10 @@ double DihedralAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &
 	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;
@@ -351,7 +351,7 @@ double CosinusAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p
 
 	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;
 }
@@ -371,15 +371,16 @@ 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);
 
+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);
 
@@ -437,9 +438,9 @@ class SphericalDots
 			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;
 		}

include/cif++/Secondary.hpp

@@ -184,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) {}
 
@@ -193,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*;
@@ -214,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_view name) const;
-
-	void property(const std::string_view name, const std::string &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)
-	{
-		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;
 
@@ -242,6 +351,8 @@ class Residue
 	int seqID() const { return mSeqID; }
 	std::string entityID() const;
 
+	EntityType entityType() const;
+
 	std::string authAsymID() const;
 	std::string authSeqID() const;
 	std::string authInsCode() const;
@@ -277,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() {}
 
@@ -285,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;
 };
@@ -354,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;
@@ -391,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_view 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(); }
 };
 
 // --------------------------------------------------------------------
@@ -439,25 +616,46 @@ 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,
@@ -469,45 +667,64 @@ class Structure
 	/// \brief Return the atom closest to point \a p with atom type \a type in a residue of type \a res_type
 	Atom getAtomByPositionAndType(Point p, std::string_view type, std::string_view res_type) const;
 
-	/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
-	const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0) const;
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	Residue &getResidue(const std::string &asymID)
+	{
+		return getResidue(asymID, 0, "");
+	}
 
-	/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
-	Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0);
+	/// \brief Get a non-poly residue for an asym with id \a asymID
+	const Residue &getResidue(const std::string &asymID) const
+	{
+		return getResidue(asymID, 0, "");
+	}
 
-	/// \brief Get a the single residue for an asym with id \a asymID
-	const Residue &getResidue(const std::string &asymID) const;
+	/// \brief Get a 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);
 
-	/// \brief Get a the single residue for an asym with id \a asymID
-	Residue &getResidue(const std::string &asymID);
+	/// \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);
+	}
 
-	// map between auth and label locations
+	/// \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);
 
-	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 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);
+	}
 
-	std::tuple<std::string, std::string, std::string, std::string> MapLabelToAuth(
-		const std::string &asymID, int seqID, const std::string &compID);
+	/// \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());
+	}
 
-	// returns chain, seqnr, icode
-	std::tuple<char, int, char> MapLabelToAuth(
-		const std::string &asymID, int seqID) const;
-
-	// 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;
-
-	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 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 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
@@ -521,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);
@@ -531,34 +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;
-	cif::Datablock &datablock() const;
+	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;
 
 	std::string insertCompound(const std::string &compoundID, bool isEntity);
 
+	std::string createEntityForBranch(Branch &branch);
+
 	void loadData();
-	void updateAtomIndex();
 
 	void loadAtomsForModel(StructureOpenOptions options);
 
-	File &mFile;
+	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

src/AtomType.cpp

@@ -159,11 +159,221 @@ const AtomTypeInfo kKnownAtoms[] =
 	{ Ts,	"Tennessine",		"Ts",	294,		true,  {	kNA,	kNA,	kNA,	165,	kNA,	kNA,	kNA } },  //	117	Ts	 Tenness­ine       
 	{ Og,	"Oganesson",		"Og",	294,		true,  {	kNA,	kNA,	kNA,	157,	kNA,	kNA,	kNA } },  //	118	Og	 Oga­nesson        
 
-	{ D,	"Deuterium",		"D",	2.014f,	false, {	53,		25,		37,		32,		kNA,	kNA,	120 } },  //	1	D	 Deuterium         
+	{ D,	"Deuterium",		"D",	2.014f,		false, {	53,		25,		37,		32,		kNA,	kNA,	120 } },  //	1	D	 Deuterium         
 };                                                                                                                                                    
 
 uint32_t kKnownAtomsCount = sizeof(kKnownAtoms) / sizeof(AtomTypeInfo);
 
+// --------------------------------------------------------------------
+// Crystal ionic radii, as taken from Wikipedia (https://en.m.wikipedia.org/wiki/Ionic_radius)
+
+const struct IonicRadii
+{
+	AtomType type;
+	float radii[11];
+} kCrystalIonicRadii[] = {
+	{ H,	{ kNA,	kNA,	208,	-4,		kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Hydrogen 
+	{ Li,	{ kNA,	kNA,	kNA,	90,		kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Lithium 
+	{ Be,	{ kNA,	kNA,	kNA,	kNA,	59,		kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Beryllium 
+	{ B,	{ kNA,	kNA,	kNA,	kNA,	kNA,	41,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Boron 
+	{ C,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	30,		kNA,	kNA,	kNA,	kNA } },	// Carbon 
+	{ N,	{ 132,	kNA,	kNA,	kNA,	kNA,	30,		kNA,	27,		kNA,	kNA,	kNA } },	// Nitrogen 
+	{ O,	{ kNA,	126,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Oxygen 
+	{ F,	{ kNA,	kNA,	119,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	22,		kNA } },	// Fluorine 
+	{ Na,	{ kNA,	kNA,	kNA,	116,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Sodium 
+	{ Mg,	{ kNA,	kNA,	kNA,	kNA,	86,		kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Magnesium 
+	{ Al,	{ kNA,	kNA,	kNA,	kNA,	kNA,	67.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Aluminium 
+	{ Si,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	54,		kNA,	kNA,	kNA,	kNA } },	// Silicon 
+	{ P,	{ kNA,	kNA,	kNA,	kNA,	kNA,	58,		kNA,	52,		kNA,	kNA,	kNA } },	// Phosphorus 
+	{ S,	{ kNA,	170,	kNA,	kNA,	kNA,	kNA,	51,		kNA,	43,		kNA,	kNA } },	// Sulfur 
+	{ Cl,	{ kNA,	kNA,	181,	kNA,	kNA,	kNA,	kNA,	26,		kNA,	41,		kNA } },	// Chlorine 
+	{ K,	{ kNA,	kNA,	kNA,	152,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Potassium 
+	{ Ca,	{ kNA,	kNA,	kNA,	kNA,	114,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Calcium 
+	{ Sc,	{ kNA,	kNA,	kNA,	kNA,	kNA,	88.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Scandium 
+	{ Ti,	{ kNA,	kNA,	kNA,	kNA,	100,	81,		74.5,	kNA,	kNA,	kNA,	kNA } },	// Titanium 
+	{ V,	{ kNA,	kNA,	kNA,	kNA,	93,		78,		72,		68,		kNA,	kNA,	kNA } },	// Vanadium 
+	{ Cr,	{ kNA,	kNA,	kNA,	kNA,	87,		75.5,	69,		63,		58,		kNA,	kNA } },	// Chromium ls 
+	// { Cr,{ kNA,	kNA,	kNA,	kNA,	94,		kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Chromium hs 
+	{ Mn,	{ kNA,	kNA,	kNA,	kNA,	81,		72,		67,		47,		39.5,	60,		kNA } },	// Manganese ls 
+	// { Mn,{ kNA,	kNA,	kNA,	kNA,	97,		78.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Manganese hs 
+	{ Fe,	{ kNA,	kNA,	kNA,	kNA,	75,		69,		72.5,	kNA,	39,		kNA,	kNA } },	// Iron ls 
+	// { Fe,{ kNA,	kNA,	kNA,	kNA,	92,		78.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Iron hs 
+	{ Co,	{ kNA,	kNA,	kNA,	kNA,	79,		68.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Cobalt ls 
+	// { Co,{ kNA,	kNA,	kNA,	kNA,	88.5,	75,		67,		kNA,	kNA,	kNA,	kNA } },	// Cobalt hs 
+	{ Ni,	{ kNA,	kNA,	kNA,	kNA,	83,		70,		62,		kNA,	kNA,	kNA,	kNA } },	// Nickel ls 
+	// { Ni,{ kNA,	kNA,	kNA,	kNA,	kNA,	74,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Nickel hs 
+	{ Cu,	{ kNA,	kNA,	kNA,	91,		87,		68,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Copper 
+	{ Zn,	{ kNA,	kNA,	kNA,	kNA,	88	,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Zinc 
+	{ Ga,	{ kNA,	kNA,	kNA,	kNA,	kNA,	76,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Gallium 
+	{ Ge,	{ kNA,	kNA,	kNA,	kNA,	87,		kNA,	67,		kNA,	kNA,	kNA,	kNA } },	// Germanium 
+	{ As,	{ kNA,	kNA,	kNA,	kNA,	kNA,	72,		kNA,	60,		kNA,	kNA,	kNA } },	// Arsenic 
+	{ Se,	{ kNA,	184,	kNA,	kNA,	kNA,	kNA,	64,		kNA,	56,		kNA,	kNA } },	// Selenium 
+	{ Br,	{ kNA,	kNA,	182,	kNA,	kNA,	73,		kNA,	45,		kNA,	53,		kNA } },	// Bromine 
+	{ Rb,	{ kNA,	kNA,	kNA,	166,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Rubidium 
+	{ Sr,	{ kNA,	kNA,	kNA,	kNA,	132,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Strontium 
+	{ Y,	{ kNA,	kNA,	kNA,	kNA,	kNA,	104,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Yttrium 
+	{ Zr,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	86,		kNA,	kNA,	kNA,	kNA } },	// Zirconium 
+	{ Nb,	{ kNA,	kNA,	kNA,	kNA,	kNA,	86,		82,		78,		kNA,	kNA,	kNA } },	// Niobium 
+	{ Mo,	{ kNA,	kNA,	kNA,	kNA,	kNA,	83,		79,		75,		73,		kNA,	kNA } },	// Molybdenum 
+	{ Tc,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	78.5,	74,		kNA,	70,		kNA } },	// Technetium 
+	{ Ru,	{ kNA,	kNA,	kNA,	kNA,	kNA,	82,		76,		70.5,	kNA,	52,		150 } },	// Ruthenium 
+	{ Rh,	{ kNA,	kNA,	kNA,	kNA,	kNA,	80.5,	74,		69,		kNA,	kNA,	kNA } },	// Rhodium 
+	{ Pd,	{ kNA,	kNA,	kNA,	73,		100,	90,		75.5,	kNA,	kNA,	kNA,	kNA } },	// Palladium 
+	{ Ag,	{ kNA,	kNA,	kNA,	129,	108,	89,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Silver 
+	{ Cd,	{ kNA,	kNA,	kNA,	kNA,	109,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Cadmium 
+	{ In,	{ kNA,	kNA,	kNA,	kNA,	kNA,	94,		kNA,	kNA,	kNA,	kNA,	kNA } },	// Indium 
+	{ Sn,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	83,		kNA,	kNA,	kNA,	kNA } },	// Tin 
+	{ Sb,	{ kNA,	kNA,	kNA,	kNA,	kNA,	90,		kNA,	74,		kNA,	kNA,	kNA } },	// Antimony 
+	{ Te,	{ kNA,	207,	kNA,	kNA,	kNA,	kNA,	111,	kNA,	70,		kNA,	kNA } },	// Tellurium 
+	{ I,	{ kNA,	kNA,	206,	kNA,	kNA,	kNA,	kNA,	109,	kNA,	67,		kNA } },	// Iodine 
+	{ Xe,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	62 } },		// Xenon 
+	{ Cs,	{ kNA,	kNA,	kNA,	167,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Caesium 
+	{ Ba,	{ kNA,	kNA,	kNA,	kNA,	149,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Barium 
+	{ La,	{ kNA,	kNA,	kNA,	kNA,	kNA,	117.2,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Lanthanum 
+	{ Ce,	{ kNA,	kNA,	kNA,	kNA,	kNA,	115,	101,	kNA,	kNA,	kNA,	kNA } },	// Cerium 
+	{ Pr,	{ kNA,	kNA,	kNA,	kNA,	kNA,	113,	99,		kNA,	kNA,	kNA,	kNA } },	// Praseodymium 
+	{ Nd,	{ kNA,	kNA,	kNA,	kNA,	143,	112.3,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Neodymium 
+	{ Pm,	{ kNA,	kNA,	kNA,	kNA,	kNA,	111,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Promethium 
+	{ Sm,	{ kNA,	kNA,	kNA,	kNA,	136,	109.8,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Samarium 
+	{ Eu,	{ kNA,	kNA,	kNA,	kNA,	131,	108.7,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Europium 
+	{ Gd,	{ kNA,	kNA,	kNA,	kNA,	kNA,	107.8,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Gadolinium 
+	{ Tb,	{ kNA,	kNA,	kNA,	kNA,	kNA,	106.3,	90,		kNA,	kNA,	kNA,	kNA } },	// Terbium 
+	{ Dy,	{ kNA,	kNA,	kNA,	kNA,	121,	105.2,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Dysprosium 
+	{ Ho,	{ kNA,	kNA,	kNA,	kNA,	kNA,	104.1,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Holmium 
+	{ Er,	{ kNA,	kNA,	kNA,	kNA,	kNA,	103,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Erbium 
+	{ Tm,	{ kNA,	kNA,	kNA,	kNA,	117,	102,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Thulium 
+	{ Yb,	{ kNA,	kNA,	kNA,	kNA,	116,	100.8,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Ytterbium 
+	{ Lu,	{ kNA,	kNA,	kNA,	kNA,	kNA,	100.1,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Lutetium 
+	{ Hf,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	85,		kNA,	kNA,	kNA,	kNA } },	// Hafnium 
+	{ Ta,	{ kNA,	kNA,	kNA,	kNA,	kNA,	86,		82,		78,		kNA,	kNA,	kNA } },	// Tantalum 
+	{ W,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	80,		76,		74,		kNA,	kNA } },	// Tungsten 
+	{ Re,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	77,		72,		69,		67,		kNA } },	// Rhenium 
+	{ Os,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	77,		71.5,	68.5,	66.5,	53 } },		// Osmium 
+	{ Ir,	{ kNA,	kNA,	kNA,	kNA,	kNA,	82,		76.5,	71,		kNA,	kNA,	kNA } },	// Iridium 
+	{ Pt,	{ kNA,	kNA,	kNA,	kNA,	94,		kNA,	76.5,	71,		kNA,	kNA,	kNA } },	// Platinum 
+	{ Au,	{ kNA,	kNA,	kNA,	151,	kNA,	99,		kNA,	71,		kNA,	kNA,	kNA } },	// Gold 
+	{ Hg,	{ kNA,	kNA,	kNA,	133,	116,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Mercury 
+	{ Tl,	{ kNA,	kNA,	kNA,	164,	kNA,	102.5,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Thallium 
+	{ Pb,	{ kNA,	kNA,	kNA,	kNA,	133,	kNA,	91.5,	kNA,	kNA,	kNA,	kNA } },	// Lead 
+	{ Bi,	{ kNA,	kNA,	kNA,	kNA,	kNA,	117,	kNA,	90,		kNA,	kNA,	kNA } },	// Bismuth 
+	{ Po,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	108,	kNA,	81,		kNA,	kNA } },	// Polonium 
+	{ At,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	76,		kNA } },	// Astatine 
+	{ Fr,	{ kNA,	kNA,	kNA,	194,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Francium 
+	{ Ra,	{ kNA,	kNA,	kNA,	kNA,	162,	kNA,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Radium 
+	{ Ac,	{ kNA,	kNA,	kNA,	kNA,	kNA,	126,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Actinium 
+	{ Th,	{ kNA,	kNA,	kNA,	kNA,	kNA,	kNA,	108,	kNA,	kNA,	kNA,	kNA } },	// Thorium 
+	{ Pa,	{ kNA,	kNA,	kNA,	kNA,	kNA,	116,	104,	92,		kNA,	kNA,	kNA } },	// Protactinium 
+	{ U,	{ kNA,	kNA,	kNA,	kNA,	kNA,	116.5,	103,	90,		87,		kNA,	kNA } },	// Uranium 
+	{ Np,	{ kNA,	kNA,	kNA,	kNA,	124,	115,	101,	89,		86,		85,		kNA } },	// Neptunium 
+	{ Pu,	{ kNA,	kNA,	kNA,	kNA,	kNA,	114,	100,	88,		85,		kNA,	kNA } },	// Plutonium 
+	{ Am,	{ kNA,	kNA,	kNA,	kNA,	140,	111.5,	99,		kNA,	kNA,	kNA,	kNA } },	// Americium 
+	{ Cm,	{ kNA,	kNA,	kNA,	kNA,	kNA,	111,	99,		kNA,	kNA,	kNA,	kNA } },	// Curium 
+	{ Bk,	{ kNA,	kNA,	kNA,	kNA,	kNA,	110,	97,		kNA,	kNA,	kNA,	kNA } },	// Berkelium 
+	{ Cf,	{ kNA,	kNA,	kNA,	kNA,	kNA,	109,	96.1,	kNA,	kNA,	kNA,	kNA } },	// Californium 
+	{ Es,	{ kNA,	kNA,	kNA,	kNA,	kNA,	92.8,	kNA,	kNA,	kNA,	kNA,	kNA } },	// Einsteinium 	
+}, kEffectiveIonicRadii[] = {
+	{ H,	{ kNA, 	kNA, 	139.9,	-18,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Hydrogen 
+	{ Li,	{ kNA, 	kNA, 	kNA, 	76,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Lithium 
+	{ Be,	{ kNA, 	kNA, 	kNA, 	kNA, 	45,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Beryllium 
+	{ B,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	27,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Boron 
+	{ C,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	16,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Carbon 
+	{ N,	{ 146,	kNA, 	kNA, 	kNA, 	kNA, 	16,		kNA, 	13,		kNA, 	kNA, 	kNA } }, 	// Nitrogen 
+	{ O,	{ kNA, 	140,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Oxygen 
+	{ F,	{ kNA, 	kNA, 	133,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	8,		kNA } }, 	// Fluorine 
+	{ Na,	{ kNA, 	kNA, 	kNA, 	102,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Sodium 
+	{ Mg,	{ kNA, 	kNA, 	kNA, 	kNA, 	72,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Magnesium 
+	{ Al,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	53.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Aluminium 
+	{ Si,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	40,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Silicon 
+	{ P,	{ 212,	kNA, 	kNA, 	kNA, 	kNA, 	44,		kNA, 	38,		kNA, 	kNA, 	kNA } }, 	// Phosphorus 
+	{ S,	{ kNA, 	184,	kNA, 	kNA, 	kNA, 	kNA, 	37,		kNA, 	29,		kNA, 	kNA } }, 	// Sulfur 
+	{ Cl,	{ kNA, 	kNA, 	181,	kNA, 	kNA, 	kNA, 	kNA, 	12,		kNA, 	27,		kNA } }, 	// Chlorine 
+	{ K,	{ kNA, 	kNA, 	kNA, 	138,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Potassium 
+	{ Ca,	{ kNA, 	kNA, 	kNA, 	kNA, 	100,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Calcium 
+	{ Sc,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	74.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Scandium 
+	{ Ti,	{ kNA, 	kNA, 	kNA, 	kNA, 	86,		67,		60.5,	kNA, 	kNA, 	kNA, 	kNA } }, 	// Titanium 
+	{ V,	{ kNA, 	kNA, 	kNA, 	kNA, 	79,		64,		58,		54,		kNA, 	kNA, 	kNA } }, 	// Vanadium 
+	{ Cr,	{ kNA, 	kNA, 	kNA, 	kNA, 	73,		61.5,	55,		49,		44,		kNA, 	kNA } }, 	// Chromium ls 
+	{ Cr,	{ kNA, 	kNA, 	kNA, 	kNA, 	80,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Chromium hs 
+	{ Mn,	{ kNA, 	kNA, 	kNA, 	kNA, 	67,		58,		53,		33,		25.5,	46,		kNA } }, 	// Manganese ls 
+	{ Mn,	{ kNA, 	kNA, 	kNA, 	kNA, 	83,		64.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Manganese hs 
+	{ Fe,	{ kNA, 	kNA, 	kNA, 	kNA, 	61,		55,		58.5,	kNA, 	25,		kNA, 	kNA } }, 	// Iron ls 
+	{ Fe,	{ kNA, 	kNA, 	kNA, 	kNA, 	78,		64.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Iron hs 
+	{ Co,	{ kNA, 	kNA, 	kNA, 	kNA, 	65,		54.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Cobalt ls 
+	{ Co,	{ kNA, 	kNA, 	kNA, 	kNA, 	74.5,	61,		53,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Cobalt hs 
+	{ Ni,	{ kNA, 	kNA, 	kNA, 	kNA, 	69,		56,		48,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Nickel ls 
+	{ Ni,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	60,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Nickel hs 
+	{ Cu,	{ kNA, 	kNA, 	kNA, 	77,		73,		54, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Copper 
+	{ Zn,	{ kNA, 	kNA, 	kNA, 	kNA, 	74,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Zinc 
+	{ Ga,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	62,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Gallium 
+	{ Ge,	{ kNA, 	kNA, 	kNA, 	kNA, 	73,		kNA, 	53,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Germanium 
+	{ As,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	58,		kNA, 	46,		kNA, 	kNA, 	kNA } }, 	// Arsenic 
+	{ Se,	{ kNA, 	198,	kNA, 	kNA, 	kNA, 	kNA, 	50,		kNA, 	42,		kNA, 	kNA } }, 	// Selenium 
+	{ Br,	{ kNA, 	kNA, 	196,	kNA, 	kNA, 	59,		kNA, 	31,		kNA, 	39,		kNA } }, 	// Bromine 
+	{ Rb,	{ kNA, 	kNA, 	kNA, 	152,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Rubidium 
+	{ Sr,	{ kNA, 	kNA, 	kNA, 	kNA, 	118,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Strontium 
+	{ Y,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	90,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Yttrium 
+	{ Zr,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	72,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Zirconium 
+	{ Nb,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	72,		68,		64,		kNA, 	kNA, 	kNA } }, 	// Niobium 
+	{ Mo,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	69,		65,		61,		59,		kNA, 	kNA } }, 	// Molybdenum 
+	{ Tc,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	64.5,	60,		kNA, 	56,		kNA } }, 	// Technetium 
+	{ Ru,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	68,		62,		56.5,	kNA, 	38,		36 } },		// Ruthenium 
+	{ Rh,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	66.5,	60,		55,		kNA, 	kNA, 	kNA } }, 	// Rhodium 
+	{ Pd,	{ kNA, 	kNA, 	kNA, 	59,		86,		76,		61.5,	kNA, 	kNA, 	kNA, 	kNA } }, 	// Palladium 
+	{ Ag,	{ kNA, 	kNA, 	kNA, 	115,	94,		75,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Silver 
+	{ Cd,	{ kNA, 	kNA, 	kNA, 	kNA, 	95,		kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Cadmium 
+	{ In,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	80,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Indium 
+	{ Sn,	{ kNA, 	kNA, 	kNA, 	kNA, 	118,	kNA, 	69,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Tin 
+	{ Sb,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	76,		kNA, 	60,		kNA, 	kNA, 	kNA } }, 	// Antimony 
+	{ Te,	{ kNA, 	221,	kNA, 	kNA, 	kNA, 	kNA, 	97,		kNA, 	56,		kNA, 	kNA } }, 	// Tellurium 
+	{ I,	{ kNA, 	kNA, 	220,	kNA, 	kNA, 	kNA, 	kNA, 	95,		kNA, 	53,		kNA } }, 	// Iodine 
+	{ Xe,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	48 } },		// Xenon 
+	{ Cs,	{ kNA, 	kNA, 	kNA, 	167,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Caesium 
+	{ Ba,	{ kNA, 	kNA, 	kNA, 	kNA, 	135,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Barium 
+	{ La,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	103.2,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Lanthanum 
+	{ Ce,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	101,	87,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Cerium 
+	{ Pr,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	99,		85,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Praseodymium 
+	{ Nd,	{ kNA, 	kNA, 	kNA, 	kNA, 	129,	98.3,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Neodymium 
+	{ Pm,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	97,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Promethium 
+	{ Sm,	{ kNA, 	kNA, 	kNA, 	kNA, 	122,	95.8,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Samarium 
+	{ Eu,	{ kNA, 	kNA, 	kNA, 	kNA, 	117,	94.7,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Europium 
+	{ Gd,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	93.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Gadolinium 
+	{ Tb,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	92.3,	76,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Terbium 
+	{ Dy,	{ kNA, 	kNA, 	kNA, 	kNA, 	107,	91.2,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Dysprosium 
+	{ Ho,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	90.1,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Holmium 
+	{ Er,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	89,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Erbium 
+	{ Tm,	{ kNA, 	kNA, 	kNA, 	kNA, 	103,	88,		kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Thulium 
+	{ Yb,	{ kNA, 	kNA, 	kNA, 	kNA, 	102,	86.8,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Ytterbium 
+	{ Lu,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	86.1,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Lutetium 
+	{ Hf,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	71,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Hafnium 
+	{ Ta,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	72,		68,		64,		kNA, 	kNA, 	kNA } }, 	// Tantalum 
+	{ W,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	66,		62,		60,		kNA, 	kNA } }, 	// Tungsten 
+	{ Re,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	63,		58,		55,		53,		kNA } }, 	// Rhenium 
+	{ Os,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	63,		57.5,	54.5,	52.5,	39 } },		// Osmium 
+	{ Ir,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	68,		62.5,	57,		kNA, 	kNA, 	kNA } }, 	// Iridium 
+	{ Pt,	{ kNA, 	kNA, 	kNA, 	kNA, 	80,		kNA, 	62.5,	57,		kNA, 	kNA, 	kNA } }, 	// Platinum 
+	{ Au,	{ kNA, 	kNA, 	kNA, 	137,	kNA, 	85,		kNA, 	57,		kNA, 	kNA, 	kNA } }, 	// Gold 
+	{ Hg,	{ kNA, 	kNA, 	kNA, 	119,	102,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Mercury 
+	{ Tl,	{ kNA, 	kNA, 	kNA, 	150,	kNA, 	88.5,	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Thallium 
+	{ Pb,	{ kNA, 	kNA, 	kNA, 	kNA, 	119,	kNA, 	77.5,	kNA, 	kNA, 	kNA, 	kNA } }, 	// Lead 
+	{ Bi,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	103,	kNA, 	76,		kNA, 	kNA, 	kNA } }, 	// Bismuth 
+	{ Po,	{ kNA, 	223,	kNA, 	kNA, 	kNA, 	kNA, 	94,		kNA, 	67,		kNA, 	kNA } }, 	// Polonium 
+	{ At,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	62,		kNA } }, 	// Astatine 
+	{ Fr,	{ kNA, 	kNA, 	kNA, 	180,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Francium 
+	{ Ra,	{ kNA, 	kNA, 	kNA, 	kNA, 	148,	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } }, 	// Radium 
+	{ Ac,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	106.5, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA } },		// Actinium 
+	{ Th,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	94,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Thorium 
+	{ Pa,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	104,	90,		78,		kNA, 	kNA, 	kNA } }, 	// Protactinium 
+	{ U,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	102.5,	89,		76,		73,		kNA, 	kNA } }, 	// Uranium 
+	{ Np,	{ kNA, 	kNA, 	kNA, 	kNA, 	110,	101,	87,		75,		72,		71,		kNA } }, 	// Neptunium 
+	{ Pu,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	100,	86,		74,		71,		kNA, 	kNA } }, 	// Plutonium 
+	{ Am,	{ kNA, 	kNA, 	kNA, 	kNA, 	126,	97.5,	85,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Americium 
+	{ Cm,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	97,		85,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Curium 
+	{ Bk,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	96,		83,		kNA, 	kNA, 	kNA, 	kNA } }, 	// Berkelium 
+	{ Cf,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	95,		82.1,	kNA, 	kNA, 	kNA, 	kNA } }, 	// Californium 
+	{ Es,	{ kNA, 	kNA, 	kNA, 	kNA, 	kNA, 	83.5,	kNA, 	kNA,	kNA,	kNA, 	kNA } }, 	// Einsteinium 	
+};
+
+
 // --------------------------------------------------------------------
 // The coefficients from Waasmaier & Kirfel (1995), Acta Cryst. A51, 416-431.
 
@@ -886,5 +1096,45 @@ auto AtomTypeTraits::elsf() const -> const SFData&
 
 	throw std::runtime_error("No scattering factor found for " + name());
 }
-	
+
+// ionic radii
+
+float AtomTypeTraits::crystal_ionic_radius(int charge) const
+{
+	float result = data::kNA;
+
+	if (charge >= -3 and charge <= 8)
+	{
+		for (auto &r : data::kCrystalIonicRadii)
+		{
+			if (r.type != mInfo->type)
+				continue;
+			
+			result = r.radii[charge < 0 ? charge + 3 : charge + 2] / 100.0f;
+			break;
+		}
+	}
+
+	return result;
+}
+
+float AtomTypeTraits::effective_ionic_radius(int charge) const
+{
+	float result = data::kNA;
+
+	if (charge >= -3 and charge <= 8)
+	{
+		for (auto &r : data::kEffectiveIonicRadii)
+		{
+			if (r.type != mInfo->type)
+				continue;
+			
+			result = r.radii[charge < 0 ? charge + 3 : charge + 2] / 100.0f;
+			break;
+		}
+	}
+
+	return result;
+}
+
 }

src/BondMap.cpp

@@ -180,7 +180,10 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &a
 
 	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())
@@ -230,7 +233,7 @@ 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;
@@ -252,54 +255,70 @@ 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;
+	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
 		{
 			lastAsymID = asymID;
 			lastSeqID = seqID;
+			lastAuthSeqID = authSeqID;
 			continue;
 		}
 
-		auto c = atomMapByAsymSeqAndAtom[make_tuple(asymID, lastSeqID, "C")];
-		auto n = atomMapByAsymSeqAndAtom[make_tuple(asymID, seqID, "N")];
+		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);
 
-		if (not(c.empty() or n.empty()))
 			bindAtoms(c, n);
+		}
+		// if (not(c.empty() or n.empty()))
 
 		lastSeqID = seqID;
+		lastAuthSeqID = authSeqID;
 	}
 
 	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) =
+		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_label_seq_id", "ptnr2_label_seq_id",
+				"ptnr1_auth_seq_id", "ptnr2_auth_seq_id");
 
-		std::string a = atomMapByAsymSeqAndAtom[make_tuple(asym1, seqId1, atomId1)];
-		std::string b = atomMapByAsymSeqAndAtom[make_tuple(asym2, seqId2, atomId2)];
+		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);
 
-		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
@@ -308,7 +327,7 @@ BondMap::BondMap(const Structure &p)
 	{
 		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;
 		}
@@ -370,15 +389,12 @@ 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)
 			{

src/Cif++.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
@@ -30,10 +30,10 @@
 #include <numeric>
 #include <regex>
 #include <set>
+#include <shared_mutex>
 #include <stack>
 #include <tuple>
 #include <unordered_map>
-#include <shared_mutex>
 
 #include <filesystem>
 
@@ -570,36 +570,6 @@ void Datablock::write(std::ostream &os, const std::vector<std::string> &order)
 
 		cat.write(os);
 	}
-
-	//	// mmcif support, sort of. First write the 'entry' Category
-	//	// and if it exists, _AND_ we have a Validator, write out the
-	//	// auditConform record.
-	//
-	//	for (auto& cat: mCategories)
-	//	{
-	//		if (cat.name() == "entry")
-	//		{
-	//			cat.write(os);
-	//
-	//			if (mValidator != nullptr)
-	//			{
-	//				Category auditConform(*this, "audit_conform", nullptr);
-	//				auditConform.emplace({
-	//					{ "dict_name", mValidator->dictName() },
-	//					{ "dict_version", mValidator->dictVersion() }
-	//				});
-	//				auditConform.write(os);
-	//			}
-	//
-	//			break;
-	//		}
-	//	}
-	//
-	//	for (auto& cat: mCategories)
-	//	{
-	//		if (cat.name() != "entry" and cat.name() != "audit_conform")
-	//			cat.write(os);
-	//	}
 }
 
 bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
@@ -636,51 +606,43 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
 	{
 		std::string nA = *catA_i;
 		ba::to_lower(nA);
-		
+
 		std::string nB = *catB_i;
 		ba::to_lower(nB);
-		
+
 		int d = nA.compare(nB);
 		if (d > 0)
 		{
-			auto cat = dbB.get(*catB_i);
-			
-			if (cat == nullptr)
-				missingA.push_back(*catB_i);
-
+			missingA.push_back(*catB_i);
 			++catB_i;
 		}
 		else if (d < 0)
 		{
-			auto cat = dbA.get(*catA_i);
-			
-			if (cat == nullptr)
-				missingB.push_back(*catA_i);
-
+			missingB.push_back(*catA_i);
 			++catA_i;
 		}
 		else
 			++catA_i, ++catB_i;
 	}
-	
+
 	while (catA_i != catA.end())
 		missingB.push_back(*catA_i++);
 
 	while (catB_i != catB.end())
 		missingA.push_back(*catB_i++);
 
-	if (not (missingA.empty() and missingB.empty()))
+	if (not(missingA.empty() and missingB.empty()))
 	{
 		if (cif::VERBOSE > 1)
 		{
 			std::cerr << "compare of datablocks failed" << std::endl;
 			if (not missingA.empty())
 				std::cerr << "Categories missing in A: " << ba::join(missingA, ", ") << std::endl
-					<< std::endl;
-		
+						  << std::endl;
+
 			if (not missingB.empty())
 				std::cerr << "Categories missing in B: " << ba::join(missingB, ", ") << std::endl
-					<< std::endl;
+						  << std::endl;
 
 			result = false;
 		}
@@ -706,7 +668,7 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
 			++catA_i;
 		else
 		{
-			if (not (*dbA.get(*catA_i) == *dbB.get(*catB_i)))
+			if (not(*dbA.get(*catA_i) == *dbB.get(*catB_i)))
 			{
 				if (cif::VERBOSE > 1)
 				{
@@ -724,10 +686,10 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
 	return result;
 }
 
-std::ostream& operator<<(std::ostream &os, const Datablock &data)
+std::ostream &operator<<(std::ostream &os, const Datablock &data)
 {
 	// whoohoo... this sucks!
-	const_cast<Datablock&>(data).write(os);
+	const_cast<Datablock &>(data).write(os);
 	return os;
 }
 
@@ -1162,7 +1124,7 @@ void CatIndex::reconstruct()
 		insert(r.mData);
 
 	// maybe reconstruction can be done quicker by using the following commented code.
-	// however, I've not had the time to think of a way to std::set the red/black flag correctly in that case.
+	// however, I've not had the time to think of a way to set the red/black flag correctly in that case.
 
 	//	std::vector<ItemRow*> rows;
 	//	transform(mCat.begin(), mCat.end(), backInserter(rows),
@@ -1254,82 +1216,15 @@ size_t CatIndex::size() const
 	return result;
 }
 
-//bool CatIndex::isValid() const
-//{
-//	bool result = true;
-//
-//	if (mRoot != nullptr)
-//	{
-//		uint32_t minBlack = numeric_limits<uint32_t>::max();
-//		uint32_t maxBlack = 0;
-//
-//		assert(not mRoot->mRed);
-//
-//		result = isValid(mRoot, false, 0, minBlack, maxBlack);
-//		assert(minBlack == maxBlack);
-//	}
-//
-//	return result;
-//}
-//
-//bool CatIndex::validate(entry* h, bool isParentRed, uint32_t blackDepth, uint32_t& minBlack, uint32_t& maxBlack) const
-//{
-//	bool result = true;
-//
-//	if (h->mRed)
-//		assert(not isParentRed);
-//	else
-//		++blackDepth;
-//
-//	if (isParentRed)
-//		assert(not h->mRed);
-//
-//	if (h->mLeft != nullptr and h->mRight != nullptr)
-//	{
-//		if (isRed(h->mLeft))
-//			assert(not isRed(h->mRight));
-//		if (isRed(h->mRight))
-//			assert(not isRed(h->mLeft));
-//	}
-//
-//	if (h->mLeft != nullptr)
-//	{
-//		assert(mComp(h->mLeft->mRow, h->mRow) < 0);
-//		validate(h->mLeft, h->mRed, blackDepth, minBlack, maxBlack);
-//	}
-//	else
-//	{
-//		if (minBlack > blackDepth)
-//			minBlack = blackDepth;
-//		if (maxBlack < blackDepth)
-//			maxBlack = blackDepth;
-//	}
-//
-//	if (h->mRight != nullptr)
-//	{
-//		assert(mComp(h->mRight->mRow, h->mRow) > 0);
-//		validate(h->mRight, h->mRight, blackDepth, minBlack, maxBlack);
-//	}
-//	else
-//	{
-//		if (minBlack > blackDepth)
-//			minBlack = blackDepth;
-//		if (maxBlack < blackDepth)
-//			maxBlack = blackDepth;
-//	}
-//}
-
 // --------------------------------------------------------------------
 
 RowSet::RowSet(Category &cat)
 	: mCat(&cat)
-// , mCond(nullptr)
 {
 }
 
 RowSet::RowSet(Category &cat, Condition &&cond)
 	: mCat(&cat)
-// , mCond(new Condition(std::forward<Condition>(cond)))
 {
 	cond.prepare(cat);
 
@@ -1343,21 +1238,17 @@ RowSet::RowSet(Category &cat, Condition &&cond)
 RowSet::RowSet(const RowSet &rhs)
 	: mCat(rhs.mCat)
 	, mItems(rhs.mItems)
-// , mCond(nullptr)
 {
 }
 
 RowSet::RowSet(RowSet &&rhs)
 	: mCat(rhs.mCat)
 	, mItems(std::move(rhs.mItems))
-// , mCond(rhs.mCond)
 {
-	// rhs.mCond = nullptr;
 }
 
 RowSet::~RowSet()
 {
-	// delete mCond;
 }
 
 RowSet &RowSet::operator=(const RowSet &rhs)
@@ -1469,7 +1360,7 @@ void Category::updateLinks()
 			auto childCat = mDb.get(link->mChildCategory);
 			if (childCat == nullptr)
 				continue;
-			mChildLinks.push_back({ childCat, link });
+			mChildLinks.push_back({childCat, link});
 		}
 
 		for (auto link : mValidator->getLinksForChild(mName))
@@ -1477,7 +1368,7 @@ void Category::updateLinks()
 			auto parentCat = mDb.get(link->mParentCategory);
 			if (parentCat == nullptr)
 				continue;
-			mParentLinks.push_back({ parentCat, link });
+			mParentLinks.push_back({parentCat, link});
 		}
 	}
 }
@@ -1497,7 +1388,7 @@ size_t Category::getColumnIndex(std::string_view name) const
 			break;
 	}
 
-	if (VERBOSE and result == mColumns.size() and mCatValidator != nullptr) // validate the name, if it is known at all (since it was not found)
+	if (VERBOSE > 0 and result == mColumns.size() and mCatValidator != nullptr) // validate the name, if it is known at all (since it was not found)
 	{
 		auto iv = mCatValidator->getValidatorForItem(name);
 		if (iv == nullptr)
@@ -1543,21 +1434,6 @@ size_t Category::addColumn(std::string_view name)
 	return result;
 }
 
-// RowSet Category::find(Condition&& cond)
-// {
-// 	RowSet result(*this);
-
-// 	cond.prepare(*this);
-
-// 	for (auto r: *this)
-// 	{
-// 		if (cond(*this, r))
-// 			result.push_back(r);
-// 	}
-
-// 	return result;
-// }
-
 void Category::reorderByIndex()
 {
 	if (mIndex != nullptr)
@@ -1601,11 +1477,13 @@ std::string Category::getUniqueID(std::function<std::string(int)> generator)
 	if (mCatValidator != nullptr and mCatValidator->mKeys.size() == 1)
 		key = mCatValidator->mKeys.front();
 
-	size_t nr = size();
+	// calling size() often is a waste of resources
+	if (mLastUniqueNr == 0)
+		mLastUniqueNr = size();
 
 	for (;;)
 	{
-		std::string result = generator(int(nr++));
+		std::string result = generator(static_cast<int>(mLastUniqueNr++));
 
 		if (exists(Key(key) == result))
 			continue;
@@ -1665,21 +1543,6 @@ Row Category::operator[](Condition &&cond)
 	return result;
 }
 
-// RowSet Category::find(Condition&& cond)
-// {
-// 	// return RowSet(*this, std::forward<Condition>(cond));
-// 	RowSet result(*this);
-
-// 	cond.prepare(*this);
-
-// 	for (auto r: *this)
-// 	{
-// 		if (cond(*this, r))
-// 			result.insert(result.end(), r);
-// 	}
-// 	return result;
-// }
-
 bool Category::exists(Condition &&cond) const
 {
 	bool result = false;
@@ -2037,11 +1900,16 @@ bool Category::hasParent(Row r, const Category &parentCat, const ValidateLink &l
 			if (mCatValidator->mMandatoryFields.count(name) and field.is_null())
 				cond = std::move(cond) and (Key(link.mParentKeys[ix]) == Empty());
 		}
-		else
+		else if (parentCat.mCatValidator->mMandatoryFields.count(link.mParentKeys[ix]))
 		{
 			const char *value = field.c_str();
 			cond = std::move(cond) and (Key(link.mParentKeys[ix]) == value);
 		}
+		else
+		{
+			const char *value = field.c_str();
+			cond = std::move(cond) and (Key(link.mParentKeys[ix]) == value or Key(link.mParentKeys[ix]) == Empty());
+		}
 	}
 
 	if (result and not cond.empty())
@@ -2338,9 +2206,27 @@ void Category::validateLinks() const
 		size_t missing = 0;
 		for (auto r : *this)
 			if (not hasParent(r, *parentCat, *link))
+			{
+				if (cif::VERBOSE > 1)
+				{
+					if (missing == 0)
+					{
+						std::cerr << "Links for " << link->mLinkGroupLabel << " are incomplete" << std::endl
+								  << "  These are the items in " << mName << " that don't have matching parent items in " << parentCat->mName << std::endl
+								  << std::endl;
+					}
+
+					for (auto k : link->mChildKeys)
+						std::cerr << k << ": " << r[k].as<std::string>() << std::endl;
+					
+					std::cerr << std::endl;
+				}
+
 				++missing;
 
-		if (missing)
+			}
+
+		if (missing and VERBOSE == 1)
 		{
 			std::cerr << "Links for " << link->mLinkGroupLabel << " are incomplete" << std::endl
 					  << "  There are " << missing << " items in " << mName << " that don't have matching parent items in " << parentCat->mName << std::endl;
@@ -2408,26 +2294,26 @@ std::set<size_t> Category::keyFieldsByIndex() const
 
 bool operator==(const Category &a, const Category &b)
 {
-	using namespace std::placeholders; 
+	using namespace std::placeholders;
 
 	bool result = true;
 
-//	set<std::string> tagsA(a.fields()), tagsB(b.fields());
-//	
-//	if (tagsA != tagsB)
-//		std::cout << "Unequal number of fields" << std::endl;
+	//	set<std::string> tagsA(a.fields()), tagsB(b.fields());
+	//
+	//	if (tagsA != tagsB)
+	//		std::cout << "Unequal number of fields" << std::endl;
 
-	auto& validator = a.getValidator();
+	auto &validator = a.getValidator();
 	auto catValidator = validator.getValidatorForCategory(a.name());
 	if (catValidator == nullptr)
 		throw std::runtime_error("missing cat validator");
-	
-	typedef std::function<int(const char*,const char*)> compType;
-	std::vector<std::tuple<std::string,compType>> tags;
+
+	typedef std::function<int(const char *, const char *)> compType;
+	std::vector<std::tuple<std::string, compType>> tags;
 	auto keys = catValidator->mKeys;
 	std::vector<size_t> keyIx;
-	
-	for (auto& tag: a.fields())
+
+	for (auto &tag : a.fields())
 	{
 		auto iv = catValidator->getValidatorForItem(tag);
 		if (iv == nullptr)
@@ -2436,24 +2322,25 @@ bool operator==(const Category &a, const Category &b)
 		if (tv == nullptr)
 			throw std::runtime_error("missing type validator");
 		tags.push_back(std::make_tuple(tag, std::bind(&cif::ValidateType::compare, tv, std::placeholders::_1, std::placeholders::_2)));
-		
-		auto pred = [tag](const std::string& s) -> bool { return cif::iequals(tag, s) == 0; };
+
+		auto pred = [tag](const std::string &s) -> bool
+		{ return cif::iequals(tag, s) == 0; };
 		if (find_if(keys.begin(), keys.end(), pred) == keys.end())
 			keyIx.push_back(tags.size() - 1);
 	}
-	
+
 	// a.reorderByIndex();
 	// b.reorderByIndex();
-	
-	auto rowEqual = [&](const cif::Row& ra, const cif::Row& rb)
+
+	auto rowEqual = [&](const cif::Row &ra, const cif::Row &rb)
 	{
 		int d = 0;
 
-		for (auto kix: keyIx)
+		for (auto kix : keyIx)
 		{
 			std::string tag;
 			compType compare;
-			
+
 			std::tie(tag, compare) = tags[kix];
 
 			d = compare(ra[tag].c_str(), rb[tag].c_str());
@@ -2465,7 +2352,7 @@ bool operator==(const Category &a, const Category &b)
 				break;
 			}
 		}
-		
+
 		return d == 0;
 	};
 
@@ -2483,9 +2370,9 @@ bool operator==(const Category &a, const Category &b)
 			else
 				return false;
 		}
-		
+
 		cif::Row ra = *ai, rb = *bi;
-		
+
 		if (not rowEqual(ra, rb))
 		{
 			if (cif::VERBOSE > 1)
@@ -2493,21 +2380,25 @@ bool operator==(const Category &a, const Category &b)
 			else
 				return false;
 		}
-		
+
 		std::vector<std::string> missingA, missingB, different;
-		
-		for (auto& tt: tags)
+
+		for (auto &tt : tags)
 		{
 			std::string tag;
 			compType compare;
-			
+
 			std::tie(tag, compare) = tt;
-			
+
 			// make it an option to compare unapplicable to empty or something
-			
-			const char* ta = ra[tag].c_str();	if (strcmp(ta, ".") == 0 or strcmp(ta, "?") == 0) ta = "";
-			const char* tb = rb[tag].c_str();	if (strcmp(tb, ".") == 0 or strcmp(tb, "?") == 0) tb = "";
-			
+
+			const char *ta = ra[tag].c_str();
+			if (strcmp(ta, ".") == 0 or strcmp(ta, "?") == 0)
+				ta = "";
+			const char *tb = rb[tag].c_str();
+			if (strcmp(tb, ".") == 0 or strcmp(tb, "?") == 0)
+				tb = "";
+
 			if (compare(ta, tb) != 0)
 			{
 				if (cif::VERBOSE > 1)
@@ -2519,7 +2410,7 @@ bool operator==(const Category &a, const Category &b)
 					return false;
 			}
 		}
-		
+
 		++ai;
 		++bi;
 	}
@@ -2527,24 +2418,12 @@ bool operator==(const Category &a, const Category &b)
 	return result;
 }
 
-// auto Category::iterator::operator++() -> iterator&
-// {
-// 	mCurrent = Row(mCurrent.data()->mNext);
-// 	return *this;
-// }
-
-// auto Category::const_iterator::operator++() -> const_iterator&
-// {
-// 	mCurrent = Row(mCurrent.data()->mNext);
-// 	return *this;
-// }
-
 namespace detail
 {
 
 	size_t writeValue(std::ostream &os, std::string value, size_t offset, size_t width)
 	{
-		if (value.find('\n') != std::string::npos or width == 0 or value.length() >= 132) // write as text field
+		if (value.find('\n') != std::string::npos or width == 0 or value.length() > 132) // write as text field
 		{
 			ba::replace_all(value, "\n;", "\n\\;");
 
@@ -2647,7 +2526,7 @@ void Category::write(std::ostream &os, const std::vector<size_t> &order, bool in
 					if (not isUnquotedString(v->mText))
 						l += 2;
 
-					if (l >= 132)
+					if (l > 132)
 						continue;
 
 					if (columnWidths[v->mColumnIndex] < l + 1)
@@ -2683,7 +2562,7 @@ void Category::write(std::ostream &os, const std::vector<size_t> &order, bool in
 				if (l < w)
 					l = w;
 
-				if (offset + l >= 132 and offset > 0)
+				if (offset + l > 132 and offset > 0)
 				{
 					os << std::endl;
 					offset = 0;
@@ -2691,7 +2570,7 @@ void Category::write(std::ostream &os, const std::vector<size_t> &order, bool in
 
 				offset = detail::writeValue(os, s, offset, w);
 
-				if (offset >= 132)
+				if (offset > 132)
 				{
 					os << std::endl;
 					offset = 0;
@@ -2900,7 +2779,7 @@ void Category::update_value(RowSet &&rows, const std::string &tag, const std::st
 				}
 
 				// cannot update this...
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Cannot update child " << childCat->mName << "." << childTag << " with value " << value << std::endl;
 			}
 
@@ -3009,21 +2888,22 @@ void Row::assign(const Item &value, bool skipUpdateLinked)
 	assign(value.name(), value.value(), skipUpdateLinked);
 }
 
-void Row::assign(std::string_view name, const std::string &value, bool skipUpdateLinked)
+void Row::assign(std::string_view name, const std::string &value, bool skipUpdateLinked, bool validate)
 {
 	try
 	{
 		auto cat = mData->mCategory;
-		assign(cat->addColumn(name), value, skipUpdateLinked);
+		assign(cat->addColumn(name), value, skipUpdateLinked, validate);
 	}
 	catch (const std::exception &ex)
 	{
-		std::cerr << "Could not assign value '" << value << "' to column _" << mData->mCategory->name() << '.' << name << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Could not assign value '" << value << "' to column _" << mData->mCategory->name() << '.' << name << std::endl;
 		throw;
 	}
 }
 
-void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
+void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked, bool validate)
 {
 	if (mData == nullptr)
 		throw std::logic_error("invalid Row, no data assigning value '" + value + "' to column with index " + std::to_string(column));
@@ -3049,7 +2929,7 @@ void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
 	std::string oldStrValue = oldValue ? oldValue : "";
 
 	// check the value
-	if (col.mValidator)
+	if (col.mValidator and validate)
 		(*col.mValidator)(value);
 
 	// If the field is part of the Key for this Category, remove it from the index
@@ -3181,7 +3061,7 @@ void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
 			auto rows_n = childCat->find(std::move(cond_n));
 			if (not rows_n.empty())
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Will not rename in child category since there are already rows that link to the parent" << std::endl;
 
 				continue;
@@ -3387,7 +3267,7 @@ void Row::swap(size_t cix, ItemRow *a, ItemRow *b)
 					}
 					else
 					{
-						if (VERBOSE)
+						if (VERBOSE > 0)
 							std::cerr << "In " << childCat->mName << " changing " << linkChildColName << ": " << r[linkChildColName].as<std::string>() << " => " << (i ? i->mText : "") << std::endl;
 						r[linkChildColName] = i ? i->mText : "";
 					}
@@ -3496,11 +3376,17 @@ File::File(const std::filesystem::path &path, bool validate)
 	}
 	catch (const std::exception &ex)
 	{
-		std::cerr << "Error while loading file " << path << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Error while loading file " << path << std::endl;
 		throw;
 	}
 }
 
+File::File(const char *data, std::size_t length)
+{
+	load(data, length);
+}
+
 File::File(File &&rhs)
 	: mHead(nullptr)
 	, mValidator(nullptr)
@@ -3564,7 +3450,8 @@ void File::load(const std::filesystem::path &p)
 	}
 	catch (const std::exception &ex)
 	{
-		std::cerr << "Error loading file " << path << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Error loading file " << path << std::endl;
 		throw;
 	}
 }
@@ -3616,6 +3503,25 @@ void File::load(std::istream &is, const std::string &datablock)
 	}
 }
 
+void File::load(const char *data, std::size_t length)
+{
+	bool gzipped = length > 2 and data[0] == static_cast<char>(0x1f) and data[1] == static_cast<char>(0x8b);
+
+	struct membuf : public std::streambuf
+	{
+		membuf(char *data, size_t length) { this->setg(data, data, data + length); }
+	} buffer(const_cast<char *>(data), length);
+
+	std::istream is(&buffer);
+
+	io::filtering_stream<io::input> in;
+	if (gzipped)
+		in.push(io::gzip_decompressor());
+	in.push(is);
+
+	load(is);
+}
+
 void File::save(std::ostream &os)
 {
 	Datablock *e = mHead;
@@ -3656,11 +3562,26 @@ Datablock &File::operator[](std::string_view name)
 	return *result;
 }
 
+Datablock &File::front()
+{
+	assert(mHead);
+	return *mHead;
+}
+
+Datablock &File::back()
+{
+	assert(mHead);
+	auto *block = mHead;
+	while (block->mNext != nullptr)
+		block = block->mNext;
+	return *block;
+}
+
 bool File::isValid()
 {
 	if (mValidator == nullptr)
 	{
-		if (VERBOSE)
+		if (VERBOSE > 0)
 			std::cerr << "No dictionary loaded explicitly, loading default" << std::endl;
 
 		loadDictionary();

src/Cif2PDB.cpp

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

src/CifParser.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
@@ -42,7 +42,7 @@ namespace cif
 const uint32_t kMaxLineLength = 132;
 
 const uint8_t kCharTraitsTable[128] = {
-//	0	1	2	3	4	5	6	7	8	9	a	b	c	d	e	f
+	//	0	1	2	3	4	5	6	7	8	9	a	b	c	d	e	f
 	14, 15, 14, 14, 14, 15, 15, 14, 15, 15, 15, 15, 15, 15, 15, 15, //	2
 	15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 10, 15, 15, 15, 15, //	3
 	15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, //	4
@@ -82,6 +82,29 @@ const char *SacParser::kValueName[] = {
 
 // --------------------------------------------------------------------
 
+bool isUnquotedString(const char *s)
+{
+	auto ss = s;
+
+	bool result = isOrdinary(*s++);
+	while (result and *s != 0)
+	{
+		result = isNonBlank(*s);
+		++s;
+	}
+
+	// but be careful it does not contain e.g. stop_
+	if (result)
+	{
+		static const std::regex reservedRx(R"((^(?:data|save)|.*(?:loop|stop|global))_.+)", std::regex_constants::icase);
+		result = not std::regex_match(ss, reservedRx);
+	}
+
+	return result;
+}
+
+// --------------------------------------------------------------------
+
 SacParser::SacParser(std::istream &is, bool init)
 	: mData(is)
 {
@@ -151,23 +174,26 @@ void SacParser::retract()
 	mTokenValue.pop_back();
 }
 
-void SacParser::restart()
+
+int SacParser::restart(int start)
 {
+	int result = 0;
+
 	while (not mTokenValue.empty())
 		retract();
 
-	switch (mStart)
+	switch (start)
 	{
 		case eStateStart:
-			mState = mStart = eStateFloat;
+			result = eStateFloat;
 			break;
 
 		case eStateFloat:
-			mState = mStart = eStateInt;
+			result = eStateInt;
 			break;
 
 		case eStateInt:
-			mState = mStart = eStateValue;
+			result = eStateValue;
 			break;
 
 		default:
@@ -175,6 +201,8 @@ void SacParser::restart()
 	}
 
 	mBol = false;
+
+	return result;
 }
 
 void SacParser::match(SacParser::CIFToken t)
@@ -191,7 +219,7 @@ SacParser::CIFToken SacParser::getNextToken()
 
 	CIFToken result = eCIFTokenUnknown;
 	int quoteChar = 0;
-	mState = mStart = eStateStart;
+	int state = eStateStart, start = eStateStart;
 	mBol = false;
 
 	mTokenValue.clear();
@@ -201,7 +229,7 @@ SacParser::CIFToken SacParser::getNextToken()
 	{
 		auto ch = getNextChar();
 
-		switch (mState)
+		switch (state)
 		{
 			case eStateStart:
 				if (ch == kEOF)
@@ -209,27 +237,23 @@ SacParser::CIFToken SacParser::getNextToken()
 				else if (ch == '\n')
 				{
 					mBol = true;
-					mState = eStateWhite;
+					state = eStateWhite;
 				}
 				else if (ch == ' ' or ch == '\t')
-					mState = eStateWhite;
+					state = eStateWhite;
 				else if (ch == '#')
-					mState = eStateComment;
-				else if (ch == '.')
-					mState = eStateDot;
+					state = eStateComment;
 				else if (ch == '_')
-					mState = eStateTag;
+					state = eStateTag;
 				else if (ch == ';' and mBol)
-					mState = eStateTextField;
+					state = eStateTextField;
 				else if (ch == '\'' or ch == '"')
 				{
 					quoteChar = ch;
-					mState = eStateQuotedString;
+					state = eStateQuotedString;
 				}
-				else if (ch == '?')
-					mState = eStateQuestionMark;
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateWhite:
@@ -237,7 +261,7 @@ SacParser::CIFToken SacParser::getNextToken()
 					result = eCIFTokenEOF;
 				else if (not isspace(ch))
 				{
-					mState = eStateStart;
+					state = eStateStart;
 					retract();
 					mTokenValue.clear();
 				}
@@ -248,7 +272,7 @@ SacParser::CIFToken SacParser::getNextToken()
 			case eStateComment:
 				if (ch == '\n')
 				{
-					mState = eStateStart;
+					state = eStateStart;
 					mBol = true;
 					mTokenValue.clear();
 				}
@@ -258,34 +282,9 @@ SacParser::CIFToken SacParser::getNextToken()
 					error("invalid character in comment");
 				break;
 
-			case eStateQuestionMark:
-				if (isNonBlank(ch))
-					mState = eStateValue;
-				else
-				{
-					retract();
-					result = eCIFTokenValue;
-					mTokenValue.clear();
-					mTokenType = eCIFValueUnknown;
-				}
-				break;
-
-			case eStateDot:
-				if (isdigit(ch))
-					mState = eStateFloat + 2;
-				else if (isspace(ch))
-				{
-					retract();
-					result = eCIFTokenValue;
-					mTokenType = eCIFValueInapplicable;
-				}
-				else
-					mState = eStateValue;
-				break;
-
 			case eStateTextField:
 				if (ch == '\n')
-					mState = eStateTextField + 1;
+					state = eStateTextField + 1;
 				else if (ch == kEOF)
 					error("unterminated textfield");
 				else if (not isAnyPrint(ch))
@@ -295,7 +294,7 @@ SacParser::CIFToken SacParser::getNextToken()
 
 			case eStateTextField + 1:
 				if (isTextLead(ch) or ch == ' ' or ch == '\t')
-					mState = eStateTextField;
+					state = eStateTextField;
 				else if (ch == ';')
 				{
 					assert(mTokenValue.length() >= 2);
@@ -313,9 +312,10 @@ SacParser::CIFToken SacParser::getNextToken()
 				if (ch == kEOF)
 					error("unterminated quoted string");
 				else if (ch == quoteChar)
-					mState = eStateQuotedStringQuote;
+					state = eStateQuotedStringQuote;
 				else if (not isAnyPrint(ch))
-					error("invalid character in quoted string");
+					std::cerr << "invalid character in quoted string '" << std::string({static_cast<char>(ch)}) << "' (" << ch << ") line: " << mLineNr << std::endl;
+					// error("invalid character in quoted string");
 				break;
 
 			case eStateQuotedStringQuote:
@@ -331,7 +331,7 @@ SacParser::CIFToken SacParser::getNextToken()
 				else if (ch == quoteChar)
 					;
 				else if (isAnyPrint(ch))
-					mState = eStateQuotedString;
+					state = eStateQuotedString;
 				else if (ch == kEOF)
 					error("unterminated quoted string");
 				else
@@ -349,12 +349,12 @@ SacParser::CIFToken SacParser::getNextToken()
 			case eStateFloat:
 				if (ch == '+' or ch == '-')
 				{
-					mState = eStateFloat + 1;
+					state = eStateFloat + 1;
 				}
 				else if (isdigit(ch))
-					mState = eStateFloat + 1;
+					state = eStateFloat + 1;
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateFloat + 1:
@@ -362,9 +362,9 @@ SacParser::CIFToken SacParser::getNextToken()
 				//					mState = eStateNumericSuffix;
 				//				else
 				if (ch == '.')
-					mState = eStateFloat + 2;
+					state = eStateFloat + 2;
 				else if (tolower(ch) == 'e')
-					mState = eStateFloat + 3;
+					state = eStateFloat + 3;
 				else if (isWhite(ch) or ch == kEOF)
 				{
 					retract();
@@ -372,16 +372,13 @@ SacParser::CIFToken SacParser::getNextToken()
 					mTokenType = eCIFValueInt;
 				}
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			// parsed '.'
 			case eStateFloat + 2:
-				//				if (ch == '(')	// numeric???
-				//					mState = eStateNumericSuffix;
-				//				else
 				if (tolower(ch) == 'e')
-					mState = eStateFloat + 3;
+					state = eStateFloat + 3;
 				else if (isWhite(ch) or ch == kEOF)
 				{
 					retract();
@@ -389,30 +386,27 @@ SacParser::CIFToken SacParser::getNextToken()
 					mTokenType = eCIFValueFloat;
 				}
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			// parsed 'e'
 			case eStateFloat + 3:
 				if (ch == '-' or ch == '+')
-					mState = eStateFloat + 4;
+					state = eStateFloat + 4;
 				else if (isdigit(ch))
-					mState = eStateFloat + 5;
+					state = eStateFloat + 5;
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateFloat + 4:
 				if (isdigit(ch))
-					mState = eStateFloat + 5;
+					state = eStateFloat + 5;
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateFloat + 5:
-				//				if (ch == '(')
-				//					mState = eStateNumericSuffix;
-				//				else
 				if (isWhite(ch) or ch == kEOF)
 				{
 					retract();
@@ -420,14 +414,14 @@ SacParser::CIFToken SacParser::getNextToken()
 					mTokenType = eCIFValueFloat;
 				}
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateInt:
 				if (isdigit(ch) or ch == '+' or ch == '-')
-					mState = eStateInt + 1;
+					state = eStateInt + 1;
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
 			case eStateInt + 1:
@@ -438,35 +432,11 @@ SacParser::CIFToken SacParser::getNextToken()
 					mTokenType = eCIFValueInt;
 				}
 				else
-					restart();
+					state = start = restart(start);
 				break;
 
-				//			case eStateNumericSuffix:
-				//				if (isdigit(ch))
-				//					mState = eStateNumericSuffix + 1;
-				//				else
-				//					restart();
-				//				break;
-				//
-				//			case eStateNumericSuffix + 1:
-				//				if (ch == ')')
-				//				{
-				//					result = eCIFTokenValue;
-				//					mTokenType = eCIFValueNumeric;
-				//				}
-				//				else if (not isdigit(ch))
-				//					restart();
-				//				break;
-
 			case eStateValue:
-				if (isNonBlank(ch))
-					mState = eStateValue + 1;
-				else
-					error("invalid character at this position");
-				break;
-
-			case eStateValue + 1:
-				if (ch == '_') // first _, check for keywords
+				if (ch == '_')
 				{
 					std::string s = toLowerCopy(mTokenValue);
 
@@ -476,23 +446,40 @@ SacParser::CIFToken SacParser::getNextToken()
 						result = eCIFTokenSTOP;
 					else if (s == "loop_")
 						result = eCIFTokenLOOP;
-					else if (s == "data_" or s == "save_")
-						mState = eStateValue + 2;
+					else if (s == "data_")
+					{
+						state = eStateDATA;
+						continue;
+					}
+					else if (s == "save_")
+					{
+						state = eStateSAVE;
+						continue;
+					}
 				}
-				else if (not isNonBlank(ch))
+
+				if (result == eCIFTokenUnknown and not isNonBlank(ch))
 				{
 					retract();
 					result = eCIFTokenValue;
-					mTokenType = eCIFValueString;
+
+					if (mTokenValue == ".")
+						mTokenType = eCIFValueInapplicable;
+					else if (mTokenValue == "?")
+					{
+						mTokenType = eCIFValueUnknown;
+						mTokenValue.clear();
+					}
 				}
 				break;
 
-			case eStateValue + 2:
+			case eStateDATA:
+			case eStateSAVE:
 				if (not isNonBlank(ch))
 				{
 					retract();
 
-					if (tolower(mTokenValue[0]) == 'd')
+					if (state == eStateDATA)
 						result = eCIFTokenDATA;
 					else
 						result = eCIFTokenSAVE;
@@ -521,6 +508,7 @@ SacParser::CIFToken SacParser::getNextToken()
 	return result;
 }
 
+
 DatablockIndex SacParser::indexDatablocks()
 {
 	DatablockIndex index;
@@ -1220,7 +1208,7 @@ void DictParser::linkItems()
 	{
 		for (auto &iv : cv.mItemValidators)
 		{
-			if (iv.mType == nullptr)
+			if (iv.mType == nullptr and cif::VERBOSE >= 0)
 				std::cerr << "Missing item_type for " << iv.mTag << std::endl;
 		}
 	}
@@ -1255,7 +1243,8 @@ void DictParser::loadDictionary()
 	}
 	catch (const std::exception &)
 	{
-		std::cerr << "Error parsing dictionary" << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Error parsing dictionary" << std::endl;
 		throw;
 	}
 

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();
 }
 
 // --------------------------------------------------------------------
@@ -675,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);
 
@@ -1237,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;
 }
@@ -1284,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

@@ -354,7 +354,7 @@ 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;
 }
 

src/Compound.cpp

@@ -193,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;
 		}
@@ -403,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";
@@ -520,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;
@@ -611,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";
@@ -645,13 +645,13 @@ CompoundFactory::CompoundFactory()
 	auto ccd = cif::loadResource("components.cif");
 	if (ccd)
 		mImpl.reset(new CCDCompoundFactoryImpl(mImpl));
-	else if (cif::VERBOSE)
+	else if (cif::VERBOSE > 0)
 		std::cerr << "CCD components.cif file was not found" << std::endl;
 
 	const char *clibd_mon = getenv("CLIBD_MON");
 	if (clibd_mon != nullptr and fs::is_directory(clibd_mon))
 		mImpl.reset(new CCP4CompoundFactoryImpl(clibd_mon));
-	else if (cif::VERBOSE)
+	else if (cif::VERBOSE > 0)
 		std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;
 
 }
@@ -695,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;
 	}
 }
@@ -715,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/PDB2Cif.cpp

@@ -268,7 +268,8 @@ int PDBRecord::vI(int columnFirst, int columnLast)
 	}
 	catch (const std::exception &ex)
 	{
-		std::cerr << "Trying to parse '" << std::string(mValue + columnFirst - 7, mValue + columnLast - 7) << '\'' << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Trying to parse '" << std::string(mValue + columnFirst - 7, mValue + columnLast - 7) << '\'' << std::endl;
 		throw;
 	}
 
@@ -337,7 +338,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
 				}
 				else if (not isspace(ch))
 				{
-					if (cif::VERBOSE)
+					if (cif::VERBOSE > 0)
 						std::cerr << "skipping invalid character in SOURCE ID: " << ch << std::endl;
 				}
 				break;
@@ -354,7 +355,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
 			case eColon:
 				if (ch == ';')
 				{
-					if (cif::VERBOSE)
+					if (cif::VERBOSE > 0)
 						std::cerr << "Empty value for SOURCE: " << id << std::endl;
 					state = eStart;
 				}
@@ -418,7 +419,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
 			case eError:
 				if (ch == ';')
 				{
-					if (cif::VERBOSE)
+					if (cif::VERBOSE > 0)
 						std::cerr << "Skipping invalid header line: '" << std::string(start, mP) << std::endl;
 					state = eStart;
 				}
@@ -656,7 +657,7 @@ class PDBFileParser
 		int mSeqNum;
 		char mIcode;
 
-		int mDbSeqNum;
+		int mDbSeqNum = 0;
 		bool mSeen = false;
 		std::set<std::string> mAlts;
 
@@ -832,7 +833,7 @@ class PDBFileParser
 		if (not mChainSeq2AsymSeq.count(key))
 		{
 			ec = error::make_error_code(error::pdbErrors::residueNotFound);
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Residue " << chainID << resSeq << iCode << " could not be mapped" << std::endl;
 		}
 		else
@@ -929,7 +930,7 @@ class PDBFileParser
 		}
 		catch (const std::exception &ex)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << ex.what() << std::endl;
 			ec = error::make_error_code(error::pdbErrors::invalidDate);
 		}
@@ -940,7 +941,10 @@ class PDBFileParser
 	std::string pdb2cifDate(std::string s)
 	{
 		std::error_code ec;
-		return pdb2cifDate(s, ec);
+		auto result = pdb2cifDate(s, ec);
+		if (ec and cif::VERBOSE > 0)
+			std::cerr << "Invalid date(" << s << "): " << ec.message() << std::endl;
+		return result;
 	}
 
 	std::string pdb2cifAuth(std::string author)
@@ -1160,7 +1164,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
 			if (is.eof())
 				break;
 
-			std::cerr << "Line number " << lineNr << " is empty!" << std::endl;
+			if (cif::VERBOSE > 0)
+				std::cerr << "Line number " << lineNr << " is empty!" << std::endl;
 
 			getline(is, lookahead);
 			++lineNr;
@@ -1278,7 +1283,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
 				}
 				catch (const std::exception &ex)
 				{
-					std::cerr << "Dropping FORMUL line (" << (lineNr - 1) << ") with invalid component number '" << value.substr(1, 3) << '\'' << std::endl;
+					if (cif::VERBOSE >= 0)
+						std::cerr << "Dropping FORMUL line (" << (lineNr - 1) << ") with invalid component number '" << value.substr(1, 3) << '\'' << std::endl;
 					continue;
 					// throw_with_nested(std::runtime_error("Invalid component number '" + value.substr(1, 3) + '\''));
 				}
@@ -1305,7 +1311,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
 			}
 			catch (const std::exception &ex)
 			{
-				std::cerr << "Error parsing FORMUL at line " << lineNr << std::endl;
+				if (cif::VERBOSE >= 0)
+					std::cerr << "Error parsing FORMUL at line " << lineNr << std::endl;
 				throw;
 			}
 		}
@@ -1412,7 +1419,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
 
 	if (not dropped.empty())
 	{
-		std::cerr << "Dropped unsupported records: " << ba::join(dropped, ", ") << std::endl;
+		if (cif::VERBOSE >= 0)
+			std::cerr << "Dropped unsupported records: " << ba::join(dropped, ", ") << std::endl;
 	}
 
 	if (mData == nullptr)
@@ -1440,7 +1448,7 @@ void PDBFileParser::Match(const std::string &expected, bool throwIfMissing)
 	{
 		if (throwIfMissing)
 			throw std::runtime_error("Expected record " + expected + " but found " + mRec->mName);
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "Expected record " + expected + " but found " + mRec->mName << std::endl;
 	}
 }
@@ -1575,7 +1583,8 @@ void PDBFileParser::ParseTitle()
 
 			if (not iequals(key, "MOL_ID") and mCompounds.empty())
 			{
-				std::cerr << "Ignoring invalid COMPND record" << std::endl;
+				if (cif::VERBOSE >= 0)
+					std::cerr << "Ignoring invalid COMPND record" << std::endl;
 				break;
 			}
 
@@ -1625,7 +1634,7 @@ void PDBFileParser::ParseTitle()
 		//			auto colon = s.find(": ");
 		//			if (colon == std::string::npos)
 		//			{
-		//				if (cif::VERBOSE)
+		//				if (cif::VERBOSE > 0)
 		//					std::cerr << "invalid source field, missing colon (" << s << ')' << std::endl;
 		//				continue;
 		//			}
@@ -1713,7 +1722,7 @@ void PDBFileParser::ParseTitle()
 	// NUMMDL
 	if (mRec->is("NUMMDL"))
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "skipping unimplemented NUMMDL record" << std::endl;
 		GetNextRecord();
 	}
@@ -1816,7 +1825,7 @@ void PDBFileParser::ParseTitle()
 	// SPRSDE
 	if (mRec->is("SPRSDE"))
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "skipping unimplemented SPRSDE record" << std::endl;
 		GetNextRecord();
 	}
@@ -2265,7 +2274,7 @@ void PDBFileParser::ParseRemarks()
 									state = eMCP;
 								else if (subtopic == "CHIRAL CENTERS")
 									state = eChC;
-								else if (cif::VERBOSE)
+								else if (cif::VERBOSE > 0)
 									throw std::runtime_error("Unknown subtopic in REMARK 500: " + subtopic);
 
 								headerSeen = false;
@@ -2342,7 +2351,7 @@ void PDBFileParser::ParseRemarks()
 									}
 									catch (const std::exception &ex)
 									{
-										if (cif::VERBOSE)
+										if (cif::VERBOSE > 0)
 											std::cerr << "Dropping REMARK 500 at line " << mRec->mLineNr << " due to invalid symmetry operation" << std::endl;
 										continue;
 									}
@@ -2675,7 +2684,7 @@ void PDBFileParser::ParseRemarks()
 							case sStart:
 								if (s == "SITE")
 									state = sID;
-								else if (cif::VERBOSE)
+								else if (cif::VERBOSE > 0)
 									throw std::runtime_error("Invalid REMARK 800 record, expected SITE");
 								break;
 
@@ -2685,7 +2694,7 @@ void PDBFileParser::ParseRemarks()
 									id = m[1].str();
 									state = sEvidence;
 								}
-								else if (cif::VERBOSE)
+								else if (cif::VERBOSE > 0)
 									throw std::runtime_error("Invalid REMARK 800 record, expected SITE_IDENTIFIER");
 								break;
 
@@ -2695,7 +2704,7 @@ void PDBFileParser::ParseRemarks()
 									evidence = m[1].str();
 									state = sDesc;
 								}
-								else if (cif::VERBOSE)
+								else if (cif::VERBOSE > 0)
 									throw std::runtime_error("Invalid REMARK 800 record, expected SITE_IDENTIFIER");
 								break;
 
@@ -2906,7 +2915,7 @@ void PDBFileParser::ParseRemark200()
 		collectionDate = pdb2cifDate(rm200("DATE OF DATA COLLECTION", diffrnNr), ec);
 		if (ec)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << ec.message() << " for pdbx_collection_date" << std::endl;
 
 			// The date field can become truncated when multiple values are available
@@ -3025,7 +3034,7 @@ void PDBFileParser::ParseRemark200()
 	else if (inRM200({"HIGHEST RESOLUTION SHELL, RANGE LOW (A)", "COMPLETENESS FOR SHELL (%)",
 	                  "R MERGE FOR SHELL (I)", "R SYM FOR SHELL (I)", "<I/SIGMA(I)> FOR SHELL", "DATA REDUNDANCY IN SHELL"}))
 	{
-		if (cif::VERBOSE)
+		if (cif::VERBOSE > 0)
 			std::cerr << "Not writing reflns_shell record since d_res_high is missing" << std::endl;
 	}
 }
@@ -3595,7 +3604,7 @@ void PDBFileParser::ConstructEntities()
 			{
 				auto &r = chain.mResiduesSeen[lastResidueIndex + 1];
 
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 				{
 					std::cerr << "Detected residues that cannot be aligned to SEQRES" << std::endl
 							  << "First residue is " << chain.mDbref.chainID << ':' << r.mSeqNum << r.mIcode << std::endl;
@@ -4005,7 +4014,7 @@ void PDBFileParser::ConstructEntities()
 					tie(asym, labelSeq, std::ignore) = MapResidue(seqadv.chainID, seqadv.seqNum, seqadv.iCode, ec);
 					if (ec)
 					{
-						if (cif::VERBOSE)
+						if (cif::VERBOSE > 0)
 							std::cerr << "dropping unmatched SEQADV record" << std::endl;
 						continue;
 					}
@@ -4319,7 +4328,7 @@ void PDBFileParser::ConstructEntities()
 		tie(asymID, seq, std::ignore) = MapResidue(chainID, seqNum, iCode, ec);
 		if (ec) // no need to write a modres if it could not be found
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "dropping unmapped MODRES record" << std::endl;
 			continue;
 		}
@@ -4415,7 +4424,7 @@ void PDBFileParser::ConstructEntities()
 		tie(asymID, seqNr, isPolymer) = MapResidue(unobs.chain, unobs.seq, unobs.iCode, ec);
 		if (ec)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "error mapping unobserved residue" << std::endl;
 			continue;
 		}
@@ -4676,7 +4685,7 @@ void PDBFileParser::ParseSecondaryStructure()
 
 		if (ec)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Could not map residue for HELIX " << vI(8, 10) << std::endl;
 		}
 		else
@@ -4791,7 +4800,7 @@ void PDBFileParser::ParseSecondaryStructure()
 
 		if (ec)
 		{
-			if (cif::VERBOSE)
+			if (cif::VERBOSE > 0)
 				std::cerr << "Dropping SHEET record " << vI(8, 10) << std::endl;
 		}
 		else
@@ -4827,7 +4836,7 @@ void PDBFileParser::ParseSecondaryStructure()
 
 				if (ec)
 				{
-					if (cif::VERBOSE)
+					if (cif::VERBOSE > 0)
 						std::cerr << "skipping unmatched pdbx_struct_sheet_hbond record" << std::endl;
 				}
 				else
@@ -4927,7 +4936,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 
 			if (ec)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Dropping SSBOND " << vI(8, 10) << std::endl;
 				continue;
 			}
@@ -4948,7 +4957,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 			}
 			catch (const std::exception &ex)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Dropping SSBOND " << vI(8, 10) << " due to invalid symmetry operation" << std::endl;
 				continue;
 			}
@@ -4993,7 +5002,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 
 		if (mRec->is("LINK  ") or mRec->is("LINKR "))
 		{
-			if (cif::VERBOSE and mRec->is("LINKR "))
+			if (cif::VERBOSE > 0 and mRec->is("LINKR "))
 				std::cerr << "Accepting non-standard LINKR record, but ignoring extra information" << std::endl;
 
 			//	 1 -  6         Record name    "LINK  "
@@ -5046,7 +5055,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 
 			if (ec)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Dropping LINK record at line " << mRec->mLineNr << std::endl;
 				continue;
 			}
@@ -5062,7 +5071,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 				}
 				catch (const std::invalid_argument &)
 				{
-					if (cif::VERBOSE)
+					if (cif::VERBOSE > 0)
 						std::cerr << "Distance value '" << distance << "' is not a valid float in LINK record" << std::endl;
 					swap(ccp4LinkID, distance); // assume this is a ccp4_link_id... oh really?
 				}
@@ -5078,7 +5087,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 			}
 			catch (const std::exception &ex)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Dropping LINK record at line " << mRec->mLineNr << " due to invalid symmetry operation" << std::endl;
 				continue;
 			}
@@ -5149,7 +5158,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
 
 			if (ec)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Dropping CISPEP record at line " << mRec->mLineNr << std::endl;
 				continue;
 			}
@@ -5215,7 +5224,7 @@ void PDBFileParser::ParseMiscellaneousFeatures()
 
 			if (ec)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "skipping struct_site_gen record" << std::endl;
 			}
 			else
@@ -5518,7 +5527,7 @@ void PDBFileParser::ParseCoordinate(int modelNr)
 		{
 			if (groupPDB == "HETATM")
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Changing atom from HETATM to ATOM at line " << mRec->mLineNr << std::endl;
 				groupPDB = "ATOM";
 			}
@@ -5527,7 +5536,7 @@ void PDBFileParser::ParseCoordinate(int modelNr)
 		{
 			if (groupPDB == "ATOM")
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Changing atom from ATOM to HETATM at line " << mRec->mLineNr << std::endl;
 				groupPDB = "HETATM";
 			}
@@ -5698,7 +5707,8 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
 		}
 		catch (const std::exception &ex)
 		{
-			std::cerr << "Error parsing REMARK 3" << std::endl;
+			if (cif::VERBOSE >= 0)
+				std::cerr << "Error parsing REMARK 3" << std::endl;
 			throw;
 		}
 		//
@@ -5750,12 +5760,12 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
 
 				if ((symm1.empty() or symm1 == "1_555") and (symm2.empty() or symm2 == "1_555"))
 					distance = static_cast<float>(mmcif::Distance(mmcif::Point{x1, y1, z1}, mmcif::Point{x2, y2, z2}));
-				else if (cif::VERBOSE)
+				else if (cif::VERBOSE > 0)
 					std::cerr << "Cannot calculate distance for link since one of the atoms is in another dimension" << std::endl;
 			}
 			catch (std::exception &ex)
 			{
-				if (cif::VERBOSE)
+				if (cif::VERBOSE > 0)
 					std::cerr << "Error finding atom for LINK distance calculation: " << ex.what() << std::endl;
 			}
 
@@ -5764,10 +5774,13 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
 	}
 	catch (const std::exception &ex)
 	{
-		std::cerr << "Error parsing PDB";
-		if (mRec != nullptr)
-			std::cerr << " at line " << mRec->mLineNr;
-		std::cerr << std::endl;
+		if (cif::VERBOSE >= 0)
+		{
+			std::cerr << "Error parsing PDB";
+			if (mRec != nullptr)
+				std::cerr << " at line " << mRec->mLineNr;
+			std::cerr << std::endl;
+		}
 		throw;
 	}
 }
@@ -5947,7 +5960,7 @@ int PDBFileParser::PDBChain::AlignResToSeqRes()
 			switch (tb(x, y))
 			{
 				case -1:
-					//					if (cif::VERBOSE)
+					//					if (cif::VERBOSE > 0)
 					//						std::cerr << "A residue found in the ATOM records "
 					//							 << "(" << ry[y].mMonID << " @ " << mDbref.chainID << ":" << ry[y].mSeqNum
 					//							 	<<  ((ry[y].mIcode == ' ' or ry[y].mIcode == 0) ? "" : std::string{ ry[y].mIcode }) << ")"
@@ -6034,7 +6047,6 @@ void ReadPDBFile(std::istream &pdbFile, cif::File &cifFile)
 
 	p.Parse(pdbFile, cifFile);
 
-	if (not cifFile.isValid())
-		// throw std::runtime_error("Resulting mmCIF file is invalid");
+	if (not cifFile.isValid() and cif::VERBOSE >= 0)
 		std::cerr << "Resulting mmCIF file is not valid!" << std::endl;
 }

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

@@ -295,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:

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/rename-compound-test.cpp

@@ -32,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

@@ -78,8 +78,8 @@ BOOST_AUTO_TEST_CASE(create_nonpoly_1)
     cif::VERBOSE = 1;
 
 	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);
 
@@ -118,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
@@ -141,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
@@ -171,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");
+}