Backporting CMakeLists.txt file from version 3

.gitignore

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

CMakeLists.txt

@@ -25,7 +25,7 @@
 cmake_minimum_required(VERSION 3.16)
 
 # set the project name
-project(cifpp VERSION 4.2.0 LANGUAGES CXX)
+project(cifpp VERSION 2.0.5 LANGUAGES CXX)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
 
@@ -70,8 +70,6 @@ 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
@@ -394,7 +392,6 @@ if(CIFPP_BUILD_TESTS)
 		# pdb2cif
 		rename-compound
 		structure
-		sugar
 		unit)
 
 	foreach(CIFPP_TEST IN LISTS CIFPP_tests)

LICENSE

@@ -1,7 +1,6 @@
-BSD-2-Clause License
+SPDX-License-Identifier: BSD-2-Clause
 
 Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
-All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
@@ -21,4 +20,4 @@ ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 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.
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

changelog

@@ -1,49 +1,3 @@
-Version 4.2.0
-- Yet another rewrite of resource loading
-
-Version 4.1.1
-- Fall back to zero charge for scattering factors if the atom
-  was not found in the table.
-- Improve code to locate resources, failing less.
-
-Version 4.1.0
-- Some interface changes for mmcif::Atom
-
-Version 4.0.1
-- Added a bunch of const methods to Datablock and Category.
-- Changed PDB writing interface to accept Datablock instead of File.
-
-Version 4.0.0
-- getResidue in mmcif::Structure now requires both a
-  sequence ID and an auth sequence ID. As a result the code was cleaned
-  up considerably.
-
-Version 3.0.5
-- mmcif::Structure redesign. It is now a wrapper around a cif::Datablock.
-
-Version 3.0.4
-- Fix in mmCIF parser, now correctly handles the unquoted
-  string ??
-
-Version 3.0.3
-- Better configuration checks, for atomic e.g.
-- Fixed a problem introduced in refactoring mmcif::Atom
-- Version string creation
-
-Version 3.0.2
-- refactored mmcif::Atom for performance reasons
-
-Version 3.0.1
-- Fixed processing of proline restraints file from CCP4, proline
-  is a peptide, really.
-- Added code to facilitate DSSP
-
-Version 3.0.0
-- Replaced many strings in the API with string_view for
-  performance reasons.
-- Upgraded mmcif::Structure
-- various other small fixes
-
 Version 2.0.5
 - Backporting updated CMakeLists.txt file
 

cmake/FindFilesystem.cmake

@@ -59,10 +59,8 @@ 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()
 

data/ccd-subset.cif

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

include/cif++/AtomType.hpp

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

include/cif++/BondMap.hpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -26,9 +26,9 @@
 
 #pragma once
 
+#include <unordered_map>
 #include <filesystem>
 #include <stdexcept>
-#include <unordered_map>
 
 #include "cif++/Structure.hpp"
 
@@ -38,40 +38,39 @@ namespace mmcif
 class BondMapException : public std::runtime_error
 {
   public:
-	BondMapException(const std::string &msg)
-		: runtime_error(msg)
-	{
-	}
+	BondMapException(const std::string& msg)
+		: runtime_error(msg) {}
 };
 
 class BondMap
 {
   public:
-	BondMap(const Structure &p);
+	BondMap(const Structure& p);
+	
+	BondMap(const BondMap&) = delete;
+	BondMap& operator=(const BondMap&) = delete;
 
-	BondMap(const BondMap &) = delete;
-	BondMap &operator=(const BondMap &) = delete;
-
-	bool operator()(const Atom &a, const Atom &b) const
+	bool operator()(const Atom& a, const Atom& b) const
 	{
 		return isBonded(index.at(a.id()), index.at(b.id()));
 	}
 
-	bool is1_4(const Atom &a, const Atom &b) const
+	bool is1_4(const Atom& a, const Atom& b) const
 	{
 		uint32_t ixa = index.at(a.id());
 		uint32_t ixb = index.at(b.id());
-
+	
 		return bond_1_4.count(key(ixa, ixb));
 	}
-
+	
 	// links coming from the struct_conn records:
-	std::vector<std::string> linked(const Atom &a) const;
+	std::vector<std::string> linked(const Atom& a) const;
 
 	// This list of atomID's is comming from either CCD or the CCP4 dictionaries loaded
-	static std::vector<std::string> atomIDsForCompound(const std::string &compoundID);
-
+	static std::vector<std::string> atomIDsForCompound(const std::string& compoundID);
+	
   private:
+
 	bool isBonded(uint32_t ai, uint32_t bi) const
 	{
 		return bond.count(key(ai, bi)) != 0;
@@ -83,19 +82,20 @@ class BondMap
 			std::swap(a, b);
 		return static_cast<uint64_t>(a) | (static_cast<uint64_t>(b) << 32);
 	}
-
-	std::tuple<uint32_t, uint32_t> dekey(uint64_t k) const
+	
+	std::tuple<uint32_t,uint32_t> dekey(uint64_t k) const
 	{
 		return std::make_tuple(
 			static_cast<uint32_t>(k >> 32),
-			static_cast<uint32_t>(k));
+			static_cast<uint32_t>(k)
+		);
 	}
-
+	
 	uint32_t dim;
-	std::unordered_map<std::string, uint32_t> index;
+	std::unordered_map<std::string,uint32_t> index;
 	std::set<uint64_t> bond, bond_1_4;
 
-	std::map<std::string, std::set<std::string>> link;
+	std::map<std::string,std::set<std::string>> link;
 };
 
-} // namespace mmcif
+}

include/cif++/Cif++.hpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -26,82 +26,78 @@
 
 #pragma once
 
+#include <string>
+
 #include <array>
-#include <cstring>
 #include <functional>
-#include <iomanip>
 #include <iostream>
 #include <list>
 #include <optional>
 #include <regex>
 #include <set>
-#include <shared_mutex>
 #include <sstream>
-#include <string>
-
-#include <boost/format.hpp>
 
 #include "cif++/CifUtils.hpp"
 
 /*
-    Simple C++ interface to CIF files.
+	Simple C++ interface to CIF files.
+	
+	Assumptions: a file contains one or more datablocks modelled by the class datablock.
+	Each datablock contains categories. These map to the original tables used to fill
+	the mmCIF file. Each Category can contain multiple Items, the columns in the table.
+	
+	Values are stored as character strings internally.
+	
+	Synopsis:
+	
+	// create a cif file
+	
+	cif::datablock e("1MVE");
+	e.append(cif::Category{"_entry", { "id", "1MVE" } });
+	
+	cif::Category atomSite("atom_site");
+	size_t nr{};
+	for (auto& myAtom: atoms)
+	{
+		atomSite.push_back({
+			{ "group_PDB", "ATOM" },
+			{ "id", ++nr },
+			{ "type_symbol", myAtom.type.str() },
+			...
+		});
+	}
+	
+	e.append(move(atomSite));
+	
+	cif::File f;
+	f.append(e);
+	
+	ofstream os("1mve.cif");
+	f.write(os);
 
-    Assumptions: a file contains one or more datablocks modelled by the class datablock.
-    Each datablock contains categories. These map to the original tables used to fill
-    the mmCIF file. Each Category can contain multiple Items, the columns in the table.
+	// read
+	f.read(ifstream{"1mve.cif"});
+	
+	auto& e = f.firstDatablock();
+	
+	cout << "ID of datablock: " << e.id() << endl;
+	
+	auto& atoms = e["atom_site"];
+	for (auto& atom: atoms)
+	{
+		cout << atom["group_PDB"] << ", "
+			 << atom["id"] << ", "
+			 ...
 
-    Values are stored as character strings internally.
+		float x, y, z;
+		cif::tie(x, y, z) = atom.get("Cartn_x", "Cartn_y", "Cartn_z");
+		...
+	}
 
-    Synopsis:
-
-    // create a cif file
-
-    cif::datablock e("1MVE");
-    e.append(cif::Category{"_entry", { "id", "1MVE" } });
-
-    cif::Category atomSite("atom_site");
-    size_t nr{};
-    for (auto& myAtom: atoms)
-    {
-        atomSite.push_back({
-            { "group_PDB", "ATOM" },
-            { "id", ++nr },
-            { "type_symbol", myAtom.type.str() },
-            ...
-        });
-    }
-
-    e.append(move(atomSite));
-
-    cif::File f;
-    f.append(e);
-
-    ofstream os("1mve.cif");
-    f.write(os);
-
-    // read
-    f.read(ifstream{"1mve.cif"});
-
-    auto& e = f.firstDatablock();
-
-    cout << "ID of datablock: " << e.id() << endl;
-
-    auto& atoms = e["atom_site"];
-    for (auto& atom: atoms)
-    {
-        cout << atom["group_PDB"] << ", "
-             << atom["id"] << ", "
-             ...
-
-        float x, y, z;
-        cif::tie(x, y, z) = atom.get("Cartn_x", "Cartn_y", "Cartn_z");
-        ...
-    }
-
-    Another way of querying a Category is by using this construct:
-
-    auto cat& = e["atom_site"];
-    auto Rows = cat.find(Key("label_asym_id") == "A" and Key("label_seq_id") == 1);
+	Another way of querying a Category is by using this construct:
+	
+	auto cat& = e["atom_site"];
+	auto Rows = cat.find(Key("label_asym_id") == "A" and Key("label_seq_id") == 1);
 
 
 */
@@ -145,27 +141,14 @@ class Item
   public:
 	Item() {}
 
-	Item(std::string_view name, char value)
-		: mName(name)
-		, mValue({value})
-	{
-	}
-
-	template <typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
-	Item(std::string_view name, const T &value, const char *fmt)
-		: mName(name)
-		, mValue((boost::format(fmt) % value).str())
-	{
-	}
-
 	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
-	Item(const std::string_view name, const T &value)
+	Item(const std::string &name, const T &value)
 		: mName(name)
 		, mValue(std::to_string(value))
 	{
 	}
 
-	Item(const std::string_view name, const std::string_view value)
+	Item(const std::string &name, const std::string &value)
 		: mName(name)
 		, mValue(value)
 	{
@@ -176,7 +159,6 @@ class Item
 		, mValue(rhs.mValue)
 	{
 	}
-
 	Item(Item &&rhs) noexcept
 		: mName(std::move(rhs.mName))
 		, mValue(std::move(rhs.mValue))
@@ -239,7 +221,7 @@ class Datablock
 	using iterator = CategoryList::iterator;
 	using const_iterator = CategoryList::const_iterator;
 
-	Datablock(const std::string_view name);
+	Datablock(const std::string &name);
 	~Datablock();
 
 	Datablock(const Datablock &) = delete;
@@ -248,53 +230,44 @@ class Datablock
 	std::string getName() const { return mName; }
 	void setName(const std::string &n) { mName = n; }
 
+	std::string firstItem(const std::string &tag) const;
+
 	iterator begin() { return mCategories.begin(); }
 	iterator end() { return mCategories.end(); }
 
 	const_iterator begin() const { return mCategories.begin(); }
 	const_iterator end() const { return mCategories.end(); }
 
-	/// \brief Access to the category with name \a name, will create it if it doesn't exist.
-	Category &operator[](std::string_view name);
+	Category &operator[](const std::string &name);
 
-	// /// \brief Access to the category with name \a name, will throw if it doesn't exist.
-	// const Category &operator[](std::string_view name) const;
-
-	/// \brief Access to the category with name \a name, will return an empty category if is doesn't exist.
-	const Category &operator[](std::string_view name) const;
-
-	std::tuple<iterator, bool> emplace(std::string_view name);
+	std::tuple<iterator, bool> emplace(const std::string &name);
 
 	bool isValid();
 	void validateLinks() const;
 
-	void setValidator(const Validator *v);
+	void setValidator(Validator *v);
 
 	// this one only looks up a Category, returns nullptr if it does not exist
-	const Category *get(std::string_view name) const;
-	Category *get(std::string_view name);
+	const Category *get(const std::string &name) const;
+	Category *get(const std::string &name);
 
 	void getTagOrder(std::vector<std::string> &tags) const;
 	void write(std::ostream &os, const std::vector<std::string> &order);
 	void write(std::ostream &os);
 
 	// convenience function, add a line to the software category
-	void add_software(const std::string_view name, const std::string &classification,
+	void add_software(const std::string &name, const std::string &classification,
 		const std::string &versionNr, const std::string &versionDate);
 
 	friend bool operator==(const Datablock &lhs, const Datablock &rhs);
 
-	friend std::ostream &operator<<(std::ostream &os, const Datablock &data);
+	friend std::ostream& operator<<(std::ostream &os, const Datablock &data);
 
   private:
-	CategoryList mCategories; // LRU
-	mutable std::shared_mutex mLock;
+	std::list<Category> mCategories;
 	std::string mName;
-	const Validator *mValidator;
+	Validator *mValidator;
 	Datablock *mNext;
-
-	// for returning empty categories in the const operator[]
-	mutable std::unique_ptr<Category> mNullCategory;
 };
 
 // --------------------------------------------------------------------
@@ -377,14 +350,14 @@ namespace detail
 	  private:
 		friend class ::cif::Row;
 
-		ItemReference(std::string_view name, size_t column, Row &row)
+		ItemReference(const char *name, size_t column, Row &row)
 			: mName(name)
 			, mColumn(column)
 			, mRow(row)
 		{
 		}
 
-		ItemReference(std::string_view name, size_t column, const Row &row)
+		ItemReference(const char *name, size_t column, const Row &row)
 			: mName(name)
 			, mColumn(column)
 			, mRow(const_cast<Row &>(row))
@@ -392,7 +365,7 @@ namespace detail
 		{
 		}
 
-		std::string_view mName;
+		const char *mName;
 		size_t mColumn;
 		Row &mRow;
 		bool mConst = false;
@@ -773,16 +746,16 @@ class Row
 		return detail::ItemReference(itemTag, column, *this);
 	}
 
-	const detail::ItemReference operator[](std::string_view itemTag) const
+	const detail::ItemReference operator[](const std::string &itemTag) const
 	{
-		size_t column = ColumnForItemTag(itemTag);
-		return detail::ItemReference(itemTag, column, *this);
+		size_t column = ColumnForItemTag(itemTag.c_str());
+		return detail::ItemReference(itemTag.c_str(), column, *this);
 	}
 
-	detail::ItemReference operator[](std::string_view itemTag)
+	detail::ItemReference operator[](const std::string &itemTag)
 	{
-		size_t column = ColumnForItemTag(itemTag);
-		return detail::ItemReference(itemTag, column, *this);
+		size_t column = ColumnForItemTag(itemTag.c_str());
+		return detail::ItemReference(itemTag.c_str(), column, *this);
 	}
 
 	template <typename... Ts, size_t N>
@@ -803,7 +776,7 @@ class Row
 	}
 
 	void assign(const std::vector<Item> &values);
-	void assign(std::string_view name, const std::string &value, bool updateLinked, bool validate = true);
+	void assign(const std::string &name, const std::string &value, bool updateLinked);
 
 	bool operator==(const Row &rhs) const
 	{
@@ -825,12 +798,12 @@ class Row
 	friend std::ostream &operator<<(std::ostream &os, const Row &row);
 
   private:
-	void assign(size_t column, const std::string &value, bool updateLinked, bool validate = true);
+	void assign(size_t column, const std::string &value, bool updateLinked);
 	void assign(const Item &i, bool updateLinked);
 
 	static void swap(size_t column, ItemRow *a, ItemRow *b);
 
-	size_t ColumnForItemTag(std::string_view itemTag) const;
+	size_t ColumnForItemTag(const char *itemTag) const;
 
 	mutable ItemRow *mData;
 	uint32_t mLineNr = 0;
@@ -1194,16 +1167,13 @@ 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)
+	Key(const std::string &itemTag)
 		: mItemTag(itemTag)
 	{
 	}
-	explicit Key(const char *itemTag)
+	Key(const char *itemTag)
 		: mItemTag(itemTag)
 	{
 	}
@@ -1636,7 +1606,7 @@ class conditional_iterator_proxy
 	using iterator = conditional_iterator_impl;
 	using reference = typename iterator::reference;
 
-	template <typename... Ns>
+	template<typename... Ns>
 	conditional_iterator_proxy(CategoryType &cat, row_iterator pos, Condition &&cond, Ns... names);
 
 	conditional_iterator_proxy(conditional_iterator_proxy &&p);
@@ -1846,12 +1816,12 @@ class Category
 	friend class Row;
 	friend class detail::ItemReference;
 
-	Category(Datablock &db, const std::string_view name, const Validator *Validator);
+	Category(Datablock &db, const std::string &name, Validator *Validator);
 	Category(const Category &) = delete;
 	Category &operator=(const Category &) = delete;
 	~Category();
 
-	const std::string &name() const { return mName; }
+	const std::string name() const { return mName; }
 
 	using iterator = iterator_impl<Row>;
 	using const_iterator = iterator_impl<const Row>;
@@ -1873,9 +1843,6 @@ class Category
 	Row front() { return Row(mHead); }
 	Row back() { return Row(mTail); }
 
-	const Row front() const { return Row(mHead); }
-	const Row back() const { return Row(mTail); }
-
 	Row operator[](Condition &&cond);
 	const Row operator[](Condition &&cond) const;
 
@@ -1910,7 +1877,7 @@ class Category
 
 	conditional_iterator_proxy<const Category> find(const_iterator pos, Condition &&cond) const
 	{
-		return conditional_iterator_proxy<const Category>{const_cast<Category &>(*this), pos, std::forward<Condition>(cond)};
+		return conditional_iterator_proxy<const Category>{const_cast<Category&>(*this), pos, std::forward<Condition>(cond)};
 	}
 
 	template <typename... Ts, typename... Ns>
@@ -1931,14 +1898,14 @@ class Category
 	conditional_iterator_proxy<Category, Ts...> find(const_iterator pos, Condition &&cond, Ns... names)
 	{
 		static_assert(sizeof...(Ts) == sizeof...(Ns), "The number of column titles should be equal to the number of types to return");
-		return {*this, pos, std::forward<Condition>(cond), std::forward<Ns>(names)...};
+		return {*this, pos, std::forward<Condition>(cond), std::forward<Ns>(names)... };
 	}
 
 	template <typename... Ts, typename... Ns>
 	conditional_iterator_proxy<const Category, Ts...> find(const_iterator pos, Condition &&cond, Ns... names) const
 	{
 		static_assert(sizeof...(Ts) == sizeof...(Ns), "The number of column titles should be equal to the number of types to return");
-		return {*this, pos, std::forward<Condition>(cond), std::forward<Ns>(names)...};
+		return {*this, pos, std::forward<Condition>(cond), std::forward<Ns>(names)... };
 	}
 
 	// --------------------------------------------------------------------
@@ -1981,13 +1948,13 @@ class Category
 	}
 
 	template <typename T>
-	T find1(Condition &&cond, const char *column) const
+	T find1(Condition &&cond, const char* column) const
 	{
 		return find1<T>(cbegin(), std::forward<Condition>(cond), column);
 	}
 
 	template <typename T>
