Chapter 01

The fly and its connectome

A connectome is a wiring diagram of a nervous system at the resolution of individual synapses: a list of neurons, and for every pair, the number of synaptic contacts one makes onto the other. It is a static anatomical object. It says where the wires run and how many contacts they make; it does not say what passes along them, or when, or with what sign.

The MaleCNS v1.0 central nervous system, rendered by the dataset's authors (FlyEM / HHMI Janelia, University of Cambridge, MRC LMB, Google Research; male-cns.janelia.org, CC BY 4.0). Forager takes from this reconstruction only the directed graph of neurons and their connections.

How a connectome is made

The tissue is fixed, stained with heavy metals and cut or milled into sections a few tens of nanometres thick, each of which is imaged by electron microscopy. The images are aligned into a single three-dimensional volume at a resolution of a few nanometres per voxel. Membranes are visible as dark lines, synapses as dense presynaptic specialisations with clustered vesicles opposite a postsynaptic density.

Segmentation follows: a convolutional network assigns each voxel to a cell, producing millions of fragments that are agglomerated into candidate neurons. Synapse detection is a second network that marks presynaptic sites and pairs each with its postsynaptic partners. Both steps make errors, and the largest share of the labour in any modern connectome is proofreading, in which people inspect merged or split fragments, correct them, and attach names and cell types to the corrected cells. The output is a table with one row per directed neuron pair and a count of the synaptic contacts between them.

MaleCNS v1.0

MaleCNS v1.0 is the connectome of the central nervous system of one adult male Drosophila melanogaster: the brain, including both optic lobes, and the ventral nerve cord, with the neck connective intact between them. It was produced by the FlyEM project team at HHMI Janelia Research Campus together with the University of Cambridge, the MRC Laboratory of Molecular Biology and Google Research, and released under a CC BY 4.0 licence. The reconstruction is fully proofread and annotated, with cell types assigned across the whole nervous system.

Forager uses the neurons that carry a type annotation and the connections between them that have at least five synaptic contacts, reduced to the largest weakly connected component. That graph has 163,972 neurons and 6,143,838 directed connections. The threshold of five contacts is a conventional cut in fly connectomics below which a connection may be a segmentation or detection artefact; the weakly connected component is taken so that no neuron in the model is unreachable from every other. A smaller graph of the 16,384 best-connected neurons was also used, and is described in chapter 3.

Why a fly

The fly is the standard animal for the cellular study of learning. Its brain has roughly 140,000 neurons, small enough to be reconstructed in full and large enough to support behaviour that can be scored. Its genetics allow any identified cell type to be silenced, activated or recorded in a living animal, and the same cell types can be found from one fly to the next. A hungry fly that smells an odour while receiving sugar will later approach that odour; one that smells it while being shocked will later avoid it. This associative olfactory learning is fast, reproducible and quantifiable, and the circuit that implements it, the mushroom body, is one of the best-characterised learning circuits in any animal. Chapter 2 describes it.

That circuit is also where the idea of a reward-modulated local synaptic rule comes from. The connectome supplies the wiring on which such a rule could act; the mushroom body literature supplies the form of the rule.

What Forager keeps and what it discards

From the reconstruction Forager keeps one thing: who connects to whom, and with how many contacts. The contact counts travel with the graph as anatomical evidence, but they were not used to set the initial synaptic weights in the runs reported here; those weights were drawn at random, with random signs, as chapter 3 explains.

Everything else is discarded. The transmitter released at each synapse, and hence whether it excites or inhibits, is not used. The shape of each neuron, the position of each synapse on it and the distances signals must travel are not used. Time is not represented beyond a few discrete settling steps. There is no membrane biophysics, no spiking, no adaptation. Neuromodulators are reduced to a single scalar broadcast to every synapse at once. Glia, which in the fly ensheath and partition the neuropils, are absent. The optic lobes and nerve cord are present only as nodes and edges; no photoreceptor or motor neuron is driven by anything resembling light or a leg.

This is not a simulation of a fly; it is a recurrent network whose sparsity pattern happens to be the measured graph of a fly, and any similarity to the animal's behaviour would have to be demonstrated, not assumed.

What remains is still an unusual object. Six million directed connections among 164,000 neurons, with the degree distribution, the modularity and the recurrent loops of a real nervous system, is a sparse matrix that no simple generative model reproduces. Whether that structure matters for a learning rule placed on top of it is the question the rest of this book asks, and chapter 10 answers for the three runs that were made.

Chapter 02

The mushroom body and dopamine

The mushroom body is the fly's associative memory centre. It attaches a value, good or bad, to a pattern of sensory input, with a circuit simple enough to draw on one page. The learning rule Forager places on the whole connectome is an abstraction of what happens at one class of synapse in this circuit.

The odour code

Odours arrive at the antennal lobe, where the axons of receptor neurons sort into about fifty glomeruli, one per receptor type. Projection neurons carry the glomerular signals forward to the calyx of the mushroom body, where they synapse onto the Kenyon cells. There are about 2,000 Kenyon cells per hemisphere, so the code is expanded roughly forty-fold on the way in.

Each Kenyon cell has on average six or seven dendritic claws, each wrapped around a single projection-neuron bouton, and the glomeruli it samples appear to be a random draw. A Kenyon cell fires only when more than half of its claws are active at once, and the anterior paired lateral neuron, a single GABAergic cell per hemisphere that receives input from all Kenyon cells and inhibits them all, raises the threshold further as the population becomes active. The result is a sparse code: a given odour drives a few per cent of the Kenyon cells, and different odours drive nearly non-overlapping sets.

Output neurons and compartments

Kenyon cell axons run in parallel bundles that form the lobes of the mushroom body. Along those lobes lie the dendrites of the mushroom body output neurons, the MBONs. Aso and colleagues counted 34 MBONs of 21 types per hemisphere in 2014, and found that their dendrites tile the lobes into 15 compartments, each compartment the private territory of one or a few MBON types. Every MBON therefore sums the same sparse Kenyon cell code, but through its own set of synapses in its own stretch of the axon bundle. The MBONs are few, and they project to premotor and neuromodulatory regions; collectively their firing biases the fly towards approach or avoidance.

Dopamine

The same compartments are innervated by dopaminergic neurons, about 130 per hemisphere in some 20 types, from two clusters. Cells of the PAM cluster mostly project to the medial lobes and are activated by reward such as sugar; cells of the PPL1 cluster mostly project to the vertical lobes and are activated by punishment such as electric shock. Each dopaminergic type sends its axons to one, or at most two, compartments, and the boundaries of its terminals coincide with the boundaries of the MBON dendrites there. Compartment by compartment, one MBON reads out the Kenyon cell code and one dopaminergic input tells it whether the current moment is good or bad.

Learning happens at the Kenyon cell to MBON synapse and requires three things to coincide: presynaptic activity in a Kenyon cell, meaning an odour is present; the postsynaptic MBON, which sets the compartment; and dopamine released into that compartment. When they coincide, the synapses from the active Kenyon cells onto that MBON are depressed. After a fly has smelt an odour while being shocked, the MBONs in the punishment-innervated compartments respond less to that odour, the balance of approach and avoidance outputs tips, and the fly avoids it. The rule is local in space, since each synapse uses only its own two neurons and the dopamine at its own location, and it is a three-factor rule, since no two of the three factors are enough.

The mushroom body circuit in outline: projection neurons drive a sparse population of Kenyon cells; Kenyon cells converge on a mushroom body output neuron; a dopaminergic neuron innervating that compartment gates the plasticity of the Kenyon-cell-to-output synapses. The three factors of the rule are the presynaptic Kenyon cell, the postsynaptic output neuron and the dopamine signal.

The dopamine-like signal in Forager

Forager borrows the third factor and nothing else of this anatomy. After the evaluator returns a value $$E_t$$ in eV per atom at step $$t$$ , the controller forms a quality

\[q_t=\tanh\!\left[\frac{E_0-E_t}{s}\right],\]

where $$E_0$$ is the reference value of the current fidelity stage and $$s=0.1$$ eV/atom is a scale. For a mixing quantity $$E_0=0$$ eV/atom, so a negative mixing enthalpy gives a positive quality, and the hyperbolic tangent saturates the signal at $\pm1$ once the answer is a few tenths of an electronvolt per atom from the reference. Each stage, the machine-learned screen and the density-functional calculation, has its own $$E_0$$ and its own baseline, and the two are never mixed.

