label and auth seq id, some improvements

include/cif++/Structure.hpp

@@ -76,6 +76,8 @@ class Atom
 
 
 		bool getAnisoU(float anisou[6]) const;
 		bool getAnisoU(float anisou[6]) const;
 
 
+		int charge() const;
+
 		void moveTo(const Point &p);
 		void moveTo(const Point &p);
 
 
 		const Compound &comp() const;
 		const Compound &comp() const;
@@ -488,6 +490,48 @@ class Polymer : public std::vector<Monomer>
 };
 };
 
 
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------
+// Sugar and Branch, to describe glycosylation sites
+
+class Branch;
+
+class Sugar : public Residue
+{
+  public:
+	Sugar(const Branch &branch, const std::string &compoundID,
+		const std::string &asymID, int authSeqID);
+
+	int num() const { return std::stoi(mAuthSeqID); }
+	std::string name() const;
+
+	// Sugar &
+
+  private:
+	const Branch &mBranch;
+};
+
+class Branch : public std::vector<Sugar>
+{
+  public:
+	Branch(Structure &structure, const std::string &asymID);
+
+	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;
+
+	Structure *mStructure;
+	std::string mAsymID;
+};
+
+// --------------------------------------------------------------------
 // file is a reference to the data stored in e.g. the cif file.
 // file is a reference to the data stored in e.g. the cif file.
 // This object is not copyable.
 // This object is not copyable.
 
 
@@ -513,9 +557,9 @@ class File : public cif::File
 	void save(const std::filesystem::path &p) override;
 	void save(const std::filesystem::path &p) override;
 
 
 	void load(std::istream &is) override;
 	void load(std::istream &is) override;
-	
-	using cif::File::save;
+
 	using cif::File::load;
 	using cif::File::load;
+	using cif::File::save;
 
 
 	cif::Datablock &data() { return front(); }
 	cif::Datablock &data() { return front(); }
 };
 };
@@ -558,9 +602,17 @@ class Structure
 	const std::list<Polymer> &polymers() const { return mPolymers; }
 	const std::list<Polymer> &polymers() const { return mPolymers; }
 	std::list<Polymer> &polymers() { return mPolymers; }
 	std::list<Polymer> &polymers() { return mPolymers; }
 
 
-	const std::vector<Residue> &nonPolymers() const { return mNonPolymers; }
-	const std::vector<Residue> &branchResidues() const { return mBranchResidues; }
+	Polymer &getPolymerByAsymID(const std::string &asymID);
+	const Polymer &getPolymerByAsymID(const std::string &asymID) const;
+
+	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(std::string id) const;
 	Atom getAtomByID(std::string id) const;
 	// Atom getAtomByLocation(Point pt, float maxDistance) const;
 	// Atom getAtomByLocation(Point pt, float maxDistance) const;
 
 
@@ -642,6 +694,9 @@ class Structure
 	/// \return				The newly create asym ID
 	/// \return				The newly create asym ID
 	std::string createNonpoly(const std::string &entity_id, std::vector<std::vector<cif::Item>> &atom_info);
 	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 \a nag_atoms
+	Branch& createBranch(std::vector<std::vector<cif::Item>> &nag_atoms);
+
 	/// \brief Remove a residue, can be monomer or nonpoly
 	/// \brief Remove a residue, can be monomer or nonpoly
 	///
 	///
 	/// \param asym_id     The asym ID
 	/// \param asym_id     The asym ID
@@ -665,7 +720,6 @@ class Structure
 	void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
 	void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
 
 
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
-	const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }
 
 
 	void cleanupEmptyCategories();
 	void cleanupEmptyCategories();
 
 
@@ -688,6 +742,8 @@ class Structure
 
 
 	std::string insertCompound(const std::string &compoundID, bool isEntity);
 	std::string insertCompound(const std::string &compoundID, bool isEntity);
 
 
+	std::string createEntityForBranch(Branch &branch);
+
 	void loadData();
 	void loadData();
 	void updateAtomIndex();
 	void updateAtomIndex();
 
 
@@ -698,7 +754,8 @@ class Structure
 	AtomView mAtoms;
 	AtomView mAtoms;
 	std::vector<size_t> mAtomIndex;
 	std::vector<size_t> mAtomIndex;
 	std::list<Polymer> mPolymers;
 	std::list<Polymer> mPolymers;
