Fix scattering factors error

- change PDB writing interface

.gitignore

@@ -13,3 +13,4 @@ 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.0.0 LANGUAGES CXX)
+project(cifpp VERSION 4.1.1 LANGUAGES CXX)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
 

changelog

@@ -1,3 +1,15 @@
+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

include/cif++/Cif++.hpp

@@ -254,8 +254,15 @@ class Datablock
 	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);
 
+	// /// \brief Access to the category with name \a name, will throw if it doesn't exist.
+	// const Category &operator[](std::string_view name) const;
+
+	/// \brief Access to the category with name \a name, will return an empty category if is doesn't exist.
+	const Category &operator[](std::string_view name) const;
+
 	std::tuple<iterator, bool> emplace(std::string_view name);
 
 	bool isValid();
@@ -285,6 +292,9 @@ class Datablock
 	std::string mName;
 	const Validator *mValidator;
 	Datablock *mNext;
+
+	// for returning empty categories in the const operator[]
+	mutable std::unique_ptr<Category> mNullCategory;
 };
 
 // --------------------------------------------------------------------
@@ -1863,6 +1873,9 @@ class Category
 	Row front() { return Row(mHead); }
 	Row back() { return Row(mTail); }
 
+	const Row front() const { return Row(mHead); }
+	const Row back() const { return Row(mTail); }
+
 	Row operator[](Condition &&cond);
 	const Row operator[](Condition &&cond) const;
 

include/cif++/Cif2PDB.hpp

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

include/cif++/Structure.hpp

@@ -80,7 +80,7 @@ class Atom
 
 		void moveTo(const Point &p);
 
-		const Compound &comp() const;
+		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);
@@ -186,7 +186,7 @@ class Atom
 	bool isSymmetryCopy() const { return impl().mSymmetryCopy; }
 	std::string symmetry() const { return impl().mSymmetryOperator; }
 
-	const Compound &comp() const { return impl().comp(); }
+	const Compound &compound() const;
 	bool isWater() const { return impl().mCompID == "HOH" or impl().mCompID == "H2O" or impl().mCompID == "WAT"; }
 	int charge() const;
 
@@ -527,6 +527,9 @@ class Sugar : public Residue
 	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;
 
@@ -534,9 +537,14 @@ class Sugar : public Residue
 	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;
+	const Branch *mBranch;
 	Atom mLink;
 };
 
@@ -623,10 +631,14 @@ class Structure
 
 	Structure(cif::Datablock &db, size_t modelNr = 1, StructureOpenOptions options = {});
 
+	Structure(Structure &&s) = default;
+
 	// Create a read-only clone of the current structure (for multithreaded calculations that move atoms)
 	Structure(const Structure &);
 
 	Structure &operator=(const Structure &) = delete;
+	Structure &operator=(Structure &&s) = default;
+
 	~Structure();
 
 	const AtomView &atoms() const { return mAtoms; }
@@ -710,7 +722,11 @@ class Structure
 	}
 
 	// Actions
-	void removeAtom(Atom &a);
+	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,
@@ -815,6 +831,9 @@ class Structure
 
 	Atom &emplace_atom(Atom &&atom);
 
+	void removeAtom(Atom &a, bool removeFromResidue);
+	void removeSugar(Sugar &sugar);
+
 	cif::Datablock &mDb;
 	size_t mModelNr;
 	AtomView mAtoms;

src/AtomType.cpp

@@ -1083,6 +1083,20 @@ auto AtomTypeTraits::wksf(int charge) const -> const SFData&
 			return sf.sf;
 	}
 
+	if (charge != 0)
+	{
+		// Oops, not found. Fall back to zero charge and see if we can use that
+
+		if (cif::VERBOSE > 0)
+			std::cerr << "No scattering factor found for " << name() << " with charge " << charge << " will try to fall back to zero charge..." << std::endl;
+
+		for (auto& sf: data::kWKSFData)
+		{
+			if (sf.symbol == mInfo->type and sf.charge == 0)
+				return sf.sf;
+		}
+	}
+
 	throw std::runtime_error("No scattering factor found for " + name() + std::to_string(charge));
 }
 	

src/Cif++.cpp

@@ -409,6 +409,22 @@ Category &Datablock::operator[](std::string_view name)
 	return *i;
 }
 
+const Category &Datablock::operator[](std::string_view name) const
+{
+	using namespace std::literals;
+
+	auto result = get(name);
+	if (result == nullptr)
+		// throw std::out_of_range("The category with name " + std::string(name) + " does not exist");
+	{
+		std::unique_lock lock(mLock);
+		if (not mNullCategory)
+			mNullCategory.reset(new Category(const_cast<Datablock&>(*this), "<null>", nullptr));
+		result = mNullCategory.get();
+	}
+	return *result;
+}
+
 Category *Datablock::get(std::string_view name)
 {
 	std::shared_lock lock(mLock);
@@ -1525,6 +1541,11 @@ void Category::drop(const std::string &field)
 	}
 }
 