The quality is centred on a running baseline $$b_t$$ ,

\[\delta_t=q_t-b_t,\qquad b_{t+1}=b_t+\beta\,\delta_t,\]

with $\beta=0.1$ . The baseline is an exponential moving average of past qualities with a time constant of about ten evaluations, so $\delta_t$ is positive when the latest answer is better than the recent run of answers and negative when it is worse. In the synaptic update of chapter 4 this quantity is written $$r-b$$ , and it multiplies the local, two-neuron term at every synapse. It is the same number at all 6,143,838 synapses in the same update.

What the signal is not

It is not a temporal-difference error. A TD critic maintains a value function over states and compares the reward received with the reward predicted, so that its error signal anticipates future outcomes. Forager's $\delta_t$ has no state, predicts nothing, and compares the present answer only with a smoothed average of past answers. It is closer to the reward-minus-baseline term of a simple policy-gradient estimator than to anything in reinforcement learning with a critic.

It is not a concentration. Dopamine in the fly is released in a particular volume, diffuses, and is cleared on a timescale that itself shapes learning. The Forager signal is a dimensionless number between $$-2$$ and $$2$$ that exists for one update and is then replaced.

It is not compartment-specific, and this is the largest departure from the biology. The fly delivers dopamine to one of 15 compartments, chosen by which dopaminergic neuron fired, so that the third factor already carries information about which synapses should change. Forager broadcasts one scalar to six million synapses at once, on a graph where the notion of a compartment does not exist. Every synapse in the network learns that the whole answer was better or worse than usual; none learns whether it contributed. Whether the two local factors, presynaptic activity and postsynaptic loss gradient, can supply that missing information is the question of chapter 4, and the outcome in these runs, reported in chapter 10, is that they did not.

Chapter 03

From wiring diagram to dynamical system

A wiring diagram is a list of pairs. To make it compute, each pair must become a number that multiplies a signal, each neuron must become a rule for combining the numbers that reach it, and the whole must be driven by something and read out somewhere. This chapter fixes those choices.

The weight matrix

The 6,143,838 directed connections of the graph become a sparse matrix $$W$$ with $$N=163{,}972$$ rows and columns. The entry $w_{ij}$ sits in row $$i$$ , column $$j$$ , and is the weight of the connection from source neuron $$j$$ to destination neuron $$i$$ ; destination rows and source columns, so that $$Wh$$ is the vector of recurrent inputs. Self-connections are removed. Where the connection table lists the same ordered pair more than once, the rows are aggregated into a single edge. The sparsity pattern of $$W$$ is thereafter fixed: no connection is created or deleted by learning, only rescaled.

The values are initialised at random, with random signs. The synapse counts that accompany each connection are anatomical evidence of how many contacts two cells make; they are not conductances, and they were not used to initialise the weights in these runs. Nor were the transmitter annotations used to set signs. What the fly contributes is which of the $N^2\approx2.7\times10^{10}$ possible entries are allowed to be non-zero, roughly one in four thousand.

The rate neuron

Each neuron carries a scalar activity $h_i\in(-1,1)$ . Given a candidate description $$x$$ , the network settles for three steps,

\[h_i^{(k+1)}=\tanh\Big(b_i+(Bx)_i+\sum_{j\to i}w_{ij}\,h_j^{(k)}\Big),\qquad k=0,1,2,\]

where $$b_i$$ is a per-neuron bias, $$(Bx)_i$$ is the input drive described below, and the sum runs over the neurons $$j$$ that synapse onto $$i$$ . The initial state $h^{(0)}$ is not zero but the state retained from the previous decision, so a residue of the last candidate's activity is present when the next is scored.

One model neuron. It receives a fixed random projection of the candidate's description, $(Bx)_i$, and recurrent input $\sum_j w_{ij}h_j$ from the neurons that synapse onto it. Its goodness is $h_i^2$. A synapse changes using only the two neurons it joins and one broadcast scalar, $r-b$.

Three steps is a small number and it should be read literally. Information about the drive at neuron $$j$$ reaches neuron $$i$$ only if there is a path of at most three synapses from $$j$$ to $$i$$ , and the contribution of a path of length three is attenuated by three factors of $$w$$ and three passes through the tanh. The state after three steps is not a fixed point of the map, and no claim is made that one exists or that it would be reached. It is a snapshot of a transient, taken at the same moment for every candidate.

The transient is kept bounded by a constraint on the incoming weights of every neuron, $\sum_j|w_{ij}|\le\rho$ with $\rho=0.85$ . Since $$|h_j|<1$$ , the recurrent input to any neuron is smaller than $$0.85$$ in magnitude at every step, and the slope of tanh is at most one, so the recurrent part of the map shrinks differences between states to at most $$0.85$$ of their size per step in the maximum norm. That is enough to guarantee the iteration cannot diverge, whatever the drive, and would make it contract to a unique fixed point if it were run to convergence. It is not run to convergence.

The drive

The candidate is described by a vector $$x$$ of twenty numbers. Fifteen are composition and geometry descriptors of the sixteen-atom cell: the cell lengths, the cell angles, the volume per atom, the mean and spread of the atomic number over the sites, statistics of the nearest-neighbour distances, the fraction of neighbouring pairs that are the same element, and the number of sites. To these are added the fraction of the evaluation budget remaining, the reward received at the previous step, two flags marking whether the candidate is being scored for the machine-learned screen or for the density-functional stage, and a constant 1 that gives the projection an offset.

The encoder $$B$$ is a matrix of $163{,}972\times20$ fixed random numbers, drawn once and never trained. Every neuron in the network, from photoreceptor to leg motor neuron, therefore receives its own fixed random linear combination of the twenty descriptors. The choice is deliberate: there is no principled way to route an alloy's cell angles into a fly's olfactory system, so the input is spread everywhere without pretending to anatomy.

The score

After settling, the goodness of neuron $$i$$ is $$h_i^2$$ , and the score of the candidate is the mean goodness over the network,

\[\overline{h^2}=\frac1N\sum_{i=1}^{N}h_i^2 .\]

A mean over 164,000 terms is a law-of-large-numbers quantity. Each $$h_i$$ is a tanh of a random projection of $$x$$ plus a bounded recurrent correction, and the projections are independent across neurons with a common scale set by $\lVert Bx\rVert$ . Averaging $$h_i^2$$ over $$i$$ therefore averages over the randomness in $$B$$ and leaves a smooth function of the drive magnitude, with fluctuations of order $N^{-1/2}\approx0.25\,\%$ . The construction predicts that candidates with longer description vectors score higher, that the recurrent term contributes a nearly candidate-independent offset, and that learning, which can move individual $$h_i$$ but not undo the averaging, has little purchase on differences between candidates. Chapter 10 confirms this empirically.

Scores become a choice through a softmax at temperature $\tau=0.06$ , mixed with uniform exploration $\epsilon=0.15$ over the $$n$$ candidates not yet evaluated,

\[p_c=(1-\epsilon)\,\frac{\exp(\overline{h^2_c}/\tau)}{\sum_{c'}\exp(\overline{h^2_{c'}}/\tau)}+\frac{\epsilon}{n},\]

and one candidate is sampled from $$p$$ . The temperature sets how large a difference in mean goodness counts as a preference: two candidates whose scores differ by $$0.06$$ are chosen in the ratio $$e:1$$ before exploration is mixed in.

Scale

On the full graph one decision means settling the network three times for each of up to 64 candidates, a sparse matrix-vector product with six million non-zeros at every step, and takes about 13 s on a single accelerator. One learning update, described in chapter 4, takes about 2.6 s. A run of fifty decisions is therefore a matter of minutes for the controller, against roughly ten minutes for each of the ten density-functional calculations it steers.

A reduced graph of 16,384 neurons was used for one of the three runs. It keeps the best-connected neurons and the measured edges among them. It is a sub-network of the graph, not a coarse-graining of it: a neuron retained in the reduction has lost the input of every neighbour that was not, and nothing stands in for the missing cells.

Chapter 04

Forward-Forward

