dials.export: Add option to specify composition of asymmetric unit so…

… a complete SHELX instruction file can be written. (dials#2623)


---------

Co-authored-by: Ben Williams <benjaminhwilliams@users.noreply.github.com>

newsfragments/2623.feature

@@ -0,0 +1 @@
+``dials.export``: Add option to specify composition of asymmetric unit for SHELX ``.ins`` file output.

src/dials/command_line/export.py

@@ -45,6 +45,11 @@
 XDS format exports a models.expt file as XDS.INP and XPARM.XDS files. If a
 reflection file is given it will be exported as a SPOT.XDS file.
 
+SHELX format exports intensity data in HKLF 4 format for use in the SHELX suite
+of programs. As this file format does not contain unit cell parameters or
+symmetry information a minimal instruction file is also written. Optionally the
+expected contents of the asymmetric unit can be added to this instruciton file.
+
 PETS format exports intensity data and diffraction data in the CIF format
 used by PETS. This is primarily intended to produce files suitable for
 dynamic diffraction refinement using Jana2020, which requires this format.
@@ -72,6 +77,11 @@
   dials.export indexed.refl format=xds
   dials.export models.expt format=xds
   dials.export models.expt indexed.refl format=xds
+
+  # Export to shelx
+  dials.export scaled.expt scaled.refl format=shelx
+  dials.export scaled.expt scaled.refl format=shelx shelx.hklout=dials.hkl
+  dials.export scaled.expt scaled.refl format=shelx composition=C3H7NO2S
 """
 
 phil_scope = parse(
@@ -249,6 +259,9 @@
     ins = dials.ins
       .type = path
       .help = "The output ins file"
+    composition = CH
+      .type = str
+      .help = "The chemical composition of the asymmetric unit"
     scale = True
       .type = bool
       .help = "Scale reflections to maximise output precision in SHELX 8.2f format"

src/dials/util/export_shelx.py

@@ -1,9 +1,11 @@
 from __future__ import annotations
 
 import logging
+import re
 from math import isclose
 
 from cctbx import crystal, miller
+from cctbx.eltbx import chemical_elements
 from iotbx.shelx.write_ins import LATT_SYMM
 
 from dials.algorithms.scaling.scaling_library import determine_best_unit_cell
@@ -96,12 +98,13 @@ def export_shelx(scaled_data, experiment_list, params):
     _write_ins(
         experiment_list,
         best_unit_cell=params.mtz.best_unit_cell,
+        composition=params.shelx.composition,
         ins_file=params.shelx.ins,
     )
     logger.info(f"Written {params.shelx.ins}")
 
 
-def _write_ins(experiment_list, best_unit_cell, ins_file):
+def _write_ins(experiment_list, best_unit_cell, composition, ins_file):
     sg = experiment_list[0].crystal.get_space_group()
     unit_cells = []
     wavelengths = []
@@ -167,3 +170,36 @@ def _write_ins(experiment_list, best_unit_cell, ins_file):
                 )
             )
         LATT_SYMM(f, sg)
+        logger.info(
+            f"Using composition {composition} to write SFAC & UNIT instructions"
+        )
+        sorted_composition = _parse_compound(composition)
+        f.write("SFAC %s\n" % " ".join(sorted_composition))
+        f.write(
+            "UNIT %s\n"
+            % " ".join(
+                [str(sorted_composition[e] * sg.order_z()) for e in sorted_composition]
+            )
+        )
+        f.write("HKLF 4\n")  # Intensities
+        f.write("END\n")
+
+
+def _parse_compound(composition):
+    elements = chemical_elements.proper_caps_list()[:94]
+    m = re.findall(r"([A-Z][a-z]?)(\d*)", composition)
+    if all(e[1] == "" for e in m):
+        # Assume user has just supplied list of elements so UNIT instruction cannot be calculated
+        result = {element: 0 for element, count in m}
+    else:
+        result = {element: int(count or 1) for element, count in m}
+    # Check for elements without scattering factors in SHELXL
+    unmatched = list(set(result).difference(elements))
+    if unmatched:
+        raise Sorry(f"Element(s) {' '.join(unmatched)} in {composition} not recognised")
+
+    if "C" in result:
+        # move C to start of list
+        elements.insert(0, elements.pop(5))
+    sorted_result = {e: result[e] for e in elements if e in result}
+    return sorted_result

tests/command_line/test_export.py

@@ -584,6 +584,39 @@ def test_shelx_ins_best_unit_cell(dials_data, tmp_path):
                 assert result == pytest.approx(cell_esds[instruction], abs=0.001)
 
 
+def test_shelx_ins_composition(dials_data, tmp_path):
+    # Call dials.export
+    result = subprocess.run(
+        [
+            shutil.which("dials.export"),
+            "intensity=scale",
+            "format=shelx",
+            "composition=C3H7NO2S",
+            dials_data("l_cysteine_4_sweeps_scaled", pathlib=True)
+            / "scaled_20_25.expt",
+            dials_data("l_cysteine_4_sweeps_scaled", pathlib=True)
+            / "scaled_20_25.refl",
+        ],
+        cwd=tmp_path,
+        capture_output=True,
+    )
+    assert not result.returncode and not result.stderr
+    assert (tmp_path / "dials.ins").is_file()
+
+    sfac_unit = {
+        "SFAC": "C H N O S",
+        "UNIT": "12 28 4 8 4",
+    }
+
+    with (tmp_path / "dials.ins").open() as fh:
+        for line in fh:
+            tokens = line.split()
+            instruction = tokens[0]
+            if instruction in sfac_unit:
+                result = " ".join(tokens[1:6])
+                assert result == sfac_unit[instruction]
+
+
 def test_export_sum_or_profile_only(dials_data, tmp_path):
     expt = dials_data("insulin_processed", pathlib=True) / "integrated.expt"
     refl = dials_data("insulin_processed", pathlib=True) / "integrated.refl"