Pairs being lost

[tdudgeon]

When reading a Buster output CIF file I'm finding different behaviour in the way the file is read.

Originally I was getting the Document object using

doc = cif.read(input)

and the _refine pairs that are defined in the file, as expected, are all present in the Document.

Then I had to change to using:

struc = gemmi.read_structure(input)
doc = struc.make_mmcif_document()

as I needed to perform operations on the Structure before generating the Document.
But this way results in quite different behaviour.
Firstly the _refine data is now a one-row loop not pairs (surprising, but explainable), but more problematical a lot of the _refine data is missing.

In the source pairs I have:

_refine.entry_id
_refine.pdbx_refine_id
_refine.ls_d_res_high
_refine.ls_d_res_low'
_refine.ls_percent_reflns_obs
_refine.ls_number_reflns_obs
_refine.ls_number_reflns_R_free
_refine.ls_R_factor_obs
_refine.ls_R_factor_R_work
_refine.ls_R_factor_R_free
_refine.pdbx_ls_cross_valid_method
_refine.pdbx_R_Free_selection_details
_refine.B_iso_mean
_refine.aniso_B[1][1]
_refine.aniso_B[2][2]
_refine.aniso_B[3][3]_refine.aniso_B[1][2]
_refine.aniso_B[1][3]
_refine.aniso_B[2][3]
_refine.pdbx_overall_SU_R_Blow_DPI
_refine.pdbx_overall_SU_R_free_Blow_DPI
_refine.overall_SU_R_Cruickshank_DPI
_refine.pdbx_overall_SU_R_free_Cruickshank_DPI
_refine.correlation_coeff_Fo_to_Fc
_refine.correlation_coeff_Fo_to_Fc_free

But in the Document generated via the Structure I only have these in the loop:

_refine.entry_id
_refine.pdbx_refine_id
_refine.ls_d_res_high
_refine.ls_d_res_low
_refine.ls_percent_reflns_obs
_refine.ls_number_reflns_obs
_refine.ls_number_reflns_R_work
_refine.ls_number_reflns_R_free
_refine.ls_R_factor_obs
_refine.ls_R_factor_R_work
_refine.ls_R_factor_R_free
_refine.aniso_B[1][1]
_refine.aniso_B[2][2]
_refine.aniso_B[3][3]
_refine.aniso_B[1][2]
_refine.aniso_B[1][3]
_refine.aniso_B[2][3]

What is happening here? Presumably a bug?

Example file:
BUSTER_model.cif.txt

Example code:
cif_reader.py.txt

[tdudgeon]

Sorry, I added the wrong code file. Try this:

handle_cif.py.txt

[tdudgeon]

My AI friend tells me that this is the problem:

The bug is a read/write asymmetry in the round-trip gemmi.read_structure() → make_mmcif_document().