Backpropagation trains a network with one forward pass, which computes the activities, and one backward pass, which carries the derivative of a single global loss back through every layer in reverse order. The backward pass needs the forward activities to have been stored, needs to know the exact function each layer computed, and needs a return path whose weights are the transposes of the forward weights. Hinton's Forward-Forward algorithm (December 2022) removes the backward pass. In its place are two forward passes that run in exactly the same way on different data with opposite objectives.

Hinton's algorithm

The positive pass is run on real data. The negative pass is run on negative data, which may be supplied from outside or generated by the network itself. Each layer has its own objective, and it is a local one: the layer's goodness should be high on positive data and low on negative data. The goodness Hinton uses in most of his experiments is the sum of the squared activities of the layer's rectified linear units,

\[G=\sum_i h_i^2 ,\]

taken before the layer is normalised. The layer is trained to classify its input as positive or negative with probability

\[p(\text{positive})=\sigma(G-\theta),\]

where $\sigma$ is the logistic function and $\theta$ is a threshold. The weights of the layer move along the derivative of this per-layer objective, which for a squared-activity goodness is simple to write down. No derivative is passed from one layer to another. Hinton describes the procedure as greedy: each layer solves its own discrimination problem given whatever the layer below sends it.

There is a difficulty with stacking. If the first hidden layer has learned to give real data a long activity vector and negative data a short one, the second layer could tell positive from negative by measuring the length of its input and would need to learn nothing. Forward-Forward therefore normalises the length of each layer's activity vector before it is passed on, using the simplest version of layer normalisation, which divides by the length without first subtracting the mean. The activity vector has a length and an orientation. The length defines the goodness of that layer; only the orientation reaches the next layer, which is forced to find new features in the relative activities.

The negative data has to come from somewhere, and its design decides what the layers learn. In the unsupervised example Hinton builds hybrid images: a smooth random mask with large regions of ones and zeros is used to combine one digit with the mask and a second digit with its complement, giving an image whose short-range correlations are those of real digits and whose long-range correlations are not. In the supervised variant the label is written into the first ten pixels of the image; positive data is an image with its correct label, negative data the same image with an incorrect one, so that anything in the image that does not correlate with the label is ignored. A third source, discussed but not the main path of the paper, is top-down prediction by the network itself.

Forward-Forward in outline. Positive and negative data are each passed forward through the layers; each layer compares its own goodness $\sum_i h_i^2$ with a threshold $\theta$ and adjusts its weights locally. The activity is normalised in length before it reaches the next layer. There is no backward pass.

What it buys and what it costs

The argument for the algorithm is not accuracy. It is that the things a backward pass requires are things a brain, or a piece of cheap analogue hardware, is unlikely to have. There is no weight transport: no second set of connections carrying the transposed weights backwards. No activities are stored for a later phase. If an unknown or stochastic black box is inserted between two layers, backpropagation cannot pass through it without a differentiable model of it, while Forward-Forward is unchanged, because nothing has to pass through it. And because the two passes can be separated in time, sensory data can be pipelined through the layers and learning can proceed while the data streams, without stopping to propagate derivatives.

The costs are equally concrete. Credit assignment is weaker: a layer is told only whether its own goodness was on the right side of its threshold, never what the layers above it would have needed. Somebody has to design the negative data, and the algorithm is only as good as the contrast it is given. On permutation-invariant MNIST Hinton reports 1.36 % test error with four layers of 2,000 units after 60 epochs, against the 1.4 % a comparable backpropagation network reaches in about 20; on CIFAR-10 with two hidden layers his own table gives 41 % test error for Forward-Forward against 37 % for backpropagation. His assessment is that the algorithm is somewhat slower, generalises somewhat less well on the problems he tried, and is unlikely to replace backpropagation where power is not the constraint.

Forager's variant

The circuit is not layered. Its units are the 163,972 neurons of the measured graph and its connections are the 6,143,838 measured directed edges, so the network is recurrent and there is no "next layer" to normalise for. Each neuron $$i$$ settles for three steps under the rate equation of chapter 2, $h_i^{(k+1)}=\tanh\big(z_i^{(k)}\big)$ with pre-activation $z_i^{(k)}=b_i+(Bx)_i+\sum_{j\to i}w_{ij}h_j^{(k)}$ , where $$x$$ is the 20-number candidate vector, $$B$$ the fixed random encoder, $$b_i$$ the bias and $w_{ij}$ the weight of the connection from $$j$$ to $$i$$ .

The goodness is taken per neuron rather than per layer, $$h_i^2$$ , and each neuron carries its own loss,

\[\ell_i=\operatorname{softplus}\!\big[-y\,(h_i^2-\theta)\big],\qquad y=\pm1,\quad \theta=0.18 .\]

Since $\operatorname{softplus}(u)=\log(1+e^{u})=-\log\sigma(-u)$ , this is the negative log of Hinton's probability $\sigma\!\big(y(h_i^2-\theta)\big)$ that the phase was correctly classified, so it is the same objective written for one neuron. Because $h_i\in(-1,1)$ the goodness of one neuron is at most 1, and $\theta=0.18$ is crossed when $$|h_i|$$ exceeds about 0.42.

The phases are not real and fake data. After at least two labelled outcomes exist at a fidelity, the best observed outcome so far is replayed as the positive phase ( $$y=+1$$ ) and the worst as the negative phase ( $$y=-1$$ ), from the stored candidate vectors. "Negative" therefore means ranked lowest among the things that have actually been evaluated at that fidelity; it does not mean physically impossible, and candidates that were never queried are never labelled.

The weight update is

\[\Delta w_{ij}\propto-(r-b)\,\frac{\partial \ell_i}{\partial z_i}\,h_j , \qquad \frac{\partial \ell_i}{\partial z_i}=-2y\,h_i\,(1-h_i^2)\,\sigma\!\big[-y(h_i^2-\theta)\big],\]

where $$z_i$$ is the pre-activation of neuron $$i$$ at the last settling step and $$h_j$$ is the presynaptic activity at that same step, held fixed. The derivative is stopped there: it does not go through the tanh of any other neuron, and it does not go back through the earlier settling steps in which $$h_j$$ was itself computed. Every factor is available at the synapse, the postsynaptic derivative on one side and the presynaptic activity on the other. The third factor $$r-b$$ , the reward minus its running baseline of chapter 2, is a scalar broadcast to every synapse and makes the rule a three-factor rule: the contrast is scaled up when the last answer was better than expected and turned over when it was worse. The learning rate is 0.03. Biases move by the same rule with $$h_j$$ replaced by 1.

After each update the incoming weights of every neuron are rescaled so that $\sum_j|w_{ij}|\le\rho=0.85$ . Without this a neuron could raise its goodness on the positive phase by growing its weights alone, which is the degenerate solution the loss invites. The bound also changes the dynamics: with the slope of tanh at most 1, the recurrent input to any neuron is at most 0.85 in magnitude, and a difference between two memory states shrinks by at least that factor at every settling step, so the settled activity is dominated by $$b_i+(Bx)_i$$ rather than by what the neuron's neighbours are doing.

"Without backpropagation" has a precise scope. MACE is a differentiable model trained by backpropagation, and its forces are derivatives of its energy; Quantum ESPRESSO minimises a functional by its own means. The claim is only that no global gradient is propagated through the fly-derived controller: no loss is differentiated across neurons or across settling steps, and nothing is differentiated through the sequence of past decisions.

The rule is local by construction, and that locality is exactly what leaves it unable, by itself, to assign credit across neurons or across steps; a neuron is told whether its own activity was on the right side of $\theta$ , never whether its contribution to the score was what made the chosen alloy good or bad. Chapter 10 shows what that looks like in the three runs.

Chapter 05

Refractory alloys on the bcc lattice

The refractory metals are the transition metals that melt hottest: tungsten at 3695 K, tantalum at 3290 K, molybdenum at 2896 K, niobium at 2750 K, with vanadium, hafnium, zirconium and titanium a step below. Four of them, Mo, Nb, Ta and W, are body-centred cubic at every temperature up to melting. Ti, Zr and Hf are hexagonal close-packed at room temperature and become bcc only above roughly 1150 K, 1140 K and 2000 K respectively; V is bcc throughout. A crystal structure shared by most of the family is what makes them mix.

Solid solutions of many principal elements

