Skip to content

L
libcifpp
Commit f3c2e591 by Maarten L. Hekkelman
Options

Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop

parents 24ab660e 4732004b
Showing with 1001 additions and 1453 deletions
CMakeLists.txt
......@@ -25,7 +25,7 @@
cmake_minimum_required(VERSION 3.16)
# set the project name
project(cifpp VERSION 3.0.1 LANGUAGES CXX)
project(cifpp VERSION 3.0.2 LANGUAGES CXX)
list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
......@@ -133,8 +133,6 @@ if(MSVC)
message(FATAL_ERROR "Unsupported or unknown processor type ${CMAKE_SYSTEM_PROCESSOR}")
endif()
set(COFF_SPEC "--coff=${COFF_TYPE}")
# for mrc, just in case
list(APPEND CMAKE_PREFIX_PATH "$ENV{LOCALAPPDATA}/mrc")
endif()
......@@ -150,12 +148,12 @@ endif()
if(WIN32 AND BUILD_SHARED_LIBS)
message("Not using resources when building shared libraries for Windows")
else()
find_program(MRC mrc)
find_package(Mrc)
if(MRC)
option(CIFPP_USE_RSRC "Use mrc to create resources" ON)
if(MRC_FOUND)
option(USE_RSRC "Use mrc to create resources" ON)
else()
message("Using resources not possible since mrc was not found")
message(WARNING "Not using resources since mrc was not found")
endif()
if(CIFPP_USE_RSRC STREQUAL "ON")
......@@ -176,9 +174,7 @@ find_package(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_opt
if(NOT MSVC AND Boost_USE_STATIC_LIBS)
find_package(ZLIB REQUIRED)
find_package(BZip2 REQUIRED)
list(APPEND CIFPP_REQUIRED_LIBRARIES ${ZLIB_LIBRARIES} ${BZip2_LIBRARIES})
list(APPEND CIFPP_REQUIRED_LIBRARIES ${ZLIB_LIBRARIES})
endif()
include(FindFilesystem)
......@@ -341,6 +337,13 @@ install(TARGETS cifpp
RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
if(MSVC)
install(
FILES $<TARGET_PDB_FILE:${PROJECT_NAME}>
DESTINATION ${CMAKE_INSTALL_LIBDIR}
OPTIONAL)
endif()
install(EXPORT cifppTargets
FILE "cifppTargets.cmake"
NAMESPACE cifpp::
......@@ -414,13 +417,6 @@ option(CIFPP_BUILD_TESTS "Build test exectuables" OFF)
if(CIFPP_BUILD_TESTS)
if(CIFPP_USE_RSRC)
add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj
COMMAND ${MRC} -o ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj ${CMAKE_SOURCE_DIR}/rsrc/mmcif_pdbx_v50.dic ${COFF_SPEC}
)
set(CIFPP_TEST_RESOURCE ${CMAKE_CURRENT_BINARY_DIR}/cifpp_test_rsrc.obj)
endif()
list(APPEND CIFPP_tests
# pdb2cif
rename-compound
......@@ -431,7 +427,7 @@ if(CIFPP_BUILD_TESTS)
set(CIFPP_TEST "${CIFPP_TEST}-test")
set(CIFPP_TEST_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/test/${CIFPP_TEST}.cpp")
add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE} ${CIFPP_TEST_RESOURCE})
add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE})
target_include_directories(${CIFPP_TEST} PRIVATE
${CMAKE_CURRENT_SOURCE_DIR}/include
......@@ -439,7 +435,11 @@ if(CIFPP_BUILD_TESTS)
)
target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp ${Boost_LIBRARIES} ${CIFPP_REQUIRED_LIBRARIES})
if(CIFPP_USE_RSRC)
mrc_target_resources(${CIFPP_TEST} ${CMAKE_SOURCE_DIR}/rsrc/mmcif_pdbx_v50.dic)
endif()
if(MSVC)
# Specify unwind semantics so that MSVC knowns how to handle exceptions
target_compile_options(${CIFPP_TEST} PRIVATE /EHsc)
......
Config.cmake.in
......@@ -4,7 +4,6 @@ include(CMakeFindDependencyMacro)
find_dependency(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_options)
if(NOT WIN32)
find_dependency(ZLIB)
find_dependency(BZip2)
endif()
INCLUDE("${CMAKE_CURRENT_LIST_DIR}/cifppTargets.cmake")
......
changelog
Version 3.0.2
- refactored mmcif::Atom for performance reasons
Version 3.0.1
- Fixed processing of proline restraints file from CCP4, proline
is a peptide, really.
- Added code to facilitate DSSP
Version 3.0.0
- Replaced many strings in the API with string_view for
performance reasons.
......
include/cif++/Cif++.hpp
......@@ -792,7 +792,7 @@ class Row
}
void assign(const std::vector<Item> &values);
void assign(std::string_view name, const std::string &value, bool updateLinked);
void assign(std::string_view name, const std::string &value, bool updateLinked, bool validate = true);
bool operator==(const Row &rhs) const
{
......@@ -814,7 +814,7 @@ class Row
friend std::ostream &operator<<(std::ostream &os, const Row &row);
private:
void assign(size_t column, const std::string &value, bool updateLinked);
void assign(size_t column, const std::string &value, bool updateLinked, bool validate = true);
void assign(const Item &i, bool updateLinked);
static void swap(size_t column, ItemRow *a, ItemRow *b);
......@@ -2152,6 +2152,7 @@ class Category
std::vector<ItemColumn> mColumns;
ItemRow *mHead;
ItemRow *mTail;
size_t mLastUniqueNr = 0;
class CatIndex *mIndex;
std::vector<Linked> mParentLinks, mChildLinks;
......
include/cif++/Secondary.hpp
......@@ -184,6 +184,15 @@ class DSSP
std::tuple<ResidueInfo,double> acceptor(int i) const;
std::tuple<ResidueInfo,double> donor(int i) const;
/// \brief Simple compare equals
bool operator==(const ResidueInfo &rhs) const
{
return mImpl == rhs.mImpl;
}
/// \brief Returns \result true if there is a bond between two residues
friend bool TestBond(ResidueInfo const &a, ResidueInfo const &b);
private:
ResidueInfo(Res* res) : mImpl(res) {}
......@@ -193,7 +202,7 @@ class DSSP
class iterator
{
public:
using iterator_category = std::input_iterator_tag;
using iterator_category = std::bidirectional_iterator_tag;
using value_type = ResidueInfo;
using difference_type = std::ptrdiff_t;
using pointer = value_type*;
......@@ -214,6 +223,14 @@ class DSSP
return tmp;
}
iterator& operator--();
iterator operator--(int)
{
auto tmp(*this);
this->operator--();
return tmp;
}
bool operator==(const iterator& rhs) const { return mCurrent.mImpl == rhs.mCurrent.mImpl; }
bool operator!=(const iterator& rhs) const { return mCurrent.mImpl != rhs.mCurrent.mImpl; }
......
include/cif++/Structure.hpp
......@@ -60,30 +60,101 @@ class File;
class Atom
{
private:
struct AtomImpl : public std::enable_shared_from_this<AtomImpl>
{
AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row);
// constructor for a symmetry copy of an atom
AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op);
AtomImpl(const AtomImpl &i) = default;
void prefetch();
int compare(const AtomImpl &b) const;
bool getAnisoU(float anisou[6]) const;
void moveTo(const Point &p);
const Compound &comp() const;
const std::string get_property(const std::string_view name) const;
void set_property(const std::string_view name, const std::string &value);
const cif::Datablock &mDb;
std::string mID;
AtomType mType;
std::string mAtomID;
std::string mCompID;
std::string mAsymID;
int mSeqID;
std::string mAltID;
std::string mAuthSeqID;
Point mLocation;
int mRefcount;
cif::Row mRow;
mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
mutable const Compound *mCompound = nullptr;
bool mSymmetryCopy = false;
bool mClone = false;
std::string mSymmetryOperator = "1_555";
};
public:
Atom();
Atom(struct AtomImpl *impl);
Atom(const Atom &rhs);
Atom() {}
Atom(std::shared_ptr<AtomImpl> impl)
: mImpl(impl) {}
Atom(const Atom &rhs)
: mImpl(rhs.mImpl) {}
Atom(cif::Datablock &db, cif::Row &row);
// a special constructor to create symmetry copies
Atom(const Atom &rhs, const Point &symmmetry_location, const std::string &symmetry_operation);
~Atom();
explicit operator bool() const { return mImpl_ != nullptr; }
explicit operator bool() const { return (bool)mImpl; }
// return a copy of this atom, with data copied instead of referenced
Atom clone() const;
Atom clone() const
{
auto copy = std::make_shared<AtomImpl>(*mImpl);
copy->mClone = true;
return Atom(copy);
}
Atom &operator=(const Atom &rhs) = default;
template <typename T>
T get_property(const std::string_view name) const;
Atom &operator=(const Atom &rhs);
void set_property(const std::string_view name, const std::string &value)
{
mImpl->set_property(name, value);
}
const std::string &id() const;
AtomType type() const;
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
void property(const std::string_view name, const T &value)
{
set_property(name, std::to_string(value));
}
Point location() const;
void location(Point p);
const std::string &id() const { return mImpl->mID; }
AtomType type() const { return mImpl->mType; }
Point location() const { return mImpl->mLocation; }
void location(Point p) { mImpl->moveTo(p); }
/// \brief Translate the position of this atom by \a t
void translate(Point t);
......@@ -91,47 +162,40 @@ class Atom
/// \brief Rotate the position of this atom by \a q
void rotate(Quaternion q);
/// \brief Translate and rotate the position of this atom by \a t and \a q
void translateAndRotate(Point t, Quaternion q);
/// \brief Translate, rotate and translate again the coordinates this atom by \a t1 , \a q and \a t2
void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
// for direct access to underlying data, be careful!
const cif::Row getRow() const;
const cif::Row getRow() const { return mImpl->mRow; }
const cif::Row getRowAniso() const;
// Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
bool isSymmetryCopy() const;
std::string symmetry() const;
// const clipper::RTop_orth& symop() const;
bool isSymmetryCopy() const { return mImpl->mSymmetryCopy; }
std::string symmetry() const { return mImpl->mSymmetryOperator; }
const Compound &comp() const;
bool isWater() const;
const Compound &comp() const { return mImpl->comp(); }
bool isWater() const { return mImpl->mCompID == "HOH" or mImpl->mCompID == "H2O" or mImpl->mCompID == "WAT"; }
int charge() const;
float uIso() const;
bool getAnisoU(float anisou[6]) const;
bool getAnisoU(float anisou[6]) const { return mImpl->getAnisoU(anisou); }
float occupancy() const;
template <typename T>
T property(const std::string_view name) const;
void property(const std::string_view name, const std::string &value);
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
void property(const std::string_view name, const T &value)
{
property(name, std::to_string(value));
}
// specifications
std::string labelAtomID() const;
std::string labelCompID() const;
std::string labelAsymID() const;
const std::string& labelAtomID() const { return mImpl->mAtomID; }
const std::string& labelCompID() const { return mImpl->mCompID; }
const std::string& labelAsymID() const { return mImpl->mAsymID; }
std::string labelEntityID() const;
int labelSeqID() const;
std::string labelAltID() const;
bool isAlternate() const;
int labelSeqID() const { return mImpl->mSeqID; }
const std::string& labelAltID() const { return mImpl->mAltID; }
bool isAlternate() const { return not mImpl->mAltID.empty(); }
std::string authAtomID() const;
std::string authCompID() const;
std::string authAsymID() const;
std::string authSeqID() const;
const std::string& authSeqID() const { return mImpl->mAuthSeqID; }
std::string pdbxAuthInsCode() const;
std::string pdbxAuthAltID() const;
......@@ -140,13 +204,6 @@ class Atom
bool operator==(const Atom &rhs) const;
// // get clipper format Atom
// clipper::Atom toClipper() const;
// Radius calculation based on integrating the density until perc of electrons is found
void calculateRadius(float resHigh, float resLow, float perc);
float radius() const;
// access data in compound for this atom
// convenience routine
......@@ -158,10 +215,10 @@ class Atom
void swap(Atom &b)
{
std::swap(mImpl_, b.mImpl_);
std::swap(mImpl, b.mImpl);
}
int compare(const Atom &b) const;
int compare(const Atom &b) const { return mImpl->compare(*b.mImpl); }
bool operator<(const Atom &rhs) const
{
......@@ -172,14 +229,31 @@ class Atom
private:
friend class Structure;
void setID(int id);
AtomImpl *impl();
const AtomImpl *impl() const;
void setID(int id);
struct AtomImpl *mImpl_;
std::shared_ptr<AtomImpl> mImpl;
};
template <>
inline std::string Atom::get_property<std::string>(const std::string_view name) const
{
return mImpl->get_property(name);
}
template <>
inline int Atom::get_property<int>(const std::string_view name) const
{
auto v = mImpl->get_property(name);
return v.empty() ? 0 : stoi(v);
}
template <>
inline float Atom::get_property<float>(const std::string_view name) const
{
return stof(mImpl->get_property(name));
}
inline void swap(mmcif::Atom &a, mmcif::Atom &b)
{
a.swap(b);
......@@ -202,19 +276,16 @@ typedef std::vector<Atom> AtomView;
class Residue
{
public:
// constructors should be private, but that's not possible for now (needed in emplace)
// constructor for waters
Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, const std::string &authSeqID);
// constructor for a residue without a sequence number
Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID);
// constructor for a residue with a sequence number
// constructor
Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, int seqID, const std::string &authSeqID);
const std::string &asymID, int seqID = 0, const std::string &authSeqID = {})
: mStructure(&structure)
, mCompoundID(compoundID)
, mAsymID(asymID)
, mSeqID(seqID)
, mAuthSeqID(authSeqID)
{
}
Residue(const Residue &rhs) = delete;
Residue &operator=(const Residue &rhs) = delete;
......@@ -227,6 +298,11 @@ class Residue
const Compound &compound() const;
const AtomView &atoms() const;
void addAtom(const Atom &atom)
{
mAtoms.push_back(atom);
}
/// \brief Unique atoms returns only the atoms without alternates and the first of each alternate atom id.
AtomView unique_atoms() const;
......@@ -277,6 +353,8 @@ class Residue
friend std::ostream &operator<<(std::ostream &os, const Residue &res);
friend Structure;
protected:
Residue() {}
......@@ -537,6 +615,12 @@ class Structure
/// \brief Rotate the coordinates of all atoms in the structure by \a q
void rotate(Quaternion t);
/// \brief Translate and rotate the coordinates of all atoms in the structure by \a t and \a q
void translateAndRotate(Point t, Quaternion q);
/// \brief Translate, rotate and translate again the coordinates of all atoms in the structure by \a t1 , \a q and \a t2
void translateRotateAndTranslate(Point t1, Quaternion q, Point t2);
const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }
......
src/BondMap.cpp
......@@ -180,7 +180,10 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &a
auto compound = mmcif::CompoundFactory::instance().create(compoundID);
if (not compound)
std::cerr << "Missing compound bond info for " << compoundID << std::endl;
{
if (cif::VERBOSE >= 0)
std::cerr << "Missing compound bond info for " << compoundID << std::endl;
}
else
{
for (auto &atom : compound->bonds())
......@@ -308,7 +311,7 @@ BondMap::BondMap(const Structure &p)
{
if (c == "HOH" or c == "H2O" or c == "WAT")
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping water in bond map calculation" << std::endl;
continue;
}
......
src/Cif++.cpp
/*-
* 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
......@@ -30,10 +30,10 @@
#include <numeric>
#include <regex>
#include <set>
#include <shared_mutex>
#include <stack>
#include <tuple>
#include <unordered_map>
#include <shared_mutex>
#include <filesystem>
......@@ -570,36 +570,6 @@ void Datablock::write(std::ostream &os, const std::vector<std::string> &order)
cat.write(os);
}
// // mmcif support, sort of. First write the 'entry' Category
// // and if it exists, _AND_ we have a Validator, write out the
// // auditConform record.
//
// for (auto& cat: mCategories)
// {
// if (cat.name() == "entry")
// {
// cat.write(os);
//
// if (mValidator != nullptr)
// {
// Category auditConform(*this, "audit_conform", nullptr);
// auditConform.emplace({
// { "dict_name", mValidator->dictName() },
// { "dict_version", mValidator->dictVersion() }
// });
// auditConform.write(os);
// }
//
// break;
// }
// }
//
// for (auto& cat: mCategories)
// {
// if (cat.name() != "entry" and cat.name() != "audit_conform")
// cat.write(os);
// }
}
bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
......@@ -636,15 +606,15 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
{
std::string nA = *catA_i;
ba::to_lower(nA);
std::string nB = *catB_i;
ba::to_lower(nB);
int d = nA.compare(nB);
if (d > 0)
{
auto cat = dbB.get(*catB_i);
if (cat == nullptr)
missingA.push_back(*catB_i);
......@@ -653,7 +623,7 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
else if (d < 0)
{
auto cat = dbA.get(*catA_i);
if (cat == nullptr)
missingB.push_back(*catA_i);
......@@ -662,25 +632,25 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
else
++catA_i, ++catB_i;
}
while (catA_i != catA.end())
missingB.push_back(*catA_i++);
while (catB_i != catB.end())
missingA.push_back(*catB_i++);
if (not (missingA.empty() and missingB.empty()))
if (not(missingA.empty() and missingB.empty()))
{
if (cif::VERBOSE > 1)
{
std::cerr << "compare of datablocks failed" << std::endl;
if (not missingA.empty())
std::cerr << "Categories missing in A: " << ba::join(missingA, ", ") << std::endl
<< std::endl;
<< std::endl;
if (not missingB.empty())
std::cerr << "Categories missing in B: " << ba::join(missingB, ", ") << std::endl
<< std::endl;
<< std::endl;
result = false;
}
......@@ -706,7 +676,7 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
++catA_i;
else
{
if (not (*dbA.get(*catA_i) == *dbB.get(*catB_i)))
if (not(*dbA.get(*catA_i) == *dbB.get(*catB_i)))
{
if (cif::VERBOSE > 1)
{
......@@ -724,10 +694,10 @@ bool operator==(const cif::Datablock &dbA, const cif::Datablock &dbB)
return result;
}
std::ostream& operator<<(std::ostream &os, const Datablock &data)
std::ostream &operator<<(std::ostream &os, const Datablock &data)
{
// whoohoo... this sucks!
const_cast<Datablock&>(data).write(os);
const_cast<Datablock &>(data).write(os);
return os;
}
......@@ -1162,7 +1132,7 @@ void CatIndex::reconstruct()
insert(r.mData);
// maybe reconstruction can be done quicker by using the following commented code.
// however, I've not had the time to think of a way to std::set the red/black flag correctly in that case.
// however, I've not had the time to think of a way to set the red/black flag correctly in that case.
// std::vector<ItemRow*> rows;
// transform(mCat.begin(), mCat.end(), backInserter(rows),
......@@ -1254,82 +1224,15 @@ size_t CatIndex::size() const
return result;
}
//bool CatIndex::isValid() const
//{
// bool result = true;
//
// if (mRoot != nullptr)
// {
// uint32_t minBlack = numeric_limits<uint32_t>::max();
// uint32_t maxBlack = 0;
//
// assert(not mRoot->mRed);
//
// result = isValid(mRoot, false, 0, minBlack, maxBlack);
// assert(minBlack == maxBlack);
// }
//
// return result;
//}
//
//bool CatIndex::validate(entry* h, bool isParentRed, uint32_t blackDepth, uint32_t& minBlack, uint32_t& maxBlack) const
//{
// bool result = true;
//
// if (h->mRed)
// assert(not isParentRed);
// else
// ++blackDepth;
//
// if (isParentRed)
// assert(not h->mRed);
//
// if (h->mLeft != nullptr and h->mRight != nullptr)
// {
// if (isRed(h->mLeft))
// assert(not isRed(h->mRight));
// if (isRed(h->mRight))
// assert(not isRed(h->mLeft));
// }
//
// if (h->mLeft != nullptr)
// {
// assert(mComp(h->mLeft->mRow, h->mRow) < 0);
// validate(h->mLeft, h->mRed, blackDepth, minBlack, maxBlack);
// }
// else
// {
// if (minBlack > blackDepth)
// minBlack = blackDepth;
// if (maxBlack < blackDepth)
// maxBlack = blackDepth;
// }
//
// if (h->mRight != nullptr)
// {
// assert(mComp(h->mRight->mRow, h->mRow) > 0);
// validate(h->mRight, h->mRight, blackDepth, minBlack, maxBlack);
// }
// else
// {
// if (minBlack > blackDepth)
// minBlack = blackDepth;
// if (maxBlack < blackDepth)
// maxBlack = blackDepth;
// }
//}
// --------------------------------------------------------------------
RowSet::RowSet(Category &cat)
: mCat(&cat)
// , mCond(nullptr)
{
}
RowSet::RowSet(Category &cat, Condition &&cond)
: mCat(&cat)
// , mCond(new Condition(std::forward<Condition>(cond)))
{
cond.prepare(cat);
......@@ -1343,21 +1246,17 @@ RowSet::RowSet(Category &cat, Condition &&cond)
RowSet::RowSet(const RowSet &rhs)
: mCat(rhs.mCat)
, mItems(rhs.mItems)
// , mCond(nullptr)
{
}
RowSet::RowSet(RowSet &&rhs)
: mCat(rhs.mCat)
, mItems(std::move(rhs.mItems))
// , mCond(rhs.mCond)
{
// rhs.mCond = nullptr;
}
RowSet::~RowSet()
{
// delete mCond;
}
RowSet &RowSet::operator=(const RowSet &rhs)
......@@ -1469,7 +1368,7 @@ void Category::updateLinks()
auto childCat = mDb.get(link->mChildCategory);
if (childCat == nullptr)
continue;
mChildLinks.push_back({ childCat, link });
mChildLinks.push_back({childCat, link});
}
for (auto link : mValidator->getLinksForChild(mName))
......@@ -1477,7 +1376,7 @@ void Category::updateLinks()
auto parentCat = mDb.get(link->mParentCategory);
if (parentCat == nullptr)
continue;
mParentLinks.push_back({ parentCat, link });
mParentLinks.push_back({parentCat, link});
}
}
}
......@@ -1497,7 +1396,7 @@ size_t Category::getColumnIndex(std::string_view name) const
break;
}
if (VERBOSE and result == mColumns.size() and mCatValidator != nullptr) // validate the name, if it is known at all (since it was not found)
if (VERBOSE > 0 and result == mColumns.size() and mCatValidator != nullptr) // validate the name, if it is known at all (since it was not found)
{
auto iv = mCatValidator->getValidatorForItem(name);
if (iv == nullptr)
......@@ -1543,21 +1442,6 @@ size_t Category::addColumn(std::string_view name)
return result;
}
// RowSet Category::find(Condition&& cond)
// {
// RowSet result(*this);
// cond.prepare(*this);
// for (auto r: *this)
// {
// if (cond(*this, r))
// result.push_back(r);
// }
// return result;
// }
void Category::reorderByIndex()
{
if (mIndex != nullptr)
......@@ -1601,11 +1485,13 @@ std::string Category::getUniqueID(std::function<std::string(int)> generator)
if (mCatValidator != nullptr and mCatValidator->mKeys.size() == 1)
key = mCatValidator->mKeys.front();
size_t nr = size();
// calling size() often is a waste of resources
if (mLastUniqueNr == 0)
mLastUniqueNr = size();
for (;;)
{
std::string result = generator(int(nr++));
std::string result = generator(static_cast<int>(mLastUniqueNr++));
if (exists(Key(key) == result))
continue;
......@@ -1665,21 +1551,6 @@ Row Category::operator[](Condition &&cond)
return result;
}
// RowSet Category::find(Condition&& cond)
// {
// // return RowSet(*this, std::forward<Condition>(cond));
// RowSet result(*this);
// cond.prepare(*this);
// for (auto r: *this)
// {
// if (cond(*this, r))
// result.insert(result.end(), r);
// }
// return result;
// }
bool Category::exists(Condition &&cond) const
{
bool result = false;
......@@ -2340,7 +2211,7 @@ void Category::validateLinks() const
if (not hasParent(r, *parentCat, *link))
++missing;
if (missing)
if (missing and VERBOSE >= 0)
{
std::cerr << "Links for " << link->mLinkGroupLabel << " are incomplete" << std::endl
<< " There are " << missing << " items in " << mName << " that don't have matching parent items in " << parentCat->mName << std::endl;
......@@ -2408,26 +2279,26 @@ std::set<size_t> Category::keyFieldsByIndex() const
bool operator==(const Category &a, const Category &b)
{
using namespace std::placeholders;
using namespace std::placeholders;
bool result = true;
// set<std::string> tagsA(a.fields()), tagsB(b.fields());
//
// if (tagsA != tagsB)
// std::cout << "Unequal number of fields" << std::endl;
// set<std::string> tagsA(a.fields()), tagsB(b.fields());
//
// if (tagsA != tagsB)
// std::cout << "Unequal number of fields" << std::endl;
auto& validator = a.getValidator();
auto &validator = a.getValidator();
auto catValidator = validator.getValidatorForCategory(a.name());
if (catValidator == nullptr)
throw std::runtime_error("missing cat validator");
typedef std::function<int(const char*,const char*)> compType;
std::vector<std::tuple<std::string,compType>> tags;
typedef std::function<int(const char *, const char *)> compType;
std::vector<std::tuple<std::string, compType>> tags;
auto keys = catValidator->mKeys;
std::vector<size_t> keyIx;
for (auto& tag: a.fields())
for (auto &tag : a.fields())
{
auto iv = catValidator->getValidatorForItem(tag);
if (iv == nullptr)
......@@ -2436,24 +2307,25 @@ bool operator==(const Category &a, const Category &b)
if (tv == nullptr)
throw std::runtime_error("missing type validator");
tags.push_back(std::make_tuple(tag, std::bind(&cif::ValidateType::compare, tv, std::placeholders::_1, std::placeholders::_2)));
auto pred = [tag](const std::string& s) -> bool { return cif::iequals(tag, s) == 0; };
auto pred = [tag](const std::string &s) -> bool
{ return cif::iequals(tag, s) == 0; };
if (find_if(keys.begin(), keys.end(), pred) == keys.end())
keyIx.push_back(tags.size() - 1);
}
// a.reorderByIndex();
// b.reorderByIndex();
auto rowEqual = [&](const cif::Row& ra, const cif::Row& rb)
auto rowEqual = [&](const cif::Row &ra, const cif::Row &rb)
{
int d = 0;
for (auto kix: keyIx)
for (auto kix : keyIx)
{
std::string tag;
compType compare;
std::tie(tag, compare) = tags[kix];
d = compare(ra[tag].c_str(), rb[tag].c_str());
......@@ -2465,7 +2337,7 @@ bool operator==(const Category &a, const Category &b)
break;
}
}
return d == 0;
};
......@@ -2483,9 +2355,9 @@ bool operator==(const Category &a, const Category &b)
else
return false;
}
cif::Row ra = *ai, rb = *bi;
if (not rowEqual(ra, rb))
{
if (cif::VERBOSE > 1)
......@@ -2493,21 +2365,25 @@ bool operator==(const Category &a, const Category &b)
else
return false;
}
std::vector<std::string> missingA, missingB, different;
for (auto& tt: tags)
for (auto &tt : tags)
{
std::string tag;
compType compare;
std::tie(tag, compare) = tt;
// make it an option to compare unapplicable to empty or something
const char* ta = ra[tag].c_str(); if (strcmp(ta, ".") == 0 or strcmp(ta, "?") == 0) ta = "";
const char* tb = rb[tag].c_str(); if (strcmp(tb, ".") == 0 or strcmp(tb, "?") == 0) tb = "";
const char *ta = ra[tag].c_str();
if (strcmp(ta, ".") == 0 or strcmp(ta, "?") == 0)
ta = "";
const char *tb = rb[tag].c_str();
if (strcmp(tb, ".") == 0 or strcmp(tb, "?") == 0)
tb = "";
if (compare(ta, tb) != 0)
{
if (cif::VERBOSE > 1)
......@@ -2519,7 +2395,7 @@ bool operator==(const Category &a, const Category &b)
return false;
}
}
++ai;
++bi;
}
......@@ -2527,18 +2403,6 @@ bool operator==(const Category &a, const Category &b)
return result;
}
// auto Category::iterator::operator++() -> iterator&
// {
// mCurrent = Row(mCurrent.data()->mNext);
// return *this;
// }
// auto Category::const_iterator::operator++() -> const_iterator&
// {
// mCurrent = Row(mCurrent.data()->mNext);
// return *this;
// }
namespace detail
{
......@@ -2900,7 +2764,7 @@ void Category::update_value(RowSet &&rows, const std::string &tag, const std::st
}
// cannot update this...
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Cannot update child " << childCat->mName << "." << childTag << " with value " << value << std::endl;
}
......@@ -3009,21 +2873,22 @@ void Row::assign(const Item &value, bool skipUpdateLinked)
assign(value.name(), value.value(), skipUpdateLinked);
}
void Row::assign(std::string_view name, const std::string &value, bool skipUpdateLinked)
void Row::assign(std::string_view name, const std::string &value, bool skipUpdateLinked, bool validate)
{
try
{
auto cat = mData->mCategory;
assign(cat->addColumn(name), value, skipUpdateLinked);
assign(cat->addColumn(name), value, skipUpdateLinked, validate);
}
catch (const std::exception &ex)
{
std::cerr << "Could not assign value '" << value << "' to column _" << mData->mCategory->name() << '.' << name << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Could not assign value '" << value << "' to column _" << mData->mCategory->name() << '.' << name << std::endl;
throw;
}
}
void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked, bool validate)
{
if (mData == nullptr)
throw std::logic_error("invalid Row, no data assigning value '" + value + "' to column with index " + std::to_string(column));
......@@ -3049,7 +2914,7 @@ void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
std::string oldStrValue = oldValue ? oldValue : "";
// check the value
if (col.mValidator)
if (col.mValidator and validate)
(*col.mValidator)(value);
// If the field is part of the Key for this Category, remove it from the index
......@@ -3181,7 +3046,7 @@ void Row::assign(size_t column, const std::string &value, bool skipUpdateLinked)
auto rows_n = childCat->find(std::move(cond_n));
if (not rows_n.empty())
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Will not rename in child category since there are already rows that link to the parent" << std::endl;
continue;
......@@ -3387,7 +3252,7 @@ void Row::swap(size_t cix, ItemRow *a, ItemRow *b)
}
else
{
if (VERBOSE)
if (VERBOSE > 0)
std::cerr << "In " << childCat->mName << " changing " << linkChildColName << ": " << r[linkChildColName].as<std::string>() << " => " << (i ? i->mText : "") << std::endl;
r[linkChildColName] = i ? i->mText : "";
}
......@@ -3496,7 +3361,8 @@ File::File(const std::filesystem::path &path, bool validate)
}
catch (const std::exception &ex)
{
std::cerr << "Error while loading file " << path << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error while loading file " << path << std::endl;
throw;
}
}
......@@ -3564,7 +3430,8 @@ void File::load(const std::filesystem::path &p)
}
catch (const std::exception &ex)
{
std::cerr << "Error loading file " << path << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error loading file " << path << std::endl;
throw;
}
}
......@@ -3660,7 +3527,7 @@ bool File::isValid()
{
if (mValidator == nullptr)
{
if (VERBOSE)
if (VERBOSE > 0)
std::cerr << "No dictionary loaded explicitly, loading default" << std::endl;
loadDictionary();
......
src/Cif2PDB.cpp
......@@ -753,7 +753,8 @@ class Ff : public FBase
}
catch (const std::exception& ex)
{
std::cerr << "Failed to write '" << s << "' as a double, this indicates an error in the code for writing PDB files" << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Failed to write '" << s << "' as a double, this indicates an error in the code for writing PDB files" << std::endl;
os << s;
}
}
......@@ -2329,7 +2330,8 @@ void WriteRemark200(std::ostream& pdbFile, Datablock& db)
}
catch (const std::exception& ex)
{
std::cerr << ex.what() << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << ex.what() << std::endl;
}
}
......@@ -2390,7 +2392,8 @@ void WriteRemark280(std::ostream& pdbFile, Datablock& db)
}
catch (const std::exception& ex)
{
std::cerr << ex.what() << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << ex.what() << std::endl;
}
}
......
src/CifParser.cpp
......@@ -288,7 +288,7 @@ SacParser::CIFToken SacParser::getNextToken()
mState = eStateTextField + 1;
else if (ch == kEOF)
error("unterminated textfield");
else if (not isAnyPrint(ch))
else if (not isAnyPrint(ch) and cif::VERBOSE >= 0)
// error("invalid character in text field '" + string({ static_cast<char>(ch) }) + "' (" + to_string((int)ch) + ")");
std::cerr << "invalid character in text field '" << std::string({static_cast<char>(ch)}) << "' (" << ch << ") line: " << mLineNr << std::endl;
break;
......@@ -1220,7 +1220,7 @@ void DictParser::linkItems()
{
for (auto &iv : cv.mItemValidators)
{
if (iv.mType == nullptr)
if (iv.mType == nullptr and cif::VERBOSE >= 0)
std::cerr << "Missing item_type for " << iv.mTag << std::endl;
}
}
......@@ -1255,7 +1255,8 @@ void DictParser::loadDictionary()
}
catch (const std::exception &)
{
std::cerr << "Error parsing dictionary" << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error parsing dictionary" << std::endl;
throw;
}
......
src/CifUtils.cpp
......@@ -1237,7 +1237,7 @@ std::filesystem::path gDataDir;
void addDataDirectory(std::filesystem::path dataDir)
{
if (VERBOSE and not fs::exists(dataDir))
if (VERBOSE > 0 and not fs::exists(dataDir))
std::cerr << "The specified data directory " << dataDir << " does not exist" << std::endl;
gDataDir = dataDir;
}
......
src/CifValidator.cpp
......@@ -354,7 +354,7 @@ void Validator::reportError(const std::string &msg, bool fatal) const
{
if (mStrict or fatal)
throw ValidationError(msg);
else if (VERBOSE)
else if (VERBOSE > 0)
std::cerr << msg << std::endl;
}
......
src/Compound.cpp
......@@ -193,7 +193,7 @@ Compound::Compound(cif::Datablock &db, const std::string &id, const std::string
bond.type = BondType::delo;
else
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Unimplemented chem_comp_bond.type " << btype << " in " << id << std::endl;
bond.type = BondType::sing;
}
......@@ -403,7 +403,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const std::filesystem::path &file, std:
if (cif::iequals(id, "gly"))
type = "peptide linking";
else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide") or cif::iequals(group, "p-peptide"))
type = "L-peptide linking";
else if (cif::iequals(group, "DNA"))
type = "DNA linking";
......@@ -520,7 +520,7 @@ Compound *CCDCompoundFactoryImpl::create(const std::string &id)
}
}
if (result == nullptr and cif::VERBOSE)
if (result == nullptr and cif::VERBOSE > 0)
std::cerr << "Could not locate compound " << id << " in the CCD components file" << std::endl;
return result;
......@@ -611,7 +611,7 @@ Compound *CCP4CompoundFactoryImpl::create(const std::string &id)
if (cif::iequals(id, "gly"))
type = "peptide linking";
else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide"))
else if (cif::iequals(group, "l-peptide") or cif::iequals(group, "L-peptide linking") or cif::iequals(group, "peptide") or cif::iequals(group, "p-peptide"))
type = "L-peptide linking";
else if (cif::iequals(group, "DNA"))
type = "DNA linking";
......@@ -645,13 +645,13 @@ CompoundFactory::CompoundFactory()
auto ccd = cif::loadResource("components.cif");
if (ccd)
mImpl.reset(new CCDCompoundFactoryImpl(mImpl));
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "CCD components.cif file was not found" << std::endl;
const char *clibd_mon = getenv("CLIBD_MON");
if (clibd_mon != nullptr and fs::is_directory(clibd_mon))
mImpl.reset(new CCP4CompoundFactoryImpl(clibd_mon));
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "CCP4 monomers library not found, CLIBD_MON is not defined" << std::endl;
}
......@@ -695,7 +695,8 @@ void CompoundFactory::setDefaultDictionary(const std::filesystem::path &inDictFi
}
catch (const std::exception &)
{
std::cerr << "Error loading dictionary " << inDictFile << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error loading dictionary " << inDictFile << std::endl;
throw;
}
}
......@@ -715,7 +716,8 @@ void CompoundFactory::pushDictionary(const std::filesystem::path &inDictFile)
}
catch (const std::exception &)
{
std::cerr << "Error loading dictionary " << inDictFile << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error loading dictionary " << inDictFile << std::endl;
throw;
}
}
......
src/PDB2Cif.cpp
......@@ -268,7 +268,8 @@ int PDBRecord::vI(int columnFirst, int columnLast)
}
catch (const std::exception &ex)
{
std::cerr << "Trying to parse '" << std::string(mValue + columnFirst - 7, mValue + columnLast - 7) << '\'' << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Trying to parse '" << std::string(mValue + columnFirst - 7, mValue + columnLast - 7) << '\'' << std::endl;
throw;
}
......@@ -337,7 +338,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
}
else if (not isspace(ch))
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping invalid character in SOURCE ID: " << ch << std::endl;
}
break;
......@@ -354,7 +355,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
case eColon:
if (ch == ';')
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Empty value for SOURCE: " << id << std::endl;
state = eStart;
}
......@@ -418,7 +419,7 @@ std::tuple<std::string, std::string> SpecificationListParser::GetNextSpecificati
case eError:
if (ch == ';')
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Skipping invalid header line: '" << std::string(start, mP) << std::endl;
state = eStart;
}
......@@ -832,7 +833,7 @@ class PDBFileParser
if (not mChainSeq2AsymSeq.count(key))
{
ec = error::make_error_code(error::pdbErrors::residueNotFound);
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Residue " << chainID << resSeq << iCode << " could not be mapped" << std::endl;
}
else
......@@ -929,7 +930,7 @@ class PDBFileParser
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << ex.what() << std::endl;
ec = error::make_error_code(error::pdbErrors::invalidDate);
}
......@@ -1160,7 +1161,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
if (is.eof())
break;
std::cerr << "Line number " << lineNr << " is empty!" << std::endl;
if (cif::VERBOSE > 0)
std::cerr << "Line number " << lineNr << " is empty!" << std::endl;
getline(is, lookahead);
++lineNr;
......@@ -1278,7 +1280,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
}
catch (const std::exception &ex)
{
std::cerr << "Dropping FORMUL line (" << (lineNr - 1) << ") with invalid component number '" << value.substr(1, 3) << '\'' << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Dropping FORMUL line (" << (lineNr - 1) << ") with invalid component number '" << value.substr(1, 3) << '\'' << std::endl;
continue;
// throw_with_nested(std::runtime_error("Invalid component number '" + value.substr(1, 3) + '\''));
}
......@@ -1305,7 +1308,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
}
catch (const std::exception &ex)
{
std::cerr << "Error parsing FORMUL at line " << lineNr << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error parsing FORMUL at line " << lineNr << std::endl;
throw;
}
}
......@@ -1412,7 +1416,8 @@ void PDBFileParser::PreParseInput(std::istream &is)
if (not dropped.empty())
{
std::cerr << "Dropped unsupported records: " << ba::join(dropped, ", ") << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Dropped unsupported records: " << ba::join(dropped, ", ") << std::endl;
}
if (mData == nullptr)
......@@ -1440,7 +1445,7 @@ void PDBFileParser::Match(const std::string &expected, bool throwIfMissing)
{
if (throwIfMissing)
throw std::runtime_error("Expected record " + expected + " but found " + mRec->mName);
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Expected record " + expected + " but found " + mRec->mName << std::endl;
}
}
......@@ -1575,7 +1580,8 @@ void PDBFileParser::ParseTitle()
if (not iequals(key, "MOL_ID") and mCompounds.empty())
{
std::cerr << "Ignoring invalid COMPND record" << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Ignoring invalid COMPND record" << std::endl;
break;
}
......@@ -1625,7 +1631,7 @@ void PDBFileParser::ParseTitle()
// auto colon = s.find(": ");
// if (colon == std::string::npos)
// {
// if (cif::VERBOSE)
// if (cif::VERBOSE > 0)
// std::cerr << "invalid source field, missing colon (" << s << ')' << std::endl;
// continue;
// }
......@@ -1713,7 +1719,7 @@ void PDBFileParser::ParseTitle()
// NUMMDL
if (mRec->is("NUMMDL"))
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping unimplemented NUMMDL record" << std::endl;
GetNextRecord();
}
......@@ -1816,7 +1822,7 @@ void PDBFileParser::ParseTitle()
// SPRSDE
if (mRec->is("SPRSDE"))
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping unimplemented SPRSDE record" << std::endl;
GetNextRecord();
}
......@@ -2265,7 +2271,7 @@ void PDBFileParser::ParseRemarks()
state = eMCP;
else if (subtopic == "CHIRAL CENTERS")
state = eChC;
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
throw std::runtime_error("Unknown subtopic in REMARK 500: " + subtopic);
headerSeen = false;
......@@ -2342,7 +2348,7 @@ void PDBFileParser::ParseRemarks()
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping REMARK 500 at line " << mRec->mLineNr << " due to invalid symmetry operation" << std::endl;
continue;
}
......@@ -2675,7 +2681,7 @@ void PDBFileParser::ParseRemarks()
case sStart:
if (s == "SITE")
state = sID;
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
throw std::runtime_error("Invalid REMARK 800 record, expected SITE");
break;
......@@ -2685,7 +2691,7 @@ void PDBFileParser::ParseRemarks()
id = m[1].str();
state = sEvidence;
}
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
throw std::runtime_error("Invalid REMARK 800 record, expected SITE_IDENTIFIER");
break;
......@@ -2695,7 +2701,7 @@ void PDBFileParser::ParseRemarks()
evidence = m[1].str();
state = sDesc;
}
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
throw std::runtime_error("Invalid REMARK 800 record, expected SITE_IDENTIFIER");
break;
......@@ -2906,7 +2912,7 @@ void PDBFileParser::ParseRemark200()
collectionDate = pdb2cifDate(rm200("DATE OF DATA COLLECTION", diffrnNr), ec);
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << ec.message() << " for pdbx_collection_date" << std::endl;
// The date field can become truncated when multiple values are available
......@@ -3025,7 +3031,7 @@ void PDBFileParser::ParseRemark200()
else if (inRM200({"HIGHEST RESOLUTION SHELL, RANGE LOW (A)", "COMPLETENESS FOR SHELL (%)",
"R MERGE FOR SHELL (I)", "R SYM FOR SHELL (I)", "<I/SIGMA(I)> FOR SHELL", "DATA REDUNDANCY IN SHELL"}))
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Not writing reflns_shell record since d_res_high is missing" << std::endl;
}
}
......@@ -3595,7 +3601,7 @@ void PDBFileParser::ConstructEntities()
{
auto &r = chain.mResiduesSeen[lastResidueIndex + 1];
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cerr << "Detected residues that cannot be aligned to SEQRES" << std::endl
<< "First residue is " << chain.mDbref.chainID << ':' << r.mSeqNum << r.mIcode << std::endl;
......@@ -4005,7 +4011,7 @@ void PDBFileParser::ConstructEntities()
tie(asym, labelSeq, std::ignore) = MapResidue(seqadv.chainID, seqadv.seqNum, seqadv.iCode, ec);
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "dropping unmatched SEQADV record" << std::endl;
continue;
}
......@@ -4319,7 +4325,7 @@ void PDBFileParser::ConstructEntities()
tie(asymID, seq, std::ignore) = MapResidue(chainID, seqNum, iCode, ec);
if (ec) // no need to write a modres if it could not be found
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "dropping unmapped MODRES record" << std::endl;
continue;
}
......@@ -4415,7 +4421,7 @@ void PDBFileParser::ConstructEntities()
tie(asymID, seqNr, isPolymer) = MapResidue(unobs.chain, unobs.seq, unobs.iCode, ec);
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "error mapping unobserved residue" << std::endl;
continue;
}
......@@ -4676,7 +4682,7 @@ void PDBFileParser::ParseSecondaryStructure()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Could not map residue for HELIX " << vI(8, 10) << std::endl;
}
else
......@@ -4791,7 +4797,7 @@ void PDBFileParser::ParseSecondaryStructure()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping SHEET record " << vI(8, 10) << std::endl;
}
else
......@@ -4827,7 +4833,7 @@ void PDBFileParser::ParseSecondaryStructure()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping unmatched pdbx_struct_sheet_hbond record" << std::endl;
}
else
......@@ -4927,7 +4933,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping SSBOND " << vI(8, 10) << std::endl;
continue;
}
......@@ -4948,7 +4954,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping SSBOND " << vI(8, 10) << " due to invalid symmetry operation" << std::endl;
continue;
}
......@@ -4993,7 +4999,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
if (mRec->is("LINK ") or mRec->is("LINKR "))
{
if (cif::VERBOSE and mRec->is("LINKR "))
if (cif::VERBOSE > 0 and mRec->is("LINKR "))
std::cerr << "Accepting non-standard LINKR record, but ignoring extra information" << std::endl;
// 1 - 6 Record name "LINK "
......@@ -5046,7 +5052,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping LINK record at line " << mRec->mLineNr << std::endl;
continue;
}
......@@ -5062,7 +5068,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
}
catch (const std::invalid_argument &)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Distance value '" << distance << "' is not a valid float in LINK record" << std::endl;
swap(ccp4LinkID, distance); // assume this is a ccp4_link_id... oh really?
}
......@@ -5078,7 +5084,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping LINK record at line " << mRec->mLineNr << " due to invalid symmetry operation" << std::endl;
continue;
}
......@@ -5149,7 +5155,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Dropping CISPEP record at line " << mRec->mLineNr << std::endl;
continue;
}
......@@ -5215,7 +5221,7 @@ void PDBFileParser::ParseMiscellaneousFeatures()
if (ec)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping struct_site_gen record" << std::endl;
}
else
......@@ -5518,7 +5524,7 @@ void PDBFileParser::ParseCoordinate(int modelNr)
{
if (groupPDB == "HETATM")
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Changing atom from HETATM to ATOM at line " << mRec->mLineNr << std::endl;
groupPDB = "ATOM";
}
......@@ -5527,7 +5533,7 @@ void PDBFileParser::ParseCoordinate(int modelNr)
{
if (groupPDB == "ATOM")
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Changing atom from ATOM to HETATM at line " << mRec->mLineNr << std::endl;
groupPDB = "HETATM";
}
......@@ -5698,7 +5704,8 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
}
catch (const std::exception &ex)
{
std::cerr << "Error parsing REMARK 3" << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error parsing REMARK 3" << std::endl;
throw;
}
//
......@@ -5750,12 +5757,12 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
if ((symm1.empty() or symm1 == "1_555") and (symm2.empty() or symm2 == "1_555"))
distance = static_cast<float>(mmcif::Distance(mmcif::Point{x1, y1, z1}, mmcif::Point{x2, y2, z2}));
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "Cannot calculate distance for link since one of the atoms is in another dimension" << std::endl;
}
catch (std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Error finding atom for LINK distance calculation: " << ex.what() << std::endl;
}
......@@ -5764,10 +5771,13 @@ void PDBFileParser::Parse(std::istream &is, cif::File &result)
}
catch (const std::exception &ex)
{
std::cerr << "Error parsing PDB";
if (mRec != nullptr)
std::cerr << " at line " << mRec->mLineNr;
std::cerr << std::endl;
if (cif::VERBOSE >= 0)
{
std::cerr << "Error parsing PDB";
if (mRec != nullptr)
std::cerr << " at line " << mRec->mLineNr;
std::cerr << std::endl;
}
throw;
}
}
......@@ -5947,7 +5957,7 @@ int PDBFileParser::PDBChain::AlignResToSeqRes()
switch (tb(x, y))
{
case -1:
// if (cif::VERBOSE)
// if (cif::VERBOSE > 0)
// std::cerr << "A residue found in the ATOM records "
// << "(" << ry[y].mMonID << " @ " << mDbref.chainID << ":" << ry[y].mSeqNum
// << ((ry[y].mIcode == ' ' or ry[y].mIcode == 0) ? "" : std::string{ ry[y].mIcode }) << ")"
......@@ -6034,7 +6044,6 @@ void ReadPDBFile(std::istream &pdbFile, cif::File &cifFile)
p.Parse(pdbFile, cifFile);
if (not cifFile.isValid())
// throw std::runtime_error("Resulting mmCIF file is invalid");
if (not cifFile.isValid() and cif::VERBOSE >= 0)
std::cerr << "Resulting mmCIF file is not valid!" << std::endl;
}
src/PDB2CifRemark3.cpp
......@@ -1320,7 +1320,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
if (line != "REFINEMENT.")
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Unexpected data in REMARK 3" << std::endl;
return false;
}
......@@ -1332,7 +1332,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
if (not std::regex_match(line, m, rxp))
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Expected valid PROGRAM line in REMARK 3" << std::endl;
return false;
}
......@@ -1367,8 +1367,9 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
}
catch(const std::exception& e)
{
std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
<< e.what() << '\n';
if (cif::VERBOSE >= 0)
std::cerr << "Error parsing REMARK 3 with " << parser->program() << std::endl
<< e.what() << '\n';
score = 0;
}
......@@ -1411,7 +1412,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
tryParser(new TNT_Remark3Parser(program, expMethod, r, db));
else if (ba::starts_with(program, "X-PLOR"))
tryParser(new XPLOR_Remark3Parser(program, expMethod, r, db));
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "Skipping unknown program (" << program << ") in REMARK 3" << std::endl;
}
......@@ -1420,7 +1421,8 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
bool guessProgram = scores.empty() or scores.front().score < 0.9f;;
if (guessProgram)
{
std::cerr << "Unknown or untrusted program in REMARK 3, trying all parsers to see if there is a match" << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Unknown or untrusted program in REMARK 3, trying all parsers to see if there is a match" << std::endl;
tryParser(new BUSTER_TNT_Remark3Parser("BUSTER-TNT", expMethod, r, db));
tryParser(new CNS_Remark3Parser("CNS", expMethod, r, db));
......@@ -1444,7 +1446,7 @@ bool Remark3Parser::parse(const std::string& expMethod, PDBRecord* r, cif::Datab
auto& best = scores.front();
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Choosing " << best.parser->program() << " version '" << best.parser->version() << "' as refinement program. Score = " << best.score << std::endl;
auto& software = db["software"];
......
src/Secondary.cpp
/*-
* 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
......@@ -26,14 +26,14 @@
// Calculate DSSP-like secondary structure information
#include <numeric>
#include <iomanip>
#include <numeric>
#include <thread>
#include <boost/algorithm/string.hpp>
#include "cif++/Structure.hpp"
#include "cif++/Secondary.hpp"
#include "cif++/Structure.hpp"
namespace ba = boost::algorithm;
......@@ -47,37 +47,37 @@ struct Res;
enum ResidueType
{
kUnknownResidue,
//
kAlanine, // A ala
kArginine, // R arg
kAsparagine, // N asn
kAsparticAcid, // D asp
kCysteine, // C cys
kGlutamicAcid, // E glu
kGlutamine, // Q gln
kGlycine, // G gly
kHistidine, // H his
kIsoleucine, // I ile
kLeucine, // L leu
kLysine, // K lys
kMethionine, // M met
kPhenylalanine, // F phe
kProline, // P pro
kSerine, // S ser
kThreonine, // T thr
kTryptophan, // W trp
kTyrosine, // Y tyr
kValine, // V val
kAlanine, // A ala
kArginine, // R arg
kAsparagine, // N asn
kAsparticAcid, // D asp
kCysteine, // C cys
kGlutamicAcid, // E glu
kGlutamine, // Q gln
kGlycine, // G gly
kHistidine, // H his
kIsoleucine, // I ile
kLeucine, // L leu
kLysine, // K lys
kMethionine, // M met
kPhenylalanine, // F phe
kProline, // P pro
kSerine, // S ser
kThreonine, // T thr
kTryptophan, // W trp
kTyrosine, // Y tyr
kValine, // V val
kResidueTypeCount
};
struct ResidueInfo
{
ResidueType type;
char code;
char name[4];
ResidueType type;
char code;
char name[4];
};
const ResidueInfo kResidueInfo[] = {
......@@ -109,8 +109,8 @@ ResidueType MapResidue(std::string inName)
ba::trim(inName);
ResidueType result = kUnknownResidue;
for (auto& ri: kResidueInfo)
for (auto &ri : kResidueInfo)
{
if (inName == ri.name)
{
......@@ -118,37 +118,39 @@ ResidueType MapResidue(std::string inName)
break;
}
}
return result;
}
struct HBond
{
Res* residue;
double energy;
Res *residue;
double energy;
};
enum BridgeType
{
btNoBridge, btParallel, btAntiParallel
btNoBridge,
btParallel,
btAntiParallel
};
struct Bridge
{
BridgeType type;
uint32_t sheet, ladder;
std::set<Bridge*> link;
std::deque<uint32_t> i, j;
std::string chainI, chainJ;
bool operator<(const Bridge& b) const { return chainI < b.chainI or (chainI == b.chainI and i.front() < b.i.front()); }
BridgeType type;
uint32_t sheet, ladder;
std::set<Bridge *> link;
std::deque<uint32_t> i, j;
std::string chainI, chainJ;
bool operator<(const Bridge &b) const { return chainI < b.chainI or (chainI == b.chainI and i.front() < b.i.front()); }
};
struct BridgeParner
{
Res* residue;
uint32_t ladder;
bool parallel;
Res *residue;
uint32_t ladder;
bool parallel;
};
// --------------------------------------------------------------------
......@@ -159,7 +161,7 @@ const float
kMinimalCADistance = 9.0f,
kMinHBondEnergy = -9.9f,
kMaxHBondEnergy = -0.5f,
kCouplingConstant = -27.888f, // = -332 * 0.42 * 0.2
kCouplingConstant = -27.888f, // = -332 * 0.42 * 0.2
kMaxPeptideBondLength = 2.5f;
const float
......@@ -172,18 +174,19 @@ const float
struct Res
{
Res(const Monomer& m, int nr, ChainBreak brk)
: mM(m), mNumber(nr)
Res(const Monomer &m, int nr, ChainBreak brk)
: mM(m)
, mNumber(nr)
, mType(MapResidue(m.compoundID()))
, mChainBreak(brk)
{
// update the box containing all atoms
mBox[0].mX = mBox[0].mY = mBox[0].mZ = std::numeric_limits<float>::max();
mBox[0].mX = mBox[0].mY = mBox[0].mZ = std::numeric_limits<float>::max();
mBox[1].mX = mBox[1].mY = mBox[1].mZ = -std::numeric_limits<float>::max();
mH = mmcif::Point{ std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max() };
mH = mmcif::Point{std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()};
for (auto& a: mM.unique_atoms())
for (auto &a : mM.unique_atoms())
{
if (a.labelAtomID() == "CA")
{
......@@ -227,44 +230,44 @@ struct Res
{
// assign the Hydrogen
mH = mN;
if (mType != kProline and mPrev != nullptr)
{
auto pc = mPrev->mC;
auto po = mPrev->mO;
float CODistance = static_cast<float>(Distance(pc, po));
mH.mX += (pc.mX - po.mX) / CODistance;
mH.mY += (pc.mY - po.mY) / CODistance;
mH.mZ += (pc.mZ - po.mZ) / CODistance;
mH.mX += (pc.mX - po.mX) / CODistance;
mH.mY += (pc.mY - po.mY) / CODistance;
mH.mZ += (pc.mZ - po.mZ) / CODistance;
}
}
void SetSecondaryStructure(SecondaryStructureType inSS) { mSecondaryStructure = inSS; }
SecondaryStructureType GetSecondaryStructure() const { return mSecondaryStructure; }
void SetBetaPartner(uint32_t n, Res& inResidue, uint32_t inLadder, bool inParallel)
void SetSecondaryStructure(SecondaryStructureType inSS) { mSecondaryStructure = inSS; }
SecondaryStructureType GetSecondaryStructure() const { return mSecondaryStructure; }
void SetBetaPartner(uint32_t n, Res &inResidue, uint32_t inLadder, bool inParallel)
{
assert(n == 0 or n == 1);
mBetaPartner[n].residue = &inResidue;
mBetaPartner[n].ladder = inLadder;
mBetaPartner[n].parallel = inParallel;
}
BridgeParner GetBetaPartner(uint32_t n) const
{
assert(n == 0 or n == 1);
return mBetaPartner[n];
}
void SetSheet(uint32_t inSheet) { mSheet = inSheet; }
uint32_t GetSheet() const { return mSheet; }
bool IsBend() const { return mBend; }
void SetBend(bool inBend) { mBend = inBend; }
void SetSheet(uint32_t inSheet) { mSheet = inSheet; }
uint32_t GetSheet() const { return mSheet; }
bool IsBend() const { return mBend; }
void SetBend(bool inBend) { mBend = inBend; }
Helix GetHelixFlag(HelixType helixType) const
{
size_t stride = static_cast<size_t>(helixType);
......@@ -292,7 +295,7 @@ struct Res
throw std::runtime_error("Only cysteine residues can form sulphur bridges");
mSSBridgeNr = inBridgeNr;
}
uint8_t GetSSBridgeNr() const
{
if (mType != kCysteine)
......@@ -300,21 +303,20 @@ struct Res
return mSSBridgeNr;
}
double CalculateSurface(const std::vector<Res>& inResidues);
double CalculateSurface(const Point& inAtom, float inRadius, const std::vector<Res*>& inNeighbours);
double CalculateSurface(const std::vector<Res> &inResidues);
double CalculateSurface(const Point &inAtom, float inRadius, const std::vector<Res *> &inNeighbours);
bool AtomIntersectsBox(const Point& atom, float inRadius) const
bool AtomIntersectsBox(const Point &atom, float inRadius) const
{
return
atom.mX + inRadius >= mBox[0].mX and
atom.mX - inRadius <= mBox[1].mX and
atom.mY + inRadius >= mBox[0].mY and
atom.mY - inRadius <= mBox[1].mY and
atom.mZ + inRadius >= mBox[0].mZ and
atom.mZ - inRadius <= mBox[1].mZ;
return atom.mX + inRadius >= mBox[0].mX and
atom.mX - inRadius <= mBox[1].mX and
atom.mY + inRadius >= mBox[0].mY and
atom.mY - inRadius <= mBox[1].mY and
atom.mZ + inRadius >= mBox[0].mZ and
atom.mZ - inRadius <= mBox[1].mZ;
}
void ExtendBox(const Point& atom, float inRadius)
void ExtendBox(const Point &atom, float inRadius)
{
if (mBox[0].mX > atom.mX - inRadius)
mBox[0].mX = atom.mX - inRadius;
......@@ -330,10 +332,10 @@ struct Res
mBox[1].mZ = atom.mZ + inRadius;
}
Res* mNext = nullptr;
Res* mPrev = nullptr;
Res *mNext = nullptr;
Res *mPrev = nullptr;
const Monomer& mM;
const Monomer &mM;
std::string mAltID;
int mNumber;
......@@ -344,14 +346,14 @@ struct Res
Point mCenter;
std::vector<Point> mSideChain;
double mAccessibility = 0;
ResidueType mType;
uint8_t mSSBridgeNr = 0;
SecondaryStructureType mSecondaryStructure = ssLoop;
HBond mHBondDonor[2] = {}, mHBondAcceptor[2] = {};
BridgeParner mBetaPartner[2] = {};
uint32_t mSheet = 0;
Helix mHelixFlags[4] = { Helix::None, Helix::None, Helix::None, Helix::None }; //
Helix mHelixFlags[4] = {Helix::None, Helix::None, Helix::None, Helix::None}; //
bool mBend = false;
ChainBreak mChainBreak = ChainBreak::None;
};
......@@ -360,19 +362,20 @@ struct Res
class Accumulator
{
public:
public:
struct candidate
{
Point location;
double radius;
double distance;
Point location;
double radius;
double distance;
bool operator<(const candidate& rhs) const
{ return distance < rhs.distance; }
bool operator<(const candidate &rhs) const
{
return distance < rhs.distance;
}
};
void operator()(const Point& a, const Point& b, double d, double r)
void operator()(const Point &a, const Point &b, double d, double r)
{
double distance = DistanceSquared(a, b);
......@@ -384,7 +387,7 @@ class Accumulator
if (distance < test and distance > 0.0001)
{
candidate c = { b - a, r * r, distance };
candidate c = {b - a, r * r, distance};
m_x.push_back(c);
push_heap(m_x.begin(), m_x.end());
......@@ -396,19 +399,18 @@ class Accumulator
sort_heap(m_x.begin(), m_x.end());
}
std::vector<candidate> m_x;
std::vector<candidate> m_x;
};
// we use a fibonacci sphere to calculate the even distribution of the dots
class MSurfaceDots
{
public:
static MSurfaceDots &Instance();
static MSurfaceDots& Instance();
size_t size() const { return mPoints.size(); }
const Point& operator[](size_t inIx) const { return mPoints[inIx]; }
double weight() const { return mWeight; }
size_t size() const { return mPoints.size(); }
const Point &operator[](size_t inIx) const { return mPoints[inIx]; }
double weight() const { return mWeight; }
private:
MSurfaceDots(int32_t inN);
......@@ -417,7 +419,7 @@ class MSurfaceDots
double mWeight;
};
MSurfaceDots& MSurfaceDots::Instance()
MSurfaceDots &MSurfaceDots::Instance()
{
const int32_t kN = 200;
......@@ -442,11 +444,11 @@ MSurfaceDots::MSurfaceDots(int32_t N)
}
}
double Res::CalculateSurface(const Point& inAtom, float inRadius, const std::vector<Res*>& inNeighbours)
double Res::CalculateSurface(const Point &inAtom, float inRadius, const std::vector<Res *> &inNeighbours)
{
Accumulator accumulate;
for (auto r: inNeighbours)
for (auto r : inNeighbours)
{
if (r->AtomIntersectsBox(inAtom, inRadius))
{
......@@ -455,7 +457,7 @@ double Res::CalculateSurface(const Point& inAtom, float inRadius, const std::vec
accumulate(inAtom, r->mC, inRadius, kRadiusC);
accumulate(inAtom, r->mO, inRadius, kRadiusO);
for (auto& atom: r->mSideChain)
for (auto &atom : r->mSideChain)
accumulate(inAtom, atom, inRadius, kRadiusSideAtom);
}
}
......@@ -465,7 +467,7 @@ double Res::CalculateSurface(const Point& inAtom, float inRadius, const std::vec
float radius = inRadius + kRadiusWater;
double surface = 0;
MSurfaceDots& surfaceDots = MSurfaceDots::Instance();
MSurfaceDots &surfaceDots = MSurfaceDots::Instance();
for (size_t i = 0; i < surfaceDots.size(); ++i)
{
......@@ -482,51 +484,51 @@ double Res::CalculateSurface(const Point& inAtom, float inRadius, const std::vec
return surface * radius * radius;
}
double Res::CalculateSurface(const std::vector<Res>& inResidues)
double Res::CalculateSurface(const std::vector<Res> &inResidues)
{
std::vector<Res*> neighbours;
std::vector<Res *> neighbours;
for (auto& r: inResidues)
for (auto &r : inResidues)
{
Point center = r.mCenter;
double radius = r.mRadius;
if (Distance(mCenter, center) < mRadius + radius)
neighbours.push_back(const_cast<Res*>(&r));
neighbours.push_back(const_cast<Res *>(&r));
}
mAccessibility = CalculateSurface(mN, kRadiusN, neighbours) +
CalculateSurface(mCAlpha, kRadiusCA, neighbours) +
CalculateSurface(mC, kRadiusC, neighbours) +
CalculateSurface(mO, kRadiusO, neighbours);
CalculateSurface(mCAlpha, kRadiusCA, neighbours) +
CalculateSurface(mC, kRadiusC, neighbours) +
CalculateSurface(mO, kRadiusO, neighbours);
for (auto& atom: mSideChain)
for (auto &atom : mSideChain)
mAccessibility += CalculateSurface(atom, kRadiusSideAtom, neighbours);
return mAccessibility;
}
void CalculateAccessibilities(std::vector<Res>& inResidues, DSSP_Statistics& stats)
void CalculateAccessibilities(std::vector<Res> &inResidues, DSSP_Statistics &stats)
{
stats.accessibleSurface = 0;
for (auto& residue: inResidues)
for (auto &residue : inResidues)
stats.accessibleSurface += residue.CalculateSurface(inResidues);
}
// --------------------------------------------------------------------
// TODO: use the angle to improve bond energy calculation.
double CalculateHBondEnergy(Res& inDonor, Res& inAcceptor)
double CalculateHBondEnergy(Res &inDonor, Res &inAcceptor)
{
double result = 0;
if (inDonor.mType != kProline)
{
double distanceHO = Distance(inDonor.mH, inAcceptor.mO);
double distanceHC = Distance(inDonor.mH, inAcceptor.mC);
double distanceNC = Distance(inDonor.mN, inAcceptor.mC);
double distanceNO = Distance(inDonor.mN, inAcceptor.mO);
if (distanceHO < kMinimalDistance or distanceHC < kMinimalDistance or distanceNC < kMinimalDistance or distanceNO < kMinimalDistance)
result = kMinHBondEnergy;
else
......@@ -550,7 +552,7 @@ double CalculateHBondEnergy(Res& inDonor, Res& inAcceptor)
{
inDonor.mHBondAcceptor[1].residue = &inAcceptor;
inDonor.mHBondAcceptor[1].energy = result;
}
}
// and acceptor
if (result < inAcceptor.mHBondDonor[0].energy)
......@@ -563,25 +565,24 @@ double CalculateHBondEnergy(Res& inDonor, Res& inAcceptor)
{
inAcceptor.mHBondDonor[1].residue = &inDonor;
inAcceptor.mHBondDonor[1].energy = result;
}
}
return result;
}
// --------------------------------------------------------------------
void CalculateHBondEnergies(std::vector<Res>& inResidues)
void CalculateHBondEnergies(std::vector<Res> &inResidues)
{
// Calculate the HBond energies
for (uint32_t i = 0; i + 1 < inResidues.size(); ++i)
{
auto& ri = inResidues[i];
auto &ri = inResidues[i];
for (uint32_t j = i + 1; j < inResidues.size(); ++j)
{
auto& rj = inResidues[j];
auto &rj = inResidues[j];
if (Distance(ri.mCAlpha, rj.mCAlpha) < kMinimalCADistance)
{
CalculateHBondEnergy(ri, rj);
......@@ -594,7 +595,7 @@ void CalculateHBondEnergies(std::vector<Res>& inResidues)
// --------------------------------------------------------------------
bool NoChainBreak(const Res* a, const Res* b)
bool NoChainBreak(const Res *a, const Res *b)
{
bool result = a->mM.asymID() == b->mM.asymID();
for (auto r = a; result and r != b; r = r->mNext)
......@@ -608,35 +609,39 @@ bool NoChainBreak(const Res* a, const Res* b)
return result;
}
bool NoChainBreak(const Res& a, const Res& b)
bool NoChainBreak(const Res &a, const Res &b)
{
return NoChainBreak(&a, &b);
}
// --------------------------------------------------------------------
bool TestBond(const Res* a, const Res* b)
bool TestBond(const Res *a, const Res *b)
{
return (a->mHBondAcceptor[0].residue == b and a->mHBondAcceptor[0].energy < kMaxHBondEnergy) or
(a->mHBondAcceptor[1].residue == b and a->mHBondAcceptor[1].energy < kMaxHBondEnergy);
}
bool TestBond(DSSP::ResidueInfo const &a, DSSP::ResidueInfo const &b)
{
return
(a->mHBondAcceptor[0].residue == b and a->mHBondAcceptor[0].energy < kMaxHBondEnergy) or
(a->mHBondAcceptor[1].residue == b and a->mHBondAcceptor[1].energy < kMaxHBondEnergy);
return a and b and TestBond(a.mImpl, b.mImpl);
}
// --------------------------------------------------------------------
BridgeType TestBridge(const Res& r1, const Res& r2)
{ // I. a d II. a d parallel
auto a = r1.mPrev; // \ /
auto b = &r1; // b e b e
auto c = r1.mNext; // / \ ..
auto d = r2.mPrev; // c f c f
auto e = &r2; //
auto f = r2.mNext; // III. a <- f IV. a f antiparallel
//
BridgeType result = btNoBridge; // b e b <-> e
//
// c -> d c d
BridgeType TestBridge(const Res &r1, const Res &r2)
{ // I. a d II. a d parallel
auto a = r1.mPrev; // \ /
auto b = &r1; // b e b e
auto c = r1.mNext; // / \ ..
auto d = r2.mPrev; // c f c f
auto e = &r2; //
auto f = r2.mNext; // III. a <- f IV. a f antiparallel
//
BridgeType result = btNoBridge; // b e b <-> e
//
// c -> d c d
if (a and c and NoChainBreak(a, c) and d and f and NoChainBreak(d, f))
{
if ((TestBond(c, e) and TestBond(e, a)) or (TestBond(f, b) and TestBond(b, d)))
......@@ -644,83 +649,80 @@ BridgeType TestBridge(const Res& r1, const Res& r2)
else if ((TestBond(c, d) and TestBond(f, a)) or (TestBond(e, b) and TestBond(b, e)))
result = btAntiParallel;
}
return result;
}
// --------------------------------------------------------------------
// return true if any of the residues in bridge a is identical to any of the residues in bridge b
bool Linked(const Bridge& a, const Bridge& b)
bool Linked(const Bridge &a, const Bridge &b)
{
return
find_first_of(a.i.begin(), a.i.end(), b.i.begin(), b.i.end()) != a.i.end() or
find_first_of(a.i.begin(), a.i.end(), b.j.begin(), b.j.end()) != a.i.end() or
find_first_of(a.j.begin(), a.j.end(), b.i.begin(), b.i.end()) != a.j.end() or
find_first_of(a.j.begin(), a.j.end(), b.j.begin(), b.j.end()) != a.j.end();
return find_first_of(a.i.begin(), a.i.end(), b.i.begin(), b.i.end()) != a.i.end() or
find_first_of(a.i.begin(), a.i.end(), b.j.begin(), b.j.end()) != a.i.end() or
find_first_of(a.j.begin(), a.j.end(), b.i.begin(), b.i.end()) != a.j.end() or
find_first_of(a.j.begin(), a.j.end(), b.j.begin(), b.j.end()) != a.j.end();
}
// --------------------------------------------------------------------
void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
void CalculateBetaSheets(std::vector<Res> &inResidues, DSSP_Statistics &stats)
{
// Calculate Bridges
std::vector<Bridge> bridges;
if (inResidues.size() > 4)
for (uint32_t i = 1; i + 4 < inResidues.size(); ++i)
{
for (uint32_t i = 1; i + 4 < inResidues.size(); ++i)
auto &ri = inResidues[i];
for (uint32_t j = i + 3; j + 1 < inResidues.size(); ++j)
{
auto& ri = inResidues[i];
for (uint32_t j = i + 3; j + 1 < inResidues.size(); ++j)
auto &rj = inResidues[j];
BridgeType type = TestBridge(ri, rj);
if (type == btNoBridge)
continue;
bool found = false;
for (Bridge &bridge : bridges)
{
auto& rj = inResidues[j];
BridgeType type = TestBridge(ri, rj);
if (type == btNoBridge)
if (type != bridge.type or i != bridge.i.back() + 1)
continue;
bool found = false;
for (Bridge& bridge : bridges)
if (type == btParallel and bridge.j.back() + 1 == j)
{
if (type != bridge.type or i != bridge.i.back() + 1)
continue;
if (type == btParallel and bridge.j.back() + 1 == j)
{
bridge.i.push_back(i);
bridge.j.push_back(j);
found = true;
break;
}
if (type == btAntiParallel and bridge.j.front() - 1 == j)
{
bridge.i.push_back(i);
bridge.j.push_front(j);
found = true;
break;
}
bridge.i.push_back(i);
bridge.j.push_back(j);
found = true;
break;
}
if (not found)
if (type == btAntiParallel and bridge.j.front() - 1 == j)
{
Bridge bridge = {};
bridge.type = type;
bridge.i.push_back(i);
bridge.chainI = ri.mM.asymID();
bridge.j.push_back(j);
bridge.chainJ = rj.mM.asymID();
bridges.push_back(bridge);
bridge.j.push_front(j);
found = true;
break;
}
}
if (not found)
{
Bridge bridge = {};
bridge.type = type;
bridge.i.push_back(i);
bridge.chainI = ri.mM.asymID();
bridge.j.push_back(j);
bridge.chainJ = rj.mM.asymID();
bridges.push_back(bridge);
}
}
}
// extend ladders
std::sort(bridges.begin(), bridges.end());
for (uint32_t i = 0; i < bridges.size(); ++i)
{
for (uint32_t j = i + 1; j < bridges.size(); ++j)
......@@ -742,7 +744,7 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
{
continue;
}
bool bulge;
if (bridges[i].type == btParallel)
bulge = ((jbj - jei < 6 and ibj - iei < 3) or (jbj - jei < 3));
......@@ -763,25 +765,25 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
}
// Sheet
std::set<Bridge*> ladderset;
for (Bridge& bridge : bridges)
std::set<Bridge *> ladderset;
for (Bridge &bridge : bridges)
{
ladderset.insert(&bridge);
size_t n = bridge.i.size();
if (n > kHistogramSize)
n = kHistogramSize;
if (bridge.type == btParallel)
stats.parallelBridgesPerLadderHistogram[n - 1] += 1;
else
stats.antiparallelBridgesPerLadderHistogram[n - 1] += 1;
}
uint32_t sheet = 1, ladder = 0;
while (not ladderset.empty())
{
std::set<Bridge*> sheetset;
std::set<Bridge *> sheetset;
sheetset.insert(*ladderset.begin());
ladderset.erase(ladderset.begin());
......@@ -789,9 +791,9 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
while (not done)
{
done = true;
for (Bridge* a : sheetset)
for (Bridge *a : sheetset)
{
for (Bridge* b : ladderset)
for (Bridge *b : ladderset)
{
if (Linked(*a, *b))
{
......@@ -806,15 +808,15 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
}
}
for (Bridge* bridge : sheetset)
for (Bridge *bridge : sheetset)
{
bridge->ladder = ladder;
bridge->sheet = sheet;
bridge->link = sheetset;
++ladder;
}
size_t nrOfLaddersPerSheet = sheetset.size();
if (nrOfLaddersPerSheet > kHistogramSize)
nrOfLaddersPerSheet = kHistogramSize;
......@@ -822,17 +824,17 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
stats.laddersPerSheetHistogram[0] += 1;
else if (nrOfLaddersPerSheet > 1)
stats.laddersPerSheetHistogram[nrOfLaddersPerSheet - 1] += 1;
++sheet;
}
for (Bridge& bridge : bridges)
for (Bridge &bridge : bridges)
{
// find out if any of the i and j set members already have
// a bridge assigned, if so, we're assigning bridge 2
uint32_t betai = 0, betaj = 0;
for (uint32_t l : bridge.i)
{
if (inResidues[l].GetBetaPartner(0).residue != nullptr)
......@@ -850,15 +852,15 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
break;
}
}
SecondaryStructureType ss = ssBetabridge;
if (bridge.i.size() > 1)
ss = ssStrand;
if (bridge.type == btParallel)
{
stats.nrOfHBondsInParallelBridges += bridge.i.back() - bridge.i.front() + 2;
std::deque<uint32_t>::iterator j = bridge.j.begin();
for (uint32_t i : bridge.i)
inResidues[i].SetBetaPartner(betai, inResidues[*j++], bridge.ladder, true);
......@@ -898,10 +900,10 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
// --------------------------------------------------------------------
void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, bool inPreferPiHelices = true)
void CalculateAlphaHelices(std::vector<Res> &inResidues, DSSP_Statistics &stats, bool inPreferPiHelices = true)
{
// Helix and Turn
for (HelixType helixType: { HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi })
for (HelixType helixType : {HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi})
{
uint32_t stride = static_cast<uint32_t>(helixType) + 3;
......@@ -915,7 +917,7 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
if (inResidues[j].GetHelixFlag(helixType) == Helix::None)
inResidues[j].SetHelixFlag(helixType, Helix::Middle);
}
if (inResidues[i].GetHelixFlag(helixType) == Helix::End)
inResidues[i].SetHelixFlag(helixType, Helix::StartAndEnd);
else
......@@ -923,8 +925,8 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
}
}
}
for (auto& r : inResidues)
for (auto &r : inResidues)
{
double kappa = r.mM.kappa();
r.SetBend(kappa != 360 and kappa > 70);
......@@ -961,7 +963,7 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
bool empty = true;
for (uint32_t j = i; empty and j <= i + 4; ++j)
empty = inResidues[j].GetSecondaryStructure() == ssLoop or inResidues[j].GetSecondaryStructure() == ssHelix_5 or
(inPreferPiHelices and inResidues[j].GetSecondaryStructure() == ssAlphahelix);
(inPreferPiHelices and inResidues[j].GetSecondaryStructure() == ssAlphahelix);
if (empty)
{
for (uint32_t j = i; j <= i + 4; ++j)
......@@ -969,19 +971,19 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
}
}
}
for (uint32_t i = 1; i + 1 < inResidues.size(); ++i)
{
if (inResidues[i].GetSecondaryStructure() == ssLoop)
{
bool isTurn = false;
for (HelixType helixType: { HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi })
for (HelixType helixType : {HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi})
{
uint32_t stride = 3 + static_cast<uint32_t>(helixType);
for (uint32_t k = 1; k < stride and not isTurn; ++k)
isTurn = (i >= k) and inResidues[i - k].IsHelixStart(helixType);
}
if (isTurn)
inResidues[i].SetSecondaryStructure(ssTurn);
else if (inResidues[i].IsBend())
......@@ -991,7 +993,7 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
std::string asym;
size_t helixLength = 0;
for (auto r: inResidues)
for (auto r : inResidues)
{
if (r.mM.asymID() != asym)
{
......@@ -1014,7 +1016,7 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
// --------------------------------------------------------------------
void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, int stretch_length)
void CalculatePPHelices(std::vector<Res> &inResidues, DSSP_Statistics &stats, int stretch_length)
{
size_t N = inResidues.size();
......@@ -1049,7 +1051,7 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
// auto phi_avg = (phi[i + 0] + phi[i + 1]) / 2;
// auto phi_sq = (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) +
// (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg);
// if (phi_sq >= 200)
// continue;
......@@ -1065,11 +1067,11 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
case Helix::None:
inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Start);
break;
case Helix::End:
inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
break;
default:
break;
}
......@@ -1099,7 +1101,7 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
// auto phi_sq = (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) +
// (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg) +
// (phi[i + 2] - phi_avg) * (phi[i + 2] - phi_avg);
// if (phi_sq >= 300)
// continue;
......@@ -1116,11 +1118,11 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
case Helix::None:
inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Start);
break;
case Helix::End:
inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::StartAndEnd);
break;
default:
break;
}
......@@ -1133,7 +1135,7 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
if (inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
if (inResidues[i + 2].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
inResidues[i + 2].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
......@@ -1150,50 +1152,52 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, in
struct DSSPImpl
{
DSSPImpl(const Structure& s, int min_poly_proline_stretch_length);
const Structure& mStructure;
const std::list<Polymer>& mPolymers;
std::vector<Res> mResidues;
std::vector<std::pair<Res*,Res*>> mSSBonds;
DSSPImpl(const Structure &s, int min_poly_proline_stretch_length);
const Structure &mStructure;
const std::list<Polymer> &mPolymers;
std::vector<Res> mResidues;
std::vector<std::pair<Res *, Res *>> mSSBonds;
int m_min_poly_proline_stretch_length;
auto findRes(const std::string& asymID, int seqID)
auto findRes(const std::string &asymID, int seqID)
{
return std::find_if(mResidues.begin(), mResidues.end(), [&](auto& r) { return r.mM.asymID() == asymID and r.mM.seqID() == seqID; });
return std::find_if(mResidues.begin(), mResidues.end(), [&](auto &r)
{ return r.mM.asymID() == asymID and r.mM.seqID() == seqID; });
}
void calculateSurface();
void calculateSecondaryStructure();
DSSP_Statistics mStats = {};
DSSP_Statistics mStats = {};
};
// --------------------------------------------------------------------
DSSPImpl::DSSPImpl(const Structure& s, int min_poly_proline_stretch_length)
DSSPImpl::DSSPImpl(const Structure &s, int min_poly_proline_stretch_length)
: mStructure(s)
, mPolymers(mStructure.polymers())
, m_min_poly_proline_stretch_length(min_poly_proline_stretch_length)
{
size_t nRes = accumulate(mPolymers.begin(), mPolymers.end(),
0ULL, [](size_t s, auto& p) { return s + p.size(); });
0ULL, [](size_t s, auto &p)
{ return s + p.size(); });
mStats.nrOfChains = static_cast<uint32_t>(mPolymers.size());
mResidues.reserve(nRes);
int resNumber = 0;
for (auto& p: mPolymers)
for (auto &p : mPolymers)
{
ChainBreak brk = ChainBreak::NewChain;
for (auto& m: p)
for (auto &m : p)
{
if (not m.isComplete())
continue;
++resNumber;
if (not mResidues.empty() and
......@@ -1219,15 +1223,15 @@ DSSPImpl::DSSPImpl(const Structure& s, int min_poly_proline_stretch_length)
{
mResidues[i].mNext = &mResidues[i + 1];
mResidues[i + 1].mPrev = &mResidues[i];
mResidues[i + 1].assignHydrogen();
}
}
void DSSPImpl::calculateSecondaryStructure()
{
auto& db = mStructure.getFile().data();
for (auto r: db["struct_conn"].find(cif::Key("conn_type_id") == "disulf"))
auto &db = mStructure.getFile().data();
for (auto r : db["struct_conn"].find(cif::Key("conn_type_id") == "disulf"))
{
std::string asym1, asym2;
int seq1, seq2;
......@@ -1236,7 +1240,7 @@ void DSSPImpl::calculateSecondaryStructure()
auto r1 = findRes(asym1, seq1);
if (r1 == mResidues.end())
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Missing (incomplete?) residue for SS bond when trying to find " << asym1 << '/' << seq1 << std::endl;
continue;
// throw std::runtime_error("Invalid file, missing residue for SS bond");
......@@ -1245,7 +1249,7 @@ void DSSPImpl::calculateSecondaryStructure()
auto r2 = findRes(asym2, seq2);
if (r2 == mResidues.end())
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Missing (incomplete?) residue for SS bond when trying to find " << asym2 << '/' << seq2 << std::endl;
continue;
// throw std::runtime_error("Invalid file, missing residue for SS bond");
......@@ -1261,29 +1265,29 @@ void DSSPImpl::calculateSecondaryStructure()
if (cif::VERBOSE > 1)
{
for (auto& r: mResidues)
for (auto &r : mResidues)
{
auto& m = r.mM;
char helix[5] = { };
for (HelixType helixType: { HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi, HelixType::rh_pp })
auto &m = r.mM;
char helix[5] = {};
for (HelixType helixType : {HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi, HelixType::rh_pp})
{
switch (r.GetHelixFlag(helixType))
{
case Helix::Start: helix[static_cast<int>(helixType)] = '>'; break;
case Helix::Middle: helix[static_cast<int>(helixType)] = helixType == HelixType::rh_pp ? 'P' : '3' + static_cast<char>(helixType); break;
case Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break;
case Helix::End: helix[static_cast<int>(helixType)] = '<'; break;
case Helix::None: helix[static_cast<int>(helixType)] = ' '; break;
case Helix::Start: helix[static_cast<int>(helixType)] = '>'; break;
case Helix::Middle: helix[static_cast<int>(helixType)] = helixType == HelixType::rh_pp ? 'P' : '3' + static_cast<char>(helixType); break;
case Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break;
case Helix::End: helix[static_cast<int>(helixType)] = '<'; break;
case Helix::None: helix[static_cast<int>(helixType)] = ' '; break;
}
}
auto id = m.asymID() + ':' + std::to_string(m.seqID()) + '/' + m.compoundID();
std::cerr << id << std::string(12 - id.length(), ' ')
<< char(r.mSecondaryStructure) << ' '
<< helix
<< std::endl;
<< char(r.mSecondaryStructure) << ' '
<< helix
<< std::endl;
}
}
......@@ -1292,11 +1296,11 @@ void DSSPImpl::calculateSecondaryStructure()
mStats.nrOfIntraChainSSBridges = 0;
uint8_t ssBondNr = 0;
for (const auto& [a, b]: mSSBonds)
for (const auto &[a, b] : mSSBonds)
{
if (a == b)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "In the SS bonds list, the residue " << a->mM << " is bonded to itself" << std::endl;
continue;
}
......@@ -1308,7 +1312,7 @@ void DSSPImpl::calculateSecondaryStructure()
}
mStats.nrOfHBonds = 0;
for (auto& r: mResidues)
for (auto &r : mResidues)
{
auto donor = r.mHBondDonor;
......@@ -1332,7 +1336,7 @@ void DSSPImpl::calculateSurface()
// --------------------------------------------------------------------
const Monomer& DSSP::ResidueInfo::residue() const
const Monomer &DSSP::ResidueInfo::residue() const
{
return mImpl->mM;
}
......@@ -1377,7 +1381,7 @@ double DSSP::ResidueInfo::accessibility() const
return mImpl->mAccessibility;
}
std::tuple<DSSP::ResidueInfo,int,bool> DSSP::ResidueInfo::bridgePartner(int i) const
std::tuple<DSSP::ResidueInfo, int, bool> DSSP::ResidueInfo::bridgePartner(int i) const
{
auto bp = mImpl->GetBetaPartner(i);
......@@ -1391,45 +1395,51 @@ int DSSP::ResidueInfo::sheet() const
return mImpl->GetSheet();
}
std::tuple<DSSP::ResidueInfo,double> DSSP::ResidueInfo::acceptor(int i) const
std::tuple<DSSP::ResidueInfo, double> DSSP::ResidueInfo::acceptor(int i) const
{
auto& a = mImpl->mHBondAcceptor[i];
return { ResidueInfo(a.residue), a.energy };
auto &a = mImpl->mHBondAcceptor[i];
return {ResidueInfo(a.residue), a.energy};
}
std::tuple<DSSP::ResidueInfo,double> DSSP::ResidueInfo::donor(int i) const
std::tuple<DSSP::ResidueInfo, double> DSSP::ResidueInfo::donor(int i) const
{
auto& d = mImpl->mHBondDonor[i];
return { ResidueInfo(d.residue), d.energy };
auto &d = mImpl->mHBondDonor[i];
return {ResidueInfo(d.residue), d.energy};
}
// --------------------------------------------------------------------
DSSP::iterator::iterator(Res* res)
DSSP::iterator::iterator(Res *res)
: mCurrent(res)
{
}
DSSP::iterator::iterator(const iterator& i)
DSSP::iterator::iterator(const iterator &i)
: mCurrent(i.mCurrent)
{
}
DSSP::iterator& DSSP::iterator::operator=(const iterator& i)
DSSP::iterator &DSSP::iterator::operator=(const iterator &i)
{
mCurrent = i.mCurrent;
return *this;
}
DSSP::iterator& DSSP::iterator::operator++()
DSSP::iterator &DSSP::iterator::operator++()
{
++mCurrent.mImpl;
return *this;
}
DSSP::iterator &DSSP::iterator::operator--()
{
--mCurrent.mImpl;
return *this;
}
// --------------------------------------------------------------------
DSSP::DSSP(const Structure& s, int min_poly_proline_stretch, bool calculateSurfaceAccessibility)
DSSP::DSSP(const Structure &s, int min_poly_proline_stretch, bool calculateSurfaceAccessibility)
: mImpl(new DSSPImpl(s, min_poly_proline_stretch))
{
if (calculateSurfaceAccessibility)
......@@ -1455,7 +1465,7 @@ DSSP::iterator DSSP::begin() const
DSSP::iterator DSSP::end() const
{
// careful now, MSVC is picky when it comes to dereferencing iterators that are at the end.
Res* res = nullptr;
Res *res = nullptr;
if (not mImpl->mResidues.empty())
{
res = mImpl->mResidues.data();
......@@ -1465,55 +1475,58 @@ DSSP::iterator DSSP::end() const
return iterator(res);
}
SecondaryStructureType DSSP::operator()(const std::string& inAsymID, int inSeqID) const
SecondaryStructureType DSSP::operator()(const std::string &inAsymID, int inSeqID) const
{
SecondaryStructureType result = ssLoop;
auto i = find_if(mImpl->mResidues.begin(), mImpl->mResidues.end(),
[&](auto& r) { return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
[&](auto &r)
{ return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
if (i != mImpl->mResidues.end())
result = i->mSecondaryStructure;
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << std::endl;
return result;
}
SecondaryStructureType DSSP::operator()(const Monomer& m) const
SecondaryStructureType DSSP::operator()(const Monomer &m) const
{
return operator()(m.asymID(), m.seqID());
}
double DSSP::accessibility(const std::string& inAsymID, int inSeqID) const
double DSSP::accessibility(const std::string &inAsymID, int inSeqID) const
{
SecondaryStructureType result = ssLoop;
auto i = find_if(mImpl->mResidues.begin(), mImpl->mResidues.end(),
[&](auto& r) { return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
[&](auto &r)
{ return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
if (i != mImpl->mResidues.end())
result = i->mSecondaryStructure;
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << std::endl;
return result;
}
double DSSP::accessibility(const Monomer& m) const
double DSSP::accessibility(const Monomer &m) const
{
return accessibility(m.asymID(), m.seqID());
}
bool DSSP::isAlphaHelixEndBeforeStart(const Monomer& m) const
bool DSSP::isAlphaHelixEndBeforeStart(const Monomer &m) const
{
return isAlphaHelixEndBeforeStart(m.asymID(), m.seqID());
}
bool DSSP::isAlphaHelixEndBeforeStart(const std::string& inAsymID, int inSeqID) const
bool DSSP::isAlphaHelixEndBeforeStart(const std::string &inAsymID, int inSeqID) const
{
auto i = find_if(mImpl->mResidues.begin(), mImpl->mResidues.end(),
[&](auto& r) { return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
[&](auto &r)
{ return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
bool result = false;
if (i != mImpl->mResidues.end() and i + 1 != mImpl->mResidues.end())
result = i->GetHelixFlag(HelixType::rh_alpha) == Helix::End and (i + 1)->GetHelixFlag(HelixType::rh_alpha) == Helix::Start;
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
std::cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << std::endl;
return result;
......@@ -1524,4 +1537,4 @@ DSSP_Statistics DSSP::GetStatistics() const
return mImpl->mStats;
}
}
} // namespace mmcif
src/Structure.cpp
/*-
* 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
......@@ -108,8 +108,8 @@ void FileImpl::load_data(const char *data, size_t length)
// And validate, otherwise lots of functionality won't work
// if (mData.getValidator() == nullptr)
mData.loadDictionary("mmcif_pdbx_v50");
if (not mData.isValid())
std::cerr << "Invalid mmCIF file" << (cif::VERBOSE ? "." : " use --verbose option to see errors") << std::endl;
if (not mData.isValid() and cif::VERBOSE >= 0)
std::cerr << "Invalid mmCIF file" << (cif::VERBOSE > 0 ? "." : " use --verbose option to see errors") << std::endl;
}
void FileImpl::load(const std::filesystem::path &path)
......@@ -140,14 +140,14 @@ void FileImpl::load(const std::filesystem::path &path)
{
try
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "unrecognized file extension, trying cif" << std::endl;
mData.load(in);
}
catch (const cif::CifParserError &e)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Not cif, trying plain old PDB" << std::endl;
// pffft...
......@@ -169,7 +169,8 @@ void FileImpl::load(const std::filesystem::path &path)
}
catch (const std::exception &ex)
{
std::cerr << "Error trying to load file " << path << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Error trying to load file " << path << std::endl;
throw;
}
......@@ -179,8 +180,8 @@ void FileImpl::load(const std::filesystem::path &path)
// And validate, otherwise lots of functionality won't work
// if (mData.getValidator() == nullptr)
mData.loadDictionary("mmcif_pdbx_v50");
if (not mData.isValid())
std::cerr << "Invalid mmCIF file" << (cif::VERBOSE ? "." : " use --verbose option to see errors") << std::endl;
if (not mData.isValid() and cif::VERBOSE >= 0)
std::cerr << "Invalid mmCIF file" << (cif::VERBOSE > 0 ? "." : " use --verbose option to see errors") << std::endl;
}
void FileImpl::save(const std::filesystem::path &path)
......@@ -202,494 +203,237 @@ void FileImpl::save(const std::filesystem::path &path)
// --------------------------------------------------------------------
// Atom
struct AtomImpl
Atom::AtomImpl::AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row)
: mDb(db)
, mID(id)
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{
AtomImpl(const AtomImpl &i)
: mDb(i.mDb)
, mID(i.mID)
, mType(i.mType)
, mAtomID(i.mAtomID)
, mCompID(i.mCompID)
, mAsymID(i.mAsymID)
, mSeqID(i.mSeqID)
, mAltID(i.mAltID)
, mLocation(i.mLocation)
, mRefcount(1)
, mRow(i.mRow)
, mCachedRefs(i.mCachedRefs)
, mCompound(i.mCompound)
, mRadius(i.mRadius)
, mSymmetryCopy(i.mSymmetryCopy)
, mClone(true)
// , mRTop(i.mRTop), mD(i.mD)
{
}
AtomImpl(cif::Datablock &db, const std::string &id)
: mDb(db)
, mID(id)
, mRefcount(1)
, mCompound(nullptr)
{
auto &cat = db["atom_site"];
mRow = cat[cif::Key("id") == mID];
prefetch();
}
AtomImpl(cif::Datablock &db, cif::Row &row)
: mDb(db)
, mID(row["id"].as<std::string>())
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{
prefetch();
}
prefetch();
}
AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row)
: mDb(db)
, mID(id)
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{
prefetch();
}
// constructor for a symmetry copy of an atom
Atom::AtomImpl::AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op)
: mDb(impl.mDb)
, mID(impl.mID)
, mType(impl.mType)
AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op)
: mDb(impl.mDb)
, mID(impl.mID)
, mType(impl.mType)
, mAtomID(impl.mAtomID)
, mCompID(impl.mCompID)
, mAsymID(impl.mAsymID)
, mSeqID(impl.mSeqID)
, mAltID(impl.mAltID)
, mLocation(loc)
, mRefcount(1)
, mRow(impl.mRow)
, mCachedRefs(impl.mCachedRefs)
, mCompound(impl.mCompound)
, mRadius(impl.mRadius)
, mSymmetryCopy(true)
, mSymmetryOperator(sym_op)
{
}
, mAtomID(impl.mAtomID)
, mCompID(impl.mCompID)
, mAsymID(impl.mAsymID)
, mSeqID(impl.mSeqID)
, mAltID(impl.mAltID)
, mAuthSeqID(impl.mAuthSeqID)
void prefetch()
{
// Prefetch some data
std::string symbol;
cif::tie(symbol, mAtomID, mCompID, mAsymID, mSeqID, mAltID) =
mRow.get("type_symbol", "label_atom_id", "label_comp_id", "label_asym_id", "label_seq_id", "label_alt_id");
, mLocation(loc)
, mRefcount(1)
, mRow(impl.mRow)
, mCachedRefs(impl.mCachedRefs)
, mCompound(impl.mCompound)
, mSymmetryCopy(true)
, mSymmetryOperator(sym_op)
{
}
if (symbol != "X")
mType = AtomTypeTraits(symbol).type();
void Atom::AtomImpl::prefetch()
{
// Prefetch some data
std::string symbol;
cif::tie(symbol, mAtomID, mCompID, mAsymID, mSeqID, mAltID, mAuthSeqID) =
mRow.get("type_symbol", "label_atom_id", "label_comp_id", "label_asym_id", "label_seq_id", "label_alt_id", "auth_seq_id");
float x, y, z;
cif::tie(x, y, z) = mRow.get("Cartn_x", "Cartn_y", "Cartn_z");
if (symbol != "X")
mType = AtomTypeTraits(symbol).type();
mLocation = Point(x, y, z);
float x, y, z;
cif::tie(x, y, z) = mRow.get("Cartn_x", "Cartn_y", "Cartn_z");
std::string compID;
cif::tie(compID) = mRow.get("label_comp_id");
mLocation = Point(x, y, z);
}
// mCompound = CompoundFactory::instance().create(compID);
}
int Atom::AtomImpl::compare(const AtomImpl &b) const
{
int d = mAsymID.compare(b.mAsymID);
if (d == 0)
d = mSeqID - b.mSeqID;
if (d == 0)
d = mAtomID.compare(b.mAtomID);
if (d == 0)
d = mAuthSeqID.compare(b.mAuthSeqID);
void reference()
{
++mRefcount;
}
return d;
}
void release()
{
if (--mRefcount <= 0)
delete this;
}
bool Atom::AtomImpl::getAnisoU(float anisou[6]) const
{
bool result = false;
bool getAnisoU(float anisou[6]) const
auto cat = mDb.get("atom_site_anisotrop");
if (cat)
{
bool result = false;
auto cat = mDb.get("atom_site_anisotrop");
if (cat)
try
{
try
{
auto r = cat->find1(cif::Key("id") == mID);
cif::tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
result = true;
}
catch(const std::exception& e)
{
}
auto r = cat->find1(cif::Key("id") == mID);
cif::tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
result = true;
}
return result;
}
void moveTo(const Point &p)
{
assert(not mSymmetryCopy);
if (mSymmetryCopy)
throw std::runtime_error("Moving symmetry copy");
if (not mClone)
catch (const std::exception &e)
{
property("Cartn_x", std::to_string(p.getX()));
property("Cartn_y", std::to_string(p.getY()));
property("Cartn_z", std::to_string(p.getZ()));
}
// boost::format kPosFmt("%.3f");
//
// mRow["Cartn_x"] = (kPosFmt % p.getX()).str();
// mRow["Cartn_y"] = (kPosFmt % p.getY()).str();
// mRow["Cartn_z"] = (kPosFmt % p.getZ()).str();
mLocation = p;
}
const Compound &comp() const
{
if (mCompound == nullptr)
{
std::string compID;
cif::tie(compID) = mRow.get("label_comp_id");
mCompound = CompoundFactory::instance().create(compID);
if (cif::VERBOSE and mCompound == nullptr)
std::cerr << "Compound not found: '" << compID << '\'' << std::endl;
}
if (mCompound == nullptr)
throw std::runtime_error("no compound");
return *mCompound;
}
return result;
}
bool isWater() const
{
// mCompound may still be null here, and besides, this check is not that exciting anyway
return mCompID == "HOH" or mCompID == "H2O" or mCompID == "WAT";
}
void Atom::AtomImpl::moveTo(const Point &p)
{
assert(not mSymmetryCopy);
if (mSymmetryCopy)
throw std::runtime_error("Moving symmetry copy");
float radius() const
if (not mClone)
{
return mRadius;
mRow.assign("Cartn_x", std::to_string(p.getX()), true, false);
mRow.assign("Cartn_y", std::to_string(p.getY()), true, false);
mRow.assign("Cartn_z", std::to_string(p.getZ()), true, false);
}
const std::string property(const std::string_view name) const
{
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag == name)
return ref.as<std::string>();
}
mCachedRefs.emplace_back(name, mRow[name]);
return std::get<1>(mCachedRefs.back()).as<std::string>();
}
mLocation = p;
}
void property(const std::string_view name, const std::string &value)
const Compound &Atom::AtomImpl::comp() const
{
if (mCompound == nullptr)
{
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag != name)
continue;
ref = value;
return;
}
std::string compID;
cif::tie(compID) = mRow.get("label_comp_id");
mCachedRefs.emplace_back(name, mRow[name]);
std::get<1>(mCachedRefs.back()) = value;
}
mCompound = CompoundFactory::instance().create(compID);
int compare(const AtomImpl &b) const
{
int d = mAsymID.compare(b.mAsymID);
if (d == 0)
d = mSeqID - b.mSeqID;
if (d == 0)
d = mAtomID.compare(b.mAtomID);
return d;
if (cif::VERBOSE > 0 and mCompound == nullptr)
std::cerr << "Compound not found: '" << compID << '\'' << std::endl;
}
void swapAtomLabels(AtomImpl &b)
{
std::swap(mAtomID, b.mAtomID);
}
const cif::Datablock &mDb;
std::string mID;
AtomType mType;
std::string mAtomID;
std::string mCompID;
std::string mAsymID;
int mSeqID;
std::string mAltID;
Point mLocation;
int mRefcount;
cif::Row mRow;
if (mCompound == nullptr)
throw std::runtime_error("no compound");
mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
mutable const Compound *mCompound = nullptr;
float mRadius = std::nanf("4");
bool mSymmetryCopy = false;
bool mClone = false;
std::string mSymmetryOperator = "1_555";
// clipper::RTop_orth mRTop;
// Point mD;
// int32_t mRTix;
};
//Atom::Atom(const File& f, const std::string& id)
// : mImpl(new AtomImpl(f, id))
//{
//}
//
Atom::Atom()
: mImpl_(nullptr)
{
return *mCompound;
}
Atom::Atom(AtomImpl *impl)
: mImpl_(impl)
const std::string Atom::AtomImpl::get_property(const std::string_view name) const
{
}
Atom::Atom(cif::Datablock &db, cif::Row &row)
: mImpl_(new AtomImpl(db, row))
{
}
AtomImpl *Atom::impl()
{
if (mImpl_ == nullptr)
throw std::runtime_error("atom is not set");
return mImpl_;
}
const AtomImpl *Atom::impl() const
{
if (mImpl_ == nullptr)
throw std::runtime_error("atom is not set");
return mImpl_;
}
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag == name)
return ref.as<std::string>();
}
Atom Atom::clone() const
{
return Atom(mImpl_ ? new AtomImpl(*mImpl_) : nullptr);
mCachedRefs.emplace_back(name, mRow[name]);
return std::get<1>(mCachedRefs.back()).as<std::string>();
}
Atom::Atom(const Atom &rhs, const Point &loc, const std::string &sym_op)
: mImpl_(new AtomImpl(*rhs.mImpl_, loc, sym_op))
void Atom::AtomImpl::set_property(const std::string_view name, const std::string &value)
{
}
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag != name)
continue;
Atom::Atom(const Atom &rhs)
: mImpl_(rhs.mImpl_)
{
if (mImpl_)
mImpl_->reference();
}
ref = value;
return;
}
Atom::~Atom()
{
if (mImpl_)
mImpl_->release();
mCachedRefs.emplace_back(name, mRow[name]);
std::get<1>(mCachedRefs.back()) = value;
}
Atom &Atom::operator=(const Atom &rhs)
Atom::Atom(cif::Datablock &db, cif::Row &row)
: Atom(std::make_shared<AtomImpl>(db, row["id"].as<std::string>(), row))
{
if (this != &rhs)
{
if (mImpl_)
mImpl_->release();
mImpl_ = rhs.mImpl_;
if (mImpl_)
mImpl_->reference();
}
return *this;
}
const cif::Row Atom::getRow() const
Atom::Atom(const Atom &rhs, const Point &loc, const std::string &sym_op)
: Atom(std::make_shared<AtomImpl>(*rhs.mImpl, loc, sym_op))
{
return mImpl_->mRow;
}
const cif::Row Atom::getRowAniso() const
{
auto &db = mImpl_->mDb;
auto &db = mImpl->mDb;
auto cat = db.get("atom_site_anisotrop");
if (not cat)
return {};
else
return cat->find1(cif::Key("id") == mImpl_->mID);
}
template <>
std::string Atom::property<std::string>(const std::string_view name) const
{
return impl()->property(name);
}
template <>
int Atom::property<int>(const std::string_view name) const
{
auto v = impl()->property(name);
return v.empty() ? 0 : stoi(v);
}
template <>
float Atom::property<float>(const std::string_view name) const
{
return stof(impl()->property(name));
}
void Atom::property(const std::string_view name, const std::string &value)
{
impl()->property(name, value);
}
const std::string &Atom::id() const
{
return impl()->mID;
}
AtomType Atom::type() const
{
return impl()->mType;
}
int Atom::charge() const
{
return property<int>("pdbx_formal_charge");
return cat->find1(cif::Key("id") == mImpl->mID);
}
float Atom::uIso() const
{
float result;
if (not property<std::string>("U_iso_or_equiv").empty())
result = property<float>("U_iso_or_equiv");
else if (not property<std::string>("B_iso_or_equiv").empty())
result = property<float>("B_iso_or_equiv") / static_cast<float>(8 * kPI * kPI);
if (not get_property<std::string>("U_iso_or_equiv").empty())
result = get_property<float>("U_iso_or_equiv");
else if (not get_property<std::string>("B_iso_or_equiv").empty())
result = get_property<float>("B_iso_or_equiv") / static_cast<float>(8 * kPI * kPI);
else
throw std::runtime_error("Missing B_iso or U_iso");
return result;
}
bool Atom::getAnisoU(float anisou[6]) const
{
return impl()->getAnisoU(anisou);
}
float Atom::occupancy() const
std::string Atom::labelID() const
{
return property<float>("occupancy");
return mImpl->mCompID + '_' + mImpl->mAsymID + '_' + std::to_string(mImpl->mSeqID) + ':' + mImpl->mAtomID;
}
std::string Atom::labelAtomID() const
std::string Atom::pdbID() const
{
return impl()->mAtomID;
return get_property<std::string>("auth_comp_id") + '_' +
get_property<std::string>("auth_asym_id") + '_' +
get_property<std::string>("auth_seq_id") +
get_property<std::string>("pdbx_PDB_ins_code");
}
std::string Atom::labelCompID() const
int Atom::charge() const
{
return impl()->mCompID;
return get_property<int>("pdbx_formal_charge");
}
std::string Atom::labelAsymID() const
float Atom::occupancy() const
{
return impl()->mAsymID;
return get_property<float>("occupancy");
}
std::string Atom::labelEntityID() const
{
return property<std::string>("label_entity_id");
}
std::string Atom::labelAltID() const
{
return impl()->mAltID;
}
bool Atom::isAlternate() const
{
return not impl()->mAltID.empty();
}
int Atom::labelSeqID() const
{
return impl()->mSeqID;
}
std::string Atom::authAsymID() const
{
return property<std::string>("auth_asym_id");
return get_property<std::string>("label_entity_id");
}
std::string Atom::authAtomID() const
{
return property<std::string>("auth_atom_id");
}
std::string Atom::pdbxAuthAltID() const
{
return property<std::string>("pdbx_auth_alt_id");
}
std::string Atom::pdbxAuthInsCode() const
{
return property<std::string>("pdbx_PDB_ins_code");
return get_property<std::string>("auth_atom_id");
}
std::string Atom::authCompID() const
{
return property<std::string>("auth_comp_id");
return get_property<std::string>("auth_comp_id");
}
std::string Atom::authSeqID() const
{
return property<std::string>("auth_seq_id");
}
std::string Atom::labelID() const
{
return property<std::string>("label_comp_id") + '_' + impl()->mAsymID + '_' + std::to_string(impl()->mSeqID) + ':' + impl()->mAtomID;
}
std::string Atom::pdbID() const
std::string Atom::authAsymID() const
{
return property<std::string>("auth_comp_id") + '_' +
property<std::string>("auth_asym_id") + '_' +
property<std::string>("auth_seq_id") +
property<std::string>("pdbx_PDB_ins_code");
return get_property<std::string>("auth_asym_id");
}
Point Atom::location() const
std::string Atom::pdbxAuthInsCode() const
{
return impl()->mLocation;
return get_property<std::string>("pdbx_PDB_ins_code");
}
void Atom::location(Point p)
std::string Atom::pdbxAuthAltID() const
{
impl()->moveTo(p);
return get_property<std::string>("pdbx_auth_alt_id");
}
void Atom::translate(Point t)
......@@ -706,80 +450,32 @@ void Atom::rotate(Quaternion q)
location(loc);
}
// Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
// {
// return Atom(new AtomImpl(*impl(), d, rt));
// }
// bool Atom::isSymmetryCopy() const
// {
// return impl()->mSymmetryCopy;
// }
// std::string Atom::symmetry() const
// {
// return clipper::Symop(impl()->mRTop).format() + "\n" + impl()->mRTop.format();
// }
bool Atom::isSymmetryCopy() const
void Atom::translateAndRotate(Point t, Quaternion q)
{
return mImpl_->mSymmetryCopy;
}
std::string Atom::symmetry() const
{
return mImpl_->mSymmetryOperator;
}
// const clipper::RTop_orth& Atom::symop() const
// {
// return impl()->mRTop;
// }
const Compound &Atom::comp() const
{
return impl()->comp();
auto loc = location();
loc += t;
loc.rotate(q);
location(loc);
}
bool Atom::isWater() const
void Atom::translateRotateAndTranslate(Point t1, Quaternion q, Point t2)
{
return impl()->isWater();
auto loc = location();
loc += t1;
loc.rotate(q);
loc += t2;
location(loc);
}
bool Atom::operator==(const Atom &rhs) const
{
return impl() == rhs.impl() or
(&impl()->mDb == &rhs.impl()->mDb and impl()->mID == rhs.impl()->mID);
}
// clipper::Atom Atom::toClipper() const
// {
// return impl()->toClipper();
// }
// void Atom::calculateRadius(float resHigh, float resLow, float perc)
// {
// AtomShape shape(*this, resHigh, resLow, false);
// impl()->mRadius = shape.radius();
// // verbose
// if (cif::VERBOSE > 1)
// cout << "Calculated radius for " << AtomTypeTraits(impl()->mType).name() << " with charge " << charge() << " is " << impl()->mRadius << std::endl;
// }
float Atom::radius() const
{
return impl()->mRadius;
}
int Atom::compare(const Atom &b) const
{
return impl() == b.impl() ? 0 : impl()->compare(*b.impl());
return mImpl == rhs.mImpl or
(&mImpl->mDb == &rhs.mImpl->mDb and mImpl->mID == rhs.mImpl->mID);
}
void Atom::setID(int id)
{
impl()->mID = std::to_string(id);
mImpl->mID = std::to_string(id);
}
std::ostream &operator<<(std::ostream &os, const Atom &atom)
......@@ -797,66 +493,6 @@ std::ostream &operator<<(std::ostream &os, const Atom &atom)
// --------------------------------------------------------------------
// residue
// First constructor used to be for waters only, but now accepts sugars as well.
Residue::Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, const std::string &authSeqID)
: mStructure(&structure)
, mCompoundID(compoundID)
, mAsymID(asymID)
, mAuthSeqID(authSeqID)
{
for (auto &a : mStructure->atoms())
{
if (a.labelAsymID() != mAsymID or
a.labelCompID() != mCompoundID)
continue;
if (compoundID == "HOH")
{
if (not mAuthSeqID.empty() and a.authSeqID() != mAuthSeqID)
continue;
}
else
{
if (mSeqID > 0 and a.labelSeqID() != mSeqID)
continue;
}
mAtoms.push_back(a);
}
assert(not mAtoms.empty());
}
Residue::Residue(const Structure &structure, const std::string &compoundID, const std::string &asymID)
: Residue(structure, compoundID, asymID, 0, {})
{
}
Residue::Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, int seqID, const std::string &authSeqID)
: mStructure(&structure)
, mCompoundID(compoundID)
, mAsymID(asymID)
, mSeqID(seqID)
, mAuthSeqID(authSeqID)
{
assert(mCompoundID != "HOH");
for (auto &a : mStructure->atoms())
{
if (mSeqID > 0 and a.labelSeqID() != mSeqID)
continue;
if (a.labelAsymID() != mAsymID or
a.labelCompID() != mCompoundID)
continue;
mAtoms.push_back(a);
}
}
Residue::Residue(Residue &&rhs)
: mStructure(rhs.mStructure)
, mCompoundID(std::move(rhs.mCompoundID))
......@@ -865,13 +501,13 @@ Residue::Residue(Residue &&rhs)
, mAuthSeqID(rhs.mAuthSeqID)
, mAtoms(std::move(rhs.mAtoms))
{
//std::cerr << "move constructor residue" << std::endl;
// std::cerr << "move constructor residue" << std::endl;
rhs.mStructure = nullptr;
}
Residue &Residue::operator=(Residue &&rhs)
{
//std::cerr << "move assignment residue" << std::endl;
// std::cerr << "move assignment residue" << std::endl;
mStructure = rhs.mStructure;
rhs.mStructure = nullptr;
mCompoundID = std::move(rhs.mCompoundID);
......@@ -885,7 +521,7 @@ Residue &Residue::operator=(Residue &&rhs)
Residue::~Residue()
{
//std::cerr << "~Residue" << std::endl;
// std::cerr << "~Residue" << std::endl;
}
std::string Residue::entityID() const
......@@ -999,7 +635,7 @@ AtomView Residue::unique_atoms() const
firstAlt = alt;
else if (alt != firstAlt)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "skipping alternate atom " << atom << std::endl;
continue;
}
......@@ -1145,11 +781,6 @@ std::ostream &operator<<(std::ostream &os, const Residue &res)
// --------------------------------------------------------------------
// monomer
//Monomer::Monomer(Monomer&& rhs)
// : Residue(std::move(rhs)), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex)
//{
//}
Monomer::Monomer(const Polymer &polymer, size_t index, int seqID, const std::string &authSeqID, const std::string &compoundID)
: Residue(*polymer.structure(), compoundID, polymer.asymID(), seqID, authSeqID)
, mPolymer(&polymer)
......@@ -1162,23 +793,11 @@ Monomer::Monomer(Monomer &&rhs)
, mPolymer(rhs.mPolymer)
, mIndex(rhs.mIndex)
{
std::cerr << "move constructor monomer" << std::endl;
// mStructure = rhs.mStructure; rhs.mStructure = nullptr;
// mCompoundID = std::move(rhs.mCompoundID);
// mAsymID = std::move(rhs.mAsymID);
// mSeqID = rhs.mSeqID;
// mAtoms = std::move(rhs.mAtoms);
//
// mPolymer = rhs.mPolymer; rhs.mPolymer = nullptr;
// mIndex = rhs.mIndex;
rhs.mPolymer = nullptr;
}
Monomer &Monomer::operator=(Monomer &&rhs)
{
std::cerr << "move assignment monomer" << std::endl;
Residue::operator=(std::move(rhs));
mPolymer = rhs.mPolymer;
rhs.mPolymer = nullptr;
......@@ -1222,7 +841,7 @@ float Monomer::phi() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << ex.what() << std::endl;
}
......@@ -1244,7 +863,7 @@ float Monomer::psi() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << ex.what() << std::endl;
}
......@@ -1268,7 +887,7 @@ float Monomer::alpha() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << ex.what() << std::endl;
}
......@@ -1296,7 +915,7 @@ float Monomer::kappa() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "When trying to calculate kappa for " << asymID() << ':' << seqID() << ": "
<< ex.what() << std::endl;
}
......@@ -1319,7 +938,7 @@ float Monomer::tco() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "When trying to calculate tco for " << asymID() << ':' << seqID() << ": "
<< ex.what() << std::endl;
}
......@@ -1338,7 +957,7 @@ float Monomer::omega() const
}
catch (const std::exception &ex)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "When trying to calculate omega for " << asymID() << ':' << seqID() << ": "
<< ex.what() << std::endl;
}
......@@ -1409,7 +1028,7 @@ float Monomer::chi(size_t nr) const
}
catch (const std::exception &e)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << e.what() << std::endl;
result = 0;
}
......@@ -1550,94 +1169,6 @@ bool Monomer::isCis(const mmcif::Monomer &a, const mmcif::Monomer &b)
// --------------------------------------------------------------------
// polymer
//
//Polymer::iterator::iterator(const Polymer& p, uint32_t index)
// : mPolymer(&p), mIndex(index), mCurrent(p, index)
//{
// auto& polySeq = mPolymer->mPolySeq;
//
// if (index < polySeq.size())
// {
// int seqID;
// std::string asymID, monID;
// cif::tie(asymID, seqID, monID) =
// polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
//
// mCurrent = Monomer(*mPolymer, index, seqID, monID, "");
// }
//}
//
//Monomer Polymer::operator[](size_t index) const
//{
// if (index >= mPolySeq.size())
// throw out_of_range("Invalid index for residue in polymer");
//
// std::string compoundID;
// int seqID;
//
// auto r = mPolySeq[index];
//
// cif::tie(seqID, compoundID) =
// r.get("seq_id", "mon_id");
//
// return Monomer(const_cast<Polymer&>(*this), index, seqID, compoundID, "");
//}
//
//Polymer::iterator::iterator(const iterator& rhs)
// : mPolymer(rhs.mPolymer), mIndex(rhs.mIndex), mCurrent(rhs.mCurrent)
//{
//}
//
//Polymer::iterator& Polymer::iterator::operator++()
//{
// auto& polySeq = mPolymer->mPolySeq;
//
// if (mIndex < polySeq.size())
// ++mIndex;
//
// if (mIndex < polySeq.size())
// {
// int seqID;
// std::string asymID, monID;
// cif::tie(asymID, seqID, monID) =
// polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
//
// mCurrent = Monomer(*mPolymer, mIndex, seqID, monID, "");
// }
//
// return *this;
//}
//Polymer::Polymer(const Structure& s, const std::string& asymID)
// : mStructure(const_cast<Structure*>(&s)), mAsymID(asymID)
// , mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID))
//{
// mEntityID = mPolySeq.front()["entity_id"].as<std::string>();
//
//#if DEBUG
// for (auto r: mPolySeq)
// assert(r["entity_id"] == mEntityID);
//#endif
//
//}
//Polymer::Polymer(Polymer&& rhs)
// : std::vector<Monomer>(std::move(rhs))
// , mStructure(rhs.mStructure)
// , mEntityID(std::move(rhs.mEntityID)), mAsymID(std::move(rhs.mAsymID)), mPolySeq(std::move(rhs.mPolySeq))
//{
// rhs.mStructure = nullptr;
//}
//
//Polymer& Polymer::operator=(Polymer&& rhs)
//{
// std::vector<Monomer>::operator=(std::move(rhs));
// mStructure = rhs.mStructure; rhs.mStructure = nullptr;
// mEntityID = std::move(rhs.mEntityID);
// mAsymID = std::move(rhs.mAsymID);
// mPolySeq = std::move(rhs.mPolySeq);
// return *this;
//}
Polymer::Polymer(const Structure &s, const std::string &entityID, const std::string &asymID)
: mStructure(const_cast<Structure *>(&s))
......@@ -1663,7 +1194,7 @@ Polymer::Polymer(const Structure &s, const std::string &entityID, const std::str
ix[seqID] = index;
emplace_back(*this, index, seqID, authSeqID, compoundID);
}
else if (cif::VERBOSE)
else if (cif::VERBOSE > 0)
{
Monomer m{*this, index, seqID, authSeqID, compoundID};
std::cerr << "Dropping alternate residue " << m << std::endl;
......@@ -1745,7 +1276,7 @@ File::~File()
delete mImpl;
}
cif::Datablock& File::createDatablock(const std::string_view name)
cif::Datablock &File::createDatablock(const std::string_view name)
{
auto db = new cif::Datablock(name);
......@@ -1808,7 +1339,7 @@ Structure::Structure(File &f, size_t modelNr, StructureOpenOptions options)
cif::tie(model_nr) = atomCat.front().get("pdbx_PDB_model_num");
if (model_nr and *model_nr != mModelNr)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "No atoms loaded for model 1, trying model " << *model_nr << std::endl;
mModelNr = *model_nr;
loadAtomsForModel(options);
......@@ -1816,7 +1347,10 @@ Structure::Structure(File &f, size_t modelNr, StructureOpenOptions options)
}
if (mAtoms.empty())
std::cerr << "Warning: no atoms loaded" << std::endl;
{
if (cif::VERBOSE >= 0)
std::cerr << "Warning: no atoms loaded" << std::endl;
}
else
loadData();
}
......@@ -1839,11 +1373,10 @@ void Structure::loadAtomsForModel(StructureOpenOptions options)
if ((options bitand StructureOpenOptions::SkipHydrogen) and type_symbol == "H")
continue;
mAtoms.emplace_back(new AtomImpl(*db, id, a));
mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(*db, id, a));
}
}
Structure::Structure(const Structure &s)
: mFile(s.mFile)
, mModelNr(s.mModelNr)
......@@ -1884,7 +1417,7 @@ void Structure::loadData()
r.get("asym_id", "mon_id", "pdb_seq_num");
if (monID == "HOH")
mNonPolymers.emplace_back(*this, monID, asymID, pdbSeqNum);
mNonPolymers.emplace_back(*this, monID, asymID, 0, pdbSeqNum);
else if (mNonPolymers.empty() or mNonPolymers.back().asymID() != asymID)
mNonPolymers.emplace_back(*this, monID, asymID);
}
......@@ -1897,7 +1430,43 @@ void Structure::loadData()
cif::tie(asymID, monID, num) =
r.get("asym_id", "mon_id", "num");
mBranchResidues.emplace_back(*this, monID, asymID, num);
mBranchResidues.emplace_back(*this, monID, asymID, 0, num);
}
// place atoms in residues
using key_type = std::tuple<std::string, std::string, int>;
using map_type = std::map<key_type, Residue *>;
map_type resMap;
for (auto &poly : mPolymers)
{
for (auto &res : poly)
resMap[{res.asymID(), res.compoundID(), res.seqID()}] = &res;
}
for (auto &res : mNonPolymers)
resMap[{res.asymID(), res.compoundID(), (res.isWater() ? std::stoi(res.mAuthSeqID) : res.seqID())}] = &res;
for (auto &res : mBranchResidues)
resMap[{res.asymID(), res.compoundID(), res.seqID()}] = &res;
for (auto &atom : mAtoms)
{
key_type k(atom.labelAsymID(), atom.labelCompID(), atom.isWater() ? std::stoi(atom.authSeqID()) : atom.labelSeqID());
auto ri = resMap.find(k);
if (ri == resMap.end())
{
if (cif::VERBOSE > 0)
std::cerr << "Missing residue for atom " << atom << std::endl;
assert(false);
continue;
}
ri->second->addAtom(atom);
}
}
......@@ -1948,7 +1517,7 @@ AtomView Structure::waters() const
for (auto &a : mAtoms)
{
if (a.property<std::string>("label_entity_id") == waterEntityID)
if (a.get_property<std::string>("label_entity_id") == waterEntityID)
result.push_back(a);
}
......@@ -2021,7 +1590,7 @@ Atom Structure::getAtomByPositionAndType(Point p, std::string_view type, std::st
if (a.labelCompID() != res_type)
continue;
if (a.labelAtomID() != type)
continue;
......@@ -2077,7 +1646,7 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin
Residue &Structure::getResidue(const std::string &asymID, const std::string &compID, int seqID)
{
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID, compID, seqID));
return const_cast<Residue &>(const_cast<Structure const &>(*this).getResidue(asymID, compID, seqID));
}
const Residue &Structure::getResidue(const std::string &asymID) const
......@@ -2095,7 +1664,7 @@ const Residue &Structure::getResidue(const std::string &asymID) const
Residue &Structure::getResidue(const std::string &asymID)
{
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID));
return const_cast<Residue &>(const_cast<Structure const &>(*this).getResidue(asymID));
}
Residue &Structure::getResidue(const mmcif::Atom &atom)
......@@ -2291,8 +1860,7 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
auto r = chemComp.find(cif::Key("id") == compoundID);
if (r.empty())
{
chemComp.emplace({
{"id", compoundID},
chemComp.emplace({{"id", compoundID},
{"name", compound->name()},
{"formula", compound->formula()},
{"formula_weight", compound->formulaWeight()},
......@@ -2308,19 +1876,17 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
{
entity_id = pdbxEntityNonpoly.find1<std::string>("comp_id"_key == compoundID, "entity_id");
}
catch(const std::exception& ex)
catch (const std::exception &ex)
{
auto &entity = db["entity"];
entity_id = entity.getUniqueID("");
entity.emplace({
{"id", entity_id},
entity.emplace({{"id", entity_id},
{"type", "non-polymer"},
{"pdbx_description", compound->name()},
{"formula_weight", compound->formulaWeight()}});
pdbxEntityNonpoly.emplace({
{"entity_id", entity_id},
pdbxEntityNonpoly.emplace({{"entity_id", entity_id},
{"name", compound->name()},
{"comp_id", compoundID}});
}
......@@ -2329,65 +1895,6 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
return entity_id;
}
// // --------------------------------------------------------------------
// Structure::residue_iterator::residue_iterator(const Structure* s, poly_iterator polyIter, size_t polyResIndex, size_t nonPolyIndex)
// : mStructure(*s), mPolyIter(polyIter), mPolyResIndex(polyResIndex), mNonPolyIndex(nonPolyIndex)
// {
// while (mPolyIter != mStructure.mPolymers.end() and mPolyResIndex == mPolyIter->size())
// ++mPolyIter;
// }
// auto Structure::residue_iterator::operator*() -> reference
// {
// if (mPolyIter != mStructure.mPolymers.end())
// return (*mPolyIter)[mPolyResIndex];
// else
// return mStructure.mNonPolymers[mNonPolyIndex];
// }
// auto Structure::residue_iterator::operator->() -> pointer
// {
// if (mPolyIter != mStructure.mPolymers.end())
// return &(*mPolyIter)[mPolyResIndex];
// else
// return &mStructure.mNonPolymers[mNonPolyIndex];
// }
// Structure::residue_iterator& Structure::residue_iterator::operator++()
// {
// if (mPolyIter != mStructure.mPolymers.end())
// {
// ++mPolyResIndex;
// if (mPolyResIndex >= mPolyIter->size())
// {
// ++mPolyIter;
// mPolyResIndex = 0;
// }
// }
// else
// ++mNonPolyIndex;
// return *this;
// }
// Structure::residue_iterator Structure::residue_iterator::operator++(int)
// {
// auto result = *this;
// operator++();
// return result;
// }
// bool Structure::residue_iterator::operator==(const Structure::residue_iterator& rhs) const
// {
// return mPolyIter == rhs.mPolyIter and mPolyResIndex == rhs.mPolyResIndex and mNonPolyIndex == rhs.mNonPolyIndex;
// }
// bool Structure::residue_iterator::operator!=(const Structure::residue_iterator& rhs) const
// {
// return mPolyIter != rhs.mPolyIter or mPolyResIndex != rhs.mPolyResIndex or mNonPolyIndex != rhs.mNonPolyIndex;
// }
// --------------------------------------------------------------------
void Structure::removeAtom(Atom &a)
......@@ -2434,7 +1941,7 @@ void Structure::swapAtoms(Atom &a1, Atom &a2)
auto l2 = r2["label_atom_id"];
l1.swap(l2);
a1.impl()->swapAtomLabels(*a2.impl());
std::swap(a1.mImpl->mAtomID, a2.mImpl->mAtomID);
auto l3 = r1["auth_atom_id"];
auto l4 = r2["auth_atom_id"];
......@@ -2524,7 +2031,8 @@ void Structure::changeResidue(Residue &res, const std::string &newCompound,
{ return a.labelAtomID() == a1; });
if (i == atoms.end())
{
std::cerr << "Missing atom for atom ID " << a1 << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "Missing atom for atom ID " << a1 << std::endl;
continue;
}
......@@ -2564,50 +2072,47 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
auto &struct_asym = db["struct_asym"];
std::string asym_id = struct_asym.getUniqueID();
struct_asym.emplace({
{ "id", asym_id },
{ "pdbx_blank_PDB_chainid_flag", "N" },
{ "pdbx_modified", "N" },
{ "entity_id", entity_id },
{ "details", "?" }
});
struct_asym.emplace({{"id", asym_id},
{"pdbx_blank_PDB_chainid_flag", "N"},
{"pdbx_modified", "N"},
{"entity_id", entity_id},
{"details", "?"}});
std::string comp_id = db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id");
auto &atom_site = db["atom_site"];
for (auto& atom: atoms)
auto &res = mNonPolymers.emplace_back(*this, comp_id, asym_id);
for (auto &atom : atoms)
{
auto atom_id = atom_site.getUniqueID("");
auto &&[row, inserted ] = atom_site.emplace({
{ "group_PDB", atom.property<std::string>("group_PDB") },
{ "id", atom_id },
{ "type_symbol", atom.property<std::string>("type_symbol") },
{ "label_atom_id", atom.property<std::string>("label_atom_id") },
{ "label_alt_id", atom.property<std::string>("label_alt_id") },
{ "label_comp_id", comp_id },
{ "label_asym_id", asym_id },
{ "label_entity_id", entity_id },
{ "label_seq_id", "." },
{ "pdbx_PDB_ins_code", "" },
{ "Cartn_x", atom.property<std::string>("Cartn_x") },
{ "Cartn_y", atom.property<std::string>("Cartn_y") },
{ "Cartn_z", atom.property<std::string>("Cartn_z") },
{ "occupancy", atom.property<std::string>("occupancy") },
{ "B_iso_or_equiv", atom.property<std::string>("B_iso_or_equiv") },
{ "pdbx_formal_charge", atom.property<std::string>("pdbx_formal_charge") },
{ "auth_seq_id", "" },
{ "auth_comp_id", comp_id },
{ "auth_asym_id", asym_id },
{ "auth_atom_id", atom.property<std::string>("label_atom_id") },
{ "pdbx_PDB_model_num", 1 }
});
mAtoms.emplace_back(new AtomImpl(db, atom_id, row));
}
mNonPolymers.emplace_back(*this, comp_id, asym_id);
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")},
{"label_alt_id", atom.get_property<std::string>("label_alt_id")},
{"label_comp_id", comp_id},
{"label_asym_id", asym_id},
{"label_entity_id", entity_id},
{"label_seq_id", "."},
{"pdbx_PDB_ins_code", ""},
{"Cartn_x", atom.get_property<std::string>("Cartn_x")},
{"Cartn_y", atom.get_property<std::string>("Cartn_y")},
{"Cartn_z", atom.get_property<std::string>("Cartn_z")},
{"occupancy", atom.get_property<std::string>("occupancy")},
{"B_iso_or_equiv", atom.get_property<std::string>("B_iso_or_equiv")},
{"pdbx_formal_charge", atom.get_property<std::string>("pdbx_formal_charge")},
{"auth_seq_id", ""},
{"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}});
auto &newAtom = mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
res.addAtom(newAtom);
}
return asym_id;
}
......@@ -2685,7 +2190,7 @@ void Structure::cleanupEmptyCategories()
{
// is this correct?
std::set<std::string> asym_ids;
for (const auto &[ asym_id ] : db["pdbx_branch_scheme"].find<std::string>("entity_id"_key == id, "asym_id"))
for (const auto &[asym_id] : db["pdbx_branch_scheme"].find<std::string>("entity_id"_key == id, "asym_id"))
asym_ids.insert(asym_id);
count = asym_ids.size();
}
......@@ -2696,14 +2201,26 @@ void Structure::cleanupEmptyCategories()
void Structure::translate(Point t)
{
for (auto& a: mAtoms)
for (auto &a : mAtoms)
a.translate(t);
}
void Structure::rotate(Quaternion q)
{
for (auto& a: mAtoms)
for (auto &a : mAtoms)
a.rotate(q);
}
void Structure::translateAndRotate(Point t, Quaternion q)
{
for (auto &a : mAtoms)
a.translateAndRotate(t, q);
}
void Structure::translateRotateAndTranslate(Point t1, Quaternion q, Point t2)
{
for (auto &a : mAtoms)
a.translateRotateAndTranslate(t1, q, t2);
}
} // namespace mmcif
src/TlsParser.cpp
......@@ -248,7 +248,7 @@ struct TLSSelectionNot : public TLSSelection
for (auto& r: residues)
r.selected = not r.selected;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "NOT" << std::endl;
DumpSelection(residues, indentLevel);
......@@ -267,7 +267,7 @@ struct TLSSelectionAll : public TLSSelection
for (auto& r: residues)
r.selected = true;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "ALL" << std::endl;
DumpSelection(residues, indentLevel);
......@@ -287,7 +287,7 @@ struct TLSSelectionChain : public TLSSelectionAll
for (auto& r: residues)
r.selected = allChains or r.chainID == m_chain;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "CHAIN " << m_chain << std::endl;
DumpSelection(residues, indentLevel);
......@@ -307,7 +307,7 @@ struct TLSSelectionResID : public TLSSelectionAll
for (auto& r: residues)
r.selected = r.seqNr == m_seq_nr and r.iCode == m_icode;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "ResID " << m_seq_nr << (m_icode ? std::string { m_icode} : "") << std::endl;
DumpSelection(residues, indentLevel);
......@@ -331,7 +331,7 @@ struct TLSSelectionRangeSeq : public TLSSelectionAll
(r.seqNr <= m_last or m_last == kResidueNrWildcard));
}
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Range " << m_first << ':' << m_last << std::endl;
DumpSelection(residues, indentLevel);
......@@ -374,7 +374,7 @@ struct TLSSelectionRangeID : public TLSSelectionAll
}
}
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Through " << m_first << ':' << m_last << std::endl;
DumpSelection(residues, indentLevel);
......@@ -407,7 +407,7 @@ struct TLSSelectionUnion : public TLSSelection
for (auto ai = a.begin(), bi = b.begin(), ri = residues.begin(); ri != residues.end(); ++ai, ++bi, ++ri)
ri->selected = ai->selected or bi->selected;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Union" << std::endl;
DumpSelection(residues, indentLevel);
......@@ -440,7 +440,7 @@ struct TLSSelectionIntersection : public TLSSelection
for (auto ai = a.begin(), bi = b.begin(), ri = residues.begin(); ri != residues.end(); ++ai, ++bi, ++ri)
ri->selected = ai->selected and bi->selected;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Intersection" << std::endl;
DumpSelection(residues, indentLevel);
......@@ -462,7 +462,7 @@ struct TLSSelectionByName : public TLSSelectionAll
for (auto& r: residues)
r.selected = r.name == m_name;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Name " << m_name << std::endl;
DumpSelection(residues, indentLevel);
......@@ -488,7 +488,7 @@ struct TLSSelectionByElement : public TLSSelectionAll
for (auto& r: residues)
r.selected = iequals(r.name, m_element);
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
{
std::cout << std::string(indentLevel * 2, ' ') << "Element " << m_element << std::endl;
DumpSelection(residues, indentLevel);
......@@ -890,7 +890,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::Parse()
Match(pt_EOLN);
if (extraParenthesis)
if (extraParenthesis and cif::VERBOSE > 0)
std::cerr << "WARNING: too many closing parenthesis in TLS selection statement" << std::endl;
return result;
......@@ -931,7 +931,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::ParseFactor()
case '(':
Match('(');
result = ParseAtomSelection();
if (m_lookahead == pt_EOLN)
if (m_lookahead == pt_EOLN and cif::VERBOSE > 0)
std::cerr << "WARNING: missing closing parenthesis in TLS selection statement" << std::endl;
else
Match(')');
......@@ -1033,7 +1033,7 @@ TLSSelectionPtr TLSSelectionParserImplPhenix::ParseFactor()
result.reset(new TLSSelectionRangeID(from, to, icode_from, icode_to));
else
{
if (cif::VERBOSE and (icode_from or icode_to))
if (cif::VERBOSE > 0 and (icode_from or icode_to))
std::cerr << "Warning, ignoring insertion codes" << std::endl;
result.reset(new TLSSelectionRangeSeq(from, to));
......@@ -1231,7 +1231,8 @@ TLSSelectionPtr TLSSelectionParserImplBuster::ParseGroup()
std::tie(chain2, seqNr2) = ParseAtom();
if (chain1 != chain2)
{
std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
if (cif::VERBOSE > 0)
std::cerr << "Warning, ranges over multiple chains detected" << std::endl;
TLSSelectionPtr sc1(new TLSSelectionChain(chain1));
TLSSelectionPtr sr1(new TLSSelectionRangeSeq(seqNr1, kResidueNrWildcard));
......@@ -1289,7 +1290,7 @@ std::tuple<std::string,int> TLSSelectionParserImplBuster::ParseAtom()
Match(':');
std::string atom = m_value_s;
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Warning: ignoring atom ID '" << atom << "' in TLS selection" << std::endl;
Match(bt_IDENT);
......@@ -1810,7 +1811,8 @@ class TLSSelectionParser
}
catch (const std::exception& ex)
{
std::cerr << "ParseError: " << ex.what() << std::endl;
if (cif::VERBOSE >= 0)
std::cerr << "ParseError: " << ex.what() << std::endl;
}
return result;
......@@ -1834,14 +1836,14 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to old BUSTER" << std::endl;
result = busterOld.Parse(selection);
}
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to PHENIX" << std::endl;
result = phenix.Parse(selection);
}
......@@ -1852,35 +1854,35 @@ TLSSelectionPtr ParseSelectionDetails(const std::string& program, const std::str
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to BUSTER" << std::endl;
result = buster.Parse(selection);
}
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to old BUSTER" << std::endl;
result = busterOld.Parse(selection);
}
}
else
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "No known program specified, trying PHENIX" << std::endl;
result = phenix.Parse(selection);
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to BUSTER" << std::endl;
result = buster.Parse(selection);
}
if (not result)
{
if (cif::VERBOSE)
if (cif::VERBOSE > 0)
std::cerr << "Falling back to old BUSTER" << std::endl;
result = busterOld.Parse(selection);
}
......
test/structure-test.cpp
......@@ -179,3 +179,23 @@ _struct_asym.details ?
<< structure.getFile().data() << std::endl;
}
}
// // --------------------------------------------------------------------
// BOOST_AUTO_TEST_CASE(test_load_1)
// {
// mmcif::File cf(gTestDir / "5v3g.cif.gz");
// mmcif::Structure s(cf);
// for (auto &poly : s.polymers())
// {
// std::cout << std::string(80, '=') << std::endl;
// for (auto &res : poly)
// {
// std::cout << res << std::endl;
// for (auto &atom : res.atoms())
// std::cout << " " << atom << std::endl;
// }
// }
// }
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment