Compare commits

...

15 Commits

Author SHA1 Message Date
Maarten L. Hekkelman
be0f885fa5 Merge branch 'trunk' into dict-for-data 2025-03-10 14:00:19 +01:00
Maarten L. Hekkelman
0459d344e9 Fixes in error reporting 2025-02-17 12:32:14 +01:00
Maarten L. Hekkelman
71e525cd76 Refactored dictionary loading 2025-02-17 09:40:36 +01:00
Maarten L. Hekkelman
1480706d8b change for mingw 2025-02-05 16:05:08 +01:00
Maarten L. Hekkelman
96655b6d80 revert 2025-01-29 17:12:59 +01:00
Maarten L. Hekkelman
eed2aa0d0d better way to include eigen3 2025-01-29 17:01:44 +01:00
Maarten L. Hekkelman
de0c078a23 Update changelog 2025-01-29 16:08:55 +01:00
Maarten L. Hekkelman
321e995a54 Add some comments 2025-01-29 16:07:03 +01:00
Maarten L. Hekkelman
da9f1f81d7 Fix eigen3 problems on github? 2025-01-29 15:57:16 +01:00
Maarten L. Hekkelman
c6d4477a24 Using eigen quaternions 2025-01-29 15:37:57 +01:00
Maarten L. Hekkelman
523b073cdc own eigen 2025-01-29 14:25:28 +01:00
Maarten L. Hekkelman
2591bee21b test for github actions, own eigen library 2025-01-29 13:54:20 +01:00
Maarten L. Hekkelman
d881ca00c9 cleanup 2025-01-29 13:54:00 +01:00
Maarten L. Hekkelman
329dbff474 replace deprecated call 2025-01-29 13:42:37 +01:00
Maarten L. Hekkelman
d84a9fe6dc Deal with missing entity.type 2025-01-29 13:41:11 +01:00
18 changed files with 361 additions and 319 deletions

View File

@@ -27,7 +27,7 @@ cmake_minimum_required(VERSION 3.23)
# set the project name
project(
libcifpp
VERSION 7.0.10
VERSION 8.0.0
LANGUAGES CXX)
list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
@@ -41,6 +41,7 @@ include(CheckCXXSourceCompiles)
include(GenerateExportHeader)
include(CTest)
include(FetchContent)
include(ExternalProject)
# FindBoost, take care of it now.
if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.30)
@@ -236,21 +237,18 @@ if(Eigen3_FOUND AND TARGET Eigen3::Eigen)
get_target_property(EIGEN_INCLUDE_DIR Eigen3::Eigen
INTERFACE_INCLUDE_DIRECTORIES)
else()
# Create a private copy of eigen3 and populate it only, no need to build
FetchContent_Declare(
my-eigen3
# Use ExternalProject since FetchContent always tries to install the result...
ExternalProject_Add(my-eigen3
GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
GIT_TAG 3.4.0)
FetchContent_GetProperties(my-eigen3)
if(NOT my-eigen3_POPULATED)
FetchContent_Populate(my-eigen3)
endif()
set(EIGEN_INCLUDE_DIR ${my-eigen3_SOURCE_DIR})
GIT_TAG 3.4.0
INSTALL_COMMAND "")
ExternalProject_Get_Property(my-eigen3 SOURCE_DIR)
set(EIGEN_INCLUDE_DIR ${SOURCE_DIR})
endif()
message(STATUS "Eigen include dir is ${EIGEN_INCLUDE_DIR}")
# Create a revision file, containing the current git version info
include(VersionString)
write_version_header(${CMAKE_CURRENT_SOURCE_DIR}/src/ LIB_NAME "LibCIFPP")
@@ -337,6 +335,10 @@ set(project_headers
add_library(cifpp)
add_library(cifpp::cifpp ALIAS cifpp)
if(TARGET my-eigen3)
add_dependencies(cifpp my-eigen3)
endif()
target_sources(cifpp
PRIVATE ${project_sources}
${CMAKE_CURRENT_SOURCE_DIR}/src/symop_table_data.hpp

View File

@@ -1,4 +1,12 @@
Version 8.0.0
- A dictionary is for a datablock and a file can have
datablocks with differing dictionaries.
Version 7.0.10
- Deal with missing _entity.type in reconstructing mmCIF files
- Replace code creating quaternions from rotation matrices
that might sometimes give incorrect results. Or at least,
the test code failed on this particular kind of code. Sometimes.
- Fix reconstruction to build pdbx_nonpoly_scheme
Version 7.0.9

View File

@@ -98,6 +98,19 @@ class datablock : public std::list<category>
m_name = name;
}
/**
* @brief Attempt to load the dictionary specified in audit_conform category
*
*/
void load_dictionary();
/**
* @brief Load the dictionary named @a dict_name
*
* @param dict_name
*/
void load_dictionary(std::string_view dict_name);
/**
* @brief Set the validator object to @a v
*

View File

@@ -118,19 +118,6 @@ class file : public std::list<datablock>
/** @endcond */
/**
* @brief Set the validator object to @a v
*/
void set_validator(const validator *v);
/**
* @brief Get the validator object
*/
const validator *get_validator() const
{
return m_validator;
}
/**
* @brief Validate the content and return true if everything was valid.
*
@@ -165,32 +152,6 @@ class file : public std::list<datablock>
*/
bool validate_links() const;
/**
* @brief Attempt to load a dictionary (validator) based on
* the contents of the *audit_conform* category, if available.
*/
void load_dictionary();
/**
* @brief Attempt to load the named dictionary @a name and
* create a validator based on it.
*
* The @a name can be the name of a single file, or even the
* stem of that filename. So, e.g. mmcif_pdbx is valid.
*
* Since libcifpp can use extensions to validators, you
* can add them to the name. So if you would like to add
* the dssp extensions you would have to write:
*
* @code{cpp}
* file.load_dictionary("mmcif_pdbx;dssp-extension");
* @endcode
*
* @param name The name of the dictionary to load
*/
void load_dictionary(std::string_view name);
/**
* @brief Return true if a datablock with the name @a name is part of this file
*/
@@ -243,6 +204,18 @@ class file : public std::list<datablock>
/** Load the data from @a is */
void load(std::istream &is);
/** Load the data from the file specified by @a p using validator @a v */
void load(const std::filesystem::path &p, const validator &v);
/** Load the data from @a is using validator @a v */
void load(std::istream &is, const validator &v);
/** Load the data from the file specified by @a p using a validator constructed from dictionary @a dict */
void load(const std::filesystem::path &p, std::string_view dict);
/** Load the data from @a is using a validator constructed from dictionary @a dict */
void load(std::istream &is, std::string_view dict);
/** Save the data to the file specified by @a p */
void save(const std::filesystem::path &p) const;
@@ -257,9 +230,6 @@ class file : public std::list<datablock>
f.save(os);
return os;
}
private:
const validator *m_validator = nullptr;
};
} // namespace cif