-	std::vector<Residue> mNonPolymers, mBranchResidues;
+	std::list<Branch> mBranches;
+	std::vector<Residue> mNonPolymers;
 };
 };
 
 
 } // namespace mmcif
 } // namespace mmcif

src/Cif++.cpp

@@ -613,20 +613,12 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
 		int d = nA.compare(nB);
 		int d = nA.compare(nB);
 		if (d > 0)
 		if (d > 0)
 		{
 		{
-			auto cat = dbB.get(*catB_i);
-
-			if (cat == nullptr)
-				missingA.push_back(*catB_i);
-
+			missingA.push_back(*catB_i);
 			++catB_i;
 			++catB_i;
 		}
 		}
 		else if (d < 0)
 		else if (d < 0)
 		{
 		{
-			auto cat = dbA.get(*catA_i);
-
-			if (cat == nullptr)
-				missingB.push_back(*catA_i);
-
+			missingB.push_back(*catA_i);
 			++catA_i;
 			++catA_i;
 		}
 		}
 		else
 		else

src/Structure.cpp

@@ -136,6 +136,20 @@ bool Atom::AtomImpl::getAnisoU(float anisou[6]) const
 	return result;
 	return result;
 }
 }
 
 
+int Atom::AtomImpl::charge() const
+{
+	auto formalCharge = mRow["pdbx_formal_charge"].as<std::optional<int>>();
+
+	if (not formalCharge.has_value() and AtomTypeTraits(mType).isMetal())
+	{
+		auto &compound = comp();
+		
+		formalCharge = compound.formalCharge();
+	}
+
+	return formalCharge.value_or(0);
+}
+
 void Atom::AtomImpl::moveTo(const Point &p)
 void Atom::AtomImpl::moveTo(const Point &p)
 {
 {
 	assert(not mSymmetryCopy);
 	assert(not mSymmetryCopy);
@@ -247,7 +261,7 @@ std::string Atom::pdbID() const
 
 
 int Atom::charge() const
 int Atom::charge() const
 {
 {
-	return get_property<int>("pdbx_formal_charge");
+	return impl().charge();
 }
 }
 
 
 float Atom::occupancy() const
 float Atom::occupancy() const
@@ -1109,6 +1123,118 @@ 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)
+{
+}
+
+// bool Sugar::hasLinkedSugarAtLeavingO(int leavingO) const
+// {
+// 	return false;
+// }
+
+// Sugar &Sugar::operator[](int leavingO)
+// {
+// 	throw std::logic_error("not implemented");
+// }
+
+// const Sugar &Sugar::operator[](int leavingO) const
+// {
+// 	throw std::logic_error("not implemented");
+// }
+
+std::string Sugar::name() const
+{
+	std::string result;
+
+	if (mCompoundID == "MAN")
+		result += "alpha-D-mannopyranose";
+	else if (mCompoundID == "BMA")
+		result += "beta-D-mannopyranose";
+	else if (mCompoundID == "NAG")
+		result += "2-acetamido-2-deoxy-beta-D-glucopyranose";
+	else if (mCompoundID == "NDG")
+		result += "2-acetamido-2-deoxy-alpha-D-glucopyranose";
+	else if (mCompoundID == "FUC")
+		result += "alpha-L-fucopyranose";
+	else if (mCompoundID == "FUL")
+		result += "beta-L-fucopyranose";
+	else
+	{
+		auto compound = CompoundFactory::instance().create(mCompoundID);
+		if (compound)
+			result += compound->name();
+		else
+			result += mCompoundID;
+	}
+
+	return result;
+}
+
+Branch::Branch(Structure &structure, const std::string &asymID)
+	: mStructure(&structure)
+	, mAsymID(asymID)
+{
+	using namespace cif::literals;
+
+	auto &db = structure.datablock();
+
+	auto &struct_asym = db["struct_asym"];
+	for (const auto &[entity_id] : struct_asym.find<std::string>("id"_key == asymID, "entity_id"))
+	{
+		auto &pdbx_entity_branch_list = db["pdbx_entity_branch_list"];
+
+		for (const auto&[comp_id, num] : pdbx_entity_branch_list.find<std::string,int>(
+			"entity_id"_key == entity_id, "comp_id", "num"
+		))
+		{
+			emplace_back(*this, comp_id, asymID, num);
+		}
+
+		break;
+	}
+}
+
+std::string Branch::name() const
+{
+	// TODO: calculate
+
+	return front().name();
+
+	// std::string result;
+
+	// for (auto i = begin(); i != end(); ++i)
+	// {
+	// 	if (i->next != sugar->c1)
+	// 		continue;
+
+	// 	auto n = entityName(i) + "-(1-" + std::to_string(i->leaving_o) + ")";
+
+	// 	if (result.empty())
+	// 		result = n;
+	// 	else
+	// 		result += "-[" + n + ']';
+	// }
+
+	// if (not result.empty() and result.back() != ']')
+	// 	result += '-';
+
+
+	// return result;
+
+
+}
+
+float Branch::weight() const
+{
+	// TODO: calculate
+	return 0;
+}
+
+// --------------------------------------------------------------------
 // File
 // File
 
 
 void File::load(const std::filesystem::path &path)
 void File::load(const std::filesystem::path &path)
@@ -1304,6 +1430,14 @@ void Structure::loadData()
 			mPolymers.emplace_back(*this, entityID, asymID);
 			mPolymers.emplace_back(*this, entityID, asymID);
 	}
 	}
 
 
