- Fix removing atoms

- Optimize isUnquotedString

include/cif++/Structure.hpp

@@ -373,6 +373,16 @@ class Residue
 
 
 	friend Structure;
 	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:
   protected:
 	Residue() {}
 	Residue() {}
 
 
@@ -706,6 +716,11 @@ class Structure
 	/// \param seq_id      The sequence ID
 	/// \param seq_id      The sequence ID
 	void removeResidue(const std::string &asym_id, int seq_id);
 	void removeResidue(const std::string &asym_id, int seq_id);
 
 
+	/// \brief Remove residue \a res
+	///
+	/// \param res         The residue to remove
+	void removeResidue(mmcif::Residue &res);
+
 	/// \brief Translate the coordinates of all atoms in the structure by \a t
 	/// \brief Translate the coordinates of all atoms in the structure by \a t
 	void translate(Point t);
 	void translate(Point t);
 
 
@@ -733,6 +748,8 @@ class Structure
 		return mDb;
 		return mDb;
 	}
 	}
 
 
+	void validateAtoms() const;
+
   private:
   private:
 	friend Polymer;
 	friend Polymer;
 	friend Residue;
 	friend Residue;

src/CifParser.cpp

@@ -96,7 +96,7 @@ bool isUnquotedString(const char *s)
 	// but be careful it does not contain e.g. stop_
 	// but be careful it does not contain e.g. stop_
 	if (result)
 	if (result)
 	{
 	{
-		const std::regex reservedRx(R"((^(?:data|save)|.*(?:loop|stop|global))_.+)", std::regex_constants::icase);
+		static const std::regex reservedRx(R"((^(?:data|save)|.*(?:loop|stop|global))_.+)", std::regex_constants::icase);
 		result = not std::regex_match(ss, reservedRx);
 		result = not std::regex_match(ss, reservedRx);
 	}
 	}
 
 

src/Structure.cpp

@@ -1779,24 +1779,30 @@ Atom& Structure::emplace_atom(Atom &&atom)
 
 
 void Structure::removeAtom(Atom &a)
 void Structure::removeAtom(Atom &a)
 {
 {
+	using namespace cif::literals;
+
 	cif::Datablock &db = datablock();
 	cif::Datablock &db = datablock();
 
 
 	auto &atomSites = db["atom_site"];
 	auto &atomSites = db["atom_site"];
+	atomSites.erase("id"_key == a.id());
 
 
-	for (auto i = atomSites.begin(); i != atomSites.end(); ++i)
+	try
 	{
 	{
-		std::string id;
-		cif::tie(id) = i->get("id");
-
-		if (id == a.id())
-		{
-			atomSites.erase(i);
-			break;
-		}
+		auto &res = getResidue(a);
+		res.mAtoms.erase(std::remove(res.mAtoms.begin(), res.mAtoms.end(), a), res.mAtoms.end());
+	}
+	catch(const std::exception& e)
+	{
+		if (cif::VERBOSE > 0)
+			std::cerr << "Error removing atom from residue: " << e.what() << std::endl;
 	}
 	}
 
 
 	assert(mAtomIndex.size() == mAtoms.size());
 	assert(mAtomIndex.size() == mAtoms.size());
 
 
+#ifndef NDEBUG
+	bool removed = false;
+#endif
+
 	int L = 0, R = mAtomIndex.size() - 1;
 	int L = 0, R = mAtomIndex.size() - 1;
 	while (L <= R)
 	while (L <= R)
 	{
 	{
@@ -1809,7 +1815,18 @@ void Structure::removeAtom(Atom &a)
 		if (d == 0)
 		if (d == 0)
 		{
 		{
 			mAtoms.erase(mAtoms.begin() + mAtomIndex[i]);
 			mAtoms.erase(mAtoms.begin() + mAtomIndex[i]);
+
+			auto ai = mAtomIndex[i];
 			mAtomIndex.erase(mAtomIndex.begin() + i);
 			mAtomIndex.erase(mAtomIndex.begin() + i);
+
+			for (auto &j : mAtomIndex)
+			{
+				if (j > ai)
+					--j;
+			}
+#ifndef NDEBUG
+			removed = true;
+#endif
 			break;
 			break;
 		}
 		}
 
 
@@ -1818,6 +1835,9 @@ void Structure::removeAtom(Atom &a)
 		else
 		else
 			R = i - 1;
 			R = i - 1;
 	}
 	}
+#ifndef NDEBUG
+assert(removed);
+#endif
 }
 }
 
 
 void Structure::swapAtoms(Atom a1, Atom a2)
 void Structure::swapAtoms(Atom a1, Atom a2)
@@ -1959,6 +1979,27 @@ void Structure::changeResidue(Residue &res, const std::string &newCompound,
 	}
 	}
 }
 }
 
 
+void Structure::removeResidue(Residue &res)
+{
+	using namespace cif::literals;
+
+	std::string asymID = res.asymID();
+
+	auto atoms = res.atoms();
+
+	for (auto atom : atoms)
+		removeAtom(atom);
+	
+	for (auto npi = mNonPolymers.begin(); npi != mNonPolymers.end(); ++npi)
+	{
+		if (not (*npi == res))
+			continue;
+		
+		mNonPolymers.erase(npi);
+		break;
+	}
+}
+
 std::string Structure::createNonPolyEntity(const std::string &comp_id)
 std::string Structure::createNonPolyEntity(const std::string &comp_id)
 {
 {
 	return insertCompound(comp_id, true);
 	return insertCompound(comp_id, true);
@@ -2431,4 +2472,47 @@ void Structure::translateRotateAndTranslate(Point t1, Quaternion q, Point t2)
 		a.translateRotateAndTranslate(t1, q, t2);
 		a.translateRotateAndTranslate(t1, q, t2);
 }
 }
 
 
+void Structure::validateAtoms() const
+{
+	// validate order
+	assert(mAtoms.size() == mAtomIndex.size());
+	for (size_t i = 0; i + i < mAtoms.size(); ++i)
+		assert(mAtoms[mAtomIndex[i]].id().compare(mAtoms[mAtomIndex[i + 1]].id()) < 0);
+
+	std::vector<Atom> atoms = mAtoms;
+
+	auto removeAtom = [&atoms](const Atom &a)
+	{
+		auto i = std::find(atoms.begin(), atoms.end(), a);
+		assert(i != atoms.end());
+		atoms.erase(i);
+	};
+	
+	for (auto &poly : mPolymers)
+	{
+		for (auto &monomer : poly)
+		{
+			for (auto &atom : monomer.atoms())
+				removeAtom(atom);
+		}
+	}
+
+	for (auto &branch : mBranches)
+	{
+		for (auto &sugar : branch)
+		{
+			for (auto &atom : sugar.atoms())
+				removeAtom(atom);
+		}
+	}
+
+	for (auto &res : mNonPolymers)
+	{
+		for (auto &atom : res.atoms())
+			removeAtom(atom);
+	}
+
+	assert(atoms.empty());
+}
+
 } // namespace mmcif
 } // namespace mmcif

test/structure-test.cpp

@@ -344,3 +344,16 @@ BOOST_AUTO_TEST_CASE(test_load_1)
 		BOOST_CHECK_EQUAL(poly.size(), pdbx_poly_seq_scheme.find("asym_id"_key == poly.asymID()).size());
 		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());
+}