mirror of
https://github.com/PDB-REDO/libcifpp.git
synced 2026-06-04 13:54:25 +08:00
Compiling, but not working yet
This commit is contained in:
Maarten L. Hekkelman
@@ -8,7 +8,7 @@ Points
|
||||
|
||||
The most basic type in use is :cpp:type:`cif::point`. It can be thought of as a point in space with three coordinates, but it is also often used as a vector in 3d space. To keep the interface simple there's no separate vector type.
|
||||
|
||||
Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`cif::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
|
||||
Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`glm::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
|
||||
|
||||
Quaternions
|
||||
-----------
|
||||
@@ -91,7 +91,7 @@ To give an idea how this works, here's a piece of code copied from one of the un
|
||||
auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
|
||||
|
||||
// The distance between these symmetry atoms should be equal to the distance in the struct_conn record
|
||||
assert(cif::distance(sa1, sa2) == dist);
|
||||
assert(glm::distance(sa1, sa2) == dist);
|
||||
|
||||
// And to show how you can obtain the closest symmetry copy of an atom near another one:
|
||||
// here we request the symmetry copy of p2 that lies closest to p1
|
||||
|
||||
@@ -33,8 +33,11 @@
|
||||
#include <cstdlib>
|
||||
#include <format>
|
||||
#include <functional>
|
||||
#include <glm/ext/quaternion_geometric.hpp>
|
||||
#include <glm/ext/quaternion_trigonometric.hpp>
|
||||
#include <glm/glm.hpp>
|
||||
#include <glm/gtc/quaternion.hpp>
|
||||
#include <glm/trigonometric.hpp>
|
||||
#include <limits>
|
||||
#include <numbers>
|
||||
#include <optional>
|
||||
@@ -205,10 +208,17 @@ std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts);
|
||||
// --------------------------------------------------------------------
|
||||
|
||||
/// \brief Return a quaternion created from angle @a angle and axis @a axis
|
||||
quaternion construct_from_angle_axis(float angle, point axis);
|
||||
constexpr quaternion construct_from_angle_axis(float angle, point axis)
|
||||
{
|
||||
glm::normalize(axis);
|
||||
return glm::angleAxis(glm::radians(angle), axis);
|
||||
}
|
||||
|
||||
/// \brief Return a tuple of an angle and an axis for quaternion @a q
|
||||
std::tuple<float, point> quaternion_to_angle_axis(quaternion q);
|
||||
constexpr std::tuple<float, point> quaternion_to_angle_axis(quaternion q)
|
||||
{
|
||||
return { glm::degrees(glm::angle(q)), glm::axis(q) };
|
||||
}
|
||||
|
||||
/// @brief Given four points and an angle, return the quaternion required to rotate
|
||||
/// point p4 along the p2-p3 axis and around point p3 to obtain the required within
|
||||
@@ -236,3 +246,17 @@ double RMSd(const std::vector<point> &a, const std::vector<point> &b);
|
||||
|
||||
|
||||
} // namespace cif
|
||||
|
||||
template <>
|
||||
struct std::formatter<glm::vec3> : std::formatter<std::string_view> // NOLINT
|
||||
{
|
||||
template <class FmtContext>
|
||||
auto format(glm::vec3 pt, FmtContext &ctx) const
|
||||
{
|
||||
std::string temp;
|
||||
std::format_to(std::back_inserter(temp), "( {}, {}, {} )", pt.x, pt.y, pt.z);
|
||||
|
||||
return std::formatter<std::string_view>::format(temp, ctx);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@@ -2288,17 +2288,15 @@ std::string structure::create_non_poly(const std::string &compound_id, bool skip
|
||||
if (skip_hydrogen and cif::atom_type_traits(a.type_symbol).symbol() == "H")
|
||||
continue;
|
||||
|
||||
auto ax = a.get_location().get_x();
|
||||
auto ay = a.get_location().get_y();
|
||||
auto az = a.get_location().get_z();
|
||||
auto aloc = a.get_location();
|
||||
|
||||
atoms.emplace_back(cif::row_initializer{
|
||||
{ "type_symbol", cif::atom_type_traits(a.type_symbol).symbol() },
|
||||
{ "label_atom_id", a.id },
|
||||
{ "auth_atom_id", a.id },
|
||||
{ "Cartn_x", ax },
|
||||
{ "Cartn_y", ay },
|
||||
{ "Cartn_z", az },
|
||||
{ "Cartn_x", aloc.x },
|
||||
{ "Cartn_y", aloc.y },
|
||||
{ "Cartn_z", aloc.z },
|
||||
{ "B_iso_or_equiv", 30.00 } });
|
||||
}
|
||||
|
||||
|
||||
@@ -257,29 +257,29 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
|
||||
|
||||
// --------------------------------------------------------------------
|
||||
|
||||
std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
|
||||
std::tuple<point, float> smallest_sphere_around_2_points(std::array<point, 2> pts)
|
||||
{
|
||||
return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
|
||||
return { (pts[0] + pts[1]) / 2.0f, distance(pts[0], pts[1]) / 2.0f };
|
||||
}
|
||||
|
||||
std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
|
||||
std::tuple<point, float> smallest_sphere_around_3_points(std::array<point, 3> pts)
|
||||
{
|
||||
// Find two bisectors
|
||||
auto vz = cross(pts[1] - pts[0], pts[2] - pts[0]);
|
||||
|
||||
auto bs1 = cross(vz, pts[1] - pts[0]);
|
||||
bs1.normalize();
|
||||
normalize(bs1);
|
||||
|
||||
auto v1 = (pts[1] - pts[0]);
|
||||
v1.normalize();
|
||||
normalize(v1);
|
||||
|
||||
auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
|
||||
|
||||
auto bs2 = cross(vz, pts[2] - pts[0]);
|
||||
bs2.normalize();
|
||||
normalize(bs2);
|
||||
|
||||
auto v2 = (pts[2] - pts[0]);
|
||||
v2.normalize();
|
||||
normalize(v2);
|
||||
|
||||
auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
|
||||
|
||||
@@ -300,7 +300,7 @@ std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point,
|
||||
return smallest_sphere_around_2_points({ pts[1], pts[2] });
|
||||
}
|
||||
|
||||
std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
|
||||
std::tuple<point, float> smallest_sphere_around_4_points(std::array<point, 4> pts)
|
||||
{
|
||||
auto t0 = -norm_squared(pts[0]);
|
||||
auto t1 = -norm_squared(pts[1]);
|
||||
@@ -427,19 +427,19 @@ bool point_in_circle(point p, std::vector<point> c)
|
||||
case 2:
|
||||
{
|
||||
auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
|
||||
return cif::distance_squared(p, center) <= radius * radius;
|
||||
return distance_squared(p, center) <= radius * radius;
|
||||
}
|
||||
|
||||
case 3:
|
||||
{
|
||||
auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
|
||||
return cif::distance_squared(p, center) <= radius * radius;
|
||||
return distance_squared(p, center) <= radius * radius;
|
||||
}
|
||||
|
||||
case 4:
|
||||
{
|
||||
auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
|
||||
return cif::distance_squared(p, center) <= radius * radius;
|
||||
return distance_squared(p, center) <= radius * radius;
|
||||
}
|
||||
|
||||
default:
|
||||
|
||||
@@ -448,7 +448,7 @@ std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b
|
||||
if (m_cell.get_a() == 0 or m_cell.get_b() == 0 or m_cell.get_c() == 0)
|
||||
throw std::runtime_error("Invalid cell, contains a dimension that is zero");
|
||||
|
||||
point result_fsb;
|
||||
point result_fsb{};
|
||||
float result_d = std::numeric_limits<float>::max();
|
||||
sym_op result_s;
|
||||
|
||||
|
||||
@@ -443,7 +443,7 @@ TEST_CASE("test_alternates_1")
|
||||
const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
|
||||
cif::file file(example.string());
|
||||
|
||||
auto &db = file.front();
|
||||
// auto &db = file.front();
|
||||
|
||||
cif::mm::structure s(file);
|
||||
|
||||
|
||||
@@ -87,7 +87,7 @@ TEST_CASE("t1")
|
||||
cif::center_points(p1);
|
||||
|
||||
for (auto &p : p2)
|
||||
p.rotate(q);
|
||||
p = q * p;
|
||||
|
||||
cif::center_points(p2);
|
||||
|
||||
@@ -98,7 +98,7 @@ TEST_CASE("t1")
|
||||
CHECK_THAT(std::fmod(360 + angle, 360), Catch::Matchers::WithinRel(std::fmod(360 - angle0, 360), 0.01));
|
||||
|
||||
for (auto &p : p1)
|
||||
p.rotate(q2);
|
||||
p = q2 * p;
|
||||
|
||||
auto rmsd = cif::RMSd(p1, p2);
|
||||
|
||||
@@ -115,7 +115,7 @@ TEST_CASE("t2")
|
||||
{ 1, 2, 0 }
|
||||
};
|
||||
|
||||
cif::point xp = cif::cross(p[1] - p[0], p[2] - p[0]);
|
||||
auto xp = glm::cross(p[1] - p[0], p[2] - p[0]);
|
||||
|
||||
auto q = cif::construct_from_angle_axis(45, xp);
|
||||
|
||||
@@ -132,16 +132,16 @@ TEST_CASE("t3")
|
||||
{ 1, 2, 0 }
|
||||
};
|
||||
|
||||
cif::point xp = cif::cross(p[1] - p[0], p[2] - p[0]);
|
||||
cif::point xp = glm::cross(p[1] - p[0], p[2] - p[0]);
|
||||
|
||||
auto q = cif::construct_from_angle_axis(45, xp);
|
||||
|
||||
auto v = p[1];
|
||||
v -= p[0];
|
||||
v.rotate(q);
|
||||
v = q * v;
|
||||
v += p[0];
|
||||
|
||||
std::cout << v << '\n';
|
||||
std::println(std::cout, "{}", v);
|
||||
|
||||
double a = cif::angle(v, p[0], p[1]);
|
||||
|
||||
@@ -165,7 +165,7 @@ TEST_CASE("dh_q_0")
|
||||
|
||||
auto q = cif::construct_from_angle_axis(90, axis);
|
||||
|
||||
p.rotate(q);
|
||||
p = q * p;
|
||||
|
||||
REQUIRE(std::abs(p.x - 1.f) < 0.01f);
|
||||
REQUIRE(std::abs(p.y - 0.f) < 0.01f);
|
||||
@@ -176,7 +176,7 @@ TEST_CASE("dh_q_0")
|
||||
|
||||
q = cif::construct_from_angle_axis(-90, axis);
|
||||
|
||||
p.rotate(q);
|
||||
p = q * p;
|
||||
|
||||
REQUIRE(std::abs(p.x - 1.f) < 0.01f);
|
||||
REQUIRE(std::abs(p.y - 1.f) < 0.01f);
|
||||
@@ -222,7 +222,9 @@ TEST_CASE("dh_q_1")
|
||||
{
|
||||
auto q = cif::construct_for_dihedral_angle(pts[0], pts[1], pts[2], pts[3], angle, 1);
|
||||
|
||||
pts[3].rotate(q, pts[2]);
|
||||
pts[3] -= pts[2];
|
||||
pts[3] = q * pts[3];
|
||||
pts[3] += pts[2];
|
||||
|
||||
auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
|
||||
CHECK_THAT(dh, Catch::Matchers::WithinRel(angle, 0.1f));
|
||||
@@ -248,7 +250,7 @@ TEST_CASE("dh_q_1")
|
||||
|
||||
cif::point p1{ 1, 1, 1 };
|
||||
cif::point p2 = p1;
|
||||
p2.rotate(q);
|
||||
p2 = q * p2;
|
||||
|
||||
cif::matrix3x3<float> rot_c({ static_cast<float>(d[0]),
|
||||
static_cast<float>(d[1]),
|
||||
@@ -318,7 +320,7 @@ TEST_CASE("dh_q_1")
|
||||
|
||||
// cif::point p1{ 1, 1, 1 };
|
||||
// cif::point p2 = p1;
|
||||
// p2.rotate(q);
|
||||
// p2 = q * p2;
|
||||
|
||||
// cif::point p3 = rot * p1;
|
||||
|
||||
@@ -453,7 +455,7 @@ TEST_CASE("symm_2bi3_1")
|
||||
auto sa1 = c.symmetry_copy(a1.get_location(), cif::sym_op(symm1));
|
||||
auto sa2 = c.symmetry_copy(a2.get_location(), cif::sym_op(symm2));
|
||||
|
||||
CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
|
||||
CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
|
||||
|
||||
auto pa1 = a1.get_location();
|
||||
|
||||
@@ -491,17 +493,20 @@ TEST_CASE("symm_2bi3_1a")
|
||||
"ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
|
||||
"pdbx_dist_value"))
|
||||
{
|
||||
cif::point p1 = atom_site.find1<float, float, float>(
|
||||
auto t1 = atom_site.find1<float, float, float>(
|
||||
"label_asym_id"_key == asym1 and "label_seq_id"_key == seqid1 and "auth_seq_id"_key == authseqid1 and "label_atom_id"_key == atomid1,
|
||||
"cartn_x", "cartn_y", "cartn_z");
|
||||
cif::point p2 = atom_site.find1<float, float, float>(
|
||||
auto t2 = atom_site.find1<float, float, float>(
|
||||
"label_asym_id"_key == asym2 and "label_seq_id"_key == seqid2 and "auth_seq_id"_key == authseqid2 and "label_atom_id"_key == atomid2,
|
||||
"cartn_x", "cartn_y", "cartn_z");
|
||||
|
||||
cif::point p1{ std::get<0>(t1), std::get<1>(t1), std::get<2>(t1) };
|
||||
cif::point p2{ std::get<0>(t2), std::get<1>(t2), std::get<2>(t2) };
|
||||
|
||||
auto sa1 = c.symmetry_copy(p1, cif::sym_op(symm1));
|
||||
auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
|
||||
|
||||
CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
|
||||
CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
|
||||
|
||||
const auto &[d, p, so] = c.closest_symmetry_copy(p1, p2);
|
||||
|
||||
@@ -582,7 +587,7 @@ TEST_CASE("smallest_sphere-1")
|
||||
for (int i = 0; i < 1000; ++i)
|
||||
{
|
||||
auto [c, r] = cif::smallest_sphere_around_points(pts);
|
||||
CHECK_THAT(cif::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
|
||||
CHECK_THAT(glm::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
|
||||
CHECK_THAT(r, Catch::Matchers::WithinAbs(7.31248331f, 0.01f));
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user