+	auto &branchScheme = category("pdbx_branch_scheme");
+
+	for (const auto &[asymID] : branchScheme.rows<std::string>("asym_id"))
+	{
+		if (mBranches.empty() or mBranches.back().asymID() != asymID)
+			mBranches.emplace_back(*this, asymID);
+	}
+
 	auto &nonPolyScheme = category("pdbx_nonpoly_scheme");
 	auto &nonPolyScheme = category("pdbx_nonpoly_scheme");
 
 
 	for (auto &r : nonPolyScheme)
 	for (auto &r : nonPolyScheme)
@@ -1318,17 +1452,6 @@ void Structure::loadData()
 			mNonPolymers.emplace_back(*this, monID, asymID);
 			mNonPolymers.emplace_back(*this, monID, asymID);
 	}
 	}
 
 
-	auto &branchScheme = category("pdbx_branch_scheme");
-
-	for (auto &r : branchScheme)
-	{
-		std::string asymID, monID, num;
-		cif::tie(asymID, monID, num) =
-			r.get("asym_id", "mon_id", "num");
-
-		mBranchResidues.emplace_back(*this, monID, asymID, 0, num);
-	}
-
 	// place atoms in residues
 	// place atoms in residues
 
 
 	using key_type = std::tuple<std::string, int>;
 	using key_type = std::tuple<std::string, int>;
@@ -1344,16 +1467,30 @@ void Structure::loadData()
 	for (auto &res : mNonPolymers)
 	for (auto &res : mNonPolymers)
 		resMap[{res.asymID(), (res.isWater() ? std::stoi(res.mAuthSeqID) : res.seqID())}] = &res;
 		resMap[{res.asymID(), (res.isWater() ? std::stoi(res.mAuthSeqID) : res.seqID())}] = &res;
 
 
-	for (auto &res : mBranchResidues)
-		resMap[{res.asymID(), res.seqID()}] = &res;
+	std::set<std::string> sugars;
+	for (auto &branch : mBranches)
+	{
+		for (auto &sugar : branch)
+		{
+			resMap[{sugar.asymID(), sugar.num()}] = &sugar;
+			sugars.insert(sugar.compoundID());
+		}
+	}
 
 
 	for (auto &atom : mAtoms)
 	for (auto &atom : mAtoms)
 	{
 	{
 		key_type k(atom.labelAsymID(), atom.isWater() ? std::stoi(atom.authSeqID()) : atom.labelSeqID());
 		key_type k(atom.labelAsymID(), atom.isWater() ? std::stoi(atom.authSeqID()) : atom.labelSeqID());
 		auto ri = resMap.find(k);
 		auto ri = resMap.find(k);
 
 
+		if (ri == resMap.end() and sugars.count(atom.labelCompID()))
+		{
+			k = { atom.labelAsymID(), std::stoi(atom.authSeqID()) };
+			ri = resMap.find(k);
+		}
+
 		if (ri == resMap.end())
 		if (ri == resMap.end())
 		{
 		{
+
 			if (cif::VERBOSE > 0)
 			if (cif::VERBOSE > 0)
 				std::cerr << "Missing residue for atom " << atom << std::endl;
 				std::cerr << "Missing residue for atom " << atom << std::endl;
 
 
@@ -1529,6 +1666,18 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin
 		}
 		}
 	}
 	}
 
 
