Skip to content

to_pdb: write REMARK 3 / 200 / 230 / 240 from Metadata #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CV-GPhL wants to merge 1 commit into
base: master
Choose a base branch
from
Open

to_pdb: write REMARK 3 / 200 / 230 / 240 from Metadata #2

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
130 changes: 130 additions & 0 deletions src/to_pdb.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,12 +197,142 @@ void write_ncs_op(const NcsOp& op, std::ostream& os) {
}
}

// Inverse of read_remark3_line() in pdb.cpp: emit REMARK 3 (refinement) from
// Metadata. gemmi reads these fields on input but, until now, dropped them when
// writing PDB. Key strings must match the parser exactly so it round-trips.
void write_remark_3(const Metadata& meta, std::ostream& os) {
if (meta.refinement.empty())
return;
char buf[88];
const RefinementInfo& ref = meta.refinement[0];
WRITELN("REMARK 3");
WRITELN("REMARK 3 REFINEMENT.");
for (const SoftwareItem& s : meta.software)
if (s.classification == SoftwareItem::Refinement) {
std::string v = s.name;
if (!s.version.empty()) { v += ' '; v += s.version; }
WRITELN("REMARK 3 PROGRAM : %s", v.c_str());
}
WRITELN("REMARK 3");
WRITELN("REMARK 3 DATA USED IN REFINEMENT.");
if (!std::isnan(ref.resolution_high))
WRITELN("REMARK 3 RESOLUTION RANGE HIGH (ANGSTROMS) : %.2f", ref.resolution_high);
if (!std::isnan(ref.resolution_low))
WRITELN("REMARK 3 RESOLUTION RANGE LOW (ANGSTROMS) : %.2f", ref.resolution_low);
if (!std::isnan(ref.completeness))
WRITELN("REMARK 3 COMPLETENESS FOR RANGE (%%) : %.1f", ref.completeness);
if (ref.reflection_count != -1)
WRITELN("REMARK 3 NUMBER OF REFLECTIONS : %d", ref.reflection_count);
WRITELN("REMARK 3");
WRITELN("REMARK 3 FIT TO DATA USED IN REFINEMENT.");
if (!ref.cross_validation_method.empty())
WRITELN("REMARK 3 CROSS-VALIDATION METHOD : %s", ref.cross_validation_method.c_str());
if (!ref.rfree_selection_method.empty())
WRITELN("REMARK 3 FREE R VALUE TEST SET SELECTION : %s", ref.rfree_selection_method.c_str());
if (!std::isnan(ref.r_all))
WRITELN("REMARK 3 R VALUE (WORKING + TEST SET) : %.4f", ref.r_all);
if (!std::isnan(ref.r_work))
WRITELN("REMARK 3 R VALUE (WORKING SET) : %.4f", ref.r_work);
if (!std::isnan(ref.r_free))
WRITELN("REMARK 3 FREE R VALUE : %.4f", ref.r_free);
if (ref.rfree_set_count != -1)
WRITELN("REMARK 3 FREE R VALUE TEST SET COUNT : %d", ref.rfree_set_count);
if (!std::isnan(ref.mean_b)) {
WRITELN("REMARK 3");
WRITELN("REMARK 3 B VALUES.");
WRITELN("REMARK 3 MEAN B VALUE (OVERALL, A**2) : %.2f", ref.mean_b);
}
}
Comment on lines +203 to +245

@gemini-code-assist gemini-code-assist Bot Jul 2, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

In joint refinement structures (such as X-ray + neutron, e.g., 5moo), there can be multiple refinement blocks in meta.refinement. Currently, write_remark_3 only writes the first refinement block (meta.refinement[0]), which silently drops the refinement details of any subsequent blocks. Looping over meta.refinement for the data-collection and fit sections resolves this issue while keeping the program header at the top.

