No longer fail pdb conversion when missing compound info

.github/workflows/build-documentation.yml

@@ -19,7 +19,7 @@ jobs:
   docs:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v1
+    - uses: actions/checkout@v4
 
     - name: Set reusable strings
       # Turn repeated input strings (such as the build output directory) into step outputs. These step outputs can be used throughout the workflow file.
@@ -47,7 +47,7 @@ jobs:
         ls -l ${{ steps.strings.outputs.build-output-dir }}/docs/sphinx
 
     - name: Upload artifact
-      uses: actions/upload-pages-artifact@v2
+      uses: actions/upload-pages-artifact@v3
       with:
         path: ${{ steps.strings.outputs.build-output-dir }}/docs/sphinx
 
@@ -62,4 +62,4 @@ jobs:
     steps:
       - name: Deploy to GitHub Pages
         id: deployment
-        uses: actions/deploy-pages@v2
+        uses: actions/deploy-pages@v4

CMakeLists.txt

@@ -27,7 +27,7 @@ cmake_minimum_required(VERSION 3.23)
 # set the project name
 project(
 	libcifpp
-	VERSION 8.0.0
+	VERSION 8.0.1
 	LANGUAGES CXX)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
@@ -247,6 +247,8 @@ else()
 	set(EIGEN_INCLUDE_DIR ${SOURCE_DIR})
 endif()
 
+message(STATUS "Eigen include dir is ${EIGEN_INCLUDE_DIR}")
+
 # Create a revision file, containing the current git version info
 include(VersionString)
 write_version_header(${CMAKE_CURRENT_SOURCE_DIR}/src/ LIB_NAME "LibCIFPP")

changelog

@@ -1,3 +1,10 @@
+Version 8.0.1
+- Fix cif::mm::structure::cleanup_empty_categories, removed too much
+- Add default value for B_iso_or_equiv in residue::create_new_atom
+- Reconstruct some branch records in bare pdbx files
+- Fix parsing PDB files (bug due to missing validator in dest. cat.)
+- Do not fail conversion of PDB files when compound info is missing
+
 Version 8.0.0
 - A dictionary is for a datablock and a file can have
   datablocks with differing dictionaries.

include/cif++/model.hpp

@@ -34,7 +34,7 @@
 #include <numeric>
 
 #if __cpp_lib_format
-#include <format>
+# include <format>
 #endif
 
 /** @file model.hpp
@@ -540,6 +540,9 @@ class residue
 	/// \brief Return the atom with atom_id @a atomID
 	atom get_atom_by_atom_id(const std::string &atomID) const;
 
+	/// \brief Return the atom with atom_id @a atomID and alternate_id @a altID
+	atom get_atom_by_atom_id(const std::string &atomID, const std::string &altID) const;
+
 	/// \brief Return the list of atoms having ID \a atomID
 	///
 	/// This includes all alternate atoms with this ID
@@ -622,6 +625,9 @@ class monomer : public residue
 	bool is_first_in_chain() const; ///< Return if this residue is the first residue in the chain
 	bool is_last_in_chain() const;  ///< Return if this residue is the last residue in the chain
 
+	const monomer &prev() const; // Return previous monomer in polymer
+	const monomer &next() const; // Return next monomer in polymer
+
 	// convenience
 	bool has_alpha() const; ///< Return if a alpha value can be calculated (depends on location in chain)
 	bool has_kappa() const; ///< Return if a kappa value can be calculated (depends on location in chain)
@@ -992,18 +998,18 @@ class structure
 	/**
 	 * @brief Change residue @a res to a new compound ID optionally
 	 * remapping atoms.
-	 * 
+	 *
 	 * A new chem_comp entry as well as an entity is created if needed and
 	 * if the list of @a remappedAtoms is not empty it is used to remap.
-	 * 
+	 *
 	 * The array in @a remappedAtoms contains tuples of strings, both
 	 * strings contain an atom_id. The first is the one in the current
 	 * residue and the second is the atom_id that should be used instead.
 	 * If the second string is empty, the atom is removed from the residue.
-	 * 
-	 * @param res 
-	 * @param newcompound 
-	 * @param remappedAtoms 
+	 *
+	 * @param res
+	 * @param newcompound
+	 * @param remappedAtoms
 	 */
 	void change_residue(residue &res, const std::string &newcompound,
 		const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);