+	for (auto &branch : mBranches)
+	{
+		if (branch.asymID() != asymID)
+			continue;
+
+		for (auto &sugar : branch)
+		{
+			if (sugar.asymID() == asymID and sugar.compoundID() == compID and sugar.num() == seqID)
+				return sugar;
+		}
+	}
+
 	if (seqID == 0)
 	if (seqID == 0)
 	{
 	{
 		for (auto &res : mNonPolymers)
 		for (auto &res : mNonPolymers)
@@ -1540,14 +1689,6 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin
 		}
 		}
 	}
 	}
 
 
-	for (auto &res : mBranchResidues)
-	{
-		if (res.asymID() != asymID or res.compoundID() != compID or res.seqID() != seqID)
-			continue;
-
-		return res;
-	}
-
 	throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
 	throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
 }
 }
 
 
@@ -1594,12 +1735,16 @@ const Residue &Structure::getResidue(const std::string &asymID, int seqID) const
 		}
 		}
 	}
 	}
 
 
-	for (auto &res : mBranchResidues)
+	for (auto &branch : mBranches)
 	{
 	{
-		if (res.asymID() != asymID or res.seqID() != seqID)
+		if (branch.asymID() != asymID)
 			continue;
 			continue;
 
 
-		return res;
+		for (auto &sugar : branch)
+		{
+			if (sugar.asymID() == asymID and sugar.num() == seqID)
+				return sugar;
+		}
 	}
 	}
 
 
 	throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
 	throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
@@ -1612,12 +1757,77 @@ Residue &Structure::getResidue(const std::string &asymID)
 
 
 Residue &Structure::getResidue(const mmcif::Atom &atom)
 Residue &Structure::getResidue(const mmcif::Atom &atom)
 {
 {
-	return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID());
+	using namespace cif::literals;
+
+	auto asymID = atom.labelAsymID();
+	auto entityID = atom.labelEntityID();
+	auto type = mDb["entity"].find1<std::string>("id"_key == entityID, "type");
+
+	if (type == "water")
+	{
+		auto authSeqID = atom.authSeqID();
+
+		for (auto &res : mNonPolymers)
+		{
+			if (res.asymID() != asymID or res.authSeqID() != authSeqID)
+				continue;
+
+			return res;
+		}
+
+		throw std::out_of_range("Could not find water " + asymID + '/' + authSeqID);
+	}
+	else if (type == "branched")
+	{
+		auto authSeqID = std::stoi(atom.authSeqID());
+
+		for (auto &branch : mBranches)
+		{
+			if (branch.asymID() != asymID)
+				continue;
+
+			for (auto &sugar : branch)
+			{
+				if (sugar.asymID() == asymID and sugar.num() == authSeqID)
+					return sugar;
+			}
+		}
+
+		throw std::out_of_range("Could not find sugar residue " + asymID + '/' + std::to_string(authSeqID));
+	}
+	else if (type == "polymer")
+	{
+		auto seqID = atom.labelSeqID();
+
+		for (auto &poly : mPolymers)
+		{
+			if (poly.asymID() != asymID)
+				continue;
+
+			for (auto &res : poly)
+			{
+				if (res.seqID() == seqID)
+					return res;
+			}
+		}
+
+		throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
+	}
+
+	for (auto &res : mNonPolymers)
+	{
+		if (res.asymID() != asymID)
+			continue;
+
+		return res;
+	}
+
+	throw std::out_of_range("Could not find residue " + asymID);
 }
 }
 
 
 const Residue &Structure::getResidue(const mmcif::Atom &atom) const
 const Residue &Structure::getResidue(const mmcif::Atom &atom) const
 {
 {
-	return getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID());
+	return const_cast<Structure*>(this)->getResidue(atom.labelAsymID(), atom.labelCompID(), atom.labelSeqID());
 }
 }
 
 
 std::tuple<char, int, char> Structure::MapLabelToAuth(
 std::tuple<char, int, char> Structure::MapLabelToAuth(
@@ -1673,6 +1883,25 @@ std::tuple<char, int, char> Structure::MapLabelToAuth(
 		}
 		}
 	}
 	}
 
 