-	T find1(const_iterator pos, Condition &&cond, const char *column) const
+	T find1(const_iterator pos, Condition &&cond, const char* column) const
 	{
 		auto h = find<T>(pos, std::forward<Condition>(cond), column);
 
@@ -2097,7 +2064,7 @@ class Category
 
 	Datablock &db() { return mDb; }
 
-	void setValidator(const Validator *v);
+	void setValidator(Validator *v);
 
 	iset fields() const;
 	iset mandatoryFields() const;
@@ -2110,8 +2077,8 @@ class Category
 	void getTagOrder(std::vector<std::string> &tags) const;
 
 	// return index for known column, or the next available column index
-	size_t getColumnIndex(std::string_view name) const;
-	bool hasColumn(std::string_view name) const;
+	size_t getColumnIndex(const std::string &name) const;
+	bool hasColumn(const std::string &name) const;
 	const std::string &getColumnName(size_t columnIndex) const;
 	std::vector<std::string> getColumnNames() const;
 
@@ -2152,27 +2119,16 @@ class Category
 	void write(std::ostream &os, const std::vector<std::string> &order);
 	void write(std::ostream &os, const std::vector<size_t> &order, bool includeEmptyColumns);
 
-	size_t addColumn(std::string_view name);
-
-	struct Linked
-	{
-		Category *linked;
-		const ValidateLink *v;
-	};
-
-	void updateLinks();
+	size_t addColumn(const std::string &name);
 
 	Datablock &mDb;
 	std::string mName;
-	const Validator *mValidator;
+	Validator *mValidator;
 	const ValidateCategory *mCatValidator = nullptr;
 	std::vector<ItemColumn> mColumns;
 	ItemRow *mHead;
 	ItemRow *mTail;
-	size_t mLastUniqueNr = 0;
 	class CatIndex *mIndex;
-
-	std::vector<Linked> mParentLinks, mChildLinks;
 };
 
 // --------------------------------------------------------------------
@@ -2186,21 +2142,17 @@ 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;
 
-	virtual ~File();
+	~File();
 
-	virtual void load(const std::filesystem::path &p);
-	virtual void save(const std::filesystem::path &p);
+	void load(const std::filesystem::path &p);
+	void save(const std::filesystem::path &p);
 
-	virtual void load(std::istream &is);
-	virtual void save(std::ostream &os);
-
-	virtual void load(const char *data, std::size_t length);
+	void load(std::istream &is);
+	void save(std::ostream &os);
 
 	/// \brief Load only the data block \a datablock from the mmCIF file
 	void load(std::istream &is, const std::string &datablock);
@@ -2210,8 +2162,7 @@ class File
 
 	void loadDictionary();                 // load the default dictionary, that is mmcifDdl in this case
 	void loadDictionary(const char *dict); // load one of the compiled in dictionaries
-
-	void setValidator(const Validator *v);
+	void loadDictionary(std::istream &is); // load dictionary from input stream
 
 	bool isValid();
 	void validateLinks() const;
@@ -2232,14 +2183,8 @@ 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);
+	Datablock *get(const std::string &name) const;
+	Datablock &operator[](const std::string &name);
 
 	struct iterator
 	{
@@ -2275,17 +2220,16 @@ class File
 	iterator begin() const;
 	iterator end() const;
 
-	Datablock &front();
-	Datablock &back();
-
 	bool empty() const { return mHead == nullptr; }
 
 	const Validator &getValidator() const;
 	void getTagOrder(std::vector<std::string> &tags) const;
 
   private:
+	void setValidator(Validator *v);
+
 	Datablock *mHead;
-	const Validator *mValidator;
+	Validator *mValidator;
 };
 
 // --------------------------------------------------------------------

include/cif++/Cif2PDB.hpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -28,12 +28,12 @@
 
 #include "cif++/Cif++.hpp"
 
-void WritePDBFile(std::ostream &os, const cif::Datablock &data);
+void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile);
 
 /// \brief Just the HEADER, COMPND, SOURCE and AUTHOR lines
-void WritePDBHeaderLines(std::ostream &os, const cif::Datablock &data);
+void WritePDBHeaderLines(std::ostream& os, cif::File& cifFile);
 
-std::string GetPDBHEADERLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
-std::string GetPDBCOMPNDLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
-std::string GetPDBSOURCELine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
-std::string GetPDBAUTHORLine(const cif::Datablock &data, std::string::size_type truncate_at = 127);
+std::string GetPDBHEADERLine(cif::File& cifFile, std::string::size_type truncate_at = 127);
+std::string GetPDBCOMPNDLine(cif::File& cifFile, std::string::size_type truncate_at = 127);
+std::string GetPDBSOURCELine(cif::File& cifFile, std::string::size_type truncate_at = 127);
+std::string GetPDBAUTHORLine(cif::File& cifFile, std::string::size_type truncate_at = 127);

include/cif++/CifParser.hpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -28,8 +28,8 @@
 
 #include "cif++/Cif++.hpp"
 
-#include <map>
 #include <stack>