A conventional alloy has one base metal and small additions. Around 2004 Cantor and Yeh independently asked what happens when four or five elements are mixed in roughly equal amounts, and found that some such mixtures form a single disordered solid solution instead of the several intermetallic compounds a phase diagram would suggest. The many-component equiatomic mixture was named a high-entropy alloy, because the configurational entropy of a random solution grows with the number of components and can, at high temperature, outweigh the enthalpy that favours ordered compounds. Senkov and colleagues applied the idea to refractory metals in 2010: MoNbTaW and MoNbTaVW cast as single-phase bcc solid solutions with yield strengths above 400 MPa at 1600 °C, a temperature at which nickel superalloys have melted.

Whether a given mixture prefers one disordered phase or a set of ordered ones is decided by the balance between mixing enthalpy and mixing entropy, by the size mismatch of the atoms, and by the presence of competing compounds such as Laves or σ phases. The first of these is a quantity a total-energy method can compute.

The body-centred cubic cell: eight corner sites, each shared by eight neighbouring cells, and one site at the body centre, two atoms per cell. Each atom has eight nearest neighbours at a distance of $\sqrt{3}\,a/2$. A candidate alloy in these runs is a 2×2×2 block of such cells with an element placed on each of its sixteen sites.

The sixty-four candidates

Every candidate is a solid solution of three, four or five elements drawn from Mo, Nb, Ta, W, V, Ti, Zr and Hf, at equal or nearly equal fractions, on the bcc lattice. The supercell is 2×2×2 conventional cells, sixteen sites. For a quaternary that is four atoms of each element; for a quinary, sixteen is not divisible by five, so the composition is near-equiatomic with one element given four sites and the others three. Sites are assigned at random with a fixed seed, once per composition, so that a candidate is one definite arrangement of atoms and not an average over arrangements.

The lattice parameter follows Vegard's rule,

\[a=\sum_i x_i\,a_i ,\]

the composition-weighted mean of the pure elements' bcc lattice parameters. For MoNbTaW it gives $$a = 3.229$$ Å. The rule is an empirical linear interpolation; deviations of a per cent or two are common in real alloys and the geometry is not relaxed before it is handed to the evaluators, except in the pressure run described in the next chapter.

The mixing enthalpy

The figure of merit for two of the three runs is the mixing enthalpy per atom against the pure bcc elements,

\[\Delta H_{\mathrm{mix}}=\frac{E_{\mathrm{alloy}}}{N}-\sum_i x_i\,\frac{E_i}{N_i},\]

where $E_{\mathrm{alloy}}$ is the total energy of the sixteen-atom cell, $$E_i$$ the energy of a pure bcc cell of element $$i$$ computed by the same evaluator at the same settings, and $$N$$ , $$N_i$$ the numbers of atoms. A negative value means the elements would rather sit mixed on the bcc lattice than segregated into pure bcc metals; a positive value means the reverse. The subtraction removes the large and method-dependent absolute energies, and because the references are computed by the evaluator that scores the alloy, a MACE mixing enthalpy and a Quantum ESPRESSO mixing enthalpy are each internally consistent even though their absolute energies are on unrelated scales.

Two conventions in this definition should be kept in view. The references for Ti, Zr and Hf are bcc, which at 0 K is not their ground state; the mixing enthalpy so defined measures mixing on the shared lattice, not stability against the elements as they actually are, and it is more negative than a formation energy against hcp references would be. And a sixteen-site cell with a single random occupancy is a small and particular sample of a disordered alloy: it carries whatever short-range correlations the seed happened to produce, and it cannot represent the long-wavelength disorder of a real solid solution. A special quasirandom structure would be chosen to match the pair correlations of the ideal random alloy; these cells were not. The energies are therefore comparable among themselves, all computed the same way, and should be read as a ranking of candidates rather than as thermodynamic data for the materials.

Chapter 06

Enthalpy, entropy, temperature and pressure

A total-energy calculation is done at zero temperature and, unless the cell is relaxed under load, at zero pressure. An alloy in service is neither. The third run therefore rewarded a quantity that depends on both: the Gibbs energy of mixing at a temperature $$T$$ and a pressure $$P$$ ,

\[\Delta G_{\mathrm{mix}}(T,P)=\Delta H_{\mathrm{mix}}(P)-T\,\Delta S_{\mathrm{conf}} ,\]

with $$T = 1500$$ K and $$P = 5$$ GPa. The two terms come from different places and carry different standing, and the run kept them apart on every recorded evaluation.

Configurational entropy

Put $$N$$ atoms of $$n$$ species on $$N$$ lattice sites at random. The number of distinguishable arrangements is the multinomial coefficient $N!/\prod_i (x_iN)!$ , and Boltzmann's $S=k_B\ln\Omega$ with Stirling's approximation gives the entropy per atom of the ideal random solution,

\[\Delta S_{\mathrm{conf}}=-k_B\sum_i x_i\ln x_i ,\]

which for an equiatomic mixture is $k_B\ln n$ . The logarithm is slow: $\ln 2 = 0.69$ , $\ln 3 = 1.10$ , $\ln 4 = 1.39$ , $\ln 5 = 1.61$ . In energy units $k_B\ln 5 = 0.139$ meV K⁻¹ per atom, so at 1500 K the entropy term of an equiatomic quinary is $T\Delta S_{\mathrm{conf}}\approx 0.21$ eV/atom, and that of a ternary about 0.14 eV/atom. Mixing enthalpies of these alloys are tens of meV/atom. At 1500 K the entropy term is the larger one by an order of magnitude, and it depends only on how many elements are mixed and in what proportions.

Left, an ordered arrangement of two species: one configuration, $S_{\mathrm{conf}}=0$. Right, the same atoms placed at random: the ideal-solution entropy $-k_B\sum_i x_i\ln x_i$. Below, the equiatomic value $\ln n$ for two to five components.

The consequence for a search is mechanical. Ranked by $\Delta G_{\mathrm{mix}}$ at 1500 K, the five-component alloys come first, then the quaternaries, then the ternaries, and only within a group does $\Delta H_{\mathrm{mix}}(P)$ decide the order. That is a property of the ideal-solution model, not a finding about the alloys. The model also leaves out everything that is not the counting of arrangements: vibrational entropy, which differs between alloys by a few tenths of $$k_B$$ per atom; electronic and magnetic entropy; short-range order, which lowers the configurational entropy below the ideal value whenever like or unlike neighbours are preferred; and the competing intermetallic phases whose formation is exactly what a high mixing entropy is supposed to suppress. A negative $\Delta G_{\mathrm{mix}}$ in this model says that the random solution is preferred to the pure bcc elements at 1500 K; it does not say the random solution is the stable phase.

Pressure

Under a hydrostatic pressure the quantity that is minimised at fixed $$T$$ and $$P$$ is the enthalpy, $$H=E+PV$$ , and the mixing enthalpy at pressure is the same difference as before with $$H$$ in place of $$E$$ ,

\[\Delta H_{\mathrm{mix}}(P)=\frac{H_{\mathrm{alloy}}(P)}{N}-\sum_i x_i\,\frac{H_i(P)}{N_i},\]

each enthalpy evaluated at the volume the crystal takes under that pressure. Five gigapascals is a modest load for these metals: bulk moduli are 160–310 GPa, so the volume shrinks by two or three per cent and the $$PV$$ term per atom is of order 0.5 eV, most of which cancels in the difference. What survives is the change in mixing energy with compression, which is small but real, and which is what a search "at pressure" is meant to see.

A supercell under hydrostatic pressure. Positions and cell vectors are relaxed until the internal stress balances $P$; the reported quantity is the enthalpy $H=E+PV$ at that volume.

The pressure was computed, not modelled. For every alloy and every pure-element reference the screening potential relaxed the atomic positions and the cell shape and volume at 5 GPa, using the FIRE optimiser on a cell filter that adds the external pressure to the stress, until the largest force fell below 0.02 eV/Å or 500 steps had passed; a relaxation that did not converge produced no result and no reward. The residual pressure the potential reported at the final geometry was recorded. Quantum ESPRESSO then ran a single self-consistent calculation at exactly that geometry and reported its own pressure there, which was recorded as a diagnostic: where the two potentials disagree about the equation of state, the DFT pressure at the MACE-relaxed cell departs from 5 GPa, and the departures ranged over a few gigapascals in the run. There was no relaxation in DFT. That is a cost decision, stated plainly: a DFT cell relaxation of a sixteen-atom refractory cell costs several times a single-point calculation, and the run's budget of ten DFT evaluations was spent on ten alloys rather than on three relaxed ones.

