New Upstream Release - libpdb-redo
Ready changes
Summary
Merged new upstream version: 3.0.6 (was: 3.0.5).
Diff
diff --git a/.gitignore b/.gitignore
index d0cf1a7..013fff9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@ data/
.gitignore
**/build/
src/revision.hpp
+Testing/
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0fccbd0..87f9354 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,7 +25,7 @@
cmake_minimum_required(VERSION 3.16)
# set the project name
-project(pdb-redo VERSION 3.0.5 LANGUAGES CXX C)
+project(pdb-redo VERSION 3.0.6 LANGUAGES CXX C)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake")
@@ -38,6 +38,7 @@ include(CheckIncludeFiles)
include(CheckLibraryExists)
include(CMakePackageConfigHelpers)
include(CheckCXXSourceCompiles)
+include(Dart)
include(GenerateExportHeader)
set(CXX_EXTENSIONS OFF)
@@ -108,15 +109,15 @@ if(MSVC)
endif()
# Libraries
+find_package(Eigen3 REQUIRED)
find_package(CCP4 REQUIRED ccp4c clipper-core clipper-ccp4 clipper-contrib)
+find_package(GSL REQUIRED)
if(NOT PDB_REDO_META)
find_package(newuoa REQUIRED)
- find_package(cifpp 5.0.4 REQUIRED)
+ find_package(cifpp 5.0.8 REQUIRED)
endif()
-find_package(GSL REQUIRED)
-
# Create a revision file, containing the current git version info
include(VersionString)
write_version_header("${PROJECT_SOURCE_DIR}/src" "LibPDBREDO")
@@ -133,7 +134,6 @@ set(project_sources
${PROJECT_SOURCE_DIR}/src/ResolutionCalculator.cpp
${PROJECT_SOURCE_DIR}/src/SkipList.cpp
${PROJECT_SOURCE_DIR}/src/Statistics.cpp
- ${PROJECT_SOURCE_DIR}/src/Symmetry-2.cpp
${PROJECT_SOURCE_DIR}/src/Minimizer.cpp
${PROJECT_SOURCE_DIR}/src/Restraints.cpp
)
@@ -149,7 +149,6 @@ set(project_headers
${PROJECT_SOURCE_DIR}/include/pdb-redo/ResolutionCalculator.hpp
${PROJECT_SOURCE_DIR}/include/pdb-redo/SkipList.hpp
${PROJECT_SOURCE_DIR}/include/pdb-redo/Statistics.hpp
- ${PROJECT_SOURCE_DIR}/include/pdb-redo/Symmetry-2.hpp
${PROJECT_SOURCE_DIR}/include/pdb-redo/Minimizer.hpp
${PROJECT_SOURCE_DIR}/include/pdb-redo/Restraints.hpp
)
@@ -161,10 +160,10 @@ generate_export_header(pdb-redo EXPORT_FILE_NAME pdb-redo/exports.hpp)
target_include_directories(pdb-redo
PUBLIC
- "$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${PROJECT_BINARY_DIR}>"
+ "$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include;${EIGEN3_INCLUDE_DIR}>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
${CCP4_INCLUDE_DIRS})
-target_link_libraries(pdb-redo PUBLIC cifpp::cifpp newuoa::newuoa ${CCP4_LIBRARIES} ${GSL_LIBRARIES})
+target_link_libraries(pdb-redo PUBLIC cifpp::cifpp newuoa::newuoa ${CCP4_LIBRARIES} GSL::gsl)
if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
target_link_options(pdb-redo PRIVATE -undefined dynamic_lookup)
diff --git a/changelog b/changelog
index 1795c22..edd474f 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+Version 3.0.6
+- Fix in calculating symmetry atom locations
+
Version 3.0.5
- Refinement work
- Fix regression in bondmap voor sugars
diff --git a/cmake/pdb-redoConfig.cmake.in b/cmake/pdb-redoConfig.cmake.in
index 8089aea..ac694da 100644
--- a/cmake/pdb-redoConfig.cmake.in
+++ b/cmake/pdb-redoConfig.cmake.in
@@ -4,9 +4,7 @@ include(CMakeFindDependencyMacro)
find_dependency(newuoa 0.1.2 REQUIRED)
find_dependency(cifpp 5.0.4 REQUIRED)
find_dependency(CCP4 REQUIRED ccp4c clipper-core clipper-ccp4 clipper-contrib)
-
-include(FindPkgConfig)
-pkg_check_modules(GSL gsl REQUIRED)
+find_dependency(GSL REQUIRED)
INCLUDE("${CMAKE_CURRENT_LIST_DIR}/pdb-redoTargets.cmake")
diff --git a/debian/changelog b/debian/changelog
index 7d9990a..9819f02 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+libpdb-redo (3.0.6-1) UNRELEASED; urgency=low
+
+ * New upstream release.
+
+ -- Debian Janitor <janitor@jelmer.uk> Tue, 06 Jun 2023 19:16:43 -0000
+
libpdb-redo (3.0.5-2) unstable; urgency=medium
* Fix libcifpp dependency.
diff --git a/include/pdb-redo/BondMap.hpp b/include/pdb-redo/BondMap.hpp
index 4af2506..f361021 100644
--- a/include/pdb-redo/BondMap.hpp
+++ b/include/pdb-redo/BondMap.hpp
@@ -49,20 +49,24 @@ class BondMapException : public std::runtime_error
class BondMap
{
public:
- BondMap(const cif::mm::structure &structure)
- : BondMap(structure.get_datablock(), structure.get_model_nr())
+ BondMap(const cif::mm::structure &structure, std::optional<std::tuple<cif::point,float>> around = {})
+ : BondMap(structure.get_datablock(), around, structure.get_model_nr())
{
-
}
- BondMap(const cif::datablock &db, size_t model_nr = 1);
+ BondMap(const cif::datablock &db, std::optional<std::tuple<cif::point,float>> around = {}, size_t model_nr = 1);
BondMap(const BondMap &) = delete;
BondMap &operator=(const BondMap &) = delete;
+ BondMap(BondMap &&);
+ BondMap &operator=(BondMap &&);
+
bool operator()(const std::string &atom_1, const std::string &atom_2) const
{
- return isBonded(index.at(atom_1), index.at(atom_2));
+ auto aix1 = index.find(atom_1);
+ auto aix2 = index.find(atom_2);
+ return aix1 != index.end() and aix2 != index.end() and isBonded(aix1->second, aix2->second);
}
bool operator()(const cif::mm::atom &atom_1, const cif::mm::atom &atom_2) const
diff --git a/include/pdb-redo/ClipperWrapper.hpp b/include/pdb-redo/ClipperWrapper.hpp
index b97d42a..7f82b27 100644
--- a/include/pdb-redo/ClipperWrapper.hpp
+++ b/include/pdb-redo/ClipperWrapper.hpp
@@ -35,16 +35,6 @@ namespace pdb_redo
clipper::Atom toClipper(const cif::mm::atom &atom);
clipper::Atom toClipper(cif::row_handle atom, cif::row_handle aniso_row);
-inline clipper::Coord_orth toClipper(const cif::point &pt)
-{
- return { pt.m_x, pt.m_y, pt.m_z };
-}
-
-inline cif::point toPoint(const clipper::Coord_orth &pt)
-{
- return { static_cast<float>(pt.x()), static_cast<float>(pt.y()), static_cast<float>(pt.z()) };
-}
-
// --------------------------------------------------------------------
clipper::Spacegroup getSpacegroup(const cif::datablock &db);
diff --git a/include/pdb-redo/Compound.hpp b/include/pdb-redo/Compound.hpp
index 7545550..5d9f0aa 100644
--- a/include/pdb-redo/Compound.hpp
+++ b/include/pdb-redo/Compound.hpp
@@ -192,10 +192,7 @@ class Compound
// std::vector<std::tuple<std::string, std::string>> mapToIsomer(const Compound &c) const;
private:
- ~Compound();
-
cif::file mCF;
-
std::string mID;
std::string mName;
std::string mGroup;
@@ -302,7 +299,6 @@ class Link
float chiralVolume(const std::string &id, const std::string &compound_id_1, const std::string &compound_id_2) const;
private:
- ~Link();
std::string mID;
std::vector<LinkBond> mBonds;
diff --git a/include/pdb-redo/DistanceMap.hpp b/include/pdb-redo/DistanceMap.hpp
index f4e5750..4dcd156 100644
--- a/include/pdb-redo/DistanceMap.hpp
+++ b/include/pdb-redo/DistanceMap.hpp
@@ -26,13 +26,9 @@
#pragma once
-#include <unordered_map>
-
-#include <clipper/clipper.h>
-
#include <cif++.hpp>
-#include <pdb-redo/ClipperWrapper.hpp>
+#include <unordered_map>
#ifdef near
#undef near
@@ -45,11 +41,10 @@ class DistanceMap
{
public:
- DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell,
- float maxDistance);
+ DistanceMap(const cif::mm::structure &p, const cif::crystal &crystal, float maxDistance);
DistanceMap(const cif::mm::structure &p, float maxDistance)
- : DistanceMap(p, getSpacegroup(p.get_datablock()), getCell(p.get_datablock()), maxDistance)
+ : DistanceMap(p, cif::crystal(p.get_datablock()), maxDistance)
{
}
@@ -60,30 +55,28 @@ class DistanceMap
std::vector<cif::mm::atom> near(const cif::mm::atom &atom, float maxDistance = 3.5f) const;
- static std::vector<clipper::RTop_orth>
- AlternativeSites(const clipper::Spacegroup &spacegroup, const clipper::Cell &cell);
-
private:
using DistKeyType = std::tuple<size_t, size_t>;
- using DistValueType = std::tuple<float, int32_t>;
+ using DistValueType = std::tuple<float, cif::sym_op, bool>;
using DistMap = std::map<DistKeyType, DistValueType>;
void AddDistancesForAtoms(const std::vector<std::tuple<size_t,cif::point>> &a,
- const std::vector<std::tuple<size_t,cif::point>> &b, DistMap &dm, int32_t rtix);
+ const std::vector<std::tuple<size_t,cif::point>> &b, DistMap &dm);
+ void AddDistancesForAtoms(const std::vector<std::tuple<size_t,cif::point>> &a,
+ const std::vector<std::tuple<size_t,cif::point>> &b, DistMap &dm, cif::sym_op symop);
+
+ cif::point offsetToOrigin(const cif::point &p) const;
const cif::mm::structure &mStructure;
- clipper::Cell cell;
- clipper::Spacegroup spacegroup;
+ cif::crystal crystal;
size_t dim;
std::unordered_map<std::string, size_t> index;
std::map<size_t, std::string> rIndex;
float mMaxDistance, mMaxDistanceSQ;
- std::vector<std::tuple<float, int32_t>> mA;
+ std::vector<std::tuple<float, cif::sym_op, bool>> mA;
std::vector<size_t> mIA, mJA;
- cif::point mD; // needed to move atoms to center
- std::vector<clipper::RTop_orth> mRtOrth;
};
} // namespace pdb_redo
diff --git a/include/pdb-redo/Minimizer.hpp b/include/pdb-redo/Minimizer.hpp
index 799c55f..6fbb020 100644
--- a/include/pdb-redo/Minimizer.hpp
+++ b/include/pdb-redo/Minimizer.hpp
@@ -95,24 +95,49 @@ class Minimizer
virtual ~Minimizer() {}
// factory method:
- static Minimizer *create(const cif::mm::polymer &poly, int first, int last,
- const BondMap &bm, const XMap &xMap, float mapWeight = 60, float plane5AtomsESD = 0.11);
+ static Minimizer *create(const cif::crystal &crystal, const cif::mm::polymer &poly, int first, int last, const XMap &xMap);
- static Minimizer *create(cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms,
- const BondMap &bm, const XMap &xMap, float mapWeight = 60, float plane5AtomsESD = 0.11)
+ static Minimizer *create(const cif::crystal &crystal, cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms, const XMap &xMap)
{
- return create(structure, atoms, bm, plane5AtomsESD, &xMap, mapWeight);
+ return create(crystal, structure, atoms, &xMap);
}
// factory method for minimizer without density:
- static Minimizer *create(cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms,
- const BondMap &bm, float plane5AtomsESD = 0.11)
+ static Minimizer *create(const cif::crystal &crystal, cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms)
{
- return create(structure, atoms, bm, plane5AtomsESD, nullptr, 0);
+ return create(crystal, structure, atoms, nullptr);
}
+ // Drop all torsion restraints
void dropTorsionRestraints();
+ // Filter based on result of callback
+ // Signature of callback should be: bool (*filter)(cif::mm::atom a1, cif::mm::atom a2, cif::mm::atom a3, cif::mm::atom a4)
+ template <typename F>
+ void filterTorsionRestraints(F &&cb)
+ {
+ auto e = std::remove_if(mTorsionRestraints.begin(), mTorsionRestraints.end(),
+ [this, cb = std::move(cb)](TorsionRestraint &r)
+ {
+ return r.mA >= mAtoms.size() or r.mB >= mAtoms.size() or r.mC >= mAtoms.size() or r.mD >= mAtoms.size() or
+ cb(mAtoms[r.mA], mAtoms[r.mB], mAtoms[r.mC], mAtoms[r.mD]);
+ });
+
+ for (auto i = e; i != mTorsionRestraints.end(); ++i)
+ mRestraints.erase(std::remove(mRestraints.begin(), mRestraints.end(), &*i), mRestraints.end());
+
+ mTorsionRestraints.erase(e, mTorsionRestraints.end());
+ }
+
+ // Set the map weight, default is 60
+ void setMapWeight(float mapWeight);
+
+ // Set the chiral volume ESD, default is 0.2
+ void setChiralVolumeESD(float chiralityESD);
+
+ // Set the planarity ESD, default is 0.11
+ void setPlanarityESD(float planarityESD);
+
void printStats();
virtual double refine(bool storeAtoms) = 0;
@@ -121,15 +146,15 @@ class Minimizer
virtual void storeAtomLocations() = 0;
protected:
- Minimizer(const cif::mm::structure &structure, const BondMap &bm, float plane5AtomsESD);
- static Minimizer *create(cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms,
- const BondMap &bm, float plane5AtomsESD, const XMap *xMap, float mapWeight);
+ Minimizer(const cif::mm::structure &structure);
+
+ static Minimizer *create(const cif::crystal &crystal, cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms, const XMap *xMap);
virtual void addResidue(const cif::mm::residue &res);
virtual void addPolySection(const cif::mm::polymer &poly, int first, int last);
virtual void addDensityMap(const XMap &xMap, float mapWeight);
- virtual void Finish();
+ virtual void Finish(const cif::crystal &crystal);
double score(const AtomLocationProvider &loc);
@@ -161,11 +186,11 @@ class Minimizer
AtomRef ref(const cif::mm::atom &atom);
+ BondMap createBondMap();
+
bool mElectronScattering = false; // TODO: use!
const cif::mm::structure &mStructure;
- const BondMap &mBonds;
- float mPlane5ESD;
std::vector<cif::mm::atom> mAtoms, mReferencedAtoms;
std::vector<size_t> mRef2AtomIndex;
diff --git a/include/pdb-redo/Restraints.hpp b/include/pdb-redo/Restraints.hpp
index b898aa3..6830a90 100644
--- a/include/pdb-redo/Restraints.hpp
+++ b/include/pdb-redo/Restraints.hpp
@@ -134,6 +134,9 @@ struct TransPeptideRestraint : public TorsionRestraint
}
};
+const double kChiralVolumeESD = 0.2; // according to coot that's a reasonable value...
+
+
struct ChiralVolumeRestraint : public Restraint
{
ChiralVolumeRestraint(AtomRef c, AtomRef a1, AtomRef a2, AtomRef a3, double volume)
@@ -150,7 +153,7 @@ struct ChiralVolumeRestraint : public Restraint
virtual void print(const AtomLocationProvider &atoms) const;
AtomRef mCentre, mA1, mA2, mA3;
- double mVolume;
+ double mVolume, mESD = kChiralVolumeESD;
};
struct PlanarityRestraint : public Restraint
diff --git a/include/pdb-redo/Statistics.hpp b/include/pdb-redo/Statistics.hpp
index 7532e88..ce71146 100644
--- a/include/pdb-redo/Statistics.hpp
+++ b/include/pdb-redo/Statistics.hpp
@@ -27,7 +27,6 @@
#pragma once
#include <pdb-redo/BondMap.hpp>
-#include "pdb-redo/DistanceMap.hpp"
#include "pdb-redo/MapMaker.hpp"
namespace pdb_redo
@@ -89,8 +88,10 @@ class StatsCollector
virtual std::vector<ResidueStatistics> collect() const;
+ virtual std::vector<ResidueStatistics> collect(const std::string &asymID) const;
+
virtual std::vector<ResidueStatistics> collect(const std::string &asymID,
- int resFirst, int resLast, bool authNameSpace = true) const;
+ int resFirst, int resLast, bool authNameSpace = false) const;
virtual ResidueStatistics collect(std::initializer_list<const cif::mm::residue *> residues) const;
@@ -151,15 +152,14 @@ class StatsCollector
class EDIAStatsCollector : public StatsCollector
{
public:
- EDIAStatsCollector(MapMaker<float> &mm,
- cif::mm::structure &structure, bool electronScattering,
- const BondMap &bondMap);
+ EDIAStatsCollector(const MapMaker<float> &mm,
+ cif::mm::structure &structure, bool electronScattering);
protected:
virtual void calculate(std::vector<AtomData> &atomData) const;
- DistanceMap mDistanceMap;
- const BondMap &mBondMap;
+ BondMap createBondMap(std::vector<AtomData> &atomData) const;
+
std::map<cif::atom_type, float> mRadii;
};
diff --git a/include/pdb-redo/Symmetry-2.hpp b/include/pdb-redo/Symmetry-2.hpp
deleted file mode 100644
index eb469db..0000000
--- a/include/pdb-redo/Symmetry-2.hpp
+++ /dev/null
@@ -1,214 +0,0 @@
-/*-
- * 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
- * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#pragma once
-
-#include <cif++.hpp>
-
-#include "pdb-redo/ClipperWrapper.hpp"
-
-namespace pdb_redo
-{
-
-// --------------------------------------------------------------------
-// Functions to use when working with symmetry stuff
-
-clipper::Coord_orth CalculateOffsetForCell(const cif::mm::structure &structure, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell);
-std::vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup &spacegroup, const clipper::Cell &cell);
-// int GetSpacegroupNumber(std::string spacegroup); // alternative for clipper's parsing code
-// std::string SpacegroupToHall(std::string spacegroup);
-
-cif::mm::atom symmetryCopy(const cif::mm::atom &atom, const cif::point &d,
- const clipper::Spacegroup &spacegroup, const clipper::Cell &cell, const clipper::RTop_orth &rt);
-
-std::string describeRToperation(const clipper::Spacegroup &spacegroup, const clipper::Cell &cell, const clipper::RTop_orth &rt);
-
-// --------------------------------------------------------------------
-// To iterate over all symmetry copies of an atom
-
-class SymmetryAtomIteratorFactory
-{
- public:
- SymmetryAtomIteratorFactory(const cif::mm::structure &structure, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell);
- SymmetryAtomIteratorFactory(const cif::mm::structure &structure, int spacegroupNr, const clipper::Cell &cell);
-
- SymmetryAtomIteratorFactory(const SymmetryAtomIteratorFactory &) = delete;
- SymmetryAtomIteratorFactory &operator=(const SymmetryAtomIteratorFactory &) = delete;
-
- template <typename ConditionFunc>
- class SymmetryAtomIterator
- {
- public:
- using iterator_category = std::forward_iterator_tag;
- using value_type = cif::mm::atom;
- using difference_type = std::ptrdiff_t;
- using pointer = value_type*;
- using reference = value_type &;
-
- SymmetryAtomIterator(const SymmetryAtomIteratorFactory &factory, const cif::mm::atom &atom, ConditionFunc &cond)
- : m_f(&factory)
- , m_i(0)
- , m_a(atom)
- , m_c(atom)
- , m_cond(cond)
- {
- while (not test() and m_i < m_f->mRtOrth.size())
- ++m_i;
- }
-
- SymmetryAtomIterator(const SymmetryAtomIteratorFactory &factory, const cif::mm::atom &atom, ConditionFunc &cond, int)
- : SymmetryAtomIterator(factory, atom, cond)
- {
- m_i = m_f->mRtOrth.size();
- }
-
- SymmetryAtomIterator(const SymmetryAtomIterator &iter)
- : m_f(iter.m_f)
- , m_i(iter.m_i)
- , m_a(iter.m_a)
- , m_c(iter.m_c)
- , m_cond(iter.m_cond)
- {
- }
-
- SymmetryAtomIterator &operator=(const SymmetryAtomIterator &iter)
- {
- if (this != &iter)
- {
- m_f = iter.m_f;
- m_i = iter.m_i;
- m_a = iter.m_a;
- m_c = iter.m_c;
- }
- return *this;
- }
-
- reference operator*() { return m_c; }
- pointer operator->() { return &m_c; }
-
- SymmetryAtomIterator operator++()
- {
- while (m_i < m_f->mRtOrth.size())
- {
- ++m_i;
- if (test())
- break;
- }
-
- return *this;
- }
-
- SymmetryAtomIterator operator++(int)
- {
- SymmetryAtomIterator result(*this);
- this->operator++();
- return result;
- }
-
- bool operator==(const SymmetryAtomIterator &iter) const
- {
- return m_f == iter.m_f and m_i == iter.m_i;
- }
-
- bool operator!=(const SymmetryAtomIterator &iter) const
- {
- return m_f != iter.m_f or m_i != iter.m_i;
- }
-
- private:
- bool test()
- {
- bool result = false;
-
- if (m_i < m_f->mRtOrth.size())
- {
- auto &rt = m_f->mRtOrth[m_i];
- auto loc = m_a.get_location();
-
- loc += m_f->mD;
- loc = toPoint(toClipper(loc).transform(rt));
- loc -= m_f->mD;
-
- if (m_cond(loc))
- {
- std::string rt_operation = describeRToperation(m_f->mSpacegroup, m_f->mCell, rt);
- m_c = cif::mm::atom(m_a, loc, rt_operation);
- result = true;
- }
- }
-
- return result;
- }
-
- const SymmetryAtomIteratorFactory *m_f;
- size_t m_i;
- cif::mm::atom m_a, m_c;
- ConditionFunc &m_cond;
- };
-
- template <typename ConditionFunc>
- class SymmetryAtomIteratorRange
- {
- public:
- SymmetryAtomIteratorRange(const SymmetryAtomIteratorFactory &f, const cif::mm::atom &a, ConditionFunc &&cond)
- : m_f(f)
- , m_a(a)
- , m_cond(std::move(cond))
- {
- }
-
- SymmetryAtomIterator<ConditionFunc> begin()
- {
- return SymmetryAtomIterator(m_f, m_a, m_cond);
- }
-
- SymmetryAtomIterator<ConditionFunc> end()
- {
- return SymmetryAtomIterator(m_f, m_a, m_cond, 1);
- }
-
- private:
- const SymmetryAtomIteratorFactory &m_f;
- cif::mm::atom m_a;
- ConditionFunc m_cond;
- };
-
- template <typename ConditionFunc>
- SymmetryAtomIteratorRange<ConditionFunc> operator()(const cif::mm::atom &a, ConditionFunc &&cond) const
- {
- return SymmetryAtomIteratorRange(*this, a, std::forward<ConditionFunc>(cond));
- }
-
- // std::string symop_mmcif(const cif::mm::atom& a) const;
-
- private:
- clipper::Spacegroup mSpacegroup;
- cif::point mD; // needed to move atoms to center
- std::vector<clipper::RTop_orth> mRtOrth;
- clipper::Cell mCell;
-};
-
-} // namespace pdb_redo
\ No newline at end of file
diff --git a/src/BondMap.cpp b/src/BondMap.cpp
index 6518deb..8d6700a 100644
--- a/src/BondMap.cpp
+++ b/src/BondMap.cpp
@@ -141,7 +141,27 @@ bool CompoundBondMap::bonded(const std::string &compoundID, const std::string &a
// --------------------------------------------------------------------
-BondMap::BondMap(const cif::datablock &db, size_t model_nr)
+BondMap::BondMap(BondMap &&bm)
+ : dim(bm.dim)
+ , index(std::move(bm.index))
+ , bond(std::move(bm.bond))
+ , bond_1_4(std::move(bm.bond_1_4))
+ , link(std::move(bm.link))
+{
+}
+
+BondMap &BondMap::operator=(BondMap &&bm)
+{
+ dim = bm.dim;
+ index.swap(bm.index);
+ bond.swap(bm.bond);
+ bond_1_4.swap(bm.bond_1_4);
+ link.swap(bm.link);
+
+ return *this;
+}
+
+BondMap::BondMap(const cif::datablock &db, std::optional<std::tuple<cif::point,float>> around, size_t model_nr)
{
using namespace cif::literals;
@@ -150,8 +170,23 @@ BondMap::BondMap(const cif::datablock &db, size_t model_nr)
// First collect the atoms from the datablock
std::vector<cif::row_handle> atoms;
+ cif::crystal crystal(db);
+
for (auto rh : db["atom_site"].find("pdbx_PDB_model_num"_key == model_nr or "pdbx_PDB_model_num"_key == cif::null))
- atoms.push_back(rh);
+ {
+ if (around)
+ {
+ const auto &[p, r] = *around;
+ const auto &[x, y, z] = rh.get<float,float,float>("Cartn_x", "Cartn_y", "Cartn_z");
+
+ const auto &[d, pt, op] = crystal.closest_symmetry_copy(p, {x, y, z});
+
+ if (d <= r)
+ atoms.push_back(rh);
+ }
+ else
+ atoms.push_back(rh);
+ }
dim = uint32_t(atoms.size());
@@ -190,7 +225,7 @@ BondMap::BondMap(const cif::datablock &db, size_t model_nr)
compounds.insert(mon_id);
}
- cif::Progress progress(compounds.size(), "Creating bond map");
+ cif::progress_bar progress_bar(compounds.size(), "Creating bond map");
// some helper indices to speed things up a bit
using atom_map_key_type = std::tuple<std::string, int, std::string, std::string>;
@@ -249,6 +284,8 @@ BondMap::BondMap(const cif::datablock &db, size_t model_nr)
for (auto c : compounds)
{
+ progress_bar.consumed(1);
+
if (c == "HOH" or c == "H2O" or c == "WAT")
{
if (cif::VERBOSE > 1)
diff --git a/src/ClipperWrapper.cpp b/src/ClipperWrapper.cpp
index ce1f747..95c141d 100644
--- a/src/ClipperWrapper.cpp
+++ b/src/ClipperWrapper.cpp
@@ -49,13 +49,18 @@ clipper::Atom toClipper(cif::row_handle atom, cif::row_handle aniso_row)
if (not atom["pdbx_formal_charge"].empty())
{
int charge = atom["pdbx_formal_charge"].as<int>();
- if (std::abs(charge) > 1)
- element += std::to_string(charge);
- if (charge < 0)
- element += '-';
- else
- element += '+';
+
+ if (cif::atom_type_traits(element).has_sf(charge) and charge != 0)
+ {
+ element += std::to_string(std::abs(charge));
+
+ if (charge < 0)
+ element += '-';
+ else
+ element += '+';
+ }
}
+
result.set_element(element);
if (not atom["U_iso_or_equiv"].empty())
@@ -86,7 +91,7 @@ clipper::Atom toClipper(const cif::mm::atom &atom)
clipper::Spacegroup getSpacegroup(const cif::datablock &db)
{
- std::string spacegroup = db["symmetry"].find1<std::string>(cif::key("entry_id") == db.name(), "space_group_name_H-M");
+ std::string spacegroup = db["symmetry"].find_first<std::string>(cif::key("entry_id") == db.name(), "space_group_name_H-M");
if (spacegroup == "P 1-")
spacegroup = "P -1";
diff --git a/src/Compound.cpp b/src/Compound.cpp
index 549ebd0..e7f229b 100644
--- a/src/Compound.cpp
+++ b/src/Compound.cpp
@@ -541,6 +541,8 @@ float Compound::chiralVolume(const std::string ¢reID) const
Link::Link(cif::datablock &db)
{
mID = db.name();
+ if (cif::starts_with(mID, "link_"))
+ mID.erase(mID.begin(), mID.begin() + 5);
auto &linkBonds = db["chem_link_bond"];
@@ -659,9 +661,7 @@ Link::Link(cif::datablock &db)
const Link &Link::create(const std::string &id)
{
- auto result = CompoundFactory::instance().getLink(id);
- if (result == nullptr)
- result = CompoundFactory::instance().createLink(id);
+ auto result = CompoundFactory::instance().createLink(id);
if (result == nullptr)
throw std::runtime_error("Link with id " + id + " not found");
@@ -669,10 +669,6 @@ const Link &Link::create(const std::string &id)
return *result;
}
-Link::~Link()
-{
-}
-
float Link::atomBondValue(const LinkAtom &atom1, const LinkAtom &atom2) const
{
auto i = find_if(mBonds.begin(), mBonds.end(),
@@ -833,8 +829,8 @@ class CompoundFactoryImpl
delete mNext;
}
- Compound *get(std::string id);
- Compound *create(std::string id);
+ const Compound *get(std::string id);
+ const Compound *create(std::string id);
const Link *getLink(std::string id);
const Link *createLink(std::string id);
@@ -885,8 +881,8 @@ class CompoundFactoryImpl
std::shared_timed_mutex mMutex;
std::string mPath;
- std::vector<Compound *> mCompounds;
- std::vector<const Link *> mLinks;
+ std::vector<std::unique_ptr<const Compound>> mCompounds;
+ std::vector<std::unique_ptr<const Link>> mLinks;
std::set<std::string> mKnownPeptides;
std::set<std::string> mKnownBases;
std::set<std::string> mMissing;
@@ -948,19 +944,19 @@ CompoundFactoryImpl::CompoundFactoryImpl(std::istream &data, CompoundFactoryImpl
}
}
-Compound *CompoundFactoryImpl::get(std::string id)
+const Compound *CompoundFactoryImpl::get(std::string id)
{
std::shared_lock lock(mMutex);
cif::to_upper(id);
- Compound *result = nullptr;
+ const Compound *result = nullptr;
- for (auto cmp : mCompounds)
+ for (auto &cmp : mCompounds)
{
if (cmp->id() == id)
{
- result = cmp;
+ result = cmp.get();
break;
}
}
@@ -971,11 +967,11 @@ Compound *CompoundFactoryImpl::get(std::string id)
return result;
}
-Compound *CompoundFactoryImpl::create(std::string id)
+const Compound *CompoundFactoryImpl::create(std::string id)
{
cif::to_upper(id);
- Compound *result = get(id);
+ const Compound *result = get(id);
if (result == nullptr and mMissing.count(id) == 0 and not mFile.empty())
{
std::unique_lock lock(mMutex);
@@ -1009,14 +1005,14 @@ Compound *CompoundFactoryImpl::create(std::string id)
mMissing.insert(id);
else
{
- mCompounds.push_back(new Compound(resFile.string(), id, name, group));
- result = mCompounds.back();
+ mCompounds.emplace_back(new Compound(resFile.string(), id, name, group));
+ result = mCompounds.back().get();
}
}
else
{
- mCompounds.push_back(new Compound(mPath, id, name, group));
- result = mCompounds.back();
+ mCompounds.emplace_back(new Compound(mPath, id, name, group));
+ result = mCompounds.back().get();
}
}
@@ -1035,11 +1031,11 @@ const Link *CompoundFactoryImpl::getLink(std::string id)
const Link *result = nullptr;
- for (auto link : mLinks)
+ for (auto &link : mLinks)
{
- if (link->id() == id)
+ if (cif::iequals(link->id(), id))
{
- result = link;
+ result = link.get();
break;
}
}
@@ -1061,10 +1057,7 @@ const Link *CompoundFactoryImpl::createLink(std::string id)
std::unique_lock lock(mMutex);
if (mFile.contains("link_" + id))
- {
- result = new Link(mFile["link_" + id]);
- mLinks.push_back(result);
- }
+ result = mLinks.emplace_back(std::make_unique<Link>(mFile["link_" + id])).get();
if (result == nullptr and mNext != nullptr)
result = mNext->createLink(id);
diff --git a/src/DistanceMap.cpp b/src/DistanceMap.cpp
index 644ea2e..4a8b63e 100644
--- a/src/DistanceMap.cpp
+++ b/src/DistanceMap.cpp
@@ -24,14 +24,12 @@
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
-#include <atomic>
-#include <mutex>
+#include "pdb-redo/DistanceMap.hpp"
#include <cif++/utilities.hpp>
-#include "pdb-redo/ClipperWrapper.hpp"
-#include "pdb-redo/DistanceMap.hpp"
-#include "pdb-redo/Symmetry-2.hpp"
+#include <atomic>
+#include <mutex>
namespace pdb_redo
{
@@ -40,38 +38,6 @@ using cif::point;
// --------------------------------------------------------------------
-std::vector<clipper::RTop_orth> DistanceMap::AlternativeSites(const clipper::Spacegroup &spacegroup,
- const clipper::Cell &cell)
-{
- std::vector<clipper::RTop_orth> result;
-
- // to make the operation at index 0 equal to identity
- result.push_back(clipper::RTop_orth::identity());
-
- for (int i = 0; i < spacegroup.num_symops(); ++i)
- {
- const auto &symop = spacegroup.symop(i);
-
- for (int u : { -1, 0, 1 })
- for (int v : { -1, 0, 1 })
- for (int w : { -1, 0, 1 })
- {
- if (i == 0 and u == 0 and v == 0 and w == 0)
- continue;
-
- auto rtop = clipper::RTop_frac(
- symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w))
- .rtop_orth(cell);
-
- result.push_back(std::move(rtop));
- }
- }
-
- return result;
-}
-
-// --------------------------------------------------------------------
-
std::tuple<point, float> calculateCenterAndRadius(const std::vector<std::tuple<size_t,point>> &atoms)
{
std::vector<point> pts;
@@ -93,11 +59,10 @@ std::tuple<point, float> calculateCenterAndRadius(const std::vector<std::tuple<s
// --------------------------------------------------------------------
-DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell,
+DistanceMap::DistanceMap(const cif::mm::structure &p, const cif::crystal &crystal,
float maxDistance)
: mStructure(p)
- , cell(cell)
- , spacegroup(spacegroup)
+ , crystal(crystal)
, dim(0)
, mMaxDistance(maxDistance)
, mMaxDistanceSQ(maxDistance * maxDistance)
@@ -146,63 +111,9 @@ DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup
pMax.m_z = pt.m_z;
};
- // correct locations so that the median of x, y and z are inside the cell
- std::vector<float> c(locations.size());
- auto median = [&]()
- {
- return dim % 1 == 0
- ? c[dim / 2]
- : (c[dim / 2 - 1] + c[dim / 2]) / 2;
- };
-
- transform(locations.begin(), locations.end(), c.begin(), [](point &l)
- { return static_cast<float>(l.m_x); });
- sort(c.begin(), c.end());
- float mx = median();
-
- transform(locations.begin(), locations.end(), c.begin(), [](point &l)
- { return static_cast<float>(l.m_y); });
- sort(c.begin(), c.end());
- float my = median();
-
- transform(locations.begin(), locations.end(), c.begin(), [](point &l)
- { return static_cast<float>(l.m_z); });
- sort(c.begin(), c.end());
- float mz = median();
-
- if (cif::VERBOSE > 1)
- std::cerr << "median position of atoms: " << point(mx, my, mz) << std::endl;
-
- auto calculateD = [&](float m, float c)
- {
- float d = 0;
- while (m + d < -(c / 2))
- d += c;
- while (m + d > (c / 2))
- d -= c;
- return d;
- };
-
- mD.m_x = calculateD(mx, static_cast<float>(cell.a()));
- mD.m_y = calculateD(my, static_cast<float>(cell.b()));
- mD.m_z = calculateD(mz, static_cast<float>(cell.c()));
-
- clipper::Coord_orth D = toClipper(mD);
-
- if (mD.m_x != 0 or mD.m_y != 0 or mD.m_z != 0)
- {
- if (cif::VERBOSE > 1)
- std::cerr << "moving coorinates by " << mD.m_x << ", " << mD.m_y << " and " << mD.m_z << std::endl;
-
- for_each(locations.begin(), locations.end(), [&](point &p)
- { p += mD; });
- }
-
pMin -= mMaxDistance; // extend bounding box
pMax += mMaxDistance;
- mRtOrth = AlternativeSites(spacegroup, cell);
-
DistMap dist;
std::vector<std::tuple<cif::point, float, std::vector<std::tuple<size_t,point>>>> residues;
@@ -217,28 +128,30 @@ DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup
rAtoms.emplace_back(i, locations[i]);
}
- AddDistancesForAtoms(rAtoms, rAtoms, dist, 0);
+ AddDistancesForAtoms(rAtoms, rAtoms, dist);
auto &&[center, radius] = calculateCenterAndRadius(rAtoms);
residues.emplace_back(center, radius, std::move(rAtoms));
}
// treat waters special
- std::string water_entity_id = db["entity"].find1<std::string>("type"_key == "water", "id");
-
- for (size_t i = 0; i < dim; ++i)
+ auto water_entity_id = db["entity"].find1<std::optional<std::string>>("type"_key == "water", "id");
+ if (water_entity_id.has_value())
{
- if (atoms[i]["label_entity_id"] == water_entity_id)
+ for (size_t i = 0; i < dim; ++i)
{
- auto pt = locations[i];
- residues.emplace_back(pt, 0.f, std::vector<std::tuple<size_t,point>>{ { i, pt } });
+ if (atoms[i]["label_entity_id"] == *water_entity_id)
+ {
+ auto pt = locations[i];
+ residues.emplace_back(pt, 0.f, std::vector<std::tuple<size_t,point>>{ { i, pt } });
+ }
}
}
// loop over pdbx_nonpoly_scheme
for (const auto &[asymID, entityID] : db["pdbx_nonpoly_scheme"].rows<std::string, std::string>("asym_id", "entity_id"))
{
- if (entityID == water_entity_id)
+ if (water_entity_id.has_value() and entityID == *water_entity_id)
continue;
std::vector<std::tuple<size_t,point>> rAtoms;
@@ -248,29 +161,29 @@ DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup
rAtoms.emplace_back(i, locations[i]);
}
- AddDistancesForAtoms(rAtoms, rAtoms, dist, 0);
+ AddDistancesForAtoms(rAtoms, rAtoms, dist);
auto &&[center, radius] = calculateCenterAndRadius(rAtoms);
residues.emplace_back(center, radius, std::move(rAtoms));
}
// loop over pdbx_branch_scheme
- for (const auto &[asym_id, pdb_seq_num] : db["pdbx_branch_scheme"].rows<std::string, std::string>("asym_id", "pdb_seq_num"))
+ for (const auto &[asym_id, pdb_seq_num] : db["pdbx_branch_scheme"].rows<std::string, std::string>("asym_id", "num"))
{
std::vector<std::tuple<size_t,point>> rAtoms;
for (size_t i = 0; i < dim; ++i)
{
- if (atoms[i]["label_asym_id"] == asym_id and atoms[i]["label_seq_id"] == pdb_seq_num)
+ if (atoms[i]["label_asym_id"] == asym_id and atoms[i]["auth_seq_id"] == pdb_seq_num)
rAtoms.emplace_back(i, locations[i]);
}
- AddDistancesForAtoms(rAtoms, rAtoms, dist, 0);
+ AddDistancesForAtoms(rAtoms, rAtoms, dist);
auto &&[center, radius] = calculateCenterAndRadius(rAtoms);
residues.emplace_back(center, radius, std::move(rAtoms));
}
- cif::Progress progress(residues.size() * (residues.size() - 1), "Creating distance map");
+ cif::progress_bar progress_bar((residues.size() * (residues.size() - 1)) / 2, "Creating distance map");
for (size_t i = 0; i + 1 < residues.size(); ++i)
{
@@ -278,43 +191,21 @@ DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup
for (size_t j = i + 1; j < residues.size(); ++j)
{
- progress.consumed(1);
+ progress_bar.consumed(1);
const auto &[centerJ, radiusJ, atomsJ] = residues[j];
// first case, no symmetry operations
-
auto d = distance(centerI, centerJ) - radiusI - radiusJ;
if (d < mMaxDistance)
{
- AddDistancesForAtoms(atomsI, atomsJ, dist, 0);
+ AddDistancesForAtoms(atomsI, atomsJ, dist);
continue;
}
- // now try all symmetry operations to see if we can move rj close to ri
-
- clipper::Coord_orth cI = toClipper(centerI);
- clipper::Coord_orth cJ = toClipper(centerJ);
-
- auto minR2 = d;
-
- int32_t kbest = 0;
- for (int32_t k = 1; k < static_cast<int32_t>(mRtOrth.size()); ++k)
- {
- auto &rt = mRtOrth[k];
-
- auto pJ = (cJ + D).transform(rt) - D;
- double r2 = std::sqrt((cI - pJ).lengthsq()) - radiusI - radiusJ;
-
- if (minR2 > r2)
- {
- minR2 = static_cast<float>(r2);
- kbest = k;
- }
- }
-
- if (minR2 < mMaxDistance)
- AddDistancesForAtoms(atomsI, atomsJ, dist, kbest);
+ const auto &[ds, p, symop] = crystal.closest_symmetry_copy(centerI, centerJ);
+ if (ds - radiusI - radiusJ < mMaxDistance)
+ AddDistancesForAtoms(atomsI, atomsJ, dist, symop);
}
}
@@ -350,30 +241,47 @@ DistanceMap::DistanceMap(const cif::mm::structure &p, const clipper::Spacegroup
// --------------------------------------------------------------------
-void DistanceMap::AddDistancesForAtoms(const std::vector<std::tuple<size_t,point>> &a, const std::vector<std::tuple<size_t,point>> &b, DistMap &dm, int32_t rtix)
+void DistanceMap::AddDistancesForAtoms(const std::vector<std::tuple<size_t,point>> &a, const std::vector<std::tuple<size_t,point>> &b, DistMap &dm)
{
for (const auto &[ixa, loc_a] : a)
{
- clipper::Coord_orth pa = toClipper(loc_a);
-
for (const auto &[ixb, loc_b] : b)
{
if (ixa == ixb)
continue;
- clipper::Coord_orth pb = toClipper(loc_b);
+ float d = cif::distance_squared(loc_a, loc_b);
+
+ if (d > mMaxDistanceSQ)
+ continue;
+
+ d = std::sqrt(d);
+
+ dm[std::make_tuple(ixa, ixb)] = std::make_tuple(d, cif::sym_op{}, false);
+ dm[std::make_tuple(ixb, ixa)] = std::make_tuple(d, cif::sym_op{}, false);
+ }
+ }
+}
+
+void DistanceMap::AddDistancesForAtoms(const std::vector<std::tuple<size_t,point>> &a, const std::vector<std::tuple<size_t,point>> &b,
+ DistMap &dm, cif::sym_op symop)
+{
+ for (const auto &[ixa, loc_a] : a)
+ {
+ for (const auto &[ixb, loc_b] : b)
+ {
+ if (ixa == ixb)
+ continue;
- if (rtix)
- pb = pb.transform(mRtOrth[rtix]);
+ float d = cif::distance_squared(loc_a, crystal.symmetry_copy(loc_b, symop));
- auto d = static_cast<float>((pa - pb).lengthsq());
if (d > mMaxDistanceSQ)
continue;
d = std::sqrt(d);
- dm[std::make_tuple(ixa, ixb)] = std::make_tuple(d, rtix);
- dm[std::make_tuple(ixb, ixa)] = std::make_tuple(d, -rtix);
+ dm[std::make_tuple(ixa, ixb)] = std::make_tuple(d, symop, false);
+ dm[std::make_tuple(ixb, ixa)] = std::make_tuple(d, symop, true);
}
}
}
@@ -388,7 +296,7 @@ float DistanceMap::operator()(const std::string &a, const std::string &b) const
}
catch (const std::out_of_range &ex)
{
- throw std::runtime_error("atom " + a + " not found in distance map");
+ throw std::out_of_range("atom " + a + " not found in distance map");
}
try
@@ -397,7 +305,7 @@ float DistanceMap::operator()(const std::string &a, const std::string &b) const
}
catch (const std::out_of_range &ex)
{
- throw std::runtime_error("atom " + b + " not found in distance map");
+ throw std::out_of_range("atom " + b + " not found in distance map");
}
// if (ixb < ixa)
@@ -449,9 +357,7 @@ std::vector<cif::mm::atom> DistanceMap::near(const cif::mm::atom &atom, float ma
for (size_t i = mIA[ixa]; i < mIA[ixa + 1]; ++i)
{
- float d;
- int32_t rti;
- std::tie(d, rti) = mA[i];
+ const auto &[d, symop, inverse] = mA[i];
if (d > maxDistance)
continue;
@@ -459,17 +365,24 @@ std::vector<cif::mm::atom> DistanceMap::near(const cif::mm::atom &atom, float ma
size_t ixb = mJA[i];
std::string b_id = rIndex.at(ixb);
- auto altb = atom_site.find1<std::string>("id"_key == b_id, "label_alt_id");
+ auto altb = atom_site.find_first<std::string>("id"_key == b_id, "label_alt_id");
if (altb != alta and not altb.empty() and not alta.empty())
continue;
auto atom_b = mStructure.get_atom_by_id(b_id);
- if (rti > 0)
- result.emplace_back(symmetryCopy(atom_b, mD, spacegroup, cell, mRtOrth.at(rti)));
- else if (rti < 0)
- result.emplace_back(symmetryCopy(atom_b, mD, spacegroup, cell, mRtOrth.at(-rti).inverse()));
+ if (symop)
+ {
+ cif::point p = atom_b.get_location();
+
+ if (inverse)
+ p = crystal.inverse_symmetry_copy(p, symop);
+ else
+ p = crystal.symmetry_copy(p, symop);
+
+ result.emplace_back(atom_b, p, symop.string());
+ }
else
result.emplace_back(atom_b);
}
diff --git a/src/MapMaker.cpp b/src/MapMaker.cpp
index 030a68f..9970d92 100644
--- a/src/MapMaker.cpp
+++ b/src/MapMaker.cpp
@@ -37,7 +37,6 @@
#include "pdb-redo/MapMaker.hpp"
#include "pdb-redo/ResolutionCalculator.hpp"
#include "pdb-redo/Statistics.hpp"
-#include "pdb-redo/Symmetry-2.hpp"
#ifdef _MSC_VER
#include <io.h>
@@ -411,33 +410,58 @@ void Map<FTYPE>::write_masked(const std::filesystem::path &f, clipper::Grid_rang
template <typename FTYPE>
Map<FTYPE> Map<FTYPE>::masked(const cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms) const
{
+ using clipper::Coord_frac;
+ using clipper::Coord_orth;
+ using clipper::Coord_map;
+ using clipper::Coord_grid;
+
Map<FTYPE> result(*this);
- auto rtops = AlternativeSites(getSpacegroup(structure.get_datablock()), getCell(structure.get_datablock()));
-
for (auto &atom : atoms)
{
float radius = cif::atom_type_traits(atom.get_type()).radius(cif::radius_type::van_der_waals);
if (std::isnan(radius))
radius = cif::atom_type_traits(atom.get_type()).radius(cif::radius_type::calculated);
+
if (std::isnan(radius)) // TODO: now what?
- radius = 200;
+ {
+ std::cerr << "Could not define radius for atom " << atom << std::endl;
+ continue;
+ }
- auto cloc = toClipper(atom.get_location());
+ float radiusSq = radius * radius;
- for (auto &rt : rtops)
- {
- auto rcloc = cloc.transform(rt);
+ auto o = Coord_orth(radius, radius, radius).coord_frac(mMap.cell());
+ o[0] = std::abs(o[0]);
+ o[1] = std::abs(o[1]);
+ o[2] = std::abs(o[2]);
+
+ auto pp = atom.get_location();
+ Coord_orth cloc = pp;
- iterateGrid(toClipper(toPoint(rcloc)), radius, result.mMap,
- [&result, radiusSq = radius * radius, a = atom.get_location()](auto iw)
+ Coord_frac fp = cloc.coord_frac(mMap.cell());
+ Coord_frac fMin = fp - o, fMax = fp + o;
+
+ // see if the box around p actually overlaps the cell
+
+ if (fMin.u() > 1 or fMax.u() < 0 or
+ fMin.v() > 1 or fMax.v() < 0 or
+ fMin.w() > 1 or fMax.w() < 0)
+ continue;
+
+ Coord_map mMin = fMin.coord_map(mMap.grid_sampling()), mMax = fMax.coord_map(mMap.grid_sampling());
+ Coord_grid gMin = mMin.floor(), gMax = mMax.ceil();
+
+ auto i0 = clipper::Xmap_base::Map_reference_coord(mMap, gMin);
+ for (auto iu = i0; iu.coord().u() <= gMax[0]; iu.next_u())
+ for (auto iv = iu; iv.coord().v() <= gMax[1]; iv.next_v())
+ for (auto iw = iv; iw.coord().w() <= gMax[2]; iw.next_w())
{
- cif::point p = toPoint(iw.coord_orth());
+ cif::point gp = iw.coord_orth();
- if (distance_squared(a, p) < radiusSq)
+ if (distance_squared(gp, pp) < radiusSq)
result.mMap[iw] = -10;
- });
- }
+ }
}
return result;
@@ -450,7 +474,7 @@ float Map<FTYPE>::z_weighted_density(const cif::mm::structure &structure, const
for (auto &atom : atoms)
{
- auto co = toClipper(atom.get_location());
+ clipper::Coord_orth co = atom.get_location();
auto a_cf = co.coord_frac(mMap.cell());
auto a_cm = a_cf.coord_map(mMap.grid_sampling());
diff --git a/src/Minimizer.cpp b/src/Minimizer.cpp
index 77a729e..45f8283 100644
--- a/src/Minimizer.cpp
+++ b/src/Minimizer.cpp
@@ -35,7 +35,6 @@
#include <regex>
#include "pdb-redo/Minimizer.hpp"
-#include "pdb-redo/Symmetry-2.hpp"
namespace fs = std::filesystem;
@@ -47,8 +46,14 @@ namespace pdb_redo
const uint32_t kRefSentinel = std::numeric_limits<uint32_t>::max();
const double
- kNonBondedContactDistanceSq = 15.0 * 15.0,
- kMaxPeptideBondLengthSq = 3.5 * 3.5;
+ kMaxNonBondedContactDistance = 10.0,
+ kMaxPeptideBondLength = 3.5,
+ kMaxPeptideBondLengthSq = kMaxPeptideBondLength * kMaxPeptideBondLength;
+
+const double
+ kDefaultMapWeight = 60,
+ kDefaultPlane5ESD = 0.11,
+ kDefaultChiralVolumeESD = 0.2;
// --------------------------------------------------------------------
@@ -73,15 +78,20 @@ std::string AtomLocationProvider::atom(AtomRef atomID) const
if (atomID >= mAtoms.size())
throw std::range_error("Unknown atom " + std::to_string(atomID));
auto &a = mAtoms[atomID];
- return std::to_string(a.get_label_seq_id()) + ' ' + a.get_label_atom_id();
+
+ auto seq_id = a.get_label_seq_id();
+
+ std::string symmetry;
+ if (a.symmetry() != "1_555")
+ symmetry = " " + a.symmetry();
+
+ return a.get_label_asym_id() + (seq_id ? std::to_string(seq_id) : a.get_auth_seq_id()) + ' ' + a.get_label_atom_id() + symmetry;
}
// --------------------------------------------------------------------
-Minimizer::Minimizer(const cif::mm::structure &structure, const BondMap &bonds, float plane5AtomsESD)
+Minimizer::Minimizer(const cif::mm::structure &structure)
: mStructure(structure)
- , mBonds(bonds)
- , mPlane5ESD(plane5AtomsESD)
{
}
@@ -295,7 +305,7 @@ void Minimizer::addPolySection(const cif::mm::polymer &poly, int first, int last
ref(r.get_atom_by_atom_id("CA"))
};
- mPlanarityRestraints.emplace_back(PlanarityRestraint{ move(atoms), mPlane5ESD });
+ mPlanarityRestraints.emplace_back(PlanarityRestraint{ move(atoms), kDefaultPlane5ESD });
}
}
catch (const std::exception &ex)
@@ -341,7 +351,7 @@ void Minimizer::addDensityMap(const XMap &xMap, float mapWeight)
mDensityRestraint.reset(new DensityRestraint(move(densityAtoms), xMap, mapWeight));
}
-void Minimizer::Finish()
+void Minimizer::Finish(const cif::crystal &crystal)
{
if (mAtoms.empty())
throw std::runtime_error("No atoms to refine");
@@ -355,9 +365,10 @@ void Minimizer::Finish()
std::set<std::tuple<AtomRef, AtomRef>> nbc;
- SymmetryAtomIteratorFactory saif(mStructure, getSpacegroup(mStructure.get_datablock()), getCell(mStructure.get_datablock()));
+ BondMap bm = createBondMap();
+ // const BondMap &bm = mBonds;
- auto add_nbc = [this, &nbc, &libAtom](const cif::mm::atom &a1, const cif::mm::atom &a2)
+ auto add_nbc = [this, &nbc, &libAtom, &bm](const cif::mm::atom &a1, const cif::mm::atom &a2)
{
AtomRef ra1 = ref(a1);
AtomRef ra2 = ref(a2);
@@ -398,7 +409,7 @@ void Minimizer::Finish()
double minDist = 2.8;
- if (mBonds.is1_4(a1, a2))
+ if (bm.is1_4(a1, a2))
{
if (cif::VERBOSE > 1)
std::cerr << "1_4 for " << a1 << " and " << a2 << std::endl;
@@ -504,17 +515,34 @@ void Minimizer::Finish()
};
// now add the non-bonded restraints
+
for (auto &a1 : mAtoms)
{
for (auto a2 : mStructure.atoms())
{
- if (a1 == a2 or mBonds(a1, a2))
+ if (a1 == a2)
continue;
- for (auto s_a2 : saif(a2, [l = a1.get_location()](const cif::point &p)
- { return distance_squared(p, l) <= kNonBondedContactDistanceSq; }))
+ if (distance_squared(a1, a2) < kMaxNonBondedContactDistance * kMaxNonBondedContactDistance)
+ {
+ if (not bm(a1, a2))
+ add_nbc(a1, a2);
+ // there used to be a continue here, but that's wrong of course
+ }
+
+ const auto &[d, p, symop] = crystal.closest_symmetry_copy(a1.get_location(), a2.get_location());
+
+ if (symop != cif::sym_op() and d < kMaxNonBondedContactDistance)
{
- add_nbc(a1, a2);
+ AtomRef ra1 = ref(a1);
+ AtomRef ra2 = ref(cif::mm::atom(a2, p, symop.string()));
+
+ if (not nbc.count(std::make_tuple(ra1, ra2)))
+ {
+ mNonBondedContactRestraints.emplace_back(ra1, ra2, 2.8, 0.02);
+ nbc.insert(std::make_tuple(ra1, ra2));
+ nbc.insert(std::make_tuple(ra2, ra1));
+ }
}
}
}
@@ -574,9 +602,29 @@ void Minimizer::dropTorsionRestraints()
mTorsionRestraints.clear();
}
+void Minimizer::setMapWeight(float mapWeight)
+{
+ mDensityRestraint->mMapWeight = mapWeight;
+}
+
+void Minimizer::setChiralVolumeESD(float chiralityESD)
+{
+ for (auto &r : mChiralVolumeRestraints)
+ r.mESD = chiralityESD;
+}
+
+void Minimizer::setPlanarityESD(float planarityESD)
+{
+ for (auto &r : mPlanarityRestraints)
+ r.mESD = planarityESD;
+}
+
AtomRef Minimizer::ref(const cif::mm::atom &atom)
{
std::string atomID = atom.id();
+ if (atom.is_symmetry_copy())
+ atomID += ':' + atom.symmetry();
+
AtomRef result;
auto k = mRefIndex.find(atomID);
@@ -598,16 +646,19 @@ void Minimizer::addLinkRestraints(const cif::mm::residue &a, const cif::mm::resi
auto c1 = cif::compound_factory::instance().create(a.get_compound_id());
auto c2 = cif::compound_factory::instance().create(b.get_compound_id());
- auto getCompoundAtom = [&](const LinkAtom &la)
- {
- return la.compID == 1 ? c1->get_atom_by_atom_id(la.atomID) : c2->get_atom_by_atom_id(la.atomID);
- };
-
assert(link.bonds().size() == 1);
bool a_is_1 = link.bonds().front().atom[0].compID == 1 ?
link.bonds().front().atom[0].atomID == atom_id_a :
link.bonds().front().atom[1].atomID == atom_id_a;
+ auto getCompoundAtom = [&](const LinkAtom &la)
+ {
+ if (la.compID == 1)
+ return a_is_1 ? c1->get_atom_by_atom_id(la.atomID) : c2->get_atom_by_atom_id(la.atomID);
+ else
+ return a_is_1 ? c2->get_atom_by_atom_id(la.atomID) : c1->get_atom_by_atom_id(la.atomID);
+ };
+
auto getAtom = [&](const LinkAtom &la)
{
if (la.compID == 1)
@@ -960,14 +1011,14 @@ void GSLDFCollector::add(AtomRef atom, double dx, double dy, double dz)
class GSLMinimizer : public Minimizer
{
public:
- GSLMinimizer(const cif::mm::structure &structure, const BondMap &bm, float plane5AtomsESD)
- : Minimizer(structure, bm, plane5AtomsESD)
+ GSLMinimizer(const cif::mm::structure &structure)
+ : Minimizer(structure)
{
}
- virtual void Finish()
+ virtual void Finish(const cif::crystal &crystal)
{
- Minimizer::Finish();
+ Minimizer::Finish(crystal);
for (auto &a : mReferencedAtoms)
mFixedLocations.push_back(a.get_location());
@@ -1007,7 +1058,8 @@ double GSLMinimizer::refine(bool storeAtoms)
fdf.n = mAtoms.size() * 3;
fdf.params = this;
- auto T = gsl_multimin_fdfminimizer_conjugate_pr;
+ // auto T = gsl_multimin_fdfminimizer_conjugate_pr;
+ auto T = gsl_multimin_fdfminimizer_vector_bfgs2;
auto x = gsl_vector_alloc(3 * mAtoms.size());
size_t ix = 0;
@@ -1021,8 +1073,10 @@ double GSLMinimizer::refine(bool storeAtoms)
m_s = gsl_multimin_fdfminimizer_alloc(T, 3 * mAtoms.size());
- float tolerance = 0.06f;
- double stepSize = 0.1 * gsl_blas_dnrm2(x);
+ // float tolerance = 0.06f;
+ // double stepSize = 0.25 * gsl_blas_dnrm2(x);
+ float tolerance = 0.1f;
+ double stepSize = 0.25;
gsl_multimin_fdfminimizer_set(m_s, &fdf, x, stepSize, tolerance);
@@ -1057,7 +1111,7 @@ double GSLMinimizer::refine(bool storeAtoms)
{
if (status != GSL_ENOPROG)
std::cerr << "Unexpected result from gsl_multimin_fdfminimizer_iterate: " << status << std::endl;
- else if (cif::VERBOSE > 0)
+ else if (cif::VERBOSE > 1)
std::cerr << "Minimizer stopped at iteration " << i << " at " << m_s->f << std::endl;
break;
}
@@ -1074,7 +1128,7 @@ double GSLMinimizer::refine(bool storeAtoms)
if (status == GSL_SUCCESS)
{
- if (cif::VERBOSE > 0)
+ if (cif::VERBOSE > 1)
std::cerr << "Minimum found at iteration " << i << " at " << m_s->f << std::endl;
break;
}
@@ -1174,20 +1228,19 @@ void GSLMinimizer::Fdf(const gsl_vector *x, double *f, gsl_vector *df)
// --------------------------------------------------------------------
-Minimizer *Minimizer::create(const cif::mm::polymer &poly, int first, int last, const BondMap &bonds,
- const XMap &xMap, float mapWeight, float plane5AtomsESD)
+Minimizer *Minimizer::create(const cif::crystal &crystal, const cif::mm::polymer &poly, int first, int last,
+ const XMap &xMap)
{
- std::unique_ptr<Minimizer> result(new GSLMinimizer(*poly.get_structure(), bonds, plane5AtomsESD));
+ std::unique_ptr<Minimizer> result(new GSLMinimizer(*poly.get_structure()));
result->addPolySection(poly, first, last);
- result->addDensityMap(xMap, mapWeight);
- result->Finish();
+ result->addDensityMap(xMap, kDefaultMapWeight);
+ result->Finish(crystal);
return result.release();
}
-Minimizer *Minimizer::create(cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms,
- const BondMap &bm, float plane5AtomsESD, const XMap *xMap, float mapWeight)
+Minimizer *Minimizer::create(const cif::crystal &crystal, cif::mm::structure &structure, const std::vector<cif::mm::atom> &atoms, const XMap *xMap)
{
- std::unique_ptr<Minimizer> result(new GSLMinimizer(structure, bm, plane5AtomsESD));
+ std::unique_ptr<Minimizer> result(new GSLMinimizer(structure));
std::vector<const cif::mm::residue *> residues;
@@ -1289,7 +1342,7 @@ Minimizer *Minimizer::create(cif::mm::structure &structure, const std::vector<ci
else
{
residues.emplace_back(&structure.get_residue(ptnr1_label_asym_id, ptnr1_label_seq_id, ptnr1_auth_seq_id));
- linked.emplace_back(rb, residues.back(), ptnr1_label_atom_id, ptnr2_label_atom_id, link_id);
+ linked.emplace_back(residues.back(), rb, ptnr1_label_atom_id, ptnr2_label_atom_id, link_id);
}
}
@@ -1305,19 +1358,52 @@ Minimizer *Minimizer::create(cif::mm::structure &structure, const std::vector<ci
try
{
result->addLinkRestraints(*a, *b, atom_a, atom_b, a->get_compound_id() + "-" + b->get_compound_id());
- }
- catch (const std::exception &e)
+ continue;
+ }
+ catch (...) {}
+
+ try
{
- result->addLinkRestraints(*b, *a, atom_a, atom_b, b->get_compound_id() + "-" + a->get_compound_id());
+ result->addLinkRestraints(*b, *a, atom_b, atom_a, b->get_compound_id() + "-" + a->get_compound_id());
+ continue;
}
+ catch (...) {}
+
+ // Last resort, if link is NAG-ASN, try pyr-ASN instead:
+
+ if (a->get_compound_id() == "NAG" and b->get_compound_id() == "ASN")
+ result->addLinkRestraints(*b, *a, atom_b, atom_a, "pyr-ASN");
+ else if (b->get_compound_id() == "NAG" and a->get_compound_id() == "ASN")
+ result->addLinkRestraints(*a, *b, atom_a, atom_b, "pyr-ASN");
+ else
+ throw std::runtime_error("Missing link information for " + a->get_compound_id() + " and " + b->get_compound_id());
}
- if (xMap != nullptr and mapWeight != 0)
- result->addDensityMap(*xMap, mapWeight);
+ if (xMap != nullptr)
+ result->addDensityMap(*xMap, kDefaultMapWeight);
- result->Finish();
+ result->Finish(crystal);
return result.release();
}
+BondMap Minimizer::createBondMap()
+{
+ std::vector<cif::point> pts;
+ for (auto a : mReferencedAtoms)
+ pts.emplace_back(a.get_location());
+
+ cif::point center = cif::centroid(pts);
+ float radius = 0;
+
+ for (auto pt : pts)
+ {
+ auto d = distance(pt, center);
+ if (radius < d)
+ radius = d;
+ }
+
+ return { mStructure.get_datablock(), std::make_tuple(center, radius + kMaxNonBondedContactDistance) };
+}
+
} // namespace pdb_redo
diff --git a/src/Restraints.cpp b/src/Restraints.cpp
index 63cda6b..6fa79c1 100644
--- a/src/Restraints.cpp
+++ b/src/Restraints.cpp
@@ -29,24 +29,18 @@
Date: dinsdag 22 mei, 2018
*/
-#include <numeric>
+#include "pdb-redo/Minimizer.hpp"
+#include "pdb-redo/Restraints.hpp"
#include <cif++.hpp>
-#include "pdb-redo/ClipperWrapper.hpp"
-#include "pdb-redo/Minimizer.hpp"
-#include "pdb-redo/Restraints.hpp"
+#include <Eigen/Eigenvalues>
+
+#include <numeric>
namespace pdb_redo
{
-using clipper::Coord_frac;
-using clipper::Coord_grid;
-using clipper::Coord_map;
-using clipper::Coord_orth;
-
-// --------------------------------------------------------------------
-
double BondRestraint::f(const AtomLocationProvider &atoms) const
{
double d = mDist - distance(atoms[mA], atoms[mB]);
@@ -292,15 +286,13 @@ void TorsionRestraint::print(const AtomLocationProvider &atoms) const
// --------------------------------------------------------------------
-const double kChiralVolumeESD = 0.2; // according to coot that's a reasonable value...
-
double ChiralVolumeRestraint::f(const AtomLocationProvider &atoms) const
{
auto chiralVolume = dot_product(atoms[mA1] - atoms[mCentre],
cross_product(atoms[mA2] - atoms[mCentre], atoms[mA3] - atoms[mCentre]));
double d = mVolume - chiralVolume;
- double result = (d * d) / (kChiralVolumeESD * kChiralVolumeESD);
+ double result = (d * d) / (mESD * mESD);
if (cif::VERBOSE > 2)
std::cerr << "chiral::f() = " << result << std::endl;
@@ -321,7 +313,7 @@ void ChiralVolumeRestraint::df(const AtomLocationProvider &atoms, DFCollector &d
auto chiralVolume = dot_product(a, cross_product(b, c));
auto d = chiralVolume - mVolume;
- auto s = 2 * d / (kChiralVolumeESD * kChiralVolumeESD);
+ auto s = 2 * d / (mESD * mESD);
df.add(mCentre, s * DPoint{
-(b.m_y * c.m_z - b.m_z * c.m_y) - (a.m_z * c.m_y - a.m_y * c.m_z) - (a.m_y * b.m_z - a.m_z * b.m_y),
@@ -335,7 +327,7 @@ void ChiralVolumeRestraint::df(const AtomLocationProvider &atoms, DFCollector &d
void ChiralVolumeRestraint::print(const AtomLocationProvider &atoms) const
{
- std::cout << "chiral volume " << atoms.atom(mA1) << " ; " << atoms.atom(mA2) << " ; " << atoms.atom(mA3) << " => " << mVolume << " / " << kChiralVolumeESD << std::endl;
+ std::cout << "chiral volume " << atoms.atom(mA1) << " ; " << atoms.atom(mA2) << " ; " << atoms.atom(mA3) << " => " << mVolume << " / " << mESD << std::endl;
}
// --------------------------------------------------------------------
@@ -348,26 +340,41 @@ void PlanarityRestraint::calculatePlaneFunction(const AtomLocationProvider &atom
center += atoms[a];
center /= mAtoms.size();
- clipper::Matrix<double> mat(3, 3);
+ double Cxx = 0, Cyy = 0, Czz = 0, Cxy = 0, Cxz = 0, Cyz = 0;
+
for (auto &a : mAtoms)
{
- mat(0, 0) += (atoms[a].m_x - center.m_x) * (atoms[a].m_x - center.m_x);
- mat(1, 1) += (atoms[a].m_y - center.m_y) * (atoms[a].m_y - center.m_y);
- mat(2, 2) += (atoms[a].m_z - center.m_z) * (atoms[a].m_z - center.m_z);
- mat(0, 1) += (atoms[a].m_x - center.m_x) * (atoms[a].m_y - center.m_y);
- mat(0, 2) += (atoms[a].m_x - center.m_x) * (atoms[a].m_z - center.m_z);
- mat(1, 2) += (atoms[a].m_y - center.m_y) * (atoms[a].m_z - center.m_z);
+ Cxx += (atoms[a].m_x - center.m_x) * (atoms[a].m_x - center.m_x);
+ Cyy += (atoms[a].m_y - center.m_y) * (atoms[a].m_y - center.m_y);
+ Czz += (atoms[a].m_z - center.m_z) * (atoms[a].m_z - center.m_z);
+ Cxy += (atoms[a].m_x - center.m_x) * (atoms[a].m_y - center.m_y);
+ Cxz += (atoms[a].m_x - center.m_x) * (atoms[a].m_z - center.m_z);
+ Cyz += (atoms[a].m_y - center.m_y) * (atoms[a].m_z - center.m_z);
}
+
+ Eigen::Matrix3d mat;
+ mat << Cxx, Cxy, Cxz,
+ Cxy, Cyy, Cyz,
+ Cxz, Cyz, Czz;
+
+ Eigen::EigenSolver<Eigen::Matrix3d> es(mat);
+
+ auto ev = es.eigenvalues();
+
+ float b_ev = std::numeric_limits<float>::max();
+ for (size_t i = 0; i < 3; ++i)
+ {
+ if (ev[i].real() > b_ev)
+ continue;
- mat(1, 0) = mat(0, 1);
- mat(2, 0) = mat(0, 2);
- mat(2, 1) = mat(1, 2);
+ b_ev = ev[i].real();
- /*auto eigen = */mat.eigen(true);
+ auto col = es.eigenvectors().col(i);
- abcd[0] = mat(0, 0);
- abcd[1] = mat(1, 0);
- abcd[2] = mat(2, 0);
+ abcd[0] = col(0).real();
+ abcd[1] = col(1).real();
+ abcd[2] = col(2).real();
+ }
double sumSq = 1e-20 + abcd[0] * abcd[0] + abcd[1] * abcd[1] + abcd[2] * abcd[2];
@@ -464,12 +471,13 @@ double NonBondedContactRestraint::f(const AtomLocationProvider &atoms) const
{
double d = mMinDist - std::sqrt(distance);
result = (d * d) / (mDistESD * mDistESD);
+
+ if (cif::VERBOSE > 2)
+ std::cerr << "non-bonded-contact::f() = " << result << " min-dist is " << mMinDist << " and dist is " << std::sqrt(distance)
+ << " a1: " << atoms.atom(mA) << " a2: " << atoms.atom(mB) << std::endl
+ << " a1: " << atoms[mA] << " a2: " << atoms[mB] << std::endl;
}
- if (cif::VERBOSE > 2)
- std::cerr << "non-bonded-contact::f() = " << result << " min-dist is " << mMinDist << " and dist is " << std::sqrt(distance)
- << " a1: " << atoms.atom(mA) << " a2: " << atoms.atom(mB) << std::endl
- << " a1: " << atoms[mA] << " a2: " << atoms[mB] << std::endl;
return result;
}
@@ -519,8 +527,8 @@ double DensityRestraint::f(const AtomLocationProvider &atoms) const
for (auto &a : mAtoms)
{
- Coord_orth p = toClipper(atoms[a.first]);
- Coord_frac pf = p.coord_frac(mXMap.cell());
+ clipper::Coord_orth p{atoms[a.first].m_x, atoms[a.first].m_y, atoms[a.first].m_z };
+ clipper::Coord_frac pf = p.coord_frac(mXMap.cell());
result += a.second * mXMap.interp<clipper::Interp_cubic>(pf);
}
@@ -538,8 +546,8 @@ void DensityRestraint::df(const AtomLocationProvider &atoms, DFCollector &df) co
for (auto &a : mAtoms)
{
- Coord_orth p = toClipper(atoms[a.first]);
- Coord_frac pf = p.coord_frac(mXMap.cell());
+ clipper::Coord_orth p{atoms[a.first].m_x, atoms[a.first].m_y, atoms[a.first].m_z };
+ clipper::Coord_frac pf = p.coord_frac(mXMap.cell());
auto pm = pf.coord_map(mXMap.grid_sampling());
clipper::Grad_map<double> grad;
diff --git a/src/Statistics.cpp b/src/Statistics.cpp
index e290cdb..f47d614 100644
--- a/src/Statistics.cpp
+++ b/src/Statistics.cpp
@@ -32,6 +32,7 @@
#include "pdb-redo/AtomShape.hpp"
#include "pdb-redo/BondMap.hpp"
#include "pdb-redo/ClipperWrapper.hpp"
+#include "pdb-redo/DistanceMap.hpp"
#include "pdb-redo/Statistics.hpp"
// --------------------------------------------------------------------
@@ -39,10 +40,6 @@
namespace pdb_redo
{
-using clipper::Coord_grid;
-using clipper::Coord_orth;
-using clipper::Xmap;
-
using cif::atom_type_traits;
// --------------------------------------------------------------------
@@ -233,13 +230,13 @@ class PointWeightFunction
struct AtomGridData
{
- AtomGridData(const Coord_grid &gp, double density)
+ AtomGridData(const clipper::Coord_grid &gp, double density)
: p(gp)
, density(density)
{
}
- Coord_grid p;
+ clipper::Coord_grid p;
double density;
};
@@ -311,7 +308,7 @@ struct AtomData
// --------------------------------------------------------------------
-std::tuple<float, float> CalculateMapStatistics(const Xmap<float> &f)
+std::tuple<float, float> CalculateMapStatistics(const clipper::Xmap<float> &f)
{
double sum = 0, sum2 = 0;
int count = 0;
@@ -548,7 +545,7 @@ void StatsCollector::initialize()
std::vector<ResidueStatistics> StatsCollector::collect() const
{
- std::vector<std::tuple<std::string, int, std::string>> residues;
+ residue_list residues;
std::vector<cif::mm::atom> atoms;
for (auto atom : mStructure.atoms())
@@ -569,6 +566,29 @@ std::vector<ResidueStatistics> StatsCollector::collect() const
return collect(residues, bbox, true);
}
+std::vector<ResidueStatistics> StatsCollector::collect(const std::string &asymID) const
+{
+ residue_list residues;
+ std::vector<cif::mm::atom> atoms;
+
+ for (auto atom_id : mStructure.get_datablock()["atom_site"].find<std::string>(cif::key("label_asym_id") == asymID, "id"))
+ {
+ auto &atom = atoms.emplace_back(mStructure.get_atom_by_id(atom_id));
+
+ auto k = std::make_tuple(atom.get_label_asym_id(), atom.get_label_seq_id(), atom.get_auth_seq_id());
+
+ if (residues.empty() or residues.back() != k)
+ {
+ residues.emplace_back(move(k));
+ atoms.emplace_back(std::move(atom));
+ }
+ }
+
+ BoundingBox bbox(mStructure, atoms, 5.0f);
+ return collect(residues, bbox, false);
+}
+
+
std::vector<ResidueStatistics> StatsCollector::collect(const std::string &asymID, int resFirst, int resLast, bool authNameSpace) const
{
residue_list residues;
@@ -951,13 +971,13 @@ void StatsCollector::sumDensity(std::vector<AtomData> &atomData,
auto radius = data.radius;
double sumDensity = 0;
- iterateGrid(toClipper(atom.get_location()), radius, Fb, [&](Xmap_base::Map_reference_coord &iw)
+ iterateGrid(atom.get_location(), radius, Fb, [&, radius_sq = radius * radius](Xmap_base::Map_reference_coord &iw)
{
- auto p = toPoint(iw.coord_orth());
+ cif::point p = iw.coord_orth();
- double d = distance(p, atom.get_location());
+ double d = distance_squared(p, atom.get_location());
- if (d <= radius)
+ if (d <= radius_sq)
{
double density = shape.calculatedDensity(p);
@@ -983,12 +1003,17 @@ void StatsCollector::collectSums(std::vector<AtomData> &atomData, GridPtDataMap
const Xmap<float> &Fb = mMapMaker.fb();
const Xmap<float> &Fd = mMapMaker.fd();
- cif::Progress progress(atomData.size(), "Stats calculation");
+ cif::progress_bar progress_bar(atomData.size(), "Stats calculation");
// Iterate over the atom data to collect the sums
for (auto &d : atomData)
{
- auto rmsScaledF = mRmsScaled.at(d.asymID);
+ auto rmsi = mRmsScaled.find(d.asymID);
+
+ if (rmsi == mRmsScaled.end())
+ continue;
+
+ auto rmsScaledF = rmsi->second;
for (auto gp : d.points)
{
@@ -1001,7 +1026,7 @@ void StatsCollector::collectSums(std::vector<AtomData> &atomData, GridPtDataMap
double e = gp.density / gpd;
double t = e * mSZ / rmsScaledF.second;
- Xmap_base::Map_reference_coord ix(Fb, gp.p);
+ clipper::Xmap_base::Map_reference_coord ix(Fb, gp.p);
double fb = Fb[ix];
double fd = Fd[ix];
@@ -1037,7 +1062,7 @@ void StatsCollector::collectSums(std::vector<AtomData> &atomData, GridPtDataMap
d.sums.edSums[1] += ed3;
}
- progress.consumed(1);
+ progress_bar.consumed(1);
}
}
@@ -1052,11 +1077,8 @@ void StatsCollector::calculate(std::vector<AtomData> &atomData) const
// --------------------------------------------------------------------
-EDIAStatsCollector::EDIAStatsCollector(MapMaker<float> &mm,
- cif::mm::structure &structure, bool electronScattering, const BondMap &bondMap)
+EDIAStatsCollector::EDIAStatsCollector(const MapMaker<float> &mm, cif::mm::structure &structure, bool electronScattering)
: StatsCollector(mm, structure, electronScattering)
- , mDistanceMap(structure, mm.spacegroup(), mm.cell(), 3.5f)
- , mBondMap(bondMap)
{
// create a atom radius map, for EDIA
@@ -1108,7 +1130,7 @@ void EDIAStatsCollector::calculate(std::vector<AtomData> &atomData) const
{
StatsCollector::calculate(atomData);
- const Xmap<float> &Fb = mMapMaker.fb();
+ const clipper::Xmap<float> &Fb = mMapMaker.fb();
// Xmap<float>& fd = mMapMaker.fd();
struct lessAtom
@@ -1120,21 +1142,19 @@ void EDIAStatsCollector::calculate(std::vector<AtomData> &atomData) const
// Calculate EDIA scores
- cif::Progress progress(atomData.size(), "EDIA calculation");
+ DistanceMap dm(mStructure, 3.5f);
+ BondMap bm = createBondMap(atomData);
+
+ cif::progress_bar progress_bar(atomData.size(), "EDIA calculation");
for (auto &data : atomData)
{
auto &atom = data.atom;
float radius = mRadii.at(atom.get_type());
- // if (cif::VERBOSE > 2)
- // std::cerr << (atomData.size() + 1) << '\t'
- // << atom_type_traits(atom.get_type()).symbol() << '\t'
- // << radius << std::endl;
-
PointWeightFunction w(atom.get_location(), radius);
- std::vector<cif::mm::atom> atomsNearBy = mDistanceMap.near(atom, 3.5f);
+ std::vector<cif::mm::atom> atomsNearBy = dm.near(atom, 3.5f);
std::vector<PointWeightFunction> wn;
for (auto a : atomsNearBy)
@@ -1142,9 +1162,9 @@ void EDIAStatsCollector::calculate(std::vector<AtomData> &atomData) const
float ediaSum[2] = {};
- iterateGrid(toClipper(atom.get_location()), radius, Fb, [&](auto iw)
+ iterateGrid(atom.get_location(), radius, Fb, [&](auto iw)
{
- cif::point p = toPoint(iw.coord_orth());
+ cif::point p = iw.coord_orth();
// EDIA calculations
auto fb = Fb[iw];
@@ -1185,7 +1205,7 @@ void EDIAStatsCollector::calculate(std::vector<AtomData> &atomData) const
{
S.insert(atomsNearBy[i]);
- if (not mBondMap(atomsNearBy[i], atom))
+ if (not bm(atomsNearBy[i], atom))
I.insert(atomsNearBy[i]);
}
}
@@ -1231,8 +1251,27 @@ void EDIAStatsCollector::calculate(std::vector<AtomData> &atomData) const
if (data.edia < 0)
data.edia = 0;
- progress.consumed(1);
+ progress_bar.consumed(1);
+ }
+}
+
+BondMap EDIAStatsCollector::createBondMap(std::vector<AtomData> &atomData) const
+{
+ std::vector<cif::point> pts;
+ for (auto a : atomData)
+ pts.emplace_back(a.atom.get_location());
+
+ cif::point center = cif::centroid(pts);
+ float radius = 0;
+
+ for (auto pt : pts)
+ {
+ auto d = distance(pt, center);
+ if (radius < d)
+ radius = d;
}
+
+ return { mStructure.get_datablock(), std::make_tuple(center, radius + 3.5f) };
}
} // namespace pdb_redo
diff --git a/src/Symmetry-2.cpp b/src/Symmetry-2.cpp
deleted file mode 100644
index 44523a2..0000000
--- a/src/Symmetry-2.cpp
+++ /dev/null
@@ -1,366 +0,0 @@
-/*-
- * 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
- * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <atomic>
-#include <mutex>
-
-#include "cif++.hpp"
-
-#include "pdb-redo/ClipperWrapper.hpp"
-#include "pdb-redo/Symmetry-2.hpp"
-
-namespace c = cif;
-
-namespace pdb_redo
-{
-
-// --------------------------------------------------------------------
-
-clipper::Coord_orth CalculateOffsetForCell(const cif::mm::structure &structure, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell)
-{
- using namespace cif::literals;
-
- auto &atoms = structure.atoms();
- size_t dim = atoms.size();
-
- std::vector<clipper::Coord_orth> locations;
- locations.reserve(dim);
-
- // bounding box
- cif::point pMin(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()),
- pMax(std::numeric_limits<float>::min(), std::numeric_limits<float>::min(), std::numeric_limits<float>::min());
-
- for (auto &atom : atoms)
- {
- cif::point pt = atom.get_location();
- locations.push_back(toClipper(pt));
-
- if (pMin.m_x > pt.m_x)
- pMin.m_x = pt.m_x;
- if (pMin.m_y > pt.m_y)
- pMin.m_y = pt.m_y;
- if (pMin.m_z > pt.m_z)
- pMin.m_z = pt.m_z;
-
- if (pMax.m_x < pt.m_x)
- pMax.m_x = pt.m_x;
- if (pMax.m_y < pt.m_y)
- pMax.m_y = pt.m_y;
- if (pMax.m_z < pt.m_z)
- pMax.m_z = pt.m_z;
- };
-
- // correct locations so that the median of x, y and z are inside the cell
- std::vector<float> c(dim);
- auto median = [&]()
- {
- return dim % 1 == 0
- ? c[dim / 2]
- : (c[dim / 2 - 1] + c[dim / 2]) / 2;
- };
-
- std::transform(locations.begin(), locations.end(), c.begin(), [](auto &l)
- { return static_cast<float>(l[0]); });
- std::sort(c.begin(), c.end());
- float mx = median();
-
- std::transform(locations.begin(), locations.end(), c.begin(), [](auto &l)
- { return static_cast<float>(l[1]); });
- std::sort(c.begin(), c.end());
- float my = median();
-
- std::transform(locations.begin(), locations.end(), c.begin(), [](auto &l)
- { return static_cast<float>(l[2]); });
- std::sort(c.begin(), c.end());
- float mz = median();
-
- if (cif::VERBOSE > 1)
- std::cerr << "median position of atoms: " << cif::point(mx, my, mz) << std::endl;
-
- auto calculateD = [&](float m, float c)
- {
- float d = 0;
- assert(c != 0);
- if (c != 0)
- {
- while (m + d < -(c / 2))
- d += c;
- while (m + d > (c / 2))
- d -= c;
- }
- return d;
- };
-
- if (cell.a() == 0 or cell.b() == 0 or cell.c() == 0)
- throw std::runtime_error("Invalid cell, contains a dimension that is zero");
-
- cif::point D;
-
- D.m_x = calculateD(mx, static_cast<float>(cell.a()));
- D.m_y = calculateD(my, static_cast<float>(cell.b()));
- D.m_z = calculateD(mz, static_cast<float>(cell.c()));
-
- if (D.m_x != 0 or D.m_y != 0 or D.m_z != 0)
- {
- if (cif::VERBOSE > 1)
- std::cerr << "moving coorinates by " << D.m_x << ", " << D.m_y << " and " << D.m_z << std::endl;
- }
-
- return toClipper(D);
-}
-
-// --------------------------------------------------------------------
-
-std::vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup &spacegroup,
- const clipper::Cell &cell)
-{
- std::vector<clipper::RTop_orth> result;
-
- // to make the operation at index 0 equal to identity
- result.push_back(clipper::RTop_orth::identity());
-
- for (int i = 0; i < spacegroup.num_symops(); ++i)
- {
- const auto &symop = spacegroup.symop(i);
-
- for (int u : {-1, 0, 1})
- for (int v : {-1, 0, 1})
- for (int w : {-1, 0, 1})
- {
- if (i == 0 and u == 0 and v == 0 and w == 0)
- continue;
-
- auto rtop = clipper::RTop_frac(
- symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w))
- .rtop_orth(cell);
-
- result.push_back(std::move(rtop));
- }
- }
-
- return result;
-}
-
-// --------------------------------------------------------------------
-// Unfortunately, clipper has a different numbering scheme than PDB
-// for rotation numbers. So we created a table to map those.
-// Perhaps a bit over the top, but hey....
-
-int32_t GetRotationalIndexNumber(int spacegroup, const clipper::RTop_frac &rt)
-{
- auto k = GetSymOpDataForRTop_frac(rt);
-
- const size_t N = cif::kSymopNrTableSize;
- int32_t L = 0, R = static_cast<int32_t>(N - 1);
- while (L <= R)
- {
- int32_t i = (L + R) / 2;
- if (cif::kSymopNrTable[i].spacegroup() < spacegroup)
- L = i + 1;
- else
- R = i - 1;
- }
-
- for (size_t i = L; i < N and cif::kSymopNrTable[i].spacegroup() == spacegroup; ++i)
- {
- if (cif::kSymopNrTable[i].symop() == k)
- return cif::kSymopNrTable[i].rotational_number();
- }
-
- throw std::runtime_error("Symmetry operation was not found in table, cannot find rotational number");
-}
-
-// -----------------------------------------------------------------------
-
-std::string SpacegroupToHall(std::string spacegroup)
-{
- int nr = cif::get_space_group_number(spacegroup);
-
- // yeah, sucks, I know, might be looping three times this way
- std::string result;
- for (size_t i = 0; i < cif::kNrOfSpaceGroups; ++i)
- {
- auto &sp = cif::kSpaceGroups[i];
- if (sp.nr == nr)
- {
- result = sp.Hall;
- break;
- }
- }
-
- if (result.empty())
- throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
-
- return result;
-}
-
-clipper::Spgr_descr GetCCP4SpacegroupDescr(int nr)
-{
- for (size_t i = 0; i < cif::kNrOfSpaceGroups; ++i)
- {
- auto &sg = cif::kSpaceGroups[i];
- if (sg.nr == nr)
- return clipper::Spgr_descr(sg.Hall, clipper::Spgr_descr::Hall);
- }
-
- throw std::runtime_error("Invalid spacegroup number: " + std::to_string(nr));
-}
-
-// --------------------------------------------------------------------
-
-std::string describeRToperation(const clipper::Spacegroup &spacegroup, const clipper::Cell &cell, const clipper::RTop_orth &rt)
-{
- auto spacegroup_nr = getSpacegroupNumber(spacegroup);
-
- if (not(rt.is_null() or rt.equals(clipper::RTop_orth::identity(), 0.0001f)))
- {
- for (int i = 0; i < spacegroup.num_symops(); ++i)
- {
- const auto &symop = spacegroup.symop(i);
-
- for (int u : {-1, 0, 1})
- for (int v : {-1, 0, 1})
- for (int w : {-1, 0, 1})
- {
- // if (i == 0 and u == 0 and v == 0 and w == 0)
- // continue;
-
- auto rtop = clipper::RTop_frac(
- symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w))
- .rtop_orth(cell);
-
- if (rtop.rot().equals(rt.rot(), 0.00001) and rtop.trn().equals(rt.trn(), 0.000001))
- {
- // gotcha
-
- auto rtop_f = rtop.rtop_frac(cell);
-
- int rnr = GetRotationalIndexNumber(spacegroup_nr, rtop_f);
-
- uint32_t t[3] = {
- static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),
- static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[1]))),
- static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[2])))};
-
- if (t[0] > 9 or t[1] > 9 or t[2] > 9)
- throw std::runtime_error("Symmetry operation has an out-of-range translation.");
-
- return std::to_string(rnr) + "_" + std::to_string(t[0]) + std::to_string(t[1]) + std::to_string(t[2]);
- }
- }
- }
- }
-
- return "1_555";
-}
-
-// --------------------------------------------------------------------
-
-cif::mm::atom symmetryCopy(const cif::mm::atom &atom, const cif::point &d,
- const clipper::Spacegroup &spacegroup, const clipper::Cell &cell, const clipper::RTop_orth &rt)
-{
- auto loc = atom.get_location();
- loc += d;
- loc = toPoint(toClipper(loc).transform(rt));
- loc -= d;
-
- std::string rt_operation = describeRToperation(spacegroup, cell, rt);
-
- return cif::mm::atom(atom, loc, rt_operation);
-}
-
-// -----------------------------------------------------------------------
-
-SymmetryAtomIteratorFactory::SymmetryAtomIteratorFactory(const cif::mm::structure &p, const clipper::Spacegroup &spacegroup, const clipper::Cell &cell)
- : mSpacegroup(spacegroup)
- , mD(toPoint(CalculateOffsetForCell(p, mSpacegroup, cell)))
- , mRtOrth(AlternativeSites(mSpacegroup, cell))
- , mCell(cell)
-{
-}
-
-SymmetryAtomIteratorFactory::SymmetryAtomIteratorFactory(const cif::mm::structure &p, int spacegroupNr, const clipper::Cell &cell)
- : SymmetryAtomIteratorFactory(p, clipper::Spacegroup{GetCCP4SpacegroupDescr(spacegroupNr)}, cell)
-{
-}
-
-// string SymmetryAtomIteratorFactory::symop_mmcif(const cif::cif::mm::atom& a) const
-// {
-// string result;
-
-// if (not a.isSymmetryCopy())
-// result = "1_555";
-// else
-// {
-// auto rtop_o = a.symop();
-
-// for (int i = 0; i < mSpacegroup.num_symops(); ++i)
-// {
-// const auto& symop = mSpacegroup.symop(i);
-
-// for (int u: { -1, 0, 1})
-// for (int v: { -1, 0, 1})
-// for (int w: { -1, 0, 1})
-// {
-// // if (i == 0 and u == 0 and v == 0 and w == 0)
-// // continue;
-
-// auto rtop = clipper::RTop_frac(
-// symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
-// ).rtop_orth(mCell);
-
-// if (rtop.rot().equals(rtop_o.rot(), 0.00001) and rtop.trn().equals(rtop_o.trn(), 0.000001))
-// {
-// // gotcha
-
-// auto rtop_f = rtop.rtop_frac(mCell);
-
-// int rnr = GetRotationalIndexNumber(mSpacegroupNr, rtop_f);
-
-// uint32_t t[3] = {
-// static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),
-// static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[1]))),
-// static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[2])))
-// };
-
-// if (t[0] > 9 or t[1] > 9 or t[2] > 9)
-// throw runtime_error("Symmetry operation has an out-of-range translation.");
-
-// result += to_string(rnr) + "_"
-// + to_string(t[0])
-// + to_string(t[1])
-// + to_string(t[2]);
-
-// return result;
-// }
-// }
-// }
-// }
-
-// return result;
-// }
-
-} // namespace pdb_redo
\ No newline at end of file
diff --git a/test/pdb-redo-unit-test.cpp b/test/pdb-redo-unit-test.cpp
index 08623c4..f07e7a5 100644
--- a/test/pdb-redo-unit-test.cpp
+++ b/test/pdb-redo-unit-test.cpp
@@ -27,6 +27,13 @@
#define BOOST_TEST_ALTERNATIVE_INIT_API
#include <boost/test/included/unit_test.hpp>
+#include "pdb-redo/AtomShape.hpp"
+#include "pdb-redo/ClipperWrapper.hpp"
+#include "pdb-redo/DistanceMap.hpp"
+#include "pdb-redo/MapMaker.hpp"
+#include "pdb-redo/Statistics.hpp"
+#include "pdb-redo/SkipList.hpp"
+
#include <clipper/clipper-ccp4.h>
#include <clipper/clipper-contrib.h>
@@ -36,11 +43,6 @@
#include <cif++.hpp>
#include <cif++/text.hpp>
-#include "pdb-redo/AtomShape.hpp"
-#include "pdb-redo/MapMaker.hpp"
-#include "pdb-redo/Statistics.hpp"
-#include "pdb-redo/SkipList.hpp"
-
namespace fs = std::filesystem;
namespace tt = boost::test_tools;
namespace utf = boost::unit_test;
@@ -66,35 +68,6 @@ bool init_unit_test()
return true;
}
-// --------------------------------------------------------------------
-
-BOOST_AUTO_TEST_CASE(symm_1)
-{
- if (not fs::exists(gTestDir / "1mdn.mtz"))
- return;
-
- const std::string spacegroup = "I 1 21 1";
-
- int nr = cif::get_space_group_number(spacegroup);
- BOOST_TEST(nr == 5005);
-
- // clipper::Spacegroup sgr(clipper::Spgr_descr{ spacegroup, clipper::Spgr_descr::HM });
-
- // BOOST_TEST(sgr.symbol_hm() == spacegroup);
-
- using clipper::CCP4MTZfile;
-
- CCP4MTZfile mtzin;
- mtzin.open_read((gTestDir / "1mdn.mtz").string());
-
- HKL_info mHKLInfo;
- mtzin.import_hkl_info(mHKLInfo);
-
- nr = pdb_redo::getSpacegroupNumber(mHKLInfo.spacegroup());
-
- BOOST_TEST(nr == 5005);
-}
-
// --------------------------------------------------------------------
// skip list test
@@ -346,9 +319,7 @@ BOOST_AUTO_TEST_CASE(stats_1)
float samplingRate = 1.5;
mm.loadMTZ(gTestDir / ".." / "examples" / "1cbs_map.mtz", samplingRate);
- BondMap bm(structure);
-
- pdb_redo::EDIAStatsCollector collector(mm, structure, false, bm);
+ pdb_redo::EDIAStatsCollector collector(mm, structure, false);
auto r = collector.collect();
auto ti = test.begin();
@@ -423,9 +394,7 @@ BOOST_AUTO_TEST_CASE(stats_2)
float samplingRate = 0.75;
mm.loadMTZ(gTestDir / ".." / "examples" / "1cbs_map.mtz", samplingRate);
- BondMap bm(structure);
-
- pdb_redo::EDIAStatsCollector collector(mm, structure, false, bm);
+ pdb_redo::EDIAStatsCollector collector(mm, structure, false);
auto r = collector.collect();
for (auto& ri: r)
@@ -444,3 +413,61 @@ BOOST_AUTO_TEST_CASE(stats_2)
}
}
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(bond_map_1)
+{
+ const fs::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+ cif::file file = cif::pdb::read(example.string());
+
+ BOOST_CHECK(file.is_valid());
+
+ auto &db = file.front();
+ auto &atom_site = db["atom_site"];
+
+ std::vector<std::string> atom_ids;
+ std::vector<cif::point> atom_locs;
+ for (const auto &[id, x, y, z] : atom_site.rows<std::string,float,float,float>("id", "cartn_x", "cartn_y", "cartn_z"))
+ {
+ atom_ids.emplace_back(id);
+ atom_locs.emplace_back(x, y, z);
+ }
+
+ BondMap bm(db);
+
+ using key_type = std::tuple<std::string,std::string>;
+ std::map<key_type,bool> bonded;
+
+ for (size_t i = 0; i + 1 < atom_ids.size(); ++i)
+ {
+ auto a = atom_ids[i];
+
+ for (size_t j = i + 1; j < atom_ids.size(); ++j)
+ {
+ auto b = atom_ids[j];
+
+ bonded.emplace(std::make_tuple(a, b), bm(a, b));
+ }
+ }
+
+ cif::point c = atom_locs.front();
+ BondMap bm2(db, std::make_tuple(c, 6.0f));
+
+ for (size_t i = 0; i + 1 < atom_ids.size(); ++i)
+ {
+ auto a = atom_ids[i];
+ auto pa = atom_locs[i];
+
+ for (size_t j = i + 1; j < atom_ids.size(); ++j)
+ {
+ auto b = atom_ids[j];
+ auto pb = atom_locs[j];
+
+ if (distance(pa, c) < 6 and distance(pb, c) < 6)
+ BOOST_TEST(bm2(a, b) == bonded[std::make_tuple(a, b)]);
+ else
+ // BOOST_CHECK_THROW(bm2(a, b), std::out_of_range);
+ BOOST_TEST(bm2(a, b) == false);
+ }
+ }
+}
\ No newline at end of file
diff --git a/test/rsr-test.cpp b/test/rsr-test.cpp
index 9f003c6..037fdda 100644
--- a/test/rsr-test.cpp
+++ b/test/rsr-test.cpp
@@ -43,6 +43,7 @@ BOOST_AUTO_TEST_CASE(refine_0)
cif::file file(example.string());
cif::mm::structure structure(file);
+ cif::crystal crystal(structure.get_datablock());
auto &chain = structure.polymers().front();
@@ -50,9 +51,7 @@ BOOST_AUTO_TEST_CASE(refine_0)
float samplingRate = 0.75;
mm.loadMTZ(gTestDir / ".." / "examples" / "1cbs_map.mtz", samplingRate);
- pdb_redo::BondMap bonds(structure);
-
- auto minimizer = pdb_redo::Minimizer::create(chain, 3, 3, bonds, mm.fb());
+ auto minimizer = pdb_redo::Minimizer::create(crystal, chain, 3, 3, mm.fb());
minimizer->printStats();
@@ -89,7 +88,7 @@ BOOST_AUTO_TEST_CASE(refine_0)
auto rmsd = std::sqrt(d_sum / atoms3.size());
std::cout << "RMSd: " << rmsd << std::endl;
- BOOST_CHECK(rmsd < 0.2);
+ BOOST_CHECK(rmsd < 0.25);
std::cout << std::string(cif::get_terminal_width(), '=') << std::endl;
}
@@ -100,8 +99,7 @@ BOOST_AUTO_TEST_CASE(refine_1)
cif::file file(example.string());
cif::mm::structure structure(file);
-
- pdb_redo::BondMap bonds(structure);
+ cif::crystal crystal(structure.get_datablock());
pdb_redo::MapMaker<float> mm;
float samplingRate = 0.75;
@@ -110,7 +108,7 @@ BOOST_AUTO_TEST_CASE(refine_1)
auto &chain = structure.polymers().front();
auto &atoms3 = chain.at(2).atoms();
- auto minimizer = pdb_redo::Minimizer::create(chain, 3, 3, bonds, mm.fb());
+ auto minimizer = pdb_redo::Minimizer::create(crystal, chain, 3, 3, mm.fb());
minimizer->printStats();
@@ -164,108 +162,109 @@ BOOST_AUTO_TEST_CASE(refine_1)
auto rmsd = std::sqrt(d_sum / atoms3.size());
std::cout << "RMSd: " << rmsd << std::endl;
- BOOST_CHECK(rmsd < 0.2);
+ BOOST_CHECK(rmsd < 0.27);
std::cout << std::string(cif::get_terminal_width(), '=') << std::endl;
}
-// BOOST_AUTO_TEST_CASE(refine_2)
-// {
-// const fs::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
-// cif::file file(example.string());
+BOOST_AUTO_TEST_CASE(refine_2)
+{
+ const fs::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
+ cif::file file(example.string());
-// cif::mm::structure structure(file);
+ cif::mm::structure structure(file);
-// const float kNearBy = 3;
+ const float kNearBy = 3;
-// pdb_redo::DistanceMap dm(structure, kNearBy);
+ pdb_redo::DistanceMap dm(structure, kNearBy);
-// // Move the REA residue somewhat
+ // Move the REA residue somewhat
-// auto &rea = structure.get_residue("B");
+ auto &rea = structure.get_residue("B");
-// std::vector<cif::mm::atom> atoms;
-// for (auto a : rea.atoms())
-// {
-// for (auto b : dm.near(a, kNearBy))
-// {
-// if (find(atoms.begin(), atoms.end(), b) != atoms.end())
-// continue;
+ std::vector<cif::mm::atom> atoms;
+ for (auto a : rea.atoms())
+ {
+ for (auto b : dm.near(a, kNearBy))
+ {
+ if (find(atoms.begin(), atoms.end(), b) != atoms.end())
+ continue;
-// atoms.push_back(b);
-// }
-// }
+ atoms.push_back(b);
+ }
+ }
-// // // translate by { 0.1, 0.1, 0.1 } and then
-// // // rotate around 1, 0, 0 for 5 degrees
+ // // translate by { 0.1, 0.1, 0.1 } and then
+ // // rotate around 1, 0, 0 for 5 degrees
-// // const float angle = 5 * (cif::kPI / 180);
-// // cif::quaternion q(
-// // std::cos(angle / 2), std::sin(angle / 2), 0, 0
-// // );
+ // const float angle = 5 * (cif::kPI / 180);
+ // cif::quaternion q(
+ // std::cos(angle / 2), std::sin(angle / 2), 0, 0
+ // );
-// // for (auto a : rea.atoms())
-// // a.translate_and_rotate({ 0.1, 0.1, 0.1 }, q);
+ // for (auto a : rea.atoms())
+ // a.translate_and_rotate({ 0.1, 0.1, 0.1 }, q);
-// for (auto a : rea.atoms())
-// {
-// auto l = a.get_location();
-// l = nudge(l, 0.5);
-// a.set_location(l);
-// }
+ for (auto a : rea.atoms())
+ {
+ auto l = a.get_location();
+ l = nudge(l, 0.5);
+ a.set_location(l);
+ }
-// cif::file refFile(example.string());
-// cif::mm::structure reference(refFile);
+ cif::file refFile(example.string());
+ cif::mm::structure reference(refFile);
-// auto &atomsRea = rea.atoms();
-// auto &refAtomsRea = reference.get_residue("B").atoms();
+ auto &atomsRea = rea.atoms();
+ auto &refAtomsRea = reference.get_residue("B").atoms();
-// BOOST_ASSERT(atomsRea.size() == refAtomsRea.size());
+ BOOST_ASSERT(atomsRea.size() == refAtomsRea.size());
-// for (size_t i = 0; i < atomsRea.size(); ++i)
-// {
-// auto a1 = atomsRea.at(i);
-// auto a2 = refAtomsRea.at(i);
+ for (size_t i = 0; i < atomsRea.size(); ++i)
+ {
+ auto a1 = atomsRea.at(i);
+ auto a2 = refAtomsRea.at(i);
-// std::cout << a1 << ": " << a1.get_location() << " => " << a2.get_location() << " distance: " << distance(a1, a2) << std::endl;
-// }
+ std::cout << a1 << ": " << a1.get_location() << " => " << a2.get_location() << " distance: " << distance(a1, a2) << std::endl;
+ }
-// std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
+ std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
-// pdb_redo::MapMaker<float> mm;
-// float samplingRate = 0.75;
-// mm.loadMTZ(gTestDir / ".." / "examples" / "1cbs_map.mtz", samplingRate);
+ pdb_redo::MapMaker<float> mm;
+ float samplingRate = 0.75;
+ mm.loadMTZ(gTestDir / ".." / "examples" / "1cbs_map.mtz", samplingRate);
-// pdb_redo::BondMap bonds(structure);
+ cif::crystal crystal(structure.get_datablock());
+ auto minimizer = pdb_redo::Minimizer::create(crystal, structure, atoms, mm.fb());
-// auto minimizer = pdb_redo::Minimizer::create(structure, atoms, bonds, mm.fb());
+ minimizer->printStats();
-// minimizer->printStats();
+ auto score = minimizer->refine(true);
-// auto score = minimizer->refine(true);
+ std::cout << "minimizer score: " << score << std::endl;
-// std::cout << "minimizer score: " << score << std::endl;
+ minimizer->printStats();
-// minimizer->printStats();
+ std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
-// std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
+ double d_sum = 0;
-// double d_sum = 0;
+ for (size_t i = 0; i < atomsRea.size(); ++i)
+ {
+ auto a1 = atomsRea.at(i);
+ auto a2 = refAtomsRea.at(i);
-// for (size_t i = 0; i < atomsRea.size(); ++i)
-// {
-// auto a1 = atomsRea.at(i);
-// auto a2 = refAtomsRea.at(i);
+ auto d = distance(a1, a2);
+ d_sum += d * d;
-// auto d = distance(a1, a2);
-// d_sum += d * d;
+ std::cout << a1 << ": " << a1.get_location() << " => " << a2.get_location() << " distance: " << d << std::endl;
+ }
-// std::cout << a1 << ": " << a1.get_location() << " => " << a2.get_location() << " distance: " << d << std::endl;
-// }
+ file.save("/tmp/rsr-test-3.cif");
-// auto rmsd = std::sqrt(d_sum / atomsRea.size());
-// std::cout << "RMSd: " << rmsd << std::endl;
-// BOOST_CHECK(rmsd < 0.2);
+ auto rmsd = std::sqrt(d_sum / atomsRea.size());
+ std::cout << "RMSd: " << rmsd << std::endl;
+ BOOST_CHECK(rmsd < 0.35);
-// std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
-// }
+ std::cout << std::string(cif::get_terminal_width(), '-') << std::endl;
+}