+const Row Category::operator[](Condition &&cond) const
+{
+	return const_cast<Category*>(this)->operator[](std::forward<Condition>(cond));
+}
+
 Row Category::operator[](Condition &&cond)
 {
 	Row result;

src/Cif2PDB.cpp

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

src/CifUtils.cpp

@@ -1204,15 +1204,18 @@ std::filesystem::path gDataDir;
 
 void addDataDirectory(std::filesystem::path dataDir)
 {
-	if (VERBOSE > 0 and not fs::exists(dataDir))
-		std::cerr << "The specified data directory " << dataDir << " does not exist" << std::endl;
-	gDataDir = dataDir;
+	std::error_code ec;
+	if (fs::exists(dataDir, ec))
+		gDataDir = dataDir;
+	else if (VERBOSE > 0)
+		std::cerr << "The specified data directory " << dataDir << " cannot be used: " << ec.message() << std::endl;
 }
 
 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());
+	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());
 
 	gLocalResources[name] = dataFile;
 }
@@ -1220,42 +1223,51 @@ void addFileResource(const std::string &name, std::filesystem::path 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 (gLocalResources.count(name.string()))
 	{
-		std::unique_ptr<std::ifstream> file(new std::ifstream(gLocalResources[name.string()], std::ios::binary));
-		if (file->is_open())
-			result.reset(file.release());
+		try
+		{
+			std::unique_ptr<std::ifstream> file(new std::ifstream(gLocalResources[name.string()], std::ios::binary));
+			if (file->is_open())
+				result.reset(file.release());
+		}
+		catch (...) {}
 	}
 
-	if (not result and not fs::exists(p) and not gDataDir.empty())
-		p = gDataDir / name;
+	if (not result and (not fs::exists(p, ec) or ec) and not gDataDir.empty())
+	{
+		auto p2 = gDataDir / p;
+		if (fs::exists(p2, ec) and not ec)
+			swap(p, p2);
+	}
 
 #if defined(CACHE_DIR)
-	if (not result and not fs::exists(p))
+	if (not result and (not fs::exists(p, ec) or ec))
 	{
 		auto p2 = fs::path(CACHE_DIR) / p;
-		if (fs::exists(p2))
+		if (fs::exists(p2, ec) and not ec)
 			swap(p, p2);
 	}
 #endif
 
 #if defined(DATA_DIR)
-	if (not result and not fs::exists(p))
+	if (not result and (not fs::exists(p, ec) or ec))
 	{
 		auto p2 = fs::path(DATA_DIR) / p;
-		if (fs::exists(p2))
+		if (fs::exists(p2, ec) and not ec)
 			swap(p, p2);
 	}
 #endif
 
 #if defined(CCP4) and CCP4
-	if (not result and not fs::exists(p))
+	if (not result and (not fs::exists(p, ec) or ec))
 	{
 		const char* CCP4_DIR = getenv("CCP4");
-		if (CCP4_DIR != nullptr and fs::exists(CCP4_DIR))
+		if (CCP4_DIR != nullptr and fs::exists(CCP4_DIR, ec) and not ec)
 		{
 			auto p2 = fs::path(DATA_DIR) / p;
 			if (fs::exists(p2))
@@ -1264,11 +1276,15 @@ std::unique_ptr<std::istream> loadResource(std::filesystem::path name)
 	}
 #endif
 
-	if (not result and fs::exists(p))
+	if (not result and fs::exists(p, ec) and not ec)
 	{
-		std::unique_ptr<std::ifstream> file(new std::ifstream(p, std::ios::binary));
-		if (file->is_open())
-			result.reset(file.release());
+		try
+		{
+			std::unique_ptr<std::ifstream> file(new std::ifstream(p, std::ios::binary));
+			if (file->is_open())
+				result.reset(file.release());
+		}
+		catch (...) {}		
 	}
 
 	if (not result and gResourceData)

src/CifValidator.cpp

@@ -389,6 +389,8 @@ const Validator &ValidatorFactory::operator[](std::string_view dictionary)
 		mValidators.emplace_back(dictionary, *data);
 	else
 	{
+		std::error_code ec;
+
 		// might be a compressed dictionary on disk
 		fs::path p = dictionary;
 		if (p.extension() == ".dic")
@@ -397,12 +399,12 @@ const Validator &ValidatorFactory::operator[](std::string_view dictionary)
 			p = p.parent_path() / (p.filename().string() + ".dic.gz");
 
 #if defined(CACHE_DIR) and defined(DATA_DIR)
-		if (not fs::exists(p))
+		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))
+				if (fs::exists(p2, ec) and not ec)
 				{
 					swap(p, p2);
 					break;
@@ -411,7 +413,7 @@ const Validator &ValidatorFactory::operator[](std::string_view dictionary)
 		}
 #endif
 
-		if (fs::exists(p))
+		if (fs::exists(p, ec) and not ec)
 		{
 			std::ifstream file(p, std::ios::binary);
 			if (not file.is_open())

src/Structure.cpp

@@ -148,20 +148,10 @@ int Atom::AtomImpl::charge() const
 
 	if (not formalCharge.has_value())
 	{
-		auto &compound = comp();
+		auto c = compound();
 
-		if (compound.atoms().size() == 1)
-			formalCharge = compound.atoms().front().charge;
-		// {
-		// 	for (auto cAtom : compound.atoms())
-		// 	{
-		// 		if (cAtom.id != mAtomID)
-		// 			continue;
-				
-		// 		formalCharge = cAtom.charge;
-		// 		break;
-		// 	}
-		// }
+		if (c != nullptr and c->atoms().size() == 1)
+			formalCharge = c->atoms().front().charge;
 	}
 
 	return formalCharge.value_or(0);
@@ -189,23 +179,16 @@ void Atom::AtomImpl::moveTo(const Point &p)
 	mLocation = p;
 }
 
-const Compound &Atom::AtomImpl::comp() const
+const Compound *Atom::AtomImpl::compound() const
 {
 	if (mCompound == nullptr)
 	{
-		std::string compID;
-		cif::tie(compID) = mRow.get("label_comp_id");
+		std::string compID = get_property("label_comp_id");
 
 		mCompound = CompoundFactory::instance().create(compID);
-
-		if (cif::VERBOSE > 0 and mCompound == nullptr)
-			std::cerr << "Compound not found: '" << compID << '\'' << std::endl;
 	}
 
-	if (mCompound == nullptr)
-		throw std::runtime_error("no compound");
-
-	return *mCompound;
+	return mCompound;
 }
 
 const std::string Atom::AtomImpl::get_property(const std::string_view name) const
@@ -216,7 +199,7 @@ const std::string Atom::AtomImpl::get_property(const std::string_view name) cons
 			return ref.as<std::string>();
 	}
 