src/model.cpp

@@ -376,6 +376,7 @@ atom residue::create_new_atom(atom_type inType, const std::string &inAtomID, poi
 		{ "auth_comp_id", m_compound_id },
 		{ "auth_seq_id", m_auth_seq_id },
 		{ "occupancy", 1.0f, 2 },
+		{ "B_iso_or_equiv", 20.0f },
 		{ "pdbx_PDB_model_num", m_structure->get_model_nr() },
 	});
 
@@ -450,6 +451,28 @@ atom residue::get_atom_by_atom_id(const std::string &atom_id) const
 	return result;
 }
 
+atom residue::get_atom_by_atom_id(const std::string &atomID, const std::string &altID) const
+{
+	if (altID.empty())
+		return get_atom_by_atom_id(atomID);
+
+	atom result;
+
+	for (auto &a : m_atoms)
+	{
+		if (auto a_alt_id = a.get_label_alt_id(); a.get_label_atom_id() == atomID and (a_alt_id.empty() or a_alt_id == altID))
+		{
+			result = a;
+			break;
+		}
+	}
+
+	if (not result and VERBOSE > 1)
+		std::cerr << "atom with atom_id " << atomID << " and alt_id " << altID << " not found in residue " << m_asym_id << ':' << m_seq_id << '\n';
+
+	return result;
+}
+
 // residue is a single entity if the atoms for the asym with m_asym_id is equal
 // to the number of atoms in this residue...  hope this is correct....
 bool residue::is_entity() const
@@ -562,6 +585,16 @@ bool monomer::is_last_in_chain() const
 	return m_index + 1 == m_polymer->size();
 }
 
