1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
// This program reads PDB files and reports discrepancies between
// symmetry from CRYST1 and symmetry from REMARK 290.
#include <cstdio>
#include <gemmi/dirwalk.hpp> // for PdbWalk
#include <gemmi/gz.hpp> // for MaybeGzipped
#include <gemmi/pdb.hpp> // for read_pdb
#include <gemmi/remarks.hpp> // for read_remark_290
void check_remark290(const std::string& path) {
using namespace gemmi;
Structure st = read_pdb(MaybeGzipped(path));
std::vector<Op> ops = read_remark_290(st.raw_remarks);
if (ops.empty()) {
if (st.cell.is_crystal())
printf("no remark 290: %s\n", path.c_str());
return;
}
const SpaceGroup* sg = find_spacegroup_by_ops(split_centering_vectors(ops));
if (sg) {
const SpaceGroup* cryst1_sg = st.find_spacegroup();
if (sg != cryst1_sg) {
printf("%s\n", path.c_str());
printf("CRYST1: %s -> %s\n", st.spacegroup_hm.c_str(),
cryst1_sg ? cryst1_sg->xhm().c_str() : "n/a");
printf("REMARK 290: %zu symops -> %s\n", ops.size(), sg->xhm().c_str());
}
} else {
printf("failed remark 290 ops to spacegroup: %s\n", path.c_str());
}
}
int main(int argc, char* argv[]) {
int counter = 0;
try {
for (int pos = 1; pos < argc; ++pos) {
for (const std::string& path : gemmi::PdbWalk(argv[pos])) {
check_remark290(path);
if (++counter % 10000 == 0) {
std::printf("[progress: %d files]\n", counter);
std::fflush(stdout);
}
}
}
} catch (std::runtime_error& err) {
std::fprintf(stderr, "Error: %s\n", err.what());
return 1;
}
return 0;
}
|