-	mCachedRefs.emplace_back(name, const_cast<cif::Row&>(mRow)[name]);
+	mCachedRefs.emplace_back(name, const_cast<cif::Row &>(mRow)[name]);
 	return std::get<1>(mCachedRefs.back()).as<std::string>();
 }
 
@@ -282,6 +265,21 @@ std::string Atom::pdbID() const
 	       get_property<std::string>("pdbx_PDB_ins_code");
 }
 
+const Compound &Atom::compound() const
+{
+	auto result = impl().compound();
+
+	if (result == nullptr)
+	{
+		if (cif::VERBOSE > 0)
+			std::cerr << "Compound not found: '" << get_property<std::string>("label_comp_id") << '\'' << std::endl;
+
+		throw std::runtime_error("no compound");
+	}
+
+	return *result;
+}
+
 int Atom::charge() const
 {
 	return impl().charge();
@@ -357,8 +355,8 @@ bool Atom::operator==(const Atom &rhs) const
 {
 	if (mImpl == rhs.mImpl)
 		return true;
-	
-	if (not (mImpl and rhs.mImpl))
+
+	if (not(mImpl and rhs.mImpl))
 		return false;
 
 	return &mImpl->mDb == &rhs.mImpl->mDb and mImpl->mID == rhs.mImpl->mID;
@@ -1130,10 +1128,28 @@ int Polymer::Distance(const Monomer &a, const Monomer &b) const
 Sugar::Sugar(const Branch &branch, const std::string &compoundID,
 	const std::string &asymID, int authSeqID)
 	: Residue(branch.structure(), compoundID, asymID, 0, std::to_string(authSeqID))
-	, mBranch(branch)
+	, mBranch(&branch)
 {
 }
 
+Sugar::Sugar(Sugar &&rhs)
+	: Residue(std::forward<Residue>(rhs))
+	, mBranch(rhs.mBranch)
+{
+
+}
+
+Sugar &Sugar::operator=(Sugar &&rhs)
+{
+	if (this != &rhs)
+	{
+		Residue::operator=(std::forward<Residue>(rhs));
+		mBranch = rhs.mBranch;
+	}
+
+	return *this;
+}
+
 // bool Sugar::hasLinkedSugarAtLeavingO(int leavingO) const
 // {
 // 	return false;
@@ -1190,16 +1206,14 @@ Branch::Branch(Structure &structure, const std::string &asymID)
 
 	for (const auto &[entity_id] : struct_asym.find<std::string>("id"_key == asymID, "entity_id"))
 	{
-		for (const auto&[comp_id, num] : branch_list.find<std::string,int>(
-			"entity_id"_key == entity_id, "comp_id", "num"
-		))
+		for (const auto &[comp_id, num] : branch_list.find<std::string, int>(
+				 "entity_id"_key == entity_id, "comp_id", "num"))
 		{
 			emplace_back(*this, comp_id, asymID, num);
 		}
 
 		for (const auto &[num1, num2, atom1, atom2] : branch_link.find<size_t, size_t, std::string, std::string>(
-			"entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"
-		))
+				 "entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"))
 		{
 			if (not cif::iequals(atom1, "c1"))
 				throw std::runtime_error("invalid pdbx_entity_branch_link");
@@ -1224,8 +1238,7 @@ void Branch::linkAtoms()
 	auto entity_id = front().entityID();
 
 	for (const auto &[num1, num2, atom1, atom2] : branch_link.find<size_t, size_t, std::string, std::string>(
-		"entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"
-	))
+			 "entity_id"_key == entity_id, "entity_branch_list_num_1", "entity_branch_list_num_2", "atom_id_1", "atom_id_2"))
 	{
 		if (not cif::iequals(atom1, "c1"))
 			throw std::runtime_error("invalid pdbx_entity_branch_link");
@@ -1242,7 +1255,6 @@ std::string Branch::name() const
 	return empty() ? "" : name(front());
 }
 
