Jungfraujoch/image_analysis/structure_refinement/SolventBuild.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <map>
#include <set>
#include <string>
#include <vector>

#include "gemmi/model.hpp"

#include "SolventSites.h"   // SolventSite, SolventMaps

// One round of the solvent build that RefineModel runs within its weight trajectory (Refinement.h): waters placed at the
// peaks of the difference map a refined model leaves, the model's waters the maps no longer support taken out, and
// waters and ions relabelled.
//
// A peak of mFo-DFc becomes a water where it has a hydrogen-bond partner and nothing too close - on a special
// position exactly on it, at the occupancy that makes its coinciding copies one water; a peak that
// cannot be a single solvent atom - one bonded to the model, beside a negative peak on a model atom, where a
// truncated side chain's missing atoms would be, one of a cluster of peaks at bonded distances, elongated, or
// denser than the model's S and P - is reported as not solvent and left alone. A water or a single-atom ion whose
// surroundings and density support another species, not ambiguously (SolventSites.h), is relabelled as it: a water
// becomes an ion, an ion water or another ion.
struct SolventEdit {
    enum class Kind { Added, Removed, NotSolvent, Relabelled };
    Kind kind = Kind::NotSolvent;
    int round = 0;
    std::string site;                // the water as modelled ("HOH A401"); empty for a peak not taken
    std::string chain;               // where an added water goes: the chain of the atom nearest to it
    gemmi::SeqId seqid;              // a removed water's residue
    char altloc = '\0';
    gemmi::Position pos;
    double difference_sigma = 0;     // mFo-DFc at the site, in its rms
    double map_sigma = 0;            // 2mFo-DFc
    double b_iso = 0;                // an added water's starting B, a removed one's
    double occupancy = 1;            // an added water's: 1 over its site's multiplicity on a special position
    double elongation = 0;           // a peak's density: the longest axis over the shortest, rms
    std::string species;             // a relabelled site's new residue name (HOH, MG, NA, K, CA, CL); its old one in `site`
    std::string reason;              // why: the partners an added water has, what a removed or refused one lacks, the
                                     // scores a relabelled one has
};

const char *SolventEditKindName(SolventEdit::Kind kind);

struct SolventRound {
    std::vector<SolventEdit> edits;
    int added = 0, removed = 0, not_solvent = 0, relabelled = 0;
    std::map<std::string, int> refused;   // the peaks neither taken nor reported, by what they lack
};

// The resolution waters are built to at most, following the ordered-solvent defaults of phenix.refine
// (Afonine et al. (2012) Acta Cryst. D68, 352-367): no picking where the data stop short of 2.8 A.
constexpr double SOLVENT_BUILD_D_MIN = 2.8;

// At most this many rounds, each judged by refining the solvent it changed (a third one gained at most 0.0004 in
// R-free where it was kept, on twelve sets, and was refused on the others).
constexpr int SOLVENT_BUILD_ROUNDS = 2;

// What the build did, over the rounds kept.
struct SolventBuildResult {
    bool ran = false;
    int rounds = 0;                  // the rounds kept
    int added = 0;                   // waters added and still in the model
    int removed = 0;                 // waters of the model as given taken out
    int relabelled = 0;              // sites relabelled, in the rounds kept
    int not_solvent = 0;             // peaks of the final model's difference map that are no solvent
    std::vector<SolventEdit> edits;  // the kept rounds' additions and removals, then the final not-solvent peaks
    std::string stop;                // why it stopped
    double seconds = 0;
};

// A round on `model` (the refined model the maps were made from, in their cell and space group): `maps` its
// 2mFo-DFc and mFo-DFc (the anomalous one unused), `sites` its solvent audit (AnalyseSolventSites on the same
// maps), `built` the waters added by earlier rounds (by name, "HOH A401"), `b_ceiling` the highest B its refinement
// allowed; `relabelled` the sites relabelled by earlier rounds (by their new name), which are not relabelled again.
SolventRound PickSolvent(const gemmi::Structure &model, const SolventMaps &maps, const std::vector<SolventSite> &sites,
                         const std::set<std::string> &built, const std::set<std::string> &relabelled, double b_ceiling,
                         int round);

// The N, O and S of a metal's inner sphere in the model's own copy (no symmetry mate): within the library's distance
// plus INNER_SHELL_MARGIN, at the coordination of those within an octahedral site's - the atoms a metal's records
// name. `own` is the metal's residue.
struct InnerSphereDonor {
    int chain = 0, residue = 0, atom = 0;   // indices into the model
    double distance = 0;
};
std::vector<InnerSphereDonor> InnerSphereDonors(const gemmi::Model &model, const gemmi::Residue &own,
                                                const gemmi::Atom &metal);

// The element of a species by its residue name (HOH: O, MG: Mg, ...).
gemmi::El SolventSpeciesElement(const std::string &species);

// A record (struct_conn metalc, named rugnux_metalcN with N on from `record`) for each of InnerSphereDonors() of the
// single-atom ion `res` of `chain`, where it is a metal.
void AddMetalRecords(gemmi::Structure &st, const gemmi::Chain &chain, const gemmi::Residue &res, int &record);

// The round's waters taken out and put in: each added one to the chain of its nearest atom, numbered on from that
// chain's last residue (its `site` is set to the name it gets). A relabelled site gets its new residue name,
// element and atom name; its metal records go, and a metal gets one for each N, O or S of its inner sphere (within
// the library's distance plus 0.33 A, as the records of the model as given are taken).
void ApplySolventEdits(gemmi::Structure &st, std::vector<SolventEdit> &edits);