+#include <map>
 
 namespace cif
 {
@@ -39,7 +39,7 @@ namespace cif
 class CifParserError : public std::runtime_error
 {
   public:
-	CifParserError(uint32_t lineNr, const std::string &message);
+	CifParserError(uint32_t lineNr, const std::string& message);
 };
 
 // --------------------------------------------------------------------
@@ -48,8 +48,7 @@ extern const uint32_t kMaxLineLength;
 
 extern const uint8_t kCharTraitsTable[128];
 
-enum CharTraitsMask : uint8_t
-{
+enum CharTraitsMask: uint8_t {
 	kOrdinaryMask = 1 << 0,
 	kNonBlankMask = 1 << 1,
 	kTextLeadMask = 1 << 2,
@@ -76,17 +75,30 @@ inline bool isTextLead(int ch)
 	return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0;
 }
 
-inline bool isAnyPrint(int ch)
+inline bool isAnyPrint(int ch)	
 {
-	return ch == '\t' or
-	       (ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
+	return ch == '\t' or 
+		(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
 }
 
-bool isUnquotedString(const char *s);
+inline bool isUnquotedString(const char* s)
+{
+	bool result = isOrdinary(*s++);
+	while (result and *s != 0)
+	{
+		result = isNonBlank(*s);
+		++s;
+	}
+	return result;
+}
 
 // --------------------------------------------------------------------
 
-using DatablockIndex = std::map<std::string, std::size_t>;
+std::tuple<std::string,std::string> splitTagName(const std::string& tag);
+
+// --------------------------------------------------------------------
+
+using DatablockIndex = std::map<std::string,std::size_t>;
 
 // --------------------------------------------------------------------
 // sac Parser, analogous to SAX Parser (simple api for xml)
@@ -94,15 +106,15 @@ using DatablockIndex = std::map<std::string, std::size_t>;
 class SacParser
 {
   public:
-	SacParser(std::istream &is, bool init = true);
+	SacParser(std::istream& is, bool init = true);
 	virtual ~SacParser() {}
 
 	enum CIFToken
 	{
 		eCIFTokenUnknown,
-
+		
 		eCIFTokenEOF,
-
+	
 		eCIFTokenDATA,
 		eCIFTokenLOOP,
 		eCIFTokenGLOBAL,
@@ -112,7 +124,7 @@ class SacParser
 		eCIFTokenValue,
 	};
 
-	static const char *kTokenName[];
+	static const char* kTokenName[];
 
 	enum CIFValueType
 	{
@@ -125,39 +137,40 @@ class SacParser
 		eCIFValueUnknown
 	};
 
-	static const char *kValueName[];
-
+	static const char* kValueName[];
+	
 	int getNextChar();
 
 	void retract();
-	int restart(int start);
-
+	void restart();
+	
 	CIFToken getNextToken();
 	void match(CIFToken token);
 
-	bool parseSingleDatablock(const std::string &datablock);
+	bool parseSingleDatablock(const std::string& datablock);
 
 	DatablockIndex indexDatablocks();
-	bool parseSingleDatablock(const std::string &datablock, const DatablockIndex &index);
+	bool parseSingleDatablock(const std::string& datablock, const DatablockIndex &index);
 
 	void parseFile();
 	void parseGlobal();
 	void parseDataBlock();
 
 	virtual void parseSaveFrame();
-
+	
 	void parseDictionary();
-
-	void error(const std::string &msg);
-
+	
+	void error(const std::string& msg);
+	
 	// production methods, these are pure virtual here
-
-	virtual void produceDatablock(const std::string &name) = 0;
-	virtual void produceCategory(const std::string &name) = 0;
+	
+	virtual void produceDatablock(const std::string& name) = 0;
+	virtual void produceCategory(const std::string& name) = 0;
 	virtual void produceRow() = 0;
-	virtual void produceItem(const std::string &category, const std::string &item, const std::string &value) = 0;
+	virtual void produceItem(const std::string& category, const std::string& item, const std::string& value) = 0;
 
   protected:
+
 	enum State
 	{
 		eStateStart,
@@ -172,21 +185,21 @@ class SacParser
 		eStateTextField,
 		eStateFloat = 100,
 		eStateInt = 110,
-		eStateValue = 300,
-		eStateDATA,
-		eStateSAVE
+//		eStateNumericSuffix = 200,
+		eStateValue = 300
 	};
 
-	std::istream &mData;
+	std::istream&			mData;
 
 	// Parser state
-	bool mValidate;
-	uint32_t mLineNr;
-	bool mBol;
-	CIFToken mLookahead;
-	std::string mTokenValue;
-	CIFValueType mTokenType;
-	std::stack<int> mBuffer;
+	bool					mValidate;
+	uint32_t					mLineNr;
+	bool					mBol;
+	int						mState, mStart;
+	CIFToken				mLookahead;
+	std::string				mTokenValue;
+	CIFValueType			mTokenType;
+	std::stack<int>			mBuffer;
 };
 
 // --------------------------------------------------------------------
@@ -194,18 +207,18 @@ class SacParser
 class Parser : public SacParser
 {
   public:
-	Parser(std::istream &is, File &f, bool init = true);
+	Parser(std::istream& is, File& f, bool init = true);
 
-	virtual void produceDatablock(const std::string &name);
-	virtual void produceCategory(const std::string &name);
+	virtual void produceDatablock(const std::string& name);
+	virtual void produceCategory(const std::string& name);
 	virtual void produceRow();
-	virtual void produceItem(const std::string &category, const std::string &item, const std::string &value);
+	virtual void produceItem(const std::string& category, const std::string& item, const std::string& value);
 
   protected:
-	File &mFile;
-	Datablock *mDataBlock;
-	Datablock::iterator mCat;
-	Row mRow;
+	File&					mFile;
+	Datablock*				mDataBlock;
+	Datablock::iterator		mCat;
+	Row						mRow;
 };
 
 // --------------------------------------------------------------------
@@ -213,21 +226,23 @@ class Parser : public SacParser
 class DictParser : public Parser
 {
   public:
-	DictParser(Validator &validator, std::istream &is);
+
+	DictParser(Validator& validator, std::istream& is);
 	~DictParser();
-
+	
 	void loadDictionary();
-
+	
   private:
-	virtual void parseSaveFrame();
 
+	virtual void parseSaveFrame();
+	
 	bool collectItemTypes();
 	void linkItems();
 
-	Validator &mValidator;
-	File mFile;
-	struct DictParserDataImpl *mImpl;
-	bool mCollectedItemTypes = false;
+	Validator&						mValidator;
+	File							mFile;
+	struct DictParserDataImpl*		mImpl;
+	bool							mCollectedItemTypes = false;
 };
 
-} // namespace cif
+}

include/cif++/CifUtils.hpp

@@ -26,7 +26,11 @@
 
 #pragma once
 
+#include <cassert>
 #include <filesystem>
+#include <iostream>
+#include <list>
+#include <memory>
 #include <set>
 #include <vector>
 
@@ -63,8 +67,8 @@ std::string get_version_nr();
 // some basic utilities: Since we're using ASCII input only, we define for optimisation
 // our own case conversion routines.
 
-bool iequals(std::string_view a, std::string_view b);
-int icompare(std::string_view a, std::string_view b);
+bool iequals(const std::string &a, const std::string &b);
+int icompare(const std::string &a, const std::string &b);
 
 bool iequals(const char *a, const char *b);
 int icompare(const char *a, const char *b);
@@ -96,7 +100,7 @@ inline char tolower(int ch)
 
 // --------------------------------------------------------------------
 
-std::tuple<std::string, std::string> splitTagName(std::string_view tag);
+std::tuple<std::string, std::string> splitTagName(const std::string &tag);
 
 // --------------------------------------------------------------------
 // generate a cif name, mainly used to generate asym_id's

include/cif++/CifValidator.hpp

@@ -36,19 +36,18 @@
 
 namespace cif
 {
-
+	
 struct ValidateCategory;
-class ValidatorFactory;
 
 // --------------------------------------------------------------------
 
 class ValidationError : public std::exception
 {
   public:
-	ValidationError(const std::string &msg);
-	ValidationError(const std::string &cat, const std::string &item,
-		const std::string &msg);
-	const char *what() const noexcept { return mMsg.c_str(); }
+	ValidationError(const std::string& msg);
+	ValidationError(const std::string& cat, const std::string& item,
+		const std::string& msg);
+	const char* what() const noexcept		{ return mMsg.c_str(); }
 	std::string mMsg;
 };
 
@@ -56,60 +55,58 @@ class ValidationError : public std::exception
 
 enum class DDL_PrimitiveType
 {
-	Char,
-	UChar,
-	Numb
+	Char, UChar, Numb
 };
 
-DDL_PrimitiveType mapToPrimitiveType(std::string_view s);
+DDL_PrimitiveType mapToPrimitiveType(const std::string& s);
 
 struct ValidateType
 {
-	std::string mName;
-	DDL_PrimitiveType mPrimitiveType;
+	std::string				mName;
+	DDL_PrimitiveType		mPrimitiveType;
 	// std::regex				mRx;
-	boost::regex mRx;
+	boost::regex			mRx;
 
-	bool operator<(const ValidateType &rhs) const
+	bool operator<(const ValidateType& rhs) const
 	{
 		return icompare(mName, rhs.mName) < 0;
 	}
 
-	// compare values based on type
-	//	int compare(const std::string& a, const std::string& b) const
-	//	{
-	//		return compare(a.c_str(), b.c_str());
-	//	}
-
-	int compare(const char *a, const char *b) const;
+	// compare values based on type	
+//	int compare(const std::string& a, const std::string& b) const
+//	{
+//		return compare(a.c_str(), b.c_str());
+//	}
+	
+	int compare(const char* a, const char* b) const;
 };
 
 struct ValidateItem
 {
-	std::string mTag;
-	bool mMandatory;
-	const ValidateType *mType;
-	cif::iset mEnums;
-	std::string mDefault;
-	bool mDefaultIsNull;
-	ValidateCategory *mCategory = nullptr;
+	std::string				mTag;
+	bool					mMandatory;
+	const ValidateType*		mType;
+	cif::iset				mEnums;
+	std::string				mDefault;
+	bool					mDefaultIsNull;
+	ValidateCategory*		mCategory = nullptr;
 
 	// ItemLinked is used for non-key links
 	struct ItemLinked
 	{
-		ValidateItem *mParent;
-		std::string mParentItem;
-		std::string mChildItem;
+		ValidateItem*		mParent;
+		std::string			mParentItem;
+		std::string			mChildItem;
 	};
 
-	std::vector<ItemLinked> mLinked;
-
-	bool operator<(const ValidateItem &rhs) const
+	std::vector<ItemLinked>	mLinked;
+	
+	bool operator<(const ValidateItem& rhs) const
 	{
 		return icompare(mTag, rhs.mTag) < 0;
 	}
 
-	bool operator==(const ValidateItem &rhs) const
+	bool operator==(const ValidateItem& rhs) const
 	{
 		return iequals(mTag, rhs.mTag);
 	}
@@ -119,22 +116,22 @@ struct ValidateItem
 
 struct ValidateCategory
 {
-	std::string mName;
-	std::vector<std::string> mKeys;
-	cif::iset mGroups;
-	cif::iset mMandatoryFields;
-	std::set<ValidateItem> mItemValidators;
+	std::string					mName;
+	std::vector<std::string>	mKeys;
+	cif::iset					mGroups;
+	cif::iset					mMandatoryFields;
+	std::set<ValidateItem>		mItemValidators;
 
-	bool operator<(const ValidateCategory &rhs) const
+	bool operator<(const ValidateCategory& rhs) const
 	{
 		return icompare(mName, rhs.mName) < 0;
 	}
 
-	void addItemValidator(ValidateItem &&v);
-
-	const ValidateItem *getValidatorForItem(std::string_view tag) const;
-
-	const std::set<ValidateItem> &itemValidators() const
+	void addItemValidator(ValidateItem&& v);
+	
+	const ValidateItem* getValidatorForItem(std::string tag) const;
+	
+	const std::set<ValidateItem>& itemValidators() const
 	{
 		return mItemValidators;
 	}
@@ -142,12 +139,12 @@ struct ValidateCategory
 
 struct ValidateLink
 {
-	int mLinkGroupID;
-	std::string mParentCategory;
-	std::vector<std::string> mParentKeys;
-	std::string mChildCategory;
-	std::vector<std::string> mChildKeys;
-	std::string mLinkGroupLabel;
+	int							mLinkGroupID;
+	std::string					mParentCategory;
+	std::vector<std::string>	mParentKeys;
+	std::string					mChildCategory;
+	std::vector<std::string>	mChildKeys;
+	std::string					mLinkGroupLabel;
 };
 
 // --------------------------------------------------------------------
@@ -155,72 +152,47 @@ struct ValidateLink
 class Validator
 {
   public:
+	friend class DictParser;
 
-	Validator(std::string_view name, std::istream &is);
+	Validator();
 	~Validator();
 
-	Validator(const Validator &rhs) = delete;
-	Validator &operator=(const Validator &rhs) = delete;
+	Validator(const Validator& rhs) = delete;
+	Validator& operator=(const Validator& rhs) = delete;
+	
+	Validator(Validator&& rhs);
+	Validator& operator=(Validator&& rhs);
+	
+	void addTypeValidator(ValidateType&& v);
+	const ValidateType* getValidatorForType(std::string typeCode) const;
 
-	Validator(Validator &&rhs);
-	Validator &operator=(Validator &&rhs);
+	void addCategoryValidator(ValidateCategory&& v);
+	const ValidateCategory* getValidatorForCategory(std::string category) const;
 
-	friend class DictParser;
-	friend class ValidatorFactory;
+	void addLinkValidator(ValidateLink&& v);
+	std::vector<const ValidateLink*> getLinksForParent(const std::string& category) const;
+	std::vector<const ValidateLink*> getLinksForChild(const std::string& category) const;
 
-	void addTypeValidator(ValidateType &&v);
-	const ValidateType *getValidatorForType(std::string_view typeCode) const;
+	void reportError(const std::string& msg, bool fatal);
+	
+	std::string dictName() const					{ return mName; }
+	void dictName(const std::string& name)			{ mName = name; }
 
-	void addCategoryValidator(ValidateCategory &&v);
-	const ValidateCategory *getValidatorForCategory(std::string_view category) const;
-
-	void addLinkValidator(ValidateLink &&v);
-	std::vector<const ValidateLink *> getLinksForParent(std::string_view category) const;
-	std::vector<const ValidateLink *> getLinksForChild(std::string_view category) const;
-
-	void reportError(const std::string &msg, bool fatal) const;
-
-	std::string dictName() const { return mName; }
-	void dictName(const std::string &name) { mName = name; }
-
-	std::string dictVersion() const { return mVersion; }
-	void dictVersion(const std::string &version) { mVersion = version; }
+	std::string dictVersion() const				{ return mVersion; }
+	void dictVersion(const std::string& version)	{ mVersion = version; }
 
   private:
 
 	// name is fully qualified here:
-	ValidateItem *getValidatorForItem(std::string_view name) const;
+	ValidateItem* getValidatorForItem(std::string name) const;
 
-	std::string mName;
-	std::string mVersion;
-	bool mStrict = false;
-	//	std::set<uint32_t>			mSubCategories;
-	std::set<ValidateType> mTypeValidators;
-	std::set<ValidateCategory> mCategoryValidators;
-	std::vector<ValidateLink> mLinkValidators;
+	std::string					mName;
+	std::string					mVersion;
+	bool						mStrict = false;
+//	std::set<uint32_t>			mSubCategories;
+	std::set<ValidateType>		mTypeValidators;
+	std::set<ValidateCategory>	mCategoryValidators;
+	std::vector<ValidateLink>	mLinkValidators;
 };
 
-// --------------------------------------------------------------------
-
-class ValidatorFactory
-{
-  public:
-
-	static ValidatorFactory &instance()
-	{
-		return sInstance;
-	}
-
-	const Validator &operator[](std::string_view dictionary);
-
-  private:
-
-	static ValidatorFactory sInstance;
-
-	ValidatorFactory();
-
-	std::mutex mMutex;
-	std::list<Validator> mValidators;
-};
-
-} // namespace cif
+}

include/cif++/Compound.hpp

@@ -104,7 +104,6 @@ class Compound
 	std::string id() const { return mID; }
 	std::string name() const { return mName; }
 	std::string type() const { return mType; }
-	std::string group() const { return mGroup; }
 	std::string formula() const { return mFormula; }
 	float formulaWeight() const { return mFormulaWeight; }
 	int formalCharge() const { return mFormalCharge; }
@@ -131,12 +130,11 @@ class Compound
 	friend class CCP4CompoundFactoryImpl;
 
 	Compound(cif::Datablock &db);
-	Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group);
+	Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type);
 
 	std::string mID;
 	std::string mName;
 	std::string mType;
-	std::string mGroup;
 	std::string mFormula;
 	float mFormulaWeight = 0;
 	int mFormalCharge = 0;

include/cif++/Matrix.hpp

@@ -0,0 +1,391 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ * 
+ * Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
+ * Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ * 
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+// --------------------------------------------------------------------
+// uBlas compatible matrix types
+
+#pragma once
+
+#include <iostream>
+#include <vector>
+
+// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
+// element m i,j is mapped to [i * n + j] and thus storage is row major
+
+template <typename T>
+class MatrixBase
+{
+  public:
+	using value_type = T;
+
+	virtual ~MatrixBase() {}
+
+	virtual uint32_t dim_m() const = 0;
+	virtual uint32_t dim_n() const = 0;
+
+	virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
+	virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
+
+	MatrixBase &operator*=(const value_type &rhs);
+
+	MatrixBase &operator-=(const value_type &rhs);
+};
+
+template <typename T>
+MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
+{
+	for (uint32_t i = 0; i < dim_m(); ++i)
+	{
+		for (uint32_t j = 0; j < dim_n(); ++j)
+		{
+			operator()(i, j) *= rhs;
+		}
+	}
+
+	return *this;
+}
+
+template <typename T>
+MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
+{
+	for (uint32_t i = 0; i < dim_m(); ++i)
+	{
+		for (uint32_t j = 0; j < dim_n(); ++j)
+		{
+			operator()(i, j) -= rhs;
+		}
+	}
+
+	return *this;
+}
+
+template <typename T>
+std::ostream &operator<<(std::ostream &lhs, const MatrixBase<T> &rhs)
+{
+	lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
+	for (uint32_t i = 0; i < rhs.dim_m(); ++i)
+	{
+		lhs << '(';
+		for (uint32_t j = 0; j < rhs.dim_n(); ++j)
+		{
+			if (j > 0)
+				lhs << ',';
+			lhs << rhs(i, j);
+		}
+		lhs << ')';
+	}
+	lhs << ')';
+
+	return lhs;
+}
+
+template <typename T>
+class Matrix : public MatrixBase<T>
+{
+  public:
+	using value_type = T;
+
+	template <typename T2>
+	Matrix(const MatrixBase<T2> &m)
+		: m_m(m.dim_m())
+		, m_n(m.dim_n())
+	{
+		m_data = new value_type[m_m * m_n];
+		for (uint32_t i = 0; i < m_m; ++i)
+		{
+			for (uint32_t j = 0; j < m_n; ++j)
+				operator()(i, j) = m(i, j);
+		}
+	}
+
+	Matrix()
+		: m_data(nullptr)
+		, m_m(0)
+		, m_n(0)
+	{
+	}
+
+	Matrix(const Matrix &m)
+		: m_m(m.m_m)
+		, m_n(m.m_n)
+	{
+		m_data = new value_type[m_m * m_n];
+		std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
+	}
+
+	Matrix &operator=(const Matrix &m)
+	{
+		value_type *t = new value_type[m.m_m * m.m_n];
+		std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
+
+		delete[] m_data;
+		m_data = t;
+		m_m = m.m_m;
+		m_n = m.m_n;
+
+		return *this;
+	}
+
+	Matrix(uint32_t m, uint32_t n, T v = T())
+		: m_m(m)
+		, m_n(n)
+	{
+		m_data = new value_type[m_m * m_n];
+		std::fill(m_data, m_data + (m_m * m_n), v);
+	}
+
+	virtual ~Matrix()
+	{
+		delete[] m_data;
+	}
+
+	virtual uint32_t dim_m() const { return m_m; }
+	virtual uint32_t dim_n() const { return m_n; }
+
+	virtual value_type operator()(uint32_t i, uint32_t j) const
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+	virtual value_type &operator()(uint32_t i, uint32_t j)
+	{
+		assert(i < m_m);
+		assert(j < m_n);
+		return m_data[i * m_n + j];
+	}
+
+	template <typename Func>
+	void each(Func f)
+	{
+		for (uint32_t i = 0; i < m_m * m_n; ++i)
+			f(m_data[i]);
+	}
+
+	template <typename U>
+	Matrix &operator/=(U v)
+	{
+		for (uint32_t i = 0; i < m_m * m_n; ++i)
+			m_data[i] /= v;
+
+		return *this;
+	}
+
+  private:
+	value_type *m_data;
+	uint32_t m_m, m_n;
+};
+
+// --------------------------------------------------------------------
+
+template <typename T>
+class SymmetricMatrix : public MatrixBase<T>
+{
+  public:
+	typedef typename MatrixBase<T>::value_type value_type;
+
+	SymmetricMatrix(uint32_t n, T v = T())
+		: m_owner(true)
+		, m_n(n)
+	{
+		uint32_t N = (m_n * (m_n + 1)) / 2;
+		m_data = new value_type[N];
+		std::fill(m_data, m_data + N, v);
+	}
+
+	SymmetricMatrix(const T *data, uint32_t n)
+		: m_owner(false)
+		, m_data(const_cast<T *>(data))
+		, m_n(n)
+	{
+	}
+
+	virtual ~SymmetricMatrix()
+	{
+		if (m_owner)
+			delete[] m_data;
+	}
+
+	virtual uint32_t dim_m() const { return m_n; }
+	virtual uint32_t dim_n() const { return m_n; }
+
+	T operator()(uint32_t i, uint32_t j) const;
+	virtual T &operator()(uint32_t i, uint32_t j);
+
+	// erase two rows, add one at the end (for neighbour joining)
+	void erase_2(uint32_t i, uint32_t j);
+
+	template <typename Func>
+	void each(Func f)
+	{
+		uint32_t N = (m_n * (m_n + 1)) / 2;
+
+		for (uint32_t i = 0; i < N; ++i)
+			f(m_data[i]);
+	}
+
+	template <typename U>
+	SymmetricMatrix &operator/=(U v)
+	{
+		uint32_t N = (m_n * (m_n + 1)) / 2;
+
+		for (uint32_t i = 0; i < N; ++i)
+			m_data[i] /= v;
+
+		return *this;
+	}
+
+  private:
+	bool m_owner;
+	value_type *m_data;
+	uint32_t m_n;
+};
+
+template <typename T>
+inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
+{
+	return i < j
+	           ? m_data[(j * (j + 1)) / 2 + i]
+	           : m_data[(i * (i + 1)) / 2 + j];
+}
+
+template <typename T>
+inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
+{
+	if (i > j)
+		std::swap(i, j);
+	assert(j < m_n);
+	return m_data[(j * (j + 1)) / 2 + i];
+}
+
+template <typename T>
+void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
+{
+	uint32_t s = 0, d = 0;
+	for (uint32_t i = 0; i < m_n; ++i)
+	{
+		for (uint32_t j = 0; j < i; ++j)
+		{
+			if (i != di and j != dj and i != dj and j != di)
+			{
+				if (s != d)
+					m_data[d] = m_data[s];
+				++d;
+			}
+
+			++s;
+		}
+	}
+
+	--m_n;
+}
+
+template <typename T>
+class IdentityMatrix : public MatrixBase<T>
+{
+  public:
+	typedef typename MatrixBase<T>::value_type value_type;
+
+	IdentityMatrix(uint32_t n)
+		: m_n(n)
+	{
+	}
+
+	virtual uint32_t dim_m() const { return m_n; }
+	virtual uint32_t dim_n() const { return m_n; }
+
+	virtual value_type operator()(uint32_t i, uint32_t j) const
+	{
+		value_type result = 0;
+		if (i == j)
+			result = 1;
+		return result;
+	}
+
+  private:
+	uint32_t m_n;
+};
+
+// --------------------------------------------------------------------
+// matrix functions
+
+template <typename T>
+Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
+{
+	Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
+
+	for (uint32_t i = 0; i < result.dim_m(); ++i)
+	{
+		for (uint32_t j = 0; j < result.dim_n(); ++j)
+		{
+			for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
+				result(i, j) += lhs(li, j) * rhs(i, rj);
+		}
+	}
+
+	return result;
+}
+
+template <typename T>
+Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
+{
+	Matrix<T> result(lhs);
+	result *= rhs;
+
+	return result;
+}
+
+template <typename T>
+Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
+{
+	Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
+
+	for (uint32_t i = 0; i < result.dim_m(); ++i)
+	{
+		for (uint32_t j = 0; j < result.dim_n(); ++j)
+		{
+			result(i, j) = lhs(i, j) - rhs(i, j);
+		}
+	}
+
+	return result;
+}
+
+template <typename T>
+Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
+{
+	Matrix<T> result(lhs.dim_m(), lhs.dim_n());
+	result -= rhs;
+	return result;
+}
+
+// template <typename T>
+// symmetric_matrix<T> hammingDistance(const MatrixBase<T> &lhs, T rhs);
+
+// template <typename T>
+// std::vector<T> sum(const MatrixBase<T> &m);

include/cif++/Point.hpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -37,7 +37,7 @@
 namespace mmcif
 {
 
-typedef boost::math::quaternion<float> Quaternion;
+typedef boost::math::quaternion<float>	Quaternion;
 
 const double
 	kPI = 3.141592653589793238462643383279502884;
@@ -51,43 +51,26 @@ const double
 //	float x, y, z;
 //	tie(x, y, z) = atom.loc();
 
-template <typename F>
+template<typename F>
 struct PointF
 {
 	typedef F FType;
 
 	FType mX, mY, mZ;
+	
+	PointF()							: mX(0), mY(0), mZ(0) {}
+	PointF(FType x, FType y, FType z)	: mX(x), mY(y), mZ(z) {}
 
-	PointF()
-		: mX(0)
-		, mY(0)
-		, mZ(0)
-	{
-	}
-	PointF(FType x, FType y, FType z)
-		: mX(x)
-		, mY(y)
-		, mZ(z)
-	{
-	}
-
-	template <typename PF>
-	PointF(const PointF<PF> &pt)
+	template<typename PF>
+	PointF(const PointF<PF>& pt)
 		: mX(static_cast<F>(pt.mX))
 		, mY(static_cast<F>(pt.mY))
-		, mZ(static_cast<F>(pt.mZ))
-	{
-	}
+		, mZ(static_cast<F>(pt.mZ)) {}
 
-#if HAVE_LIBCLIPPER
-	PointF(const clipper::Coord_orth &pt)
-		: mX(pt[0])
-		, mY(pt[1])
-		, mZ(pt[2])
-	{
-	}
+#if HAVE_LIBCLIPPER	
+	PointF(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {}
 
-	PointF &operator=(const clipper::Coord_orth &rhs)
+	PointF& operator=(const clipper::Coord_orth& rhs)
 	{
 		mX = rhs[0];
 		mY = rhs[1];
@@ -96,72 +79,72 @@ struct PointF
 	}
 #endif
 
-	template <typename PF>
-	PointF &operator=(const PointF<PF> &rhs)
+	template<typename PF>
+	PointF& operator=(const PointF<PF>& rhs)
 	{
 		mX = static_cast<F>(rhs.mX);
 		mY = static_cast<F>(rhs.mY);
 		mZ = static_cast<F>(rhs.mZ);
 		return *this;
 	}
+	
+	FType& getX()			{ return mX; }
+	FType getX() const		{ return mX; }
+	void setX(FType x)		{ mX = x; }
 
-	FType &getX() { return mX; }
-	FType getX() const { return mX; }
-	void setX(FType x) { mX = x; }
+	FType& getY()			{ return mY; }
+	FType getY() const		{ return mY; }
+	void setY(FType y)		{ mY = y; }
 
-	FType &getY() { return mY; }
-	FType getY() const { return mY; }
-	void setY(FType y) { mY = y; }
-
-	FType &getZ() { return mZ; }
-	FType getZ() const { return mZ; }
-	void setZ(FType z) { mZ = z; }
-
-	PointF &operator+=(const PointF &rhs)
+	FType& getZ()			{ return mZ; }
+	FType getZ() const		{ return mZ; }
+	void setZ(FType z)		{ mZ = z; }
+	
+	PointF& operator+=(const PointF& rhs)
 	{
 		mX += rhs.mX;
 		mY += rhs.mY;
 		mZ += rhs.mZ;
-
+		
 		return *this;
 	}
-
-	PointF &operator+=(FType d)
+	
+	PointF& operator+=(FType d)
 	{
 		mX += d;
 		mY += d;
 		mZ += d;
-
+		
 		return *this;
 	}
 
-	PointF &operator-=(const PointF &rhs)
+	PointF& operator-=(const PointF& rhs)
 	{
 		mX -= rhs.mX;
 		mY -= rhs.mY;
 		mZ -= rhs.mZ;
-
+		
 		return *this;
 	}
 
-	PointF &operator-=(FType d)
+	PointF& operator-=(FType d)
 	{
 		mX -= d;
 		mY -= d;
 		mZ -= d;
-
+		
 		return *this;
 	}
 
-	PointF &operator*=(FType rhs)
+	PointF& operator*=(FType rhs)
 	{
 		mX *= rhs;
 		mY *= rhs;
 		mZ *= rhs;
 		return *this;
 	}
-
-	PointF &operator/=(FType rhs)
+	
+	PointF& operator/=(FType rhs)
 	{
 		mX /= rhs;
 		mY /= rhs;
@@ -179,18 +162,18 @@ struct PointF
 		}
 		return length;
 	}
-
-	void rotate(const boost::math::quaternion<FType> &q)
+	
+	void rotate(const boost::math::quaternion<FType>& q)
 	{
 		boost::math::quaternion<FType> p(0, mX, mY, mZ);
-
+		
 		p = q * p * boost::math::conj(q);
-
+	
 		mX = p.R_component_2();
 		mY = p.R_component_3();
 		mZ = p.R_component_4();
 	}
-
+	
 #if HAVE_LIBCLIPPER
 	operator clipper::Coord_orth() const
 	{
@@ -198,21 +181,21 @@ struct PointF
 	}
 #endif
 
-	operator std::tuple<const FType &, const FType &, const FType &>() const
+	operator std::tuple<const FType&, const FType&, const FType&>() const
 	{
 		return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
 	}
 
-	operator std::tuple<FType &, FType &, FType &>()
+	operator std::tuple<FType&,FType&,FType&>()
 	{
 		return std::make_tuple(std::ref(mX), std::ref(mY), std::ref(mZ));
 	}
-
-	bool operator==(const PointF &rhs) const
+	
+	bool operator==(const PointF& rhs) const
 	{
 		return mX == rhs.mX and mY == rhs.mY and mZ == rhs.mZ;
 	}
-
+	
 	// consider point as a vector... perhaps I should rename Point?
 	FType lengthsq() const
 	{
@@ -221,52 +204,52 @@ struct PointF
 
 	FType length() const
 	{
-		return std::sqrt(mX * mX + mY * mY + mZ * mZ);
+		return sqrt(mX * mX + mY * mY + mZ * mZ);
 	}
 };
 
 typedef PointF<float> Point;
 typedef PointF<double> DPoint;
 
-template <typename F>
-inline std::ostream &operator<<(std::ostream &os, const PointF<F> &pt)
+template<typename F>
+inline std::ostream& operator<<(std::ostream& os, const PointF<F>& pt)
 {
 	os << '(' << pt.mX << ',' << pt.mY << ',' << pt.mZ << ')';
-	return os;
+	return os; 
 }
 
-template <typename F>
-inline PointF<F> operator+(const PointF<F> &lhs, const PointF<F> &rhs)
+template<typename F>
+inline PointF<F> operator+(const PointF<F>& lhs, const PointF<F>& rhs)
 {
 	return PointF<F>(lhs.mX + rhs.mX, lhs.mY + rhs.mY, lhs.mZ + rhs.mZ);
 }
 
-template <typename F>
-inline PointF<F> operator-(const PointF<F> &lhs, const PointF<F> &rhs)
+template<typename F>
+inline PointF<F> operator-(const PointF<F>& lhs, const PointF<F>& rhs)
 {
 	return PointF<F>(lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ);
 }
 
-template <typename F>
-inline PointF<F> operator-(const PointF<F> &pt)
+template<typename F>
+inline PointF<F> operator-(const PointF<F>& pt)
 {
 	return PointF<F>(-pt.mX, -pt.mY, -pt.mZ);
 }
 
-template <typename F>
-inline PointF<F> operator*(const PointF<F> &pt, F f)
+template<typename F>
+inline PointF<F> operator*(const PointF<F>& pt, F f)
 {
 	return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
 }
 
-template <typename F>
-inline PointF<F> operator*(F f, const PointF<F> &pt)
+template<typename F>
+inline PointF<F> operator*(F f, const PointF<F>& pt)
 {
 	return PointF<F>(pt.mX * f, pt.mY * f, pt.mZ * f);
 }
 
-template <typename F>
-inline PointF<F> operator/(const PointF<F> &pt, F f)
+template<typename F>
+inline PointF<F> operator/(const PointF<F>& pt, F f)
 {
 	return PointF<F>(pt.mX / f, pt.mY / f, pt.mZ / f);
 }
@@ -274,96 +257,97 @@ inline PointF<F> operator/(const PointF<F> &pt, F f)
 // --------------------------------------------------------------------
 // several standard 3d operations
 
-template <typename F>
-inline double DistanceSquared(const PointF<F> &a, const PointF<F> &b)
+template<typename F>
+inline double DistanceSquared(const PointF<F>& a, const PointF<F>& b)
 {
-	return (a.mX - b.mX) * (a.mX - b.mX) +
-	       (a.mY - b.mY) * (a.mY - b.mY) +
-	       (a.mZ - b.mZ) * (a.mZ - b.mZ);
+	return
+		(a.mX - b.mX) * (a.mX - b.mX) +
+		(a.mY - b.mY) * (a.mY - b.mY) +
+		(a.mZ - b.mZ) * (a.mZ - b.mZ);
 }
 
-template <typename F>
-inline double Distance(const PointF<F> &a, const PointF<F> &b)
+template<typename F>
+inline double Distance(const PointF<F>& a, const PointF<F>& b)
 {
-	return std::sqrt(
+	return sqrt(
 		(a.mX - b.mX) * (a.mX - b.mX) +
 		(a.mY - b.mY) * (a.mY - b.mY) +
 		(a.mZ - b.mZ) * (a.mZ - b.mZ));
 }
 
-template <typename F>
-inline F DotProduct(const PointF<F> &a, const PointF<F> &b)
+template<typename F>
+inline F DotProduct(const PointF<F>& a, const PointF<F>& b)
 {
 	return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ;
 }
 
-template <typename F>
-inline PointF<F> CrossProduct(const PointF<F> &a, const PointF<F> &b)
+template<typename F>
+inline PointF<F> CrossProduct(const PointF<F>& a, const PointF<F>& b)
 {
 	return PointF<F>(a.mY * b.mZ - b.mY * a.mZ,
-		a.mZ * b.mX - b.mZ * a.mX,
-		a.mX * b.mY - b.mX * a.mY);
+				  a.mZ * b.mX - b.mZ * a.mX,
+				  a.mX * b.mY - b.mX * a.mY);
 }
 
-template <typename F>
-double Angle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3)
+template<typename F>
+double Angle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3)
 {
 	PointF<F> v1 = p1 - p2;
 	PointF<F> v2 = p3 - p2;
-
+	
 	return std::acos(DotProduct(v1, v2) / (v1.length() * v2.length())) * 180 / kPI;
 }
 
-template <typename F>
-double DihedralAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
+template<typename F>
+double DihedralAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4)
 {
-	PointF<F> v12 = p1 - p2; // vector from p2 to p1
-	PointF<F> v43 = p4 - p3; // vector from p3 to p4
-
-	PointF<F> z = p2 - p3; // vector from p3 to p2
-
+	PointF<F> v12 = p1 - p2;	// vector from p2 to p1
+	PointF<F> v43 = p4 - p3;	// vector from p3 to p4
+	
+	PointF<F> z = p2 - p3;		// vector from p3 to p2
+	
 	PointF<F> p = CrossProduct(z, v12);
 	PointF<F> x = CrossProduct(z, v43);
 	PointF<F> y = CrossProduct(z, x);
-
+	
 	double u = DotProduct(x, x);
 	double v = DotProduct(y, y);
-
+	
 	double result = 360;
 	if (u > 0 and v > 0)
 	{
-		u = DotProduct(p, x) / std::sqrt(u);
-		v = DotProduct(p, y) / std::sqrt(v);
+		u = DotProduct(p, x) / sqrt(u);
+		v = DotProduct(p, y) / sqrt(v);
 		if (u != 0 or v != 0)
-			result = std::atan2(v, u) * 180 / kPI;
+			result = atan2(v, u) * 180 / kPI;
 	}
-
+	
 	return result;
 }
 