+const monomer &monomer::prev() const
+{
+	return m_polymer->at(m_index - 1);
+}
+
+const monomer &monomer::next() const
+{
+	return m_polymer->at(m_index + 1);
+}
+
 bool monomer::has_alpha() const
 {
 	return m_index >= 1 and m_index + 2 < m_polymer->size();
@@ -578,7 +611,7 @@ float monomer::phi() const
 
 	if (m_index > 0)
 	{
-		auto &prev = m_polymer->operator[](m_index - 1);
+		auto &prev = m_polymer->at(m_index - 1);
 		if (prev.m_seq_id + 1 == m_seq_id)
 		{
 			auto a1 = prev.C();
@@ -600,7 +633,7 @@ float monomer::psi() const
 
 	if (m_index + 1 < m_polymer->size())
 	{
-		auto &next = m_polymer->operator[](m_index + 1);
+		auto &next = m_polymer->at(m_index + 1);
 		if (m_seq_id + 1 == next.m_seq_id)
 		{
 			auto a1 = N();
@@ -624,9 +657,9 @@ float monomer::alpha() const
 	{
 		if (m_index >= 1 and m_index + 2 < m_polymer->size())
 		{
-			auto &prev = m_polymer->operator[](m_index - 1);
-			auto &next = m_polymer->operator[](m_index + 1);
-			auto &nextNext = m_polymer->operator[](m_index + 2);
+			auto &prev = m_polymer->at(m_index - 1);
+			auto &next = m_polymer->at(m_index + 1);
+			auto &nextNext = m_polymer->at(m_index + 2);
 
 			result = static_cast<float>(dihedral_angle(prev.CAlpha().get_location(), CAlpha().get_location(), next.CAlpha().get_location(), nextNext.CAlpha().get_location()));
 		}
@@ -648,8 +681,8 @@ float monomer::kappa() const
 	{
 		if (m_index >= 2 and m_index + 2 < m_polymer->size())
 		{
-			auto &prevPrev = m_polymer->operator[](m_index - 2);
-			auto &nextNext = m_polymer->operator[](m_index + 2);
+			auto &prevPrev = m_polymer->at(m_index - 2);
+			auto &nextNext = m_polymer->at(m_index + 2);
 
 			if (prevPrev.m_seq_id + 4 == nextNext.m_seq_id)
 			{
@@ -677,7 +710,7 @@ float monomer::tco() const
 	{
 		if (m_index > 0)
 		{
-			auto &prev = m_polymer->operator[](m_index - 1);
+			auto &prev = m_polymer->at(m_index - 1);
 			if (prev.m_seq_id + 1 == m_seq_id)
 				result = static_cast<float>(cosinus_angle(C().get_location(), O().get_location(), prev.C().get_location(), prev.O().get_location()));
 		}
@@ -699,7 +732,7 @@ float monomer::omega() const
 	try
 	{
 		if (not is_last_in_chain())
-			result = omega(*this, m_polymer->operator[](m_index + 1));
+			result = omega(*this, m_polymer->at(m_index + 1));
 	}
 	catch (const std::exception &ex)
 	{
@@ -795,7 +828,7 @@ bool monomer::is_cis() const
 
 	if (m_index + 1 < m_polymer->size())
 	{
-		auto &next = m_polymer->operator[](m_index + 1);
+		auto &next = m_polymer->at(m_index + 1);
 
 		result = monomer::is_cis(*this, next);
 	}
@@ -2717,9 +2750,12 @@ std::string structure::create_entity_for_branch(branch &branch)
 
 void structure::cleanup_empty_categories()
 {
+
 	using namespace literals;
 
 	auto &atomSite = m_db["atom_site"];
+	auto &pdbxPolySeqScheme = m_db["pdbx_poly_seq_scheme"];
+	auto &entityPolySeq = m_db["entity_poly_seq"];
 
 	// Remove chem_comp's for which there are no atoms at all
 	auto &chem_comp = m_db["chem_comp"];
@@ -2728,8 +2764,12 @@ void structure::cleanup_empty_categories()
 	for (auto chemComp : chem_comp)
 	{
 		std::string compID = chemComp["id"].as<std::string>();
-		if (atomSite.contains("label_comp_id"_key == compID or "auth_comp_id"_key == compID))
+		if (atomSite.contains("label_comp_id"_key == compID or "auth_comp_id"_key == compID) or
+			pdbxPolySeqScheme.contains("mon_id"_key == compID or "auth_mon_id"_key == compID or "pdb_mon_id"_key == compID) or
+			entityPolySeq.contains("mon_id"_key == compID))
+		{
 			continue;
+		}
 
 		obsoleteChemComps.push_back(chemComp);
 	}
@@ -2908,6 +2948,14 @@ static int compare_numbers(std::string_view a, std::string_view b)
 	return result;
 }
 
+int compare_cif_id(const std::string &a, const std::string &b)
+{
+	int d = a.length() - b.length();
+	if (d == 0)
+		d = a.compare(b);
+	return d;
+}
+
 void structure::reorder_atoms()
 {
 	auto &atom_site = m_db["atom_site"];
@@ -2919,7 +2967,7 @@ void structure::reorder_atoms()
 			// First by model number
 			d = a.get<int>("pdbx_PDB_model_num") - b.get<int>("pdbx_PDB_model_num");
 			if (d == 0)
-				d = a.get<std::string>("label_asym_id").compare(b.get<std::string>("label_asym_id"));
+				d = compare_cif_id(a.get<std::string>("label_asym_id"), b.get<std::string>("label_asym_id"));
 			if (d == 0)
 			{
 				auto na = a.get<std::optional<int>>("label_seq_id");

src/pdb/pdb2cif.cpp

@@ -3334,15 +3334,9 @@ void PDBFileParser::ParseRemark350()
 
 						std::string type = mat == std::vector<double>{ 1, 0, 0, 0, 1, 0, 0, 0, 1 } and vec == std::vector<double>{ 0, 0, 0 } ? "identity operation" : "crystal symmetry operation";
 
-						// if (type == "identity operation")
-						// {
-
-						// }
-						// else
-						try
-						{
-							// clang-format off
-							getCategory("pdbx_struct_oper_list")->emplace({
+						auto pdbx_struct_oper_list = getCategory("pdbx_struct_oper_list");
+						if (not pdbx_struct_oper_list->contains(cif::key("id") == operID))
+							getCategory("pdbx_struct_oper_list")->emplace({ // clang-format off
 								{ "id", operID },
 								{ "type", type },
 								// { "name", "" },
@@ -3360,12 +3354,7 @@ void PDBFileParser::ParseRemark350()
 								{ "matrix[3][3]", cif::format("%12.10f", mat[8]).str() },
 								{ "vector[3]", cif::format("%12.10f", vec[2]).str() }
 							});
-							// clang-format on
-						}
-						catch (duplicate_key_error &ex)
-						{
-							// so what?
-						}
+																			// clang-format on
 
 						mat.clear();
 						vec.clear();
@@ -4300,6 +4289,8 @@ void PDBFileParser::ConstructEntities()
 			type = "polypeptide(L)";
 		else if (mightBeDNA and not mightBePolyPeptide)
 			type = "polyribonucleotide";
+		else
+			type = "other";
 
 		// clang-format off
 		getCategory("entity_poly")->emplace({
@@ -4505,7 +4496,7 @@ void PDBFileParser::ConstructEntities()
 	int modResID = 1;
 	std::set<std::string> modResSet;
 	for (auto rec = FindRecord("MODRES"); rec != nullptr and rec->is("MODRES");
-		 rec = rec->mNext)                     //	 1 -  6        Record name   "MODRES"
+		rec = rec->mNext)                      //	 1 -  6        Record name   "MODRES"
 	{                                          //	 8 - 11        IDcode        idCode      ID code of this datablock.
 		std::string resName = rec->vS(13, 15); //	13 - 15        Residue name  resName     Residue name used in this datablock.
 		char chainID = rec->vC(17);            //	17             Character     chainID     Chain identifier.
@@ -5627,7 +5618,7 @@ void PDBFileParser::ParseCoordinateTransformation()
 			igiven = vC(60) == '1';   //	60             Integer       iGiven        1 if coordinates for the  representations
 			                          //	                                           which  are approximately related by the
 			GetNextRecord();          //	                                           transformations  of the molecule are
-		}                             //	                                           contained in the datablock. Otherwise, blank.
+		} //	                                           contained in the datablock. Otherwise, blank.
 
 		// clang-format off
 		getCategory("struct_ncs_oper")->emplace({
@@ -6413,7 +6404,10 @@ file read(std::istream &is)
 		// apart from the letter 'd', the test has changed into the following:
 
 		if (std::isalpha(ch) and std::toupper(ch) != 'D')
+		{
 			read_pdb_file(is, result);
+			reconstruct_pdbx(result);
+		}
 		else
 		{
 			try
@@ -6424,11 +6418,18 @@ file read(std::istream &is)
 			{
 				std::throw_with_nested(std::runtime_error("Since the file did not start with a valid PDB HEADER line mmCIF was assumed, but that failed."));
 			}
-		}
 
-		// Since we're using the cif::pdb way of reading the file, the data may need
-		// reconstruction
-		reconstruct_pdbx(result);
+			// Try to see if we can create an mm::structure out of this data.
+			// If that fails, we need to reconstruct a PDBx file out of it.
+			try
+			{
+				cif::mm::structure s(result);
+			}
+			catch (const std::exception &e)
+			{
+				reconstruct_pdbx(result);
+			}
+		}
 	}
 
 	// Must be a PDB like file, right?

src/pdb/pdb2cif_remark_3.cpp

@@ -1478,6 +1478,8 @@ bool Remark3Parser::parse(const std::string &expMethod, PDBRecord *r, cif::datab
 
 		best.parser->fixup();
 
+		auto &validator = cif::validator_factory::instance().get("mmcif_pdbx.dic");
+
 		for (auto &cat1 : best.parser->mDb)
 		{
 			if (cat1.empty())
@@ -1496,8 +1498,15 @@ bool Remark3Parser::parse(const std::string &expMethod, PDBRecord *r, cif::datab
 					auto r1 = cat1.front();
 					auto r2 = cat2.front();
 
-					for (auto item : cat1.key_items())
-						r2[item] = r1[item].text();
+					auto cv = cat1.get_cat_validator();
+					if (cv == nullptr)
+						cv = validator.get_validator_for_category(cat1.name());
+					
+					if (cv == nullptr)
+						continue;
+
+					for (auto &iv : cv->m_item_validators)
+						r2[iv.m_item_name] = r1[iv.m_item_name].text();
 				}
 			}
 			else

src/pdb/reconstruct.cpp

@@ -104,7 +104,7 @@ void checkEntities(datablock &db)
 
 		float formula_weight = 0;
 
-		if (type.empty())	// yes, that happens
+		if (type.empty()) // yes, that happens
 		{
 			const auto comp_id = db["atom_site"].find_first<std::string>("label_entity_id"_key == entity_id, "label_comp_id");
 			auto compound = cf.create(comp_id);
@@ -125,10 +125,10 @@ void checkEntities(datablock &db)
 
 			if (type.empty())
 				throw std::runtime_error("Entity without type and cannot determine what it should be");
-			
+
 			entity["type"] = type;
 		}
-		
+
 		if (type == "polymer")
 		{
 			int n = 0;
@@ -136,10 +136,10 @@ void checkEntities(datablock &db)
 			for (std::string comp_id : db["pdbx_poly_seq_scheme"].find<std::string>("entity_id"_key == entity_id, "mon_id"))
 			{
 				auto compound = cf.create(comp_id);
-				assert(compound);
-				if (not compound)
-					throw std::runtime_error("missing information for compound " + comp_id);
-				formula_weight += compound->formula_weight();
+				if (compound)
+					formula_weight += compound->formula_weight();
+				else if (cif::VERBOSE > 0)
+					std::clog << "missing information for compound " + comp_id << '\n';
 				++n;
 			}
 
@@ -154,10 +154,10 @@ void checkEntities(datablock &db)
 			for (std::string comp_id : db["pdbx_entity_branch_list"].find<std::string>("entity_id"_key == entity_id, "comp_id"))
 			{
 				auto compound = cf.create(comp_id);
-				assert(compound);
-				if (not compound)
-					throw std::runtime_error("missing information for compound " + comp_id);
-				formula_weight += compound->formula_weight();
+				if (compound)
+					formula_weight += compound->formula_weight();
+				else if (cif::VERBOSE > 0)
+					std::clog << "missing information for compound " + comp_id << '\n';
 				++n;
 			}
 
@@ -205,11 +205,11 @@ void createEntityIDs(datablock &db)
 	std::vector<residue_key_type> waters;
 
 	for (residue_key_type k : atom_site.rows<std::optional<std::string>,
-							  std::optional<int>,
-							  std::optional<std::string>,
-							  std::optional<std::string>,
-							  std::optional<int>,
-							  std::optional<std::string>>(
+			 std::optional<int>,
+			 std::optional<std::string>,
+			 std::optional<std::string>,
+			 std::optional<int>,
+			 std::optional<std::string>>(
 			 "auth_asym_id", "auth_seq_id", "auth_comp_id",
 			 "label_asym_id", "label_seq_id", "label_comp_id"))
 	{
@@ -491,7 +491,7 @@ void checkAtomRecords(datablock &db)
 		auto chem_comp_entry = chem_comp.find_first("id"_key == comp_id);
 
 		std::optional<bool> non_std;
-		if  (cf.is_monomer(comp_id))
+		if (cf.is_monomer(comp_id))
 			non_std = cf.is_std_monomer(comp_id);
 
 		if (not chem_comp_entry)
@@ -776,6 +776,7 @@ void createEntity(datablock &db)
 void createEntityPoly(datablock &db)
 {
 	using namespace literals;
+	using namespace std::literals;
 
 	auto &cf = compound_factory::instance();
 
@@ -802,7 +803,7 @@ void createEntityPoly(datablock &db)
 			auto c = cf.create(comp_id);
 
 			std::string letter;
-			char letter_can;
+			char letter_can{};
 
 			// TODO: Perhaps we should improve this...
 			if (type != "other")
@@ -914,7 +915,7 @@ void createEntityPoly(datablock &db)
 
 		entity_poly.emplace({ //
 			{ "entity_id", entity_id },
-			{ "type", type },
+			{ "type", type.empty() ? "other"s : type },
 			{ "nstd_linkage", non_std_linkage },
 			{ "nstd_monomer", non_std_monomer },
 			{ "pdbx_seq_one_letter_code", entity_seq },
@@ -1190,6 +1191,57 @@ void createPdbxNonpolyScheme(datablock &db)
 	}
 }
 
+void createPdbxBranchScheme(datablock &db)
+{
+	using namespace literals;
+
+	createPdbxEntityNonpoly(db);
+
+	auto &entity = db["entity"];
+	auto &pdbx_branch_scheme = db["pdbx_branch_scheme"];
+	auto &pdbx_entity_branch_list = db["pdbx_entity_branch_list"];
+	auto &atom_site = db["atom_site"];
+
+	for (const auto entity_id : entity.find<std::string>("type"_key == "branched", "id"))
+	{
+		for (const auto &[comp_id, asym_id, auth_seq_id] : atom_site.find<std::string, std::string, std::optional<int>>("label_entity_id"_key == entity_id, "label_comp_id", "label_asym_id", "auth_seq_id"))
+		{
+			if (not auth_seq_id.has_value())
+				throw std::runtime_error("Missing auth_seq_id on sugar atom");
+
+			int num = *auth_seq_id;
+
+			if (not pdbx_entity_branch_list.contains("entity_id"_key == entity_id and "num"_key == num))
+			{
+				pdbx_entity_branch_list.emplace({
+					// clang-format off
+
+					{ "entity_id", entity_id },
+					{ "comp_id", comp_id },
+					{ "num", num },
+
+					// clang-format on
+				});
+			}
+
+			if (not pdbx_branch_scheme.contains("entity_id"_key == entity_id and "asym_id"_key == asym_id and "num"_key == num))
+			{
+				pdbx_branch_scheme.emplace({
+					// clang-format off
+					{ "asym_id", asym_id },
+					{ "entity_id", entity_id },
+					{ "mon_id", comp_id },
+					{ "num", num },
+					{ "pdb_asym_id", asym_id },
+					{ "pdb_mon_id", comp_id },
+					{ "pdb_seq_num", num }
+					// clang-format on
+				});
+			}
+		}
+	}
+}
+
 bool reconstruct_pdbx(file &file)
 {
 	if (file.empty())
@@ -1480,9 +1532,6 @@ bool reconstruct_pdbx(file &file, const validator &validator)
 	if (auto cat = db.get("entity"); cat == nullptr or cat->empty())
 		createEntity(db);
 
-	// fill in missing formula_weight, e.g.
-	checkEntities(db);
-
 	if (auto cat = db.get("pdbx_poly_seq_scheme"); cat == nullptr or cat->empty())
 		createPdbxPolySeqScheme(db);
 
@@ -1491,6 +1540,12 @@ bool reconstruct_pdbx(file &file, const validator &validator)
 	
 	createPdbxNonpolyScheme(db);
 
+	// Create a minimal set of branch records
+	createPdbxBranchScheme(db);
+
+	// fill in missing formula_weight, e.g.
+	checkEntities(db);
+
 	// skip unknown categories for now
 	bool valid = true;
 	for (auto &cat : db)