void write_remark_3(const Metadata& meta, std::ostream& os) {
  if (meta.refinement.empty())
    return;
  char buf[88];
  WRITELN("REMARK   3");
  WRITELN("REMARK   3  REFINEMENT.");
  for (const SoftwareItem& s : meta.software)
    if (s.classification == SoftwareItem::Refinement) {
      std::string v = s.name;
      if (!s.version.empty()) { v += ' '; v += s.version; }
      WRITELN("REMARK   3   PROGRAM     : %s", v.c_str());
    }
  for (const RefinementInfo& ref : meta.refinement) {
    WRITELN("REMARK   3");
    WRITELN("REMARK   3  DATA USED IN REFINEMENT.");
    if (!std::isnan(ref.resolution_high))
      WRITELN("REMARK   3   RESOLUTION RANGE HIGH (ANGSTROMS) : %.2f", ref.resolution_high);
    if (!std::isnan(ref.resolution_low))
      WRITELN("REMARK   3   RESOLUTION RANGE LOW  (ANGSTROMS) : %.2f", ref.resolution_low);
    if (!std::isnan(ref.completeness))
      WRITELN("REMARK   3   COMPLETENESS FOR RANGE        (%%) : %.1f", ref.completeness);
    if (ref.reflection_count != -1)
      WRITELN("REMARK   3   NUMBER OF REFLECTIONS : %d", ref.reflection_count);
    WRITELN("REMARK   3");
    WRITELN("REMARK   3  FIT TO DATA USED IN REFINEMENT.");
    if (!ref.cross_validation_method.empty())
      WRITELN("REMARK   3   CROSS-VALIDATION METHOD : %s", ref.cross_validation_method.c_str());
    if (!ref.rfree_selection_method.empty())
      WRITELN("REMARK   3   FREE R VALUE TEST SET SELECTION : %s", ref.rfree_selection_method.c_str());
    if (!std::isnan(ref.r_all))
      WRITELN("REMARK   3   R VALUE     (WORKING + TEST SET) : %.4f", ref.r_all);
    if (!std::isnan(ref.r_work))
      WRITELN("REMARK   3   R VALUE            (WORKING SET) : %.4f", ref.r_work);
    if (!std::isnan(ref.r_free))
      WRITELN("REMARK   3   FREE R VALUE : %.4f", ref.r_free);
    if (ref.rfree_set_count != -1)
      WRITELN("REMARK   3   FREE R VALUE TEST SET COUNT : %d", ref.rfree_set_count);
    if (!std::isnan(ref.mean_b)) {
      WRITELN("REMARK   3");
      WRITELN("REMARK   3  B VALUES.");
      WRITELN("REMARK   3   MEAN B VALUE      (OVERALL, A**2) : %.2f", ref.mean_b);
    }
  }
}