-template <typename F>
-double CosinusAngle(const PointF<F> &p1, const PointF<F> &p2, const PointF<F> &p3, const PointF<F> &p4)
+template<typename F>
+double CosinusAngle(const PointF<F>& p1, const PointF<F>& p2, const PointF<F>& p3, const PointF<F>& p4)
 {
 	PointF<F> v12 = p1 - p2;
 	PointF<F> v34 = p3 - p4;
-
+	
 	double result = 0;
-
+	
 	double x = DotProduct(v12, v12) * DotProduct(v34, v34);
 	if (x > 0)
-		result = DotProduct(v12, v34) / std::sqrt(x);
-
+		result = DotProduct(v12, v34) / sqrt(x);
+	
 	return result;
 }
 
-template <typename F>
+template<typename F>
 auto DistancePointToLine(const PointF<F> &l1, const PointF<F> &l2, const PointF<F> &p)
 {
-	auto line = l2 - l1;
-	auto p_to_l1 = p - l1;
-	auto p_to_l2 = p - l2;
-	auto cross = CrossProduct(p_to_l1, p_to_l2);
-	return cross.length() / line.length();
+	auto line       = l2 - l1;
+    auto p_to_l1    = p - l1;
+    auto p_to_l2    = p - l2;
+    auto cross      = CrossProduct(p_to_l1, p_to_l2);
+    return cross.length() / line.length();
 }
 
 // --------------------------------------------------------------------
@@ -371,86 +355,74 @@ 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.
 
-Point Nudge(Point p, float offset);
+template<typename F>
+PointF<F> Nudge(PointF<F> p, F 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);
-
-/// \brief Returns how the two sets of points \a a and \b b can be aligned
-///        
-/// \param a	The first set of points
-/// \param b    The second set of points
-/// \result     The quaternion which should be applied to the points in \a a to
-///             obtain the best superposition.
-Quaternion AlignPoints(const std::vector<Point> &a, const std::vector<Point> &b);
-
-/// \brief The RMSd for the points in \a a and \a b
-double RMSd(const std::vector<Point> &a, const std::vector<Point> &b);
+std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q);
+Point Centroid(std::vector<Point>& Points);
+Point CenterPoints(std::vector<Point>& Points);
+Quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
+double RMSd(const std::vector<Point>& a, const std::vector<Point>& b);
 
 // --------------------------------------------------------------------
 // Helper class to generate evenly divided Points on a sphere
 // we use a fibonacci sphere to calculate even distribution of the dots
 
-template <int N>
+template<int N>
 class SphericalDots
 {
   public:
-	enum
-	{
-		P = 2 * N + 1
-	};
-	typedef typename std::array<Point, P> array_type;
-	typedef typename array_type::const_iterator iterator;
+	enum { P = 2 * N + 1 };
+	typedef typename std::array<Point,P>	array_type;
+	typedef typename array_type::const_iterator	iterator;
 
-	static SphericalDots &instance()
+	static SphericalDots& instance()
 	{
 		static SphericalDots sInstance;
 		return sInstance;
 	}
+	
+	size_t size() const							{ return mPoints.size(); }
+	const Point operator[](uint32_t inIx) const	{ return mPoints[inIx]; }
+	iterator begin() const						{ return mPoints.begin(); }
+	iterator end() const						{ return mPoints.end(); }
 
-	size_t size() const { return mPoints.size(); }
-	const Point operator[](uint32_t inIx) const { return mPoints[inIx]; }
-	iterator begin() const { return mPoints.begin(); }
-	iterator end() const { return mPoints.end(); }
-
-	double weight() const { return mWeight; }
+	double weight() const						{ return mWeight; }
 
 	SphericalDots()
 	{
-
+				
 		const double
 			kGoldenRatio = (1 + std::sqrt(5.0)) / 2;
-
+		
 		mWeight = (4 * kPI) / P;
-
+		
 		auto p = mPoints.begin();
-
+		
 		for (int32_t i = -N; i <= N; ++i)
 		{
 			double lat = std::asin((2.0 * i) / P);
 			double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio;
-
-			p->mX = std::sin(lon) * std::cos(lat);
-			p->mY = std::cos(lon) * std::cos(lat);
-			p->mZ = std::sin(lat);
+			
+			p->mX = sin(lon) * cos(lat);
+			p->mY = cos(lon) * cos(lat);
+			p->mZ =            sin(lat);
 
 			++p;
 		}
 	}
 
   private:
-	array_type mPoints;
-	double mWeight;
+
+	array_type				mPoints;
+	double					mWeight;
 };
 
 typedef SphericalDots<50> SphericalDots_50;
 
-} // namespace mmcif
+}

include/cif++/Secondary.hpp

@@ -137,22 +137,6 @@ class DSSP
 	  public:
 		friend class iterator;
 
-		ResidueInfo()
-			: mImpl(nullptr)
-		{
-		}
-
-		ResidueInfo(const ResidueInfo &rhs)
-			: mImpl(rhs.mImpl)
-		{
-		}
-
-		ResidueInfo& operator=(const ResidueInfo &rhs)
-		{
-			mImpl = rhs.mImpl;
-			return *this;
-		}
-
 		explicit operator bool() const		{ return not empty(); }
 		bool empty() const					{ return mImpl == nullptr; }
 
@@ -184,15 +168,6 @@ 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) {}
 
@@ -202,7 +177,7 @@ class DSSP
 	class iterator
 	{
 	  public:
-		using iterator_category = std::bidirectional_iterator_tag;
+		using iterator_category = std::input_iterator_tag;
 		using value_type = ResidueInfo;
 		using difference_type = std::ptrdiff_t;
 		using pointer = value_type*;
@@ -223,14 +198,6 @@ 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,112 +60,30 @@ 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 *compound() 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(std::shared_ptr<AtomImpl> impl)
-		: mImpl(impl)
-	{
-	}
-
-	Atom(const Atom &rhs)
-		: mImpl(rhs.mImpl)
-	{
-	}
+	Atom();
+	Atom(struct AtomImpl *impl);
+	Atom(const Atom &rhs);
 
 	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);
 
-	explicit operator bool() const { return (bool)mImpl; }
+	~Atom();
+
+	explicit operator bool() const { return mImpl_ != nullptr; }
 
 	// return a copy of this atom, with data copied instead of referenced
-	Atom clone() const
-	{
-		auto copy = std::make_shared<AtomImpl>(*mImpl);
-		copy->mClone = true;
-		return Atom(copy);
-	}
+	Atom clone() const;
 
-	Atom &operator=(const Atom &rhs) = default;
+	Atom &operator=(const Atom &rhs);
 
-	template <typename T>
-	T get_property(const std::string_view name) const;
+	const std::string &id() const;
+	AtomType type() const;
 
-	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);
-	}
+	Point location() const;
+	void location(Point p);
 
 	/// \brief Translate the position of this atom by \a t
 	void translate(Point t);
@@ -173,40 +91,47 @@ 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 { return impl().mRow; }
+	const cif::Row getRow() const;
 	const cif::Row getRowAniso() const;
 
-	bool isSymmetryCopy() const { return impl().mSymmetryCopy; }
-	std::string symmetry() const { return impl().mSymmetryOperator; }
+	// Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
+	bool isSymmetryCopy() const;
+	std::string symmetry() const;
+	// const clipper::RTop_orth& symop() const;
 
-	const Compound &compound() const;
-	bool isWater() const { return impl().mCompID == "HOH" or impl().mCompID == "H2O" or impl().mCompID == "WAT"; }
+	const Compound &comp() const;
+	bool isWater() const;
 	int charge() const;
 
 	float uIso() const;
-	bool getAnisoU(float anisou[6]) const { return impl().getAnisoU(anisou); }
+	bool getAnisoU(float anisou[6]) const;
 	float occupancy() const;
 
+	template <typename T>
+	T property(const std::string &name) const;
+
+	void property(const std::string &name, const std::string &value);
+
+	template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
+	void property(const std::string &name, const T &value)
+	{
+		property(name, std::to_string(value));
+	}
+
 	// specifications
-	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 labelAtomID() const;
+	std::string labelCompID() const;
+	std::string labelAsymID() const;
 	std::string labelEntityID() const;
-	int labelSeqID() const { return impl().mSeqID; }
-	const std::string &labelAltID() const { return impl().mAltID; }
-	bool isAlternate() const { return not impl().mAltID.empty(); }
+	int labelSeqID() const;
+	std::string labelAltID() const;
+	bool isAlternate() const;
 
 	std::string authAtomID() const;
 	std::string authCompID() const;
 	std::string authAsymID() const;
-	const std::string &authSeqID() const { return impl().mAuthSeqID; }
+	std::string authSeqID() const;
 	std::string pdbxAuthInsCode() const;
 	std::string pdbxAuthAltID() const;
 
@@ -214,10 +139,13 @@ 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;
-	bool operator!=(const Atom &rhs) const
-	{
-		return not operator==(rhs);
-	}
+
+	// // get clipper format Atom
+	// clipper::Atom toClipper() const;
+
+	// Radius calculation based on integrating the density until perc of electrons is found
+	void calculateRadius(float resHigh, float resLow, float perc);
+	float radius() const;
 
 	// access data in compound for this atom
 
@@ -230,10 +158,10 @@ class Atom
 
 	void swap(Atom &b)
 	{
-		std::swap(mImpl, b.mImpl);
+		std::swap(mImpl_, b.mImpl_);
 	}
 
-	int compare(const Atom &b) const { return impl().compare(*b.mImpl); }
+	int compare(const Atom &b) const;
 
 	bool operator<(const Atom &rhs) const
 	{
@@ -242,45 +170,16 @@ 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);
 
-	const AtomImpl &impl() const
-	{
-		if (not mImpl)
-			throw std::runtime_error("Uninitialized atom, not found?");
-		return *mImpl;
-	}
+	AtomImpl *impl();
+	const AtomImpl *impl() const;
 
-	std::shared_ptr<AtomImpl> mImpl;
+	struct 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);
@@ -300,26 +199,22 @@ typedef std::vector<Atom> AtomView;
 
 // --------------------------------------------------------------------
 
-enum class EntityType
-{
-	Polymer, NonPolymer, Macrolide, Water, Branched
-};
-
-// --------------------------------------------------------------------
-
 class Residue
 {
   public:
-	// constructor
+	// constructors should be private, but that's not possible for now (needed in emplace)
+
+	// constructor for waters
 	Residue(const Structure &structure, const std::string &compoundID,
-		const std::string &asymID, int seqID, const std::string &authSeqID)
-		: mStructure(&structure)
-		, mCompoundID(compoundID)
-		, mAsymID(asymID)
-		, mSeqID(seqID)
-		, mAuthSeqID(authSeqID)
-	{
-	}
+		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);
 
 	Residue(const Residue &rhs) = delete;
 	Residue &operator=(const Residue &rhs) = delete;
@@ -330,12 +225,8 @@ 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;
 
@@ -345,14 +236,10 @@ class Residue
 	Atom atomByID(const std::string &atomID) const;
 
 	const std::string &compoundID() const { return mCompoundID; }
-	void setCompoundID(const std::string &id) { mCompoundID = id; }
-
 	const std::string &asymID() const { return mAsymID; }
 	int seqID() const { return mSeqID; }
 	std::string entityID() const;
 
-	EntityType entityType() const;
-
 	std::string authAsymID() const;
 	std::string authSeqID() const;
 	std::string authInsCode() const;
@@ -388,18 +275,6 @@ 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() {}
 
@@ -408,6 +283,9 @@ 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;
 };
@@ -474,11 +352,6 @@ 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;
@@ -516,95 +389,35 @@ 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);
-
-	Sugar(Sugar &&rhs);
-	Sugar &operator=(Sugar &&rhs);
-
-	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; }
-
-	size_t getLinkNr() const
-	{
-		return mLink ? std::stoi(mLink.authSeqID()) : 0;
-	}
-
-  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 cif::File
+class File : public std::enable_shared_from_this<File>
 {
   public:
-	File() {}
-
-	File(const std::filesystem::path &path)
-	{
-		load(path);
-	}
-
-	File(const char *data, size_t length)
-	{
-		load(data, length);
-	}
+	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(const File &) = delete;
 	File &operator=(const File &) = delete;
 
-	void load(const std::filesystem::path &p) override;
-	void save(const std::filesystem::path &p) override;
+	cif::Datablock& createDatablock(const std::string &name);
 
-	void load(std::istream &is) override;
+	void load(const std::filesystem::path &path);
+	void save(const std::filesystem::path &path);
 
-	using cif::File::load;
-	using cif::File::save;
+	Structure *model(size_t nr = 1);
 
-	cif::Datablock &data() { return front(); }
+	struct FileImpl &impl() const { return *mImpl; }
+
+	cif::Datablock &data();
+	cif::File &file();
+
+  private:
+	struct FileImpl *mImpl;
 };
 
 // --------------------------------------------------------------------
@@ -624,123 +437,60 @@ inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
 class Structure
 {
   public:
-	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 = {});
-
-	Structure(Structure &&s) = default;
+	Structure(File &p, size_t modelNr = 1, StructureOpenOptions options = {});
+	Structure &operator=(const Structure &) = delete;
+	~Structure();
 
 	// Create a read-only clone of the current structure (for multithreaded calculations that move atoms)
 	Structure(const Structure &);
 
-	Structure &operator=(const Structure &) = delete;
-	// Structure &operator=(Structure &&s) = default;
-
-	~Structure();
+	File &getFile() const;
 
 	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; }
 
-	Polymer &getPolymerByAsymID(const std::string &asymID);
-
-	const Polymer &getPolymerByAsymID(const std::string &asymID) const
-	{
-		return const_cast<Structure *>(this)->getPolymerByAsymID(asymID);
-	}
-
-	const std::list<Branch> &branches() const { return 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;
+	const std::vector<Residue> &branchResidues() const { return mBranchResidues; }
+
+	Atom getAtomByID(std::string id) const;
 	// Atom getAtomByLocation(Point pt, float maxDistance) const;
 
 	Atom getAtomByLabel(const std::string &atomID, const std::string &asymID,
 		const std::string &compID, int seqID, const std::string &altID = "");
 
-	/// \brief Return the atom closest to point \a p
-	Atom getAtomByPosition(Point p) const;
+	/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
+	const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0) const;
 
-	/// \brief Return the atom closest to point \a p 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;
+	// map between auth and label locations
 
-	/// \brief Get a non-poly residue for an asym with id \a asymID
-	Residue &getResidue(const std::string &asymID)
-	{
-		return getResidue(asymID, 0, "");
-	}
+	std::tuple<std::string, int, std::string> MapAuthToLabel(const std::string &asymID,
+		const std::string &seqID, const std::string &compID, const std::string &insCode = "");
 
