Skip to content

Commit

Permalink
XdsAscii: read also merged files from XSCALE
Browse files Browse the repository at this point in the history 
xds2mtz converts it to 5-column merged MTZ
  • Loading branch information
@wojdyr
wojdyr committed Nov 19, 2024
1 parent 32b97aa commit 76f405e
Show file tree
Hide file tree
Showing 3 changed files with 156 additions and 128 deletions.
2 changes: 1 addition & 1 deletion include/gemmi/xds_ascii.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ inline bool likely_in_house_source(double wavelength) {
struct GEMMI_DLL XdsAscii {
struct Refl {
Miller hkl;
int iset = 1;
double iobs;
double sigma;
double xd;
Expand All @@ -35,7 +36,6 @@ struct GEMMI_DLL XdsAscii {
double peak;
double corr; // is it always integer?
double maxc;
int iset = 1;

// ZD can be negative for a few reflections
int frame() const { return (int) std::floor(zd + 1); }
Expand Down
241 changes: 129 additions & 112 deletions prog/xds2mtz.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,134 +156,151 @@ int GEMMI_MAIN(int argc, char **argv) {
mtz.datasets.push_back({iset.id, pxd[0], pxd[1], pxd[2], cell, wavelength});
}
}
mtz.add_column("M/ISYM", 'Y', 0, -1, false);
mtz.add_column("BATCH", 'B', 0, -1, false);
if (xds.read_columns >= 8) {
mtz.add_column("M/ISYM", 'Y', 0, -1, false);
mtz.add_column("BATCH", 'B', 0, -1, false);
}
mtz.add_column("I", 'J', 0, -1, false);
mtz.add_column("SIGI", 'Q', 0, -1, false);
mtz.add_column("XDET", 'R', 0, -1, false);
mtz.add_column("YDET", 'R', 0, -1, false);
mtz.add_column("ROT", 'R', 0, -1, false);
if (xds.read_columns >= 11) {
mtz.add_column("FRACTIONCALC", 'R', 0, -1, false);
mtz.add_column("LP", 'R', 0, -1, false);
mtz.add_column("CORR", 'R', 0, -1, false);
if (xds.read_columns > 11)
mtz.add_column("MAXC", 'I', 0, -1, false);
}
mtz.add_column("FLAG", 'I', 0, -1, false);
mtz.nreflections = (int) xds.data.size();
mtz.data.resize(mtz.columns.size() * xds.data.size());
gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
int max_frame = 0;
for (const gemmi::XdsAscii::Refl& refl : xds.data)
max_frame = std::max(max_frame, refl.frame());
int iset_offset = (max_frame + 11000) / 10000 * 10000;
// iset,frame -> batch
std::map<std::pair<int,int>, int> frames;
size_t k = 0;
for (const gemmi::XdsAscii::Refl& refl : xds.data) {
auto hkl = refl.hkl;
int isym = hkl_mover.move_to_asu(hkl);
for (size_t j = 0; j != 3; ++j)
mtz.data[k++] = (float) hkl[j];
mtz.data[k++] = (float) isym;
int frame = refl.frame();
int batch = frame + iset_offset * std::max(refl.iset - 1, 0);
frames.emplace(std::make_pair(refl.iset, frame), batch);
mtz.data[k++] = (float) batch;
mtz.data[k++] = (float) refl.iobs; // I
mtz.data[k++] = (float) std::fabs(refl.sigma); // SIGI
mtz.data[k++] = (float) refl.xd;
mtz.data[k++] = (float) refl.yd;
mtz.data[k++] = (float) xds.rot_angle(refl); // ROT
if (xds.read_columns >= 8) {
mtz.add_column("XDET", 'R', 0, -1, false);
mtz.add_column("YDET", 'R', 0, -1, false);
mtz.add_column("ROT", 'R', 0, -1, false);
if (xds.read_columns >= 11) {
mtz.data[k++] = float(0.01 * refl.peak); // FRACTIONCALC
mtz.data[k++] = (float) refl.rlp;
mtz.data[k++] = float(0.01 * refl.corr);
mtz.add_column("FRACTIONCALC", 'R', 0, -1, false);
mtz.add_column("LP", 'R', 0, -1, false);
mtz.add_column("CORR", 'R', 0, -1, false);
if (xds.read_columns > 11)
mtz.data[k++] = (float) refl.maxc;
mtz.add_column("MAXC", 'I', 0, -1, false);
}
mtz.data[k++] = refl.sigma < 0 ? 64.f : 0.f; // FLAG
mtz.add_column("FLAG", 'I', 0, -1, false);
}
// Prepare a similar batch header as Pointless.
gemmi::Mtz::Batch batch;
batch.set_dataset_id(1);
mtz.nreflections = (int) xds.data.size();
mtz.data.resize(mtz.columns.size() * xds.data.size());
if (xds.read_columns >= 8) { // unmerged file
gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
int max_frame = 0;
for (const gemmi::XdsAscii::Refl& refl : xds.data)
max_frame = std::max(max_frame, refl.frame());
int iset_offset = (max_frame + 11000) / 10000 * 10000;
// iset,frame -> batch
std::map<std::pair<int,int>, int> frames;
size_t k = 0;
for (const gemmi::XdsAscii::Refl& refl : xds.data) {
auto hkl = refl.hkl;
int isym = hkl_mover.move_to_asu(hkl);
for (size_t j = 0; j != 3; ++j)
mtz.data[k++] = (float) hkl[j];
mtz.data[k++] = (float) isym;
int frame = refl.frame();
int batch = frame + iset_offset * std::max(refl.iset - 1, 0);
frames.emplace(std::make_pair(refl.iset, frame), batch);
mtz.data[k++] = (float) batch;
mtz.data[k++] = (float) refl.iobs; // I
mtz.data[k++] = (float) std::fabs(refl.sigma); // SIGI
mtz.data[k++] = (float) refl.xd;
mtz.data[k++] = (float) refl.yd;
mtz.data[k++] = (float) xds.rot_angle(refl); // ROT
if (xds.read_columns >= 11) {
mtz.data[k++] = float(0.01 * refl.peak); // FRACTIONCALC
mtz.data[k++] = (float) refl.rlp;
mtz.data[k++] = float(0.01 * refl.corr);
if (xds.read_columns > 11)
mtz.data[k++] = (float) refl.maxc;
}
mtz.data[k++] = refl.sigma < 0 ? 64.f : 0.f; // FLAG
}
// Prepare a similar batch header as Pointless.
gemmi::Mtz::Batch batch;
batch.set_dataset_id(1);

// We don't set lbcell (refinement flags for unit cell),
// because it's probably not used by any program anyway.
// batch.ints[4] to [9] = left unset
// We don't set lbcell (refinement flags for unit cell),
// because it's probably not used by any program anyway.
// batch.ints[4] to [9] = left unset

// We also skip jumpax, which is defined as:
// reciprocal axis closest to principle goniostat axis E1
// batch.ints[11] = left unset
// We also skip jumpax, which is defined as:
// reciprocal axis closest to principle goniostat axis E1
// batch.ints[11] = left unset

// We assume one crystal, one goniostat axis, one detector.
batch.ints[12] = 1; // ncryst
batch.ints[14] = 2; // ldtype 3D
batch.ints[15] = 1; // jsaxs - goniostat scan axis number
batch.ints[17] = 1; // ngonax - number of goniostat axes
batch.ints[19] = 1; // ndet
// We assume one crystal, one goniostat axis, one detector.
batch.ints[12] = 1; // ncryst
batch.ints[14] = 2; // ldtype 3D
batch.ints[15] = 1; // jsaxs - goniostat scan axis number
batch.ints[17] = 1; // ngonax - number of goniostat axes
batch.ints[19] = 1; // ndet

batch.set_cell(mtz.cell); // batch.floats[0] to [5]
gemmi::Vec3 s0(-1, 0, 0); // will be re-set if we have geometry info
try {
gemmi::Mat33 Q = xds.calculate_conversion_from_cambridge().inverse();
s0 = -Q.multiply(xds.get_s0_direction());
// Orientation matrix U. It is calculated differently in Pointless,
// so the results are slightly different (due to limited precision
// of numbers in XDS file).
gemmi::Mat33 U = Q.multiply(xds.get_orientation());
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
batch.floats[6 + 3*i + j] = (float) U[j][i];
} catch (std::runtime_error& e) {
// if some of the headers are absent, U is not set
// and s0 is set to "idealised" s0
std::fprintf(stderr, "Note: %s, orientation matrix U not set.\n", e.what());
}
batch.floats[21] = float(xds.reflecting_range_esd); // crydat(0)
// In the so-called "Cambridge" frame (as used by Mosflm),
// the principal rotation axis is along z
// and the incident beam S0 is along x.
// Therefore, we set the rotation axis phi (scanax) to:
batch.floats[38+2] = 1.f; // scanax = [0, 0, 1]
batch.floats[47] = float(xds.oscillation_range); // phi range
// E1,E2,E3 vectors are the goniostat axes in Cambridge laboratory frame.
// E1 is set to rotation axis. E2 and E3 are not set for ngonax==1.
batch.floats[59+2] = 1.f; // e1 = scanax
batch.set_cell(mtz.cell); // batch.floats[0] to [5]
gemmi::Vec3 s0(-1, 0, 0); // will be re-set if we have geometry info
try {
gemmi::Mat33 Q = xds.calculate_conversion_from_cambridge().inverse();
s0 = -Q.multiply(xds.get_s0_direction());
// Orientation matrix U. It is calculated differently in Pointless,
// so the results are slightly different (due to limited precision
// of numbers in XDS file).
gemmi::Mat33 U = Q.multiply(xds.get_orientation());
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
batch.floats[6 + 3*i + j] = (float) U[j][i];
} catch (std::runtime_error& e) {
// if some of the headers are absent, U is not set
// and s0 is set to "idealised" s0
std::fprintf(stderr, "Note: %s, orientation matrix U not set.\n", e.what());
}
batch.floats[21] = float(xds.reflecting_range_esd); // crydat(0)
// In the so-called "Cambridge" frame (as used by Mosflm),
// the principal rotation axis is along z
// and the incident beam S0 is along x.
// Therefore, we set the rotation axis phi (scanax) to:
batch.floats[38+2] = 1.f; // scanax = [0, 0, 1]
batch.floats[47] = float(xds.oscillation_range); // phi range
// E1,E2,E3 vectors are the goniostat axes in Cambridge laboratory frame.
// E1 is set to rotation axis. E2 and E3 are not set for ngonax==1.
batch.floats[59+2] = 1.f; // e1 = scanax

// Idealised source vector is -x ([-1 0 0]), antiparallel to beam.
batch.floats[80+0] = -1.f; // source[0]
// Idealised source vector is -x ([-1 0 0]), antiparallel to beam.
batch.floats[80+0] = -1.f; // source[0]

// s0 source vector (including tilts), CMtz::MTZBAT::so in libccp4
batch.floats[83] = (float) s0.x;
batch.floats[84] = (float) s0.y;
batch.floats[85] = (float) s0.z;
// s0 source vector (including tilts), CMtz::MTZBAT::so in libccp4
batch.floats[83] = (float) s0.x;
batch.floats[84] = (float) s0.y;
batch.floats[85] = (float) s0.z;

batch.set_wavelength((float)xds.wavelength); // batch.floats[86]
// or batch.set_wavelength(iset.wavelength);
// Detector geometry.
batch.floats[111] = (float) xds.detector_distance; // dx[0]
batch.floats[113] = 1.f; // detlm[0][0][0]
batch.floats[114] = (float) xds.nx;
batch.floats[115] = 1.f;
batch.floats[116] = (float) xds.ny;
batch.axes.push_back("PHI"); // gonlab[0]
batch.set_wavelength((float)xds.wavelength); // batch.floats[86]
// or batch.set_wavelength(iset.wavelength);
// Detector geometry.
batch.floats[111] = (float) xds.detector_distance; // dx[0]
batch.floats[113] = 1.f; // detlm[0][0][0]
batch.floats[114] = (float) xds.nx;
batch.floats[115] = 1.f;
batch.floats[116] = (float) xds.ny;
batch.axes.push_back("PHI"); // gonlab[0]

for (auto& t : frames) {
batch.set_dataset_id(t.first.first);
batch.number = t.second;
int frame = t.first.second;
double phistt = xds.starting_angle +
xds.oscillation_range * (frame - xds.starting_frame);
batch.floats[36] = float(phistt);
batch.floats[37] = float(phistt + xds.oscillation_range); // phiend
mtz.batches.push_back(batch);
for (auto& t : frames) {
batch.set_dataset_id(t.first.first);
batch.number = t.second;
int frame = t.first.second;
double phistt = xds.starting_angle +
xds.oscillation_range * (frame - xds.starting_frame);
batch.floats[36] = float(phistt);
batch.floats[37] = float(phistt + xds.oscillation_range); // phiend
mtz.batches.push_back(batch);
}
mtz.sort(5);
} else { // merged
if (verbose)
std::fprintf(stderr, "Preparing merged MTZ file...\n");
size_t k = 0;
for (const gemmi::XdsAscii::Refl& refl : xds.data) {
for (size_t j = 0; j != 3; ++j)
mtz.data[k++] = (float) refl.hkl[j];
mtz.data[k++] = (float) refl.iobs; // I
mtz.data[k++] = (float) std::fabs(refl.sigma); // SIGI
}
mtz.sort(3);
}
mtz.sort(5);
if (verbose)
std::fprintf(stderr, "Writing %d reflections to %s ...\n",
mtz.nreflections, output_path);
std::fprintf(stderr, "Writing %d reflections to (%smerged) %s ...\n",
mtz.nreflections, (xds.read_columns >= 8 ? "un" : ""), output_path);
mtz.write_to_file(output_path);
} catch (std::exception& e) {
std::fprintf(stderr, "ERROR: %s\n", e.what());
Expand Down
41 changes: 26 additions & 15 deletions src/xds_ascii.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,10 +117,15 @@ void XdsAscii::read_stream(LineReaderBase&& line_reader, const std::string& sour
read_columns = 12;
char line[256];
size_t len0 = line_reader.copy_line(line, 255);
if (len0 == 0)
fail("empty file");
int iset_col = 0;
if (len0 == 0 || !(starts_with(line, "!FORMAT=XDS_ASCII MERGE=FALSE") ||
(starts_with(line, "!OUTPUT_FILE=INTEGRATE.HKL"))))
fail("not an unmerged XDS_ASCII nor INTEGRATE.HKL file: " + source_path);
const char xds_ascii_header[] = "!FORMAT=XDS_ASCII MERGE=";
char xds_ascii_type = '\0';
if (starts_with(line, xds_ascii_header))
xds_ascii_type = line[sizeof(xds_ascii_header)-1];
if (!xds_ascii_type && !starts_with(line, "!OUTPUT_FILE=INTEGRATE.HKL"))
fail("not an XDS_ASCII nor INTEGRATE.HKL file: " + source_path);
const char* rhs;
while (size_t len = line_reader.copy_line(line, 255)) {
if (line[0] == '!') {
Expand Down Expand Up @@ -195,7 +200,9 @@ void XdsAscii::read_stream(LineReaderBase&& line_reader, const std::string& sour
} else if (starts_with_ptr(line+1, "NUMBER_OF_ITEMS_IN_EACH_DATA_RECORD=", &rhs)) {
int num = simple_atoi(rhs);
// INTEGRATE.HKL has read_columns=12, as set above
if (generated_by == "XSCALE")
if (xds_ascii_type == 'T') // merged file
read_columns = 5;
else if (generated_by == "XSCALE")
read_columns = 8;
else if (generated_by == "CORRECT")
read_columns = 11;
Expand Down Expand Up @@ -239,20 +246,24 @@ void XdsAscii::read_stream(LineReaderBase&& line_reader, const std::string& sour
r.hkl[i] = simple_atoi(p, &p);
auto result = fast_from_chars(p, line+len, r.iobs); // 4
result = fast_from_chars(result.ptr, line+len, r.sigma); // 5
result = fast_from_chars(result.ptr, line+len, r.xd); // 6
result = fast_from_chars(result.ptr, line+len, r.yd); // 7
result = fast_from_chars(result.ptr, line+len, r.zd); // 8
if (read_columns >= 11) {
result = fast_from_chars(result.ptr, line+len, r.rlp); // 9
result = fast_from_chars(result.ptr, line+len, r.peak); // 10
result = fast_from_chars(result.ptr, line+len, r.corr); // 11
if (read_columns > 11) {
result = fast_from_chars(result.ptr, line+len, r.maxc); // 12
if (read_columns >= 8) {
result = fast_from_chars(result.ptr, line+len, r.xd); // 6
result = fast_from_chars(result.ptr, line+len, r.yd); // 7
result = fast_from_chars(result.ptr, line+len, r.zd); // 8
if (read_columns >= 11) {
result = fast_from_chars(result.ptr, line+len, r.rlp); // 9
result = fast_from_chars(result.ptr, line+len, r.peak); // 10
result = fast_from_chars(result.ptr, line+len, r.corr); // 11
if (read_columns >= 12) {
result = fast_from_chars(result.ptr, line+len, r.maxc); // 12
} else {
r.maxc = 0; // 12
}
} else {
r.maxc = 0;
r.rlp = r.peak = r.corr = r.maxc = 0; // 9-11
}
} else {
r.rlp = r.peak = r.corr = r.maxc = 0;
r.xd = r.yd = r.zd = 0; // 6-8
}
if (result.ec != std::errc())
fail("failed to parse data line:\n", line);
Expand Down

0 comments on commit 76f405e

Please sign in to comment.