-
 std::string Branch::name(const Sugar &s) const
 {
 	using namespace cif::literals;
@@ -1261,19 +1273,18 @@ std::string Branch::name(const Sugar &s) const
 
 	if (not result.empty() and result.back() != ']')
 		result += '-';
-	
+
 	return result + s.name();
 }
 
 float Branch::weight() const
 {
 	return std::accumulate(begin(), end(), 0.f, [](float sum, const Sugar &s)
-	{
+		{
 		auto compound = mmcif::CompoundFactory::instance().create(s.compoundID());
 		if (compound)
 			sum += compound->formulaWeight();
-		return sum;
-	});
+		return sum; });
 }
 
 // --------------------------------------------------------------------
@@ -1380,7 +1391,7 @@ void File::save(const std::filesystem::path &path)
 	out.push(outFile);
 
 	if (file.extension() == ".pdb")
-		WritePDBFile(out, *this);
+		WritePDBFile(out, data());
 	else
 		cif::File::save(out);
 }
@@ -1528,7 +1539,7 @@ void Structure::loadData()
 			{
 				if (res.asymID() != atom.labelAsymID())
 					continue;
-				
+
 				res.addAtom(atom);
 				break;
 			}
@@ -1614,7 +1625,7 @@ AtomView Structure::waters() const
 Atom Structure::getAtomByID(const std::string &id) const
 {
 	assert(mAtoms.size() == mAtomIndex.size());
-	
+
 	int L = 0, R = mAtoms.size() - 1;
 	while (L <= R)
 	{
@@ -1711,7 +1722,7 @@ Polymer &Structure::getPolymerByAsymID(const std::string &asymID)
 	{
 		if (poly.asymID() != asymID)
 			continue;
-		
+
 		return poly;
 	}
 
@@ -1794,7 +1805,7 @@ Residue &Structure::getResidue(const std::string &asymID, const std::string &com
 	throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID) + '-' + authSeqID);
 }
 
-Branch& Structure::getBranchByAsymID(const std::string &asymID)
+Branch &Structure::getBranchByAsymID(const std::string &asymID)
 {
 	for (auto &branch : mBranches)
 	{
@@ -1856,7 +1867,7 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
 
 // --------------------------------------------------------------------
 
-Atom& Structure::emplace_atom(Atom &&atom)
+Atom &Structure::emplace_atom(Atom &&atom)
 {
 	int L = 0, R = mAtomIndex.size() - 1;
 	while (L <= R)
@@ -1881,7 +1892,7 @@ Atom& Structure::emplace_atom(Atom &&atom)
 	return mAtoms.emplace_back(std::move(atom));
 }
 
-void Structure::removeAtom(Atom &a)
+void Structure::removeAtom(Atom &a, bool removeFromResidue)
 {
 	using namespace cif::literals;
 
@@ -1890,15 +1901,18 @@ void Structure::removeAtom(Atom &a)
 	auto &atomSites = db["atom_site"];
 	atomSites.erase("id"_key == a.id());
 
-	try
+	if (removeFromResidue)
 	{
-		auto &res = getResidue(a);
-		res.mAtoms.erase(std::remove(res.mAtoms.begin(), res.mAtoms.end(), a), res.mAtoms.end());
-	}
-	catch (const std::exception &ex)
-	{
-		if (cif::VERBOSE > 0)
-			std::cerr << "Error removing atom from residue: " << ex.what() << std::endl;
+		try
+		{
+			auto &res = getResidue(a);
+			res.mAtoms.erase(std::remove(res.mAtoms.begin(), res.mAtoms.end(), a), res.mAtoms.end());
+		}
+		catch (const std::exception &ex)
+		{
+			if (cif::VERBOSE > 0)
+				std::cerr << "Error removing atom from residue: " << ex.what() << std::endl;
+		}
 	}
 
 	assert(mAtomIndex.size() == mAtoms.size());
@@ -1940,7 +1954,7 @@ void Structure::removeAtom(Atom &a)
 			R = i - 1;
 	}
 #ifndef NDEBUG
-assert(removed);
+	assert(removed);
 #endif
 }
 
@@ -1964,7 +1978,7 @@ void Structure::swapAtoms(Atom a1, Atom a2)
 		auto l4 = r2["auth_atom_id"];
 		l3.swap(l4);
 	}
-	catch (const std::exception& ex)
+	catch (const std::exception &ex)
 	{
 		std::throw_with_nested(std::runtime_error("Failed to swap atoms"));
 	}
@@ -2085,19 +2099,16 @@ void Structure::removeResidue(Residue &res)
 	cif::Datablock &db = datablock();
 
 	auto atoms = res.atoms();
-	for (auto atom : atoms)
-		removeAtom(atom);
 
 	switch (res.entityType())
 	{
 		case EntityType::Polymer:
 		{
-			Monomer &monomer = dynamic_cast<Monomer&>(res);
+			Monomer &monomer = dynamic_cast<Monomer &>(res);
 
 			db["pdbx_poly_seq_scheme"].erase(
 				"asym_id"_key == res.asymID() and
-				"seq_id"_key == res.seqID()
-			);
+				"seq_id"_key == res.seqID());
 
 			for (auto &poly : mPolymers)
 				poly.erase(std::remove(poly.begin(), poly.end(), monomer), poly.end());
@@ -2114,14 +2125,106 @@ void Structure::removeResidue(Residue &res)
 			db["pdbx_nonpoly_scheme"].erase("asym_id"_key == res.asymID());
 			mNonPolymers.erase(std::remove(mNonPolymers.begin(), mNonPolymers.end(), res), mNonPolymers.end());
 			break;
-		
+
 		case EntityType::Branched:
-			throw std::runtime_error("Don't remove a sugar using removeResidue...");
+		{
+			Sugar &sugar = dynamic_cast<Sugar&>(res);
+
+			removeSugar(sugar);
+
+			atoms.clear();
+			break;
+		}
 
 		case EntityType::Macrolide:
 			// TODO: Fix this?
 			throw std::runtime_error("no support for macrolides yet");
 	}