-	/// \brief Get a non-poly residue for an asym with id \a asymID
-	const Residue &getResidue(const std::string &asymID) const
-	{
-		return getResidue(asymID, 0, "");
-	}
+	std::tuple<std::string, std::string, std::string, std::string> MapLabelToAuth(
+		const std::string &asymID, int seqID, const std::string &compID);
 
-	/// \brief Get a residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
-	Residue &getResidue(const std::string &asymID, int seqID, const std::string &authSeqID);
+	// returns chain, seqnr, icode
+	std::tuple<char, int, char> MapLabelToAuth(
+		const std::string &asymID, int seqID) const;
 
-	/// \brief Get a the single residue for an asym with id \a asymID seq id \a seqID and authSeqID \a authSeqID
-	const Residue &getResidue(const std::string &asymID, int seqID, const std::string &authSeqID) const
-	{
-		return const_cast<Structure *>(this)->getResidue(asymID, seqID, authSeqID);
-	}
+	// returns chain,seqnr,comp,iCode
+	std::tuple<std::string, int, std::string, std::string> MapLabelToPDB(
+		const std::string &asymID, int seqID, const std::string &compID,
+		const std::string &authSeqID) const;
 
-	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
-	Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID);
-
-	/// \brief Get a residue for an asym with id \a asymID, compound id \a compID, seq id \a seqID and authSeqID \a authSeqID
-	const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID, const std::string &authSeqID) const
-	{
-		return const_cast<Structure *>(this)->getResidue(asymID, compID, seqID, authSeqID);
-	}
-
-	/// \brief Get a the residue for atom \a atom
-	Residue &getResidue(const mmcif::Atom &atom)
-	{
-		return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID(), atom.authSeqID());
-	}
-
-	/// \brief Get a the residue for atom \a atom
-	const Residue &getResidue(const mmcif::Atom &atom) const
-	{
-		return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID(), atom.authSeqID());
-	}
+	std::tuple<std::string, int, std::string> MapPDBToLabel(
+		const std::string &asymID, int seqID, const std::string &compID, const std::string &iCode) const;
 
 	// Actions
-	void removeAtom(Atom &a)
-	{
-		removeAtom(a, true);
-	}
-
-	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,
+	void removeAtom(Atom &a);
+	void swapAtoms(Atom &a1, Atom &a2); // swap the labels for these atoms
+	void moveAtom(Atom &a, Point p);    // move atom to a new location
+	void changeResidue(const Residue &res, const std::string &newCompound,
 		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
@@ -754,33 +504,9 @@ class Structure
 	/// \return				The newly create asym ID
 	std::string createNonpoly(const std::string &entity_id, const std::vector<mmcif::Atom> &atoms);
 
-	/// \brief Create a new NonPolymer struct_asym with atoms constructed from info in \a atom_info, returns asym_id.
-	/// This method creates new atom records filled with info from the info.
-	///
-	/// \param entity_id	The entity ID of the new nonpoly
-	/// \param atoms		The array of sets of 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 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 Translate the coordinates of all atoms in the structure by \a t
 	void translate(Point t);
@@ -788,59 +514,33 @@ class Structure
 	/// \brief Rotate the coordinates of all atoms in the structure by \a q
 	void rotate(Quaternion t);
 
-	/// \brief Translate and rotate the coordinates of all atoms in the structure by \a t and \a q
-	void translateAndRotate(Point t, Quaternion q);
-
-	/// \brief Translate, rotate and translate again the coordinates of all atoms in the structure by \a t1 , \a q and \a t2
-	void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
-
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
+	const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }
 
 	void cleanupEmptyCategories();
 
-	/// \brief Direct access to underlying data
-	cif::Category &category(std::string_view name) const
-	{
-		return mDb[name];
-	}
-
-	cif::Datablock &datablock() const
-	{
-		return mDb;
-	}
-
-	void validateAtoms() const;
-
   private:
 	friend Polymer;
 	friend Residue;
+	// friend residue_view;
+	// friend residue_iterator;
+
+	cif::Category &category(const char *name) const;
+	cif::Datablock &datablock() const;
 
 	std::string insertCompound(const std::string &compoundID, bool isEntity);
 
-	std::string createEntityForBranch(Branch &branch);
-
 	void loadData();
+	void updateAtomIndex();
 
 	void loadAtomsForModel(StructureOpenOptions options);
 
-	template<typename... Args>
-	Atom& emplace_atom(Args ...args)
-	{
-		return emplace_atom(Atom{std::forward<Args>(args)...});
-	}
-
-	Atom &emplace_atom(Atom &&atom);
-
-	void removeAtom(Atom &a, bool removeFromResidue);
-	void removeSugar(Sugar &sugar);
-
-	cif::Datablock &mDb;
+	File &mFile;
 	size_t mModelNr;
 	AtomView mAtoms;
 	std::vector<size_t> mAtomIndex;
 	std::list<Polymer> mPolymers;
-	std::list<Branch> mBranches;
-	std::vector<Residue> mNonPolymers;
+	std::vector<Residue> mNonPolymers, mBranchResidues;
 };
 
 } // namespace mmcif

include/cif++/Symmetry.hpp

@@ -37,11 +37,6 @@ namespace mmcif
 
 // --------------------------------------------------------------------
 
-enum class SpacegroupName
-{
-	full, xHM, Hall
-};
-
 struct Spacegroup
 {
 	const char* name;
@@ -138,7 +133,6 @@ CIFPP_EXPORT extern const std::size_t kSymopNrTableSize;
 
 // --------------------------------------------------------------------
 
-int GetSpacegroupNumber(std::string spacegroup);						// alternative for clipper's parsing code, using SpacegroupName::full
-int GetSpacegroupNumber(std::string spacegroup, SpacegroupName type);	// alternative for clipper's parsing code
+int GetSpacegroupNumber(std::string spacegroup);	// alternative for clipper's parsing code
 
 }

src/BondMap.cpp

@@ -1,17 +1,17 @@
 /*-
  * SPDX-License-Identifier: BSD-2-Clause
- *
+ * 
  * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
+ * 
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions are met:
- *
+ * 
  * 1. Redistributions of source code must retain the above copyright notice, this
  *    list of conditions and the following disclaimer
  * 2. Redistributions in binary form must reproduce the above copyright notice,
  *    this list of conditions and the following disclaimer in the documentation
  *    and/or other materials provided with the distribution.
- *
+ * 
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -24,14 +24,14 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include <algorithm>
 #include <fstream>
+#include <algorithm>
 #include <mutex>
 
-#include "cif++/BondMap.hpp"
 #include "cif++/Cif++.hpp"
-#include "cif++/CifUtils.hpp"
 #include "cif++/Compound.hpp"
+#include "cif++/CifUtils.hpp"
+#include "cif++/BondMap.hpp"
 
 namespace mmcif
 {
@@ -39,79 +39,223 @@ namespace mmcif
 namespace
 {
 
-	union IDType
+union IDType
+{
+	IDType() : id_n(0){}
+	IDType(const IDType& rhs) : id_n(rhs.id_n) {}
+	IDType(const std::string& s)
+		: IDType()
 	{
-		IDType()
-			: id_n(0)
-		{
-		}
-		IDType(const IDType &rhs)
-			: id_n(rhs.id_n)
-		{
-		}
-		IDType(const std::string &s)
-			: IDType()
-		{
-			assert(s.length() <= 4);
-			if (s.length() > 4)
-				throw BondMapException("Atom ID '" + s + "' is too long");
-			std::copy(s.begin(), s.end(), id_s);
-		}
+		assert(s.length() <= 4);
+		if (s.length() > 4)
+			throw BondMapException("Atom ID '" + s + "' is too long");
+		std::copy(s.begin(), s.end(), id_s);
+	}
 
-		IDType &operator=(const IDType &rhs)
-		{
-			id_n = rhs.id_n;
-			return *this;
-		}
+	IDType& operator=(const IDType& rhs)
+	{
+		id_n = rhs.id_n;
+		return *this;
+	}
 
-		IDType &operator=(const std::string &s)
-		{
-			id_n = 0;
-			assert(s.length() <= 4);
-			if (s.length() > 4)
-				throw BondMapException("Atom ID '" + s + "' is too long");
-			std::copy(s.begin(), s.end(), id_s);
-			return *this;
-		}
+	IDType& operator=(const std::string& s)
+	{
+		id_n = 0;
+		assert(s.length() <= 4);
+		if (s.length() > 4)
+			throw BondMapException("Atom ID '" + s + "' is too long");
+		std::copy(s.begin(), s.end(), id_s);
+		return *this;
+	}
 
-		bool operator<(const IDType &rhs) const
-		{
-			return id_n < rhs.id_n;
-		}
+	bool operator<(const IDType& rhs) const
+	{
+		return id_n < rhs.id_n;
+	}
 
-		bool operator<=(const IDType &rhs) const
-		{
-			return id_n <= rhs.id_n;
-		}
+	bool operator<=(const IDType& rhs) const
+	{
+		return id_n <= rhs.id_n;
+	}
 
-		bool operator==(const IDType &rhs) const
-		{
-			return id_n == rhs.id_n;
-		}
+	bool operator==(const IDType& rhs) const
+	{
+		return id_n == rhs.id_n;
+	}
 
-		bool operator!=(const IDType &rhs) const
-		{
-			return id_n != rhs.id_n;
-		}
+	bool operator!=(const IDType& rhs) const
+	{
+		return id_n != rhs.id_n;
+	}
 
-		char id_s[4];
-		uint32_t id_n;
-	};
+	char		id_s[4];
+	uint32_t	id_n;
+};
 
-	static_assert(sizeof(IDType) == 4, "atom_id_type should be 4 bytes");
-} // namespace
+static_assert(sizeof(IDType) == 4, "atom_id_type should be 4 bytes");	
+}
+
+// // --------------------------------------------------------------------
+
+// void createBondInfoFile(const fs::path& components, const fs::path& infofile)
+// {
+// 	std::ofstream outfile(infofile.string() + ".tmp", std::ios::binary);
+// 	if (not outfile.is_open())
+// 		throw BondMapException("Could not create bond info file " + infofile.string() + ".tmp");
+
+// 	cif::File infile(components);
+
+// 	std::set<atom_id_type> atomIDs;
+// 	std::vector<atom_id_type> compoundIDs;
+
+// 	for (auto& db: infile)
+// 	{
+// 		auto chem_comp_bond = db.get("chem_comp_bond");
+// 		if (not chem_comp_bond)
+// 		{
+// 			if (cif::VERBOSE > 1)
+// 				std::cerr << "Missing chem_comp_bond category in data block " << db.getName() << std::endl;
+// 			continue;
+// 		}
+
+// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
+// 		{
+// 			atomIDs.insert(atom_id_1);
+// 			atomIDs.insert(atom_id_2);
+// 		}
+
+// 		compoundIDs.push_back({ db.getName() });
+// 	}
+
+// 	if (cif::VERBOSE)
+// 		std::cout << "Number of unique atom names is " << atomIDs.size() << std::endl
+// 				  << "Number of unique residue names is " << compoundIDs.size() << std::endl;
+
+// 	CompoundBondInfoFileHeader header = {};
+// 	header.indexEntries = compoundIDs.size();
+// 	header.atomEntries = atomIDs.size();
+
+// 	outfile << header;
+	
+// 	for (auto atomID: atomIDs)
+// 		outfile << atomID;
+
+// 	auto dataOffset = outfile.tellp();
+
+// 	std::vector<CompoundBondInfo> entries;
+// 	entries.reserve(compoundIDs.size());
+
+// 	std::map<atom_id_type, uint16_t> atomIDMap;
+// 	for (auto& atomID: atomIDs)
+// 		atomIDMap[atomID] = atomIDMap.size();
+
+// 	for (auto& db: infile)
+// 	{
+// 		auto chem_comp_bond = db.get("chem_comp_bond");
+// 		if (not chem_comp_bond)
+// 			continue;
+
+// 		std::set<uint16_t> bondedAtoms;
+
+// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
+// 		{
+// 			bondedAtoms.insert(atomIDMap[atom_id_1]);
+// 			bondedAtoms.insert(atomIDMap[atom_id_2]);
+// 		}
+
+// 		std::map<uint16_t, int32_t> bondedAtomMap;
+// 		for (auto id: bondedAtoms)
+// 			bondedAtomMap[id] = static_cast<int32_t>(bondedAtomMap.size());
+		
+// 		CompoundBondInfo info = {
+// 			db.getName(),
+// 			static_cast<uint32_t>(bondedAtomMap.size()),
+// 			outfile.tellp() - dataOffset
+// 		};
+
+// 		entries.push_back(info);
+
+// 		// An now first write the array of atom ID's in this compound
+// 		for (uint16_t id: bondedAtoms)
+// 			write(outfile, id);
+
+// 		// And then the symmetric matrix with bonds
+// 		size_t N = bondedAtoms.size();
+// 		size_t M = (N * (N - 1)) / 2;
+
+// 		size_t K = M / 8;
+// 		if (M % 8)
+// 			K += 1;
+		
+// 		std::vector<uint8_t> m(K);
+
+// 		for (const auto& [atom_id_1, atom_id_2]: chem_comp_bond->rows<std::string,std::string>({"atom_id_1", "atom_id_2"}))
+// 		{
+// 			auto a = bondedAtomMap[atomIDMap[atom_id_1]];
+// 			auto b = bondedAtomMap[atomIDMap[atom_id_2]];
+
+// 			assert(a != b);
+// 			assert((int)b < (int)N);
+
+// 			if (a > b)
+// 				std::swap(a, b);
+			
+// 			size_t ix = ((b - 1) * b) / 2 + a;
+// 			assert(ix < M);
+
+// 			auto Bix = ix / 8;
+// 			auto bix = ix % 8;
+
+// 			m[Bix] |= 1 << bix;
+// 		}
+
+// 		outfile.write(reinterpret_cast<char*>(m.data()), m.size());
+// 	}
+
+// 	header.dataSize = outfile.tellp() - dataOffset;
+
+// 	std::sort(entries.begin(), entries.end(), [](CompoundBondInfo& a, CompoundBondInfo& b)
+// 	{
+// 		return a.id < b.id;
+// 	});
+
+// 	for (auto& info: entries)
+// 		outfile << info;
+
+// 	outfile.seekp(0);
+// 	outfile << header;
+
+// 	// compress
+// 	outfile.close();
+
+// 	std::ifstream in(infofile.string() + ".tmp", std::ios::binary);
+// 	std::ofstream out(infofile, std::ios::binary);
+
+// 	{
+// 		io::filtering_stream<io::output> os;
+// 		os.push(io::gzip_compressor());
+// 		os.push(out);
+// 		io::copy(in, os);
+// 	}
+
+// 	in.close();
+// 	out.close();
+
+// 	fs::remove(infofile.string() + ".tmp");
+// }
 
 // --------------------------------------------------------------------
 
 struct CompoundBondInfo
 {
-	IDType mID;
-	std::set<std::tuple<uint32_t, uint32_t>> mBonded;
+	IDType	mID;
+	std::set<std::tuple<uint32_t,uint32_t>> mBonded;
 
 	bool bonded(uint32_t a1, uint32_t a2) const
 	{
-		return mBonded.count({a1, a2}) > 0;
+		return mBonded.count({ a1, a2 }) > 0;
 	}
+
 };
 
 // --------------------------------------------------------------------
@@ -119,18 +263,20 @@ struct CompoundBondInfo
 class CompoundBondMap
 {
   public:
+
 	static CompoundBondMap &instance()
 	{
 		static std::unique_ptr<CompoundBondMap> s_instance(new CompoundBondMap);
 		return *s_instance;
 	}
 
-	bool bonded(const std::string &compoundID, const std::string &atomID1, const std::string &atomID2);
+	bool bonded(const std::string& compoundID, const std::string& atomID1, const std::string& atomID2);
 
   private:
+
 	CompoundBondMap() {}
 
-	uint32_t getAtomID(const std::string &atomID)
+	uint32_t getAtomID(const std::string& atomID)
 	{
 		IDType id(atomID);
 
@@ -144,16 +290,16 @@ class CompoundBondMap
 		}
 		else
 			result = i->second;
-
+		
 		return result;
 	}
 
-	std::map<IDType, uint32_t> mAtomIDIndex;
+	std::map<IDType,uint32_t> mAtomIDIndex;
 	std::vector<CompoundBondInfo> mCompounds;
 	std::mutex mMutex;
 };
 
-bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &atomID1, const std::string &atomID2)
+bool CompoundBondMap::bonded(const std::string &compoundID, const std::string& atomID1, const std::string& atomID2)
 {
 	std::lock_guard lock(mMutex);
 
@@ -164,35 +310,32 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &a
 	uint32_t a2 = getAtomID(atomID2);
 	if (a1 > a2)
 		std::swap(a1, a2);
-
-	for (auto &bi : mCompounds)
+	
+	for (auto &bi: mCompounds)
 	{
 		if (bi.mID != id)
 			continue;
-
+		
 		return bi.bonded(a1, a2);
 	}
 
 	bool result = false;
 
 	// not found in our cache, calculate
-	CompoundBondInfo bondInfo{id};
+	CompoundBondInfo bondInfo{ id };
 
 	auto compound = mmcif::CompoundFactory::instance().create(compoundID);
 	if (not compound)
-	{
-		if (cif::VERBOSE >= 0)
-			std::cerr << "Missing compound bond info for " << compoundID << std::endl;
-	}
+		std::cerr << "Missing compound bond info for " << compoundID << std::endl;
 	else
 	{
-		for (auto &atom : compound->bonds())
+		for (auto &atom: compound->bonds())
 		{
 			uint32_t ca1 = getAtomID(atom.atomID[0]);
 			uint32_t ca2 = getAtomID(atom.atomID[1]);
 			if (ca1 > ca2)
 				std::swap(ca1, ca2);
-
+			
 			bondInfo.mBonded.insert({ca1, ca2});
 			result = result or (a1 == ca1 and a2 == ca2);
 		}
@@ -205,27 +348,27 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &a
 
 // --------------------------------------------------------------------
 
-BondMap::BondMap(const Structure &p)
+BondMap::BondMap(const Structure& p)
 {
-	auto &compoundBondInfo = CompoundBondMap::instance();
+	auto& compoundBondInfo = CompoundBondMap::instance();
 
 	auto atoms = p.atoms();
 	dim = uint32_t(atoms.size());
 
-	//	bond = std::vector<bool>(dim * (dim - 1), false);
+//	bond = std::vector<bool>(dim * (dim - 1), false);
 
-	for (auto &atom : atoms)
+	for (auto& atom: atoms)
 		index[atom.id()] = uint32_t(index.size());
-
-	auto bindAtoms = [this](const std::string &a, const std::string &b)
+	
+	auto bindAtoms = [this](const std::string& a, const std::string& b)
 	{
 		uint32_t ixa = index[a];
 		uint32_t ixb = index[b];
-
+		
 		bond.insert(key(ixa, ixb));
 	};
 
-	auto linkAtoms = [this, &bindAtoms](const std::string &a, const std::string &b)
+	auto linkAtoms = [this,&bindAtoms](const std::string& a, const std::string& b)
 	{
 		bindAtoms(a, b);
 
@@ -233,20 +376,20 @@ BondMap::BondMap(const Structure &p)
 		link[b].insert(a);
 	};
 
-	cif::Datablock &db = p.datablock();
+	cif::Datablock& db = p.getFile().data();
 
 	// collect all compounds first
 	std::set<std::string> compounds;
-	for (auto c : db["chem_comp"])
+	for (auto c: db["chem_comp"])
 		compounds.insert(c["id"].as<std::string>());
-
+	
 	// make sure we also have all residues in the polyseq
-	for (auto m : db["entity_poly_seq"])
+	for (auto m: db["entity_poly_seq"])
 	{
 		std::string c = m["mon_id"].as<std::string>();
 		if (compounds.count(c))
 			continue;
-
+		
 		if (cif::VERBOSE > 1)
 			std::cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << std::endl;
 		compounds.insert(c);
@@ -255,84 +398,68 @@ BondMap::BondMap(const Structure &p)
 	cif::Progress progress(compounds.size(), "Creating bond map");
 
 	// some helper indices to speed things up a bit
-	std::map<std::tuple<std::string, int, std::string, std::string>, std::string> atomMapByAsymSeqAndAtom;
-	for (auto &a : p.atoms())
+	std::map<std::tuple<std::string,int,std::string>,std::string> atomMapByAsymSeqAndAtom;
+	for (auto& a: p.atoms())
 	{
-		auto key = make_tuple(a.labelAsymID(), a.labelSeqID(), a.labelAtomID(), a.authSeqID());
+		auto key = make_tuple(a.labelAsymID(), a.labelSeqID(), a.labelAtomID());
 		atomMapByAsymSeqAndAtom[key] = a.id();
 	}
-
+	
 	// first link all residues in a polyseq
-
-	std::string lastAsymID, lastAuthSeqID;
+	
+	std::string lastAsymID;
 	int lastSeqID = 0;
-	for (const auto &[asymID, seqID, authSeqID] : db["pdbx_poly_seq_scheme"].rows<std::string,int,std::string>("asym_id", "seq_id", "pdb_seq_num"))
+	for (auto r: db["pdbx_poly_seq_scheme"])
 	{
-		if (asymID != lastAsymID) // first in a new sequece
+		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"])
+	for (auto l: db["struct_conn"])
 	{
 		std::string asym1, asym2, atomId1, atomId2;
 		int seqId1 = 0, seqId2 = 0;
-		std::string authSeqId1, authSeqId2;
-
-		cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2, authSeqId1, authSeqId2) =
+		cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2) =
 			l.get("ptnr1_label_asym_id", "ptnr2_label_asym_id",
-				"ptnr1_label_atom_id", "ptnr2_label_atom_id",
-				"ptnr1_label_seq_id", "ptnr2_label_seq_id",
-				"ptnr1_auth_seq_id", "ptnr2_auth_seq_id");
-
-		auto ka = make_tuple(asym1, seqId1, atomId1, authSeqId1);
-		auto kb = make_tuple(asym2, seqId2, atomId2, authSeqId2);
-
-		if (atomMapByAsymSeqAndAtom.count(ka) and atomMapByAsymSeqAndAtom.count(kb))
-		{
-			auto a = atomMapByAsymSeqAndAtom.at(ka);
-			auto b = atomMapByAsymSeqAndAtom.at(kb);
+				  "ptnr1_label_atom_id", "ptnr2_label_atom_id",
+				  "ptnr1_label_seq_id", "ptnr2_label_seq_id");
 
+		std::string a = atomMapByAsymSeqAndAtom[make_tuple(asym1, seqId1, atomId1)];
+		std::string b = atomMapByAsymSeqAndAtom[make_tuple(asym2, seqId2, atomId2)];
+			
+		if (not (a.empty() or b.empty()))
 			linkAtoms(a, b);
-		}
-
-		// std::string a = atomMapByAsymSeqAndAtom.at(make_tuple(asym1, seqId1, atomId1, authSeqId1));
-		// std::string b = atomMapByAsymSeqAndAtom.at(make_tuple(asym2, seqId2, atomId2, authSeqId2));
-		// if (not(a.empty() or b.empty()))
-		// 	linkAtoms(a, b);
 	}
 
 	// then link all atoms in the compounds
-
-	for (auto c : compounds)
+	
+	for (auto c: compounds)
 	{
 		if (c == "HOH" or c == "H2O" or c == "WAT")
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "skipping water in bond map calculation" << std::endl;
 			continue;
 		}
 
-		auto bonded = [c, &compoundBondInfo](const Atom &a, const Atom &b)
+		auto bonded = [c, &compoundBondInfo](const Atom& a, const Atom& b)
 		{
 			auto label_a = a.labelAtomID();
 			auto label_b = b.labelAtomID();
@@ -341,17 +468,16 @@ BondMap::BondMap(const Structure &p)
 		};
 
 		// loop over poly_seq_scheme
-		for (auto r : db["pdbx_poly_seq_scheme"].find(cif::Key("mon_id") == c))
+		for (auto r: db["pdbx_poly_seq_scheme"].find(cif::Key("mon_id") == c))
 		{
 			std::string asymID;
 			int seqID;
 			cif::tie(asymID, seqID) = r.get("asym_id", "seq_id");
-
+			
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[&](auto &a)
-				{ return a.labelAsymID() == asymID and a.labelSeqID() == seqID; });
-
+				[&](auto& a) { return a.labelAsymID() == asymID and a.labelSeqID() == seqID; });
+			
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -363,16 +489,15 @@ BondMap::BondMap(const Structure &p)
 		}
 
 		// loop over pdbx_nonpoly_scheme
-		for (auto r : db["pdbx_nonpoly_scheme"].find(cif::Key("mon_id") == c))
+		for (auto r: db["pdbx_nonpoly_scheme"].find(cif::Key("mon_id") == c))
 		{
 			std::string asymID;
 			cif::tie(asymID) = r.get("asym_id");
-
+			
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[&](auto &a)
-				{ return a.labelAsymID() == asymID; });
-
+				[&](auto& a) { return a.labelAsymID() == asymID; });
+			
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -381,7 +506,7 @@ BondMap::BondMap(const Structure &p)
 					{
 						uint32_t ixa = index[rAtoms[i].id()];
 						uint32_t ixb = index[rAtoms[j].id()];
-
+						
 						bond.insert(key(ixa, ixb));
 					}
 				}
@@ -389,13 +514,15 @@ BondMap::BondMap(const Structure &p)
 		}
 
 		// loop over pdbx_branch_scheme
-		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"))
+		for (auto r: db["pdbx_branch_scheme"].find(cif::Key("mon_id") == c))
 		{
+			std::string asymID;
+			cif::tie(asymID) = r.get("asym_id");
+			
 			std::vector<Atom> rAtoms;
 			copy_if(atoms.begin(), atoms.end(), back_inserter(rAtoms),
-				[id = asym_id, nr = pdb_seq_num](const Atom &a)
-				{ return a.labelAsymID() == id and a.authSeqID() == nr; });
-
+				[&](auto& a) { return a.labelAsymID() == asymID; });
+			
 			for (uint32_t i = 0; i + 1 < rAtoms.size(); ++i)
 			{
 				for (uint32_t j = i + 1; j < rAtoms.size(); ++j)
@@ -404,31 +531,31 @@ BondMap::BondMap(const Structure &p)
 					{
 						uint32_t ixa = index[rAtoms[i].id()];
 						uint32_t ixb = index[rAtoms[j].id()];
-
+						
 						bond.insert(key(ixa, ixb));
 					}
 				}
 			}
 		}
 	}
-
+	
 	// start by creating an index for single bonds
-
-	std::multimap<uint32_t, uint32_t> b1_2;
-	for (auto &bk : bond)
+	
+	std::multimap<uint32_t,uint32_t> b1_2;
+	for (auto& bk: bond)
 	{
 		uint32_t a, b;
 		std::tie(a, b) = dekey(bk);
-
-		b1_2.insert({a, b});
-		b1_2.insert({b, a});
+		
+		b1_2.insert({ a, b });
+		b1_2.insert({ b, a });
 	}
-
-	std::multimap<uint32_t, uint32_t> b1_3;
+	
+	std::multimap<uint32_t,uint32_t> b1_3;
 	for (uint32_t i = 0; i < dim; ++i)
 	{
 		auto a = b1_2.equal_range(i);
-
+		
 		std::vector<uint32_t> s;
 		for (auto j = a.first; j != a.second; ++j)
 			s.push_back(j->second);
@@ -439,12 +566,12 @@ BondMap::BondMap(const Structure &p)
 			{
 				uint32_t x = s[si1];
 				uint32_t y = s[si2];
-
+				
 				if (isBonded(x, y))
 					continue;
-
-				b1_3.insert({x, y});
-				b1_3.insert({y, x});
+				
+				b1_3.insert({ x, y });
+				b1_3.insert({ y, x });
 			}
 		}
 	}
@@ -453,48 +580,48 @@ BondMap::BondMap(const Structure &p)
 	{
 		auto a1 = b1_2.equal_range(i);
 		auto a2 = b1_3.equal_range(i);
-
+		
 		for (auto ai1 = a1.first; ai1 != a1.second; ++ai1)
 		{
 			for (auto ai2 = a2.first; ai2 != a2.second; ++ai2)
 			{
 				uint32_t b1 = ai1->second;
 				uint32_t b2 = ai2->second;
-
+				
 				if (isBonded(b1, b2))
 					continue;
-
+				
 				bond_1_4.insert(key(b1, b2));
 			}
 		}
 	}
 }
 
-std::vector<std::string> BondMap::linked(const Atom &a) const
+std::vector<std::string> BondMap::linked(const Atom& a) const
 {
 	auto i = link.find(a.id());
-
+	
 	std::vector<std::string> result;
-
+	
 	if (i != link.end())
 		result = std::vector<std::string>(i->second.begin(), i->second.end());
 
 	return result;
 }
 
-std::vector<std::string> BondMap::atomIDsForCompound(const std::string &compoundID)
+std::vector<std::string> BondMap::atomIDsForCompound(const std::string& compoundID)
 {
 	std::vector<std::string> result;
 
-	auto *compound = mmcif::CompoundFactory::instance().create(compoundID);
+	auto* compound = mmcif::CompoundFactory::instance().create(compoundID);
 
 	if (compound == nullptr)
 		throw BondMapException("Missing bond information for compound " + compoundID);
 
-	for (auto &compAtom : compound->atoms())
+	for (auto& compAtom: compound->atoms())
 		result.push_back(compAtom.id);
 
 	return result;
 }
 
-} // namespace mmcif
+}

src/Cif2PDB.cpp

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

src/CifUtils.cpp

@@ -25,7 +25,10 @@
  */
 
 #include <atomic>