+	if (not found)
+	{
+		auto r = db["pdbx_"].find(cif::Key("asym_id") == asymID);
+		if (r.size() == 1)
+		{
+			std::string pdb_strand_id, pdb_ins_code;
+			int pdb_seq_num;
+
+			cif::tie(pdb_strand_id, pdb_seq_num, pdb_ins_code) =
+				r.front().get("pdb_strand_id", "pdb_seq_num", "pdb_ins_code");
+
+			result = std::make_tuple(pdb_strand_id.front(), pdb_seq_num,
+				pdb_ins_code.empty() ? ' ' : pdb_ins_code.front());
+
+			found = true;
+		}
+	}
+
+
 	return result;
 	return result;
 }
 }
 
 
@@ -2019,7 +2248,8 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
 	{
 	{
 		auto atom_id = atom_site.getUniqueID("");
 		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},
 			{"id", atom_id},
 			{"type_symbol", atom.get_property<std::string>("type_symbol")},
 			{"type_symbol", atom.get_property<std::string>("type_symbol")},
 			{"label_atom_id", atom.get_property<std::string>("label_atom_id")},
 			{"label_atom_id", atom.get_property<std::string>("label_atom_id")},
@@ -2035,11 +2265,12 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
 			{"occupancy", atom.get_property<std::string>("occupancy")},
 			{"occupancy", atom.get_property<std::string>("occupancy")},
 			{"B_iso_or_equiv", atom.get_property<std::string>("B_iso_or_equiv")},
 			{"B_iso_or_equiv", atom.get_property<std::string>("B_iso_or_equiv")},
 			{"pdbx_formal_charge", atom.get_property<std::string>("pdbx_formal_charge")},
 			{"pdbx_formal_charge", atom.get_property<std::string>("pdbx_formal_charge")},
-			{"auth_seq_id", ""},
+			{"auth_seq_id", 1},
 			{"auth_comp_id", comp_id},
 			{"auth_comp_id", comp_id},
 			{"auth_asym_id", asym_id},
 			{"auth_asym_id", asym_id},
 			{"auth_atom_id", atom.get_property<std::string>("label_atom_id")},
 			{"auth_atom_id", atom.get_property<std::string>("label_atom_id")},
-			{"pdbx_PDB_model_num", 1}});
+			{"pdbx_PDB_model_num", 1}
+		});
 
 
 		auto &newAtom = mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
 		auto &newAtom = mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
 		res.addAtom(newAtom);
 		res.addAtom(newAtom);
@@ -2071,7 +2302,8 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 	auto &struct_asym = db["struct_asym"];
 	auto &struct_asym = db["struct_asym"];
 	std::string asym_id = struct_asym.getUniqueID();
 	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_blank_PDB_chainid_flag", "N"},
 		{"pdbx_modified", "N"},
 		{"pdbx_modified", "N"},
 		{"entity_id", entity_id},
 		{"entity_id", entity_id},
@@ -2101,7 +2333,7 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 		appendUnlessSet(atom, { "label_entity_id", entity_id} );
 		appendUnlessSet(atom, { "label_entity_id", entity_id} );
 		appendUnlessSet(atom, { "auth_comp_id", comp_id} );
 		appendUnlessSet(atom, { "auth_comp_id", comp_id} );
 		appendUnlessSet(atom, { "auth_asym_id", asym_id} );
 		appendUnlessSet(atom, { "auth_asym_id", asym_id} );
-		appendUnlessSet(atom, { "auth_seq_id", ""} );
+		appendUnlessSet(atom, { "auth_seq_id", 1} );
 		appendUnlessSet(atom, { "pdbx_PDB_model_num", 1} );
 		appendUnlessSet(atom, { "pdbx_PDB_model_num", 1} );
 		appendUnlessSet(atom, { "label_alt_id", ""} );
 		appendUnlessSet(atom, { "label_alt_id", ""} );
 
 
@@ -2129,6 +2361,142 @@ std::string Structure::createNonpoly(const std::string &entity_id, std::vector<s
 	return asym_id;
 	return asym_id;
 }
 }
 
 
