diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 4ba4cee51..764d41642 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.166 +* `rugnux --model` writes the model as it was placed against the data to `_model.cif`, in the cell and space group of the reflection files beside it. * `rugnux --model` places the model against the data as a rigid body before scoring it, keeping the move only if R-free improves, and writes sigma_A-weighted 2mFo-DFc and mFo-DFc maps in place of the unweighted 2Fo-Fc and Fo-Fc. * `rugnux --model` lets the model decide the enantiomorph and the indexing only where its R-work beats that of the same model in random orientations; a model the data reject leaves the reflection files byte for byte what a run with no model writes, and a model that claims neither is not put to that test at all. * The rugnux results report says whether the model was accepted, rejected or never tested and what it was allowed to decide as `MODEL_FIT`, `MODEL_FIT_SIGMA` and `MODEL_DECISIONS_TAKEN`, `MODEL_VALIDATION= PERFORMED` is written on success, and `SPACE_GROUP_ENANTIOMORPH= DETERMINED_FROM_MODEL` becomes `ASSUMED_FROM_MODEL`; `REPORT_VERSION` is 6. diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index a4db67fbc..bceaaa2e0 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -221,4 +221,6 @@ The refinement walks a coarse-to-fine ladder, 6 Å → 4.5 Å → 3.5 Å, each z The refinement sees **only the working reflections**. The step is then committed only if it lowers **R-free**, computed at full resolution on the free set it never saw; otherwise the model is put back exactly where it was read and the maps are the ones it would have given. The whole addition costs about two seconds. Note that the origin is a gauge in some space groups — free in all three directions in $P1$, and along the unique axis in a polar group — so those components of the translation are undetermined; nothing is done about that beyond the Levenberg–Marquardt damping and the R-free gate, which between them make an undetermined direction harmless rather than unstable. -The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_DEG`, `RIGID_BODY_SHIFT_A`), together with the R-free before it (`R_FREE_BEFORE_RIGID_BODY`), so what the placement bought is visible. The **placed coordinates are not written out** — the deliverable is the map — so a model overlaid on these maps is the file as it was supplied, out by the reported shift. +The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_DEG`, `RIGID_BODY_SHIFT_A`), together with the R-free before it (`R_FREE_BEFORE_RIGID_BODY`), so what the placement bought is visible. + +The coordinates **as placed** are written as `_model.cif`, so there is a coordinate file that describes the maps: the input's chains, residues, ligands, waters, B-factors, occupancies and anisotropic $U$s, moved, in the same cell and space group as `.mtz`. Taking the frame from the written reflections rather than from the input model matters — §14.6 may have relabelled them to the model's enantiomorph, which is neither the data's original group nor, necessarily, the model's. It is written **whenever the maps are**, not only where the rigid-body step was committed: the model is re-fractionalized into the data cell and may be relabelled whatever the placement decided, so an unmoved model is still not the input file; and a model the null of §14.5 rejected is scored, placed and mapped like any other, which is precisely the case where the density is worth looking at. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 886140ca6..aa32caee3 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -116,8 +116,18 @@ anisotropic *B* — and reports **R-work / R-free** and the mean 2Fo-Fc density It also writes `_2fofc.ccp4`, `_fofc.ccp4` and `_maps.mtz` next to the merged reflections — and, where the merge kept the Bijvoet split (a rotation merge always does), `_anom.ccp4`, the anomalous difference map whose strongest sites the report names -(`ANOMALOUS_SITE_01`…`10`). Nothing about the model is refined; it is only re-fractionalized into the data -cell, so a deposited model with a slightly different cell still lines up. +(`ANOMALOUS_SITE_01`…`10`). The structure itself is not refined; the model is re-fractionalized into +the data cell and then placed as **one rigid body**, so a deposited model from a crystal that is not +quite isomorphous still sits where the density is. + +Because the model moves, the input file no longer describes these maps, so the model **as placed** is +written as `_model.cif` — the input's chains, residues, ligands, waters, B-factors, +occupancies and anisotropic *U*s, at the coordinates the maps were computed from, in the same unit +cell and space group as `.mtz` beside it. That is the file to open with the maps, and the one +to hand to REFMAC5 or `phenix.refine`; its starting R-free is the `R_FREE=` the report quotes. +(`.cif` is the merged **reflections** — hence the separate name.) It is written whenever the +maps are, including for a model the data rejected: the rejection is a result, and it is exactly the +case where someone wants to look at the model in the density. The model may be **PDB or mmCIF**, gzipped or not, and the format is taken from the file's own content rather than from its name — a model downloaded as `.cif`, `.pdb`, `.ent` or with no useful diff --git a/gemmi_gph/CMakeLists.txt b/gemmi_gph/CMakeLists.txt index 323f3e88f..c8de1d26e 100644 --- a/gemmi_gph/CMakeLists.txt +++ b/gemmi_gph/CMakeLists.txt @@ -6,10 +6,14 @@ # takes a .cif as readily as a .pdb. They bring GEMMI's own copy of the PEGTL parser with them # (gemmi/third_party/tao, MIT, Colin Hirsch and Daniel Frey) - header-only, and included from # gemmi/cif.hpp exactly where upstream GEMMI expects it. +# +# to_mmcif/to_pdb are the writing side, which `rugnux --model` needs to put the model it has placed +# against the data on disk. to_mmcif.cpp calls use_hetatm() out of to_pdb.cpp, so the two come together. ADD_LIBRARY(gemmi STATIC symmetry.cpp gz.cpp mtz.cpp sprintf.cpp xds_ascii.cpp pdb.cpp resinfo.cpp polyheur.cpp calculate.cpp eig3.cpp ccp4.cpp read_cif.cpp mmcif.cpp json.cpp mmread_gz.cpp + to_mmcif.cpp to_pdb.cpp gemmi/cellred.hpp gemmi/symmetry.hpp gemmi/fail.hpp diff --git a/gemmi_gph/to_mmcif.cpp b/gemmi_gph/to_mmcif.cpp new file mode 100644 index 000000000..e97fe3ed7 --- /dev/null +++ b/gemmi_gph/to_mmcif.cpp @@ -0,0 +1,1321 @@ +// Copyright 2017-2023 Global Phasing Ltd. + +#include + +#include +#include // for isnan +#include +#include +#include // std::pair + +#include // no_sign_atoi +#include +#include // for entity_type_to_string, ... +#include // for pdbx_one_letter_code, ... +#include // for use_hetatm + +namespace gemmi { + +namespace { + +inline std::string pdbx_icode(const SeqId& seqid) { + return std::string(1, seqid.has_icode() ? seqid.icode : '?'); +} +inline std::string pdbx_icode(const ResidueId& rid) { + return pdbx_icode(rid.seqid); +} + +inline std::string subchain_or_dot(const Residue& res) { + return res.subchain.empty() ? "." : cif::quote(res.subchain); +} + +inline std::string number_or_dot(double d) { + return std::isnan(d) ? "." : to_str(d); +} +inline std::string number_or_qmark(double d) { + return std::isnan(d) ? "?" : to_str(d); +} + +// for use with non-negative Metadata fields that use -1 for N/A +inline std::string int_or_dot(int n) { + return n == -1 ? "." : std::to_string(n); +} +inline std::string int_or_qmark(int n) { + return n == -1 ? "?" : std::to_string(n); +} + +inline std::string string_or_dot(const std::string& s) { + return s.empty() ? "." : cif::quote(s); +} +inline std::string string_or_qmark(const std::string& s) { + return s.empty() ? "?" : cif::quote(s); +} + +// Quote chain name or entity id if necessary. It is necessary +// only if the chain name is missing, which was OK in the past. +// Here we use '' rather than . or ?. +inline std::string qchain(const std::string& s) { + return cif::quote(s); +} + + +void add_cif_atoms(const Structure& st, cif::Block& block, + bool use_group_pdb, bool auth_all) { + // atom list + cif::Loop& atom_loop = block.init_mmcif_loop("_atom_site.", { + "id", + "type_symbol", + "label_atom_id", + "label_alt_id", + "label_comp_id", + "label_asym_id", + "label_entity_id", + "label_seq_id", + "pdbx_PDB_ins_code", + "Cartn_x", + "Cartn_y", + "Cartn_z", + "occupancy", + "B_iso_or_equiv", + "pdbx_formal_charge", + "auth_atom_id", // optional (tags[15] is removed if !auth_all) + "auth_comp_id", // optional (tags[16] is removed if !auth_all) + "auth_seq_id", + "auth_asym_id", + "pdbx_PDB_model_num"}); + if (!auth_all) + atom_loop.tags.erase(atom_loop.tags.begin() + 15, atom_loop.tags.begin() + 17); + if (use_group_pdb) + atom_loop.tags.emplace(atom_loop.tags.begin(), "_atom_site.group_PDB"); + bool has_calc_flag = false; + bool has_tls_group_id = false; + size_t atom_site_count = 0; + for (const Model& model : st.models) + for (const Chain& chain : model.chains) + for (const Residue& res : chain.residues) + for (const Atom& atom : res.atoms) { + ++atom_site_count; + if (atom.calc_flag != CalcFlag::NotSet && + atom.calc_flag != CalcFlag::NoHydrogen) + has_calc_flag = true; + if (atom.tls_group_id >= 0) + has_tls_group_id = true; + } + if (has_calc_flag) + atom_loop.tags.emplace_back("_atom_site.calc_flag"); + if (has_tls_group_id) + atom_loop.tags.emplace_back("_atom_site.pdbx_tls_group_id"); + if (st.has_d_fraction) + atom_loop.tags.emplace_back("_atom_site.ccp4_deuterium_fraction"); + + std::vector& vv = atom_loop.values; + vv.reserve(atom_site_count * atom_loop.tags.size()); + std::vector> aniso; + int serial = 0; + for (const Model& model : st.models) { + for (const Chain& chain : model.chains) { + for (const Residue& res : chain.residues) { + bool as_het = use_hetatm(res); + std::string label_seq_id = res.label_seq.str('.'); + std::string auth_seq_id = res.seqid.num.str(); + std::string entity_id; + if (const Entity* ent = gemmi::find_entity_of_subchain(res.subchain, st.entities)) + entity_id = cif::quote(ent->name); + else + entity_id = string_or_dot(res.entity_id); + for (const Atom& atom : res.atoms) { + if (use_group_pdb) + vv.emplace_back(as_het ? "HETATM" : "ATOM"); + vv.emplace_back(std::to_string(++serial)); + vv.emplace_back(atom.element.uname()); + vv.emplace_back(cif::quote(atom.name)); + vv.emplace_back(1, atom.altloc_or('.')); + vv.emplace_back(cif::quote(res.name)); + vv.emplace_back(subchain_or_dot(res)); + vv.emplace_back(entity_id); + vv.emplace_back(label_seq_id); + vv.emplace_back(pdbx_icode(res)); + vv.emplace_back(to_str(atom.pos.x)); + vv.emplace_back(to_str(atom.pos.y)); + vv.emplace_back(to_str(atom.pos.z)); + vv.emplace_back(to_str(atom.occ)); + vv.emplace_back(to_str(atom.b_iso)); + vv.emplace_back(atom.charge == 0 ? "?" : std::to_string(atom.charge)); + if (auth_all) { + size_t atom_name_idx = vv.size() - 13; + vv.emplace_back(vv[atom_name_idx]); // auth_atom_id = label_atom_id + vv.emplace_back(vv[atom_name_idx + 2]); // auth_comp_id = label_comp_id + } + vv.emplace_back(auth_seq_id); + vv.emplace_back(qchain(chain.name)); + vv.emplace_back(std::to_string(model.num)); + if (has_calc_flag) + vv.emplace_back(&".\0.\0d\0c\0dum"[2 * (int) atom.calc_flag]); + if (has_tls_group_id) + vv.emplace_back(int_or_qmark(atom.tls_group_id)); + if (st.has_d_fraction) + vv.emplace_back(to_str(atom.fraction)); + if (atom.aniso.nonzero()) + aniso.emplace_back(serial, model.num, &atom); + } + } + } + } + if (aniso.empty()) { + block.find_mmcif_category("_atom_site_anisotrop.").erase(); + } else { + cif::Loop& aniso_loop = block.init_mmcif_loop("_atom_site_anisotrop.", { + "id", "type_symbol", "U[1][1]", "U[2][2]", + "U[3][3]", "U[1][2]", "U[1][3]", "U[2][3]"}); + if (st.models.size() > 1) + aniso_loop.tags.push_back("_atom_site_anisotrop.pdbx_PDB_model_num"); + std::vector& aniso_val = aniso_loop.values; + aniso_val.reserve(aniso_loop.tags.size() * aniso.size()); + for (const auto& a : aniso) { + aniso_val.emplace_back(std::to_string(std::get<0>(a))); + const Atom* atom = std::get<2>(a); + aniso_val.emplace_back(atom->element.uname()); + aniso_val.emplace_back(to_str(atom->aniso.u11)); + aniso_val.emplace_back(to_str(atom->aniso.u22)); + aniso_val.emplace_back(to_str(atom->aniso.u33)); + aniso_val.emplace_back(to_str(atom->aniso.u12)); + aniso_val.emplace_back(to_str(atom->aniso.u13)); + aniso_val.emplace_back(to_str(atom->aniso.u23)); + if (st.models.size() > 1) + aniso_loop.values.push_back(std::to_string(std::get<1>(a))); + } + } +} + +// the names are: monomeric, dimeric, ...meric, 21-meric, 22-meric, ... +int xmeric_to_number(const std::string& oligomeric) { + static const char names[20][10] = { + "mono", "di", "tri", "tetra", "penta", + "hexa", "hepta", "octa", "nona", "deca", + "undeca", "dodeca", "trideca", "tetradeca", "pentadeca", + "hexadeca", "heptadeca", "octadeca", "nonadeca", "eicosa" + }; + size_t len = oligomeric.length(); + const char* p = oligomeric.c_str(); + for (int i = 0; i != 20; ++i) + if (len == std::strlen(names[i]) + 5 && strncmp(p, names[i], len-5) == 0) + return i + 1; + return no_sign_atoi(p); +} + +void write_assemblies(const Structure& st, cif::Block& block) { + block.items.reserve(block.items.size() + 4); // avoid re-allocation + cif::Loop& a_loop = block.init_mmcif_loop("_pdbx_struct_assembly.", + {"id", "details", "method_details", + "oligomeric_details", "oligomeric_count"}); + cif::Loop& prop_loop = block.init_mmcif_loop("_pdbx_struct_assembly_prop.", + {"biol_id", "type", "value"}); + cif::Loop& gen_loop = block.init_mmcif_loop("_pdbx_struct_assembly_gen.", + {"assembly_id", "oper_expression", "asym_id_list"}); + cif::Loop& oper_loop = block.init_mmcif_loop("_pdbx_struct_oper_list.", + {"id", "type", + "matrix[1][1]", "matrix[1][2]", "matrix[1][3]", "vector[1]", + "matrix[2][1]", "matrix[2][2]", "matrix[2][3]", "vector[2]", + "matrix[3][1]", "matrix[3][2]", "matrix[3][3]", "vector[3]"}); + std::vector distinct_oper; + for (const Assembly& as : st.assemblies) { + std::string how_defined = "?"; + if (as.author_determined && as.software_determined) + how_defined = "author_and_software_defined_assembly"; + else if (as.author_determined) + how_defined = "author_defined_assembly"; + else if (as.software_determined) + how_defined = "software_defined_assembly"; + else if (as.special_kind == Assembly::SpecialKind::CompleteIcosahedral) + how_defined = "'complete icosahedral assembly'"; + else if (as.special_kind == Assembly::SpecialKind::RepresentativeHelical) + how_defined = "'representative helical assembly'"; + else if (as.special_kind == Assembly::SpecialKind::CompletePoint) + how_defined = "'complete point assembly'"; + std::string oligomer = to_lower(as.oligomeric_details); + int nmer = as.oligomeric_count != 0 ? as.oligomeric_count + : xmeric_to_number(oligomer); + // _pdbx_struct_assembly + a_loop.add_row({as.name, + how_defined, + string_or_qmark(as.software_name), + string_or_qmark(oligomer), + nmer == 0 ? "?" : std::to_string(nmer)}); + + // _pdbx_struct_assembly_prop + if (!std::isnan(as.absa)) + prop_loop.add_row({as.name, "'ABSA (A^2)'", to_str(as.absa)}); + if (!std::isnan(as.ssa)) + prop_loop.add_row({as.name, "'SSA (A^2)'", to_str(as.ssa)}); + if (!std::isnan(as.more)) + prop_loop.add_row({as.name, "MORE", to_str(as.more)}); + + // _pdbx_struct_assembly_gen and _pdbx_struct_oper_list + for (const Assembly::Gen& gen : as.generators) { + std::string subchain_str; + for (const std::string& name : gen.subchains) + string_append_sep(subchain_str, ',', name); + if (subchain_str.empty()) // chain names to subchain names + for (const Chain& chain : st.models[0].chains) + if (in_vector(chain.name, gen.chains)) + for (const auto& sub : chain.subchains()) + string_append_sep(subchain_str, ',', sub.front().subchain); + std::string oper_str; + for (const Assembly::Operator& oper : gen.operators) { + size_t k = 0; + for (; k != distinct_oper.size(); ++k) + if (distinct_oper[k]->transform.approx(oper.transform, 1e-9)) + break; + string_append_sep(oper_str, ',', std::to_string(k+1)); + if (k != distinct_oper.size()) + continue; + distinct_oper.emplace_back(&oper); + oper_loop.values.emplace_back(std::to_string(k+1)); + if (!oper.type.empty()) { + oper_loop.values.emplace_back(cif::quote(oper.type)); + } else if (oper.transform.is_identity()) { + oper_loop.values.emplace_back("'identity operation'"); + } else if (as.author_determined || as.software_determined) { + oper_loop.values.emplace_back("'crystal symmetry operation'"); + } else { + oper_loop.values.emplace_back("."); + } + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) + oper_loop.values.emplace_back(to_str(oper.transform.mat[i][j])); + oper_loop.values.emplace_back(to_str(oper.transform.vec.at(i))); + } + } + gen_loop.add_row({as.name, + oper_str.empty() ? "." : oper_str, + subchain_str.empty() ? "?" : subchain_str}); + } + } +} + +void write_cell_parameters(const UnitCell& cell, cif::ItemSpan& span) { + span.set_pair("_cell.length_a", to_str(cell.a)); + span.set_pair("_cell.length_b", to_str(cell.b)); + span.set_pair("_cell.length_c", to_str(cell.c)); + span.set_pair("_cell.angle_alpha", to_str(cell.alpha)); + span.set_pair("_cell.angle_beta", to_str(cell.beta)); + span.set_pair("_cell.angle_gamma", to_str(cell.gamma)); +} + +bool is_valid_block_name(const std::string& name) { + return !name.empty() && + std::all_of(name.begin(), name.end(), [](char c){ return c >= '!' && c <= '~'; }); +} + +int get_number_obs(const BasicRefinementInfo& ref) { + int nobs = ref.reflection_count; + if (nobs == -1 && ref.rfree_set_count >= 0 && ref.work_set_count >= 0) + nobs = ref.work_set_count + ref.rfree_set_count; + return nobs; +} + +int get_number_work(const BasicRefinementInfo& ref) { + int nwork = ref.work_set_count; + if (nwork == -1 && ref.rfree_set_count >= 0 && ref.reflection_count >= 0) + nwork = ref.reflection_count - ref.rfree_set_count; + return nwork; +} + +} // anonymous namespace + +void write_ncs_oper(const Structure& st, cif::Block& block) { + // _struct_ncs_oper (MTRIX) + if (st.ncs.empty()) + return; + cif::Loop& ncs_oper = block.init_mmcif_loop("_struct_ncs_oper.", + {"id", "code", + "matrix[1][1]", "matrix[1][2]", "matrix[1][3]", "vector[1]", + "matrix[2][1]", "matrix[2][2]", "matrix[2][3]", "vector[2]", + "matrix[3][1]", "matrix[3][2]", "matrix[3][3]", "vector[3]"}); + auto add_op = [&ncs_oper](const NcsOp& op) { + ncs_oper.values.emplace_back(op.id); + ncs_oper.values.emplace_back(op.given ? "given" : "generate"); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) + ncs_oper.values.emplace_back(to_str(op.tr.mat[i][j])); + ncs_oper.values.emplace_back(to_str(op.tr.vec.at(i))); + } + }; + auto identity = st.info.find("_struct_ncs_oper.id"); + if (identity != st.info.end() && + !in_vector_f([&](const NcsOp& op) { return op.id == identity->second; }, st.ncs)) + add_op(NcsOp{identity->second, true, {}}); + for (const NcsOp& op : st.ncs) + add_op(op); +} + +void write_struct_conn(const Structure& st, cif::Block& block) { + // example: + // disulf1 disulf A CYS 3 SG ? 3 ? 1_555 A CYS 18 SG ? 18 ? 1_555 ? 2.045 + std::array type_ids{}; + bool use_ccp4_link_id = false; + for (const Connection& con : st.connections) + if (!con.link_id.empty()) + use_ccp4_link_id = true; + cif::Loop& conn_loop = block.init_mmcif_loop("_struct_conn.", + {"id", "conn_type_id", + "ptnr1_label_asym_id", "ptnr1_label_comp_id", "ptnr1_label_seq_id", + "ptnr1_label_atom_id", "pdbx_ptnr1_label_alt_id", "ptnr1_auth_asym_id", + "ptnr1_auth_seq_id", "pdbx_ptnr1_PDB_ins_code", "ptnr1_symmetry", + "ptnr2_label_asym_id", "ptnr2_label_comp_id", "ptnr2_label_seq_id", + "ptnr2_label_atom_id", "pdbx_ptnr2_label_alt_id", "ptnr2_auth_asym_id", + "ptnr2_auth_seq_id", "pdbx_ptnr2_PDB_ins_code", "ptnr2_symmetry", + "details", "pdbx_dist_value"}); + if (use_ccp4_link_id) + conn_loop.tags.push_back("_struct_conn.ccp4_link_id"); + for (const Connection& con : st.connections) { + const_CRA cra1 = st.models[0].find_cra(con.partner1, true); + const_CRA cra2 = st.models[0].find_cra(con.partner2, true); + if (!cra1.residue || !cra2.residue) + continue; + const Atom* at1 = cra1.atom; + const Atom* at2 = cra2.atom; + std::string im_pdb_symbol = "?", im_dist_str = "?"; + if (at1 && at2) { + NearestImage im = st.cell.find_nearest_image(at1->pos, at2->pos, con.asu); + im_pdb_symbol = im.symmetry_code(true); + im_dist_str = to_str_prec<4>(im.dist()); + } + auto& v = conn_loop.values; + v.emplace_back(string_or_qmark(con.name)); // id + v.emplace_back(connection_type_to_string(con.type)); // conn_type_id + v.emplace_back(subchain_or_dot(*cra1.residue)); // ptnr1_label_asym_id + v.emplace_back(cra1.residue->name); // ptnr1_label_comp_id + v.emplace_back(cra1.residue->label_seq.str('.')); // ptnr1_label_seq_id + v.emplace_back(at1 ? cif::quote(at1->name) : "?"); // ptnr1_label_atom_id + v.emplace_back(1, at1 ? at1->altloc_or('?') : '?'); // pdbx_ptnr1_label_alt_id + v.emplace_back(qchain(con.partner1.chain_name)); // ptnr1_auth_asym_id + v.emplace_back(cra1.residue->seqid.num.str()); // ptnr1_auth_seq_id + v.emplace_back(pdbx_icode(con.partner1.res_id)); // ptnr1_PDB_ins_code + v.emplace_back("1_555"); // ptnr1_symmetry + v.emplace_back(subchain_or_dot(*cra2.residue)); // ptnr2_label_asym_id + v.emplace_back(cra2.residue->name); // ptnr2_label_comp_id + v.emplace_back(cra2.residue->label_seq.str('.')); // ptnr2_label_seq_id + v.emplace_back(at2 ? cif::quote(at2->name) : "?"); // ptnr2_label_atom_id + v.emplace_back(1, at2 ? at2->altloc_or('?') : '?'); // pdbx_ptnr2_label_alt_id + v.emplace_back(qchain(con.partner2.chain_name)); // ptnr2_auth_asym_id + v.emplace_back(cra2.residue->seqid.num.str()); // ptnr2_auth_seq_id + v.emplace_back(pdbx_icode(con.partner2.res_id)); // ptnr2_PDB_ins_code + v.emplace_back(im_pdb_symbol); // ptnr2_symmetry + v.emplace_back("?"); // details + v.emplace_back(im_dist_str); // pdbx_dist_value + if (use_ccp4_link_id) + v.emplace_back(string_or_qmark(con.link_id)); // ccp4_link_id + type_ids[int(con.type)] = true; + } + + cif::Loop& type_loop = block.init_mmcif_loop("_struct_conn_type.", {"id"}); + for (int i = 0; i < (int)type_ids.size() - 1; ++i) + if (type_ids[i]) + type_loop.add_row({connection_type_to_string((Connection::Type)i)}); +} + +void write_cispeps(const Structure& st, cif::Block& block) { + cif::Loop* prot_cis_loop = nullptr; + int pdbx_id = 0; + for (const CisPep& cispep : st.cispeps) { + const Model* model = &st.models[0]; + if (st.models.size() > 1) { + model = st.find_model(cispep.model_num); + if (!model) + continue; + } + const_CRA cra1 = model->find_cra(cispep.partner_c, true); + const_CRA cra2 = model->find_cra(cispep.partner_n, true); + if (!cra1.residue || !cra2.residue) + continue; + if (!prot_cis_loop) + prot_cis_loop = &block.init_mmcif_loop("_struct_mon_prot_cis.", + {"pdbx_id", "pdbx_PDB_model_num", + "label_asym_id", "label_seq_id", "label_comp_id", + "auth_asym_id", "auth_seq_id", "pdbx_PDB_ins_code", + "pdbx_label_asym_id_2", "pdbx_label_seq_id_2", "pdbx_label_comp_id_2", + "pdbx_auth_asym_id_2", "pdbx_auth_seq_id_2", "pdbx_PDB_ins_code_2", + "label_alt_id", "pdbx_omega_angle"}); + auto& v = prot_cis_loop->values; + v.emplace_back(std::to_string(++pdbx_id)); // pdbx_id + v.emplace_back(std::to_string(model->num)); // pdbx_PDB_model_num + v.emplace_back(subchain_or_dot(*cra1.residue)); // label_asym_id + v.emplace_back(cra1.residue->label_seq.str('.')); // label_seq_id + v.emplace_back(cra1.residue->name); // label_comp_id + v.emplace_back(qchain(cispep.partner_c.chain_name)); // auth_asym_id + v.emplace_back(cispep.partner_c.res_id.seqid.num.str()); // auth_seq_id + v.emplace_back(pdbx_icode(cispep.partner_c.res_id)); // pdbx_PDB_ins_code + v.emplace_back(subchain_or_dot(*cra2.residue)); // pdbx_label_asym_id_2 + v.emplace_back(cra2.residue->label_seq.str('.')); // pdbx_label_seq_id_2 + v.emplace_back(cra2.residue->name); // pdbx_label_comp_id_2 + v.emplace_back(qchain(cispep.partner_n.chain_name)); // pdbx_auth_asym_id_2 + v.emplace_back(cispep.partner_n.res_id.seqid.num.str()); // pdbx_auth_seq_id_2 + v.emplace_back(pdbx_icode(cispep.partner_n.res_id)); // pdbx_PDB_ins_code_2 + v.emplace_back(1, cispep.only_altloc ? cispep.only_altloc : '.'); + v.emplace_back(number_or_qmark(cispep.reported_angle)); + } +} + +void update_mmcif_block(const Structure& st, cif::Block& block, MmcifOutputGroups groups) { + if (st.models.empty()) + return; + + if (groups.block_name) + block.name = is_valid_block_name(st.name) ? st.name : "model"; + + auto e_id = st.info.find("_entry.id"); + std::string id = cif::quote(e_id != st.info.end() ? e_id->second : block.name); + if (groups.entry) + block.set_pair("_entry.id", id); + else if (const std::string* val = block.find_value("_entry.id")) + id = *val; + + if (groups.database_status) { + auto initial_date = st.info.find("_pdbx_database_status.recvd_initial_deposition_date"); + if (initial_date != st.info.end() && !initial_date->second.empty()) { + cif::ItemSpan span(block.items, "_pdbx_database_status."); + span.set_pair("_pdbx_database_status.entry_id", id); + span.set_pair(initial_date->first, initial_date->second); + } + } + + if (groups.author && !st.meta.authors.empty()) { + cif::Loop& loop = block.init_mmcif_loop("_audit_author.", {"pdbx_ordinal", "name"}); + int n = 0; + for (const std::string& author : st.meta.authors) + loop.add_row({std::to_string(++n), cif::quote(author)}); + } + + if (groups.cell) { + cif::ItemSpan cell_span(block.items, "_cell."); + cell_span.set_pair("_cell.entry_id", id); + write_cell_parameters(st.cell, cell_span); + auto z_pdb = st.info.find("_cell.Z_PDB"); + if (z_pdb != st.info.end()) + cell_span.set_pair(z_pdb->first, z_pdb->second); + } + + if (groups.symmetry) { + cif::ItemSpan span(block.items, "_symmetry."); + span.set_pair("_symmetry.entry_id", id); + span.set_pair("_symmetry.space_group_name_H-M", + cif::quote(st.spacegroup_hm)); + if (const SpaceGroup* sg = st.find_spacegroup()) + span.set_pair("_symmetry.Int_Tables_number", std::to_string(sg->number)); + } + + if (groups.entity) { + cif::Loop& entity_loop = block.init_mmcif_loop("_entity.", {"id", "type"}); + for (const Entity& ent : st.entities) + entity_loop.add_row({qchain(ent.name), + entity_type_to_string(ent.entity_type)}); + } + + std::map subs_to_strands; + if (groups.entity_poly || groups.struct_ref) + subs_to_strands = st.models[0].subchain_to_chain(); + + if (groups.entity_poly) { + // If the _entity_poly category is included when depositing to the PDB, + // it must contain entity_id, type, pdbx_seq_one_letter_code + // and pdbx_strand_id. The last one is not documented as required, + // but OneDep shows error when it's not included. + cif::Loop& ent_poly_loop = block.init_mmcif_loop("_entity_poly.", + {"entity_id", "type", "pdbx_strand_id", "pdbx_seq_one_letter_code"}); + for (const Entity& ent : st.entities) + if (ent.entity_type == EntityType::Polymer) { + if (ent.polymer_type == PolymerType::Unknown) + continue; // not sure what to do here + ResidueKind kind = sequence_kind(ent.polymer_type); + std::string seq1 = pdbx_one_letter_code(ent.full_sequence, kind); + std::string strand_ids; + for (const std::string& sub : ent.subchains) { + auto strand_id = subs_to_strands.find(sub); + if (strand_id != subs_to_strands.end()) { + if (!strand_ids.empty()) + strand_ids += ','; + strand_ids += strand_id->second; + } + } + ent_poly_loop.add_row({qchain(ent.name), + polymer_type_to_string(ent.polymer_type), + string_or_qmark(strand_ids), + string_or_qmark(seq1)}); + } + } + + if (groups.struct_ref) { // _struct_ref, _struct_ref_seq + block.items.reserve(block.items.size() + 2); // avoid re-allocation + cif::Loop& ref_loop = block.init_mmcif_loop("_struct_ref.", + {"id", "entity_id", "db_name", "db_code", + "pdbx_db_accession", "pdbx_db_isoform"}); + cif::Loop& seq_loop = block.init_mmcif_loop("_struct_ref_seq.", { + "align_id", "ref_id", "pdbx_strand_id", "pdbx_PDB_id_code", + "seq_align_beg", "seq_align_end", "pdbx_db_accession", + "db_align_beg", "db_align_end", + "pdbx_auth_seq_align_beg", "pdbx_seq_align_beg_ins_code", + "pdbx_auth_seq_align_end", "pdbx_seq_align_end_ins_code"}); + int counter = 0; + int counter2 = 0; + for (const Entity& ent : st.entities) + for (const Entity::DbRef& dbref : ent.dbrefs) { + ref_loop.add_row({std::to_string(++counter), + qchain(ent.name), + string_or_dot(dbref.db_name), + string_or_dot(dbref.id_code), + string_or_qmark(dbref.accession_code), + string_or_qmark(dbref.isoform)}); + for (const std::string& subchain : ent.subchains) { + auto strand_id = subs_to_strands.find(subchain); + if (strand_id == subs_to_strands.end()) + continue; + // DbRef::label_seq_begin/end (_struct_ref_seq.seq_align_beg/end) is + // not filled in when reading PDB file, so we check it here. + Residue::OptionalNum label_begin = dbref.label_seq_begin; + Residue::OptionalNum label_end = dbref.label_seq_end; + if (!label_begin || !label_end) { + ConstResidueSpan span = st.models[0].get_subchain(subchain); + try { + label_begin = span.auth_seq_id_to_label(dbref.seq_begin); + label_end = span.auth_seq_id_to_label(dbref.seq_end); + } catch (const std::out_of_range&) {} + } + SeqId begin = dbref.seq_begin; + SeqId end = dbref.seq_end; + if (!begin.num || !end.num) { + if (const Chain* chain = st.models[0].find_chain(strand_id->second)) + if (ConstResidueGroup polymer = chain->get_polymer()) { + begin = polymer.label_seq_id_to_auth(dbref.label_seq_begin); + end = polymer.label_seq_id_to_auth(dbref.label_seq_end); + } + } + seq_loop.add_row({std::to_string(++counter2), + std::to_string(counter), + strand_id->second, // pdbx_strand_id + id, + label_begin.str(), + label_end.str(), + string_or_qmark(dbref.accession_code), + dbref.db_begin.num.str(), + dbref.db_end.num.str(), + begin.num.str(), + pdbx_icode(begin), + end.num.str(), + pdbx_icode(end)}); + } + } + } + + if (groups.chem_comp) { + std::set resnames; + for (const Model& model : st.models) + for (const Chain& chain : model.chains) + for (const Residue& res : chain.residues) + resnames.insert(res.name); + for (const Entity& ent : st.entities) + for (const std::string& item : ent.full_sequence) + resnames.insert(Entity::first_mon(item)); + cif::Loop& chem_comp_loop = block.init_mmcif_loop("_chem_comp.", {"id", "type"}); + if (!st.shortened_ccd_codes.empty()) + chem_comp_loop.tags.push_back("_chem_comp.three_letter_code"); + for (const std::string& name : resnames) { + chem_comp_loop.values.push_back(cif::quote(name)); + chem_comp_loop.values.push_back("."); + if (!st.shortened_ccd_codes.empty()) { + chem_comp_loop.values.push_back(cif::quote(name)); + for (const auto& old_new : st.shortened_ccd_codes) + if (old_new.second == name) + chem_comp_loop.values.back() = old_new.first; + } + } + } + + if (groups.exptl) { + // _exptl + if (!st.meta.experiments.empty()) { + cif::Loop& loop = block.init_mmcif_loop("_exptl.", + {"entry_id", "method", "crystals_number"}); + for (const ExperimentInfo& exper : st.meta.experiments) + loop.add_row({id, cif::quote(exper.method), + int_or_qmark(exper.number_of_crystals)}); + } else { + auto exptl_method = st.info.find("_exptl.method"); + if (exptl_method != st.info.end()) { + cif::Loop& loop = block.init_mmcif_loop("_exptl.", {"entry_id", "method"}); + for (const std::string& m : gemmi::split_str(exptl_method->second, "; ")) + loop.add_row({id, cif::quote(m)}); + } + } + + // _exptl_crystal + if (!st.meta.crystals.empty()) { + cif::Loop& loop = block.init_mmcif_loop("_exptl_crystal.", + {"id", "description"}); + for (const CrystalInfo& cryst : st.meta.crystals) + loop.add_row({cryst.id, string_or_qmark(cryst.description)}); + } + + // _exptl_crystal_grow + if (std::any_of(st.meta.crystals.begin(), st.meta.crystals.end(), + [](const CrystalInfo& c) { return !c.ph_range.empty() || !std::isnan(c.ph); })) { + cif::Loop& grow_loop = block.init_mmcif_loop("_exptl_crystal_grow.", + {"crystal_id", "pH", "pdbx_pH_range"}); + for (const CrystalInfo& crystal : st.meta.crystals) + grow_loop.add_row({cif::quote(crystal.id), + number_or_qmark(crystal.ph), + string_or_qmark(crystal.ph_range)}); + } + } + + if (groups.diffrn && + std::any_of(st.meta.crystals.begin(), st.meta.crystals.end(), + [](const CrystalInfo& c) { return !c.diffractions.empty(); })) { + + cif::Loop& loop = block.init_mmcif_loop("_diffrn.", {"id", "crystal_id", "ambient_temp"}); + for (const CrystalInfo& cryst : st.meta.crystals) + for (const DiffractionInfo& diffr : cryst.diffractions) + loop.add_row({diffr.id, cryst.id, number_or_qmark(diffr.temperature)}); + // _diffrn_detector + cif::Loop& det_loop = block.init_mmcif_loop("_diffrn_detector.", + {"diffrn_id", + "pdbx_collection_date", + "detector", + "type", + "details"}); + for (const CrystalInfo& cryst : st.meta.crystals) + for (const DiffractionInfo& diffr : cryst.diffractions) + det_loop.add_row({diffr.id, + string_or_qmark(diffr.collection_date), + string_or_qmark(diffr.detector), + string_or_qmark(diffr.detector_make), + string_or_qmark(diffr.optics)}); + + // _diffrn_radiation + cif::Loop& rad_loop = block.init_mmcif_loop("_diffrn_radiation.", + {"diffrn_id", + "pdbx_scattering_type", + "pdbx_monochromatic_or_laue_m_l", + "monochromator"}); + for (const CrystalInfo& cryst : st.meta.crystals) + for (const DiffractionInfo& diffr : cryst.diffractions) + rad_loop.add_row({diffr.id, + string_or_qmark(diffr.scattering_type), + std::string(1, diffr.mono_or_laue ? diffr.mono_or_laue : '?'), + string_or_qmark(diffr.monochromator)}); + // _diffrn_source + cif::Loop& source_loop = block.init_mmcif_loop("_diffrn_source.", + {"diffrn_id", + "source", + "type", + "pdbx_synchrotron_site", + "pdbx_synchrotron_beamline", + "pdbx_wavelength_list"}); + for (const CrystalInfo& crystal : st.meta.crystals) + for (const DiffractionInfo& diffr : crystal.diffractions) + source_loop.add_row({diffr.id, + string_or_qmark(diffr.source), + string_or_qmark(diffr.source_type), + string_or_qmark(diffr.synchrotron), + string_or_qmark(diffr.beamline), + string_or_qmark(diffr.wavelengths)}); + } + + if (groups.reflns && !st.meta.experiments.empty()) { + // _reflns + cif::Loop& loop = block.init_mmcif_loop("_reflns.", { + "entry_id", + "pdbx_ordinal", + "pdbx_diffrn_id", + "number_obs", + "d_resolution_high", + "d_resolution_low", + "percent_possible_obs", + "pdbx_redundancy", + "pdbx_Rmerge_I_obs", + "pdbx_Rsym_value", + "pdbx_netI_over_sigmaI", + /*"B_iso_Wilson_estimate"*/}); + int n = 0; + for (const ExperimentInfo& exper : st.meta.experiments) + loop.add_row({id, + std::to_string(++n), + string_or_dot(join_str(exper.diffraction_ids, ",")), + int_or_qmark(exper.unique_reflections), + number_or_qmark(exper.reflections.resolution_high), + number_or_qmark(exper.reflections.resolution_low), + number_or_qmark(exper.reflections.completeness), + number_or_qmark(exper.reflections.redundancy), + number_or_qmark(exper.reflections.r_merge), + number_or_qmark(exper.reflections.r_sym), + number_or_qmark(exper.reflections.mean_I_over_sigma), + /*number_or_qmark(exper.b_wilson)*/}); + // _reflns_shell + cif::Loop* shell_loop = nullptr; + n = 0; + for (const ExperimentInfo& exper : st.meta.experiments) { + std::string diffrn_id = + string_or_dot(join_str(exper.diffraction_ids, ",")); + for (const ReflectionsInfo& shell : exper.shells) { + if (!shell_loop) + shell_loop = &block.init_mmcif_loop("_reflns_shell.", { + "pdbx_ordinal", + "pdbx_diffrn_id", + "d_res_high", + "d_res_low", + "percent_possible_all", + "pdbx_redundancy", + "Rmerge_I_obs", + "pdbx_Rsym_value", + "meanI_over_sigI_obs"}); + + shell_loop->add_row({std::to_string(++n), + diffrn_id, + number_or_qmark(shell.resolution_high), + number_or_qmark(shell.resolution_low), + number_or_qmark(shell.completeness), + number_or_qmark(shell.redundancy), + number_or_qmark(shell.r_merge), + number_or_qmark(shell.r_sym), + number_or_qmark(shell.mean_I_over_sigma)}); + } + } + } + + if (groups.refine && !st.meta.refinement.empty()) { + block.items.reserve(block.items.size() + 4); + cif::Loop& loop = block.init_mmcif_loop("_refine.", { + "entry_id", + "pdbx_refine_id", + "ls_d_res_high", + "ls_d_res_low", + "ls_percent_reflns_obs", + "ls_number_reflns_obs", + "ls_number_reflns_R_work"}); + cif::Loop& analyze_loop = block.init_mmcif_loop("_refine_analyze.", { + "entry_id", + "pdbx_refine_id", + "Luzzati_coordinate_error_obs"}); + cif::Loop& restr_loop = block.init_mmcif_loop("_refine_ls_restr.", { + "pdbx_refine_id", "type", + "number", "weight", "pdbx_restraint_function", "dev_ideal"}); + // _refine_ls_shell + std::vector shell_tags = { + "pdbx_refine_id", + "d_res_high", + "d_res_low", + "percent_reflns_obs", + "number_reflns_obs", + "number_reflns_R_work", + "number_reflns_R_free", + "R_factor_obs", + "R_factor_R_work", + "R_factor_R_free"}; + bool has_shell_fsc = false; + bool has_shell_ffcc = false; + bool has_shell_iicc = false; + for (const RefinementInfo& ref : st.meta.refinement) + for (const BasicRefinementInfo& bin : ref.bins) { + if (!std::isnan(bin.fsc_work) || !std::isnan(bin.fsc_free)) + has_shell_fsc = true; + if (!std::isnan(bin.cc_fo_fc_work) || !std::isnan(bin.cc_fo_fc_free)) + has_shell_ffcc = true; + if (!std::isnan(bin.cc_intensity_work) || !std::isnan(bin.cc_intensity_free)) + has_shell_iicc = true; + } + if (has_shell_fsc) { + shell_tags.push_back("pdbx_fsc_work"); + shell_tags.push_back("pdbx_fsc_free"); + } + if (has_shell_ffcc) { + shell_tags.push_back("correlation_coeff_Fo_to_Fc"); + shell_tags.push_back("correlation_coeff_Fo_to_Fc_free"); + } + if (has_shell_iicc) { + shell_tags.push_back("correlation_coeff_I_to_Fcsqd_work"); + shell_tags.push_back("correlation_coeff_I_to_Fcsqd_free"); + } + cif::Loop& shell_loop = block.init_mmcif_loop("_refine_ls_shell.", shell_tags); + + for (size_t i = 0; i != st.meta.refinement.size(); ++i) { + const RefinementInfo& ref = st.meta.refinement[i]; + loop.add_values({id, + cif::quote(ref.id), + number_or_dot(ref.resolution_high), + number_or_dot(ref.resolution_low), + number_or_dot(ref.completeness), + int_or_dot(get_number_obs(ref)), + int_or_qmark(get_number_work(ref))}); + auto add = [&](const std::string& tag, const std::string& val) { + if (i == 0) + loop.tags.push_back("_refine." + tag); + loop.values.push_back(val); + }; + if (st.meta.has(&RefinementInfo::rfree_set_count)) + add("ls_number_reflns_R_free", int_or_dot(ref.rfree_set_count)); + if (st.meta.has(&RefinementInfo::r_all)) + add("ls_R_factor_obs", number_or_qmark(ref.r_all)); + if (st.meta.has(&RefinementInfo::r_work)) + add("ls_R_factor_R_work", number_or_qmark(ref.r_work)); + if (st.meta.has(&RefinementInfo::r_free)) + add("ls_R_factor_R_free", number_or_qmark(ref.r_free)); + if (st.meta.has(&RefinementInfo::cross_validation_method)) + add("pdbx_ls_cross_valid_method", + string_or_qmark(ref.cross_validation_method)); + if (st.meta.has(&RefinementInfo::rfree_selection_method)) + add("pdbx_R_Free_selection_details", + string_or_qmark(ref.rfree_selection_method)); + if (st.meta.has(&RefinementInfo::mean_b)) + add("B_iso_mean", number_or_qmark(ref.mean_b)); + if (st.meta.has(&RefinementInfo::aniso_b)) { + if (i == 0) + for (const char* index : {"[1][1]", "[2][2]", "[3][3]", "[1][2]", "[1][3]", "[2][3]"}) + loop.tags.push_back(std::string("_refine.aniso_B") + index); + for (double d : ref.aniso_b.elements_pdb()) + loop.values.push_back(number_or_qmark(d)); + } + if (st.meta.has(&RefinementInfo::dpi_blow_r)) + add("pdbx_overall_SU_R_Blow_DPI", + number_or_qmark(ref.dpi_blow_r)); + if (st.meta.has(&RefinementInfo::dpi_blow_rfree)) + add("pdbx_overall_SU_R_free_Blow_DPI", + number_or_qmark(ref.dpi_blow_rfree)); + if (st.meta.has(&RefinementInfo::dpi_cruickshank_r)) + add("overall_SU_R_Cruickshank_DPI", + number_or_qmark(ref.dpi_cruickshank_r)); + if (st.meta.has(&RefinementInfo::dpi_cruickshank_rfree)) + add("pdbx_overall_SU_R_free_Cruickshank_DPI", + number_or_qmark(ref.dpi_cruickshank_rfree)); + if (st.meta.has(&RefinementInfo::cc_fo_fc_work)) + add("correlation_coeff_Fo_to_Fc", number_or_qmark(ref.cc_fo_fc_work)); + if (st.meta.has(&RefinementInfo::cc_fo_fc_free)) + add("correlation_coeff_Fo_to_Fc_free", + number_or_qmark(ref.cc_fo_fc_free)); + if (st.meta.has(&RefinementInfo::fsc_work)) + add("pdbx_average_fsc_work", number_or_qmark(ref.fsc_work)); + if (st.meta.has(&RefinementInfo::fsc_free)) + add("pdbx_average_fsc_free", number_or_qmark(ref.fsc_free)); + if (st.meta.has(&RefinementInfo::cc_intensity_work)) + add("correlation_coeff_I_to_Fcsqd_work", number_or_qmark(ref.cc_intensity_work)); + if (st.meta.has(&RefinementInfo::cc_intensity_free)) + add("correlation_coeff_I_to_Fcsqd_free", number_or_qmark(ref.cc_intensity_free)); + if (!st.meta.solved_by.empty()) + add("pdbx_method_to_determine_struct", string_or_qmark(st.meta.solved_by)); + if (!st.meta.starting_model.empty()) + add("pdbx_starting_model", string_or_qmark(st.meta.starting_model)); + if (!std::isnan(ref.luzzati_error)) + analyze_loop.add_row({id, + cif::quote(ref.id), + number_or_qmark(ref.luzzati_error)}); + for (const RefinementInfo::Restr& restr : ref.restr_stats) + restr_loop.add_row({cif::quote(ref.id), + cif::quote(restr.name), + int_or_qmark(restr.count), + number_or_qmark(restr.weight), + string_or_qmark(restr.function), + number_or_qmark(restr.dev_ideal)}); + for (const BasicRefinementInfo& bin : ref.bins) { + shell_loop.add_values({cif::quote(ref.id), + number_or_dot(bin.resolution_high), + number_or_qmark(bin.resolution_low), + number_or_qmark(bin.completeness), + int_or_qmark(get_number_obs(bin)), + int_or_qmark(get_number_work(bin)), + int_or_qmark(bin.rfree_set_count), + number_or_qmark(bin.r_all), + number_or_qmark(bin.r_work), + number_or_qmark(bin.r_free)}); + if (has_shell_fsc) + shell_loop.add_values({number_or_qmark(bin.fsc_work), + number_or_qmark(bin.fsc_free)}); + if (has_shell_ffcc) + shell_loop.add_values({number_or_qmark(bin.cc_fo_fc_work), + number_or_qmark(bin.cc_fo_fc_free)}); + if (has_shell_iicc) + shell_loop.add_values({number_or_qmark(bin.cc_intensity_work), + number_or_qmark(bin.cc_intensity_free)}); + } + } + assert(shell_loop.values.size() % shell_loop.tags.size() == 0); + assert(loop.values.size() % loop.tags.size() == 0); + } + + if (groups.title_keywords) { + auto title = st.info.find("_struct.title"); + if (title != st.info.end()) { + cif::ItemSpan span(block.items, "_struct."); + span.set_pair("_struct.entry_id", id); + span.set_pair(title->first, cif::quote(title->second)); + } + auto pdbx_keywords = st.info.find("_struct_keywords.pdbx_keywords"); + auto keywords = st.info.find("_struct_keywords.text"); + cif::ItemSpan span(block.items, "_struct_keywords."); + if (pdbx_keywords != st.info.end() || keywords != st.info.end()) + span.set_pair("_struct_keywords.entry_id", id); + if (pdbx_keywords != st.info.end()) + span.set_pair(pdbx_keywords->first, cif::quote(pdbx_keywords->second)); + if (keywords != st.info.end()) + span.set_pair(keywords->first, cif::quote(keywords->second)); + } + + if (groups.ncs) + write_ncs_oper(st, block); + + if (groups.struct_asym) { + cif::Loop& asym_loop = block.init_mmcif_loop("_struct_asym.", + {"id", "entity_id"}); + for (const Chain& chain : st.models[0].chains) + for (ConstResidueSpan& sub : chain.subchains()) { + const std::string& sub_id = sub.subchain_id(); + if (!sub_id.empty()) { + const Entity* ent = find_entity_of_subchain(sub_id, st.entities); + asym_loop.add_row({sub_id, (ent ? qchain(ent->name) : "?")}); + } + } + } + + bool nontrivial_origx = st.has_origx && !st.origx.is_identity(); + if (groups.origx && nontrivial_origx) { // _database_PDB_matrix (ORIGX) + cif::ItemSpan span(block.items, "_database_PDB_matrix."); + span.set_pair("_database_PDB_matrix.entry_id", id); + std::string tag_mat = "_database_PDB_matrix.origx[0][0]"; + std::string tag_vec = "_database_PDB_matrix.origx_vector[0]"; + for (int i = 0; i < 3; ++i) { + tag_mat[27] += 1; // origx[0] -> origx[1] -> origx[2] + tag_vec[34] += 1; + for (int j = 0; j < 3; ++j) { + tag_mat[30] = '1' + j; + span.set_pair(tag_mat, to_str(st.origx.mat[i][j])); + } + span.set_pair(tag_vec, to_str(st.origx.vec.at(i))); + } + } + + if (groups.struct_conf && !st.helices.empty()) { + cif::Loop& struct_conf_loop = block.init_mmcif_loop("_struct_conf.", + {"conf_type_id", "id", + "beg_auth_asym_id", "beg_label_asym_id", "beg_label_comp_id", + "beg_label_seq_id", "beg_auth_seq_id", "pdbx_beg_PDB_ins_code", + "end_auth_asym_id", "end_label_asym_id", "end_label_comp_id", + "end_label_seq_id", "end_auth_seq_id", "pdbx_end_PDB_ins_code", + "pdbx_PDB_helix_class", "pdbx_PDB_helix_length"}); + int count = 0; + for (const Helix& helix : st.helices) { + const_CRA cra1 = st.models[0].find_cra(helix.start); + const_CRA cra2 = st.models[0].find_cra(helix.end); + if (!cra1.residue || !cra2.residue) + continue; + struct_conf_loop.add_row({ + "HELX_P", // conf_type_id + "H" + std::to_string(++count), // id + qchain(cra1.chain->name), // beg_auth_asym_id + subchain_or_dot(*cra1.residue), // beg_label_asym_id + cra1.residue->name, // beg_label_comp_id + cra1.residue->label_seq.str(), // beg_label_seq_id + cra1.residue->seqid.num.str(), // beg_auth_seq_id + pdbx_icode(*cra1.residue), // beg_PDB_ins_code + qchain(cra2.chain->name), // end_auth_asym_id + subchain_or_dot(*cra2.residue), // end_label_asym_id + cra2.residue->name, // end_label_comp_id + cra2.residue->label_seq.str(), // end_label_seq_id + cra2.residue->seqid.num.str(), // end_auth_seq_id + pdbx_icode(*cra2.residue), // end_PDB_ins_code + std::to_string((int)helix.pdb_helix_class), // pdbx_PDB_helix_class + int_or_qmark(helix.length) // pdbx_PDB_helix_length + }); + } + if (count != 0) + block.set_pair("_struct_conf_type.id", "HELX_P"); + } + + // _struct_sheet* + if (groups.struct_sheet && !st.sheets.empty()) { + cif::Loop& sheet_loop = block.init_mmcif_loop("_struct_sheet.", + {"id", "number_strands"}); + for (const Sheet& sheet : st.sheets) + sheet_loop.add_row({string_or_dot(sheet.name), + std::to_string(sheet.strands.size())}); + + cif::Loop& order_loop = block.init_mmcif_loop("_struct_sheet_order.", + {"sheet_id", "range_id_1", "range_id_2", "sense"}); + for (const Sheet& sheet : st.sheets) + for (size_t i = 1; i < sheet.strands.size(); ++i) { + const Sheet::Strand& strand = sheet.strands[i]; + if (strand.sense != 0) + order_loop.add_row({string_or_dot(sheet.name), + std::to_string(i), std::to_string(i+1), + strand.sense > 0 ? "parallel" : "anti-parallel"}); + } + + cif::Loop& range_loop = block.init_mmcif_loop("_struct_sheet_range.", + {"sheet_id", "id", + "beg_auth_asym_id", "beg_label_asym_id", "beg_label_comp_id", + "beg_label_seq_id", "beg_auth_seq_id", "pdbx_beg_PDB_ins_code", + "end_auth_asym_id", "end_label_asym_id", "end_label_comp_id", + "end_label_seq_id", "end_auth_seq_id", "pdbx_end_PDB_ins_code"}); + for (const Sheet& sheet : st.sheets) + for (size_t i = 0; i < sheet.strands.size(); ++i) { + const Sheet::Strand& strand = sheet.strands[i]; + const_CRA cra1 = st.models[0].find_cra(strand.start); + const_CRA cra2 = st.models[0].find_cra(strand.end); + if (!cra1.residue || !cra2.residue) + continue; + range_loop.add_row({ + string_or_dot(sheet.name), // sheet_id + std::to_string(i+1), // id + qchain(cra1.chain->name), // beg_auth_asym_id + subchain_or_dot(*cra1.residue), // beg_label_asym_id + cra1.residue->name, // beg_label_comp_id + cra1.residue->label_seq.str(), // beg_label_seq_id + cra1.residue->seqid.num.str(), // beg_auth_seq_id + pdbx_icode(*cra1.residue), // beg_PDB_ins_code + qchain(cra2.chain->name), // end_auth_asym_id + subchain_or_dot(*cra2.residue), // end_label_asym_id + cra2.residue->name, // end_label_comp_id + cra2.residue->label_seq.str(), // end_label_seq_id + cra2.residue->seqid.num.str(), // end_auth_seq_id + pdbx_icode(*cra2.residue) // end_PDB_ins_code + }); + } + + cif::Loop& hbond_loop = block.init_mmcif_loop("_pdbx_struct_sheet_hbond.", + {"sheet_id", "range_id_1", "range_id_2", + "range_1_auth_asym_id", "range_1_label_asym_id", + "range_1_label_comp_id", "range_1_label_seq_id", "range_1_auth_seq_id", + "range_1_PDB_ins_code", "range_1_label_atom_id", + "range_2_auth_asym_id", "range_2_label_asym_id", + "range_2_label_comp_id", "range_2_label_seq_id", "range_2_auth_seq_id", + "range_2_PDB_ins_code", "range_2_label_atom_id"}); + for (const Sheet& sheet : st.sheets) + for (size_t i = 1; i < sheet.strands.size(); ++i) { + const Sheet::Strand& strand = sheet.strands[i]; + if (strand.hbond_atom2.atom_name.empty()) + continue; + // hbond_atomN is not a full atom "address": altloc is missing + const_CRA cra1 = st.models[0].find_cra(strand.hbond_atom1); + const_CRA cra2 = st.models[0].find_cra(strand.hbond_atom2); + if (!cra1.residue || !cra2.residue) + continue; + hbond_loop.add_row({ + string_or_dot(sheet.name), // sheet_id + std::to_string(i), // range_id_1 + std::to_string(i+1), // range_id_2 + qchain(cra1.chain->name), // range_1_auth_asym_id + subchain_or_dot(*cra1.residue), // range_1_label_asym_id + cra1.residue->name, // range_1_label_comp_id + cra1.residue->label_seq.str(), // range_1_label_seq_id + cra1.residue->seqid.num.str(), // range_1_auth_seq_id + pdbx_icode(*cra1.residue), // range_1_PDB_ins_code + cif::quote(strand.hbond_atom1.atom_name), // range_1_label_atom_id + qchain(cra2.chain->name), // range_2_auth_asym_id + subchain_or_dot(*cra2.residue), // range_2_label_asym_id + cra2.residue->name, // range_2_label_comp_id + cra2.residue->label_seq.str(), // range_2_label_seq_id + cra2.residue->seqid.num.str(), // range_2_auth_seq_id + pdbx_icode(*cra2.residue), // range_2_PDB_ins_code + cif::quote(strand.hbond_atom2.atom_name) // range_2_label_atom_id + }); + } + } + + // _pdbx_struct_assembly* and _struct_biol are REMARK 300/350 in PDB + if (groups.struct_biol && !st.meta.remark_300_detail.empty()) { + cif::ItemSpan span(block.items, "_struct_biol."); + span.set_pair("_struct_biol.id", "1"); + span.set_pair("_struct_biol.details", cif::quote(st.meta.remark_300_detail)); + } + if (groups.assembly && !st.assemblies.empty()) + write_assemblies(st, block); + + if (groups.conn) + write_struct_conn(st, block); + + if (groups.cis) // _struct_mon_prot_cis + write_cispeps(st, block); + + // _pdbx_struct_mod_residue (MODRES) + if (groups.modres && !st.mod_residues.empty()) { + bool use_ccp4_mod_id = false; + for (const ModRes& modres : st.mod_residues) + if (!modres.mod_id.empty()) + use_ccp4_mod_id = true; + cif::Loop& loop = block.init_mmcif_loop("_pdbx_struct_mod_residue.", + {"id", "auth_asym_id", "auth_seq_id", "PDB_ins_code", "auth_comp_id", + "label_comp_id", "parent_comp_id", "details"}); + if (use_ccp4_mod_id) + loop.tags.push_back("_pdbx_struct_mod_residue.ccp4_mod_id"); + int counter = 0; + for (const ModRes& modres : st.mod_residues) { + loop.add_values({std::to_string(++counter), + qchain(modres.chain_name), + modres.res_id.seqid.num.str(), + pdbx_icode(modres.res_id), + string_or_dot(modres.res_id.name), + string_or_qmark(modres.res_id.name), + string_or_qmark(modres.parent_comp_id), + string_or_qmark(modres.details)}); + if (use_ccp4_mod_id) + loop.values.push_back(string_or_qmark(modres.mod_id)); + } + } + + // _atom_sites (SCALE) + if (groups.scale && (nontrivial_origx || st.cell.explicit_matrices)) { + cif::ItemSpan span(block.items, "_atom_sites."); + span.set_pair("_atom_sites.entry_id", id); + std::string prefix = "_atom_sites.fract_transf_"; + for (int i = 0; i < 3; ++i) { + std::string idx = "[" + std::to_string(i + 1) + "]"; + const auto& frac = st.cell.frac; + std::string matrix_idx = prefix + "matrix"; + matrix_idx += idx; + span.set_pair(matrix_idx + "[1]", to_str(frac.mat[i][0])); + span.set_pair(matrix_idx + "[2]", to_str(frac.mat[i][1])); + span.set_pair(matrix_idx + "[3]", to_str(frac.mat[i][2])); + span.set_pair(cat(prefix, "vector", idx), to_str(frac.vec.at(i))); + } + } + + // _atom_type + if (groups.atom_type) { + std::array types{}; + for (const Model& model : st.models) + for (const Chain& chain : model.chains) + for (const Residue& res : chain.residues) + for (const Atom& atom : res.atoms) + types[atom.element.ordinal()] = true; + cif::Loop& atom_type_loop = block.init_mmcif_loop("_atom_type.", {"symbol"}); + for (int i = 0; i < (int)El::END; ++i) + if (types[i]) + atom_type_loop.add_row({Element((El)i).uname()}); + } + + if (groups.entity_poly_seq) { + cif::Loop& poly_loop = block.init_mmcif_loop("_entity_poly_seq.", + {"entity_id", "num", "mon_id", "hetero"}); + for (const Entity& ent : st.entities) + if (ent.entity_type == EntityType::Polymer) { + // SEQRES from PDB doesn't record microheterogeneity. + std::string hetero_no = ent.reflects_microhetero ? "n" : "?"; + for (size_t i = 0; i != ent.full_sequence.size(); ++i) { + const std::string& mon_ids = ent.full_sequence[i]; + std::string num = std::to_string(i+1); + size_t start = 0, end; + while ((end = mon_ids.find(',', start)) != std::string::npos) { + poly_loop.add_row({qchain(ent.name), num, + mon_ids.substr(start, end-start), "y"}); + start = end + 1; + } + poly_loop.add_row({qchain(ent.name), num, mon_ids.substr(start), + start == 0 ? hetero_no : "y"}); + } + } + } + + if (groups.atoms) + add_cif_atoms(st, block, groups.group_pdb, groups.auth_all); + + if (groups.tls && st.meta.get_tls_groups() != nullptr) { + // pdbx_refine_id doesn't make sense here, but it's required + // by the mmCIF spec. In joint refinement, TLS constraints can't be + // specific to a dataset, because they constrain the shared model. + cif::Loop& loop = block.init_mmcif_loop("_pdbx_refine_tls.", { + "id", "pdbx_refine_id", + "origin_x", "origin_y", "origin_z", + "T[1][1]", "T[2][2]", "T[3][3]", "T[1][2]", "T[1][3]", "T[2][3]", + "L[1][1]", "L[2][2]", "L[3][3]", "L[1][2]", "L[1][3]", "L[2][3]", + "S[1][1]", "S[1][2]", "S[1][3]", + "S[2][1]", "S[2][2]", "S[2][3]", + "S[3][1]", "S[3][2]", "S[3][3]"}); + for (const RefinementInfo& ref : st.meta.refinement) + for (const TlsGroup& tls : ref.tls_groups) { + const SMat33& T = tls.T; + const SMat33& L = tls.L; + const Mat33& S = tls.S; + auto q = number_or_qmark; + loop.add_row({string_or_dot(tls.id), cif::quote(ref.id), + q(tls.origin.x), q(tls.origin.y), q(tls.origin.z), + q(T.u11), q(T.u22), q(T.u33), q(T.u12), q(T.u13), q(T.u23), + q(L.u11), q(L.u22), q(L.u33), q(L.u12), q(L.u13), q(L.u23), + q(S[0][0]), q(S[0][1]), q(S[0][2]), + q(S[1][0]), q(S[1][1]), q(S[1][2]), + q(S[2][0]), q(S[2][1]), q(S[2][2])}); + } + cif::Loop& group_loop = block.init_mmcif_loop("_pdbx_refine_tls_group.", { + "id", "refine_tls_id", "pdbx_refine_id", + "beg_auth_asym_id", "beg_auth_seq_id", "beg_PDB_ins_code", + "end_auth_asym_id", "end_auth_seq_id", "end_PDB_ins_code", + "selection_details"}); + int counter = 1; + for (const RefinementInfo& ref : st.meta.refinement) + for (const TlsGroup& tls : ref.tls_groups) + for (const TlsGroup::Selection& sel : tls.selections) + group_loop.add_row({std::to_string(counter++), + string_or_dot(tls.id), + cif::quote(ref.id), + string_or_qmark(sel.chain), + sel.res_begin.num.str(), + pdbx_icode(sel.res_begin), + string_or_qmark(sel.chain), + sel.res_end.num.str(), + pdbx_icode(sel.res_end), + string_or_qmark(sel.details)}); + } + + if (groups.software && !st.meta.software.empty()) { + bool write_all_fields = false; + for (const SoftwareItem& item : st.meta.software) + if (!item.date.empty() || !item.description.empty() || + !item.contact_author.empty() || !item.contact_author_email.empty()) + write_all_fields = true; + cif::Loop& loop = block.init_mmcif_loop("_software.", + {"pdbx_ordinal", "classification", "name", "version"}); + if (write_all_fields) + loop.tags.insert(loop.tags.end(), + {"_software.date", "_software.description", + "_software.contact_author", "_software.contact_author_email"}); + int ordinal = 0; + for (const SoftwareItem& item : st.meta.software) { + loop.add_values({ + std::to_string(++ordinal), + cif::quote(software_classification_to_string(item.classification)), + cif::quote(item.name), + string_or_dot(item.version)}); + if (write_all_fields) + loop.add_values({ + string_or_qmark(item.date), + string_or_qmark(item.description), + string_or_qmark(item.contact_author), + string_or_qmark(item.contact_author_email)}); + } + } +} + +cif::Document make_mmcif_document(const Structure& st, MmcifOutputGroups groups) { + cif::Document doc; + doc.blocks.resize(1); + update_mmcif_block(st, doc.blocks[0], groups); + return doc; +} + +cif::Block make_mmcif_block(const Structure& st, MmcifOutputGroups groups) { + cif::Block block; + update_mmcif_block(st, block, groups); + return block; +} + +cif::Block make_mmcif_headers(const Structure& st) { + MmcifOutputGroups groups(true); + groups.atoms = false; + return make_mmcif_block(st, groups); +} + +void add_minimal_mmcif_data(const Structure& st, cif::Block& block) { + cif::ItemSpan cell_span(block.items, "_cell."); + write_cell_parameters(st.cell, cell_span); + block.set_pair("_symmetry.space_group_name_H-M", cif::quote(st.spacegroup_hm)); + write_ncs_oper(st, block); + add_cif_atoms(st, block, /*use_group_pdb=*/false, /*auth_all=*/false); +} + +} // namespace gemmi diff --git a/gemmi_gph/to_pdb.cpp b/gemmi_gph/to_pdb.cpp new file mode 100644 index 000000000..04d411511 --- /dev/null +++ b/gemmi_gph/to_pdb.cpp @@ -0,0 +1,722 @@ +// Copyright 2017-2023 Global Phasing Ltd. + +#include + +#include +#include // for isdigit +#include // for memset, memcpy +#include +#include +#include // for ostringstream + +#include // for fail +#include +#include // for find_tabulated_residue +#include + +namespace gemmi { + +#define WRITE(...) do { \ + snprintf_z(buf, 82, __VA_ARGS__); \ + buf[80] = '\n'; \ + os.write(buf, 81); \ + } while(0) + +#define WRITEU(...) do { \ + snprintf_z(buf, 82, __VA_ARGS__); \ + buf[80] = '\n'; \ + for (int i_ = 0; i_ != 80; i_++) \ + if (buf[i_] >= 'a' && buf[i_] <= 'z') buf[i_] -= 0x20; \ + os.write(buf, 81); \ + } while(0) + +#define WRITELN(...) do { \ + int length__ = snprintf_z(buf, 82, __VA_ARGS__); \ + if (length__ < 80) \ + std::memset(buf + length__, ' ', 80 - length__); \ + buf[80] = '\n'; \ + os.write(buf, 81); \ + } while(0) + +bool use_hetatm(const Residue& res) { + if (res.het_flag == 'H') + return true; + if (res.het_flag == 'A') + return false; + if (res.entity_type == EntityType::Branched || + res.entity_type == EntityType::NonPolymer || + res.entity_type == EntityType::Water) + return true; + return !find_tabulated_residue(res.name).is_standard(); +} + +namespace { + +// works for non-negative values only +void base36_encode(char* buffer, int width, int value) { + const char base36[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + buffer[width] = '\0'; + do { + buffer[--width] = base36[value % 36]; + value /= 36; + } while (value != 0 && width != 0); + while (width != 0) + buffer[--width] = ' '; +} + +// based on http://cci.lbl.gov/hybrid_36/ +std::array encode_serial_in_hybrid36(int serial) { + std::array str; + assert(serial >= 0); + if (serial < 100000) + to_chars_z(str.data(), str.data() + 8, serial); + else + base36_encode(str.data(), 5, serial + (10 * 36 * 36 * 36 * 36 - 100000)); + return str; +} + +std::array write_seq_id(const SeqId& seqid) { + std::array str; + char* ptr = str.data(); + if (*seqid.num > -1000 && *seqid.num < 10000) { + ptr = to_chars_z(ptr, ptr + 5, *seqid.num); + } else if (seqid.num.has_value()) { + base36_encode(ptr, 4, *seqid.num + (10 * 36 * 36 * 36 - 10000)); + ptr += 4; + } + *ptr++ = seqid.icode; + *ptr = '\0'; + return str; +} + +const char* find_last_break(const char *str, int max_len) { + int last_break = 0; + for (int i = 0; i < max_len; i++) { + if (str[i] == '\0') + return str + i; + if (str[i] == ' ' || str[i] == '-') + last_break = i + 1; + } + return str + (last_break != 0 ? last_break : max_len); +} + +// Write record with possible continuation lines, with the format: +// 1-6 record name, 8-10 continuation, 11-lastcol string. +void write_multiline(std::ostream& os, const char* record_name, + const std::string& text, int lastcol) { + if (text.empty()) + return; + char buf[88]; // a few bytes extra, just in case + const char *start = text.c_str(); + const char *end = find_last_break(start, lastcol-10); + WRITEU("%-6s %-70.*s", record_name, static_cast(end-start), start); + for (int n = 2; n < 1000 && *end != '\0'; ++n) { + start = end; + end = find_last_break(start, lastcol-11); + int len = int(end - start); + WRITEU("%-6s %3d %-69.*s", record_name, n, len, start); + } +} + +void write_ncs_op(const NcsOp& op, std::ostream& os) { + char buf[88]; + for (int i = 0; i < 3; ++i) { + WRITE("MTRIX%d %3.3s%10.6f%10.6f%10.6f %14.5f %-21c", i+1, + op.id.c_str(), op.tr.mat[i][0], op.tr.mat[i][1], op.tr.mat[i][2], + op.tr.vec.at(i), op.given ? '1' : ' '); + } +} + +void write_remarks(const Structure& st, std::ostream& os) { + char buf[88]; + if (st.resolution > 0) { + WRITE("%-80s", "REMARK 2"); + WRITE("REMARK 2 RESOLUTION. %7.2f %-49s", st.resolution, "ANGSTROMS."); + } + if (!st.assemblies.empty()) { + const char* preface[] = { + "REMARK 350", + "REMARK 350 COORDINATES FOR A COMPLETE MULTIMER REPRESENTING THE KNOWN", + "REMARK 350 BIOLOGICALLY SIGNIFICANT OLIGOMERIZATION STATE OF THE", + "REMARK 350 MOLECULE CAN BE GENERATED BY APPLYING BIOMT TRANSFORMATIONS", + "REMARK 350 GIVEN BELOW. BOTH NON-CRYSTALLOGRAPHIC AND", + "REMARK 350 CRYSTALLOGRAPHIC OPERATIONS ARE GIVEN." + }; + for (const char* line : preface) + WRITE("%-80s", line); + int counter = 0; + for (const Assembly& assem : st.assemblies) { + WRITE("%-80s", "REMARK 350"); + WRITE("REMARK 350 BIOMOLECULE: %-56d", ++counter); + if (assem.author_determined) + WRITEU("REMARK 350 AUTHOR DETERMINED BIOLOGICAL UNIT: %-34s", + assem.oligomeric_details.c_str()); + if (assem.software_determined) { + WRITEU("REMARK 350 SOFTWARE DETERMINED QUATERNARY STRUCTURE: %-27s", + assem.oligomeric_details.c_str()); + if (!assem.software_name.empty()) + WRITEU("REMARK 350 SOFTWARE USED: %-54s", + assem.software_name.c_str()); + if (!std::isnan(assem.absa)) + WRITELN("REMARK 350 TOTAL BURIED SURFACE AREA: %.0f ANGSTROM**2", + assem.absa); + if (!std::isnan(assem.ssa)) + WRITELN("REMARK 350 SURFACE AREA OF THE COMPLEX: %.0f ANGSTROM**2", + assem.ssa); + if (!std::isnan(assem.more)) + WRITELN("REMARK 350 CHANGE IN SOLVENT FREE ENERGY: %.1f KCAL/MOL", + assem.more); + } + int oper_cnt = 0; + for (const Assembly::Gen& gen : assem.generators) { + std::string chains_str; + if (!gen.chains.empty()) { + chains_str = join_str(gen.chains, ", "); + } else { + // subchains -> chains + std::vector chains; + for (const Chain& ch : st.first_model().chains) + if (!ch.residues.empty() && !in_vector(ch.name, chains) && + in_vector(ch.residues[0].subchain, gen.subchains)) + chains.push_back(ch.name); + chains_str = join_str(chains, ", "); + } + size_t end = chains_str.length(); + if (end >= 30) + end = chains_str.rfind(' ', 29); + WRITELN("REMARK 350 APPLY THE FOLLOWING TO CHAINS: %s", + chains_str.substr(0, end).c_str()); + while (end < chains_str.length()) { + size_t begin = end + 1; + end = chains_str.length(); + if (end - begin >= 30) + end = chains_str.rfind(' ', begin + 29); + WRITELN("REMARK 350 AND CHAINS: %s", + chains_str.substr(begin, end - begin).c_str()); + } + for (const Assembly::Operator& oper : gen.operators) { + ++oper_cnt; + const Transform& tr = oper.transform; + for (int i = 0; i < 3; ++i) + WRITE("REMARK 350 " + "BIOMT%d %3d%10.6f%10.6f%10.6f %14.5f ", + i+1, oper_cnt, + tr.mat[i][0], tr.mat[i][1], tr.mat[i][2], tr.vec.at(i)); + } + } + } + } +} + +void write_chain_atoms(const Chain& chain, std::ostream& os, + int& serial, PdbWriteOptions opt) { + char buf[88]; + buf[0] = '\0'; + if (chain.name.length() > 2) + fail("long chain name: " + chain.name); + for (const Residue& res : chain.residues) { + bool as_het = use_hetatm(res); + for (const Atom& a : res.atoms) { + serial = opt.preserve_serial ? a.serial : serial + 1; + // 1- 6 6s record name + // 7-11 5d integer serial + // 12 1 - + // 13-16 4s atom name (from 13 only if 4-char or 2-char symbol) + // 17 1c altloc + // 18-20 3s residue name + // 21 1 - + // 22 1s chain + // 23-26 4d integer residue sequence number + // 27 1c insertion code + // 28-30 3 - + // 31-38 8f x (8.3) + // 39-46 8f y + // 47-54 8f z + // 55-60 6f occupancy (6.2) + // 61-66 6f temperature factor (6.2) + // 67-76 6 - + // 73-76 segment identifier, left-justified (non-standard) + // 77-78 2s element symbol, right-justified + // 79-80 2s charge + int written_bytes = snprintf_z(buf, 82, + "%-6s%5s %-4.4s%c%3.3s" + "%2s%5s %8.3f%8.3f%8.3f", + as_het ? "HETATM" : "ATOM", + encode_serial_in_hybrid36(serial).data(), + a.padded_name().c_str(), + a.altloc ? std::toupper(a.altloc) : ' ', + res.name.c_str(), + chain.name.c_str(), + write_seq_id(res.seqid).data(), + // We want to avoid negative zero and round the numbers up + // if they originally had one digit more and that digit was 5. + a.pos.x > -5e-4 && a.pos.x < 0 ? 0 : a.pos.x + 1e-10, + a.pos.y > -5e-4 && a.pos.y < 0 ? 0 : a.pos.y + 1e-10, + a.pos.z > -5e-4 && a.pos.z < 0 ? 0 : a.pos.z + 1e-10); + if GEMMI_UNLIKELY(written_bytes > 54) { + // The only items expected to overflow above are the coordinates, + // if the integer part of the number exceeds 5 characters. + // This happens when something goes wrong and the model is far from + // the origin. Such a model should be shifted; it can't be written it + // in a spec-conforming format: Real(8.3). Here we overwrite the last + // digits - trimming is better than overflowing the line. + snprintf_z(buf+38, 82-38, "%8.3f", a.pos.y); + snprintf_z(buf+46, 82-46, "%8.3f", a.pos.z); + } + snprintf_z(buf+54, 82-54, + "%6.2f%6.2f %-4.4s%2s%c%c", + // Occupancy is stored as single prec, but we know it's <= 1, + // so no precision is lost even if it had 6 digits after dot. + a.occ + 1e-6, + // B is harder to get rounded right. It is stored as float, + // and may be given with more than single precision in mmCIF + // If it was originally %.5f (5TIS) we need to add 0.5 * 10^-5. + std::min(a.b_iso + 0.5e-5, 999.99), + res.segment.c_str(), + a.element.uname(), + // Charge is written as 1+ or 2-, etc, or just empty space. + // Sometimes PDB files have explicit 0s (5M05); we ignore them. + a.charge ? a.charge > 0 ? '0'+a.charge : '0'-a.charge : ' ', + a.charge ? a.charge > 0 ? '+' : '-' : ' '); + buf[80] = '\n'; + os.write(buf, 81); + if (a.aniso.nonzero()) { + // re-using part of the buffer + std::memcpy(buf, "ANISOU", 6); + const double eps = 1e-6; + snprintf_z(buf+28, 43, "%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f", + a.aniso.u11*1e4 + eps, a.aniso.u22*1e4 + eps, + a.aniso.u33*1e4 + eps, a.aniso.u12*1e4 + eps, + a.aniso.u13*1e4 + eps, a.aniso.u23*1e4 + eps); + buf[28+42] = ' '; + buf[80] = '\n'; + os.write(buf, 81); + } + } + if (opt.ter_records && buf[0] != '\0' && + (opt.ter_ignores_type ? &res == &chain.residues.back() + : (res.entity_type == EntityType::Polymer && + (&res == &chain.residues.back() || + (&res + 1)->entity_type != EntityType::Polymer)))) { + if (opt.numbered_ter) { + // re-using part of the buffer in the middle, e.g.: + // TER 4153 LYS B 286 + snprintf_z(buf, 82, "TER %5s", + encode_serial_in_hybrid36(++serial).data()); + std::memset(buf+11, ' ', 6); + std::memset(buf+28, ' ', 52); + buf[80] = '\n'; + os.write(buf, 81); + } else { + WRITE("%-80s", "TER"); + } + } + } +} + +} // anonymous namespace + +void write_pdb(const Structure& st, std::ostream& os, PdbWriteOptions opt) { + // check if structure can be written as pdb + for (const gemmi::Model& model : st.models) + for (const gemmi::Chain& chain : model.chains) + if (chain.name.size() > 2) + gemmi::fail("chain name too long for the PDB format: " + chain.name); + + const std::string& entry_id = st.get_info("_entry.id"); + const char* entry_id_4 = entry_id.size() <= 4 ? entry_id.c_str() : ""; + char buf[88]; + + // HEADER + if (!opt.minimal_file) { + const char* months = "JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC???"; + const std::string& date = + st.get_info("_pdbx_database_status.recvd_initial_deposition_date"); + std::string pdb_date; + if (date.size() == 10) { + unsigned month_idx = 10 * (date[5] - '0') + date[6] - '0' - 1; + std::string month(months + 3 * std::min(month_idx, 13u), 3); + pdb_date = date.substr(8, 2) + "-" + month + "-" + date.substr(2, 2); + } + // "classification" in PDB == _struct_keywords.pdbx_keywords in mmCIF + const std::string& keywords = st.get_info("_struct_keywords.pdbx_keywords"); + if (!pdb_date.empty() || !keywords.empty() || !entry_id.empty()) + WRITEU("HEADER %-40.40s%-9s %-18.18s", + keywords.c_str(), pdb_date.c_str(), entry_id.c_str()); + } + + // TITLE, KEYWDS, EXPDTA, NUMMDL, REMARK + if (!opt.minimal_file) { + write_multiline(os, "TITLE", st.get_info("_struct.title"), 80); + write_multiline(os, "KEYWDS", st.get_info("_struct_keywords.text"), 79); + std::string expdta = st.get_info("_exptl.method"); + if (expdta.empty()) + expdta = join_str(st.meta.experiments, "; ", + [](const ExperimentInfo& e) { return e.method; }); + write_multiline(os, "EXPDTA", expdta, 79); + if (st.models.size() > 1) + WRITE("NUMMDL %-6zu %63s", st.models.size(), ""); + + if (!st.raw_remarks.empty()) { + for (const std::string& line : st.raw_remarks) { + os << line; + if (line.empty() || line.back() != '\n') + os << '\n'; + } + } else { + write_remarks(st, os); + } + } + + // DBREF[12], SEQRES + if (opt.seqres_records && !st.models.empty()) { + std::vector entity_list; + entity_list.reserve(st.models[0].chains.size()); + for (const Chain& ch : st.models[0].chains) { + const Entity* entity = st.get_entity_of(ch.get_polymer()); + // If the input pdb file has no TER records the subchains and entities + // are not setup automatically. In such case it is possible to call + // setup_entities() to use heuristic to split polymer and ligands + // and assign entities. But if it was not called, we may still find + // the original SEQRES in the entity named after the chain name. + if (!entity && st.input_format == CoorFormat::Pdb && + !ch.residues.empty() && ch.residues[0].subchain.empty()) { + entity = st.get_entity(ch.name); + if (entity && !entity->subchains.empty()) + entity = nullptr; + } + entity_list.push_back(entity); + } + // DBREF / DBREF1 / DBREF2 + for (size_t i = 0; i != entity_list.size(); ++i) + if (const Entity* entity = entity_list[i]) { + const Chain& ch = st.models[0].chains[i]; + for (const Entity::DbRef& dbref : entity->dbrefs) { + bool short_record = *dbref.db_end.num < 100000 && + dbref.accession_code.size() < 9 && + dbref.id_code.size() < 13; + SeqId begin = dbref.seq_begin; + SeqId end = dbref.seq_end; + if (!begin.num || !end.num) + if (ConstResidueGroup polymer = ch.get_polymer()) { + begin = polymer.label_seq_id_to_auth(dbref.label_seq_begin); + end = polymer.label_seq_id_to_auth(dbref.label_seq_end); + } + snprintf_z(buf, 82, "DBREF %4s%2s %5s %5s %-6s ", + entry_id_4, ch.name.c_str(), + write_seq_id(begin).data(), + write_seq_id(end).data(), + dbref.db_name.c_str()); + if (dbref.db_name == "PDB" && dbref.id_code == entry_id) { + // PDB uses self-reference for fragments that don't have real + // reference. No idea why. In such case the same begin/end is used. + } else { + begin = dbref.db_begin; + end = dbref.db_end; + } + if (short_record) { + snprintf_z(buf+33, 82-33, "%-8s %-12s %5d%c %5d%c \n", + dbref.accession_code.c_str(), dbref.id_code.c_str(), + *begin.num, begin.icode, *end.num, end.icode); + } else { + buf[5] = '1'; // -> DBREF1 + snprintf_z(buf+33, 82-33, " %-33s\n", + dbref.id_code.c_str()); + } + buf[80] = '\n'; + os.write(buf, 81); + if (!short_record) + WRITE("DBREF2 %4s%2s %-22s %10d %10d ", + entry_id_4, ch.name.c_str(), + dbref.accession_code.c_str(), *begin.num, *end.num); + } + } + // SEQRES + for (size_t i = 0; i != entity_list.size(); ++i) + if (const Entity* entity = entity_list[i]) { + const Chain& ch = st.models[0].chains[i]; + int row = 0; + int col = 0; + for (const std::string& monomers : entity->full_sequence) { + if (col == 0) + snprintf_z(buf, 82, "SEQRES%4d%2s%5zu %62s\n", + ++row, ch.name.c_str(), + entity->full_sequence.size(), ""); + size_t end = std::min(monomers.find(','), std::min(monomers.length(), (size_t)3)); + std::memcpy(buf + 18 + 4*col + 4-end, monomers.c_str(), end); + if (++col == 13) { + os.write(buf, 81); + col = 0; + } + } + if (col != 0) { + buf[80] = '\n'; + os.write(buf, 81); + } + } + } + + if (!opt.minimal_file) { + // MODRES + for (const ModRes& modres : st.mod_residues) { + snprintf_z(buf, 82, "MODRES %4s %3.3s%2s %5s %3s %-41.41s %-8.8s\n", + entry_id_4, + modres.res_id.name.c_str(), + modres.chain_name.c_str(), + write_seq_id(modres.res_id.seqid).data(), + modres.parent_comp_id.c_str(), + modres.details.c_str(), + modres.mod_id.c_str()); + buf[80] = '\n'; + for (int i_ = 29; i_ != 80; i_++) + if (buf[i_] >= 'a' && buf[i_] <= 'z') buf[i_] -= 0x20; + os.write(buf, 81); + } + + // HET + if (!st.models.empty()) { + auto full_ccd_code = [&](const std::string& name) -> std::string { + for (const auto& mapping : st.shortened_ccd_codes) + if (mapping.second == name) + return "real CCD code: " + mapping.first; + return ""; + }; + for (const Chain& chain : st.models[0].chains) + for (const Residue& res : chain.residues) + if (use_hetatm(res) && !find_tabulated_residue(res.name).is_water()) { + WRITE("HET %3s %2s%5s %6zu %-51s\n", + res.name.c_str(), + chain.name.c_str(), + write_seq_id(res.seqid).data(), + res.atoms.size(), + res.name[0] == '~' ? full_ccd_code(res.name).c_str() : ""); + } + } + } + + // HETNAM - but it's used only for tilde-hetnam extension + for (const auto& mapping : st.shortened_ccd_codes) + WRITE("HETNAM %3s %55s %-9s\n", mapping.second.c_str(), "", mapping.first.c_str()); + + // HELIX + if (!opt.minimal_file && !st.helices.empty()) { + int counter = 0; + for (const Helix& helix : st.helices) { + if (++counter == 10000) + counter = 0; + // According to the PDB spec serial number can be from 1 to 999. + // Here, we allow for up to 9999 helices by using columns 7 and 11. + snprintf_z(buf, 82, "HELIX %4d%4d %3.3s%2s %5s %3.3s%2s %5s%2d %35d \n", + counter, counter, + helix.start.res_id.name.c_str(), helix.start.chain_name.c_str(), + write_seq_id(helix.start.res_id.seqid).data(), + helix.end.res_id.name.c_str(), helix.end.chain_name.c_str(), + write_seq_id(helix.end.res_id.seqid).data(), + (int) helix.pdb_helix_class, helix.length); + if (helix.length < 0) // make 72-76 blank if the length is not given + std::memset(buf+71, ' ', 5); + buf[80] = '\n'; + os.write(buf, 81); + } + } + + // SHEET + if (!opt.minimal_file && !st.sheets.empty()) { + for (const Sheet& sheet : st.sheets) { + int strand_counter = 0; + for (const Sheet::Strand& strand : sheet.strands) { + const AtomAddress& a2 = strand.hbond_atom2; + const AtomAddress& a1 = strand.hbond_atom1; + // H-bond atom names are expected to be O and N + WRITE("SHEET%5d %3.3s%2zu %3.3s%2s%5s %3.3s%2s%5s%2d %-3s%3.3s%2s%5s " + " %-3s%3.3s%2s%5s ", + ++strand_counter, sheet.name.c_str(), sheet.strands.size(), + strand.start.res_id.name.c_str(), strand.start.chain_name.c_str(), + write_seq_id(strand.start.res_id.seqid).data(), + strand.end.res_id.name.c_str(), strand.end.chain_name.c_str(), + write_seq_id(strand.end.res_id.seqid).data(), strand.sense, + a2.atom_name.c_str(), a2.res_id.name.c_str(), + a2.chain_name.c_str(), + a2.res_id.seqid.num ? write_seq_id(a2.res_id.seqid).data() : "", + a1.atom_name.c_str(), a1.res_id.name.c_str(), + a1.chain_name.c_str(), + a1.res_id.seqid.num ? write_seq_id(a1.res_id.seqid).data() : ""); + } + } + } + + if (!st.models.empty()) { + // SSBOND (note: uses only the first model and primary conformation) + if (opt.ssbond_records) { + int counter = 0; + for (const Connection& con : st.connections) + if (con.type == Connection::Disulf) { + const_CRA cra1 = st.models[0].find_cra(con.partner1, true); + const_CRA cra2 = st.models[0].find_cra(con.partner2, true); + if (!cra1.atom || !cra2.atom) + continue; + NearestImage im = st.cell.find_nearest_image(cra1.atom->pos, + cra2.atom->pos, con.asu); + if (++counter == 10000) + counter = 0; + WRITE("SSBOND%4d %3.3s%2s %5s %3.3s%2s %5s %28s %6s %5.2f ", + counter, + cra1.residue->name.c_str(), cra1.chain->name.c_str(), + write_seq_id(cra1.residue->seqid).data(), + cra2.residue->name.c_str(), cra2.chain->name.c_str(), + write_seq_id(cra2.residue->seqid).data(), + "1555", im.symmetry_code(false).c_str(), im.dist()); + } + } + + // LINK (note: uses only the first model and primary conformation) + if (opt.link_records) { + for (const Connection& con : st.connections) + if (con.type == Connection::Covale || con.type == Connection::MetalC || + con.type == Connection::Unknown) { + const_CRA cra1 = st.models[0].find_cra(con.partner1, true); + const_CRA cra2 = st.models[0].find_cra(con.partner2, true); + // In special cases (LINKR gap) atoms are not there. + if (!cra1.residue || !cra2.residue) + continue; + std::string im_pdb_symbol, im_dist_str; + bool im_same_asu = true; + if (cra1.atom && cra2.atom) { + NearestImage im = st.cell.find_nearest_image(cra1.atom->pos, + cra2.atom->pos, con.asu); + im_pdb_symbol = im.symmetry_code(false); + im_dist_str = to_str_prec<2>(im.dist()); + im_same_asu = im.same_asu(); + } + // Pdb spec: "sym1 and sym2 are right justified and are given as + // blank when the identity operator (and no cell translation) is + // to be applied to the atom." But all files from wwPDB have + // 1555 not blank, so here we also write 1555, + // except for LINKR (Refmac variant of LINK). + snprintf_z(buf, 82, "LINK %-4s%c%3.3s%2s%5s " + " %-4s%c%3.3s%2s%5s %6s %6s %5s \n", + cra1.atom ? cra1.atom->padded_name().c_str() : "", + cra1.atom && cra1.atom->altloc ? std::toupper(cra1.atom->altloc) : ' ', + cra1.residue->name.c_str(), + con.partner1.chain_name.c_str(), + write_seq_id(cra1.residue->seqid).data(), + cra2.atom ? cra2.atom->padded_name().c_str() : "", + cra2.atom && cra2.atom->altloc ? std::toupper(cra2.atom->altloc) : ' ', + cra2.residue->name.c_str(), + con.partner2.chain_name.c_str(), + write_seq_id(cra2.residue->seqid).data(), + "1555", im_pdb_symbol.c_str(), im_dist_str.c_str()); + if (im_same_asu && !con.link_id.empty()) + std::memset(buf+58, ' ', 14); // erase symmetry + if ((opt.use_linkr ||opt.use_link_id) && !con.link_id.empty()) { + if (opt.use_linkr) + buf[4] = 'R'; // LINK -> LINKR + // overwrite distance with link_id + if (opt.use_linkr || opt.use_link_id) + snprintf_z(buf+72, 82-72, "%-8s\n", con.link_id.c_str()); + } + buf[80] = '\n'; + os.write(buf, 81); + } + } + + // CISPEP + if (opt.cispep_records) { + int counter = 0; + for (const CisPep& cispep : st.cispeps) { + WRITE("CISPEP%4d %3.3s%2s %5s %3.3s%2s %5s %9d %12.2f %20s", + ++counter, + cispep.partner_c.res_id.name.c_str(), + cispep.partner_c.chain_name.c_str(), + write_seq_id(cispep.partner_c.res_id.seqid).data(), + cispep.partner_n.res_id.name.c_str(), + cispep.partner_n.chain_name.c_str(), + write_seq_id(cispep.partner_n.res_id.seqid).data(), + st.models.size() > 1 ? cispep.model_num : 0, + std::isnan(cispep.reported_angle) ? 0. : cispep.reported_angle, + ""); + if (counter == 9999) + counter = 0; + } + } + } + + // CRYST1 + if (opt.cryst1_record) { + WRITE("CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4s ", + st.cell.a, st.cell.b, st.cell.c, st.cell.alpha, st.cell.beta, st.cell.gamma, + st.spacegroup_hm.empty() ? "P 1" : st.spacegroup_hm.c_str(), + st.get_info("_cell.Z_PDB").c_str()); + } + + // ORIGX1/2/3 + if (!opt.minimal_file && st.has_origx && !st.origx.is_identity()) { + for (int i = 0; i < 3; ++i) + WRITE("ORIGX%d %13.6f%10.6f%10.6f %14.5f %24s", i+1, + st.origx.mat[i][0], st.origx.mat[i][1], st.origx.mat[i][2], + st.origx.vec.at(i), ""); + } + + // SCALE1/2/3 + if (!opt.minimal_file && st.cell.explicit_matrices) { + for (int i = 0; i < 3; ++i) + // We add a small number to avoid negative 0. + WRITE("SCALE%d %13.6f%10.6f%10.6f %14.5f %24s", i+1, + st.cell.frac.mat[i][0] + 1e-15, st.cell.frac.mat[i][1] + 1e-15, + st.cell.frac.mat[i][2] + 1e-15, st.cell.frac.vec.at(i) + 1e-15, ""); + } + + // MTRIX1/2/3 + if (!st.ncs.empty()) { + auto identity = st.info.find("_struct_ncs_oper.id"); + if (identity != st.info.end() && + !in_vector_f([&](const NcsOp& op) { return op.id == identity->second; }, st.ncs)) + write_ncs_op(NcsOp{identity->second, true, {}}, os); + for (const NcsOp& op : st.ncs) + write_ncs_op(op, os); + } + + // MODEL, ATOM, HETATM, TER, ENDMDL + if (opt.atom_records) { + for (const Model& model : st.models) { + int serial = 0; + if (st.models.size() > 1) + WRITE("MODEL %8d %65s", model.num, ""); + for (const Chain& chain : model.chains) + write_chain_atoms(chain, os, serial, opt); + if (st.models.size() > 1) + WRITE("%-80s", "ENDMDL"); + } + } + + // CONECT + if (opt.conect_records) { + auto num_str = [](const std::vector& nums, size_t i) { + return i < nums.size() ? encode_serial_in_hybrid36(nums[i]) : std::array{}; + }; + for (const auto& num_pair : st.conect_map) + for (size_t i = 0; i < num_pair.second.size(); i += 4) + WRITE("CONECT%5s%5s%5s%5s%5s%50s", + encode_serial_in_hybrid36(num_pair.first).data(), + encode_serial_in_hybrid36(num_pair.second[i]).data(), + num_str(num_pair.second, i+1).data(), + num_str(num_pair.second, i+2).data(), + num_str(num_pair.second, i+3).data(), + ""); + } + + // END + if (opt.end_record) + WRITE("%-80s", "END"); +} + +std::string make_pdb_string(const Structure& st, PdbWriteOptions opt) { + std::ostringstream os; + write_pdb(st, os, opt); + return os.str(); +} + +#undef WRITE +#undef WRITEU + +} // namespace gemmi diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 458254fe9..a8d0e7c24 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -18,6 +18,8 @@ ADD_LIBRARY(Rugnux STATIC ResultReport.h SpotWidth.cpp SpotWidth.h + WriteModel.cpp + WriteModel.h ) TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi) diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 2cf79b608..0471551c1 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -765,6 +765,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.ok = true; result.maps_prefix = output_prefix; + // The placed coordinates, for the caller to write out beside the maps once the frame is settled. + result.placed_model = std::make_shared(st); logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) " "[overall + anisotropic B + bulk solvent]", result.r_work, result.n_work, result.r_free, result.n_free); diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index cd73f1c4a..0dcc61e94 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include #include @@ -13,6 +14,7 @@ #include "../common/UnitCell.h" class Logger; +namespace gemmi { struct Structure; } // Result of validating merged intensities against an atomic model. struct ModelValidationResult { @@ -99,6 +101,12 @@ struct ModelValidationResult { // identically, and the reflections themselves must not move. See the note in ValidateAgainstModel. bool adopted_model_enantiomorph = false; int model_space_group_number = 0; + + // The model as it was scored: re-fractionalized into the data cell and moved by the rigid-body + // step, everything else exactly as it was read. Null where the validation did not run. The frame + // it is written in is only settled by AdoptModelFrame below, so it is the caller that puts it on + // disk, with WritePlacedModel (WriteModel.h). + std::shared_ptr placed_model; }; // Given merged intensities and an atomic model (PDB or mmCIF, gzipped or not - the format is taken diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 8f8368f3f..875d30809 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4,6 +4,7 @@ #include #include "Rugnux.h" #include "ModelValidation.h" +#include "WriteModel.h" #include "SpotWidth.h" #include @@ -4833,6 +4834,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b r.h = h[0]; r.k = h[1]; r.l = h[2]; } } + // The coordinates that go with the maps just written, in the frame the reflection files + // beside them carry. Written whenever the validation ran, not only where the rigid-body + // step moved the model: it is also re-fractionalized into the data cell and may be + // relabelled to the model's enantiomorph, so "did not move" is not "identical to the input + // file". A rejected model is scored, placed and mapped like any other, and that is the case + // where someone most wants to see the model in the density, so it gets its file too. + if (validation.placed_model) + WritePlacedModel(*validation.placed_model, *result.consensus_cell, + experiment_.GetSpaceGroupOrP1(), config_.output_prefix, logger); } // Not on the geometry pre-pass. Pass 1 exists to choose the space group and post-refine the diff --git a/rugnux/WriteModel.cpp b/rugnux/WriteModel.cpp new file mode 100644 index 000000000..b9e36afbf --- /dev/null +++ b/rugnux/WriteModel.cpp @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "WriteModel.h" + +#include + +#include // make_mmcif_document +#include // write_cif_to_stream +#include // setup_entities + +#include "../common/Logger.h" + +void WritePlacedModel(const gemmi::Structure &placed, + const UnitCell &cell, + const gemmi::SpaceGroup &space_group, + const std::string &output_prefix, + Logger &logger) { + const std::string path = output_prefix + "_model.cif"; + + gemmi::Structure st = placed; + // Not the model's own cell and group but the data's, taken from the same two values the reflection + // files are written from. The coordinates already sit in this cell; what is set here is the label + // the file carries, which must match the .mtz beside it - the enantiomorph in particular, since + // --model can adopt the model's and that is neither the data's original label nor, necessarily, + // the one the input model arrived with. + st.cell = cell; + st.spacegroup_hm = space_group.xhm(); + st.setup_cell_images(); + // Fills in entity types and label_asym_id for a model read from a PDB, which carries neither. + // Both are no-ops where the input already had them, i.e. for an mmCIF input. + gemmi::setup_entities(st); + + std::ofstream os(path); + gemmi::cif::write_cif_to_stream(os, gemmi::make_mmcif_document(st)); + if (!os) { + logger.Error("Model validation: cannot write the placed model to {}", path); + return; + } + logger.Info("Model validation: the model as placed against these data written to {} " + "(cell and space group {} as in the reflection files)", path, space_group.short_name()); +} diff --git a/rugnux/WriteModel.h b/rugnux/WriteModel.h new file mode 100644 index 000000000..3600daabc --- /dev/null +++ b/rugnux/WriteModel.h @@ -0,0 +1,27 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "gemmi/model.hpp" +#include "gemmi/symmetry.hpp" + +#include "../common/UnitCell.h" + +class Logger; + +// Write `placed` - the atomic model as ValidateAgainstModel left it, re-fractionalized into the data +// cell and moved by the rigid-body step - to `output_prefix`_model.cif. Everything the input carried +// comes with it: chains, residues, ligands, waters, B-factors, occupancies, anisotropic Us. +// +// `cell` and `space_group` must be the ones the reflection files beside it are written in, i.e. the +// unit cell and DiffractionExperiment::GetSpaceGroupOrP1() that WriteReflections is given after +// AdoptModelFrame has settled the enantiomorph. A coordinate file in a different frame from the .mtz +// next to it is worse than no coordinate file at all. +void WritePlacedModel(const gemmi::Structure &placed, + const UnitCell &cell, + const gemmi::SpaceGroup &space_group, + const std::string &output_prefix, + Logger &logger); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 48cca5485..b6a8a1e25 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -44,6 +44,7 @@ #include "Rugnux.h" #include "RugnuxDefaults.h" #include "ModelValidation.h" +#include "WriteModel.h" #include "ResultReport.h" // Spots kept per image (the strongest ones) and handed to indexing. Offline reprocessing is not @@ -1810,6 +1811,15 @@ static int RunRugnux(int argc, char **argv) { r.h = h[0]; r.k = h[1]; r.l = h[2]; } } + // The coordinates that go with the maps just written, in the frame the reflection files + // beside them carry. Written whenever the validation ran, not only where the rigid-body + // step moved the model: it is also re-fractionalized into the data cell and may be + // relabelled to the model's enantiomorph, so "did not move" is not "identical to the input + // file". A rejected model is scored, placed and mapped like any other, and that is the case + // where someone most wants to see the model in the density, so it gets its file too. + if (validation.placed_model) + WritePlacedModel(*validation.placed_model, *experiment.GetUnitCell(), + experiment.GetSpaceGroupOrP1(), output_prefix, logger); } // Unmerged observations, from the integrated observations rather than the merged ones: the diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 1c96726c7..0fffadf6e 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -15,6 +15,7 @@ #include "../rugnux/ModelValidation.h" #include "../rugnux/RigidBodyRefine.h" #include "../rugnux/SigmaA.h" +#include "../rugnux/WriteModel.h" namespace { // A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with @@ -70,6 +71,16 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00 return pdb + "END\n"; } + // The same two atoms as kPdb, plus an anisotropic U, a partial occupancy and a water: a written + // model has to carry what the input carried, not only its coordinates. + const char *kPdbRich = + "CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n" + "ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n" + "ANISOU 1 CA GLY A 1 1000 1200 1400 100 200 300 C\n" + "ATOM 2 CB GLY A 1 12.000 14.000 16.000 0.60 25.00 C\n" + "HETATM 3 O HOH A 101 20.000 22.000 24.000 1.00 30.00 O\n" + "END\n"; + std::string WriteTemp(const std::string &name, const char *content) { std::ofstream f(name); f << content; @@ -211,3 +222,56 @@ TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValida CHECK(unweighted.weight[0].m == 1.0); CHECK(unweighted.weight[0].d == 1.0); } + +// The model rugnux scored has to reach disk, or a user overlaying their input model on rugnux's maps +// is wrong by the whole rigid-body shift. Two things have to hold: the file carries what the input +// carried, and it carries the cell and space group the reflection files beside it are written in - +// which, once --model has adopted the model's enantiomorph, is neither the data's original label nor +// necessarily the input model's. +TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]") { + Logger logger("WriteModel_KeepsTheContentAndTakesTheGivenFrame"); + + const auto input = WriteTemp("write_model_test_input.pdb", kPdbRich); + gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect); + + // A frame that is neither the model's (P 1, 40/50/60) nor anything derived from it: the tetragonal + // lysozyme cell and one enantiomorph of its group, standing in for what AdoptModelFrame settled. + const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2"); + REQUIRE(sg != nullptr); + + WritePlacedModel(st, data_cell, *sg, "write_model_test", logger); + const std::string written = "write_model_test_model.cif"; + REQUIRE(std::filesystem::exists(written)); + + const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect); + CHECK(back.cell.a == Catch::Approx(79.0)); + CHECK(back.cell.c == Catch::Approx(38.0)); + REQUIRE(back.find_spacegroup() != nullptr); + CHECK(back.find_spacegroup()->number == 96); + + REQUIRE(back.models.size() == 1); + REQUIRE(back.models[0].chains.size() == 1); + const gemmi::Chain &chain = back.models[0].chains[0]; + REQUIRE(chain.residues.size() == 2); // GLY and the water + REQUIRE(chain.residues[0].atoms.size() == 2); + REQUIRE(chain.residues[1].atoms.size() == 1); + + // The coordinates are written as they stand - WritePlacedModel changes the label, not the atoms. + const gemmi::Atom &ca = chain.residues[0].atoms[0]; + const gemmi::Atom &cb = chain.residues[0].atoms[1]; + CHECK(ca.pos.x == Catch::Approx(10.0)); + CHECK(ca.pos.z == Catch::Approx(14.0)); + CHECK(cb.pos.y == Catch::Approx(14.0)); + // B-factors, occupancies and the anisotropic U survive the round trip. + CHECK(ca.b_iso == Catch::Approx(20.0)); + CHECK(cb.b_iso == Catch::Approx(25.0)); + CHECK(cb.occ == Catch::Approx(0.60)); + CHECK(ca.aniso.nonzero()); + CHECK(ca.aniso.u11 == Catch::Approx(0.1000).margin(1e-4)); + CHECK(ca.aniso.u12 == Catch::Approx(0.0100).margin(1e-4)); + CHECK(chain.residues[1].name == "HOH"); + + std::filesystem::remove(input); + std::filesystem::remove(written); +}