+
+	for (auto atom : atoms)
+		removeAtom(atom, false);
+}
+
+void Structure::removeSugar(Sugar &sugar)
+{
+	using namespace cif::literals;
+
+	std::string asym_id = sugar.asymID();
+	Branch &branch = getBranchByAsymID(asym_id);
+	auto si = std::find(branch.begin(), branch.end(), sugar);
+	if (si == branch.end())
+		throw std::runtime_error("Sugar not part of branch");
+	size_t six = si - branch.begin();
+
+	if (six == 0)	// first sugar, means the death of this branch
+		removeBranch(branch);
+	else
+	{
+		std::set<size_t> dix;
+		std::stack<size_t> test;
+		test.push(sugar.num());
+
+		while (not test.empty())
+		{
+			auto tix = test.top();
+			test.pop();
+
+			if (dix.count(tix))
+				continue;
+			
+			dix.insert(tix);
+
+			for (auto atom : branch[tix - 1].atoms())
+				removeAtom(atom, false);
+
+			for (auto &s : branch)
+			{
+				if (s.getLinkNr() == tix)
+					test.push(s.num());
+			}
+		}
+
+		branch.erase(remove_if(branch.begin(), branch.end(), [dix](const Sugar &s) { return dix.count(s.num()); }), branch.end());
+
+		cif::Datablock &db = datablock();
+
+		auto entity_id = createEntityForBranch(branch);
+
+		// Update the entity id of the asym
+		auto &struct_asym = db["struct_asym"];
+		auto r = struct_asym.find1("id"_key == asym_id);
+		r["entity_id"] = entity_id;
+
+		for (auto &sugar : branch)
+		{
+			for (auto atom : sugar.atoms())
+				atom.set_property("label_entity_id", entity_id);
+		}
+
+		auto &pdbx_branch_scheme = db["pdbx_branch_scheme"];
+		pdbx_branch_scheme.erase("asym_id"_key == asym_id);
+
+		for (auto &sugar : branch)
+		{
+			pdbx_branch_scheme.emplace({
+				{"asym_id", asym_id},
+				{"entity_id", entity_id},
+				{"num", sugar.num()},
+				{"mon_id", sugar.compoundID()},
+
+				{"pdb_asym_id", asym_id},
+				{"pdb_seq_num", sugar.num()},
+				{"pdb_mon_id", sugar.compoundID()},
+
+				// TODO: need fix, collect from nag_atoms?
+				{"auth_asym_id", asym_id},
+				{"auth_mon_id", sugar.compoundID()},
+				{"auth_seq_num", sugar.num()},
+
+				{"hetero", "n"}
+			});
+		}
+	}
 }
 
 void Structure::removeBranch(Branch &branch)
@@ -2172,8 +2275,7 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
 	{
 		auto atom_id = atom_site.getUniqueID("");
 
-		auto &&[row, inserted] = atom_site.emplace({
-			{"group_PDB", atom.get_property<std::string>("group_PDB")},
+		auto &&[row, inserted] = atom_site.emplace({{"group_PDB", atom.get_property<std::string>("group_PDB")},
 			{"id", atom_id},
 			{"type_symbol", atom.get_property<std::string>("type_symbol")},
 			{"label_atom_id", atom.get_property<std::string>("label_atom_id")},
@@ -2193,8 +2295,7 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
 			{"auth_comp_id", comp_id},
 			{"auth_asym_id", asym_id},
 			{"auth_atom_id", atom.get_property<std::string>("label_atom_id")},
-			{"pdbx_PDB_model_num", 1}
-		});
+			{"pdbx_PDB_model_num", 1}});
 
 		auto &newAtom = emplace_atom(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
 		res.addAtom(newAtom);
@@ -2203,16 +2304,16 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
 	auto &pdbx_nonpoly_scheme = db["pdbx_nonpoly_scheme"];
 	int ndb_nr = pdbx_nonpoly_scheme.find("asym_id"_key == asym_id and "entity_id"_key == entity_id).size() + 1;
 	pdbx_nonpoly_scheme.emplace({
-		{ "asym_id", asym_id },
-		{ "entity_id", entity_id },
-		{ "mon_id", comp_id },
-		{ "ndb_seq_num", ndb_nr },
-		{ "pdb_seq_num", res.authSeqID()  },
-		{ "auth_seq_num", res.authSeqID() },
-		{ "pdb_mon_id", comp_id },
-		{ "auth_mon_id", comp_id },
-		{ "pdb_strand_id", asym_id },
-		{ "pdb_ins_code", "." },
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"mon_id", comp_id},
+		{"ndb_seq_num", ndb_nr},
+		{"pdb_seq_num", res.authSeqID()},
+		{"auth_seq_num", res.authSeqID()},
+		{"pdb_mon_id", comp_id},
+		{"auth_mon_id", comp_id},
+		{"pdb_strand_id", asym_id},
+		{"pdb_ins_code", "."},
 	});
 
 	return asym_id;
@@ -2226,8 +2327,7 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 	auto &struct_asym = db["struct_asym"];
 	std::string asym_id = struct_asym.getUniqueID();
 
-	struct_asym.emplace({
-		{"id", asym_id},
+	struct_asym.emplace({{"id", asym_id},
 		{"pdbx_blank_PDB_chainid_flag", "N"},
 		{"pdbx_modified", "N"},
 		{"entity_id", entity_id},
@@ -2241,7 +2341,8 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 
 	auto appendUnlessSet = [](std::vector<cif::Item> &ai, cif::Item &&i)
 	{
-		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci) { return ci.name() == name; }) == ai.end())
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci)
+				{ return ci.name() == name; }) == ai.end())
 			ai.emplace_back(std::move(i));
 	};
 