+#include <chrono>
 #include <cmath>
+#include <cstdio>
+#include <filesystem>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -34,6 +37,7 @@
 #include <regex>
 #include <sstream>
 #include <thread>
+#include <tuple>
 
 #if defined(_MSC_VER)
 #define TERM_WIDTH 80
@@ -45,7 +49,6 @@
 #include <boost/algorithm/string.hpp>
 
 #include "cif++/CifUtils.hpp"
-
 #include "revision.hpp"
 
 namespace ba = boost::algorithm;
@@ -91,12 +94,11 @@ const uint8_t kCharToLowerMap[256] =
 
 // --------------------------------------------------------------------
 
-bool iequals(std::string_view a, std::string_view b)
+bool iequals(const std::string &a, const std::string &b)
 {
 	bool result = a.length() == b.length();
-	for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end(); ++ai, ++bi)
-		result = kCharToLowerMap[uint8_t(*ai)] == kCharToLowerMap[uint8_t(*bi)];
-		// result = tolower(*ai) == tolower(*bi);
+	for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end() and bi != b.end(); ++ai, ++bi)
+		result = tolower(*ai) == tolower(*bi);
 	return result;
 }
 
@@ -109,7 +111,7 @@ bool iequals(const char *a, const char *b)
 	return result and *a == *b;
 }
 
-int icompare(std::string_view a, std::string_view b)
+int icompare(const std::string &a, const std::string &b)
 {
 	int d = 0;
 	auto ai = a.begin(), bi = b.begin();
@@ -162,7 +164,7 @@ std::string toLowerCopy(const std::string &s)
 
 // --------------------------------------------------------------------
 
-std::tuple<std::string, std::string> splitTagName(std::string_view tag)
+std::tuple<std::string, std::string> splitTagName(const std::string &tag)
 {
 	if (tag.empty())
 		throw std::runtime_error("empty tag");
@@ -643,6 +645,7 @@ void ProgressImpl::PrintProgress()
 
 	msg.append("| ");
 
+	//	const char		kSpinner[] = { '|', '/', '-', '\\' };
 	const char kSpinner[] = {' ', '.', 'o', 'O', '0', 'O', 'o', '.'};
 	const size_t kSpinnerCount = sizeof(kSpinner);
 
@@ -1199,97 +1202,63 @@ namespace cif
 
 // --------------------------------------------------------------------
 
-class ResourcePool
+std::map<std::string,std::filesystem::path> gLocalResources;
+std::filesystem::path gDataDir;
+
+void addDataDirectory(std::filesystem::path dataDir)
 {
-  public:
-	static ResourcePool &instance()
-	{
-		static std::unique_ptr<ResourcePool> sInstance(new ResourcePool);
-		return *sInstance;
-	}
-
-	void pushDir(fs::path dir)
-	{
-		std::error_code ec;
-
-		if (fs::exists(dir, ec) and not ec)
-			mDirs.push_front(dir);
-	}
-
-	void pushDir(const char *path)
-	{
-		if (path != nullptr)
-			pushDir(fs::path(path));
-	}
-
-	void pushAlias(const std::string &name, std::filesystem::path dataFile)
-	{
-		std::error_code ec;
-		if (not fs::exists(dataFile, ec) or ec)
-			throw std::runtime_error("Attempt to add a file resource for " + name + " that cannot be used (" + dataFile.string() + ") :" + ec.message());
-
-		mLocalResources[name] = dataFile;
-	}
-
-	std::unique_ptr<std::istream> load(fs::path name);
-
-  private:
-	ResourcePool();
-
-	std::unique_ptr<std::ifstream> open(fs::path &p)
-	{
-		std::unique_ptr<std::ifstream> result;
-
-		try
-		{
-			if (fs::exists(p))
-			{
-				std::unique_ptr<std::ifstream> file(new std::ifstream(p, std::ios::binary));
-				if (file->is_open())
-					result.reset(file.release());
-			}
-		}
-		catch (...) {}
-
-		return result;
-	}
-
-	std::map<std::string,std::filesystem::path> mLocalResources;
-	std::deque<fs::path> mDirs;
-};
-
-ResourcePool::ResourcePool()
-{
-#if defined(DATA_DIR)
-	pushDir(DATA_DIR);
-#endif
-
-	pushDir(getenv("LIBCIFPP_DATA_DIR"));
-
-	auto ccp4 = getenv("CCP4");
-	if (ccp4 != nullptr)
-		pushDir(fs::path(ccp4) / "share" / "libcifpp");
-
-#if defined(CACHE_DIR)
-	pushDir(CACHE_DIR);
-#endif
+	if (VERBOSE and not fs::exists(dataDir))
+		std::cerr << "The specified data directory " << dataDir << " does not exist" << std::endl;
+	gDataDir = dataDir;
 }
 
-std::unique_ptr<std::istream> ResourcePool::load(fs::path name)
+void addFileResource(const std::string &name, std::filesystem::path dataFile)
+{
+	if (not fs::exists(dataFile))
+		throw std::runtime_error("Attempt to add a file resource for " + name + " that does not exist: " + dataFile.string());
+
+	gLocalResources[name] = dataFile;
+}
+
+std::unique_ptr<std::istream> loadResource(std::filesystem::path name)
 {
 	std::unique_ptr<std::istream> result;
-	std::error_code ec;
 
 	fs::path p = name;
 
-	if (mLocalResources.count(name.string()))
-		result = open(mLocalResources[name.string()]);
-
-	for (auto di = mDirs.begin(); not result and di != mDirs.end(); ++di)
+	if (gLocalResources.count(name.string()))
 	{
-		auto p2 = *di / p;
-		if (fs::exists(p2, ec) and not ec)
-			result = open(p2);
+		std::unique_ptr<std::ifstream> file(new std::ifstream(gLocalResources[name.string()], std::ios::binary));
+		if (file->is_open())
+			result.reset(file.release());
+	}
+
+	if (not result and not fs::exists(p) and not gDataDir.empty())
+		p = gDataDir / name;
+
+#if defined(CACHE_DIR)
+	if (not result and not fs::exists(p))
+	{
+		auto p2 = fs::path(CACHE_DIR) / p;
+		if (fs::exists(p2))
+			swap(p, p2);
+	}
+#endif
+
+#if defined(DATA_DIR)
+	if (not result and not fs::exists(p))
+	{
+		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));
+		if (file->is_open())
+			result.reset(file.release());
 	}
 
 	if (not result and gResourceData)
@@ -1299,24 +1268,7 @@ std::unique_ptr<std::istream> ResourcePool::load(fs::path name)
 			result.reset(new mrsrc::istream(rsrc));
 	}
 
-	return result;	
-}
-
-// --------------------------------------------------------------------
-
-void addDataDirectory(std::filesystem::path dataDir)
-{
-	ResourcePool::instance().pushDir(dataDir);
-}
-
-void addFileResource(const std::string &name, std::filesystem::path dataFile)
-{
-	ResourcePool::instance().pushAlias(name, dataFile);
-}
-
-std::unique_ptr<std::istream> loadResource(std::filesystem::path name)
-{
-	return ResourcePool::instance().load(name);
+	return result;
 }
 
 } // namespace cif

src/CifValidator.cpp

@@ -24,39 +24,32 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include <fstream>
-#include <filesystem>
-
 #include <boost/algorithm/string.hpp>
-#include <boost/iostreams/filtering_stream.hpp>
-#include <boost/iostreams/filter/gzip.hpp>
 
 #include "cif++/Cif++.hpp"
 #include "cif++/CifParser.hpp"
 #include "cif++/CifValidator.hpp"
 
 namespace ba = boost::algorithm;