Temperature

Temperature enters only through $T\Delta S_{\mathrm{conf}}$ . There is no thermal expansion, no phonon free energy, no electronic excitation and no melting in the model; the geometry at 1500 K is the geometry at 0 K under 5 GPa. So the objective of the third run is best read as the zero-temperature mixing enthalpy at pressure, corrected by the ideal entropy of mixing at 1500 K, which is the simplest thermodynamic model that distinguishes a ternary from a quinary at all.

Chapter 07

Machine-learned interatomic potentials: MACE

A density-functional calculation of a sixteen-atom refractory cell takes about ten minutes on six processor cores. A search that wants to look at forty candidates before spending that kind of time needs something that answers in seconds and is usually right. An interatomic potential is a function from atomic positions and species to an energy, fitted to reproduce a reference method; a machine-learned one is fitted with a flexible model to a large set of reference calculations rather than with a physically motivated formula to a few.

Locality

Every such potential rests on one assumption: the energy is a sum of atomic contributions, each of which depends only on the atom's surroundings within a cutoff radius,

\[E=\sum_i E_i\big(\{\mathbf r_j-\mathbf r_i : |\mathbf r_j-\mathbf r_i|<r_c\}\big).\]

Metals screen charge over a few ångströms, so for total energies the assumption is good with $$r_c$$ of 5–6 Å, which in a bcc refractory metal includes the first four or five neighbour shells. Long-range electrostatics, which the assumption excludes, matter little in a metal. The function $$E_i$$ must be invariant under translation, under permutation of identical atoms, and under rotation and reflection of the environment; how those invariances are built in is what distinguishes one family of potentials from another.

An atom's local environment. Only neighbours inside the cutoff $r_c$ enter its energy $E_i$. Two-body terms see distances; higher body-order terms see angles between pairs of neighbours and, at order four and beyond, dihedral relations. The total energy is the sum over atoms.

Body order and equivariance

The atomic cluster expansion writes $$E_i$$ as a sum of terms of increasing body order: a sum over single neighbours of functions of distance, a sum over pairs of neighbours of functions of two distances and one angle, and so on. Each term is built from products of one-particle basis functions, radial functions times spherical harmonics, and made invariant by contracting with the Clebsch–Gordan coefficients that couple angular momenta. Written this way the expansion is complete, and the body order controls how much of the many-body physics the model can represent.

MACE, introduced by Batatia, Kovács, Simm, Ortner and Csányi in 2022, uses that construction inside a message-passing neural network. Each atom carries a feature vector that is not a scalar but a collection of components transforming as spherical tensors of various ranks under rotation: the features are equivariant, so a rotated environment gives correspondingly rotated features rather than different ones. In each of two message-passing layers an atom's features are updated from those of its neighbours through a many-body message of body order four, built by the cluster-expansion product and contraction. Two layers with a 5–6 Å cutoff give each atom an effective receptive field of about twice that, and the high body order per layer means very few layers are needed. The energy is read out from the invariant part of the final features, and forces follow by differentiating the energy with respect to positions, which is why the model is trained on forces as well as energies and is smooth by construction.

A foundation potential

Fitting a potential used to mean choosing a chemical system and computing a training set for it. MACE-MP-0, published by Batatia and about fifty co-authors in 2023, was instead trained once on the Materials Project trajectory set: some 1.6 million configurations from the relaxations of 146 000 inorganic crystals, computed with the PBE functional, covering 89 elements. It is a single model meant to be used, without refitting, on any inorganic material, with an accuracy that is respectable across the periodic table and excellent nowhere in particular. The medium-sized model, run in double precision, evaluated a sixteen-atom refractory cell in 0–3 s in these runs, including the relaxation at pressure in the third run.

Its errors are systematic in a way that matters for a screen. The Materials Project trajectories consist largely of structures near equilibrium and of compounds rather than random solid solutions, and the model is known to soften the potential-energy surface: it tends to underestimate energies of strained or unusual environments relative to relaxed ones. In these runs the mixing enthalpies it returned were more negative than the DFT values by 60 meV/atom on average at zero pressure and by about 100 meV/atom at 5 GPa, with the ordering of candidates agreeing at a Spearman correlation of 0.77–0.78. A bias that is nearly constant across candidates does not disturb a ranking; the scatter about it does, and at 70–110 meV/atom mean absolute error the scatter is of the same size as the spread of mixing enthalpies being ranked.

The screen is therefore used for what it is good at: sorting sixty-four candidates into a shortlist of the most promising fourteen, at a cost that is negligible against one DFT calculation. It is not used as a source of numbers. Every mixing enthalpy that is reported as a result comes from Quantum ESPRESSO, and the MACE value for the same crystal is kept beside it as a record of what the screen predicted.

Chapter 08

Density functional theory in Quantum ESPRESSO

The energy of a crystal is, in principle, the lowest eigenvalue of a Schrödinger equation for all its electrons in the field of its nuclei. For sixteen atoms with several valence electrons each that wavefunction has too many coordinates to store, let alone to optimise. Density functional theory removes the difficulty by changing the variable.

Hohenberg, Kohn and Sham

Hohenberg and Kohn proved in 1964 that the ground-state energy of an electron gas in an external potential $v_{\mathrm{ext}}(\mathbf r)$ is a functional of the electron density $n(\mathbf r)$ alone, and that the true density minimises it. Kohn and Sham showed the next year how to use this: introduce a fictitious system of non-interacting electrons with the same density, described by single-particle orbitals $\psi_i$ that satisfy

\[\Big[-\frac{\hbar^2}{2m}\nabla^2+v_{\mathrm{eff}}(\mathbf r)\Big]\psi_i(\mathbf r)=\varepsilon_i\,\psi_i(\mathbf r),\qquad v_{\mathrm{eff}}=v_{\mathrm{ext}}+v_{\mathrm{H}}[n]+v_{\mathrm{xc}}[n],\]

where $v_{\mathrm{H}}$ is the classical electrostatic potential of the density and $v_{\mathrm{xc}}$ , the exchange–correlation potential, collects everything the non-interacting picture leaves out. The density is $n=\sum_i f_i|\psi_i|^2$ over occupied orbitals, and because $v_{\mathrm{eff}}$ depends on $$n$$ the equations are solved self-consistently: guess a density, build the potential, solve for the orbitals, form a new density, mix it with the old and repeat until the change is below a tolerance. The theory is exact if $v_{\mathrm{xc}}$ is; in practice it is approximated. These runs used the generalised-gradient functional of Perdew, Burke and Ernzerhof, which makes the exchange–correlation energy per electron a function of the local density and its gradient. For transition metals PBE gives lattice parameters within about a per cent and mixing energies to within a few tens of meV/atom of experiment where comparison is possible, and it is the functional the screening potential was trained against, so the two evaluators share a reference.

Plane waves and pseudopotentials

In a periodic crystal the orbitals are Bloch functions and can be expanded in plane waves, $\psi_{n\mathbf k}=\sum_{\mathbf G}c_{n\mathbf k}(\mathbf G)\,e^{i(\mathbf k+\mathbf G)\cdot\mathbf r}$ , over reciprocal-lattice vectors $\mathbf G$ . The basis is complete and unbiased and is truncated by one number, the kinetic-energy cutoff $\hbar^2|\mathbf k+\mathbf G|^2/2m\le E_{\mathrm{cut}}$ . Its weakness is the core: the tightly bound inner electrons and the nodes of valence orbitals near the nucleus would need enormous cutoffs to represent.

The pseudopotential construction. Inside the core radius $r_c$ the all-electron valence wavefunction oscillates and the Coulomb potential $-Z/r$ diverges; both are replaced by smooth functions that reproduce the all-electron wavefunction and its scattering properties exactly beyond $r_c$. The core electrons are frozen and removed from the calculation.