@@ -2249,17 +2350,17 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 	{
 		auto atom_id = atom_site.getUniqueID("");
 
-		appendUnlessSet(atom, { "group_PDB", "HETATM"} );
-		appendUnlessSet(atom, { "id", atom_id} );
-		appendUnlessSet(atom, { "label_comp_id", comp_id} );
-		appendUnlessSet(atom, { "label_asym_id", asym_id} );
-		appendUnlessSet(atom, { "label_seq_id", ""} );
-		appendUnlessSet(atom, { "label_entity_id", entity_id} );
-		appendUnlessSet(atom, { "auth_comp_id", comp_id} );
-		appendUnlessSet(atom, { "auth_asym_id", asym_id} );
-		appendUnlessSet(atom, { "auth_seq_id", 1} );
-		appendUnlessSet(atom, { "pdbx_PDB_model_num", 1} );
-		appendUnlessSet(atom, { "label_alt_id", ""} );
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_comp_id", comp_id});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_seq_id", ""});
+		appendUnlessSet(atom, {"label_entity_id", entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", comp_id});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"auth_seq_id", 1});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
 
 		auto &&[row, inserted] = atom_site.emplace(atom.begin(), atom.end());
 
@@ -2270,22 +2371,22 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 	auto &pdbx_nonpoly_scheme = db["pdbx_nonpoly_scheme"];
 	int ndb_nr = pdbx_nonpoly_scheme.find("asym_id"_key == asym_id and "entity_id"_key == entity_id).size() + 1;
 	pdbx_nonpoly_scheme.emplace({
-		{ "asym_id", asym_id },
-		{ "entity_id", entity_id },
-		{ "mon_id", comp_id },
-		{ "ndb_seq_num", ndb_nr },
-		{ "pdb_seq_num", res.authSeqID()  },
-		{ "auth_seq_num", res.authSeqID() },
-		{ "pdb_mon_id", comp_id },
-		{ "auth_mon_id", comp_id },
-		{ "pdb_strand_id", asym_id },
-		{ "pdb_ins_code", "." },
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"mon_id", comp_id},
+		{"ndb_seq_num", ndb_nr},
+		{"pdb_seq_num", res.authSeqID()},
+		{"auth_seq_num", res.authSeqID()},
+		{"pdb_mon_id", comp_id},
+		{"auth_mon_id", comp_id},
+		{"pdb_strand_id", asym_id},
+		{"pdb_ins_code", "."},
 	});
 
 	return asym_id;
 }
 
-Branch& Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
+Branch &Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
 {
 	// sanity check
 	for (auto &nag_atom : nag_atoms)
@@ -2311,7 +2412,8 @@ Branch& Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
 
 	auto appendUnlessSet = [](std::vector<cif::Item> &ai, cif::Item &&i)
 	{
-		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci) { return ci.name() == name; }) == ai.end())
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci)
+				{ return ci.name() == name; }) == ai.end())
 			ai.emplace_back(std::move(i));
 	};
 