-namespace fs = std::filesystem;
-namespace io = boost::iostreams;
 
 extern int VERBOSE;
 
 namespace cif
 {
 
-ValidationError::ValidationError(const std::string &msg)
+ValidationError::ValidationError(const std::string& msg)
 	: mMsg(msg)
 {
 }
 
-ValidationError::ValidationError(const std::string &cat, const std::string &item, const std::string &msg)
+ValidationError::ValidationError(const std::string& cat, const std::string& item, const std::string& msg)
 	: mMsg("When validating _" + cat + '.' + item + ": " + msg)
 {
 }
 
 // --------------------------------------------------------------------
 
-DDL_PrimitiveType mapToPrimitiveType(std::string_view s)
+DDL_PrimitiveType mapToPrimitiveType(const std::string& s)
 {
 	DDL_PrimitiveType result;
 	if (iequals(s, "char"))
@@ -72,10 +65,10 @@ DDL_PrimitiveType mapToPrimitiveType(std::string_view s)
 
 // --------------------------------------------------------------------
 
-int ValidateType::compare(const char *a, const char *b) const
+int ValidateType::compare(const char* a, const char* b) const
 {
 	int result = 0;
-
+	
 	if (*a == 0)
 		result = *b == 0 ? 0 : -1;
 	else if (*b == 0)
@@ -90,7 +83,7 @@ int ValidateType::compare(const char *a, const char *b) const
 				{
 					double da = strtod(a, nullptr);
 					double db = strtod(b, nullptr);
-
+					
 					auto d = da - db;
 					if (std::abs(d) > std::numeric_limits<double>::epsilon())
 					{
@@ -101,13 +94,13 @@ int ValidateType::compare(const char *a, const char *b) const
 					}
 					break;
 				}
-
+				
 				case DDL_PrimitiveType::UChar:
 				case DDL_PrimitiveType::Char:
 				{
 					// CIF is guaranteed to have ascii only, therefore this primitive code will do
 					// also, we're collapsing spaces
-
+					
 					auto ai = a, bi = b;
 					for (;;)
 					{
@@ -122,7 +115,7 @@ int ValidateType::compare(const char *a, const char *b) const
 							result = 1;
 							break;
 						}
-
+						
 						char ca = *ai;
 						char cb = *bi;
 
@@ -131,12 +124,12 @@ int ValidateType::compare(const char *a, const char *b) const
 							ca = tolower(ca);
 							cb = tolower(cb);
 						}
-
+						
 						result = ca - cb;
-
+						
 						if (result != 0)
 							break;
-
+						
 						if (ca == ' ')
 						{
 							while (ai[1] == ' ')
@@ -144,21 +137,21 @@ int ValidateType::compare(const char *a, const char *b) const
 							while (bi[1] == ' ')
 								++bi;
 						}
-
+						
 						++ai;
 						++bi;
 					}
-
+					
 					break;
 				}
 			}
 		}
-		catch (const std::invalid_argument &ex)
+		catch (const std::invalid_argument& ex)
 		{
 			result = 1;
 		}
 	}
-
+	
 	return result;
 }
 
@@ -172,13 +165,13 @@ int ValidateType::compare(const char *a, const char *b) const
 //
 //	if (mType == nullptr and parent != nullptr)
 //		mType = parent->mType;
-//
+//		
 //	if (parent != nullptr)
 //	{
 //		mLinked.push_back({parent, parentItem, childItem});
 //
 //		parent->mChildren.insert(this);
-////
+////	
 ////		if (mCategory->mKeys == std::vector<std::string>{mTag})
 ////			parent->mForeignKeys.insert(this);
 //	}
@@ -201,7 +194,7 @@ void ValidateItem::operator()(std::string value) const
 
 // --------------------------------------------------------------------
 
-void ValidateCategory::addItemValidator(ValidateItem &&v)
+void ValidateCategory::addItemValidator(ValidateItem&& v)
 {
 	if (v.mMandatory)
 		mMandatoryFields.insert(v.mTag);
@@ -213,10 +206,10 @@ void ValidateCategory::addItemValidator(ValidateItem &&v)
 		std::cout << "Could not add validator for item " << v.mTag << " to category " << mName << std::endl;
 }
 
-const ValidateItem *ValidateCategory::getValidatorForItem(std::string_view tag) const
+const ValidateItem* ValidateCategory::getValidatorForItem(std::string tag) const
 {
-	const ValidateItem *result = nullptr;
-	auto i = mItemValidators.find(ValidateItem{std::string(tag)});
+	const ValidateItem* result = nullptr;
+	auto i = mItemValidators.find(ValidateItem{tag});
 	if (i != mItemValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -226,29 +219,26 @@ const ValidateItem *ValidateCategory::getValidatorForItem(std::string_view tag)
 
 // --------------------------------------------------------------------
 
-Validator::Validator(std::string_view name, std::istream &is)
-	: mName(name)
+Validator::Validator()
 {
-	DictParser p(*this, is);
-	p.loadDictionary();
 }
 
 Validator::~Validator()
 {
 }
 
-void Validator::addTypeValidator(ValidateType &&v)
+void Validator::addTypeValidator(ValidateType&& v)
 {
 	auto r = mTypeValidators.insert(std::move(v));
 	if (not r.second and VERBOSE > 4)
 		std::cout << "Could not add validator for type " << v.mName << std::endl;
 }
 
-const ValidateType *Validator::getValidatorForType(std::string_view typeCode) const
+const ValidateType* Validator::getValidatorForType(std::string typeCode) const
 {
-	const ValidateType *result = nullptr;
-
-	auto i = mTypeValidators.find(ValidateType{std::string(typeCode), DDL_PrimitiveType::Char, boost::regex()});
+	const ValidateType* result = nullptr;
+	
+	auto i = mTypeValidators.find(ValidateType{ typeCode, DDL_PrimitiveType::Char, boost::regex() });
 	if (i != mTypeValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -256,17 +246,17 @@ const ValidateType *Validator::getValidatorForType(std::string_view typeCode) co
 	return result;
 }
 
-void Validator::addCategoryValidator(ValidateCategory &&v)
+void Validator::addCategoryValidator(ValidateCategory&& v)
 {
 	auto r = mCategoryValidators.insert(std::move(v));
 	if (not r.second and VERBOSE > 4)
 		std::cout << "Could not add validator for category " << v.mName << std::endl;
 }
 
-const ValidateCategory *Validator::getValidatorForCategory(std::string_view category) const
+const ValidateCategory* Validator::getValidatorForCategory(std::string category) const
 {
-	const ValidateCategory *result = nullptr;
-	auto i = mCategoryValidators.find(ValidateCategory{std::string(category)});
+	const ValidateCategory* result = nullptr;
+	auto i = mCategoryValidators.find(ValidateCategory{category});
 	if (i != mCategoryValidators.end())
 		result = &*i;
 	else if (VERBOSE > 4)
@@ -274,16 +264,16 @@ const ValidateCategory *Validator::getValidatorForCategory(std::string_view cate
 	return result;
 }
 
-ValidateItem *Validator::getValidatorForItem(std::string_view tag) const
+ValidateItem* Validator::getValidatorForItem(std::string tag) const
 {
-	ValidateItem *result = nullptr;
-
+	ValidateItem* result = nullptr;
+	
 	std::string cat, item;
 	std::tie(cat, item) = splitTagName(tag);
 
-	auto *cv = getValidatorForCategory(cat);
+	auto* cv = getValidatorForCategory(cat);
 	if (cv != nullptr)
-		result = const_cast<ValidateItem *>(cv->getValidatorForItem(item));
+		result = const_cast<ValidateItem*>(cv->getValidatorForItem(item));
 
 	if (result == nullptr and VERBOSE > 4)
 		std::cout << "No validator for item " << tag << std::endl;
@@ -291,15 +281,15 @@ ValidateItem *Validator::getValidatorForItem(std::string_view tag) const
 	return result;
 }
 
-void Validator::addLinkValidator(ValidateLink &&v)
+void Validator::addLinkValidator(ValidateLink&& v)
 {
 	assert(v.mParentKeys.size() == v.mChildKeys.size());
 	if (v.mParentKeys.size() != v.mChildKeys.size())
 		throw std::runtime_error("unequal number of keys for parent and child in link");
-
+	
 	auto pcv = getValidatorForCategory(v.mParentCategory);
 	auto ccv = getValidatorForCategory(v.mChildCategory);
-
+	
 	if (pcv == nullptr)
 		throw std::runtime_error("unknown parent category " + v.mParentCategory);
 
@@ -309,130 +299,53 @@ void Validator::addLinkValidator(ValidateLink &&v)
 	for (size_t i = 0; i < v.mParentKeys.size(); ++i)
 	{
 		auto piv = pcv->getValidatorForItem(v.mParentKeys[i]);
-
+	
 		if (piv == nullptr)
 			throw std::runtime_error("unknown parent tag _" + v.mParentCategory + '.' + v.mParentKeys[i]);
 
 		auto civ = ccv->getValidatorForItem(v.mChildKeys[i]);
 		if (civ == nullptr)
 			throw std::runtime_error("unknown child tag _" + v.mChildCategory + '.' + v.mChildKeys[i]);
-
+		
 		if (civ->mType == nullptr and piv->mType != nullptr)
-			const_cast<ValidateItem *>(civ)->mType = piv->mType;
-	}
-
+			const_cast<ValidateItem*>(civ)->mType = piv->mType;
+	}		
+	
 	mLinkValidators.emplace_back(std::move(v));
 }
 
-std::vector<const ValidateLink *> Validator::getLinksForParent(std::string_view category) const
+std::vector<const ValidateLink*> Validator::getLinksForParent(const std::string& category) const
 {
-	std::vector<const ValidateLink *> result;
+	std::vector<const ValidateLink*> result;
 
-	for (auto &l : mLinkValidators)
+	for (auto& l: mLinkValidators)
 	{
 		if (l.mParentCategory == category)
 			result.push_back(&l);
 	}
-
+	
 	return result;
 }
 
-std::vector<const ValidateLink *> Validator::getLinksForChild(std::string_view category) const
+std::vector<const ValidateLink*> Validator::getLinksForChild(const std::string& category) const
 {
-	std::vector<const ValidateLink *> result;
+	std::vector<const ValidateLink*> result;
 
-	for (auto &l : mLinkValidators)
+	for (auto& l: mLinkValidators)
 	{
 		if (l.mChildCategory == category)
 			result.push_back(&l);
 	}
-
+	
 	return result;
 }
 
-void Validator::reportError(const std::string &msg, bool fatal) const
+void Validator::reportError(const std::string& msg, bool fatal)
 {
 	if (mStrict or fatal)
 		throw ValidationError(msg);
-	else if (VERBOSE > 0)
+	else if (VERBOSE)
 		std::cerr << msg << std::endl;
 }
 
-// --------------------------------------------------------------------
-
-ValidatorFactory ValidatorFactory::sInstance;
-
-ValidatorFactory::ValidatorFactory()
-{
 }
-
-const Validator &ValidatorFactory::operator[](std::string_view dictionary)
-{
-	std::lock_guard lock(mMutex);
-
-	for (auto &validator : mValidators)
-	{
-		if (iequals(validator.mName, dictionary))
-			return validator;
-	}
-
-	// not found, add it
-
-	// too bad clang version 10 did not have a constructor for fs::path that accepts a std::string_view
-	fs::path dict_name(dictionary.data(), dictionary.data() + dictionary.length());
-
-	auto data = loadResource(dict_name);
-
-	if (not data and dict_name.extension().string() != ".dic")
-		data = loadResource(dict_name.parent_path() / (dict_name.filename().string() + ".dic"));
-
-	if (data)
-		mValidators.emplace_back(dictionary, *data);
-	else
-	{
-		std::error_code ec;
-
-		// might be a compressed dictionary on disk
-		fs::path p = dict_name;
-		if (p.extension() == ".dic")
-			p = p.parent_path() / (p.filename().string() + ".gz");
-		else
-			p = p.parent_path() / (p.filename().string() + ".dic.gz");
-
-#if defined(CACHE_DIR) and defined(DATA_DIR)
-		if (not fs::exists(p, ec) or ec)
-		{
-			for (const char *dir : {CACHE_DIR, DATA_DIR})
-			{
-				auto p2 = fs::path(dir) / p;
-				if (fs::exists(p2, ec) and not ec)
-				{
-					swap(p, p2);
-					break;
-				}
-			}
-		}
-#endif
-
-		if (fs::exists(p, ec) and not ec)
-		{
-			std::ifstream file(p, std::ios::binary);
-			if (not file.is_open())
-				throw std::runtime_error("Could not open dictionary (" + p.string() + ")");
-
-			io::filtering_stream<io::input> in;
-			in.push(io::gzip_decompressor());
-			in.push(file);
-
-			mValidators.emplace_back(dictionary, in);
-		}
-		else
-			throw std::runtime_error("Dictionary not found or defined (" + dict_name.string() + ")");
-	}
-
-	assert(iequals(mValidators.back().mName, dictionary));
-
-	return mValidators.back();
-}
-
-} // namespace cif

src/Compound.cpp

@@ -128,8 +128,6 @@ Compound::Compound(cif::Datablock &db)
 	// The name should not contain newline characters since that triggers validation errors later on
 	ba::replace_all(mName, "\n", "");
 
-	mGroup = "non-polymer";
-
 	auto &chemCompAtom = db["chem_comp_atom"];
 	for (auto row : chemCompAtom)
 	{
@@ -153,11 +151,10 @@ Compound::Compound(cif::Datablock &db)
 	}
 }
 
-Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group)
+Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type)
 	: mID(id)
 	, mName(name)
 	, mType(type)
-	, mGroup(group)
 {
 	auto &chemCompAtom = db["chem_comp_atom"];
 	for (auto row : chemCompAtom)
@@ -193,7 +190,7 @@ Compound::Compound(cif::Datablock &db, const std::string &id, const std::string
 			bond.type = BondType::delo;
 		else
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Unimplemented chem_comp_bond.type " << btype << " in " << id << std::endl;
 			bond.type = BondType::sing;
 		}
@@ -274,7 +271,7 @@ class CompoundFactoryImpl : public std::enable_shared_from_this<CompoundFactoryI
   public:
 	CompoundFactoryImpl(std::shared_ptr<CompoundFactoryImpl> next);
 
-	CompoundFactoryImpl(const fs::path &file, std::shared_ptr<CompoundFactoryImpl> next);
+	CompoundFactoryImpl(const std::filesystem::path &file, std::shared_ptr<CompoundFactoryImpl> next);
 
 	virtual ~CompoundFactoryImpl()
 	{
@@ -368,7 +365,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(std::shared_ptr<CompoundFactoryImpl> ne
 		mKnownBases.insert(key);
 }
 
-CompoundFactoryImpl::CompoundFactoryImpl(const fs::path &file, std::shared_ptr<CompoundFactoryImpl> next)
+CompoundFactoryImpl::CompoundFactoryImpl(const std::filesystem::path &file, std::shared_ptr<CompoundFactoryImpl> next)
 	: mNext(next)
 {
 	cif::File cifFile(file);
@@ -403,7 +400,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const fs::path &file, std::shared_ptr<C
 
 			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") or cif::iequals(group, "p-peptide"))
+			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
 				type = "L-peptide linking";
 			else if (cif::iequals(group, "DNA"))
 				type = "DNA linking";
@@ -414,7 +411,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const fs::path &file, std::shared_ptr<C
 
 			auto &db = cifFile["comp_" + id];
 
-			mCompounds.push_back(new Compound(db, id, name, type, group));
+			mCompounds.push_back(new Compound(db, id, name, type));
 		}
 	}
 	else
@@ -520,7 +517,7 @@ Compound *CCDCompoundFactoryImpl::create(const std::string &id)
 		}
 	}
 
-	if (result == nullptr and cif::VERBOSE > 0)
+	if (result == nullptr and cif::VERBOSE)
 		std::cerr << "Could not locate compound " << id << " in the CCD components file" << std::endl;
 
 	return result;
@@ -611,7 +608,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") or cif::iequals(group, "p-peptide"))
+			else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
 				type = "L-peptide linking";
 			else if (cif::iequals(group, "DNA"))
 				type = "DNA linking";
@@ -620,7 +617,7 @@ Compound *CCP4CompoundFactoryImpl::create(const std::string &id)
 			else
 				type = "non-polymer";
 
-			mCompounds.push_back(new Compound(db, id, name, type, group));
+			mCompounds.push_back(new Compound(db, id, name, type));
 			result = mCompounds.back();
 		}
 	}
@@ -645,13 +642,13 @@ CompoundFactory::CompoundFactory()
 	auto ccd = cif::loadResource("components.cif");
 	if (ccd)
 		mImpl.reset(new CCDCompoundFactoryImpl(mImpl));
-	else if (cif::VERBOSE > 0)
+	else if (cif::VERBOSE)
 		std::cerr << "CCD components.cif file was not found" << std::endl;
 
 	const char *clibd_mon = getenv("CLIBD_MON");
 	if (clibd_mon != nullptr and fs::is_directory(clibd_mon))
 		mImpl.reset(new CCP4CompoundFactoryImpl(clibd_mon));
-	else if (cif::VERBOSE > 0)
+	else if (cif::VERBOSE)
 		std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;
 
 }
@@ -684,7 +681,7 @@ void CompoundFactory::clear()
 		sInstance.reset();
 }
 
-void CompoundFactory::setDefaultDictionary(const fs::path &inDictFile)
+void CompoundFactory::setDefaultDictionary(const std::filesystem::path &inDictFile)
 {
 	if (not fs::exists(inDictFile))
 		throw std::runtime_error("file not found: " + inDictFile.string());
@@ -695,13 +692,12 @@ void CompoundFactory::setDefaultDictionary(const fs::path &inDictFile)
 	}
 	catch (const std::exception &)
 	{
-		if (cif::VERBOSE >= 0)
-			std::cerr << "Error loading dictionary " << inDictFile << std::endl;
+		std::cerr << "Error loading dictionary " << inDictFile << std::endl;
 		throw;
 	}
 }
 
-void CompoundFactory::pushDictionary(const fs::path &inDictFile)
+void CompoundFactory::pushDictionary(const std::filesystem::path &inDictFile)
 {
 	if (not fs::exists(inDictFile))
 		throw std::runtime_error("file not found: " + inDictFile.string());
@@ -716,8 +712,7 @@ void CompoundFactory::pushDictionary(const fs::path &inDictFile)
 	}
 	catch (const std::exception &)
 	{
-		if (cif::VERBOSE >= 0)
-			std::cerr << "Error loading dictionary " << inDictFile << std::endl;
+		std::cerr << "Error loading dictionary " << inDictFile << std::endl;
 		throw;
 	}
 }

src/PDB2CifRemark3.cpp

@@ -1320,7 +1320,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 
 	if (line != "REFINEMENT.")
 	{
-		if (cif::VERBOSE > 0)
+		if (cif::VERBOSE)
 			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 > 0)
+		if (cif::VERBOSE)
 			std::cerr << "Expected valid PROGRAM line in REMARK 3" << std::endl;
 		return false;
 	}
@@ -1367,9 +1367,8 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 		}
 		catch(const std::exception& e)
 		{
-			if (cif::VERBOSE >= 0)
-				std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
-						  << e.what() << '\n';
+			std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
+					  << e.what() << '\n';
 			score = 0;
 		}
 
@@ -1412,7 +1411,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 > 0)
+		else if (cif::VERBOSE)
 			std::cerr << "Skipping unknown program (" << program << ") in REMARK 3" << std::endl;
 	}
 
@@ -1421,8 +1420,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 	bool guessProgram = scores.empty() or scores.front().score < 0.9f;;
 	if (guessProgram)
 	{
-		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;
+		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));
@@ -1446,7 +1444,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
 
 		auto& best = scores.front();
 