View File

@@ -39,6 +39,8 @@
namespace cif
{
class validator;
// --------------------------------------------------------------------
/** Exception that is thrown when the mmCIF file contains a parsing error */
@@ -307,6 +309,14 @@ class sac_parser
class parser : public sac_parser
{
public:
/// \brief constructor, generates data into @a file from @a is using validator @a v
parser(std::istream &is, file &file, const validator *v)
: sac_parser(is)
, m_file(file)
, m_validator(v)
{
}
/// \brief constructor, generates data into @a file from @a is
parser(std::istream &is, file &file)
: sac_parser(is)
@@ -327,6 +337,7 @@ class parser : public sac_parser
file &m_file;
datablock *m_datablock = nullptr;
category *m_category = nullptr;
const validator *m_validator = nullptr;
row_handle m_row;
/** @endcond */

View File

@@ -718,7 +718,7 @@ bool category::is_valid() const
auto iv = m_cat_validator->get_validator_for_item(col.m_name);
if (iv == nullptr)
{
m_validator->report_error(validation_error::unknown_item, col.m_name, m_name, false);
m_validator->report_error(validation_error::unknown_item, m_name, col.m_name, false);
result = false;
}

View File

@@ -38,6 +38,35 @@ datablock::datablock(const datablock &db)
cat.update_links(*this);
}
void datablock::load_dictionary()
{
if (auto *audit_conform = get("audit_conform"); audit_conform and not audit_conform->empty())
{
std::string name = audit_conform->front().get<std::string>("dict_name");
if (name == "mmcif_pdbx_v50")
name = "mmcif_pdbx.dic"; // we had a bug here in libcifpp...
if (not name.empty())
{
try
{
load_dictionary(name);
}
catch (const std::exception &ex)
{
if (VERBOSE)
std::cerr << "Failed to load dictionary " << std::quoted(name) << ": " << ex.what() << '\n';
}
}
}
}
void datablock::load_dictionary(std::string_view name)
{
set_validator(&validator_factory::instance()[name]);
}
void datablock::set_validator(const validator *v)
{
m_validator = v;
@@ -62,7 +91,7 @@ const validator *datablock::get_validator() const
bool datablock::is_valid() const
{
if (m_validator == nullptr)
throw std::runtime_error("Validator not specified");
throw std::runtime_error("Validator not specified for datablock data_" + name());
bool result = true;
for (auto &cat : *this)
@@ -74,7 +103,7 @@ bool datablock::is_valid() const
bool datablock::is_valid()
{
if (m_validator == nullptr)
throw std::runtime_error("Validator not specified");
throw std::runtime_error("Validator not specified for datablock data_" + name());
bool result = true;
for (auto &cat : *this)

View File

@@ -30,40 +30,8 @@
namespace cif
{
// --------------------------------------------------------------------
// TODO: This is wrong. A validator should be assigned to datablocks,
// not to a file. Since audit_conform is a category specifying the
// content of a datablock. Not the entire file.
void file::set_validator(const validator *v)
{
m_validator = v;
for (bool first = true; auto &db : *this)
{
try
{
db.set_validator(v);
}
catch (const std::exception &e)
{
if (first)
throw;
// Accept failure on secondary datablocks
// now that many mmCIF files have invalid
// restraint data concatenated.
std::cerr << e.what() << '\n';
}
first = false;
}
}
bool file::is_valid() const
{
if (m_validator == nullptr)
std::runtime_error("No validator loaded explicitly, cannot continue");
bool result = true;
for (auto &d : *this)
result = d.is_valid() and result;
@@ -76,14 +44,6 @@ bool file::is_valid() const
bool file::is_valid()
{
if (m_validator == nullptr)
{
if (VERBOSE > 0)
std::cerr << "No dictionary loaded explicitly, loading default\n";
load_dictionary();
}
bool result = not empty();
for (auto &d : *this)
@@ -97,9 +57,6 @@ bool file::is_valid()
bool file::validate_links() const
{
if (m_validator == nullptr)
std::runtime_error("No validator loaded explicitly, cannot continue");
bool result = true;
for (auto &db : *this)
@@ -108,41 +65,41 @@ bool file::validate_links() const
return result;
}
void file::load_dictionary()
{
if (not empty())
{
auto *audit_conform = front().get("audit_conform");
if (audit_conform and not audit_conform->empty())
{
std::string name = audit_conform->front().get<std::string>("dict_name");
// void file::load_dictionary()
// {
// if (not empty())
// {
// auto *audit_conform = front().get("audit_conform");
// if (audit_conform and not audit_conform->empty())
// {
// std::string name = audit_conform->front().get<std::string>("dict_name");
if (name == "mmcif_pdbx_v50")
name = "mmcif_pdbx.dic"; // we had a bug here in libcifpp...
// if (name == "mmcif_pdbx_v50")
// name = "mmcif_pdbx.dic"; // we had a bug here in libcifpp...
if (not name.empty())
{
try
{
load_dictionary(name);
}
catch (const std::exception &ex)
{
if (VERBOSE)
std::cerr << "Failed to load dictionary " << std::quoted(name) << ": " << ex.what() << '\n';
}
}
}
}
// if (not name.empty())
// {
// try
// {
// load_dictionary(name);
// }
// catch (const std::exception &ex)
// {
// if (VERBOSE)
// std::cerr << "Failed to load dictionary " << std::quoted(name) << ": " << ex.what() << '\n';
// }
// }
// }
// }
// if (not m_validator)
// load_dictionary("mmcif_pdbx.dic"); // TODO: maybe incorrect? Perhaps improve?
}
// // if (not m_validator)
// // load_dictionary("mmcif_pdbx.dic"); // TODO: maybe incorrect? Perhaps improve?
// }
void file::load_dictionary(std::string_view name)
{
set_validator(&validator_factory::instance()[name]);
}
// void file::load_dictionary(std::string_view name)
// {
// set_validator(&validator_factory::instance()[name]);
// }
bool file::contains(std::string_view name) const
{
@@ -187,10 +144,7 @@ std::tuple<file::iterator, bool> file::emplace(std::string_view name)
}
if (is_new)
{
i = insert(end(), { name });
i->set_validator(m_validator);
}
assert(i != end());
return std::make_tuple(i, is_new);
@@ -212,18 +166,44 @@ void file::load(const std::filesystem::path &p)
}
}
void file::load(std::istream &is)
void file::load(const std::filesystem::path &p, std::string_view dict)
{
auto saved = m_validator;
set_validator(nullptr);
load(p, validator_factory::instance().operator[](dict));
}
void file::load(std::istream &is, std::string_view dict)
{
load(is, validator_factory::instance().operator[](dict));
}
void file::load(const std::filesystem::path &p, const validator &v)
{
gzio::ifstream in(p);
if (not in.is_open())
throw std::runtime_error("Could not open file '" + p.string() + '\'');
try
{
load(in, v);
}
catch (const std::exception &)
{
throw_with_nested(std::runtime_error("Error reading file '" + p.string() + '\''));
}
}
void file::load(std::istream &is, const validator &v)
{
parser p(is, *this);
p.parse_file();
for (auto &db : *this)
db.set_validator(&v);
}
if (saved != nullptr)
set_validator(saved);
else
load_dictionary();
void file::load(std::istream &is)
{
parser p(is, *this);
p.parse_file();
}
void file::save(const std::filesystem::path &p) const

View File

@@ -837,6 +837,9 @@ void parser::produce_datablock(std::string_view name)
const auto &[iter, ignore] = m_file.emplace(name);
m_datablock = &(*iter);
if (m_validator)
m_datablock->set_validator(m_validator);
}
void parser::produce_category(std::string_view name)

View File

@@ -5909,7 +5909,8 @@ void PDBFileParser::Parse(std::istream &is, cif::file &result)
{
try
{
mDatablock.set_validator(result.get_validator());
if (mDatablock.get_validator() == nullptr)
mDatablock.load_dictionary();
PreParseInput(is);
@@ -6373,10 +6374,11 @@ void read_pdb_file(std::istream &pdbFile, cif::file &cifFile)
{
PDBFileParser p;
cifFile.load_dictionary("mmcif_pdbx.dic");
p.Parse(pdbFile, cifFile);
if (not cifFile.empty() and cifFile.front().get_validator() == nullptr)
cifFile.front().load_dictionary("mmcif_pdbx.dic");
if (not cifFile.is_valid() and cif::VERBOSE >= 0)
std::cerr << "Resulting mmCIF file is not valid!\n";
}
@@ -6421,8 +6423,8 @@ file read(std::istream &is)
}
// Must be a PDB like file, right?
if (result.get_validator() == nullptr)
result.load_dictionary("mmcif_pdbx.dic");
if (not result.empty() and result.front().get_validator() == nullptr)
result.front().load_dictionary("mmcif_pdbx.dic");
return result;
}

View File

@@ -100,10 +100,35 @@ void checkEntities(datablock &db)
for (auto entity : db["entity"].find("formula_weight"_key == null or "formula_weight"_key == 0))
{
const auto &[entity_id, type] = entity.get<std::string, std::string>("id", "type");
auto &&[entity_id, type] = entity.get<std::string, std::string>("id", "type");
float formula_weight = 0;
if (type.empty()) // yes, that happens
{
const auto comp_id = db["atom_site"].find_first<std::string>("label_entity_id"_key == entity_id, "label_comp_id");
auto compound = cf.create(comp_id);
if (compound != nullptr)
{
if (compound->is_base() or compound->is_peptide())
type = "polymer";
else if (compound->is_water())
type = "water";
else
{
if (db["pdbx_entity_branch_link"].contains("entity_id"_key == entity_id))
type = "branched";
else
type = "non-polymer";
}
}
if (type.empty())
throw std::runtime_error("Entity without type and cannot determine what it should be");
entity["type"] = type;
}
if (type == "polymer")
{
int n = 0;
@@ -1012,6 +1037,10 @@ void comparePolySeqSchemes(datablock &db)
auto &ndb_poly_seq_scheme = db["ndb_poly_seq_scheme"];
auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
// Don't bother if ndb_poly_seq_scheme is empty
if (ndb_poly_seq_scheme.empty())
return;
// Since often ndb_poly_seq_scheme only contains an id and mon_id item
// we assume that it should match the accompanying pdbx_poly_seq
@@ -1427,7 +1456,7 @@ bool reconstruct_pdbx(file &file, std::string_view dictionary)
db["chem_comp"].reorder_by_index();
file.load_dictionary(dictionary);
db.load_dictionary(dictionary);
if (db.get("atom_site_anisotrop"))
checkAtomAnisotropRecords(db);

View File

@@ -306,8 +306,8 @@ bool is_valid_pdbx_file(const file &file, std::string_view dictionary, std::erro
if (not seq_can.has_value())
{
if (cif::VERBOSE > 0)
std::clog << "Warning: entity_poly has no sequence for entity_id " << entity_id << '\n';
if (cif::VERBOSE > 1)
std::clog << "Warning: entity_poly has no canonical sequence for entity_id " << entity_id << '\n';
}
else
{

View File

@@ -32,7 +32,7 @@
#include "symop_table_data.hpp"
#include <Eigen/Eigenvalues>
#include <Eigen/Eigen>
namespace cif
{
@@ -103,9 +103,9 @@ sym_op::sym_op(std::string_view s)
auto b = s.data();
auto e = b + s.length();
int rnri = 256; // default to unexisting number
int rnri = 256; // default to unexisting number
auto r = std::from_chars(b, e, rnri);
m_nr = static_cast<uint8_t>(rnri);
m_ta = r.ptr[1] - '0';
m_tb = r.ptr[2] - '0';
@@ -121,7 +121,7 @@ std::string sym_op::string() const
auto r = std::to_chars(b, b + sizeof(b), m_nr);
if ((bool)r.ec or r.ptr > b + 4)
throw std::runtime_error("Could not write out symmetry operation to string");
*r.ptr++ = '_';
*r.ptr++ = '0' + m_ta;
*r.ptr++ = '0' + m_tb;
@@ -163,41 +163,16 @@ transformation::transformation(const matrix3x3<float> &r, const cif::point &t)
void transformation::try_create_quaternion()
{
float Qxx = m_rotation(0, 0);
float Qxy = m_rotation(0, 1);
float Qxz = m_rotation(0, 2);
float Qyx = m_rotation(1, 0);
float Qyy = m_rotation(1, 1);
float Qyz = m_rotation(1, 2);
float Qzx = m_rotation(2, 0);
float Qzy = m_rotation(2, 1);
float Qzz = m_rotation(2, 2);
Eigen::Matrix3f rot;
Eigen::Matrix4f em;
rot << m_rotation(0, 0), m_rotation(0, 1), m_rotation(0, 2),
m_rotation(1, 0), m_rotation(1, 1), m_rotation(1, 2),
m_rotation(2, 0), m_rotation(2, 1), m_rotation(2, 2);
em << Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz;
Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
auto ev = es.eigenvalues();
for (std::size_t j = 0; j < 4; ++j)
if (rot * rot.transpose() == Eigen::Matrix3f::Identity() and rot.determinant() == 1)
{
if (std::abs(ev[j].real() - 1) > 0.01)
continue;
auto col = es.eigenvectors().col(j);
m_q = normalize(cif::quaternion{
static_cast<float>(col(3).real()),
static_cast<float>(col(0).real()),
static_cast<float>(col(1).real()),
static_cast<float>(col(2).real()) });
break;
Eigen::Quaternionf qe(rot);
m_q = normalize(cif::quaternion{ qe.w(), qe.x(), qe.y(), qe.z() });
}
}
@@ -297,7 +272,7 @@ point spacegroup::operator()(const point &pt, const cell &c, sym_op symop) const
{
if (symop.m_nr < 1 or symop.m_nr > size())
throw std::out_of_range("symmetry operator number out of range");
transformation t = at(symop.m_nr - 1);
t.m_translation.m_x += symop.m_ta - 5;
@@ -316,7 +291,7 @@ point spacegroup::inverse(const point &pt, const cell &c, sym_op symop) const
{
if (symop.m_nr < 1 or symop.m_nr > size())
throw std::out_of_range("symmetry operator number out of range");
transformation t = at(symop.m_nr - 1);
t.m_translation.m_x += symop.m_ta - 5;
@@ -450,13 +425,13 @@ int get_space_group_number(const datablock &db)
if (_symmetry.size() != 1)
throw std::runtime_error("Could not find a unique symmetry in this mmCIF file");
return _symmetry.front().get<int>("Int_Tables_number");
}
// --------------------------------------------------------------------
std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b) const
std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b) const
{
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");
@@ -491,7 +466,7 @@ std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b)
while (fsb.m_x + 0.5f < fa.m_x)
{
fsb.m_x += 1;
s.m_ta += 1;
s.m_ta += 1;
}
while (fsb.m_y - 0.5f > fa.m_y)
@@ -503,7 +478,7 @@ std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b)
while (fsb.m_y + 0.5f < fa.m_y)
{
fsb.m_y += 1;
s.m_tb += 1;
s.m_tb += 1;
}
while (fsb.m_z - 0.5f > fa.m_z)
@@ -515,7 +490,7 @@ std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b)
while (fsb.m_z + 0.5f < fa.m_z)
{
fsb.m_z += 1;
s.m_tc += 1;
s.m_tc += 1;
}
auto p = orthogonal(fsb, m_cell);

View File

@@ -63,9 +63,9 @@ std::string get_version_nr()
// --------------------------------------------------------------------
#ifdef _WIN32
#if defined(_WIN32) or defined(__MINGW32__)
}
#include <Windows.h>
#include <windows.h>
#include <libloaderapi.h>
#include <wincon.h>

View File

@@ -401,7 +401,7 @@ void validator::report_error(std::error_code ec, bool fatal) const
{
if (m_strict or fatal)
throw validation_exception(ec);
else
else if (VERBOSE > 0)
std::cerr << ec.message() << '\n';
}
@@ -412,7 +412,7 @@ void validator::report_error(std::error_code ec, std::string_view category,
if (m_strict or fatal)
throw ex;
else
else if (VERBOSE > 0)
std::cerr << ex.what() << '\n';
}

View File

@@ -53,8 +53,8 @@ TEST_CASE("create_nonpoly_1")
cif::VERBOSE = 1;
cif::file file;
file.load_dictionary("mmcif_pdbx.dic");
file.emplace("TEST"); // create a datablock
auto &&[dbi, ignore] = file.emplace("TEST"); // create a datablock
dbi->load_dictionary("mmcif_pdbx.dic");
cif::mm::structure structure(file);
@@ -82,7 +82,7 @@ _atom_site.pdbx_formal_charge
# that's enough to test with
)"_cf;
atoms.load_dictionary("mmcif_pdbx.dic");
atoms.front().load_dictionary("mmcif_pdbx.dic");
auto &hem_data = atoms["HEM"];
auto &atom_site = hem_data["atom_site"];
@@ -159,7 +159,7 @@ _struct_asym.details ?
_atom_type.symbol C
)"_cf;
expected.load_dictionary("mmcif_pdbx.dic");
expected.front().load_dictionary("mmcif_pdbx.dic");
if (not(expected.front() == structure.get_datablock()))
{
@@ -177,8 +177,8 @@ TEST_CASE("create_nonpoly_2")
cif::VERBOSE = 1;
cif::file file;
file.load_dictionary("mmcif_pdbx.dic");
file.emplace("TEST"); // create a datablock
auto &&[dbi, ignore] = file.emplace("TEST"); // create a datablock
dbi->load_dictionary("mmcif_pdbx.dic");
cif::mm::structure structure(file);
@@ -270,7 +270,7 @@ _struct_asym.details ?
_atom_type.symbol C
)"_cf;
expected.load_dictionary("mmcif_pdbx.dic");
expected.front().load_dictionary("mmcif_pdbx.dic");
REQUIRE(expected.front() == structure.get_datablock());
@@ -354,7 +354,7 @@ _struct_asym.details ?
#
)"_cf;
data.load_dictionary("mmcif_pdbx.dic");
data.front().load_dictionary("mmcif_pdbx.dic");
cif::mm::structure s(data);

View File

@@ -148,9 +148,8 @@ TEST_CASE("dh_q_0")
cif::point axis(1, 0, 0);
cif::point p(1, 1, 0);
cif::point t[3] =
{
cif::point t[3] = {
{ 0, 1, 0 },
{ 0, 0, 0 },
{ 1, 0, 0 }
@@ -180,7 +179,6 @@ TEST_CASE("dh_q_0")
a = cif::dihedral_angle(t[0], t[1], t[2], p);
REQUIRE(std::abs(a - 0.f) < 0.01f);
}
TEST_CASE("dh_q_1")
@@ -228,62 +226,103 @@ TEST_CASE("dh_q_1")
// --------------------------------------------------------------------
TEST_CASE("m2q_0, *utf::tolerance(0.001f)")
// TEST_CASE("m2q_0")
// {
// for (std::size_t i = 0; i < cif::kSymopNrTableSize; ++i)
// {
// auto d = cif::kSymopNrTable[i].symop().data();
// cif::matrix3x3<float> rot;
// float Qxx = rot(0, 0) = d[0];
// float Qxy = rot(0, 1) = d[1];
// float Qxz = rot(0, 2) = d[2];
// float Qyx = rot(1, 0) = d[3];
// float Qyy = rot(1, 1) = d[4];
// float Qyz = rot(1, 2) = d[5];
// float Qzx = rot(2, 0) = d[6];
// float Qzy = rot(2, 1) = d[7];
// float Qzz = rot(2, 2) = d[8];
// Eigen::Matrix4f em;
// em << Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
// Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
// Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
// Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz;
// Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
// auto ev = es.eigenvalues();
// std::size_t bestJ = 0;
// float bestEV = -1;
// for (std::size_t j = 0; j < 4; ++j)
// {
// if (bestEV < ev[j].real())
// {
// bestEV = ev[j].real();
// bestJ = j;
// }
// }
// if (std::abs(bestEV - 1) > 0.01)
// continue; // not a rotation matrix
// auto col = es.eigenvectors().col(bestJ);
// auto q = normalize(cif::quaternion{
// static_cast<float>(col(3).real()),
// static_cast<float>(col(0).real()),
// static_cast<float>(col(1).real()),
// static_cast<float>(col(2).real()) });
// cif::point p1{ 1, 1, 1 };
// cif::point p2 = p1;
// p2.rotate(q);
// cif::point p3 = rot * p1;
// REQUIRE_THAT(p2.m_x, Catch::Matchers::WithinRel(p3.m_x, 0.01f));
// REQUIRE_THAT(p2.m_y, Catch::Matchers::WithinRel(p3.m_y, 0.01f));
// REQUIRE_THAT(p2.m_z, Catch::Matchers::WithinRel(p3.m_z, 0.01f));
// }
// }
TEST_CASE("m2q_0a")
{
for (std::size_t i = 0; i < cif::kSymopNrTableSize; ++i)
{
auto d = cif::kSymopNrTable[i].symop().data();
cif::matrix3x3<float> rot;
float Qxx = rot(0, 0) = d[0];
float Qxy = rot(0, 1) = d[1];
float Qxz = rot(0, 2) = d[2];
float Qyx = rot(1, 0) = d[3];
float Qyy = rot(1, 1) = d[4];
float Qyz = rot(1, 2) = d[5];
float Qzx = rot(2, 0) = d[6];
float Qzy = rot(2, 1) = d[7];
float Qzz = rot(2, 2) = d[8];
Eigen::Matrix3f rot;
rot << d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8];
Eigen::Matrix4f em;
// check to see if this matrix contains a true rotation
if (rot * rot.transpose() != Eigen::Matrix3f::Identity() or rot.determinant() != 1)
continue;
em << Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz;
Eigen::Quaternionf qe(rot);
Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
auto q = normalize(cif::quaternion{ qe.w(), qe.x(), qe.y(), qe.z() });
auto ev = es.eigenvalues();
std::size_t bestJ = 0;
float bestEV = -1;
for (std::size_t j = 0; j < 4; ++j)
{
if (bestEV < ev[j].real())
{
bestEV = ev[j].real();
bestJ = j;
}
}
if (std::abs(bestEV - 1) > 0.01)
continue; // not a rotation matrix
auto col = es.eigenvectors().col(bestJ);
auto q = normalize(cif::quaternion{
static_cast<float>(col(3).real()),
static_cast<float>(col(0).real()),
static_cast<float>(col(1).real()),
static_cast<float>(col(2).real()) });
cif::point p1{ 1, 1, 1 };
cif::point p2 = p1;
p2.rotate(q);
cif::point p3 = rot * p1;
cif::matrix3x3<float> rot_c({
rot_c(0, 0) = d[0],
rot_c(0, 1) = d[1],
rot_c(0, 2) = d[2],
rot_c(1, 0) = d[3],
rot_c(1, 1) = d[4],
rot_c(1, 2) = d[5],
rot_c(2, 0) = d[6],
rot_c(2, 1) = d[7],
rot_c(2, 2) = d[8]
});
cif::point p3 = rot_c * p1;
REQUIRE_THAT(p2.m_x, Catch::Matchers::WithinRel(p3.m_x, 0.01f));
REQUIRE_THAT(p2.m_y, Catch::Matchers::WithinRel(p3.m_y, 0.01f));
@@ -291,7 +330,7 @@ TEST_CASE("m2q_0, *utf::tolerance(0.001f)")
}
}
// "TEST_CASE(m2q_1, *utf::tolerance(0.001f)")
// "TEST_CASE(m2q_1")
// {
// for (std::size_t i = 0; i < cif::kSymopNrTableSize; ++i)
// {
@@ -337,7 +376,7 @@ TEST_CASE("m2q_0, *utf::tolerance(0.001f)")
// static_cast<float>(em(bestJ, 0)),
// static_cast<float>(em(bestJ, 1)),
// static_cast<float>(em(bestJ, 2)) });
// cif::point p1{ 1, 1, 1 };
// cif::point p2 = p1;
// p2.rotate(q);
@@ -390,17 +429,17 @@ TEST_CASE("symm_3")
REQUIRE(sg.get_name() == "P 21 21 2");
}
TEST_CASE("symm_4, *utf::tolerance(0.1f)")
TEST_CASE("symm_4")
{
using namespace cif::literals;
// based on 2b8h
auto sg = cif::spacegroup(154); // p 32 2 1
auto c = cif::cell(107.516, 107.516, 338.487, 90.00, 90.00, 120.00);
cif::point a{ -8.688, 79.351, 10.439 }; // O6 NAG A 500
cif::point b{ -35.356, 33.693, -3.236 }; // CG2 THR D 400
cif::point sb( -6.916, 79.34, 3.236); // 4_565 copy of b
cif::point a{ -8.688, 79.351, 10.439 }; // O6 NAG A 500
cif::point b{ -35.356, 33.693, -3.236 }; // CG2 THR D 400
cif::point sb(-6.916, 79.34, 3.236); // 4_565 copy of b
REQUIRE_THAT(distance(a, sg(a, c, "1_455"_symop)), Catch::Matchers::WithinRel(static_cast<float>(c.get_a()), 0.01f));
REQUIRE_THAT(distance(a, sg(a, c, "1_545"_symop)), Catch::Matchers::WithinRel(static_cast<float>(c.get_b()), 0.01f));
@@ -411,12 +450,12 @@ TEST_CASE("symm_4, *utf::tolerance(0.1f)")
REQUIRE_THAT(sb.m_y, Catch::Matchers::WithinRel(sb2.m_y, 0.01f));
REQUIRE_THAT(sb.m_z, Catch::Matchers::WithinRel(sb2.m_z, 0.01f));
REQUIRE_THAT(distance(a, sb2), Catch::Matchers::WithinRel(7.42f, 0.01f));
REQUIRE_THAT(distance(a, sb2), Catch::Matchers::WithinRel(7.42f, 0.01f));
}
// --------------------------------------------------------------------
TEST_CASE("symm_4wvp_1, *utf::tolerance(0.1f)")
TEST_CASE("symm_4wvp_1")
{
using namespace cif::literals;
@@ -427,7 +466,7 @@ TEST_CASE("symm_4wvp_1, *utf::tolerance(0.1f)")
cif::crystal c(db);
cif::point p{ -78.722, 98.528, 11.994 };
cif::point p{ -78.722, 98.528, 11.994 };
auto a = s.get_residue("A", 10, "").get_atom_by_atom_id("O");
auto sp1 = c.symmetry_copy(a.get_location(), "2_565"_symop);
@@ -442,10 +481,9 @@ TEST_CASE("symm_4wvp_1, *utf::tolerance(0.1f)")
REQUIRE_THAT(sp2.m_x, Catch::Matchers::WithinAbs(p.m_x, 0.5f));
REQUIRE_THAT(sp2.m_y, Catch::Matchers::WithinAbs(p.m_y, 0.5f));
REQUIRE_THAT(sp2.m_z, Catch::Matchers::WithinAbs(p.m_z, 0.5f));
}
TEST_CASE("symm_2bi3_1, *utf::tolerance(0.1f)")
TEST_CASE("symm_2bi3_1")
{
cif::file f(gTestDir / "2bi3.cif.gz");
@@ -455,18 +493,15 @@ TEST_CASE("symm_2bi3_1, *utf::tolerance(0.1f)")
cif::crystal c(db);
auto struct_conn = db["struct_conn"];
for (const auto &[
asym1, seqid1, authseqid1, atomid1, symm1,
asym2, seqid2, authseqid2, atomid2, symm2,
dist] : struct_conn.find<
std::string,int,std::string,std::string,std::string,
std::string,int,std::string,std::string,std::string,
float>(
cif::key("ptnr1_symmetry") != "1_555" or cif::key("ptnr2_symmetry") != "1_555",
"ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr1_auth_seq_id", "ptnr1_label_atom_id", "ptnr1_symmetry",
"ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
"pdbx_dist_value"
))
for (const auto &[asym1, seqid1, authseqid1, atomid1, symm1,
asym2, seqid2, authseqid2, atomid2, symm2,
dist] : struct_conn.find<std::string, int, std::string, std::string, std::string,
std::string, int, std::string, std::string, std::string,
float>(
cif::key("ptnr1_symmetry") != "1_555" or cif::key("ptnr2_symmetry") != "1_555",
"ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr1_auth_seq_id", "ptnr1_label_atom_id", "ptnr1_symmetry",
"ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
"pdbx_dist_value"))
{
auto &r1 = s.get_residue(asym1, seqid1, authseqid1);
auto &r2 = s.get_residue(asym2, seqid2, authseqid2);
@@ -492,7 +527,7 @@ TEST_CASE("symm_2bi3_1, *utf::tolerance(0.1f)")
}
}
TEST_CASE("symm_2bi3_1a, *utf::tolerance(0.1f)")
TEST_CASE("symm_2bi3_1a")
{
using namespace cif::literals;
@@ -504,23 +539,20 @@ TEST_CASE("symm_2bi3_1a, *utf::tolerance(0.1f)")
auto struct_conn = db["struct_conn"];
auto atom_site = db["atom_site"];
for (const auto &[
asym1, seqid1, authseqid1, atomid1, symm1,
asym2, seqid2, authseqid2, atomid2, symm2,
dist] : struct_conn.find<
std::string,std::optional<int>,std::string,std::string,std::string,
std::string,std::optional<int>,std::string,std::string,std::string,
float>(
cif::key("ptnr1_symmetry") != "1_555" or cif::key("ptnr2_symmetry") != "1_555",
"ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr1_auth_seq_id", "ptnr1_label_atom_id", "ptnr1_symmetry",
"ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
"pdbx_dist_value"
))
for (const auto &[asym1, seqid1, authseqid1, atomid1, symm1,
asym2, seqid2, authseqid2, atomid2, symm2,
dist] : struct_conn.find<std::string, std::optional<int>, std::string, std::string, std::string,
std::string, std::optional<int>, std::string, std::string, std::string,
float>(
cif::key("ptnr1_symmetry") != "1_555" or cif::key("ptnr2_symmetry") != "1_555",
"ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr1_auth_seq_id", "ptnr1_label_atom_id", "ptnr1_symmetry",
"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>(
cif::point p1 = 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>(
cif::point p2 = 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");
@@ -540,7 +572,7 @@ TEST_CASE("symm_2bi3_1a, *utf::tolerance(0.1f)")
}
}
TEST_CASE("symm_3bwh_1, *utf::tolerance(0.1f)")
TEST_CASE("symm_3bwh_1")
{
cif::file f(gTestDir / "3bwh.cif.gz");
@@ -555,15 +587,15 @@ TEST_CASE("symm_3bwh_1, *utf::tolerance(0.1f)")
{
if (a1 == a2)
continue;
const auto&[ d, p, so ] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
const auto &[d, p, so] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
REQUIRE_THAT(d, Catch::Matchers::WithinAbs(distance(a1.get_location(), p), 0.5f));
}
}
}
TEST_CASE("volume_3bwh_1, *utf::tolerance(0.1f)")
TEST_CASE("volume_3bwh_1")
{
cif::file f(gTestDir / "1juh.cif.gz");
@@ -573,4 +605,3 @@ TEST_CASE("volume_3bwh_1, *utf::tolerance(0.1f)")
REQUIRE_THAT(c.get_cell().get_volume(), Catch::Matchers::WithinRel(741009.625f, 0.01f));
}

View File

@@ -774,7 +774,6 @@ save__cat_2.desc
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -805,7 +804,7 @@ _cat_2.desc
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
SECTION("one")
{
@@ -927,7 +926,6 @@ save__cat_1.c
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -950,7 +948,7 @@ mies Mies
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
@@ -1089,7 +1087,6 @@ save__cat_2.desc
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -1123,7 +1120,7 @@ _cat_2.desc
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -1292,7 +1289,6 @@ save__cat_2.parent_id3
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -1336,7 +1332,7 @@ _cat_2.parent_id3
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -1513,7 +1509,6 @@ cat_2 3 cat_2:cat_1:3
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -1550,7 +1545,7 @@ _cat_2.parent_id3
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -1753,7 +1748,6 @@ cat_2 1 cat_2:cat_1:1
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -1790,7 +1784,7 @@ _cat_2.parent_id_2
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
// auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -2149,7 +2143,6 @@ cat_2 1 '_cat_2.num' '_cat_3.num' cat_3
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -2191,7 +2184,7 @@ _cat_3.num
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -2434,7 +2427,6 @@ cat_2 1 '_cat_2.num' '_cat_3.num' cat_3
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -2476,7 +2468,7 @@ _cat_3.num
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];
auto &cat2 = f.front()["cat_2"];
@@ -3031,7 +3023,6 @@ save__cat_1.name
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -3054,7 +3045,7 @@ _cat_1.name
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
REQUIRE(f.is_valid());
@@ -3226,7 +3217,6 @@ save__cat_1.name
auto &validator = cif::validator_factory::instance().construct_validator("test_dict.dic", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -3253,7 +3243,7 @@ _cat_1.name
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
REQUIRE(f.is_valid());
@@ -3333,7 +3323,6 @@ save__cat_1.id_2
auto validator = cif::parse_dictionary("test", is_dict);
cif::file f;
f.set_validator(&validator);
// --------------------------------------------------------------------
@@ -3358,7 +3347,7 @@ _cat_1.id_2
} data_buffer(const_cast<char *>(data), sizeof(data) - 1);
std::istream is_data(&data_buffer);
f.load(is_data);
f.load(is_data, validator);
auto &cat1 = f.front()["cat_1"];