@@ -2319,17 +2421,17 @@ Branch& Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
 	{
 		auto atom_id = atom_site.getUniqueID("");
 
-		appendUnlessSet(atom, { "group_PDB", "HETATM"} );
-		appendUnlessSet(atom, { "id", atom_id} );
-		appendUnlessSet(atom, { "label_comp_id", "NAG"} );
-		appendUnlessSet(atom, { "label_asym_id", asym_id} );
-		appendUnlessSet(atom, { "label_seq_id", "."} );
-		appendUnlessSet(atom, { "label_entity_id", tmp_entity_id} );
-		appendUnlessSet(atom, { "auth_comp_id", "NAG"} );
-		appendUnlessSet(atom, { "auth_asym_id", asym_id} );
-		appendUnlessSet(atom, { "auth_seq_id", 1} );
-		appendUnlessSet(atom, { "pdbx_PDB_model_num", 1} );
-		appendUnlessSet(atom, { "label_alt_id", ""} );
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_comp_id", "NAG"});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_seq_id", "."});
+		appendUnlessSet(atom, {"label_entity_id", tmp_entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", "NAG"});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"auth_seq_id", 1});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
 
 		auto &&[row, inserted] = atom_site.emplace(atom.begin(), atom.end());
 
@@ -2345,33 +2447,34 @@ Branch& Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
 		{"pdbx_blank_PDB_chainid_flag", "N"},
 		{"pdbx_modified", "N"},
 		{"entity_id", entity_id},
-		{"details", "?"}});
+		{"details", "?"}
+	});
 
 	for (auto &a : sugar.atoms())
 		a.set_property("label_entity_id", entity_id);
 
 	db["pdbx_branch_scheme"].emplace({
-		{ "asym_id", asym_id },
-		{ "entity_id", entity_id },
-		{ "num", 1 },
-		{ "mon_id", "NAG" },
+		{"asym_id", asym_id},
+		{"entity_id", entity_id},
+		{"num", 1},
+		{"mon_id", "NAG"},
 
-		{ "pdb_asym_id", asym_id },
-		{ "pdb_seq_num", 1 },
-		{ "pdb_mon_id", "NAG" },
+		{"pdb_asym_id", asym_id},
+		{"pdb_seq_num", 1},
+		{"pdb_mon_id", "NAG"},
 
 		// TODO: need fix, collect from nag_atoms?
-		{ "auth_asym_id", asym_id },
-		{ "auth_mon_id", "NAG" },
-		{ "auth_seq_num", 1 },
+		{"auth_asym_id", asym_id},
+		{"auth_mon_id", "NAG"},
+		{"auth_seq_num", 1},
 
-		{ "hetero", "n" }
+		{"hetero", "n"}
 	});
 
 	return branch;
 }
 
-Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vector<cif::Item>> &atom_info,
+Branch &Structure::extendBranch(const std::string &asym_id, std::vector<std::vector<cif::Item>> &atom_info,
 	int link_sugar, const std::string &link_atom)
 {
 	// sanity check
@@ -2383,7 +2486,7 @@ Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vec
 		{
 			if (info.name() != "label_comp_id")
 				continue;
-			
+
 			if (compoundID.empty())
 				compoundID = info.value();
 			else if (info.value() != compoundID)
@@ -2402,11 +2505,13 @@ Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vec
 
 	auto appendUnlessSet = [](std::vector<cif::Item> &ai, cif::Item &&i)
 	{
-		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci) { return ci.name() == name; }) == ai.end())
+		if (find_if(ai.begin(), ai.end(), [name = i.name()](cif::Item &ci)
+				{ return ci.name() == name; }) == ai.end())
 			ai.emplace_back(std::move(i));
 	};
 
-	auto bi = std::find_if(mBranches.begin(), mBranches.end(), [asym_id](Branch &b) { return b.asymID() == asym_id; });
+	auto bi = std::find_if(mBranches.begin(), mBranches.end(), [asym_id](Branch &b)
+		{ return b.asymID() == asym_id; });
 	if (bi == mBranches.end())
 		throw std::logic_error("Create a branch first!");
 