The pseudopotential replaces the nucleus and core electrons by an effective potential that is smooth inside a core radius and identical to the true one outside it, chosen so that the valence wavefunctions it produces agree with the all-electron ones beyond $$r_c$$ . The projector-augmented-wave method of Blöchl (1994) keeps the all-electron wavefunction formally, through a linear transformation between smooth pseudo-wavefunctions and the true ones, and so recovers all-electron accuracy at pseudopotential cost. These runs used PAW datasets from Dal Corso's PSlibrary 1.0.0 with the semicore $$s$$ and $$p$$ states of each metal treated as valence: for Mo, for instance, 4s, 4p, 4d and 5s, fourteen electrons per atom. Semicore states matter for refractory metals under pressure and in alloys, where the outer core is not inert. The wavefunction cutoff was 50 Ry and the density cutoff 400 Ry, the eightfold ratio that PAW and ultrasoft datasets need for the augmentation charges.

Sampling, smearing and convergence

Quantities like the density and the energy are integrals over the Brillouin zone, approximated by sums over a mesh of $\mathbf k$ -points. A 3×3×3 Monkhorst–Pack mesh was used for the 16-atom cell, equivalent in density to a 6×6×6 mesh for the two-atom conventional cell. In a metal the occupied and empty states meet at the Fermi surface, and a sharp step in occupation makes the sum converge slowly and erratically with mesh density. The remedy is to smear the occupations: the Marzari–Vanderbilt cold-smearing function was used with a width of 0.02 Ry (0.27 eV), a choice that damps the oscillations while contributing an error in the energy that is of second order in the width and, for this scheme, nearly free of the spurious negative occupations of a Gaussian.

Self-consistency was declared when the estimated energy error fell below $10^{-7}$ Ry, with Broyden mixing of the density at a fraction 0.3 per iteration, a conservative value for metals in which charge sloshing between iterations is common. Each calculation ran on six MPI ranks and took a median of about ten minutes. The total energy in rydberg was converted to electronvolts per atom, $1\ \mathrm{Ry}=13.6057$ eV, and the pure-element references described in chapter 5 were computed with identical settings before the first alloy.

What a converged SCF is

A self-consistent calculation that converges yields the PBE energy of one arrangement of atoms in one cell at one set of numerical parameters. It is not a statement that the cutoffs and the $\mathbf k$ -mesh are converged, that the sixteen-site cell is large enough, that the geometry is the equilibrium one, or that the alloy forms. Each of those is a separate study: cutoff and mesh by systematic tightening until differences fall below a chosen tolerance, cell size by comparing supercells, geometry by relaxation, formation by comparison against competing phases. None was done inside these runs, whose purpose was to see whether a controller could learn to pick calculations, and the mixing enthalpies reported in chapter 10 should be read with that in mind: comparable among themselves, computed to a stated recipe, and not yet materials data.

Chapter 09

The closed loop: two fidelities and an acquisition policy

The task the circuit is set is budgeted sequential candidate acquisition with immediate feedback. There is a finite pool of candidates, here 64 alloys, and a fixed number of queries. Each query names one candidate, sends it to an evaluator, and receives its label at once, before the next choice is made. The objective is to find good candidates within the budget: not to model the whole pool, not to predict the labels of the candidates never queried, but to have spent the queries on the alloys with the lowest mixing energy. The natural measure is regret, the gap between the best candidate found and the best that existed.

The funnel

One run passes through two evaluators of very different cost. Of the 64 candidates, 40 are screened with MACE-MP-0, at a few hundredths of a second to a few seconds each. The 14 best screened alloys form a shortlist, and 10 of those are computed with Quantum ESPRESSO at about ten minutes each. The first evaluation of each stage is fixed by protocol and is always the reference composition MoNbTaW, so that every run and every control begins from the same labelled point; the circuit chooses the remaining 39 screens and the remaining 9 DFT calculations.

Before the first acquisition, each evaluator computes one pure bcc crystal per element at its own settings. These eight references per evaluator define the mixing quantity, $\Delta H_{\text{mix}}=E_{\text{alloy}}/N-\sum_i x_iE_i$ in eV per atom, so that a MACE mixing enthalpy is measured against MACE elements and a DFT mixing enthalpy against DFT elements. MACE and QE energies are never subtracted from one another and never enter the same reward. Each stage has its own reference $$E_0$$ , which is 0 eV/atom for a mixing quantity, its own running baseline $$b$$ and its own replay memory of best and worst outcomes. The DFT stage starts learning from its own first labels, not from the screen's.

The funnel of one run: 64 candidate alloys, 40 screened with MACE, the 14 best screened shortlisted, 10 computed with density functional theory.

One decision, in order

Every remaining candidate at the current stage is encoded into its 20-number vector $$x$$ : the 15 composition and geometry descriptors, the fraction of budget remaining, the previous reward, the two fidelity flags and a constant. The circuit is settled for three steps from the same memory state for each candidate in turn, so that the order in which candidates are enumerated cannot change the state any of them is scored against. Each settling gives a score, the mean goodness $\overline{h^2}$ over all neurons. The scores pass through a softmax at temperature 0.06, the result is mixed with 15 % uniform exploration, and one candidate is sampled; its selection probability is recorded.

The evaluator then runs. Its answer $$E_t$$ in eV/atom becomes the reward $q_t=\tanh[(E_0-E_t)/s]$ with $$s=0.1$$ eV/atom, so that a mixing enthalpy of $$-0.1$$ eV/atom earns a quality of about 0.76 and a positive one is penalised. The baseline is subtracted, $\delta_t=q_t-b_t$ , and the baseline moves, $b_{t+1}=b_t+0.1\,\delta_t$ . The labelled outcome joins the stage's memory; the best and worst in that memory are replayed as positive and negative phases, and every synapse and bias is updated by the local rule of chapter 4, scaled by $\delta_t$ . Finally the activity the circuit settled to for the chosen candidate is committed as the new memory state, and the next decision starts from it. On the whole-brain graph a decision over the 64 candidates took about 13 s and an update about 2.6 s.

One decision of the loop: the controller scores the remaining candidates and samples one; MACE screens it; the shortlist of the best screened alloys goes to Quantum ESPRESSO; the evaluator's answer becomes a reward that returns to every synapse.

What kind of acquisition this is

Three standard methods sit next to this one. A multi-armed bandit treats each candidate as an independent arm with no features, so nothing learned about one alloy transfers to another; here the candidates share descriptors, and the whole point of the circuit is to score an unseen alloy from them. Bayesian optimisation fits a surrogate with a predictive uncertainty, usually a Gaussian process, and chooses the candidate that maximises an acquisition function such as expected improvement or an upper confidence bound; the surrogate says what it expects and how sure it is, and the acquisition function trades the two. Active learning queries where the model is least sure, in order to improve the model rather than to find the optimum.

The circuit is the acquisition function directly. It has no surrogate, no predictive distribution and no uncertainty estimate of any kind; it produces a score, and the softmax turns the score into a probability. Exploration comes only from the temperature and the 15 % uniform mixture, which are the same for a candidate the circuit has never seen anything like and for one it has. This is a real difference from Bayesian optimisation, not a rephrasing of it, and it is one of the things chapter 11 puts on the list of what to change.

The shortlist makes the loop a form of multi-fidelity acquisition, because a cheap evaluator is used to decide where the expensive one is spent. But it is a fixed form. The rule "the 14 best screened go forward, and 10 of them are computed" is set in advance; the circuit never decides whether the next query should be a screen or a DFT calculation, and it never decides to spend DFT on an alloy the screen ranked low. Where MACE and DFT disagree, and they do by 60 to 100 meV/atom in bias with a Spearman rank correlation near 0.77, the shortlist inherits the screen's errors.

Caching and accounting

An identical geometry sent to the same evaluator at the same settings is not recomputed; the stored answer is returned. This matters when several runs or several controls select the same alloy, and it means no control can obtain a label it did not choose to acquire. A cached answer still consumes a query. If it did not, a policy that happened to pick alloys some earlier run had already computed would appear more sample-efficient than it was, and comparisons between policies would be unfair. Wall-clock time is a separate quantity and is not charged again for a hit.

A calculation that fails, for instance a self-consistent-field loop that does not converge, consumes its query and returns no label. The candidate is not marked negative; it is simply not in the memory, and the budget is smaller by one. In the three runs made, no evaluation failed.

Chapter 10

What three runs showed

