The paper in one minute
They replace hand-tuned atomic partial charges with a probability distribution learned from quantum-level simulations, and make that affordable by swapping the slow simulations for Gaussian-process surrogates.
Classical force fields model molecules as balls (the atoms) joined by springs (the bonds), and give every atom a fixed point charge: a small electric charge placed at the atom's centre, set once and never updated during the simulation (more in Background, idea 3). The charges are usually picked by fitting and hand-tuning, which gives one "best" set with no error bars. This paper treats the charges as unknowns and infers them with Bayes' theorem. The reference data come from ab initio molecular dynamics (AIMD): simulations where the forces are computed from quantum mechanics ("ab initio" means from first principles). Each small molecule is simulated surrounded by explicit water molecules, so the effect of the solvent is built in. The inference uses Markov chain Monte Carlo (MCMC), an algorithm that tries hundreds of thousands of charge guesses and keeps them in proportion to how well they fit. Running a classical simulation for every guess would be far too slow, so they train Gaussian-process surrogates: cheap statistical models that predict the simulation outputs (radial distribution functions, hydrogen-bond counts, ion-pair distances) directly from the charges. They do this for 18 molecular fragments that make up proteins, nucleic acids and lipids.
The learned charges match the reference better than CHARMM36, a standard protein force field (in its "nbfix" variant, explained below). They reproduce experimental solution densities within 1%, and their error bars follow chemical intuition. As a test on a real protein, they build cardiac troponin C, the calcium sensor of heart muscle, from the fragment results and compute how strongly it binds a calcium ion (Ca²⁺). The result lands close to experiment, and the spread of plausible charges shows how much that number can move.
Background from zero
Seven ideas carry the whole paper. If you can explain each of these in two sentences, you can follow any discussion of it.
1. Molecular dynamics (MD)
Put every atom of a system (a protein, its water, some ions) in a box, compute the force on each atom, move all atoms a tiny step (about 1 femtosecond, fs, which is 10⁻¹⁵ s) using Newton's laws, repeat millions of times. The output is a trajectory: a movie of atomic positions. Averages over that movie give structure (which atoms sit near which), thermodynamics (density, binding free energies) and dynamics.
2. A force field is the energy function
MD needs the force on every atom. It gets them from one function, the potential energy \(U\): a single number (in kJ/mol) for how much energy is stored in the current arrangement of all the atoms. Stretch a bond or push two atoms too close and \(U\) goes up. The force on each atom points downhill, in the direction that lowers \(U\) fastest: \(\mathbf{F}_i = -\partial U / \partial \mathbf{r}_i\), where \(\mathbf{r}_i\) is that atom's position.
A force field is the recipe for \(U\) plus its fitted numbers (the parameters). CHARMM36 (Chemistry at HARvard Macromolecular Mechanics, version 36) is a standard one for proteins. Its recipe adds up five kinds of terms:
A sum sign \(\sum\) means "add this up over every item written underneath it". Term by term:
- Bonds. Every covalent bond is a spring. \(b\) is its current length, \(b_0\) its preferred length and \(k_b\) its stiffness. The energy grows with the square of the stretch, like a spring in school physics.
- Angles. The same idea for the angle \(\vartheta\) between two bonds that share an atom, such as H–O–H in water. (This \(\vartheta\) is unrelated to the \(\theta\) used later for the charges.)
- Dihedrals. Twisting around a bond, measured by the angle \(\phi\) between four atoms in a row. The cosine gives \(n\) preferred twist positions, such as the staggered positions around a C–C bond.
- Coulomb. \(\sum_{i \lt j}\) means "every pair of atoms \(i\) and \(j\), each pair counted once". In practice pairs already joined through one or two bonds are left out, because the bonded terms cover them. Each pair adds the product of its two charges divided by the distance \(r_{ij}\) between them. Same signs give positive energy, so the atoms push apart; opposite signs give negative energy, so they pull together. It only fades as \(1/r\), so it reaches far. \(\varepsilon_0\) is a physical constant (the vacuum permittivity). This is the only term the paper changes, through the charges \(q_i\).
- Lennard-Jones (LJ). Also summed over every pair. The \(r^{-12}\) part is a steep wall that stops atoms overlapping. The \(-r^{-6}\) part is the weak, short-range attraction that all atoms feel (dispersion, part of the van der Waals forces). \(\sigma_{ij}\) is roughly the contact distance, where the energy crosses zero, and \(\varepsilon_{ij}\) is the depth of the energy well, i.e. how sticky the pair is.
Terms 1–3 are the bonded terms: they act only along the chemical bonds. Terms 4–5 are the non-bonded terms: they act between all pairs, including between the molecule and the surrounding water. CHARMM36 adds a few small correction terms not shown here. This paper learns only the charges \(q_i\) and keeps every other parameter at its CHARMM36 value.
3. Partial charges cannot be measured
A partial charge is a fraction of an electron's charge assigned to each atom (in units of the elementary charge \(e\), the charge of one proton, about 1.6 × 10⁻¹⁹ coulombs: a unit of electric charge, not energy, so −0.6 e means 60% of one electron's charge), a crude summary of where the electron cloud sits. The oxygen in a carboxylate might carry about −0.6 e. No experiment measures it directly, and many different charge sets produce the same observable. The paper calls this a many-to-one mapping, which makes the problem non-unique. The usual approach is to fit charges to a quantum calculation of the molecule in vacuum and then patch them with ad hoc corrections for being in water.
What "fixed point charge" means. Point: the charge sits at a single spot, the centre of the atom, instead of being spread out like a real electron cloud. Fixed: the value is written once in the parameter file and never changes during the simulation, whatever is nearby. Real electrons shift when, say, a Ca²⁺ ion comes close (that is polarization, idea 5), and fixed charges cannot. Example: in the TIP3P water model used with CHARMM, the oxygen carries −0.834 e and each hydrogen +0.417 e. They add up to zero because water is neutral. Every pair of charged atoms then attracts or repels by Coulomb's law, the Coulomb term in the formula above.
4. AIMD versus FFMD
AIMD (ab initio MD) computes forces from the electrons at every step, here with density functional theory (DFT). It captures polarization and chemistry but costs so much that systems stay tiny (128 water molecules in a box 1.6 nanometres, nm, wide) and short (tens of picoseconds). FFMD (force-field MD) uses the cheap classical formula above, so it can run whole proteins for microseconds, but it is only as good as its parameters. The strategy: use AIMD on small fragments as the teacher, learn classical parameters from it, then reuse those fragment parameters in large systems.
5. Polarization and the electronic continuum correction (ECC)
Real electron clouds shift in response to their surroundings. Fixed charges cannot do that, so charged groups in water interact too strongly in standard force fields. ECC is a mean-field fix: scale the charges of ionic groups by roughly \(1/\sqrt{\varepsilon_{el}}\), where \(\varepsilon_{el} \approx 1.78\) is water's electronic dielectric constant, giving about 0.75. The Jungwirth and Martinez-Seara groups are major developers of this idea (prosECCo75 is their version of CHARMM36 with charges scaled by 0.75). In this paper the total molecular charge is scaled by 0.8, with matching water and ion models from their earlier work.
6. Quantities of interest (QoIs)
These are the simulation outputs used to compare classical and AIMD runs:
- Radial distribution function (RDF), g(r): how likely you are to find a water oxygen at distance r from a given solute atom, relative to bulk water. Peaks are hydration shells. It tends to 1 at large r.
- Hydrogen-bond counts: a geometric rule between a donor (the O or N holding the hydrogen) and an acceptor (the O, N or S receiving it): distance < 3.5 Å (1 ångström = 0.1 nm) and donor–H–acceptor angle > 150°.
- Ion–solute distance distributions: with a counterion (an ion of opposite charge) held in place by a restraint, either in direct contact or with a water layer in between (solvent-shared).
7. Bayesian inference, MCMC and surrogates
Bayesian inference turns "which parameter values explain the data?" into a probability distribution. You start with a prior (which values are plausible before seeing data), multiply by the likelihood (how well each value reproduces the data), and get the posterior (which values are plausible after seeing the data). The width of the posterior is the error bar. Its single most probable point is called the MAP (maximum a posteriori) estimate.
MCMC (Markov chain Monte Carlo) is how you get the posterior in practice when there are several parameters. It is a guided random walk through parameter space: propose a small change, accept it more often if the posterior is higher there, repeat. After enough steps, the visited points are samples from the posterior. It needs a very large number of likelihood evaluations.
A surrogate (or emulator) is a cheap model trained to imitate an expensive one. Here a Gaussian process (GP), a flexible statistical model that interpolates smoothly between training examples and reports how uncertain it is, stands in for running a full simulation at every MCMC step.
Where your background already fits
You are new to molecular simulation, but the core of this paper is statistics and software, and you have both. These are honest connections you can make in the interview.
Four standard tools, and where each one shows up in the method:
- Emulator and Bayesian calibration (the whole method). An expensive simulation is replaced by a cheap statistical model of it, an emulator (this paper calls it a surrogate), and the simulation's unknown inputs are then inferred from reference data. Statisticians call this Bayesian calibration of a computer model; the classic paper is Kennedy and O'Hagan (2001). That label is my framing, not the paper's: the paper does not cite Kennedy and O'Hagan, though it does say the QoIs are "emulated" and it cites earlier work titled "Bayesian Calibration of Force-Fields from Experimental Data" (Dutta et al., 2018, on water). Useful if an interviewer asks how this relates to methods you know.
- Latin hypercube sampling (stage a, choosing the training runs). Before any surrogate exists, the method needs several thousand trial charge sets to simulate. The paper picks them by Latin hypercube sampling from SciPy, inside chemically sensible bounds. In the code it is the
RandomParamsGeneratorinbff/domain/specs.py, which draws Latin hypercube points and throws away any that break the charge rules. "Space-filling design" is the design-of-experiments name for this goal: spread the runs evenly over the space, because a Gaussian process predicts well only near its training points. A gap in the training set means poor predictions there, and the sampler could wander into it. How it works is in the stage a box. - Leave-one-out (stage b, tuning the surrogate). The GP's own settings (length scales, signal size, noise) are chosen to maximize a leave-one-out (LOO) score: for each training run, how well the GP predicts it from all the other runs. That is the paper's Eq. 4, plus log-normal priors on the settings. For a GP this needs no refitting: all N "leave one out" predictions come from one matrix inverse. In the code it is
loo_log_likelihoodinbff/bayes/likelihoods.py. Strictly it is a leave-one-out marginal likelihood, a close cousin of leave-one-out cross-validation. - PyTorch on graphics cards (stages b and c). GP fitting, the surrogate predictions and the newer MCMC sampler are written in PyTorch, so they run on graphics processing units (GPUs).
RDFs are noisy 1-D curves, like a filtered EEG trace, except the horizontal axis is distance r instead of time. The newer code (v0.3.0 onward, not in the paper) uses a small signal-processing routine on them.
The problem it solves. The likelihood has to decide how much one curve counts. The paper counts every curve as one observation. Counting every grid point (say 51) is overconfident, because neighbouring points move together. The truth is in between: roughly the number of distinct features the curve has. The code calls that number n_eff, the effective number of observations, and estimates it from the reference curve. Details in the likelihood deep dive.
How it estimates it (estimate_curve_n_eff in bff/bayes/effective_observations.py, about 20 lines of SciPy):
- Smooth. A Savitzky–Golay filter slides a window of 15 grid points along the curve, fits a cubic polynomial inside each window, and keeps the fitted value at the centre. Unlike a plain moving average, it keeps peak heights and widths, which is why it is popular for EEG and fNIRS too.
- Estimate the noise. Take the median of the absolute differences between the raw and the smoothed curve. The median makes it robust to a few large spikes, like the median absolute deviation you may have used for artefact thresholds.
- Find features.
scipy.signal.find_peaksfinds peaks (hydration shells) and, run on the flipped curve, troughs (the gaps between shells). Only features with a prominence (height above their surroundings) of at least 5 times the noise count. That is the same as rejecting EEG peaks that don't stand clear of the noise floor. - Score them. Each feature adds \(\max(1, \ln(\text{prominence}/\text{tolerance}))\): at least 1, a little more if it is very pronounced. The tolerance is a scale you set in the learn configuration (for example
tolerance: 0.1). The total is at least 1. If one QoI holds several curves, say one RDF per solute atom, their counts are added.
For the example RDF in the deep dive, this gives about 6, so the curve counts as 6 observations instead of 1 or 51. Two differences from EEG work worth knowing: it runs once, on the reference curve before sampling, not on every prediction; and the noise is statistical noise from a finite simulation, not sensor noise. Because the window is 15 points rather than a fixed distance, the result also depends on the grid spacing, a point you could raise in the interview.
bfflearn), is a command-line tool configured with YAML files (a plain-text config format). It drives GROMACS (the classical MD program) and CP2K (the quantum-chemistry program), submits jobs to SLURM (a cluster job scheduler), analyses trajectories with the MDAnalysis Python library and fits surrogates in PyTorch. An intern will very likely touch this workflow code.The method, stage by stage
The letters match panels a–c of the paper's Figure 1, plus the validation step.
aData acquisition
- Each molecule is placed in a ~1.6 nm cubic box with 128 water molecules of the TIP4P/2005 model (a standard rigid water model with four interaction sites), plus counterions if the molecule is charged. Starting parameters are CHARMM36 topologies (files listing atoms, bonds and parameters) made with the CHARMM-GUI web tool. The box is equilibrated (run until it settles) with classical MD.
- Reference: equilibrated snapshots seed AIMD runs in CP2K. The DFT approximation used is the revPBE functional plus the D3 correction for dispersion (weak van der Waals attraction), at 300 K with a 0.5 fs time step.
- Training set: several thousand charge vectors \(\theta\) are drawn by Latin hypercube sampling (explained in the box below) inside chemically sensible bounds. One atom type per molecule is "implicit": its charge is whatever keeps the total molecular charge fixed.
- For every \(\theta\), GROMACS runs three setups: solute alone, solute with a restrained ion in contact, and solute with a restrained ion one water layer away.
- QoIs are stacked into a matrix \(Y\) (one row per charge set). The AIMD QoIs form the target vector \(y\).
Why three setups: water structure alone does not pin down how the molecule interacts with ions, and ion binding is the downstream use case.
bSurrogate modelling
- The bottleneck: MCMC needs tens to hundreds of thousands of likelihood evaluations. Each would need a full MD run. Not feasible.
- The fix: learn a cheap function from charges to QoIs. They use a local Gaussian process (LGP), from Shanks et al. (Journal of Chemical Theory and Computation, 2024). Each output point (for example the RDF value at one radius \(r_i\)) is modelled by its own independent GP over the charge vector. All points of one QoI share a kernel, so the big matrix inverse is computed once per QoI and prediction becomes a matrix multiply on the GPU.
- Kernel (the function that says how similar two charge sets are, and so how much their outputs should agree): squared-exponential, also called radial basis function (RBF), with one length scale per charge. A short length scale means the output is sensitive to that charge.
- Hyperparameters, the GP's own settings (kernel variance, length scales, noise), are fitted by maximizing a leave-one-out log marginal likelihood with log-normal priors, using gradient descent in PyTorch, on an 80/20 split.
- RDF trick: subtract a sigmoid baseline (it rises from 0 to 1 around 3 Å) so the GP only learns deviations from generic RDF shape. Equivalent to giving the GP that sigmoid as its prior mean.
- They checked whether averaging over hyperparameter uncertainty (a 100-member committee) changed the charge posterior. It did not, so a single fitted GP is used.
cBayesian inference
- Prior on charges: normal, centred at the middle of each charge's allowed range, standard deviation = range / 5. Zero probability outside the bounds or if total charge is violated. This keeps the sampler away from edges where the surrogate extrapolates. A uniform prior gave essentially the same answer, so the data dominate.
- Nuisance parameters \(n_k\): one noise scale per QoI, which soaks up surrogate error plus reference noise. Log-normal prior, learned together with the charges. In effect, the model decides how much to trust each QoI instead of you hand-picking loss weights.
- Likelihood: Gaussian residuals between predicted and AIMD QoIs, independent between QoIs, same noise at every grid point of a QoI, each whole curve counting as one observation.
- Sampler:
emcee, an affine-invariant ensemble sampler (Goodman–Weare stretch move). Walkers = 5 × number of dimensions. Run until the chain is at least 100 integrated autocorrelation times (τ, roughly how many steps apart two samples must be to count as independent) long and the estimate is stable to 1%, up to 100,000 iterations. Discard 2τmax as burn-in and thin by 0.5τmin.
dValidation
- Fit a multivariate normal (Laplace approximation) to the posterior, draw 10 charge sets, and run new classical MD with them. No surrogate involved at this point.
- Compare the resulting QoIs to AIMD with the normalized mean absolute error (NMAE), and do the same for CHARMM36-nbfix: CHARMM36 plus NBFIX ("non-bonded fix") terms, hand-set Lennard-Jones corrections for specific atom pairs that weaken overly strong ion binding.
- Separate checks against experiment: solution densities and the troponin Ca²⁺ binding free energy.
The equations, annotated
Five equations carry the method. You do not need to derive them, but be able to say what each symbol is.
Under each symbol, in grey, is how to say it out loud, and each equation has a "Read aloud" line you can use if you talk through it in the interview. Subscripts are usually just read after the letter: \(n_k\) is "n k" (or, more formally, "n sub k"). The symbols table in the glossary lists the Greek letters too.
Read aloud: "p of theta and n given y equals p of y given theta and n, times p of theta, times p of n, all over p of y."
- \(p(a \mid b)\)"p of a given b"
- the probability of \(a\) if \(b\) is true; the bar \(\mid\) is read "given"
- \(\theta\)"theta"
- the partial charges being learned, as one list. Each entry is one partial charge, written \(q\) ("q") elsewhere on this page, so for acetate \(\theta\) holds a \(q\) for the oxygens, one for the carboxyl carbon, and so on.
- \(n\)"n"
- nuisance noise parameters, one per QoI
- \(y\)"y"
- reference QoIs from AIMD
- \(p(\theta), p(n)\)"p of theta", "p of n"
- the priors: what was plausible before seeing the data
- \(p(y)\)"p of y"
- evidence; a constant here, so MCMC can ignore it
Read aloud: "p of H p given E over p of H d given E equals p of E given H p over p of E given H d, times p of H p over p of H d." Or in words: posterior odds equal the likelihood ratio times the prior odds.
- \(H_p, H_d\)"H p", "H d"
- the prosecution and defence hypotheses
- \(E\)"E"
- the evidence
Same theorem. MCMC works with ratios too: each accept/reject step compares posterior values at two points, and the unknown \(p(y)\) cancels.
Read aloud: "p of y given theta and n equals the product, over k from 1 to capital K, of one over two pi n k squared, to the power n obs k over two, times exp of minus one over two n k squared, times the norm of y hat k of theta minus y k, squared."
- \(\prod_{k=1}^{K}\)"the product over k from 1 to capital K"
- multiply the bracketed term for every QoI. \(\prod\) is a capital Greek pi used as a multiplication sign, like \(\sum\) is used for adding.
- \(k\), \(K\)"k", "capital K"
- which QoI (an RDF, an H-bond count, a distance distribution); \(K\) is how many QoIs there are
- \(\hat{y}_k(\theta)\)"y hat k of theta"
- surrogate prediction for that QoI at charges \(\theta\). The hat means "estimated" or "predicted".
- \(y_k\)"y k"
- the AIMD reference value of that QoI
- \(n_k\)"n k"
- noise standard deviation for QoI \(k\), learned
- \(n^{\text{obs}}_k\)"n obs k"
- number of independent observations; set to 1 per QoI in the paper
- \(\pi\)"pi"
- the number 3.14159…, from the formula of the bell curve
- \(\exp(x)\)"exp of x" or "e to the x"
- the exponential function, \(e^x\) with \(e \approx 2.718\) (not the elementary charge)
- \(\lVert x \rVert^2\)"norm of x, squared"
- the sum of the squares of the entries of \(x\); here the squared misses added over the grid
Read it as a weighted sum of squared errors: small \(n_k\) means "trust this QoI a lot". The \((2\pi n_k^2)^{-1/2}\) factor stops the sampler from making every \(n_k\) huge to excuse bad fits.
Read aloud: "y hat k of theta equals K theta X, times the inverse of K X X plus sigma k squared I, times capital Y k." The superscript \((k)\) on each \(K\) just says "for QoI k" and is usually left out when speaking.
- \(X\)"capital X"
- the N training charge sets
- \(Y^{(k)}\)"capital Y k"
- their simulated values of QoI \(k\) (N × grid points)
- \(K_{\theta,X}\)"K theta X"
- similarity of the new \(\theta\) to each training point
- \(K_{X,X}\)"K X X"
- similarity of every training point to every other one (the N × N kernel matrix)
- \(\sigma_k^2 I\)"sigma k squared I"
- the GP's own noise term: a small number \(\sigma_k^2\) added along the diagonal. \(I\) is the identity matrix (1 on the diagonal, 0 elsewhere). Not the same as the likelihood noise \(n_k\).
- \([\,\cdot\,]^{-1}\)"the inverse of …"
- the matrix inverse. \([\,\cdot\,]^{-1} Y\) is computed once and cached, which is why each prediction is cheap
Intuition: the prediction is a similarity-weighted blend of training simulations.
Read aloud: "K i j equals alpha k squared times exp of minus one half times the sum, over d from 1 to capital D, of theta d i minus theta d j, over l k d, all squared."
- \(K_{ij}\)"K i j"
- how similar training charge sets \(i\) and \(j\) are; one entry of the kernel matrix
- \(\sum_{d=1}^{D}\)"the sum over d from 1 to capital D"
- add up over every charge \(d\); \(D\) is how many charges are learned (up to 10 here)
- \(\theta_{d,i}\)"theta d i"
- charge number \(d\) in training set number \(i\)
- \(\alpha_k^2\)"alpha k squared"
- vertical scale of the function
- \(l_{k,d}\)"l k d" (often "ell k d")
- length scale for charge \(d\): how far you can move it before the output changes
Read aloud: "N-M-A-E k equals the one-norm of y hat k minus y k, over the one-norm of y hat k plus the one-norm of y k." NMAE (normalized mean absolute error) is spelled out letter by letter.
- \(\lVert x \rVert_1\)"one-norm of x"
- the sum of the absolute values of the entries of \(x\): add up the misses, ignoring their sign
0 = perfect match, 1 = complete disagreement, reported as a percentage. Normalizing lets RDFs, counts and distributions share one scale.
Toy posterior: why some charges are pinned down
One charge, one QoI. The QoI changes linearly with the charge at a rate you set (sensitivity). The reference says the best charge is −0.62 e. The prior is the paper's recipe: bounds −1 to 0, normal centred at −0.5 with standard deviation 0.2. Watch the posterior width.
This is the mechanism behind Figure 4: oxygens that touch water change the RDFs a lot, so the data pin them down. Buried carbons barely affect any QoI, so their posterior stays close to the prior. Wide intervals are useful: they mark where you can put excess charge without hurting the fit, which is exactly how the troponin model was neutralized.
Results
The 18 fragments
These are the building blocks of backbones, side chains, nucleic acid phosphates and lipid head groups.
Agreement with AIMD (Figure 3b)
| QoI | Typical NMAE | Note |
|---|---|---|
| RDFs (hydration structure) | < 5% | for most species |
| Hydrogen-bond counts | 10–20% | the rigid water model limits fitting RDFs and H-bonds at the same time |
| Ion–solute distances | < 20% | somewhat larger for some anions |
- Better than CHARMM36-nbfix for nearly all species and QoIs. Biggest gains for charged molecules, especially anions. Neutral molecules improve a little.
- Fitting each QoI alone gives better per-QoI numbers (Supplementary Fig. S3). The joint fit is a compromise, and the remaining error reflects that trade-off.
Transferability: solution densities (Figure 3c)
Training boxes hold one solute, so solute–solute contacts were never seen. Still, simulated densities of concentrated solutions (sodium acetate, ethanol, ammonium chloride, potassium phosphates and formate, guanidinium chloride, acetamide) are within 1% of experiment, up to roughly 50% solute by mass.
Chemical intuition from the posterior (Figure 4)
- Narrow 95% intervals for oxygens exposed to water; wide for buried carbons and phosphorus. Exposed atoms control solvation structure.
- Learned charges span about −0.8 e to +0.9 e, narrower than CHARMM36-nbfix (−1.0 e to +1.5 e), partly because of the 0.8 scaling.
- Carboxylate and phosphate terminal oxygens cluster at −0.65 to −0.55 e; sulfate oxygens about −0.4 e (charge spread over three equivalent oxygens); carbonyl oxygens −0.7 to −0.4 e; hydroxyl oxygens about −0.7 e; ether oxygens −0.3 to −0.1 e.
- Nitrogens −0.9 to −0.4 e. Hydrogens on O or N clearly positive, aliphatic ones weakly positive. Carbonyl carbons strongly positive.
- Each molecule was fitted independently, yet the same chemical groups converge on similar charges. The authors present this as evidence the method learns chemistry rather than noise.
Proof of principle: Ca²⁺ binding to cardiac troponin C
Biology: cardiac troponin C (cTnC) is the calcium sensor of heart muscle. Its regulatory domain at the N-terminal end of the chain (N-cTnC; crystal structure 1AP4 in the Protein Data Bank) binds Ca²⁺ in an EF-hand loop. That is a helix–loop–helix motif, named after the E and F helices of parvalbumin, where it was first described. Calcium binding there starts contraction. The carboxylate groups (COO⁻) of aspartate (Asp) and glutamate (Glu) side chains grab the ion.
Building the protein from fragments: backbone parameters stay fixed. Charged side chains come from fragment posteriors: Asp and Glu from acetate, arginine (Arg) from guanidinium, lysine (Lys) from ethylammonium, all with 0.8 scaling. To keep the protein neutral, up to ±0.1 e is added to the atom with the widest posterior, usually a buried carbon. Several charge sets (the MAP, i.e. the most probable set, plus 10 random posterior samples) give several protein models.
Free energy method: alchemical double decoupling. You gradually switch off the ion's interactions once in the binding site (held by flat-bottom restraints) and once in bulk water, and close the thermodynamic cycle with corrections for the restraints and for the net charge change in a periodic box. The switching is controlled by a parameter λ, run in 21 steps called λ windows (11 for electrostatics, 10 for van der Waals). Free energy differences between neighbouring windows are estimated with the Bennett acceptance ratio (BAR).
ΔGbind is the binding free energy in kJ/mol. More negative means tighter binding.
| Model | ΔGbind (kJ/mol) | Reading |
|---|---|---|
| Experiment | −28.6 | target |
| This work, posterior samples | within 10 | of experiment; best sets within ~1 |
| This work, MAP charges | ≈ 6 weaker | underbinds: best charges for acetate in water are not best in the protein |
| CHARMM36 | −105.7 | far too strong; no polarization correction |
| CHARMM36-nbfix | −31.7 | close, but through error compensation: pair-specific Lennard-Jones tweaks offset overly strong electrostatics |
| prosECCo75 | −13.5 | too weak |
- ΔGbind correlates with the carboxyl oxygen charge, the atom touching the calcium (Figure 5g).
- Charge sets whose CH2COO⁻ dipole moment is 2.7–4.4 debye (D, the unit of dipole moment, a measure of charge separation) reproduce experiment best. CHARMM36 has 8.9 D and prosECCo75 6.6 D.
- The main point: a single best-fit charge set would have given one number (the MAP, off by ~6). The posterior gives a family of plausible models and shows how sensitive the prediction is to the charges. The authors call the MAP result a basic limitation of fixed-charge models: they cannot adapt to a new chemical context.
Limitations
The first six are stated in the paper. The rest are fair points you could raise yourself, which shows you read critically.
What changed in the code since the paper
The paper snapshot is tag v0.0.1. The repository has moved to v0.4.1 (25 August 2026). Mentioning one or two of these shows you looked beyond the PDF. I read the README, changelog, commit log and file list, not the full source, so skim the docs yourself before the interview.
| Change | Why it matters |
|---|---|
| Reference MD via a machine-learned interatomic potential (MLIP): CP2K labels snapshots, an MLIP is trained outside BFF, then runs the reference MD | MLIPs give close to DFT accuracy much more cheaply, so reference runs can be longer and larger. Addresses the "small, short AIMD" limitation. |
| Effective observation counts for curves (v0.3.0) | Replaces "one observation per curve" with an estimate from the number of resolved peaks. Directly relaxes a likelihood assumption from the paper. See the likelihood deep dive. |
| QoI-attributed marginals | Shows which observable supports which region of the posterior. |
| Lennard-Jones σ and ε, and dihedral force constants, as learnable parameters | The paper's "exploratory LJ work" is now a feature. |
| Hierarchical residue- and system-level charge constraints (v0.2.0) | Needed to build whole residues and proteins from fragments. |
| Own PyTorch MCMC sampler with rank-normalized split-R̂ convergence checks | The paper used emcee with autocorrelation time. R̂ ("R-hat") compares independent chains; values close to 1 mean they agree, which suggests convergence. |
Worked examples in the repo: acetate (full workflow), arbitrary tabular data (notebook), and a neon Mie potential learned from published RDFs. Running the arbitrary-data or neon notebook before the interview would give you something concrete to talk about.
Paper, current code, and pull request 18
Three documents sit next to each other, and they do different jobs. The paper learns charges and reports a posterior. Version 0.4.1 is that method as maintained code. Pull request 18 proposes extra modules for a cheaper campaign, and its description links to the draft at bff.think2earn.com. The pull request is open, with no reviews. It was opened on 26 August 2026, the day after the v0.4.1 release, and it only adds files: 376 lines, nothing deleted.
- The paper (Kostál et al., J. Chem. Theory Comput. 2026) learns partial charges for 18 biomolecular fragments from ab initio molecular dynamics in water, then checks solution densities and calcium binding to cardiac troponin C. This whole guide is about that paper.
- The current code is
mainat tagv0.4.1. Same loop as the paper, plus the changes in what changed in the code: a local Gaussian process in PyTorch, a Gaussian likelihood, and the group's own sampler. Molecular dynamics runs locally or through SLURM. There is nobff/backends/directory onmain. - The pull request adds four modules beside that loop, plus three unit tests. The systems in its writeup are two small solutes in water: acetate, with two carboxyl charges and the methyl charge fixed so the ion stays at −1, and phenol with out-of-plane virtual sites (charge points that are not atoms, placed off the ring to stand in for the electron cloud).
| Question | Paper | Code v0.4.1 | Pull request 18 |
|---|---|---|---|
| Which simulations are run? | Several thousand charge sets, drawn once by Latin hypercube sampling, then all simulated before any inference. | The same design. RandomParamsGenerator draws the hypercube and drops any set that breaks the charge rules. |
The text describes active learning on phenol: 40 simulations compared with grids of 250 and 1,500. The diff contains no code that chooses the next simulation. |
| What is the surrogate? | One local Gaussian process per quantity of interest. Each bin of a curve is its own output, and the bins share one squared-exponential kernel, so there is one set of length scales and one matrix inverse per quantity of interest. Stage b of the method says this. | LocalGaussianProcess in bff/bayes/gaussian_process.py, kernel gaussian_kernel. One predict call returns the whole output vector. A committee of hyperparameters was tried in the paper and did not move the posterior, so one fitted model is the default. |
FPCASurrogate compresses one curve by a singular value decomposition into M shape-modes (default 5) and fits M separate scikit-learn Matérn Gaussian processes, one per coefficient, each with its own length scales. It never calls LocalGaussianProcess, and the learn command never calls it. |
| Does model uncertainty decide where to simulate? | No. The training set is fixed first. During sampling, surrogate error is absorbed into the noise level nk, not used as a reject rule. | No gate. Same likelihood path. | PlausibilityGateV2 keeps a proposal only when the absolute miss, plus κ times the Gaussian-process standard deviation, stays inside τ times the target noise (defaults κ = 2, τ = 3). Dividing by that standard deviation, a Mahalanobis ratio, would do the opposite: a huge error bar makes any miss look small, so an unvisited corner looks acceptable. Neither the paper nor v0.4.1 implements that ratio. The new test rejects a prediction that already matches the target when the error bar is huge. It does not try a charge set far from the training data. |
| Logarithms and Jensen's inequality | The likelihood is Gaussian on the observable itself: a sum of squared misses. No logarithm of the observable. | Still Gaussian on the observable. The sum of squares became a mean squared error times neff. Here neff is 1, the number of grid points, or a count of peaks in the reference curve. | JensenLogObservableCorrection is for a model trained on the log of the observable. The average of a log sits below the log of the average, so a short trajectory biases that model low. The module turns the simulation variance into a noise level on the log, with its own neff = (number of frames) / (2τcorr + 1), and converts back by exp(μ + σ²/2). The learn command never trains on logs, and this neff is not the peak count above. |
| Reusing an old trajectory | The Bennett acceptance ratio (BAR) estimates the troponin binding free energy between neighbouring λ windows. It does not replace a training simulation. | The training loop still runs a fresh simulation for every training charge set. | No reweighting code. The linked draft describes the multistate Bennett acceptance ratio (MBAR): reweight frames you already saved to estimate a nearby parameter set, while the ensembles still overlap, and quotes about 0.0015 seconds for one such step. That code is not in this pull request. |
| Where does the molecular dynamics run? | Classical molecular dynamics, GROMACS in the released workflow, with cluster jobs. | job_scheduler is local or slurm only. |
ModalGromacsBackend is meant to fan GROMACS out on Modal, a serverless cloud. As committed, it imports new_optimizations.modal_bff_campaign, a module that is not in the BayesicForceFields repository, and the charge vector it is given is not copied into the job. No existing workflow calls it. The pull request text also mentions 16 concurrent runs on local NVMe under SLURM. That path is not in the diff. |
| How is the posterior sampled? | emcee, with the autocorrelation time as the convergence check. |
The group's PyTorch sampler, with rank-normalized split-R̂. | No sampler in the diff. The linked draft says the No-U-Turn Sampler in NumPyro, on a surrogate that can be differentiated. The new models are scikit-learn objects, so they are not that PyTorch surrogate. |
| What was actually checked? | 18 fragments, densities within about 1%, and the troponin free energy. | Tests for the local Gaussian process, the likelihood, and the peak-counting neff, plus the acetate example in the repository. | Three tests on made-up arrays: a synthetic curve, two hand-set error bars through the gate, and the log-noise formula. The phenol table in the pull request (40 simulations within 0.14% of a 1,500-point grid, 3.2 minutes, a 37.5-fold cut) and the draft's acetate timing (50 simulations, about 3 minutes) are claims in the writeup. The tests do not produce those numbers. |
The useful way to say this in the interview: the paper and v0.4.1 answer "which charges fit the reference, and how wide is the posterior?" Pull request 18 answers "could a later campaign simulate fewer points, compress each curve, and run the batch somewhere other than a SLURM queue?" Those modules are a proposal sitting next to the pipeline. They do not replace the local Gaussian process, the likelihood, or the job runner, and the 37.5-fold and 3-minute figures belong to the writeup, not to a test in the diff.
One line in the pull request is easy to repeat and easy to get wrong. It says that modelling an RDF as 200 bins means 200 Gaussian processes and 200 separate hyperparameter fits, and that five modes are a 97.5% cut. In the paper and in v0.4.1 the bins share one kernel per quantity of interest, so the hyperparameter count does not grow with the number of bins. Five Matérn models are fewer outputs than a full curve. They are not a 97.5% cut of the model the code fits today.
Deep dive: the likelihood, and a first intern project
Terms used in this section (32)
- Likelihood
- A score for a trial charge set: how probable the reference data would be if these charges were right. Higher means a better fit.
- QoI (quantity of interest)
- One simulation output compared with the reference: an RDF, an H-bond count, an ion-distance distribution.
- Curve
- A QoI stored as many numbers along a distance r, like an RDF.
- Grid point (bin)
- One of the fixed r values where a curve is stored, for example every 0.1 Å. The grid spacing is the gap between them.
- Reference curve
- The QoI from the AIMD (or MLIP) reference simulation: the target to match.
- Prediction
- The surrogate's estimate of a QoI for a trial charge set.
- Residual
- Prediction minus reference at one grid point.
- Sum of squares
- The squared residuals added up over the grid. What the paper uses.
- MSE (mean squared error)
- The average squared residual: the sum of squares divided by the number of grid points. What the code uses now.
- Noise level (n_k in the paper, σ_k "sigma k" in the code)
- How big a mismatch to expect even with the right charges.
- Nuisance parameter
- A parameter the model needs but you don't care about at the end; here the noise level.
- Noise model
- The assumptions about the mismatch: its size, whether the size changes along the curve, and whether errors at different points or in different QoIs move together.
- Homoskedastic / heteroskedastic
- Noise of the same size everywhere / noise whose size varies, for example along r.
- Correlated errors
- Errors that move together, like neighbouring grid points of a smooth curve.
- Covariance matrix (Σ_k, "capital sigma k")
- A table of how much every pair of grid points varies together. Its diagonal is the noise variance at each point.
- Observation count (n_obs, "n obs")
- How many independent measurements the likelihood treats a QoI as. The paper uses 1.
- Effective observation count (n_eff, "n eff")
- An estimate of how many independent pieces of information a curve really holds. Fewer than its grid points, because neighbours are correlated.
- Curve weight
- How hard a QoI pulls on the charges. Grows with n_eff and shrinks with the square of the noise level.
- Peak / trough
- A local maximum / minimum of a curve. In an RDF, peaks are hydration shells and troughs are the gaps between them.
- Prominence
- How far a peak rises above its surroundings, or a trough dips below them.
- Tolerance
- The scale, in the curve's units, used when counting peaks. A feature much larger than the tolerance counts as more than one observation.
- Savitzky–Golay filter
- Smoothing by fitting a small polynomial in a sliding window. The code uses a window of up to 15 grid points.
- Surrogate uncertainty (predictive variance)
- The Gaussian process's own estimate of how unsure its prediction is at each point.
- "Folded into"
- Lumped together into one number instead of being modelled separately.
- Block averaging
- Split a trajectory into chunks, compute the QoI in each chunk, and use the scatter between chunks to estimate the noise.
- Shrinkage (Ledoit–Wolf)
- Blending a noisy covariance estimate with a simpler one so it becomes stable and can be inverted.
- Cholesky factor, triangular solve
- A standard way to compute with a covariance matrix: compute the factor once, and each later use is cheap.
- Calibration (coverage)
- Whether stated uncertainty is honest: 95% intervals should contain the truth about 95% of the time.
- Synthetic data
- Made-up reference data generated from charges you chose, so you know the right answer.
- MLIP (machine-learned interatomic potential)
- A cheaper stand-in for DFT, trained on DFT data, used for longer reference runs.
- Physics-informed GP (ref. 22)
- A Gaussian process built with physical knowledge of liquid structure. It gives a learned covariance for structure data.
- Posterior / marginal
- The distribution of plausible charges after seeing the data / the same for one charge on its own.
The likelihood is the part of the method that scores a trial set of charges against the reference data. Its formula, Eq. 7, decides how much each simulation output (each RDF, the H-bond count, the ion distances) influences which charges come out. The next subsection explains what that means.
The paper lists the simplifications in this formula itself. Comparing the code released with the paper (tag v0.0.1) with the current code (v0.4.1) shows where the group has gone since:
- Already changed: how much weight each curve gets. The paper counts every curve as one observation. The code now estimates an effective number of observations per curve.
- Still as in the paper: one noise level per QoI, noise treated as independent along each curve and between QoIs, Gaussian noise, and the surrogate's own error folded into that one noise level.
All the open items live in one function, gaussian_log_likelihood_by_qoi in bff/bayes/likelihoods.py, which already has tests. That is why they make a good first project. The side-by-side code is in paper versus code below.
What "pulling on the charges" means
MCMC wanders through possible charge sets. At each step it compares two candidates and moves to the one with the higher posterior score more often. That score is the log prior plus the log likelihood, and the log likelihood is a sum with one term per QoI. So every QoI casts its own vote: charge sets that reproduce its reference curve raise the score, and charge sets that miss it lower the score. After many steps, the walkers crowd where the combined score is high, and that crowd is the posterior.
When two QoIs disagree, the result is a compromise, like a tug of war. Each QoI pulls the charge toward the value it likes best. How hard it pulls depends on how sharply its score drops as you move away from that value, which is set by three things:
- its noise level \(n_k\): less noise, sharper drop, stronger pull;
- its observation count \(n^{\text{obs}}_k\): more observations, stronger pull;
- its sensitivity: how much the QoI actually changes when this charge changes. An RDF around an exposed oxygen is very sensitive to that oxygen's charge.
Try it: change the noise levels, or give the RDF more observations, and watch the posterior move.
In this toy, both QoIs are scaled so a noise of 0.02 means 0.02 e of charge. The H-bond count is a single number, so it always counts as one observation. Raising the RDF's observation count is what the newer code's \(n^{\text{eff}}\) does, and it shifts the balance toward the RDF.
What Eq. 7 does, step by step
Read aloud: "log p of y given theta and n equals the sum over k of: minus the norm of y hat k of theta minus y k, squared, over two n k squared, minus n obs k times log n k; plus a constant." "Const" is a constant that does not depend on the charges or the noise.
Read it from the inside out. The numbers below use a made-up five-point piece of an RDF so you can follow every step.
| r (Å) | Reference \(y\) | Prediction \(\hat{y}\) | Difference | Squared |
|---|---|---|---|---|
| 2.6 | 0.00 | 0.00 | 0.00 | 0.000 |
| 2.8 | 0.90 | 1.00 | +0.10 | 0.010 |
| 3.0 | 2.10 | 2.40 | +0.30 | 0.090 |
| 3.2 | 1.30 | 1.20 | −0.10 | 0.010 |
| 3.4 | 0.90 | 0.90 | 0.00 | 0.000 |
| Sum of squares \(\lVert \hat{y}_k - y_k \rVert^2\) | 0.110 | |||
- Predict. Pick a trial charge set \(\theta\). Instead of running a new MD simulation, the surrogate (the Gaussian process from stage b) returns the curve those charges would produce, as a number at each grid point r. This is a matrix multiplication, so it is fast enough to repeat hundreds of thousands of times.
- Compare point by point. Subtract the reference curve from the prediction at every r. Positive means the prediction is too high there, negative too low (the Difference column).
- Square and add. Squaring makes every miss positive, so misses in opposite directions cannot cancel out. It also punishes big misses much more than small ones: the 0.30 miss at 3.0 Å costs nine times as much as a 0.10 miss. Adding the squares gives one number for the whole curve, 0.110 here. The squared form is not arbitrary; it is what assuming Gaussian noise produces.
- Divide by the noise. \(n_k\) is how much this QoI would scatter even with perfect charges: statistical noise from the short reference simulation plus the surrogate's own error. Dividing by \(2n_k^2\) turns the raw mismatch into "how many noise widths off", like a z-score. The units cancel, which is what lets an RDF (no units) and an ion distance (in Å) be added together. A QoI believed to be noisy gets a small penalty for the same mismatch, so it pulls less.
- Pay for the noise you claim. Without the \(-n^{\text{obs}}_k \log n_k\) term, the sampler could make any fit look good by claiming huge noise. The term comes from the Gaussian's normalization: a wider bell curve is lower everywhere, so claiming more noise lowers the best score you can reach. For the example, with \(n^{\text{obs}}_k = 1\):
The score is highest when \(n_k\) matches the real size of the misses, here about 0.33, the square root of 0.110. Claim less noise and the misses are punished harshly; claim more and the noise price takes over. This is how the sampler learns \(n_k\) from the data.
Noise \(n_k\) Mismatch term \(-0.110/(2n_k^2)\) Noise price \(-\log n_k\) Total 0.10 −5.50 +2.30 −3.20 0.33 −0.51 +1.11 +0.60 (best) 1.00 −0.06 0.00 −0.06 3.00 −0.01 −1.10 −1.11 - Add up the QoIs. Repeat steps 1–5 for every QoI (every RDF, the H-bond count, each ion distance) and add the results. Adding logarithms is the same as multiplying probabilities, and multiplying the probabilities of different pieces of evidence is only valid if they are independent: assumption (i) below. It is the same rule as multiplying likelihood ratios for independent pieces of evidence in forensics.
- Add the prior and compare. The posterior score is this sum plus the log prior. MCMC only ever compares the scores of two charge sets. A trial set whose prediction misses by a sum of squares of 0.02 instead of 0.110 scores +1.02 instead of +0.60 at the same noise, which makes it about 1.5 times as probable (\(e^{0.42}\)). So the walkers drift toward it.
The assumptions, in plain words
| Assumption | What it means | Why it's not quite true |
|---|---|---|
| (i) QoIs are conditionally independent | Knowing the error in the RDF tells you nothing about the error in the H-bond count. | Both come from the same short AIMD trajectory. If that trajectory happened to sample unusually ordered water, the RDF peak and the H-bond count are off in the same direction. |
| (ii) Residuals are independent along a curve | The error at r = 3.00 Å is unrelated to the error at r = 3.05 Å. | Neighbouring points of a smooth curve, computed from the same frames, move together. Counting them as separate evidence overstates the information. |
| (iii) Homoskedastic noise | "Same scatter everywhere": one noise level \(n_k\) for every point of the curve. | At short range g(r) is exactly 0 because atoms cannot overlap, so there is no noise there; the noise is largest around the peaks. The paper names this case. |
| (iv) Gaussian residuals | Errors follow a bell curve. | Reasonable for averages over many frames (central limit theorem); weaker for small counts or quantities near a hard bound. |
| (v) \(n^{\text{obs}}_k = 1\) | Each whole curve counts as one observation. | A curve with three clear hydration shells carries more information than a single number, yet gets the same weight. |
A nuance the paper stresses: correlations between QoIs that come from the charges are still captured. Changing one charge moves both the RDF and the H-bond count, and both predictions depend on the same \(\theta\). The independence assumption is only about the noise.
Why the observation count matters most
The noise level \(n_k\) is learned together with the charges. If you let it settle at its best value for a given \(\theta\), each QoI's term collapses to something simple, written with the mean squared error (MSE), the average of the squared mismatches over the curve:
Read aloud: "The maximum over n k of log L k of theta equals minus n obs k over two, times log M-S-E k of theta, plus a constant." \(L_k\) ("L k") is QoI \(k\)'s likelihood term; "max over n k" means "with n k set to its best value".
Two consequences: the units drop out, so halving the error of any QoI is worth the same whether it is an RDF or a count; and the only thing that sets how hard a QoI pulls is its observation count.
| Observations per curve | Effect |
|---|---|
| 1 (the paper) | Every QoI weighs the same. Can be too cautious: error bars wider than the data justify. |
| Number of grid points (e.g. 200) | Treats every point as independent evidence. Overconfident: error bars too narrow, because neighbouring points are strongly correlated. |
| Effective count \(n^{\text{eff}}\) (the code now) | A middle ground: roughly how many independent features the curve really has. |
Rule of thumb: when the data dominate, posterior width shrinks roughly like \(1/\sqrt{n}\). Going from 1 to 200 observations makes the error bars about 14 times narrower, which is why this choice matters.
Paper versus code: what changed and what is still open
Here is the likelihood as released with the paper (tag v0.0.1, bff/bayes/likelihoods.py, trimmed; comments are mine):
N = model.observations # 1 for every QoI in the paper
diff = y_true[qoi] - y_trial
ssq = torch.sum(diff**2, dim=1) # sum of squares over the grid
log_like_valid += -0.5 * ssq / sigma_exp**2 - N * torch.log(sigma_exp)That line is Eq. 7 exactly: \(\sigma\) is \(n_k\), ssq is \(\lVert \hat{y}_k - y_k \rVert^2\) and N is \(n^{\text{obs}}_k\). In the current version (v0.4.1), gaussian_log_likelihood_by_qoi computes:
diff = problem.observations[qoi] - y_trial
mse = torch.mean(diff**2, dim=1)
n_eff = float(model.n_eff)
contribution[mask] = (
-0.5 * n_eff * mse / sigma**2
- n_eff * torch.log(sigma)
)Read aloud: "log L k equals minus n eff k over two, times M-S-E k of theta over sigma k squared, minus n eff k times log sigma k." The code's \(\sigma_k\) ("sigma k") is the paper's \(n_k\).
- Mean instead of sum. It takes the mean squared error over the grid and then multiplies by \(n^{\text{eff}}\). Sampling a curve more finely therefore doesn't change its weight. There is a test for exactly this:
test_gaussian_log_likelihood_is_invariant_to_duplicated_bins. - One noise level per QoI, \(\sigma_k\), either learned (sampled on a log scale) or fixed from the dataset. So (ii) independent residuals and (iii) homoskedastic noise still hold.
- Terms are added across QoIs, so (i) independence between QoIs still holds.
- Invalid charges, outside the bounds or breaking the total-charge rule, get a log-likelihood of \(-\infty\), i.e. zero probability.
| Simplification | Paper code (v0.0.1) | Code now (v0.4.1) | Status |
|---|---|---|---|
| One observation per curve | N = model.observations, set to 1 | n_eff set by hand, as the number of grid points, or by peak counting (effective_observations.py, chosen in workflows/learn/main.py) | Partly solved. \(n^{\text{eff}}\) is a heuristic, not derived from the measured noise. |
| Weight depends on grid spacing | Sum of squares: with a fixed noise level, a finer grid means more weight | Mean of squares times \(n^{\text{eff}}\); tested by test_gaussian_log_likelihood_is_invariant_to_duplicated_bins | Solved. |
| (iii) Same noise all along the curve | One sigma per QoI | One sigma per QoI ("one scalar nuisance per QoI" in _split_parameters_and_sigmas) | Open. |
| (ii) Independent residuals along a curve | Plain sum of squares (no covariance) | Plain mean of squares (no covariance) | Open. \(n^{\text{eff}}\) only compensates roughly. |
| (i) Independent noise between QoIs | Terms added over QoIs | Terms added over QoIs in gaussian_log_likelihood. They are now also kept per QoI, which powers the QoI-attributed marginal plots. | Open. |
| Surrogate error folded into the noise | The learned sigma absorbs it | Same. The GP's own predictive variance is not used in the likelihood. | Open. |
| (iv) Gaussian noise | Gaussian | Gaussian | Open, lower priority. |
So the group has started on how much weight each curve gets, but the noise model itself, its shape along the curve and its correlations, is still the paper's. That gap is small enough for one person, sits in one file with existing tests, and has a measurable outcome.
Reading tip: diffs between tags are a quick way into an unfamiliar codebase. In the cloned repo, git diff v0.0.1 v0.4.1 -- bff/bayes/likelihoods.py shows this whole story. It also shows smaller fixes: in v0.0.1, if no charge constraint was passed, mask was never set and the function would fail; v0.4.1 handles that case in _valid_parameter_mask.
The effective count is chosen per QoI in the learn configuration (applied in bff/workflows/learn/main.py), using exactly one of three modes:
models:
rdf:
model_path: ../05-lgp/models/rdf.lgp
tolerance: 0.1 # estimate n_eff from the curve's peaks
pmf:
model_path: ../05-lgp/models/pmf.lgp
n_eff: 5 # set it by hand
# or: independent_observations: true → n_eff = number of grid pointsThe peak-counting estimate (bff/bayes/effective_observations.py) goes like this for each reference curve:
- Smooth the curve with a Savitzky–Golay filter (a cubic fitted in a sliding window of up to 15 points).
- Estimate the noise as the median distance between the raw and the smoothed curve.
- Find peaks and troughs whose prominence, their height above the surroundings, is at least 5 times that noise.
- Each one adds \(\max(1, \ln(\text{prominence} / \text{tolerance}))\), where the tolerance is a scale you set in the curve's units. A very pronounced feature counts a bit more than one. The total is at least 1.
It's a sensible heuristic, and the part of the code closest to your EEG signal-processing work. It does not use how the reference noise actually behaves, which is where a project could start. The package also draws QoI-attributed marginals (plot_qoi_marginals in bff/plotting.py): for each region of a charge's posterior, which QoI's likelihood favours it most.
A first intern project: measure the noise instead of assuming it
Goal: estimate the reference noise directly from the reference trajectory, use it in the likelihood, and check whether the learned charges and their error bars change, and whether they are better calibrated.
- Reproduce a baseline. Run the learn stage of the acetate example three times:
n_eff: 1,tolerance, andindependent_observations: true. Compare the marginals. This teaches you the code and shows how much the weighting matters. - Measure the noise. Cut the reference trajectory into B consecutive blocks, each long enough to be roughly independent (block averaging). Compute every QoI in each block. The scatter between blocks at each r gives the noise there. How neighbouring r values move together gives the correlations along the curve. RDF and H-bond values from the same block give the correlation between QoIs. Divide by B to get the uncertainty of the full-trajectory average.
- Handle the small-sample problem. With a handful of blocks and around a hundred grid points, the raw covariance matrix cannot be inverted. Shrinkage, for example Ledoit–Wolf, blends it with a simple diagonal matrix so it can.
- Add a noise-model option next to the current one in
likelihoods.py, keeping today's behaviour as the default: heteroskedastic (a noise level \(\sigma_k(r)\) for each grid point) and full covariance (below). Write tests in the style oftests/bayes/test_likelihoods.py. - Check calibration with a known answer. Pick a charge set, generate "reference" QoIs from it, run the inference, and record whether the 95% intervals contain the true charges. Repeat many times: a calibrated model contains the truth about 95% of the time. Compare the noise models.
- Report marginals and QoI-attribution plots for each noise model, and whether charges that matter downstream move, such as the carboxylate oxygen for troponin.
Read aloud: "log L k equals minus one half r k transpose, capital sigma k inverse, r k, minus one half log det capital sigma k, plus a constant, where r k is y hat k of theta minus y k." \(r_k\) ("r k") is the list of residuals; \(^{\top}\) ("transpose") turns that column into a row so the product gives one number; "det" is the determinant, a single number that measures the overall size of the matrix's noise.
\(\Sigma_k\) is the covariance matrix of the noise. Its diagonal holds the noise variance at each r (heteroskedastic); its off-diagonal entries say how neighbouring points move together. A diagonal \(\Sigma_k\) with equal entries gives back today's likelihood. With real correlations in \(\Sigma_k\), a hand-set \(n^{\text{eff}}\) is no longer needed, because the correlations already say how much independent information the curve holds. A learned scale factor on \(\Sigma_k\) can still absorb surrogate error.
Things to watch out for
Mentioning these shows you've thought it through.
- Short reference runs. Few blocks means the noise estimate is itself uncertain. Shrinkage helps; the longer MLIP reference runs of the newer workflow help more.
- Surrogate error is noise too. The GP reports its own predictive variance at each point, which could be added to \(\Sigma_k\) instead of being folded into one learned number.
- The duplicated-bins property breaks. A copied grid point makes \(\Sigma_k\) singular, so the grid has to be fixed or thinned before building \(\Sigma_k\).
- Cost. Every MCMC step needs a solve with \(\Sigma_k\). Computing its Cholesky factor once, as the code already does for the GP kernel, keeps the extra cost small.
- Prior work. The paper points to ref. 22 (Sullivan, Cervenka, Shanks and Hoepfner, J. Phys. Chem. B 2025), where a physics-informed Gaussian process gives a learned covariance for liquid-structure data. Brennon Shanks is on both papers. Skim it before the interview if you can.
How to pitch it in the interview
Say it in three short parts: what the paper did, what the code does now, and what you would add. Then ask questions, so it becomes a conversation rather than a speech.
1. What the paper did. "In the paper, the likelihood is Eq. 7: for each QoI, the sum of squared differences between the surrogate's prediction and the AIMD reference, divided by one learned noise level. Each whole curve counts as one observation. The paper itself says this assumes the noise is the same all along a curve and independent between points and between QoIs."
2. What the code does now. "The code released with the paper computes exactly that: a sum of squares, with the observation count set to 1. In version 0.4.1 it takes the mean squared error instead and multiplies it by an effective observation count, n_eff, which you can set by hand, set to the number of grid points, or estimate by counting peaks in the reference curve. So a curve's weight no longer depends on the grid spacing, and richer curves can count for more. What hasn't changed is the noise model: still one noise level per QoI, no correlations along a curve or between QoIs, and the surrogate's own uncertainty folded into that noise level."
3. What I would add. "As a first project I'd measure the reference noise directly, by block-averaging the reference trajectory, and use it in the likelihood: first a noise level that varies along the curve, then a full covariance matrix. I'd keep today's behaviour as the default and check calibration on synthetic data with a known answer: do the 95% intervals contain the true charges 95% of the time? That's close to how likelihood-ratio methods are validated in forensic science, which is my background."
Grid points. A curve like an RDF is stored as numbers at fixed distances r, for example every 0.1 Å from 2 to 7 Å, which gives 51 numbers. Those r values are the grid points (also called bins), and the grid spacing is the 0.1 Å gap between them.
estimate_curve_n_eff, tolerance 0.1) finds the four circled features and adds 2.2 + 1.9 + 1.2 + 1.0 ≈ 6.3, so this 51-point curve counts as about 6 observations rather than 51 or 1.Observation count and curve weight. The likelihood has to decide how much each QoI counts. The observation count says how many independent measurements a QoI is treated as. The curve's weight, meaning how hard it pulls on the charges, grows with that count and shrinks with the square of its noise level. The paper uses a count of 1 for every curve.
Why the weight "no longer depends on the grid spacing". The paper adds up the squared errors over all grid points. Store the same curve every 0.05 Å and there are twice as many terms, so with a fixed noise level the curve would count twice as much without any new information. The current code averages instead (mean squared error), so the grid spacing drops out.
Effective observation count, n_eff ("n eff"). A 51-point RDF is not 51 independent facts, because neighbouring points move together. What it really tells you is closer to a handful of features: where the first peak is and how high, where the trough is, and so on. n_eff is an estimate of that number. You can set it by hand; set it to the number of grid points (51 here, which is overconfident); or let the code count peaks.
Peaks in the reference curve. The code smooths the reference RDF, finds the clear peaks (hydration shells) and troughs (the gaps between shells), and adds at least 1 for each, more for large ones. The figure above shows the result: about 6.3.
"Richer curves can count for more". A curve with several sharp shells gets a higher n_eff than a nearly flat one, so it pulls harder on the charges.
The noise model. The assumptions about the leftover mismatch between prediction and reference: how big it is, whether its size changes along r, and whether errors at nearby points, or in different QoIs, move together. The code still uses the simplest version: one size for the whole curve and no correlations.
The surrogate's uncertainty "folded into" the noise level. A Gaussian process predicts a value and also says how unsure it is at each point, for example more unsure for charges near the edge of the training range. The code doesn't use that per-point uncertainty. Instead, the single learned noise level has to cover both the reference's noise and the surrogate's error, averaged over the whole curve. "Folded into" means lumped together into that one number.
A side finding you could ask about. Because the smoothing window is a fixed 15 points, the peak count depends on the grid spacing. The same made-up curve gives n_eff ≈ 6.3 on a 0.1 Å grid and ≈ 9.0 on a 0.05 Å grid (or with every point duplicated). The likelihood itself is built not to depend on grid spacing, so it may be worth asking whether the window should be set in ångströms instead. This comes from my own quick test on a made-up curve, not from their data.
Questions to ask while you pitch
Unfamiliar word? See the terms used in this section.
- "How did you choose the peak-counting rule for n_eff, and have you compared it with a noise estimate taken from the reference data?"
- "How much do the posteriors change between n_eff = 1, the peak-counting estimate, and one observation per grid point?" (n_eff, "n eff", is the effective number of observations: how many independent pieces of information a curve is counted as. n_eff = 1 is the paper's setting.)
- "In a quick test on a made-up RDF, the peak-counting n_eff changed with grid spacing, from about 6.3 at 0.1 Å to about 9.0 at 0.05 Å, because the smoothing window is 15 points. Is that intended, or should the window be set in ångströms?"
- "Are the reference runs from the MLIP long enough to block-average? Roughly how many independent blocks would a typical run give?"
- "Would you keep one learned noise scale per QoI on top of a measured covariance, so surrogate error is still absorbed?"
- "Ref. 22 learns a covariance for structure data with a physics-informed GP. Is plugging that into this likelihood already planned?"
- "Is anyone working on this now? If so, which part could I take on?"
Follow-ups they might ask, with short answers
- "Why not just set n_eff to the number of grid points?" Neighbouring points of an RDF are strongly correlated, so they aren't independent evidence. Counting them all would make the error bars far too narrow.
- "Won't a full covariance be slow inside MCMC?" The covariance stays fixed during sampling, so its Cholesky factor is computed once. Each step then needs one triangular solve per QoI, which is small next to the surrogate prediction.
- "How would you know it's better?" Two checks: calibration on synthetic data with known charges, and the usual validation of rerunning real MD at posterior samples and comparing with the reference.
- "What if the reference run is too short for a good covariance?" Use shrinkage toward a simple diagonal matrix, or a smooth parametric form, for example a noise level that follows the RDF plus a correlation length along r, fitted to the few blocks available.
This plan comes from reading the code at v0.4.1, not from running it. Ask the group what is already in progress.
Likely questions, with answers you can adapt
Answers are written in plain first person. Rephrase them in your own words rather than memorizing.
Explain the paper in two minutes.
Problem: fixed-charge force fields need partial charges, which aren't measurable and aren't unique, and the usual fitting gives one set with no error bars.
Method: they learn a posterior over charges using AIMD in explicit water as reference. Classical MD runs at thousands of trial charge sets train Gaussian-process surrogates, and MCMC samples the posterior through the surrogates.
Results: 18 biological fragments, better agreement with AIMD than CHARMM36-nbfix, densities within 1% of experiment, chemically sensible error bars, and a troponin Ca²⁺ binding free energy near experiment.
Why it matters: you get uncertainty, you know which charges you can change safely, and you can propagate that uncertainty into predictions like binding free energies.
Uncertainty: each charge comes with an error bar, so you know how well the data pin it down.
Charges you can change safely: a charge with a wide error bar can move within that range without making the match to AIMD worse. That matters when fragments are joined into a protein: the pieces must add up to a whole-number total charge, so a little charge has to be added or removed somewhere. The paper puts that correction (up to ±0.1 e) on the atom with the widest posterior, usually a buried carbon, because changing it barely affects anything. Changing a tightly pinned charge, like an exposed oxygen, would spoil the hydration structure.
Propagating uncertainty into predictions: instead of computing the calcium binding free energy once, with one "best" charge set, they compute it with several charge sets drawn from the posterior (the MAP plus 10 samples). The spread of the results shows how much the prediction depends on charges the data can't fully pin down. Here the results fell within 10 kJ/mol of experiment, the best within about 1, while the single best-fit set was off by about 6. It is the same idea as reporting a lab measurement with its uncertainty, like a blood-alcohol result given with an interval instead of a single number.
Why Bayesian inference instead of minimizing a loss?
Because many charge sets fit about equally well, a single optimum hides that. The posterior shows which charges are pinned down and which are flexible. That is practical: you know where to put extra charge when joining fragments, and you can run a downstream calculation over several plausible parameter sets.
Two more points. Priors regularize when reference data are scarce and expensive. And the nuisance noise parameters learn how much weight each QoI gets, which replaces hand-tuned loss weights.
Why use AIMD as the reference instead of experiments?
There are far more parameters than experimental observables, and charges are not observable at all. AIMD gives detailed, atom-resolved structure (an RDF for each solute atom) for each fragment in explicit water, including polarization, without ad hoc gas-to-liquid corrections. Experiments are then used as independent checks: solution densities and the troponin binding free energy. The authors also say the framework can take experimental data, such as scattering patterns, directly in the likelihood.
Why do you need a surrogate model at all?
MCMC needs 10⁴ to 10⁵ likelihood evaluations, and each would need an MD simulation. The surrogate turns that into a matrix multiply on the GPU. You pay the simulation cost once, up front, and those training runs are independent, so they parallelize perfectly on a cluster.
Why a local Gaussian process rather than a neural network or a full GP?
The data set is small (thousands of runs), the map from charges to structure is smooth, and you want something cheap and well behaved with interpretable hyperparameters. The length scales directly tell you how sensitive each output is to each charge.
The surrogate has to predict a whole curve, say an RDF at 100 radii, from one charge set.
A full GP models the whole curve at once, including how the value at one radius relates to the value at the next. Its covariance matrix has one row and one column for every combination of training run and grid point. With 2,000 training runs and 100 grid points that is 200,000 rows and columns: 4 × 10¹⁰ numbers, about 320 GB of memory, and it has to be inverted. Done directly, that is not practical. (Tricks exist that exploit the matrix's structure, but they add complexity.)
A local GP treats each grid point as its own small question: "given the charges, what is g(r) at this one radius?" All grid points of one QoI share the same kernel settings, so there is only one 2,000 × 2,000 matrix (4 million numbers, about 32 MB) to invert, once. Predicting the whole curve is then one matrix multiplication. "Local" means one output point at a time, not local in charge space.
What you give up: the local GP does not model how errors at neighbouring radii are correlated. That is fine for the mean prediction, which is what the likelihood uses, but it is one more reason the likelihood's noise model is a simplification (see the likelihood deep dive).
A neural network would typically need more data and give less trustworthy uncertainty in this regime.
What do the nuisance parameters do?
One per QoI. Each is a noise standard deviation: how far the prediction may miss the reference even when the charges are right. It covers two sources:
- Surrogate error. The Gaussian process is an approximation. For a given charge set, its predicted RDF differs a little from what a real MD run with those charges would give. It was trained on a finite number of runs, and each training run is itself a finite simulation with some statistical noise.
- Noise in the reference. The AIMD reference is an average over a short simulation of a small box (128 water molecules). Run it again from a different starting point and you would get a slightly different RDF. That wobble is noise in the target itself.
Because of both, even perfect charges would not match the reference exactly, and \(n_k\) is the expected size of that leftover mismatch. It is sampled along with the charges and then ignored (marginalized) when reading the charge posterior. Its practical effect is automatic weighting: a QoI that can't be fitted well gets a larger noise scale and less influence.
How do you know the MCMC has converged?
The paper tracks the integrated autocorrelation time τ, requires the chain to be at least 100τ long with the τ estimate stable to 1%, discards 2τ as burn-in and thins. The newer code uses rank-normalized split-R̂, which compares chains against each other.
MCMC produces its samples one after another. A walker starts at some charge set, proposes a nearby one, accepts or rejects it, and repeats. The list of charge sets it visits, step by step, is its chain. emcee runs many walkers side by side (5 per parameter here), each with its own chain, and the posterior is estimated from all the visited points together.
Each step starts from the previous one, so consecutive points are similar. The autocorrelation time τ ("tau") is roughly how many steps the chain needs to "forget" where it was, so a chain of length L holds about L/τ independent samples. The start of a chain still reflects the random starting point and is thrown away (burn-in). Thinning keeps only every k-th point, since close neighbours add little information.
R̂ ("R-hat") compares several chains. If they have all converged, each chain looks like the others: the spread within a chain matches the spread between chains, and R̂ is close to 1. "Rank-normalized split" is a more robust version that splits each chain in half and works with ranks instead of raw values.
The stronger check is external: draw charge sets from the posterior, run real MD without the surrogate, and confirm they reproduce AIMD. That tests the surrogate and sampler together.
What is ECC and why scale charges by 0.8?
Fixed charges can't polarize, so charged groups in water interact too strongly. ECC treats electronic polarization as a uniform background dielectric and scales ionic charges by about 1/√εel, roughly 0.75 for water. The paper uses 0.8 from the group's recent work, with matching water and ion models. The CHARMM36 result of −105.7 kJ/mol against −28.6 experimental shows what happens without it.
CHARMM36-nbfix got −31.7 kJ/mol, close to experiment. Why isn't that good enough?
It is right for the wrong reason. NBFIX adds pair-specific Lennard-Jones corrections that push ions apart to compensate for electrostatics that are too strong. Two errors cancel for this property, but there is no guarantee they cancel for other properties or other ion pairs. From forensics: a correct conclusion from an invalid method doesn't validate the method.
Why did the MAP charge set underbind calcium?
The MAP is the best set for acetate in bulk water, but a carboxymethyl side chain inside a binding pocket is a different chemical environment. Fixed charges can't adapt to that. The posterior still contains sets that match experiment within about 1 kJ/mol, and the spread tells you how much the prediction depends on the charges. That is the argument for keeping the full distribution.
What are the main limitations?
Scaling to many parameters (shown up to 10), dependence on modelling choices in the likelihood and surrogate, reliance on the quality of one DFT reference, simplified noise assumptions, and learning only charges. See the limitations section. It's good to mention that the newer code already addresses some of these (effective observations, LJ parameters, MLIP references).
How would you scale this to 30 or more parameters?
Why it's hard: the surrogate needs training runs that cover the space of possible charges, and that space grows explosively with the number of parameters. Trying just 3 values per charge in every combination takes 3¹⁰ ≈ 59,000 runs for 10 charges, but 3³⁰ ≈ 2 × 10¹⁴ for 30. Latin hypercube sampling needs far fewer runs than a full grid, but the same "curse of dimensionality" applies. MCMC also gets slower as dimensions are added.
Five ideas, each attacking a different part of the problem:
- Active learning: spend simulations where they matter. Most combinations of charges are implausible, and a first rough posterior shows which ones. Train on a modest first batch, run MCMC, then run new simulations only where the posterior is high, retrain, and repeat. Like searching a field with a metal detector: after the first beep you stop walking the whole field.
- Sensitivity screening: drop the charges that don't matter. The GP's length scales say how much each output changes with each charge. A very long length scale means the outputs barely notice that charge, like the buried carbon in the toy demo. Fix those at a sensible value and sample only the rest.
- Parameter sharing: fewer independent numbers. Chemically equivalent atoms (the two oxygens of a carboxylate, the three hydrogens of a methyl group) already share one charge through their atom type. A step further is a hierarchical prior: the same kind of atom in different molecules draws its charge from one shared distribution, so what is learned about the carboxylate oxygen in acetate also informs the one in formate.
- Gradient-based sampling: smarter moves. emcee's stretch move only looks at score values, and in many dimensions most random proposals land somewhere worse and get rejected. Hamiltonian Monte Carlo (HMC) also uses the slope of the score, which PyTorch can compute because the GP is differentiable, to make long, well-aimed moves.
- Other surrogates for bigger data. An exact GP's memory grows with the square of the number of training runs and its compute with the cube, so at around 100,000 runs it becomes impractical. Sparse GPs summarize the data with a few hundred representative "inducing points". Neural networks scale to large data sets, and training several of them (an ensemble) and comparing their predictions gives an uncertainty estimate.
How would you add Lennard-Jones parameters?
What Lennard-Jones (LJ) parameters are. Each atom type gets two numbers that describe its non-electric interactions, the fifth term in Background, idea 2 (with the curve):
- σ ("sigma"): the atom's size, roughly the distance at which two atoms start pushing each other away. In GROMACS files it is in nanometres.
- ε ("epsilon"): how sticky the atom is, the depth of the weak attraction well, as an energy in kJ/mol.
For a pair of different atoms, CHARMM combines the two: the pair's size is the average of the two sizes, and its ε is the geometric mean (the square root of the product).
Adding them to the learning. Put σ and ε into θ ("theta", the list of parameters being learned) next to the charges q, and keep the same pipeline: sample them, simulate, train the surrogate, run MCMC.
The catch is identifiability: whether the data can tell the parameters apart. A bigger σ keeps two atoms further apart, which weakens their electrostatic attraction at contact, much like a smaller charge would. So many (charge, σ) combinations fit the data about equally well. The posterior then shows a long, tilted ridge instead of a compact peak, and neither parameter is pinned down on its own. The fix is to add QoIs that respond differently to each, like densities, ion-pair distributions and hydration structure together. The code now supports LJ parameters, so this is a good question to ask them back.
What would you change or test in the paper?
Pick one and explain it well. For example: the likelihood counts each RDF as one observation, which might under- or over-weight curves with many features. The newer code's effective observation counts address exactly that, and I'd like to see how much the posteriors change. Another: add a model-discrepancy term for the DFT reference, so the posterior also reflects that revPBE-D3 is not the truth.
What would you like to work on first?
Use the pitch from the likelihood deep dive: measure the reference noise by block averaging, try heteroskedastic and full-covariance likelihoods, and test calibration on synthetic data with a known answer.
Why it suits a first project: it's self-contained (one module plus tests), statistical rather than chemical, has a measurable outcome, and addresses a limitation the paper names itself.
Show that you know where the code is today, so the idea builds on their work instead of repeating it. Something like: "I know the current code already changed the paper's likelihood, by averaging the squared errors and adding an effective observation count. What I'd add is the part that's still the paper's: the noise model itself."
The paper's code is tag v0.0.1; the current release is v0.4.1. In the likelihood, the part this project touches (bff/bayes/likelihoods.py):
- Mean instead of sum. The paper adds up the squared errors over the grid; the code takes their average (mean squared error). A curve's weight no longer depends on how finely it is stored.
- Effective observation count n_eff instead of 1. Each QoI's term is multiplied by n_eff, set by hand, set to the number of grid points, or estimated by counting peaks in the reference curve (v0.3.0).
- Terms kept per QoI. The code keeps each QoI's contribution separately, which feeds the new QoI-attributed marginal plots (which observable supports which charge values).
- A small fix. In v0.0.1 the function failed if no charge constraint was given; v0.4.1 handles that case.
Still as in the paper, and so open for your project: one noise level per QoI, constant along the curve; no correlations along a curve or between QoIs; the GP's own uncertainty not used; Gaussian noise. Full side-by-side code in paper versus code.
Elsewhere in the pipeline (details in what changed in the code): reference MD can now come from a machine-learned interatomic potential (MLIP) instead of AIMD; Lennard-Jones σ and ε and dihedral force constants can be learned, not just charges; charge constraints can be set per residue and per system, for building proteins; and the code has its own PyTorch MCMC sampler that checks convergence with R̂ ("R-hat") instead of using emcee with autocorrelation times.
Why this group, and what would you bring?
Keep it true and specific. Something like: Bayesian reasoning is what I used in forensic evidence evaluation, and this paper applies it to a problem where it clearly helps. I'm new to molecular simulation and I'll be learning GROMACS and the chemistry, but I can contribute from day one on the statistics (GPs, MCMC, validation) and on the software side (Python and PyTorch pipelines, Linux, cluster jobs, testing). I've also done signal processing on EEG and fNIRS data, which is close to working with noisy RDF curves.
Questions to ask them
- Now that reference MD can come from an MLIP, how do you check that MLIP errors don't leak into the learned charges?
- How much did the effective observation counts change the posteriors compared with one observation per curve?
- When you learn charges and Lennard-Jones parameters together, do strong correlations show up, and how do you choose QoIs to separate them?
- Would you add experimental data, such as scattering patterns or densities, into the likelihood alongside AIMD?
- What biological system comes after troponin?
- Is anyone already working on the noise model in the likelihood, for example heteroskedastic or correlated noise? I'd be interested in that as a first project.
- What does a typical intern project look like here, and which part of the pipeline would I work on?
Flashcards and quiz
Facts and numbers worth having ready. Cards you mark as known are skipped until you reset. Progress is saved in this browser only.
Multiple-choice quiz
Flashcards test recall; these test whether you can tell the right idea from a near miss. Each question has four options of similar length, and every wrong option is a mistake that is easy to make after one read of the paper: a number from elsewhere in the paper, a cause and effect reversed, or something the code does now but the paper didn't. The options are shuffled every time, so position and length give nothing away. After each answer you see why the tempting options are wrong. Progress is saved in this browser only.
Glossary
Symbols and how to say them
| Symbol | Say it | Meaning on this page |
|---|---|---|
| p(a | b) | "p of a given b" | Probability of a if b is true. The bar is read "given". |
| ŷ | "y hat" | A predicted value (the hat means "estimated"); y without a hat is the reference. |
| q* | "q star" | The best or combined value of a charge in a worked example. |
| θ | "theta" | The list of parameters being learned (here the charges). |
| ϑ | "theta" (curly form) | A bond angle, in the force-field formula. |
| q | "q" | A partial charge. |
| e | "e" | The elementary charge, the charge of one proton (a unit of charge, not energy). Charges are given in units of e. |
| σ | "sigma" | Lennard-Jones size. In the code, also the noise level of a QoI. |
| ε | "epsilon" | Lennard-Jones well depth (stickiness). |
| ε₀ | "epsilon nought" | The vacuum permittivity, a physical constant in the Coulomb term. |
| ε_el | "epsilon e-l" | Electronic dielectric constant of water (about 1.78), used by ECC. |
| φ, δ | "phi", "delta" | Dihedral (twist) angle and its phase. |
| n_k | "n k" or "n sub k" | Noise level of QoI k in the paper (a nuisance parameter). |
| n_obs, n_eff | "n obs", "n eff" | Number of observations; effective number of observations. |
| Σ | "capital sigma" | A sum sign. With a subscript (Σ_k), a covariance matrix. |
| α, l | "alpha", "ell" | GP kernel scale and length scale. |
| τ | "tau" | Integrated autocorrelation time of an MCMC chain. |
| R̂ | "R-hat" | Convergence check that compares chains. |
| λ | "lambda" | The on/off switch in alchemical free-energy calculations. |
| ΔG | "delta G" | Free energy change, such as the binding free energy. |
| ∂ | "partial" or "del" | Partial derivative, as in F = −∂U/∂r. |
| ‖x‖² | "norm of x, squared" | Sum of the squares of the entries of x. |
| ‖x‖₁ | "one-norm of x" | Sum of the absolute values of the entries of x (used in NMAE). |
| ∏ | "product" or "capital pi" | A multiplication sign over many terms, like Σ for a sum. |
| π | "pi" | The number 3.14159… |
| exp(x) | "exp of x" or "e to the x" | The exponential function (this e is 2.718…, not the elementary charge). |
| I | "I" or "the identity" | The identity matrix: 1 on the diagonal, 0 elsewhere. |
| x⊤, A−1, det A | "x transpose", "A inverse", "det A" | Row version of a column; matrix inverse; determinant (one number for a matrix's overall size). |
| Å | "ångström" (ONG-strum) | Length unit, 0.1 nm. |