Merge branch 'af-mr4' into rc175: --af-model likelihood rescoring, PAE domain search, twin comparison on the subgroup merge, refined tNCS rule
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
+15
-1
@@ -473,7 +473,21 @@ complex crystals" (2003), Protein Sci. 12, 1865-1871
|
||||
[doi:10.1110/ps.0350503](https://doi.org/10.1110/ps.0350503). A copy refined beside the others has the
|
||||
others' placements read off by least-squares superposition: W. Kabsch, "A solution for the best rotation
|
||||
to relate two sets of vectors" (1976), Acta Cryst. A32, 922-923
|
||||
[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873).
|
||||
[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873). The placements are judged by
|
||||
their log-likelihood gain, from the Rice and Woolfson distributions of the observed normalised amplitudes
|
||||
given the model's, with sigma_A from the model's expected coordinate error: R. J. Read, "Pushing the
|
||||
boundaries of molecular replacement with maximum likelihood" (2001), Acta Cryst. D57, 1373-1382
|
||||
[doi:10.1107/S0907444901012471](https://doi.org/10.1107/S0907444901012471); and the gain is read as a
|
||||
Z-score over the same orientation at other positions, as the likelihood-enhanced translation function's
|
||||
is: A. J. McCoy, R. W. Grosse-Kunstleve, L. C. Storoni and R. J. Read, "Likelihood-enhanced fast
|
||||
translation functions" (2005), Acta Cryst. D61, 458-464
|
||||
[doi:10.1107/S0907444905001617](https://doi.org/10.1107/S0907444905001617). The logarithm of the Bessel
|
||||
function I0 in it is evaluated with the polynomial approximations of M. Abramowitz and I. A. Stegun,
|
||||
"Handbook of Mathematical Functions" (1964), National Bureau of Standards, 9.8.1-9.8.2. A predicted
|
||||
model is cut into rigid domains by its predicted aligned error, as Oeffner et al. (above) do before
|
||||
molecular replacement; the error itself is AlphaFold's: J. Jumper, R. Evans, A. Pritzel et al., "Highly
|
||||
accurate protein structure prediction with AlphaFold" (2021), Nature 596, 583-589
|
||||
[doi:10.1038/s41586-021-03819-2](https://doi.org/10.1038/s41586-021-03819-2).
|
||||
|
||||
## Software and computing methods
|
||||
|
||||
|
||||
+39
-14
@@ -980,14 +980,30 @@ merge cannot tell apart, and the best three placements per group that pack throu
|
||||
solvent and R-free (to `SCORE_D_MIN`, 3 A or the data's limit). Further copies follow one at a time, up
|
||||
to the largest copy number with a Matthews probability of 1% (`COPIES_LIKELY`: the most probable): the
|
||||
translation function again with the copies already placed as a known partial structure, its peaks
|
||||
ranked by the correlation of all the copies with the data, each new copy refined beside the others held
|
||||
fixed, and kept where its `TF_Z` reaches 8, it packs and R-free drops by 0.01 (`COPIES`, and
|
||||
`MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched side by side; every group the best
|
||||
one leads is searched once more with the best one's copies' orientations tried first, and the better
|
||||
search kept.
|
||||
ranked by the likelihood gain of the new copy beside them (peaks that lower the likelihood dropped),
|
||||
each new copy refined beside the others held fixed, and kept where it stands out (below), it packs and
|
||||
R-free drops by 0.01 (`COPIES`, and `MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched
|
||||
side by side; with tNCS, every group that placed fewer copies than the best one is searched once more
|
||||
with the best one's copies' orientations tried first, and the better search kept.
|
||||
|
||||
The table has one row per group and merge. A row is `PLACED` where one of its copies' translation-
|
||||
function Z-scores (`TF_Z`, the highest) reaches 8, at most 5% of the C-alpha atoms clash with another
|
||||
The likelihood is the log-likelihood gain of the Rice/Woolfson distributions of E_obs given E_calc, with
|
||||
sigma_A falling off with the coordinate error the pLDDTs predict and its scale refined against the data,
|
||||
summed to 2.5 A (or the data's limit, or the 20000 lowest-resolution reflections in a large cell). A copy
|
||||
stands out where its translation-function Z-score reaches 8, or - only where the copies placed are as
|
||||
many as the cell most probably holds - where its likelihood after the rigid body reaches a Z-score of 8
|
||||
over the same copy in the same orientation at random positions (`LLG_Z`; taken only where `TF_Z` is
|
||||
below 8). That second way in is what places a single copy in a small cell, whose few reflections keep
|
||||
the translation function's own Z-score low.
|
||||
|
||||
Where the AlphaFold DB's predicted-aligned-error file sits beside the model under its own name
|
||||
(`AF-<id>-F1-model_v<n>.cif` and `AF-<id>-F1-predicted_aligned_error_v<n>.json`), the prediction is cut
|
||||
into rigid domains by average linkage on the PAE (domains a mean 10 A apart, of 40 residues or more;
|
||||
`AF_PAE_DOMAINS=`). Where the whole model is placed nowhere, each domain is searched by itself, side by
|
||||
side; the domains placed have to agree, and a domain alone can only confirm the answer's group
|
||||
(`AGREES`), never change it (`UNDECIDED` otherwise). `MODEL_PRIOR_SEARCH_MODEL=` names what the decision
|
||||
was made with.
|
||||
|
||||
The table has one row per group and merge. A row is `PLACED` where one of its copies stands out, at most 5% of the C-alpha atoms clash with another
|
||||
copy or a symmetry mate (`CLASH`), and one copy fits the cell (`MATTHEWS`: Vm below 1.5 A^3/Da). The
|
||||
null - the same model refitted from random orientations, each placed by the translation function as the
|
||||
search placed the real one (all copies turned together about their centroid where there are several) -
|
||||
@@ -995,9 +1011,12 @@ is built for the candidate the choice falls on; its R-work has to beat the null'
|
||||
(`NULL_SIGMA`: 3 sd allowing for the null being nine draws, as `MODEL_FIT_SIGMA`), and one that does not
|
||||
is marked `NULL` and the choice falls on the next. A placement clashing over twice the limit is not
|
||||
refined. Translational NCS is looked for in the native Patterson (a peak off the lattice of a fifth of
|
||||
the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's and
|
||||
whose search placed fewer copies leaves the choice `UNDECIDED`, since one copy beside a translated twin
|
||||
of itself fits the data no better than at random. The rotation and translation functions, the rigid
|
||||
the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's (or
|
||||
the chosen group's enantiomorph) and whose search placed fewer copies leaves the choice `UNDECIDED`,
|
||||
since one copy beside a translated twin of itself fits the data no better than at random - unless the
|
||||
chosen copies, carried into that group as they are and refined there, fit worse by 0.03 in R-free:
|
||||
where the translation does make up the difference between the groups, the copies build the same crystal
|
||||
in both. The rotation and translation functions, the rigid
|
||||
body and the structure factors the R-factors are taken on run on the
|
||||
GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` is then:
|
||||
|
||||
@@ -1012,10 +1031,16 @@ GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` i
|
||||
`ENANTIOMORPH_VETOED` and the label stays;
|
||||
- `SPACE_GROUP` - another group with the answer's rotations; written as `<prefix>_mr.mtz`/`.cif` from
|
||||
the answer's merge, the reflections that group extinguishes left out;
|
||||
- `TWIN_SUBGROUP` - no group of the answer's merge can hold the model at all (`MATTHEWS`/`CLASH` in
|
||||
every row) and a group of the twin hypothesis' subgroup can; written as `<prefix>_mr.mtz`/`.cif`
|
||||
from the twin hypothesis' merge, with its law and fraction. The twin hypothesis' merge is searched
|
||||
only in this case;
|
||||
- `TWIN_SUBGROUP` - the twin hypothesis' subgroup: either no group of the answer's merge can hold the
|
||||
model at all (`MATTHEWS`/`CLASH` in every row) and a group of the subgroup can, or the twin comparison
|
||||
chose it. The comparison runs where the answer's merge chose a group and the twin hypothesis' merge
|
||||
is there: on that merge, against its data as twinned by its law (the fraction refined on the working
|
||||
set), the chosen placement with all its symmetry images - the group's crystal read in the subgroup -
|
||||
is scored against the model searched for afresh in the subgroup; the subgroup is taken where its own
|
||||
placement leads by 0.03 in R-free and passes the subgroup's own choice (its null, its lead over the
|
||||
subgroup's other groups). Written as `<prefix>_mr.mtz`/`.cif` from the twin hypothesis' merge, with
|
||||
its law and fraction. With `--af-model`, the twin hypothesis' merge is made wherever a metric re-ask
|
||||
promoted the indexed lattice's group, so that the comparison can run;
|
||||
- `UNDECIDED` or `NOT_PLACED` - nothing is written.
|
||||
|
||||
With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `<prefix>_mr_model.pdb`/`.cif`.
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -23,6 +23,7 @@
|
||||
|
||||
#include <complex>
|
||||
#include <memory>
|
||||
#include <optional>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
@@ -43,10 +44,25 @@ struct SearchModel {
|
||||
int residues_kept = 0;
|
||||
double mass_da = 0; // of the residues kept, 110 Da each
|
||||
double radius_A = 0; // the farthest atom from the centroid
|
||||
double error_A = 0; // rms of the coordinate errors the pLDDTs predict, over the atoms kept
|
||||
std::vector<int> residue_index; // of each residue kept, in model order: its position in the prediction
|
||||
std::string description; // "the whole model", or which domain
|
||||
std::string failure_reason;
|
||||
};
|
||||
SearchModel PrepareSearchModel(const std::string &path);
|
||||
|
||||
// The rigid domains of a prediction from its predicted aligned error (PAE, the AlphaFold DB's JSON: a matrix
|
||||
// "predicted_aligned_error" over the residues of the whole prediction): the kept residues clustered by average
|
||||
// linkage on the symmetrised PAE until every two clusters are a mean PAE_DOMAIN_A apart; clusters of at least
|
||||
// DOMAIN_MIN_RESIDUES, largest first, each a list of positions in the prediction. Empty where the file cannot
|
||||
// be read or does not match the model.
|
||||
std::vector<std::vector<int>> PaeDomains(const std::string &pae_path, const SearchModel &model);
|
||||
// The search model of one domain: the residues of `model` at those positions, recentred.
|
||||
SearchModel DomainModel(const SearchModel &model, const std::vector<int> &residues);
|
||||
// The PAE file the AlphaFold DB puts beside a model: AF-<id>-F1-model_v<n>.cif ->
|
||||
// AF-<id>-F1-predicted_aligned_error_v<n>.json in the same directory; empty where the name does not follow it.
|
||||
std::string PaePathFor(const std::string &model_path);
|
||||
|
||||
// Rotations covering SO(3) with neighbours about `step_deg` apart: the identity, then shells of
|
||||
// rotation angle kappa = step, 2 step, ... 180 deg, each with its axes on a golden-spiral lattice of
|
||||
// the sphere as dense as 2 sin(kappa/2) - the distance a change of axis moves a rotation - needs.
|
||||
@@ -75,11 +91,13 @@ public:
|
||||
const gemmi::FPhiGrid<float> &Grid() const { return grid_; }
|
||||
double Box() const { return box_; }
|
||||
double Blur() const { return blur_; }
|
||||
double DMin() const { return d_min_; }
|
||||
|
||||
private:
|
||||
gemmi::FPhiGrid<float> grid_; // half-l, on a P1 box `box_` A wide
|
||||
double box_ = 0;
|
||||
double blur_ = 0;
|
||||
double d_min_ = 0;
|
||||
};
|
||||
|
||||
// The observed reflections the search scores against: E^2 in resolution bins, with each unique
|
||||
@@ -156,15 +174,62 @@ TranslationPeak TranslationFunction(const MolecularTransform &ft, const Translat
|
||||
const gemmi::Mat33 &rotation,
|
||||
const std::vector<std::complex<double>> *known = nullptr);
|
||||
|
||||
// The correlation of E^2 with |F|^2 of the known copies (KnownStructureFactors) and one more copy turned by
|
||||
// `rotation` with its centroid at `t`, over the terms of `setup`, each weighted by its multiplicity and each
|
||||
// bin's |F|^2 over its mean: the score the translation function's peaks for a further copy are ranked by.
|
||||
// Unlike the translation function, which is linear in |F|^2, it does not reward a copy that only repeats
|
||||
// the known ones - their own copy moved by one of the group's origin shifts doubles F on half the
|
||||
// reflections and cancels it on the other half, which the overlap with E^2 counts as a gain.
|
||||
double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell,
|
||||
const std::vector<std::complex<double>> &known, const gemmi::Mat33 &rotation,
|
||||
const gemmi::Fractional &t);
|
||||
// --- The likelihood a placement is judged by ---
|
||||
// The terms it is summed over: the reflections of a SearchReflections not absent in a group, each with
|
||||
// its E_obs, whether it is centric, its epsilon and 1/d^2, and for every operator j of the group h R_j (as
|
||||
// a Cartesian vector and as indices) and exp(2 pi i h.t_j), as in TranslationSetup.
|
||||
struct LikelihoodSetup {
|
||||
int ops = 0;
|
||||
std::vector<gemmi::Vec3> s; // terms x ops
|
||||
std::vector<gemmi::Miller> k; // terms x ops
|
||||
std::vector<std::complex<double>> phase; // terms x ops
|
||||
std::vector<double> eo, epsilon, s2; // terms
|
||||
std::vector<char> centric; // terms
|
||||
std::vector<int> bin; // terms
|
||||
int bins = 0;
|
||||
int max_k[3] = {0, 0, 0}; // the largest |index| of h R_j along each axis
|
||||
size_t Terms() const { return eo.size(); }
|
||||
};
|
||||
LikelihoodSetup MakeLikelihoodSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
|
||||
const gemmi::SpaceGroup &sg);
|
||||
|
||||
// g_j = G_j(h) exp(2 pi i h.t_j) of one copy in orientation `rotation`, terms x ops: the copy with its
|
||||
// centroid at t adds sum_j g_j exp(2 pi i (h R_j).t) to F(h). G_j(h) = F1(h R_j), F1 the transform of
|
||||
// the copy in `cell` with its centroid at the origin - one FFT of the cell's grid to `d_min`, read off at
|
||||
// the lattice points themselves, so exact at any resolution. Where `ft` is given and reaches `d_min`, read
|
||||
// off it instead, interpolated as the rotation and translation functions read it.
|
||||
std::vector<std::complex<double>> CopyFactors(const SearchModel &model, const LikelihoodSetup &setup,
|
||||
const gemmi::UnitCell &cell, const gemmi::Mat33 &rotation,
|
||||
double d_min, size_t nthreads, const MolecularTransform *ft = nullptr);
|
||||
// `f` += that copy's F with its centroid at the fractional position t.
|
||||
void AddCopy(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &g, const gemmi::Fractional &t,
|
||||
std::vector<std::complex<double>> &f);
|
||||
|
||||
// The log-likelihood gain of the model structure factors `f` (one per term) over a random structure,
|
||||
// with sigma_A = scale exp(-2 pi^2 error^2 / (3 d^2)): E_calc are |f| normalised in the resolution bins.
|
||||
double LogLikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &f, double error_A,
|
||||
double scale);
|
||||
// The same with the scale refined to its maximum (returned in `scale`).
|
||||
double RefinedLogLikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &f,
|
||||
double error_A, double *scale = nullptr);
|
||||
|
||||
// How the likelihood of one more copy in orientation `rotation` at `t`, beside the `known` structure factors
|
||||
// (zero for none), stands out over the same copy in the same orientation at random positions: the gain
|
||||
// with the scale refined at t, and its Z-score over the gains at those positions with that scale held.
|
||||
struct LikelihoodPeak {
|
||||
double llg = 0; // of all the copies, the known and this one
|
||||
double gain = 0; // over the known ones alone, at their own refined scale
|
||||
double z = 0;
|
||||
double scale = 0;
|
||||
};
|
||||
// `g`: CopyFactors of the copy in that orientation; `known_llg`: RefinedLogLikelihoodGain of `known` alone.
|
||||
LikelihoodPeak LikelihoodTranslationZ(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &known,
|
||||
double known_llg, const std::vector<std::complex<double>> &g,
|
||||
const gemmi::Fractional &t, double error_A, size_t nthreads);
|
||||
// The gain alone: of the known structure factors and the copy over the known ones (`llg`: of all of them).
|
||||
double LikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &known, double known_llg,
|
||||
const std::vector<std::complex<double>> &g, const gemmi::Fractional &t, double error_A,
|
||||
double *llg = nullptr, double *scale = nullptr);
|
||||
|
||||
// The model placed: rotated about its centroid by `rotation`, the centroid at `t`, in `cell` and `sg`.
|
||||
gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
|
||||
@@ -177,6 +242,28 @@ double ClashFraction(const gemmi::Structure &placed, size_t copies = 1);
|
||||
// Copies of the model as one structure: their chains in order, renamed so that no two share a name.
|
||||
gemmi::Structure AssembleCopies(const std::vector<gemmi::Structure> &copies);
|
||||
|
||||
// The copies of an assembly (AssembleCopies), its chains split evenly in order.
|
||||
std::vector<gemmi::Structure> SplitCopies(const gemmi::Structure &assembly, int copies);
|
||||
|
||||
// The setting of a subgroup's merge (cell `h`, group `h_sg`) in the cell `g` of a group `g_sg`: `basis` has h's
|
||||
// axes as columns, lattice vectors of g (with g's centring) in g's fractional coordinates, right-handed and
|
||||
// matching h's lengths and angles to 1%; `origin` is where g's origin lies in h's fractional coordinates, so
|
||||
// that the operators of g with h's rotations are h's own: x_h = basis^-1 x_g + origin. Empty where there is no
|
||||
// such setting.
|
||||
struct SubgroupSetting {
|
||||
gemmi::Mat33 basis;
|
||||
gemmi::Vec3 origin;
|
||||
};
|
||||
std::optional<SubgroupSetting> SubgroupBasis(const gemmi::UnitCell &g, const gemmi::SpaceGroup &g_sg,
|
||||
const gemmi::UnitCell &h, const gemmi::SpaceGroup &h_sg);
|
||||
|
||||
// The `copies` of a placement in group `g_sg` and every symmetry image of them that the subgroup `h_sg` does
|
||||
// not make itself, as independent copies in the subgroup's cell `h`, in its setting (SubgroupBasis); their
|
||||
// number in `copies_out`.
|
||||
gemmi::Structure ExpandIntoSubgroup(const gemmi::Structure &placed, int copies, const gemmi::SpaceGroup &g_sg,
|
||||
const SubgroupSetting &setting, const gemmi::UnitCell &h,
|
||||
const gemmi::SpaceGroup &h_sg, int *copies_out);
|
||||
|
||||
// The probabilities of 1, 2, ... copies of a model of `mass_da` in the asymmetric unit of a cell of
|
||||
// `asu_volume` A^3, from the distribution of the Matthews coefficient over protein crystals; one entry
|
||||
// per copy number whose Vm is at least 1.5 A^3/Da.
|
||||
@@ -209,6 +296,7 @@ struct MRMerge {
|
||||
const std::vector<MergedReflection> *merged = nullptr;
|
||||
UnitCell cell;
|
||||
gemmi::SpaceGroup space_group; // the group it was merged in
|
||||
std::optional<gemmi::Op> twin_law; // the twin hypothesis' merge: the law it carries
|
||||
};
|
||||
|
||||
// One space group's best placement: its first copy and every further one the search added.
|
||||
@@ -222,10 +310,18 @@ struct MRCandidate {
|
||||
double rf_z = 0; // of the first copy's orientation, over the rotation function
|
||||
double tf_z = 0; // the highest of the copies'
|
||||
std::vector<double> copy_tf_z; // each copy's, the later ones beside the earlier
|
||||
// The likelihood (LikelihoodTranslationZ), after the rigid body, taken only for a copy whose translation
|
||||
// function does not stand out (NaN otherwise): of all the copies with the last such one, and the highest
|
||||
// of their Z-scores, and each copy's.
|
||||
double llg = NAN;
|
||||
double llg_z = NAN;
|
||||
std::vector<double> copy_llg_z;
|
||||
double clash_fraction = 0; // of all the copies
|
||||
double r_work = NAN, r_free = NAN; // after the rigid body, to d_min_score
|
||||
double k_sol = NAN, b_sol = NAN; // the bulk solvent R was computed with
|
||||
double null_sigma = NAN; // R-work above its random-orientation null, in sd
|
||||
// On the twin hypothesis' merge, against the data as twinned by its law (the twin comparison).
|
||||
double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN;
|
||||
std::string verdict; // PLACED, or why not
|
||||
std::shared_ptr<gemmi::Structure> placed; // every copy
|
||||
};
|
||||
@@ -236,6 +332,8 @@ struct MolecularReplacementResult {
|
||||
std::string model_path;
|
||||
int residues_read = 0, residues_kept = 0;
|
||||
double mass_da = 0;
|
||||
std::string search_model; // the whole model, or the PAE domain the decision was made with
|
||||
int pae_domains = 0; // the rigid domains of the prediction's PAE; 0 without a PAE file
|
||||
double d_min_search = 0, d_max_search = 0;
|
||||
double d_min_score = 0; // the limit R-work / R-free are computed to
|
||||
std::vector<MRCandidate> candidates;
|
||||
@@ -261,6 +359,7 @@ struct MolecularReplacementResult {
|
||||
} hand_check;
|
||||
// Seconds per step.
|
||||
double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0;
|
||||
double t_domains = 0; // the PAE domains' own searches, wall clock (they run side by side)
|
||||
};
|
||||
|
||||
// Search `merges` (the answer's first) for the predicted model at `model_path`, in every space group with
|
||||
|
||||
@@ -5,6 +5,7 @@
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <mutex>
|
||||
#include <string>
|
||||
|
||||
#include <cufft.h>
|
||||
@@ -260,6 +261,25 @@ __global__ void translation_stats(const float *map, size_t n, double *sum, doubl
|
||||
}
|
||||
}
|
||||
|
||||
// cuFFT's r2c output of a map stored u fastest ([w][v][u/2+1], u halved) taken to MapToFPhi's half-l
|
||||
// layout (u fastest, w halved), conjugated and scaled as gemmi does.
|
||||
__global__ void half_u_to_half_w(const cufftComplex *in, int nu, int nv, int nw, float norm, float2 *out) {
|
||||
const int nwh = nw / 2 + 1, nuh = nu / 2 + 1;
|
||||
const size_t n = static_cast<size_t>(nu) * nv * nwh;
|
||||
const size_t i = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
|
||||
if (i >= n)
|
||||
return;
|
||||
const int u = static_cast<int>(i % nu), v = static_cast<int>((i / nu) % nv), l = static_cast<int>(i / (static_cast<size_t>(nu) * nv));
|
||||
cufftComplex x;
|
||||
if (u < nuh) {
|
||||
x = in[u + static_cast<size_t>(nuh) * (v + static_cast<size_t>(nv) * l)];
|
||||
} else { // F(-h) = conj F(h)
|
||||
x = in[(nu - u) + static_cast<size_t>(nuh) * ((nv - v) % nv + static_cast<size_t>(nv) * ((nw - l) % nw))];
|
||||
x.y = -x.y;
|
||||
}
|
||||
out[i] = make_float2(x.x * norm, -x.y * norm);
|
||||
}
|
||||
|
||||
Rot transposed(const gemmi::Mat33 &m) {
|
||||
Rot r;
|
||||
for (int i = 0; i < 3; i++)
|
||||
@@ -292,6 +312,51 @@ std::vector<int> bin_counts(const std::vector<int> &bin, int bins) {
|
||||
|
||||
} // namespace
|
||||
|
||||
bool MapToFPhiGPUFits(const gemmi::Grid<float> &map) {
|
||||
const size_t points = static_cast<size_t>(map.nu) * map.nv * map.nw;
|
||||
const size_t half_u = static_cast<size_t>(map.nu / 2 + 1) * map.nv * map.nw;
|
||||
const size_t half_w = static_cast<size_t>(map.nu) * map.nv * (map.nw / 2 + 1);
|
||||
const size_t bytes = points * sizeof(float) + (half_u + half_w) * sizeof(cufftComplex);
|
||||
size_t free_bytes = 0, total = 0;
|
||||
cuda_err(cudaMemGetInfo(&free_bytes, &total));
|
||||
return bytes <= total / 4;
|
||||
}
|
||||
|
||||
gemmi::FPhiGrid<float> MapToFPhiGPU(const gemmi::Grid<float> &map) {
|
||||
// One at a time: concurrent searches (the PAE domains') would each hold a large grid on the card.
|
||||
static std::mutex m;
|
||||
std::lock_guard lock(m);
|
||||
gemmi::FPhiGrid<float> hkl;
|
||||
hkl.unit_cell = map.unit_cell;
|
||||
hkl.spacegroup = map.spacegroup;
|
||||
hkl.axis_order = map.axis_order;
|
||||
hkl.half_l = true;
|
||||
hkl.set_size_without_checking(map.nu, map.nv, map.nw / 2 + 1);
|
||||
const float norm = static_cast<float>(map.unit_cell.volume / map.point_count());
|
||||
CudaStream stream;
|
||||
auto d_map = upload(map.data, stream);
|
||||
CudaDevicePtr<cufftComplex> d_half(static_cast<size_t>(map.nu / 2 + 1) * map.nv * map.nw);
|
||||
CudaDevicePtr<float2> d_out(hkl.data.size());
|
||||
cufftHandle plan;
|
||||
cufft_err(cufftPlan3d(&plan, map.nw, map.nv, map.nu, CUFFT_R2C));
|
||||
try {
|
||||
cufft_err(cufftSetStream(plan, stream));
|
||||
cufft_err(cufftExecR2C(plan, d_map.get(), d_half.get()));
|
||||
half_u_to_half_w<<<static_cast<unsigned>((hkl.data.size() + THREADS - 1) / THREADS), THREADS, 0, stream>>>(
|
||||
d_half, map.nu, map.nv, map.nw, norm, d_out);
|
||||
cuda_err(cudaGetLastError());
|
||||
static_assert(sizeof(std::complex<float>) == sizeof(float2));
|
||||
cuda_err(cudaMemcpyAsync(hkl.data.data(), d_out.get(), hkl.data.size() * sizeof(float2),
|
||||
cudaMemcpyDeviceToHost, stream));
|
||||
cuda_err(cudaStreamSynchronize(stream));
|
||||
} catch (...) {
|
||||
cufftDestroy(plan);
|
||||
throw;
|
||||
}
|
||||
cufftDestroy(plan);
|
||||
return hkl;
|
||||
}
|
||||
|
||||
struct MolecularReplacementGPUState {
|
||||
CudaStream stream;
|
||||
CudaDevicePtr<float2> grid;
|
||||
|
||||
@@ -15,6 +15,13 @@
|
||||
|
||||
#include "MolecularReplacement.h"
|
||||
|
||||
// MapToFPhi (ModelFFT.h) with cuFFT on the calling thread's current device: the same transform, the same
|
||||
// half-l layout and conventions, to float rounding. For the molecular transform of a large model, whose box is
|
||||
// hundreds of millions of points. MapToFPhiGPUFits: whether its device buffers take at most a quarter of the
|
||||
// card's total memory - decided on the total, so the same input always takes the same path.
|
||||
bool MapToFPhiGPUFits(const gemmi::Grid<float> &map);
|
||||
gemmi::FPhiGrid<float> MapToFPhiGPU(const gemmi::Grid<float> &map);
|
||||
|
||||
struct MolecularReplacementGPUState;
|
||||
struct MolecularReplacementGPUGroup;
|
||||
|
||||
|
||||
+28
-15
@@ -2587,6 +2587,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
Add(s, KeyEnum("MODEL_PRIOR", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
|
||||
Add(s, KeyInt("AF_RESIDUES_READ", mp.residues_read));
|
||||
Add(s, KeyInt("AF_RESIDUES_KEPT", mp.residues_kept));
|
||||
Add(s, KeyInt("AF_PAE_DOMAINS", mp.pae_domains));
|
||||
Add(s, KeyText("MODEL_PRIOR_SEARCH_MODEL", mp.search_model));
|
||||
Add(s, KeyReal("SEARCH_D_MIN", mp.d_min_search, "{:.2f}"));
|
||||
Add(s, KeyReal("SCORE_D_MIN", mp.d_min_score, "{:.2f}"));
|
||||
Add(s, KeyText("MODEL_PRIOR_DECISION", mp.decision));
|
||||
@@ -2612,22 +2614,23 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
Add(s, KeyText("MODEL_PRIOR_FILE_PREFIX", result.model_prior_file_prefix));
|
||||
Add(s, KeyEnum("MODEL_PRIOR_SEARCH_DEVICE", mp.gpu_search ? "GPU" : "CPU", {"GPU", "CPU"}));
|
||||
Add(s, KeyReal("MODEL_PRIOR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation
|
||||
+ mp.t_rescore + mp.t_null + mp.hand_check.seconds, "{:.2f}"));
|
||||
+ mp.t_rescore + mp.t_null + mp.t_domains + mp.hand_check.seconds, "{:.2f}"));
|
||||
Add(s, Blank());
|
||||
ReportEntry e;
|
||||
e.kind = ReportEntry::Kind::Table;
|
||||
e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "CLASH",
|
||||
"R_WORK", "R_FREE", "NULL_SIGMA", "VERDICT"};
|
||||
e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z CLASH R_WORK R_FREE NULL_SIGMA VERDICT\n"
|
||||
" ------------------- ----------- ----- ------ ------ ----- ----- ----- ------ ------ ---------- ----------";
|
||||
e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "LLG_Z", "CLASH",
|
||||
"R_WORK", "R_FREE", "R_FREE_TWIN", "NULL_SIGMA", "VERDICT"};
|
||||
e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z LLG_Z CLASH R_WORK R_FREE R_FREE_TWIN NULL_SIGMA VERDICT\n"
|
||||
" ------------------- ----------- ----- ------ ------ ----- ----- ----- ----- ------ ------ ----------- ---------- ----------";
|
||||
for (const auto &c : mp.candidates) {
|
||||
const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE";
|
||||
const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); };
|
||||
e.table.text_rows.push_back(fmt::format(
|
||||
" {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}",
|
||||
" {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:>5s} {:4.1f}% {:>6s} {:>6s} {:>11s} {:>10s} {}",
|
||||
merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z,
|
||||
std::isfinite(c.llg_z) ? fmt::format("{:.1f}", c.llg_z) : std::string("-"),
|
||||
100 * c.clash_fraction, r(c.r_work),
|
||||
r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"),
|
||||
r(c.r_free), r(c.r_free_twin), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"),
|
||||
c.verdict));
|
||||
std::vector<ReportValue> row;
|
||||
row.push_back(KeyText("", merge).value);
|
||||
@@ -2637,9 +2640,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
row.push_back(KeyInt("", c.copies_likely).value);
|
||||
row.push_back(KeyReal("", c.rf_z, "{:.1f}").value);
|
||||
row.push_back(KeyReal("", c.tf_z, "{:.1f}").value);
|
||||
row.push_back(KeyReal("", c.llg_z, "{:.1f}").value);
|
||||
row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value);
|
||||
row.push_back(KeyReal("", c.r_work, "{:.3f}").value);
|
||||
row.push_back(KeyReal("", c.r_free, "{:.3f}").value);
|
||||
row.push_back(KeyReal("", c.r_free_twin, "{:.3f}").value);
|
||||
row.push_back(KeyReal("", c.null_sigma, "{:+.1f}").value);
|
||||
row.push_back(KeyText("", c.verdict).value);
|
||||
e.table.cells.push_back(std::move(row));
|
||||
@@ -2651,16 +2656,24 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
" cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n"
|
||||
" and the best placements taken through the rigid body; then further copies, one at a time, up to\n"
|
||||
" the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n"
|
||||
" each with the translation function beside the copies already placed, kept where its TF_Z is 8\n"
|
||||
" or more, it packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES, to\n"
|
||||
" SCORE_D_MIN; TF_Z is the highest of the copies'. PLACED needs TF_Z >= 8, at most 5% of the\n"
|
||||
" each with the translation function beside the copies already placed, its peaks ranked by the\n"
|
||||
" likelihood gain of the copy (sigma_A from the pLDDT's coordinate error, its scale refined), kept\n"
|
||||
" where the copy stands out, packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES,\n"
|
||||
" to SCORE_D_MIN; TF_Z is the highest of the copies'; LLG_Z the highest Z-score of a copy's refined\n"
|
||||
" likelihood gain over the same copy at random positions, where its TF_Z was below 8. A copy stands\n"
|
||||
" out where TF_Z or LLG_Z is 8 or more. PLACED needs a copy that stands out, at most 5% of the\n"
|
||||
" C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n"
|
||||
" (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n"
|
||||
" better than the same copies refitted from random orientations (NULL where it is not; the next\n"
|
||||
" candidate is then taken). The model chooses within a merge where its placed group leads every\n"
|
||||
" other that packs by 0.03 in R-free, and the twin hypothesis' merge - searched only then - where\n"
|
||||
" no group of the answer's merge can hold it at all. A group one copy cannot fit is checked for\n"
|
||||
" packing only, not refined.\n"
|
||||
" other that packs by 0.03 in R-free, and the twin hypothesis' merge where no group of the\n"
|
||||
" answer's merge can hold it at all, or where, on that merge and against its data as twinned by\n"
|
||||
" its law, the model searched for in the subgroup leads the chosen placement read as a crystal\n"
|
||||
" there by 0.03 in R-free (R_FREE_TWIN). A group one copy cannot fit is checked for\n"
|
||||
" packing only, not refined. Where the AlphaFold DB's predicted aligned error sits beside the\n"
|
||||
" model, the prediction is cut into rigid domains (AF_PAE_DOMAINS), residues in none left out;\n"
|
||||
" where the whole model is placed nowhere, each domain is searched by itself, and the domains\n"
|
||||
" placed have to agree (MODEL_PRIOR_SEARCH_MODEL: what the decision was made with).\n"
|
||||
" The answer's files take the enantiomorph from it - only where model validation of the placed\n"
|
||||
" model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n"
|
||||
" atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n"
|
||||
@@ -2669,8 +2682,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
result.model_prior_file_prefix.empty() ? std::string("<prefix>_mr") : result.model_prior_file_prefix)));
|
||||
Add(s, Prose(fmt::format(
|
||||
" Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n"
|
||||
" R-free {:.2f}, null {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, mp.t_rotation,
|
||||
mp.t_translation, mp.t_rescore, mp.t_null, mp.hand_check.seconds), true));
|
||||
" R-free {:.2f}, null {:.2f}, PAE domains {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform,
|
||||
mp.t_rotation, mp.t_translation, mp.t_rescore, mp.t_null, mp.t_domains, mp.hand_check.seconds), true));
|
||||
}
|
||||
doc.sections.push_back(std::move(s));
|
||||
}
|
||||
|
||||
@@ -2874,16 +2874,22 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
// too, with the law - so a refinement program can take the twin hypothesis without a re-run.
|
||||
// The answer itself is not touched by any of this.
|
||||
const bool pre_twin_like = result.pre_promotion_twinning && result.pre_promotion_twinning->twinning_suspected;
|
||||
// A predicted model can test the twin hypothesis itself (--af-model's twin comparison), so with one the
|
||||
// merge in the indexed lattice's group is made wherever a metric re-ask promoted from it.
|
||||
const bool model_tests_twin = !config_.af_model_path.empty() && twin_metric_subgroup.has_value();
|
||||
if (!geometry_prepass && result.space_group && sg_search.point_group_order > 1
|
||||
&& (twin_overruled_subgroup || pre_twin_like)) {
|
||||
&& (twin_overruled_subgroup || pre_twin_like || model_tests_twin)) {
|
||||
twin_choice = ChooseTwinSubgroup(*result.space_group, sg_search,
|
||||
twin_overruled_subgroup ? twin_overruled_subgroup
|
||||
: twin_metric_subgroup);
|
||||
twin_choice_group_number = result.space_group->number;
|
||||
twin_choice_reason = twin_overruled_subgroup
|
||||
? fmt::format("a twin gate refused it ({}) and was overruled", twin_overruled_reason)
|
||||
: fmt::format("the merge read twin-like before the space-group search ({})",
|
||||
TwinningVerdictLine(*result.pre_promotion_twinning));
|
||||
: pre_twin_like
|
||||
? fmt::format("the merge read twin-like before the space-group search ({})",
|
||||
TwinningVerdictLine(*result.pre_promotion_twinning))
|
||||
: std::string("it was promoted from the indexed lattice's group, which the predicted model "
|
||||
"can test");
|
||||
if (twin_choice)
|
||||
logger.Info("Twin hypothesis: {} was adopted, but {} - the merge in {} is written too, "
|
||||
"with twin law {}", result.space_group->xhm(), twin_choice_reason,
|
||||
@@ -3992,7 +3998,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
TimingMark mark(logger, "model prior");
|
||||
std::vector<MRMerge> merges{{"MERGE", &sm.merged, *result.consensus_cell, experiment_.GetSpaceGroupOrP1()}};
|
||||
if (twin_file)
|
||||
merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group});
|
||||
merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group,
|
||||
twin_file->law});
|
||||
MolecularReplacementResult mr = RunMolecularReplacement(config_.af_model_path, merges, logger,
|
||||
static_cast<size_t>(config_.nthreads));
|
||||
if (!mr.ok)
|
||||
|
||||
@@ -202,7 +202,7 @@ void print_usage() {
|
||||
std::cout << " -z, --reference <file> Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl;
|
||||
std::cout << " --reference-column <label> Reference column to use, by MTZ label; an mmCIF is read under the labels gemmi cif2mtz gives it (IMEAN, FP, I(+), FC, F-model, ...) (default: auto - F-model, else IMEAN/I, else FP/FOBS/F)" << std::endl;
|
||||
std::cout << " --model <file> After merging, validate vs this atomic model: R-free, sigma_A-weighted 2mFo-DFc/mFo-DFc maps, the model as placed (written beside them), and - where the merge kept the Bijvoet split - an anomalous difference map with the strongest anomalous scatterers named by the atom they sit on. PDB or mmCIF, gzipped or not, recognised by content rather than by extension" << std::endl;
|
||||
std::cout << " --af-model <file> Experimental: place this predicted model (AlphaFold mmCIF, pLDDT in the B column; residues below pLDDT 70 removed) by molecular replacement, one copy in the asymmetric unit, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. Reported in the MODEL PRIOR section; it may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as <prefix>_mr.mtz/.cif with the placed model <prefix>_mr_model.pdb. --mode mx only" << std::endl;
|
||||
std::cout << " --af-model <file> Experimental: place this predicted model (AlphaFold mmCIF, pLDDT in the B column; residues below pLDDT 70 removed; where the AlphaFold DB's predicted-aligned-error JSON sits beside it under its own name, the model is cut into rigid domains) by molecular replacement, as many copies as the cell holds, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. Reported in the MODEL PRIOR section; it may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as <prefix>_mr.mtz/.cif with the placed model <prefix>_mr_model.pdb. --mode mx only" << std::endl;
|
||||
std::cout << " --export-unmerged Write <prefix>_unmerged.mtz, an unmerged MTZ (POINTLESS column layout) of the integrated observations, for aimless / pointless / careless. On by default, whenever there is an output prefix. On a rotation run the partials of each reflection are summed into one full, written at the batch its rocking curve is centred on, with the summed rocking-curve fraction in FRACTIONCALC; an event that caught less of its rocking curve than --min-captured-fraction (or --min-partiality) is not written, exactly as in the merge, since a full is what the reading program is told these rows are. Intensities carry the deterministic per-reflection corrections and nothing else - Lorentz and polarization in the LP column, the sensor's angle-dependent efficiency in QE and the flight path's attenuation in FLIGHT, so raw counts are I / LP * QE * FLIGHT - while the partiality is not divided out and the per-image scale is not applied at all, since those programs scale the data themselves. Written in --mode mx and --mode scale, and with --no-merge" << std::endl;
|
||||
std::cout << " --no-export-unmerged Do not write <prefix>_unmerged.mtz. It is the largest file a run produces, so a run whose observations are not going to another scaling program - a regression battery, say - can skip writing it" << std::endl;
|
||||
std::cout << " --export-unmerged-partials Write <prefix>_unmerged_partials.mtz, the same observations with each partial as its own row (one batch per image) flagged for the reading program to sum, instead of summed here. Off by default and independent of --export-unmerged; both can be written in one run" << std::endl;
|
||||
|
||||
@@ -4,6 +4,8 @@
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <fstream>
|
||||
#include <random>
|
||||
|
||||
#include <gemmi/dencalc.hpp>
|
||||
@@ -202,6 +204,120 @@ TEST_CASE("MolecularReplacement_SecondCopyBesideKnownFirst", "[MolecularReplacem
|
||||
PlaceSearchModel(model, r2, t1, cell, sg)}), 2) > 0.5);
|
||||
}
|
||||
|
||||
// The likelihood of the true placement stands out over the same orientation elsewhere, and the gain is zero
|
||||
// where sigma_A is.
|
||||
TEST_CASE("MolecularReplacement_LikelihoodPlacement", "[MolecularReplacement]") {
|
||||
const SearchModel model = synthetic_model();
|
||||
const gemmi::UnitCell cell(42, 50, 58, 90, 90, 90);
|
||||
const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
|
||||
const gemmi::Mat33 truth = axis_angle({1, 2, 3}, 50);
|
||||
const gemmi::Fractional t(0.17, 0.31, 0.09);
|
||||
const auto merged = observed(PlaceSearchModel(model, truth, t, cell, sg));
|
||||
const LikelihoodSetup setup = MakeLikelihoodSetup(MakeSearchReflections(merged, cell, sg, 3.5, 12.0), cell, sg);
|
||||
REQUIRE(setup.Terms() > 500);
|
||||
const std::vector<std::complex<double>> none(setup.Terms(), 0.0);
|
||||
const auto g = CopyFactors(model, setup, cell, truth, 3.5, 4);
|
||||
|
||||
const LikelihoodPeak peak = LikelihoodTranslationZ(setup, none, 0.0, g, t, 0.3, 4);
|
||||
CHECK(peak.llg > 100.0);
|
||||
CHECK(peak.z > 8.0);
|
||||
CHECK(peak.scale > 0.8);
|
||||
const double wrong = LikelihoodGain(setup, none, 0.0, g, gemmi::Fractional(0.4, 0.05, 0.3), 0.3);
|
||||
CHECK(wrong < peak.llg / 4);
|
||||
std::vector<std::complex<double>> f = none;
|
||||
AddCopy(setup, g, t, f);
|
||||
CHECK(LogLikelihoodGain(setup, f, 0.3, 0.0) == Catch::Approx(0.0).margin(1e-9));
|
||||
}
|
||||
|
||||
TEST_CASE("MolecularReplacement_PaeDomains", "[MolecularReplacement]") {
|
||||
CHECK(PaePathFor("/x/AF-P12345-F1-model_v6.cif") == "/x/AF-P12345-F1-predicted_aligned_error_v6.json");
|
||||
CHECK(PaePathFor("AF-P12345-F1-model_v4.pdb.gz") == "AF-P12345-F1-predicted_aligned_error_v4.json");
|
||||
CHECK(PaePathFor("/x/model.cif").empty());
|
||||
|
||||
// A prediction of 120 residues: two confident blocks, 1-50 and 61-120, placed relative to each other
|
||||
// with no confidence, and a loop 51-60 the trim removed.
|
||||
SearchModel m;
|
||||
gemmi::Model model(1);
|
||||
gemmi::Chain chain("A");
|
||||
for (int i = 0; i < 120; i++) {
|
||||
if (i >= 50 && i < 60)
|
||||
continue;
|
||||
gemmi::Residue r;
|
||||
r.name = "ALA";
|
||||
r.seqid = gemmi::SeqId(i + 1, ' ');
|
||||
gemmi::Atom a;
|
||||
a.name = "CA";
|
||||
a.element = gemmi::El::C;
|
||||
a.pos = gemmi::Position(3.8 * i, i < 50 ? 0.0 : 20.0, 0.0);
|
||||
a.b_iso = 10.0f;
|
||||
a.occ = 1.0f;
|
||||
r.atoms.push_back(a);
|
||||
chain.residues.push_back(r);
|
||||
m.residue_index.push_back(i);
|
||||
}
|
||||
model.chains.push_back(chain);
|
||||
m.st.models.push_back(model);
|
||||
m.residues_read = 120;
|
||||
m.residues_kept = 110;
|
||||
const std::string path = "MolecularReplacement_PaeDomains.json";
|
||||
{
|
||||
std::ofstream out(path);
|
||||
out << "[{\"predicted_aligned_error\":[";
|
||||
for (int i = 0; i < 120; i++) {
|
||||
out << (i ? "," : "") << "[";
|
||||
for (int j = 0; j < 120; j++)
|
||||
out << (j ? "," : "") << ((i < 55) == (j < 55) ? 2 : 25);
|
||||
out << "]";
|
||||
}
|
||||
out << "],\"max_predicted_aligned_error\":31.75}]";
|
||||
}
|
||||
const auto domains = PaeDomains(path, m);
|
||||
REQUIRE(domains.size() == 2);
|
||||
CHECK(domains[0].size() == 60);
|
||||
CHECK(domains[0].front() == 60);
|
||||
CHECK(domains[1].size() == 50);
|
||||
const SearchModel d = DomainModel(m, domains[1]);
|
||||
CHECK(d.residues_kept == 50);
|
||||
CHECK(d.description == "PAE domain of 50 residues (1-50)");
|
||||
std::remove(path.c_str());
|
||||
}
|
||||
|
||||
// A placement of a group carried into a subgroup's merge: in P 21 21 21 read as P 1 21 1 on the same cell each
|
||||
// copy becomes two, the second its image under a two-fold the subgroup does not have; in C 2 2 21 read as
|
||||
// P 1 21 1 on the primitive monoclinic cell half the volume holds the same two.
|
||||
TEST_CASE("MolecularReplacement_ExpandIntoSubgroup", "[MolecularReplacement]") {
|
||||
const SearchModel model = synthetic_model();
|
||||
{
|
||||
const gemmi::UnitCell cell(42, 50, 58, 90, 90, 90);
|
||||
const gemmi::SpaceGroup &g = gemmi::get_spacegroup_by_name("P 21 21 21");
|
||||
const gemmi::SpaceGroup &h = gemmi::get_spacegroup_by_name("P 1 21 1");
|
||||
const auto setting = SubgroupBasis(cell, g, cell, h);
|
||||
REQUIRE(setting);
|
||||
const gemmi::Structure placed = PlaceSearchModel(model, axis_angle({1, 2, 3}, 50), gemmi::Fractional(0.17, 0.31, 0.09), cell, g);
|
||||
int n = 0;
|
||||
const gemmi::Structure expanded = ExpandIntoSubgroup(placed, 1, g, *setting, cell, h, &n);
|
||||
CHECK(n == 2);
|
||||
// The same crystal: the copies and their P 1 21 1 mates are the P 21 21 21 crystal, which packs.
|
||||
CHECK(ClashFraction(expanded, 2) < 0.05);
|
||||
}
|
||||
{
|
||||
const gemmi::UnitCell g_cell(56, 158, 144, 90, 90, 90);
|
||||
const gemmi::SpaceGroup &g = gemmi::get_spacegroup_by_name("C 2 2 21");
|
||||
const double half = std::sqrt(56.0 * 56.0 + 158.0 * 158.0) / 2;
|
||||
const double beta = 180.0 - std::acos(28.0 / half) * 180.0 / M_PI;
|
||||
const gemmi::UnitCell h_cell(half, 144, 56, 90, beta, 90);
|
||||
const gemmi::SpaceGroup &h = gemmi::get_spacegroup_by_name("P 1 21 1");
|
||||
const auto setting = SubgroupBasis(g_cell, g, h_cell, h);
|
||||
REQUIRE(setting);
|
||||
CHECK(std::fabs(setting->basis.determinant() - 0.5) < 1e-9);
|
||||
const gemmi::Structure placed = PlaceSearchModel(model, axis_angle({1, 2, 3}, 50), gemmi::Fractional(0.17, 0.31, 0.09), g_cell, g);
|
||||
int n = 0;
|
||||
const gemmi::Structure expanded = ExpandIntoSubgroup(placed, 1, g, *setting, h_cell, h, &n);
|
||||
CHECK(n == 2);
|
||||
CHECK(ClashFraction(expanded, 2) < 0.05);
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("MolecularReplacement_MatthewsProbabilities", "[MolecularReplacement]") {
|
||||
// Vm 2.5 A^3/Da for one copy: one copy is by far the most probable, four impossible (Vm 0.6).
|
||||
const auto p = MatthewsProbabilities(2.5 * 20000, 20000);
|
||||
@@ -312,12 +428,28 @@ TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") {
|
||||
for (size_t i = 0; i < orientations.size(); i++) {
|
||||
const TranslationPeak cpu = TranslationFunction(ft, setup, orientations[i]);
|
||||
CHECK(std::fabs(tf_gpu[i].z - cpu.z) < 1e-3);
|
||||
if (i == 0) { // the solution's peak is unique; noise peaks can tie between grid points
|
||||
CHECK(tf_gpu[i].t.x == Catch::Approx(cpu.t.x));
|
||||
CHECK(tf_gpu[i].t.y == Catch::Approx(cpu.t.y));
|
||||
CHECK(tf_gpu[i].t.z == Catch::Approx(cpu.t.z));
|
||||
}
|
||||
// The solution's peak is the same up to the origins P 21 21 21 allows, which tie exactly; noise peaks
|
||||
// can tie between grid points.
|
||||
if (i == 0)
|
||||
for (double d : {tf_gpu[i].t.x - cpu.t.x, tf_gpu[i].t.y - cpu.t.y, tf_gpu[i].t.z - cpu.t.z})
|
||||
CHECK(std::fabs(2 * d - std::round(2 * d)) < 1e-6);
|
||||
}
|
||||
// The molecular transform's FFT.
|
||||
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
|
||||
dc.d_min = 4.0;
|
||||
dc.rate = 1.5;
|
||||
gemmi::Structure st = PlaceSearchModel(model, truth, t, cell, gemmi::get_spacegroup_p1());
|
||||
dc.set_grid_cell_and_spacegroup(st);
|
||||
dc.put_model_density_on_grid(st.models[0]);
|
||||
REQUIRE(MapToFPhiGPUFits(dc.grid));
|
||||
const auto f_cpu = MapToFPhi(dc.grid), f_gpu = MapToFPhiGPU(dc.grid);
|
||||
REQUIRE(f_gpu.data.size() == f_cpu.data.size());
|
||||
double worst = 0, largest = 0;
|
||||
for (size_t i = 0; i < f_cpu.data.size(); i++) {
|
||||
worst = std::max(worst, static_cast<double>(std::abs(f_gpu.data[i] - f_cpu.data[i])));
|
||||
largest = std::max(largest, static_cast<double>(std::abs(f_cpu.data[i])));
|
||||
}
|
||||
CHECK(worst < 1e-4 * largest);
|
||||
// Beside a known copy.
|
||||
const TranslationSetup setup2 = MakeTranslationSetup(refl, cell, sg, true);
|
||||
const auto known = KnownStructureFactors(ft, setup2, cell, {rotations[5]}, {gemmi::Fractional(0.3, 0.1, 0.4)});
|
||||
|
||||
Reference in New Issue
Block a user