Where the data is lost: src/mmcif.cpp, read_refinement_info() (line 667)

  for (auto row : block.find("_refine.", {"pdbx_refine_id",           // 0
                                          "?ls_d_res_high",           // 1
                                          "?ls_d_res_low",            // 2
                                          "?ls_percent_reflns_obs",   // 3
                                          "?ls_number_reflns_obs",    // 4
                                          "?ls_number_reflns_R_work", // 5
                                          "?ls_number_reflns_R_free", // 6
                                          "?ls_R_factor_obs",         // 7
                                          "?ls_R_factor_R_work",      // 8
                                          "?ls_R_factor_R_free"})) {  // 9

This only parses 10 fields. The RefinementInfo struct holds many more, and to_mmcif.cpp (lines 991–1030) conditionally writes
them all — but only when st.meta.has(field) returns true, i.e., when the value is non-NAN / non-empty. Since the reader never
populates them, has() always returns false and they are silently dropped.

Missing field mappings (CIF tag → struct member, all absent from reader):

  │                CIF tag                 │  RefinementInfo member  │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ pdbx_ls_cross_valid_method             │ cross_validation_method │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ pdbx_R_Free_selection_details          │ rfree_selection_method  │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ B_iso_mean                             │ mean_b                  │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ aniso_B[1][1] … aniso_B[2][3]          │ aniso_b (SMat33)        │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ pdbx_overall_SU_R_Blow_DPI             │ dpi_blow_r              │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ pdbx_overall_SU_R_free_Blow_DPI        │ dpi_blow_rfree          │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ overall_SU_R_Cruickshank_DPI           │ dpi_cruickshank_r       │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ pdbx_overall_SU_R_free_Cruickshank_DPI │ dpi_cruickshank_rfree   │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ correlation_coeff_Fo_to_Fc             │ cc_fo_fc_work           │
  ├────────────────────────────────────────┼─────────────────────────┤
  │ correlation_coeff_Fo_to_Fc_free        │ cc_fo_fc_free           │
  └────────────────────────────────────────┴─────────────────────────┘

All of these are present in your BUSTER_model.cif (lines 71–85) and are parsed correctly by cif.read() (which returns the raw
document without going through RefinementInfo). They're only lost when the round-trip goes through the Structure object.

The fix

Extend read_refinement_info to request these as optional columns (indices 10–24) and populate the struct fields:

  for (auto row : block.find("_refine.", {
          "pdbx_refine_id",                           // 0
          "?ls_d_res_high",                           // 1
          "?ls_d_res_low",                            // 2
          "?ls_percent_reflns_obs",                   // 3
          "?ls_number_reflns_obs",                    // 4
          "?ls_number_reflns_R_work",                 // 5
          "?ls_number_reflns_R_free",                 // 6
          "?ls_R_factor_obs",                         // 7
          "?ls_R_factor_R_work",                      // 8
          "?ls_R_factor_R_free",                      // 9
          "?pdbx_ls_cross_valid_method",              // 10
          "?pdbx_R_Free_selection_details",           // 11
          "?B_iso_mean",                              // 12
          "?aniso_B[1][1]", "?aniso_B[2][2]",        // 13, 14
          "?aniso_B[3][3]", "?aniso_B[1][2]",        // 15, 16
          "?aniso_B[1][3]", "?aniso_B[2][3]",        // 17, 18
          "?pdbx_overall_SU_R_Blow_DPI",             // 19
          "?pdbx_overall_SU_R_free_Blow_DPI",        // 20
          "?overall_SU_R_Cruickshank_DPI",           // 21
          "?pdbx_overall_SU_R_free_Cruickshank_DPI", // 22
          "?correlation_coeff_Fo_to_Fc",             // 23
          "?correlation_coeff_Fo_to_Fc_free"})) {    // 24
      // ... existing code for indices 0-9 ...
      copy_string(row, 10, ref.cross_validation_method);
      copy_string(row, 11, ref.rfree_selection_method);
      copy_double(row, 12, ref.mean_b);
      if (row.has2(13))
          ref.aniso_b = get_smat33<double>(row, 13); // reads 13..18
      copy_double(row, 19, ref.dpi_blow_r);
      copy_double(row, 20, ref.dpi_blow_rfree);
      copy_double(row, 21, ref.dpi_cruickshank_r);
      copy_double(row, 22, ref.dpi_cruickshank_rfree);
      copy_double(row, 23, ref.cc_fo_fc_work);
      copy_double(row, 24, ref.cc_fo_fc_free);
  }

The existing get_smat33(row, n) helper (line 31 in mmcif.cpp) reads 6 consecutive columns from index n, which fits perfectly
for aniso_b at indices 13–18.

[keitaroyam]

I think we can add reading/writing of more refine-related metadata. For now, you could alternatively do this, unless any metadata are modified:

import gemmi
st = gemmi.read_structure("BUSTER_model.cif")
doc = gemmi.cif.read("BUSTER_model.cif")
# ... do operations as needed
groups = gemmi.MmcifOutputGroups(False, atoms=True)
st.update_mmcif_block(doc[0], groups)
doc.write_file("modified.cif")

In this way, only _atom_site and _atom_site_anisotrop will be updated within the cif file.

Please refer to the documentation regarding MmcifOutputGroups: https://gemmi.readthedocs.io/en/stable/mol.html#id2