+Branch& Structure::createBranch(std::vector<std::vector<cif::Item>> &nag_atoms)
+{
+	// sanity check
+	for (auto &nag_atom : nag_atoms)
+	{
+		for (auto info : nag_atom)
+		{
+			if (info.name() == "label_comp_id" and info.value() != "NAG")
+				throw std::logic_error("The first sugar in a branch should be a NAG");
+		}
+	}
+
+	using namespace cif::literals;
+
+	cif::Datablock &db = datablock();
+	auto &struct_asym = db["struct_asym"];
+	std::string asym_id = struct_asym.getUniqueID();
+
+	auto &branch = mBranches.emplace_back(*this, asym_id);
+	auto &sugar = branch.emplace_back(branch, "NAG", asym_id, 1);
+
+	auto entity_id = createEntityForBranch(branch);
+
+	struct_asym.emplace({
+		{"id", asym_id},
+		{"pdbx_blank_PDB_chainid_flag", "N"},
+		{"pdbx_modified", "N"},
+		{"entity_id", entity_id},
+		{"details", "?"}});
+
+	auto &atom_site = db["atom_site"];
+
+	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())
+			ai.emplace_back(std::move(i));
+	};
+
+	for (auto &atom : 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", 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());
+
+		auto &newAtom = mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
+		sugar.addAtom(newAtom);
+	}
+
+	db["pdbx_branch_scheme"].emplace({
+		{ "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" },
+	});
+
+	return branch;
+
+	// 		db["pdbx_entity_branch"].emplace({
+	// 			{ "entity_id", entityID },
+	// 			{ "type", "oligosaccharide" }
+	// 		});
+
+	// 		for (std::size_t i = 0; i < st.size(); ++i)
+	// 		{
+	// 			db["pdbx_entity_branch_list"].emplace({
+	// 				{ "entity_id", entityID },
+	// 				{ "comp_id", st[i].compID },
+	// 				{ "num", i + 1 },
+	// 				{ "hetero", "n" }
+	// 			});
+	// 		}
+
+	// 		for (std::size_t i = 0; i + 1 < st.size(); ++i)
+	// 		{
+	// 			auto linkID = db["pdbx_entity_branch_link"].getUniqueID("");
+
+	// 			db["pdbx_entity_branch_link"].emplace({
+	// 				{ "link_id", linkID },
+	// 				{ "atom_id_1", "C1" }
+
+	// 			});
+	// 		}
+	// 	}
+
+
+
+	// return asym_id;
+}
+
+std::string Structure::createEntityForBranch(Branch &branch)
+{
+	using namespace cif::literals;
+
+	std::string entityName = branch.name(), entityID;
+
+	auto &entity = mDb["entity"];
+
+	try
+	{
+		entityID = entity.find1<std::string>("type"_key == "branched" and "pdbx_description"_key == entityName, "id");
+	}
+	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() }
+		});
+	}
+
+	return entityID;
+}
+
 void Structure::cleanupEmptyCategories()
 void Structure::cleanupEmptyCategories()
 {
 {
 	using namespace cif::literals;
 	using namespace cif::literals;

test/structure-test.cpp

@@ -118,6 +118,17 @@ HETATM C  CHD . ? -4.342  36.262 -3.536  1.00 8.00  ?
 
 
 	auto expected = R"(
 	auto expected = R"(
 data_TEST
 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     0 
+_pdbx_nonpoly_scheme.auth_seq_num    0 
+_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_
 loop_
 _atom_site.id
 _atom_site.id
@@ -141,10 +152,10 @@ _atom_site.auth_seq_id
 _atom_site.auth_comp_id
 _atom_site.auth_comp_id
 _atom_site.auth_atom_id
 _atom_site.auth_atom_id
 _atom_site.pdbx_PDB_model_num
 _atom_site.pdbx_PDB_model_num
-1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? ? HEM CHA 1
-2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? ? HEM CHB 1
-3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? ? HEM CHC 1
-4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? ? HEM CHD 1
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? 1 HEM CHA 1
+2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? 1 HEM CHB 1
+3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? 1 HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? 1 HEM CHD 1
 #
 #
 _chem_comp.id               HEM
 _chem_comp.id               HEM
 _chem_comp.type             NON-POLYMER
 _chem_comp.type             NON-POLYMER