// Inverse of read_remark_200_230_240() in pdb.cpp: emit data-collection
// REMARKs. Uses REMARK 200/230/240 for x-ray/neutron/electron respectively
// (the header record number is what the parser keys scattering type on).
void write_remark_200_230_240(const Metadata& meta, std::ostream& os) {
char buf[88];
for (size_t i = 0; i < meta.experiments.size() && i < meta.crystals.size(); ++i) {
const ExperimentInfo& e = meta.experiments[i];
const CrystalInfo& c = meta.crystals[i];
if (c.diffractions.empty())
continue;
const DiffractionInfo& d = c.diffractions[0];
Comment on lines +250 to +257

@gemini-code-assist gemini-code-assist Bot Jul 2, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

In structures sourced from mmCIF, meta.experiments and meta.crystals are not guaranteed to be aligned by index or have the same size. Using a simple index-based loop (meta.experiments[i] and meta.crystals[i]) can lead to mismatched data or out-of-bounds access. Instead, map the experiment to the correct crystal and diffraction info using diffraction_ids and id matching, with a safe fallback to the first available crystal if no match is found.

void write_remark_200_230_240(const Metadata& meta, std::ostream& os) {
  char buf[88];
  for (const ExperimentInfo& e : meta.experiments) {
    const CrystalInfo* matched_c = nullptr;
    const DiffractionInfo* matched_d = nullptr;
    for (const std::string& diff_id : e.diffraction_ids) {
      for (const CrystalInfo& c : meta.crystals) {
        for (const DiffractionInfo& d : c.diffractions) {
          if (d.id == diff_id) {
            matched_c = &c;
            matched_d = &d;
            break;
          }
        }
        if (matched_d) break;
      }
      if (matched_d) break;
    }
    if (!matched_d && !meta.crystals.empty()) {
      for (const CrystalInfo& c : meta.crystals) {
        if (!c.diffractions.empty()) {
          matched_c = &c;
          matched_d = &c.diffractions[0];
          break;
        }
      }
    }
    if (!matched_d)
      continue;
    const CrystalInfo& c = *matched_c;
    const DiffractionInfo& d = *matched_d;

const char* R = "REMARK 200";
if (d.scattering_type == "neutron") R = "REMARK 230";
else if (d.scattering_type == "electron") R = "REMARK 240";
WRITELN("%s", R);
WRITELN("%s EXPERIMENTAL DETAILS", R);
if (!e.method.empty())
WRITELN("%s EXPERIMENT TYPE : %s", R, e.method.c_str());
if (!std::isnan(d.temperature))
WRITELN("%s TEMPERATURE (KELVIN) : %g", R, d.temperature);
if (e.number_of_crystals != -1)
WRITELN("%s NUMBER OF CRYSTALS USED : %d", R, e.number_of_crystals);
if (!std::isnan(c.ph))
WRITELN("%s PH : %g", R, c.ph);
else if (!c.ph_range.empty())
WRITELN("%s PH : %s", R, c.ph_range.c_str());
if (d.source == "SYNCHROTRON") {
WRITELN("%s SYNCHROTRON (Y/N) : Y", R);
if (!d.synchrotron.empty())
WRITELN("%s RADIATION SOURCE : %s", R, d.synchrotron.c_str());
if (!d.beamline.empty())
WRITELN("%s BEAMLINE : %s", R, d.beamline.c_str());
} else if (!d.source.empty() && d.scattering_type != "neutron") {
WRITELN("%s RADIATION SOURCE : %s", R, d.source.c_str());
if (!d.source_type.empty())
WRITELN("%s X-RAY GENERATOR MODEL : %s", R, d.source_type.c_str());
} else if (!d.source.empty()) {
WRITELN("%s NEUTRON SOURCE : %s", R, d.source.c_str());
}
if (d.mono_or_laue)
WRITELN("%s MONOCHROMATIC OR LAUE (M/L) : %c", R, d.mono_or_laue);
if (!d.wavelengths.empty())
WRITELN("%s WAVELENGTH OR RANGE (A) : %s", R, d.wavelengths.c_str());
if (!d.monochromator.empty())
WRITELN("%s MONOCHROMATOR : %s", R, d.monochromator.c_str());
if (!d.detector.empty())
WRITELN("%s DETECTOR TYPE : %s", R, d.detector.c_str());
if (!d.detector_make.empty())
WRITELN("%s DETECTOR MANUFACTURER : %s", R, d.detector_make.c_str());
for (const SoftwareItem& s : meta.software) {
if (s.classification == SoftwareItem::DataReduction)
WRITELN("%s INTENSITY-INTEGRATION SOFTWARE : %s", R, s.name.c_str());
else if (s.classification == SoftwareItem::DataScaling)
WRITELN("%s DATA SCALING SOFTWARE : %s", R, s.name.c_str());
else if (s.classification == SoftwareItem::Phasing)
WRITELN("%s SOFTWARE USED : %s", R, s.name.c_str());
}
if (!meta.solved_by.empty())
WRITELN("%s METHOD USED TO DETERMINE THE STRUCTURE : %s", R, meta.solved_by.c_str());
if (!meta.starting_model.empty())
WRITELN("%s STARTING MODEL : %s", R, meta.starting_model.c_str());
const ReflectionsInfo& rf = e.reflections;
if (e.unique_reflections != -1)
WRITELN("%s NUMBER OF UNIQUE REFLECTIONS : %d", R, e.unique_reflections);
if (!std::isnan(rf.resolution_high))
WRITELN("%s RESOLUTION RANGE HIGH (A) : %.3f", R, rf.resolution_high);
if (!std::isnan(rf.resolution_low))
WRITELN("%s RESOLUTION RANGE LOW (A) : %.3f", R, rf.resolution_low);
if (!std::isnan(rf.completeness))
WRITELN("%s COMPLETENESS FOR RANGE (%%) : %.1f", R, rf.completeness);
if (!std::isnan(rf.redundancy))
WRITELN("%s DATA REDUNDANCY : %.3f", R, rf.redundancy);
if (!std::isnan(rf.r_merge))
WRITELN("%s R MERGE (I) : %.5f", R, rf.r_merge);
if (!std::isnan(rf.r_sym))
WRITELN("%s R SYM (I) : %.5f", R, rf.r_sym);
if (!std::isnan(rf.mean_I_over_sigma))
WRITELN("%s <I/SIGMA(I)> FOR THE DATA SET : %.4f", R, rf.mean_I_over_sigma);
}
}

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.");
}
write_remark_3(st.meta, os);
write_remark_200_230_240(st.meta, os);
if (!st.assemblies.empty()) {
const char* preface[] = {
"REMARK 350",
Expand Down
Loading