-		if (cif::VERBOSE > 0)
+		if (cif::VERBOSE)
 			std::cerr << "Choosing " << best.parser->program() << " version '" << best.parser->version() << "' as refinement program. Score = " << best.score << std::endl;
 
 		auto& software = db["software"];

src/Point.cpp

@@ -28,248 +28,11 @@
 #include <valarray>
 
 #include "cif++/Point.hpp"
+#include "cif++/Matrix.hpp"
 
 namespace mmcif
 {
 
-// --------------------------------------------------------------------
-// We're using expression templates here
-
-template <typename M>
-class MatrixExpression
-{
-  public:
-	uint32_t dim_m() const { return static_cast<const M&>(*this).dim_m(); }
-	uint32_t dim_n() const { return static_cast<const M&>(*this).dim_n(); }
-
-	double &operator()(uint32_t i, uint32_t j)
-	{
-		return static_cast<M&>(*this).operator()(i, j);
-	}
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return static_cast<const M&>(*this).operator()(i, j);
-	}
-};
-
-// --------------------------------------------------------------------
-// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
-// element m i,j is mapped to [i * n + j] and thus storage is row major
-
-class Matrix : public MatrixExpression<Matrix>
-{
-  public:
-	template <typename M2>
-	Matrix(const MatrixExpression<M2> &m)
-		: m_m(m.dim_m())
-		, m_n(m.dim_n())
-		, m_data(m_m * m_n)
-	{
-		for (uint32_t i = 0; i < m_m; ++i)
-		{
-			for (uint32_t j = 0; j < m_n; ++j)
-				operator()(i, j) = m(i, j);
-		}
-	}
-
-	Matrix(size_t m, size_t n, double v = 0)
-		: m_m(m)
-		, m_n(n)
-		, m_data(m_m * m_n)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	Matrix() = default;
-	Matrix(Matrix &&m) = default;
-	Matrix(const Matrix &m) = default;
-	Matrix &operator=(Matrix &&m) = default;
-	Matrix &operator=(const Matrix &m) = default;
-
-	uint32_t dim_m() const { return m_m; }
-	uint32_t dim_n() const { return m_n; }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-	double &operator()(uint32_t i, uint32_t j)
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-  private:
-	uint32_t m_m = 0, m_n = 0;
-	std::vector<double> m_data;
-};
-
-// --------------------------------------------------------------------
-
-class SymmetricMatrix : public MatrixExpression<SymmetricMatrix>
-{
-  public:
-	SymmetricMatrix(uint32_t n, double v = 0)
-		: m_n(n)
-		, m_data((m_n * (m_n + 1)) / 2)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	SymmetricMatrix() = default;
-	SymmetricMatrix(SymmetricMatrix &&m) = default;
-	SymmetricMatrix(const SymmetricMatrix &m) = default;
-	SymmetricMatrix &operator=(SymmetricMatrix &&m) = default;
-	SymmetricMatrix &operator=(const SymmetricMatrix &m) = default;
-
-	uint32_t dim_m() const { return m_n; }
-	uint32_t dim_n() const { return m_n; }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return i < j
-				? m_data[(j * (j + 1)) / 2 + i]
-				: m_data[(i * (i + 1)) / 2 + j];
-	}
-
-	double &operator()(uint32_t i, uint32_t j)
-	{
-		if (i > j)
-			std::swap(i, j);
-		assert(j < m_n);
-		return m_data[(j * (j + 1)) / 2 + i];
-	}
-
-  private:
-	uint32_t m_n;
-	std::vector<double> m_data;
-};
-
-class IdentityMatrix : public MatrixExpression<IdentityMatrix>
-{
-  public:
-	IdentityMatrix(uint32_t n)
-		: m_n(n)
-	{
-	}
-
-	uint32_t dim_m() const { return m_n; }
-	uint32_t dim_n() const { return m_n; }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return i == j ? 1 : 0;
-	}
-
-  private:
-	uint32_t m_n;
-};
-
-// --------------------------------------------------------------------
-// matrix functions, implemented as expression templates
-
-template<typename M1, typename M2>
-class MatrixSubtraction : public MatrixExpression<MatrixSubtraction<M1, M2>>
-{
-  public:
-	MatrixSubtraction(const M1 &m1, const M2 &m2)
-		: m_m1(m1), m_m2(m2)
-	{
-		assert(m_m1.dim_m() == m_m2.dim_m());
-		assert(m_m1.dim_n() == m_m2.dim_n());
-	}
-
-	uint32_t dim_m() const { return m_m1.dim_m(); }
-	uint32_t dim_n() const { return m_m1.dim_n(); }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return m_m1(i, j) - m_m2(i, j);
-	}
-
-  private:
-	const M1 &m_m1;
-	const M2 &m_m2;
-};
-
-template<typename M1, typename M2>
-MatrixSubtraction<M1, M2> operator-(const MatrixExpression<M1> &m1, const MatrixExpression<M2> &m2)
-{
-	return MatrixSubtraction(*static_cast<const M1*>(&m1), *static_cast<const M2*>(&m2));
-}
-
-template<typename M>
-class MatrixMultiplication : public MatrixExpression<MatrixMultiplication<M>>
-{
-  public:
-	MatrixMultiplication(const M &m, double v)
-		: m_m(m), m_v(v)
-	{
-	}
-
-	uint32_t dim_m() const { return m_m.dim_m(); }
-	uint32_t dim_n() const { return m_m.dim_n(); }
-
-	double operator()(uint32_t i, uint32_t j) const
-	{
-		return m_m(i, j) * m_v;
-	}
-
-  private:
-	const M &m_m;
-	double m_v;
-};
-
-template<typename M>
-MatrixMultiplication<M> operator*(const MatrixExpression<M> &m, double v)
-{
-	return MatrixMultiplication(*static_cast<const M*>(&m), v);
-}
-
-// --------------------------------------------------------------------
-
-template<class M1>
-Matrix Cofactors(const M1& m)
-{
-	Matrix cf(m.dim_m(), m.dim_m());
-
-    const size_t ixs[4][3] =
-    {
-        { 1, 2, 3 },
-        { 0, 2, 3 },
-        { 0, 1, 3 },
-        { 0, 1, 2 }
-    };
-
-    for (size_t x = 0; x < 4; ++x)
-    {
-        const size_t* ix = ixs[x];
-
-        for (size_t y = 0; y < 4; ++y)
-        {
-            const size_t* iy = ixs[y];
-
-            cf(x, y) =
-                m(ix[0], iy[0]) * m(ix[1], iy[1]) * m(ix[2], iy[2]) +
-                m(ix[0], iy[1]) * m(ix[1], iy[2]) * m(ix[2], iy[0]) +
-                m(ix[0], iy[2]) * m(ix[1], iy[0]) * m(ix[2], iy[1]) -
-                m(ix[0], iy[2]) * m(ix[1], iy[1]) * m(ix[2], iy[0]) -
-                m(ix[0], iy[1]) * m(ix[1], iy[0]) * m(ix[2], iy[2]) -
-                m(ix[0], iy[0]) * m(ix[1], iy[2]) * m(ix[2], iy[1]);
-			
-			if ((x + y) % 2 == 1)
-				cf(x, y) *= -1;
-        }
-    }
-
-	return cf;
-}
-
 // --------------------------------------------------------------------
 
 Quaternion Normalize(Quaternion q)
@@ -295,27 +58,13 @@ 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 * std::acos(q.R_component_1());
+	double angle = 2 * acos(q.R_component_1());
 	angle = angle * 180 / kPI;
 
 	// axis:
@@ -353,14 +102,14 @@ Point CenterPoints(std::vector<Point>& Points)
 	return t;
 }
 
-Point Centroid(const std::vector<Point>& pts)
+Point Centroid(std::vector<Point>& Points)
 {
 	Point result;
 	
-	for (auto &pt : pts)
+	for (Point& pt : Points)
 		result += pt;
 	
-	result /= static_cast<float>(pts.size());
+	result /= static_cast<float>(Points.size());
 	
 	return result;
 }
@@ -426,7 +175,7 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
 Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
 {
 	// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
-	Matrix M(3, 3, 0);
+	Matrix<double> M(3, 3, 0);
 
 	for (uint32_t i = 0; i < pa.size(); ++i)
 	{
@@ -439,7 +188,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 	}
 	
 	// Now calculate N, a symmetric 4x4 Matrix
-	SymmetricMatrix N(4);
+	SymmetricMatrix<double> N(4);
 	
 	N(0, 0) =  M(0, 0) + M(1, 1) + M(2, 2);
 	N(0, 1) =  M(1, 2) - M(2, 1);
@@ -476,7 +225,6 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 					M(1, 2) * M(2, 0) * M(0, 1) +
 					M(2, 1) * M(1, 0) * M(0, 2));
 	
-	// E is the determinant of N:
 	double E = 
 		(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
 		(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
@@ -486,22 +234,47 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 		(N(0,2) * N(1,3) - N(2,1) * N(0,3)) * (N(0,2) * N(1,3) - N(2,1) * N(0,3));
 	
 	// solve quartic
-	double lambda = LargestDepressedQuarticSolution(C, D, E);
+	double lm = LargestDepressedQuarticSolution(C, D, E);
 	
 	// calculate t = (N - λI)
-	Matrix t = N - IdentityMatrix(4) * lambda;
+	Matrix<double> li = IdentityMatrix<double>(4) * lm;
+	Matrix<double> t = N - li;
 	
 	// calculate a Matrix of cofactors for t
-	Matrix cf = Cofactors(t);
+	Matrix<double> cf(4, 4);
 
-	int maxR = 0;
-	for (int r = 1; r < 4; ++r)
+	const uint32_t ixs[4][3] =
 	{
-		if (std::abs(cf(r, 0)) > std::abs(cf(maxR, 0)))
+		{ 1, 2, 3 },
+		{ 0, 2, 3 },
+		{ 0, 1, 3 },
+		{ 0, 1, 2 }
+	};
+
+	uint32_t maxR = 0;
+	for (uint32_t r = 0; r < 4; ++r)
+	{
+		const uint32_t* ir = ixs[r];
+		
+		for (uint32_t c = 0; c < 4; ++c)
+		{
+			const uint32_t* ic = ixs[c];
+
+			cf(r, c) =
+				t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
+				t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
+				t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
+				t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
+				t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
+				t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
+		}
+		
+		if (r > maxR and cf(r, 0) > cf(maxR, 0))
 			maxR = r;
 	}
 	
-	Quaternion q(cf(maxR, 0), cf(maxR, 1), cf(maxR, 2), cf(maxR, 3));
+	// NOTE the negation of the y here, why? Maybe I swapped r/c above?
+	Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
 	q = Normalize(q);
 	
 	return q;

src/Symmetry.cpp

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

src/TlsParser.cpp

@@ -248,7 +248,7 @@ struct TLSSelectionNot : public TLSSelection
 		for (auto& r: residues)
 			r.selected = not r.selected;
 
-		if (cif::VERBOSE > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 > 0)
+		if (cif::VERBOSE)
 		{
 			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 and cif::VERBOSE > 0)
+	if (extraParenthesis)
 		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 and cif::VERBOSE > 0)
+			if (m_lookahead == pt_EOLN)
 				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 > 0 and (icode_from or icode_to))
+					if (cif::VERBOSE and (icode_from or icode_to))
 						std::cerr << "Warning, ignoring insertion codes" << std::endl;
 					
 					result.reset(new TLSSelectionRangeSeq(from, to));
@@ -1231,8 +1231,7 @@ TLSSelectionPtr TLSSelectionParserImplBuster::ParseGroup()
 				std::tie(chain2, seqNr2) = ParseAtom();
 				if (chain1 != chain2)
 				{
-					if (cif::VERBOSE > 0)
-						std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
+					std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
 					
 					TLSSelectionPtr sc1(new TLSSelectionChain(chain1));
 					TLSSelectionPtr sr1(new TLSSelectionRangeSeq(seqNr1, kResidueNrWildcard));
@@ -1290,7 +1289,7 @@ std::tuple<std::string,int> TLSSelectionParserImplBuster::ParseAtom()
 			Match(':');
 			std::string atom = m_value_s;
 			
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Warning: ignoring atom ID '" << atom << "' in TLS selection" << std::endl;
 			
 			Match(bt_IDENT);
@@ -1811,8 +1810,7 @@ class TLSSelectionParser
 		}
 		catch (const std::exception& ex)
 		{
-			if (cif::VERBOSE >= 0)
-				std::cerr << "ParseError: " << ex.what() << std::endl;
+			std::cerr << "ParseError: " << ex.what() << std::endl;
 		}
 		
 		return result;
@@ -1836,14 +1834,14 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str
 
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}
 		
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to PHENIX" << std::endl;
 			result = phenix.Parse(selection);
 		}
@@ -1854,35 +1852,35 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str
 
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to BUSTER" << std::endl;
 			result = buster.Parse(selection);
 		}
 
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}
 	}
 	else
 	{
-		if (cif::VERBOSE > 0)
+		if (cif::VERBOSE)
 			std::cerr << "No known program specified, trying PHENIX" << std::endl;
 
 		result = phenix.Parse(selection);
 
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to BUSTER" << std::endl;
 			result = buster.Parse(selection);
 		}
 
 		if (not result)
 		{
-			if (cif::VERBOSE > 0)
+			if (cif::VERBOSE)
 				std::cerr << "Falling back to old BUSTER" << std::endl;
 			result = busterOld.Parse(selection);
 		}

src/revision.hpp.in

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

test/pdb2cif-test.cpp

@@ -18,6 +18,7 @@ int main(int argc, char* argv[])
 	desc.add_options()
 		("input,i",		po::value<std::string>(),	"Input file")
 		("help,h",									"Display help message")
+		("version",									"Print version")
 		("verbose,v",								"Verbose output")
 		("debug,d",		po::value<int>(),			"Debug level (for even more verbose output)");
 
@@ -28,6 +29,12 @@ int main(int argc, char* argv[])
 	po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
 	po::notify(vm);
 
+	if (vm.count("version"))
+	{
+		std::cout << argv[0] << " version " PACKAGE_VERSION << std::endl;
+		exit(0);
+	}
+
 	if (vm.count("help") or vm.count("input") == 0)
 	{
 		std::cerr << desc << std::endl;

test/rename-compound-test.cpp

@@ -20,11 +20,8 @@ int main(int argc, char* argv[])
 		if (argc == 3)
 			testdir = argv[2];
 
-		if (std::filesystem::exists(testdir / ".." / "data" / "ccd-subset.cif"))
-			cif::addFileResource("components.cif", testdir / ".." / "data" / "ccd-subset.cif");
-
-		if (std::filesystem::exists(testdir / ".." / "rsrc" / "mmcif_pdbx_v50.dic"))
-			cif::addFileResource("mmcif_pdbx_v50.dic", testdir / ".." / "rsrc" / "mmcif_pdbx_v50.dic");
+		if (std::filesystem::exists(testdir / ".."/"data"/"ccd-subset.cif"))
+			cif::addFileResource("components.cif", testdir / ".."/"data"/"ccd-subset.cif");
 
 		mmcif::CompoundFactory::instance().pushDictionary(testdir / "REA.cif");
 		mmcif::CompoundFactory::instance().pushDictionary(testdir / "RXA.cif");
@@ -32,12 +29,12 @@ int main(int argc, char* argv[])
 		mmcif::File f(testdir / ".."/"examples"/"1cbs.cif.gz");
 		mmcif::Structure structure(f);
 
-		auto &res = structure.getResidue("B");
+		auto &res = structure.getResidue("B", "REA");
 		structure.changeResidue(res, "RXA", {});
 
 		structure.cleanupEmptyCategories();
 
-		f.save(std::cout);
+		f.file().save(std::cout);
 	}
 	catch (const std::exception& e)
 	{

test/structure-test.cpp

@@ -77,9 +77,12 @@ BOOST_AUTO_TEST_CASE(create_nonpoly_1)
 {
     cif::VERBOSE = 1;
 
+	// do this now, avoids the need for installing
+	cif::addFileResource("mmcif_pdbx_v50.dic", "../rsrc/mmcif_pdbx_v50.dic");
+
 	mmcif::File file;
-	file.loadDictionary("mmcif_pdbx_v50.dic");
-	file.emplace("TEST");	// create a datablock
+	file.file().loadDictionary("mmcif_pdbx_v50.dic");
+	file.createDatablock("TEST");	// create a datablock
 	
 	mmcif::Structure structure(file);
 
@@ -118,17 +121,6 @@ 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
@@ -152,10 +144,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 ? 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
+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
 #
 _chem_comp.id               HEM
 _chem_comp.type             NON-POLYMER
@@ -182,178 +174,11 @@ _struct_asym.details                       ?
 
 	expected.loadDictionary("mmcif_pdbx_v50.dic");
 
-	if (not (expected.firstDatablock() == structure.datablock()))
+	if (not (expected.firstDatablock() == structure.getFile().data()))
 	{
 		BOOST_TEST(false);
 		std::cout << expected.firstDatablock() << std::endl
 				<< std::endl
-				<< structure.datablock() << std::endl;
+				<< structure.getFile().data() << std::endl;
 	}
 }
-
-// // --------------------------------------------------------------------
-
-// BOOST_AUTO_TEST_CASE(test_load_1)
-// {
-// 	mmcif::File cf(gTestDir / "5v3g.cif.gz");
-// 	mmcif::Structure s(cf);
-
-// 	for (auto &poly : s.polymers())
-// 	{
-// 		std::cout << std::string(80, '=') << std::endl;
-// 		for (auto &res : poly)
-// 		{
-// 			std::cout << res << std::endl;
-
-// 			for (auto &atom : res.atoms())
-// 				std::cout << "  " << atom << std::endl;
-// 		}
-// 	}
-// }
-
-
-// --------------------------------------------------------------------
-
-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

@@ -1,199 +0,0 @@
-/*-
- * 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");
-
-	BOOST_CHECK_NO_THROW(mmcif::Structure s2(file));
-}
-
-// --------------------------------------------------------------------
-
-BOOST_AUTO_TEST_CASE(delete_sugar_1)
-{
-	using namespace cif::literals;
-
-	const std::filesystem::path example(gTestDir / "1juh.cif.gz");
-	mmcif::File file(example.string());
-	mmcif::Structure s(file);
-
-	// Get branch for H
-	auto &bG = s.getBranchByAsymID("G");
-
-	BOOST_CHECK_EQUAL(bG.size(), 4);
-
-	s.removeResidue(bG[1]);
-
-	BOOST_CHECK_EQUAL(bG.size(), 1);
-
-	auto &bN = s.getBranchByAsymID("G");
-
-	BOOST_CHECK_EQUAL(bN.name(), "2-acetamido-2-deoxy-beta-D-glucopyranose");
-	BOOST_CHECK_EQUAL(bN.size(), 1);
-
-	file.save(gTestDir / "test-create_sugar_3.cif");
-
-	BOOST_CHECK_NO_THROW(mmcif::Structure s2(file));
-}