diff --git a/include/cif++/compound.hpp b/include/cif++/compound.hpp index 16de77f..bd39710 100644 --- a/include/cif++/compound.hpp +++ b/include/cif++/compound.hpp @@ -179,7 +179,7 @@ class compound friend class compound_factory_impl; friend class local_compound_factory_impl; - compound(cif::datablock &db); + compound(datablock &db); std::string m_id; std::string m_name; @@ -270,11 +270,15 @@ class compound_factory return is_std_base(res_name) or is_std_peptide(res_name); } + /// Return whether @a res_name is water bool is_water(std::string_view res_name) const { return res_name == "HOH" or res_name == "H2O" or res_name == "WAT"; } + /// Return whether @a res_name already exists, without creating it. + bool exists(std::string_view res_name) const; + /// \brief Create the compound object for \a id /// /// This will create the compound instance for \a id if it doesn't exist already. @@ -331,14 +335,14 @@ class compound_factory class compound_source { public: - compound_source(const cif::file &file) + compound_source(const file &file) { - cif::compound_factory::instance().push_dictionary(file); + compound_factory::instance().push_dictionary(file); } ~compound_source() { - cif::compound_factory::instance().pop_dictionary(); + compound_factory::instance().pop_dictionary(); } }; diff --git a/include/cif++/datablock.hpp b/include/cif++/datablock.hpp index bdaae69..c29d1d1 100644 --- a/include/cif++/datablock.hpp +++ b/include/cif++/datablock.hpp @@ -43,7 +43,7 @@ namespace cif /** * @brief A datablock is a list of category objects with some additional features - * + * */ class datablock : public std::list @@ -53,7 +53,7 @@ class datablock : public std::list /** * @brief Construct a new datablock object with name @a name - * + * * @param name The name for the new datablock */ datablock(std::string_view name) @@ -80,7 +80,7 @@ class datablock : public std::list { std::swap(a.m_name, b.m_name); std::swap(a.m_validator, b.m_validator); - std::swap(static_cast&>(a), static_cast&>(b)); + std::swap(static_cast &>(a), static_cast &>(b)); } // -------------------------------------------------------------------- @@ -92,7 +92,7 @@ class datablock : public std::list /** * @brief Set the name of this datablock to @a name - * + * * @param name The new name */ void set_name(std::string_view name) @@ -102,27 +102,27 @@ class datablock : public std::list /** * @brief Attempt to load the dictionary specified in audit_conform category - * + * */ void load_dictionary(); /** * @brief Set the validator object to @a v - * + * * @param v The new validator object, may be null */ void set_validator(const validator *v); /** * @brief Get the validator object - * + * * @return const validator* The validator or nullptr if there is none */ const validator *get_validator() const; /** * @brief Validates the content of this datablock and all its content - * + * * @return true If the content is valid * @return false If the content is not valid */ @@ -131,7 +131,7 @@ class datablock : public std::list /** * @brief Validates all contained data for valid links between parents and children * as defined in the validator - * + * * @return true If all links are valid * @return false If all links are not valid */ @@ -140,7 +140,7 @@ class datablock : public std::list /** * @brief Strip removes all categories and items that are invalid according * to the assigned validator. Will also add a valid audit_conform block. - * + * * @return true if the remaining datablock is valid */ bool strip(); @@ -150,7 +150,7 @@ class datablock : public std::list /** * @brief Return the category named @a name, will create a new and empty * category named @a name if it does not exist. - * + * * @param name The name of the category to return * @return category& Reference to the named category */ @@ -159,7 +159,7 @@ class datablock : public std::list /** * @brief Return the const category named @a name, will return a reference * to a static empty category if it was not found. - * + * * @param name The name of the category to return * @return category& Reference to the named category */ @@ -168,7 +168,7 @@ class datablock : public std::list /** * @brief Return a pointer to the category named @a name or nullptr if * it does not exist. - * + * * @param name The name of the category * @return category* Pointer to the category found or nullptr */ @@ -177,13 +177,12 @@ class datablock : public std::list /** * @brief Return a pointer to the category named @a name or nullptr if * it does not exist. - * + * * @param name The name of the category * @return category* Pointer to the category found or nullptr */ const category *get(std::string_view name) const; - /** * @brief Return true if this datablock contains a non-empty category */ @@ -197,7 +196,7 @@ class datablock : public std::list * new one if it is not found. The result is a tuple of an iterator * pointing to the category and a boolean indicating whether the category * was created or not. - * + * * @param name The name for the category * @return std::tuple A tuple containing an iterator pointing * at the category and a boolean indicating whether the category was newly diff --git a/src/compound.cpp b/src/compound.cpp index 7fa0324..2ca7b92 100644 --- a/src/compound.cpp +++ b/src/compound.cpp @@ -24,14 +24,38 @@ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ -#include "cif++.hpp" +#include "cif++/compound.hpp" // for compound_atom, compound_bond, compoun... -#include -#include -#include -#include -#include -#include +#include "cif++/atom_type.hpp" // for atom_type_traits +#include "cif++/category.hpp" // for category +#include "cif++/datablock.hpp" // for datablock +#include "cif++/file.hpp" // for file +#include "cif++/item.hpp" // for item +#include "cif++/iterator.hpp" // for iterator_proxy +#include "cif++/parser.hpp" // for parser +#include "cif++/point.hpp" // for distance, point +#include "cif++/row.hpp" // for tie, row_initializer, tie_wrap +#include "cif++/text.hpp" // for iequals, replace_all, iset +#include "cif++/utilities.hpp" // for load_resource, VERBOSE, colour_type + +#include // for find_if +#include // for size_t +#include // for exception, throw_with_nested +#include // for path, exists +#include // for char_traits, basic_ostream, operator<< +#include // for operator<<, quoted +#include // for clog, cout, cerr +#include // for numeric_limits +#include // for _List_iterator +#include // for allocator, map, _Rb_tree_iterator +#include // for shared_ptr, unique_ptr, __shared_ptr_... +#include // for optional +#include // for shared_lock, shared_timed_mutex +#include // for runtime_error, invalid_argument, out_... +#include // for basic_string, string, operator==, ope... +#include // for string_view, basic_string_view +#include // for pair, exchange, move +#include // for vector namespace fs = std::filesystem; @@ -140,7 +164,7 @@ compound::compound(cif::datablock &db) cif::tie(m_id, m_name, m_type, m_formula, m_formula_weight, m_formal_charge, one_letter_code, m_parent_id) = chemComp.front().get("id", "name", "type", "formula", "formula_weight", "pdbx_formal_charge", "one_letter_code", "mon_nstd_parent_comp_id"); - + if (one_letter_code.length() == 1) m_one_letter_code = one_letter_code.front(); @@ -159,7 +183,7 @@ compound::compound(cif::datablock &db) if (stereo_config.empty()) atom.stereo_config = stereo_config_type::N; else - atom.stereo_config = parse_stereo_config_from_string(stereo_config); + atom.stereo_config = parse_stereo_config_from_string(stereo_config); m_atoms.push_back(std::move(atom)); } @@ -172,7 +196,7 @@ compound::compound(cif::datablock &db) if (valueOrder.empty()) bond.type = bond_type::sing; else - bond.type = parse_bond_type_from_string(valueOrder); + bond.type = parse_bond_type_from_string(valueOrder); m_bonds.push_back(std::move(bond)); } } @@ -231,12 +255,12 @@ float compound::bond_length(const std::string &atomId_1, const std::string &atom bool compound::is_peptide() const { - return iequals(m_type, "l-peptide linking") or iequals(m_type, "peptide linking"); + return iequals(m_type, "l-peptide linking") or iequals(m_type, "peptide linking"); } bool compound::is_base() const { - return iequals(m_type, "dna linking") or iequals(m_type, "rna linking"); + return iequals(m_type, "dna linking") or iequals(m_type, "rna linking"); } // -------------------------------------------------------------------- @@ -294,12 +318,31 @@ class compound_factory_impl : public std::enable_shared_from_thisid() == id; }) != m_compounds.end()) + return true; + + return m_next and m_next->exists_self(id); + } + + bool exists(std::string_view id) + { + std::shared_lock lock(mMutex); + + return exists_self(std::string{ id }); + } + compound *get(std::string id) { std::shared_lock lock(mMutex); compound *result = nullptr; - + for (auto impl = shared_from_this(); impl; impl = impl->m_next) { result = impl->create(id); @@ -363,7 +406,9 @@ compound *compound_factory_impl::create(const std::string &id) if (m_missing.contains(id)) return nullptr; - if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end()) + if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) + { return c->id() == id; }); + i != m_compounds.end()) return *i; compound *result = nullptr; @@ -454,7 +499,6 @@ class local_compound_factory_impl : public compound_factory_impl compound *create(const std::string &id) override; private: - compound *construct_compound(const datablock &db, const std::string &id, const std::string &name, const std::string &three_letter_code, const std::string &group); cif::file m_local_file; @@ -465,7 +509,9 @@ compound *local_compound_factory_impl::create(const std::string &id) if (m_missing.contains(id)) return nullptr; - if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end()) + if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) + { return c->id() == id; }); + i != m_compounds.end()) return *i; compound *result = nullptr; @@ -510,12 +556,10 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb, float formula_weight = 0; int formal_charge = 0; - std::map formula_data; + std::map formula_data; for (std::size_t ord = 1; const auto &[atom_id, type_symbol, type, charge, x, y, z, xi, yi, zi] : - rdb["chem_comp_atom"].rows, std::optional, std::optional, - std::optional, std::optional, std::optional>( + rdb["chem_comp_atom"].rows, std::optional, std::optional, std::optional, std::optional, std::optional>( "atom_id", "type_symbol", "type", "charge", "model_Cartn_x", "model_Cartn_y", "model_Cartn_z", "pdbx_model_Cartn_x_ideal", "pdbx_model_Cartn_y_ideal", "pdbx_model_Cartn_z_ideal")) @@ -525,16 +569,14 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb, formula_data[type_symbol] += 1; - db["chem_comp_atom"].emplace({ - { "comp_id", id }, + db["chem_comp_atom"].emplace({ { "comp_id", id }, { "atom_id", atom_id }, { "type_symbol", type_symbol }, { "charge", charge }, - { "model_Cartn_x", x.has_value() ? x : xi, 3 }, - { "model_Cartn_y", y.has_value() ? y : yi, 3 }, - { "model_Cartn_z", z.has_value() ? z : zi, 3 }, - { "pdbx_ordinal", ord++ } - }); + { "model_Cartn_x", x.has_value() ? x : xi, 3 }, + { "model_Cartn_y", y.has_value() ? y : yi, 3 }, + { "model_Cartn_z", z.has_value() ? z : zi, 3 }, + { "pdbx_ordinal", ord++ } }); formal_charge += charge; } @@ -551,21 +593,19 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb, else if (cif::iequals(type, "triple") or cif::iequals(type, "trip")) value_order = "TRIP"; - db["chem_comp_bond"].emplace({ - { "comp_id", id }, + db["chem_comp_bond"].emplace({ { "comp_id", id }, { "atom_id_1", atom_id_1 }, { "atom_id_2", atom_id_2 }, { "value_order", value_order }, { "pdbx_aromatic_flag", aromatic }, // TODO: fetch stereo_config info from chem_comp_chir - { "pdbx_ordinal", ord++ } - }); + { "pdbx_ordinal", ord++ } }); } db.emplace_back(rdb["pdbx_chem_comp_descriptor"]); std::string formula; - for (bool first = true; const auto &[symbol, count]: formula_data) + for (bool first = true; const auto &[symbol, count] : formula_data) { if (std::exchange(first, false)) formula += ' '; @@ -584,15 +624,13 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb, else type = "NON-POLYMER"; - db["chem_comp"].emplace({ - { "id", id }, + db["chem_comp"].emplace({ { "id", id }, { "name", name }, { "type", type }, { "formula", formula }, { "pdbx_formal_charge", formal_charge }, { "formula_weight", formula_weight }, - { "three_letter_code", three_letter_code } - }); + { "three_letter_code", three_letter_code } }); std::shared_lock lock(mMutex); @@ -698,6 +736,11 @@ void compound_factory::pop_dictionary() m_impl = m_impl->next(); } +bool compound_factory::exists(std::string_view id) const +{ + return m_impl and m_impl->exists(id); +} + const compound *compound_factory::create(std::string_view id) { auto result = m_impl ? m_impl->get(std::string{ id }) : nullptr; @@ -722,7 +765,7 @@ bool compound_factory::is_peptide(std::string_view res_name) const bool result = is_std_peptide(res_name); if (not result and m_impl) { - auto compound = const_cast(*this).create(res_name); + auto compound = const_cast(*this).create(res_name); result = compound != nullptr and compound->is_peptide(); } return result; @@ -734,7 +777,7 @@ bool compound_factory::is_base(std::string_view res_name) const bool result = is_std_base(res_name); if (not result and m_impl) { - auto compound = const_cast(*this).create(res_name); + auto compound = const_cast(*this).create(res_name); result = compound != nullptr and compound->is_base(); } return result; diff --git a/src/pdb/pdb2cif.cpp b/src/pdb/pdb2cif.cpp index 18b05e4..3a83e7c 100644 --- a/src/pdb/pdb2cif.cpp +++ b/src/pdb/pdb2cif.cpp @@ -3142,7 +3142,6 @@ void PDBFileParser::ParseRemark350() std::map values; std::vector asymIdList; std::smatch m; - cif::row_handle genR; std::vector mat, vec; diff --git a/src/point.cpp b/src/point.cpp index 323de8a..7d76013 100644 --- a/src/point.cpp +++ b/src/point.cpp @@ -25,17 +25,17 @@ */ #include "cif++/point.hpp" -#include "cif++/matrix.hpp" -#include -#include +#include "cif++/matrix.hpp" // for matrix_subtraction, matrix_cofactors + +#include // for uniform_real_distribution, normal_distri... namespace cif { // -------------------------------------------------------------------- -template +template quaternion_type normalize(quaternion_type q) { std::valarray t(4); @@ -126,10 +126,9 @@ quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4, p4 -= p3; p3 -= p3; - quaternion q; auto axis = -p2; - float dh = dihedral_angle(p1, p2, p3, p4); + return construct_from_angle_axis(angle - dh, axis); } @@ -293,9 +292,9 @@ quaternion align_points(const std::vector &pa, const std::vector & } quaternion q( - static_cast(cf(maxR, 0)), - static_cast(cf(maxR, 1)), - static_cast(cf(maxR, 2)), + static_cast(cf(maxR, 0)), + static_cast(cf(maxR, 1)), + static_cast(cf(maxR, 2)), static_cast(cf(maxR, 3))); q = normalize(q);