diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 23e1b4600..971f91405 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -459,11 +459,14 @@ gemmi::Fractional centroid(const gemmi::Structure &st) { // about its centroid, in standard deviations of the latter: the null model validation builds for a // model that claims something (ModelValidation.cpp). An orientation the crystal cannot tell from the // placement's - within the rigid body's reach of one of its symmetry classes - is drawn again. +// With several copies each is turned about its own centroid, all by the same rotation: the copies keep +// their arrangement, and with it the fit that tNCS gives any model placed parallel at its vector, so that +// only the model's own orientation is random. // `translate`: where given, each replicate is placed where the translation function puts the model in its // orientation (the centroids, fractional, of the model in each of the orientations given), as the search // placed the real one, instead of at the real placement's centroid; the search chose that position among // all others for the real orientation, and a null that is not given the same choice reads high. -double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg, +double null_sigma(const gemmi::Structure &placed, int copies, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg, const std::vector &merged, double d_min, RigidBodyGPUPool *pool, size_t nthreads, ScoreEngines *scorers, const std::function(const std::vector &)> &translate) { @@ -497,15 +500,23 @@ double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientatio workers.push_back(std::async(std::launch::async, [&] { for (int i = next++; i < NULL_REPLICATES; i = next++) { gemmi::Structure rep = placed; - std::vector pos = ModelPositions(rep.models[0]); - gemmi::Vec3 c; - for (const gemmi::Position &p : pos) - c += p; - c *= 1.0 / static_cast(pos.size()); - const gemmi::Vec3 to = translate ? gemmi::Vec3(placed.cell.orthogonalize(where[i])) : c; - for (gemmi::Position &p : pos) - p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + to); - SetModelPositions(rep.models[0], pos); + std::vector &chains = rep.models[0].chains; + for (int k = 0; k < copies; k++) { + // Copy k: the chains SplitCopies gives it. + std::vector atoms; + for (size_t ch = 0; ch < chains.size(); ch++) + if (static_cast(ch * copies / chains.size()) == k) + for (gemmi::Residue &r : chains[ch].residues) + for (gemmi::Atom &a : r.atoms) + atoms.push_back(&a); + gemmi::Vec3 c; + for (const gemmi::Atom *a : atoms) + c += a->pos; + c *= 1.0 / static_cast(atoms.size()); + const gemmi::Vec3 to = translate ? gemmi::Vec3(placed.cell.orthogonalize(where[i])) : c; + for (gemmi::Atom *a : atoms) + a->pos = gemmi::Position(draws[i].multiply(gemmi::Vec3(a->pos) - c) + to); + } r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads, scorers).r_work; } })); @@ -2440,7 +2451,7 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec out.push_back(p.t); return out; }; - c.null_sigma = null_sigma(*c.placed, orientations[best], c.space_group, *merges[merge].merged, + c.null_sigma = null_sigma(*c.placed, c.copies, orientations[best], c.space_group, *merges[merge].merged, result.d_min_score, engines.get(), nthreads, &scorers, translate); result.t_null += seconds_since(t_null); if (c.null_sigma >= NULL_SIGMA_MIN)