Three runs were made on the same design space of 64 near-equiatomic bcc alloys, each with 40 MACE evaluations and 10 Quantum ESPRESSO calculations, and no evaluation failed. The first used the 16,384-neuron reduction of the connectome and rewarded the mixing enthalpy at 0 K on unrelaxed cells. The second used the whole brain, 163,972 neurons and 6,143,838 edges, with the same objective. The third used the whole brain and rewarded the Gibbs energy of mixing at 1500 K and 5 GPa, with every cell relaxed at pressure by the screen before the DFT calculation.

The results are reported here in the order in which they deserve to be trusted: the density functional numbers first, then the accuracy of the screening model measured against them, then what the controller did.

The alloys

At 0 K, ten alloys reached Quantum ESPRESSO. Four have a negative same-structure mixing enthalpy: MoNbTaW and MoTaVW at −58 meV/atom, MoTaV at −44 meV/atom and MoNbTaVW at −36 meV/atom. Every alloy containing Zr or Hf is positive, between +62 and +91 meV/atom, and TiVW sits at +19 meV/atom. This is the expected chemistry. The group 5 and 6 metals have similar atomic radii and mix nearly ideally on the bcc lattice; Zr and Hf are larger and are not bcc at 0 K, and a rigid bcc cell with a compromise lattice parameter makes them pay for both facts.

Quantum ESPRESSO results for the two whole-brain runs. Left, the mixing enthalpy $\Delta H_{mix}$ at 0 K on unrelaxed cells. Right, the Gibbs mixing energy $\Delta G_{mix}$ at 1500 K and 5 GPa, split into the enthalpy at pressure and the entropy term $-T\Delta S_{conf}$.

At 1500 K and 5 GPa the quinary MoNbTiVW leads at $\Delta G_{mix}=-247$ meV/atom. Its enthalpy at pressure is only −40 meV/atom; the other −207 meV/atom is $T\Delta S_{conf}$ for five equiatomic components. MoNbTaW and MoNbTaVW follow at −233 meV/atom, MoTaVW at −223 meV/atom. The entropy term, 142 to 207 meV/atom depending on the number of components, is three to five times any enthalpy in the table, so the ranking is by number of components first and by enthalpy second. Chapter 6 says why that is a property of the ideal-solution model rather than a finding about these alloys.

Relaxation at pressure changed the two families of alloys differently. The size-mismatched alloys gained a great deal: MoNbWZr went from +91 meV/atom unrelaxed to +51 meV/atom relaxed, HfNbTiW from +63 to +14 meV/atom. The well-matched alloys hardly moved: MoNbTaW from −58 to −54 meV/atom. Both are physically sensible; local relaxation is what a misfit atom needs and a well-fitted atom does not. Every relaxed cell remained bcc, with a median cell strain of 0.8 % and no atom displaced by more than 0.31 Å. The DFT pressure at the screen-relaxed cells ranged from 1.5 to 5.0 GPa against the 5 GPa target, which says that MACE's equation of state for the Mo- and W-rich alloys is softer than PBE's.

None of these numbers is a property of a material. Each is a converged single-point energy of one 16-atom occupancy of a bcc supercell, without vibrations, without competing phases and without the many configurations a real solid solution samples.

The screening model

On the twenty alloys that both evaluators saw at identical geometries, MACE-MP-0 reproduces the DFT ordering reasonably and the DFT values poorly. The Spearman rank correlation is 0.77 on the unrelaxed cells and 0.78 on the cells relaxed at 5 GPa. The mean absolute error is 70 and 114 meV/atom. The bias, MACE minus DFT, is −60 and −97 meV/atom: the model is systematically too negative. It agrees with DFT on the sign of $\Delta H_{mix}$ for 8 of 10 unrelaxed alloys and for 5 of 10 relaxed ones.

MACE-MP-0 against Quantum ESPRESSO on the same crystals, in eV/atom. The rank order is fair; the values are systematically too negative and the slope is wrong. A least-squares line through the relaxed points gives DFT ≈ 0.28 × MACE + 28 meV/atom.

A straight line through the relaxed points has slope 0.28, so the model's dynamic range on these cells is three to four times too large. The likely cause is distributional. MACE-MP-0 was trained on relaxation trajectories of mostly ordered compounds near their equilibria; a random 16-atom bcc cell with four or five species at a compromise lattice parameter is far from that distribution, and the model over-rewards the local relaxations it finds there. The pure-element references relax to geometries the model knows well, so the error lands on the alloy side of $\Delta H_{mix}$ .

For this project the consequence is exact. A shortlist of the best ten by MACE will contain most of DFT's best ten, so MACE is fit to order candidates. Rewarding a controller with MACE's $\Delta H_{mix}$ , as the screening stage does, teaches it a landscape that is three times too steep and wrong in sign for half of the borderline alloys.

The controller

The three runs used graphs ten times apart in size and received quite different reward streams: the 1500 K run's reward-minus-baseline signal was almost entirely positive, the two 0 K runs' mostly negative. The probability with which the controller picked each alloy is nevertheless the same curve in all three, to plotting precision.

Top, the probability with which the circuit picked each chosen alloy, for the three runs; the curves coincide and equal $1/n_{remaining}$, the probability of uniform sampling over the alloys not yet evaluated. Bottom, the reward-minus-baseline signal $r-b$ each run received.

That curve is $1/n_{remaining}$ , where $n_{remaining}$ is the number of alloys not yet evaluated. It is the shape of uniform sampling with a fixed random seed. The forty alloys screened were the same set in all three runs, and the first twelve were the same alloys in the same order.

Re-scoring all 64 candidates with the saved initial and final weights makes the same point without reference to the seed. A candidate's initial score, the mean goodness $\overline{h^2}$ after settling, correlates with the norm of its 20-number description, $\lVert x\rVert$ , at $$r=1.000$$ . Whatever the 6.1 million recurrent synapses added to the score, it was a function of the input's length and nothing else. Learning then moved the weights a long way, $\lVert\Delta w\rVert_2=111$ on the whole-brain 0 K run and 66 on the 1500 K run, and moved every selection probability by at most 0.0012 and 0.0015 respectively. The whole policy lives between 0.013 and 0.019 per alloy, where uniform over 64 is 0.0156.

Selection probability over the 64 candidates with the initial weights and with the final weights of the whole-brain runs. Uniform sampling would be 0.0156 for each. The largest change in any probability is 0.0015.

Said plainly: in these three runs the fly brain did not choose the alloys. The seed did, with a slight lean towards Hf- and Zr-rich compositions because they have the largest feature norms. MACE and Quantum ESPRESSO did all of the physics, and the DFT rankings above are exactly as valid as they would be for forty alloys drawn at random, which is what they are.

Why

The cause is mechanical, not statistical, and chapter 3 anticipates it. The mean goodness over 164,000 neurons is a law-of-large-numbers quantity. Each neuron receives $$(Bx)_i$$ , a fixed random combination of the 20 inputs, and the mean of $\tanh^2$ over so many independent random projections depends on the magnitude of the drive, $\lVert Bx\rVert\propto\lVert x\rVert$ , and on almost nothing else. The recurrent term does not change this: it adds to each neuron another sum of many random contributions whose statistics are again set by the overall level of activity.

A bounded local update on synapses whose incoming weights are held to a fixed absolute sum can raise every neuron's goodness together, and it did; the 0 K run's scores all rose by about 0.125. It cannot open a gap between two candidates that the input did not already separate, because no term in the update knows which candidate it is looking at except through the same drive. Chapter 4 says the rule is local by construction and does not, by itself, assign credit across neurons. This is what that looks like when it is run at scale.

The features were not normalised, so their norm carried information about composition, mainly through the mean atomic number and the lattice parameter. With the description scaled to unit length the initial scores would have lost their main source of variation, and whatever remained would have had to come from the direction of $$x$$ rather than its length. Either way the circuit was not in the causal path between the reward and the next choice. Chapter 11 sets out what has to be controlled, and in what order, before any statement about the wiring can be tested.

Chapter 11

How to tell whether the wiring matters

The question the project asks is whether a network shaped like a fly's brain, trained by a local rule, chooses calculations better than an arbitrary network of the same size would. Three runs with one topology and one learning rule cannot answer it, however they turn out, because they vary nothing. The answer has the form of a comparison, and the comparison has to be designed so that only the thing in question differs.

Matched controls