@@ -2420,14 +2525,15 @@ Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vec
 	{
 		auto atom_id = atom_site.getUniqueID("");
 
-		appendUnlessSet(atom, { "group_PDB", "HETATM"} );
-		appendUnlessSet(atom, { "id", atom_id} );
-		appendUnlessSet(atom, { "label_comp_id", compoundID} );
-		appendUnlessSet(atom, { "label_entity_id", tmp_entity_id} );
-		appendUnlessSet(atom, { "auth_comp_id", compoundID} );
-		appendUnlessSet(atom, { "auth_asym_id", asym_id} );
-		appendUnlessSet(atom, { "pdbx_PDB_model_num", 1} );
-		appendUnlessSet(atom, { "label_alt_id", ""} );
+		appendUnlessSet(atom, {"group_PDB", "HETATM"});
+		appendUnlessSet(atom, {"id", atom_id});
+		appendUnlessSet(atom, {"label_asym_id", asym_id});
+		appendUnlessSet(atom, {"label_comp_id", compoundID});
+		appendUnlessSet(atom, {"label_entity_id", tmp_entity_id});
+		appendUnlessSet(atom, {"auth_comp_id", compoundID});
+		appendUnlessSet(atom, {"auth_asym_id", asym_id});
+		appendUnlessSet(atom, {"pdbx_PDB_model_num", 1});
+		appendUnlessSet(atom, {"label_alt_id", ""});
 
 		auto &&[row, inserted] = atom_site.emplace(atom.begin(), atom.end());
 
@@ -2439,6 +2545,11 @@ Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vec
 
 	auto entity_id = createEntityForBranch(branch);
 
+	// Update the entity id of the asym
+	auto &struct_asym = db["struct_asym"];
+	auto r = struct_asym.find1("id"_key == asym_id);
+	r["entity_id"] = entity_id;
+
 	for (auto &sugar : branch)
 	{
 		for (auto atom : sugar.atoms())
@@ -2451,24 +2562,24 @@ Branch& Structure::extendBranch(const std::string &asym_id, std::vector<std::vec
 	for (auto &sugar : branch)
 	{
 		pdbx_branch_scheme.emplace({
-			{ "asym_id", asym_id },
-			{ "entity_id", entity_id },
-			{ "num", sugar.num() },
-			{ "mon_id", sugar.compoundID() },
+			{"asym_id", asym_id},
+			{"entity_id", entity_id},
+			{"num", sugar.num()},
+			{"mon_id", sugar.compoundID()},
 
-			{ "pdb_asym_id", asym_id },
-			{ "pdb_seq_num", sugar.num() },
-			{ "pdb_mon_id", sugar.compoundID() },
+			{"pdb_asym_id", asym_id},
+			{"pdb_seq_num", sugar.num()},
+			{"pdb_mon_id", sugar.compoundID()},
 
 			// TODO: need fix, collect from nag_atoms?
-			{ "auth_asym_id", asym_id },
-			{ "auth_mon_id", sugar.compoundID() },
-			{ "auth_seq_num", sugar.num() },
+			{"auth_asym_id", asym_id},
+			{"auth_mon_id", sugar.compoundID()},
+			{"auth_seq_num", sugar.num()},
 
-			{ "hetero", "n" }
+			{"hetero", "n"}
 		});
 	}
-	
+
 	return branch;
 }
 
@@ -2484,61 +2595,54 @@ std::string Structure::createEntityForBranch(Branch &branch)
 	{
 		entityID = entity.find1<std::string>("type"_key == "branched" and "pdbx_description"_key == entityName, "id");
 	}
-	catch(const std::exception& e)
+	catch (const std::exception &e)
 	{
 		entityID = entity.getUniqueID("");
 
 		if (cif::VERBOSE)
 			std::cout << "Creating new entity " << entityID << " for branched sugar " << entityName << std::endl;
 
-		entity.emplace({
-			{ "id", entityID },
-			{ "type", "branched" },
-			{ "src_method", "man" },
-			{ "pdbx_description", entityName },
-			{ "formula_weight", branch.weight() }
-		});
-	}
+		entity.emplace({{"id", entityID},
+			{"type", "branched"},
+			{"src_method", "man"},
+			{"pdbx_description", entityName},
+			{"formula_weight", branch.weight()}});
+		auto &pdbx_entity_branch_list = mDb["pdbx_entity_branch_list"];
+		for (auto &sugar : branch)
+		{
+			pdbx_entity_branch_list.emplace({
+				{"entity_id", entityID},
+				{"comp_id", sugar.compoundID()},
+				{"num", sugar.num()},
+				{"hetero", "n"}
+			});
+		}
 
-	auto &pdbx_entity_branch_list = mDb["pdbx_entity_branch_list"];
-	pdbx_entity_branch_list.erase("entity_id"_key == entityID);
+		auto &pdbx_entity_branch_link = mDb["pdbx_entity_branch_link"];
+		for (auto &s1 : branch)
+		{
+			auto l2 = s1.getLink();
 
-	for (auto &sugar : branch)
-	{
-		pdbx_entity_branch_list.emplace({
-			{ "entity_id", entityID },
-			{ "comp_id", sugar.compoundID() },
-			{ "num", sugar.num() },
-			{ "hetero", "n" }
-		});
-	}
+			if (not l2)
+				continue;
 
-	auto &pdbx_entity_branch_link = mDb["pdbx_entity_branch_link"];
-	pdbx_entity_branch_link.erase("entity_id"_key == entityID);
+			auto &s2 = branch.at(std::stoi(l2.authSeqID()) - 1);
+			auto l1 = s2.atomByID("C1");
 
-	for (auto &s1 : branch)
-	{
-		auto l1 = s1.getLink();
-
-		if (not l1)
-			continue;
-		
-		auto &s2 = branch.at(std::stoi(l1.authSeqID()) - 1);
-		auto l2 = s2.atomByID("C1");
-
-		pdbx_entity_branch_link.emplace({
-			{ "link_id", pdbx_entity_branch_link.getUniqueID("") },
-			{ "entity_id", entityID },
-			{ "entity_branch_list_num_1", s2.authSeqID() },
-			{ "comp_id_1", s2.compoundID() },
-			{ "atom_id_1", l2.labelAtomID() },
-			{ "leaving_atom_id_1", "O1" },
-			{ "entity_branch_list_num_2", s1.authSeqID() },
-			{ "comp_id_2", s1.compoundID() },
-			{ "atom_id_2", l1.labelAtomID() },
-			{ "leaving_atom_id_2", "H" + l1.labelAtomID() },
-			{ "value_order", "sing" }
-		});
+			pdbx_entity_branch_link.emplace({
+				{"link_id", pdbx_entity_branch_link.getUniqueID("")},
+				{"entity_id", entityID},
+				{"entity_branch_list_num_1", s1.authSeqID()},
+				{"comp_id_1", s1.compoundID()},
+				{"atom_id_1", l1.labelAtomID()},
+				{"leaving_atom_id_1", "O1"},
+				{"entity_branch_list_num_2", s2.authSeqID()},
+				{"comp_id_2", s2.compoundID()},
+				{"atom_id_2", l2.labelAtomID()},
+				{"leaving_atom_id_2", "H" + l2.labelAtomID()},
+				{"value_order", "sing"}
+			});
+		}
 	}
 
 	return entityID;
@@ -2665,7 +2769,7 @@ void Structure::validateAtoms() const
 		assert(i != atoms.end());
 		atoms.erase(i);
 	};
-	
+
 	for (auto &poly : mPolymers)
 	{
 		for (auto &monomer : poly)

test/sugar-test.cpp

@@ -165,4 +165,35 @@ BOOST_AUTO_TEST_CASE(create_sugar_2)
 	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));
 }