Two properties are entangled in the controller: the topology, which comes from the connectome, and the learning rule, which comes from Forward-Forward. Each can be replaced while the other is held fixed, which gives a small matrix of conditions. The topology axis has three entries. The measured graph. A degree-preserving rewiring of it, in which directed edges are swapped in pairs so that every neuron keeps its exact in-degree and out-degree and the multiset of edge weights is unchanged, but who connects to whom is scrambled; if the measured wiring carries useful structure, this control loses it while keeping everything a degree sequence determines. And a generic sparse random graph with the same numbers of nodes and edges, which keeps only the size and density. The learning axis has two entries: the local Forward-Forward rule of chapter 4, and a control that differentiates the same per-neuron objective through the three settling steps, so that derivatives do flow between neurons within one decision. Everything else, the encoder $$B$$ , the initial weights, the softmax temperature, the exploration rate, the candidates, the budgets and the evaluators, is identical across the six cells.

Four further controls bracket the matrix. Random selection, which chooses candidates uniformly, measures what luck alone achieves under the same budget. A frozen circuit, with the same topology and initial weights but no learning, separates what the architecture does at initialisation from what learning adds. An encoder-only scorer, a ridge regression from the twenty features to the observed rewards with no recurrent network at all, tests whether the features suffice by themselves. And shuffled rewards, in which the reward is assigned to the wrong candidate, tests whether learning depends on the association between an action and its outcome, as it must if it is learning anything.

What chapter 10 changes

The runs described in the previous chapter found that the initial score was a monotone function of the feature norm, with a correlation of 1.000, and that learning moved the weights without moving the selection probabilities. This reorders the controls. Before the topology can be tested, the encoder has to be fixed: the features must be standardised so that their length carries no information about their content, and the encoder-only control must be run first, because if a ridge regression from standardised features already picks well, the recurrent network has to beat that and not uniform sampling. The frozen circuit and the shuffled-reward control then say whether what remains is learning. Only after those does the rewiring comparison mean anything.

An effect, if there is one, has to survive all of these. A controller that beats random selection but not the frozen circuit has a useful initialisation, not a useful rule; one that beats the frozen circuit but not the rewired graph has a useful rule on an arbitrary graph; only one that beats the rewired graph has learned something the wiring provided.

Measuring a search

The natural outcome of a budgeted search is the best value found, but with ten DFT evaluations out of sixty-four candidates the best value found is dominated by which fourteen reached the shortlist, and a single number per run is a poor statistic. Regret is better: at each step the difference between the best value found so far and the best value in the pool, summed over the budget, which rewards finding good candidates early and penalises finding them late. On a pool this small the true best can be computed once and the regret is exact. For the comparison to hold up, the measure must be chosen before the runs, seeds must be many, the rewiring must be repeated with independent swaps, and the uncertainty must be reported over independent seeds and tasks rather than over the fifty correlated decisions within one run. Three runs, one seed each, are a pilot; a claim would need tens.

What has to change

The list is short and concrete. Standardise the twenty features. Use the synapse counts, which are anatomical evidence, to set the magnitudes of the initial weights, so that the measured graph enters as more than a mask. Lengthen the budgets, or make the candidates cheaper, so that a rule with a learning rate of 0.03 has enough decisions to change a policy; fifty decisions were not enough. Give the controller an estimate of what it does not know, since an acquisition function without uncertainty can only exploit, and a search that only exploits is a ranking of the initial scores. And extend credit beyond one step: a search is a sequence, and the value of a screen lies in what it does to the shortlist, which the immediate reward does not see. Each of these is an experiment with a control, and the machinery to run them is what the three runs actually tested.

Chapter 12

Notes and sources

The chapters lean on the following. Journal references are given by volume and page; the preprints by arXiv number.

The fly

MaleCNS v1.0, the connectome of an adult male Drosophila melanogaster central nervous system, FlyEM project at HHMI Janelia Research Campus with the University of Cambridge, the MRC Laboratory of Molecular Biology and Google Research. Data and renderings at male-cns.janelia.org, CC BY 4.0.
Y. Aso et al., "The neuronal architecture of the mushroom body provides a logic for associative learning", eLife 3, e04577 (2014); and Y. Aso et al., "Mushroom body output neurons encode valence and guide memory-based action selection in Drosophila", eLife 3, e04580 (2014).
N. Frémaux and W. Gerstner, "Neuromodulated spike-timing-dependent plasticity, and theory of three-factor learning rules", Frontiers in Neural Circuits 9, 85 (2016).

Learning

G. Hinton, "The Forward-Forward Algorithm: Some Preliminary Investigations", arXiv:2212.13345 (2022).
R. J. Williams, "Simple statistical gradient-following algorithms for connectionist reinforcement learning", Machine Learning 8, 229 (1992), for the reward-minus-baseline form of the learning signal.

Alloys

B. Cantor, I. T. H. Chang, P. Knight and A. J. B. Vincent, "Microstructural development in equiatomic multicomponent alloys", Materials Science and Engineering A 375–377, 213 (2004); J.-W. Yeh et al., "Nanostructured high-entropy alloys with multiple principal elements: novel alloy design concepts and outcomes", Advanced Engineering Materials 6, 299 (2004).
O. N. Senkov, G. B. Wilks, D. B. Miracle, C. P. Chuang and P. K. Liaw, "Refractory high-entropy alloys", Intermetallics 18, 1758 (2010); O. N. Senkov, G. B. Wilks, J. M. Scott and D. B. Miracle, "Mechanical properties of Nb₂₅Mo₂₅Ta₂₅W₂₅ and V₂₀Nb₂₀Mo₂₀Ta₂₀W₂₀ refractory high entropy alloys", Intermetallics 19, 698 (2011).
L. Vegard, "Die Konstitution der Mischkristalle und die Raumfüllung der Atome", Zeitschrift für Physik 5, 17 (1921).
A. Zunger, S.-H. Wei, L. G. Ferreira and J. E. Bernard, "Special quasirandom structures", Physical Review Letters 65, 353 (1990), for what the sixteen-site cells are not.

Interatomic potentials

R. Drautz, "Atomic cluster expansion for accurate and transferable interatomic potentials", Physical Review B 99, 014104 (2019).
I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner and G. Csányi, "MACE: Higher order equivariant message passing neural networks for fast and accurate force fields", Advances in Neural Information Processing Systems 35 (2022).
I. Batatia et al., "A foundation model for atomistic materials chemistry", The Journal of Chemical Physics 163, 184110 (2025); preprint arXiv:2401.00096 (2024). Describes MACE-MP-0 and the Materials Project trajectory training set.
E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, "Structural relaxation made simple", Physical Review Letters 97, 170201 (2006), the FIRE optimiser; A. H. Larsen et al., "The atomic simulation environment—a Python library for working with atoms", Journal of Physics: Condensed Matter 29, 273002 (2017).

Density functional theory

P. Hohenberg and W. Kohn, Physical Review 136, B864 (1964); W. Kohn and L. J. Sham, Physical Review 140, A1133 (1965).
J. P. Perdew, K. Burke and M. Ernzerhof, "Generalized gradient approximation made simple", Physical Review Letters 77, 3865 (1996).
P. E. Blöchl, "Projector augmented-wave method", Physical Review B 50, 17953 (1994).
A. Dal Corso, "Pseudopotentials periodic table: From H to Pu", Computational Materials Science 95, 337 (2014), the PSlibrary.
H. J. Monkhorst and J. D. Pack, "Special points for Brillouin-zone integrations", Physical Review B 13, 5188 (1976).
N. Marzari, D. Vanderbilt, A. De Vita and M. C. Payne, "Thermal contraction and disordering of the Al(110) surface", Physical Review Letters 82, 3296 (1999), the cold-smearing scheme.
P. Giannozzi et al., "QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials", Journal of Physics: Condensed Matter 21, 395502 (2009); P. Giannozzi et al., "Advanced capabilities for materials modelling with Quantum ESPRESSO", Journal of Physics: Condensed Matter 29, 465901 (2017).

Drawings

The engraved diagrams were drawn in MetaPost with fiziko, a library by Sergey Slyusarev (GPL-3.0-or-later; its manual is CC BY-SA 4.0); the drawings themselves are the author's.
The crystal renders in the run pages were made with OVITO: A. Stukowski, "Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool", Modelling and Simulation in Materials Science and Engineering 18, 015012 (2010).
The brain rendering in chapter 1 is the MaleCNS authors' own, reproduced under CC BY 4.0.