[Studio picture · a steel ball resting in a watch glass]

[04 / MODEL 3]

The computed proton wells

Each hydrogen sits in a pocket between metal atoms. The shape of that pocket decides how easily the hydrogens could move together, so it was computed from the electrons up for nickel hydride and palladium deuteride.

[IN SHORT]

Each hydrogen sits in a pocket between metal atoms. The ordering model depends on the shape of that pocket more than on any other input: how stiff it is, how fast its walls steepen, and how hard one hydrogen pulls on its neighbours.

All three were computed with density functional theory for nickel hydride (NiH) and palladium deuteride (PdD) with one hydrogen per metal atom.

Nickel holds its hydrogen in a stiff pocket. Palladium holds its deuterium in a soft one, and its neighbours pull on it strongly.

Stretch the metal and the pocket softens fast. In PdD one percent of uniform stretch takes the stiffness below zero.

g2_input g3_skew/partA g3_loading

1.05 eV/Å² stiffness of the NiH pocket 0.181 eV/Å² stiffness of the much softer PdD pocket +1% uniform stretch that turns the PdD stiffness negative 35.0 meV bare harmonic quantum of one D in its well

[04.1 · THE POCKET]

A ball in a bowl

Picture each hydrogen as a ball resting in a bowl. The metal atoms around it and the electrons between them make the bowl. Three properties of the bowl decide how the hydrogens behave as a crowd.

[01 · STIFFNESS]

How steep the bowl is

Stiffness measures how fast the energy rises as the hydrogen moves off centre. A stiff pocket is a narrow, steep bowl; a soft one is wide and shallow. The unit is eV/Å², energy per distance squared.
[02 · STEEPENING]

How fast the walls rise

A perfect bowl rises with the square of the distance. These pockets rise faster further out, and fastest when the hydrogen heads straight for a metal atom. Physicists call this the anharmonic part of the well.
[03 · COUPLING]

How hard it pulls its neighbours

When one hydrogen moves, the electrons around it shift and the neighbouring hydrogens feel a tug. That tug is the coupling. It is what could make many hydrogens move in step.

These three numbers are the inputs the ordering model is most sensitive to, so they were computed instead of guessed. The tool is density functional theory (DFT), a standard quantum-mechanical method that solves for the electrons in a crystal. Here it runs in the GPAW code with the PBE functional. The metal atoms are held in place, a hydrogen is moved by a small step, and the energy and forces are read off.

[04.2 · TWO CALCULATIONS]

Moving together and moving alone

Two kinds of cell answer two questions. In a two-atom cell every hydrogen in the crystal moves at once. In a 64-atom cell one hydrogen moves and the other 31 hold still, so the forces on them can be read.

[FIG 1]Two ways to move the hydrogens

Two-atom cell · every hydrogen moves the same way

This gives the collective stiffness: how hard it is to move every hydrogen at once.

64-atom cell · one moves, the others are held

This gives the local stiffness of one hydrogen and the pull it exerts on each neighbour.

Schematic. The bowls are the pockets and the springs stand for the coupling through the electrons. Moving everything together leaves the springs unchanged, so the collective stiffness is lower than the local one by the sum of the pulls.

[THE BOOKKEEPING]

collective = local − sum of the pulls

For PdD in the 64-atom cell with every deuterium site filled, the numbers in eV/Å² are:

One D moves alone 0.591

Pull of all neighbours −0.389

All D together, 64-atom cell 0.202

All D together, two-atom cell 0.181

The neighbours' pull cancels most of the local stiffness, so moving together is far easier than moving alone. The 64-atom cell was built and sampled independently and lands close to the two-atom cell.

[04.3 · THE SHAPE]

Stiff nickel and soft palladium

When every hydrogen moves together, the NiH pocket has a stiffness of 1.050 eV/Å² and the PdD pocket has 0.181. Both steepen sharply further out. Move the slider to stretch the metal and watch the PdD bowl go flat.

[FIG 2]Energy of the pocket against displacement

Direction

View

Uniform stretch 0.00%

NiH stiffness 1.050 eV/Å²

PdD stiffness 0.181 eV/Å²

full wellharmonic part onlycomputed point

What to notice: the gap between each solid curve and its dashed harmonic part is the steepening of the walls, which is largest along [100]. Stretch by +1% and the PdD curve loses its curvature at the centre. Curves are V(u) = ½K₂u² + ¼K₄u⁴ fitted to the computed points (two-atom cell, all hydrogens moving together). Points are computed at 0 and +1% (PdD only); between them the stiffness follows a curve through the computed values at −1, 0 and +1%, and the PdD quartic term is interpolated linearly. Source: g2_input

[ALONG 100]

The hydrogen heads straight for a metal atom. The wall rises steeply: the quartic coefficient K₄ is 16.2 eV/Å⁴ in NiH and 11.2 in PdD.

[ALONG 111]

The hydrogen heads for the gap between three metal atoms. The wall rises gently: K₄ is 1.35 eV/Å⁴ in NiH and 3.29 in PdD.

[04.4 · VIBRATION]

Calculated against measured

A hydrogen in its pocket rattles. The stiffness sets the energy of one quantum of that rattle, which experiments measure as a peak in the vibration spectrum.

Counting only the bottom of the well and moving every hydrogen together, the calculation gives 66.0 meV for NiH and 19.4 meV for PdD. That is 26 and 48 percent below the measured 89 and 37 meV. The wells are strongly anharmonic: their walls steepen so fast that the bottom alone does not set the vibration energy.

One deuterium moving alone in the 64-atom cell while its neighbours are held still has a bare harmonic quantum of 35.0 meV. That value counts only the curvature at the bottom of its own well, with no anharmonic walls and no neighbour effects. It characterises the well itself. The measured peak also carries the effects left out here.

[FIG 3]Vibration energy, meV

NiH

All H together, bottom of the well 66.0

Measured peak 89

PdD

All D together, bottom of the well 19.4

One D alone, bare harmonic, 64-atom cell 35.0

Measured peak 37

What to notice: the bottom of the collective well alone falls well short of the measured peaks. The one-D bar is the bare harmonic value of a single D in its well; the measured peak also includes the anharmonic walls and the neighbours. Calculated values are harmonic energies ħω = ħ√(K₂/m) with the proton mass for NiH and the deuteron mass for PdD. Sources: g2_input g3_loading

[04.5 · STRETCH]

Stretch the metal and the well flattens

Pull the metal atoms apart and the pocket softens. NiH loses about 22 percent of its stiffness per percent of stretch. PdD would lose 113 percent, so at +1% its stiffness is already −0.013 eV/Å².

[FIG 4]Collective stiffness against strain

Strain

Scale

NiHPdDcomputedPdD at the functional's own spacing

What to notice: in absolute terms the two metals lose stiffness at a similar rate; on the relative scale PdD falls far faster because it starts from so little. The PdD line crosses zero just below +1% uniform stretch. Strain is the fractional change in length. Curves pass through the computed points at −1, 0 and +1%. Source: g2_input g3_skew/partA

[FRACTIONAL CHANGE IN STIFFNESS PER UNIT STRAIN]

Material	Uniform	Along the motion	Across the motion
NiH	−22.0	−23.9	+1.2
PdD	−113	−80.9	−15.9

Direction matters. Stretching along the line the hydrogen moves on softens the pocket most. Stretching across it leaves NiH almost unchanged and still softens PdD.

Negative does not mean free. A negative stiffness says the bottom of the bowl is flat or has a small hump in the middle. The steep walls and the hydrogen's own zero-point motion still hold it in place.

[THE MAIN UNCERTAINTY]

Which lattice spacing to use

The PdD numbers use a lattice spacing of 4.084 Å. That is the published extrapolation to full loading of measurements made at 77 K over loadings 0.80 to 0.98 (Schirber and Morosin, Phys. Rev. B 12, 117, 1975). There is no direct measurement at x = 1 behind it. When the spacing is left free, the same calculation settles at 4.134 Å, which is 1.2 percent larger. There the collective stiffness is −0.055 eV/Å² instead of +0.181. Because one percent of stretch decides the sign, the choice of spacing matters directly. Strain thresholds counted from 4.134 Å would be about 1.2 percentage points smaller. The fitted strain law, which never saw this point, predicts −0.058. A smaller uncertainty sits in the 64-atom cell. The corrected single-D well there gives a collective stiffness of 0.202 eV/Å², 11 percent above the converged two-atom value of 0.181. The cell's 3×3×3 k-point mesh was not converged further.

4.084 Å · published extrapolation to full loading: +0.181
4.134 Å · the functional's own: −0.055

[04.6 · SHEAR AND FILMS]

Shear it or grow it as a film

There are two more ways to deform PdD. Both were computed in the two-atom cell with every deuterium moving together. Shear softens the pocket steadily. In-plane strain on a thin film trades stiffness between the two directions.

[FIG 5]PdD under shear

What to notice: the lowest stiffness falls from 0.181 to 0.035 eV/Å² at 8% shear along [100] (0.030 along [110]) and would reach zero near 9 to 10%. Shear is the engineering shear strain of the cell. Source: g3_skew/partA

[FIG 6]PdD film with its height free to relax

motion in the planemotion across the layer

What to notice: squeeze the film and motion across the layer goes soft; stretch it and motion in the plane goes soft. At every strain in this range one direction is below zero. With no imposed strain the film still grows 3.2% taller because the functional prefers a larger spacing than 4.084 Å. Source: g3_skew/partA

[04.7 · LOADING]

Fewer than one hydrogen per metal atom

Real samples are seldom full. The 64-atom PdD cell was run again with some deuterium sites empty. With 4 of the 32 empty (loading x = 0.875) the gaps were either spread out or bunched together. With 8 of the 32 empty (x = 0.75) they were bunched. An empty site changes the well of the deuterium next to it.

[FIG 7]On-site stiffness of PdD against empty neighbours

x = 1x = 0.875x = 0.75full or spreadbunchedfitted law

What to notice: each empty neighbour lowers the stiffness, and the smaller lattice of a lower loading raises it. At x = 0.875 with two empty neighbours there are two points, one spread and one bunched. The law gives them the same value and the calculation does not. Each point is one D moved while every other atom is held in place. Source: g3_loading

[01 · LESS HYDROGEN MEANS A STIFFER SITE]

Lower loading shrinks the lattice: 4.084 Å at x = 1, 4.062 at 0.875 and 4.039 at 0.75. A D that keeps all its nearest neighbours gets stiffer as it does: 0.591, 0.703 and 0.732 eV/Å².

[02 · AN EMPTY NEIGHBOUR SOFTENS AND PUSHES]

A D beside empty sites is softer (0.529 eV/Å² with two, 0.477 with four) and feels a steady force toward the vacancy (0.070 and 0.136 eV/Å), so it does not sit at the centre of its site. The atoms were held on their ideal positions, so the size of that shift is only a linear estimate: about 0.13 Å in the spread cell.

[03 · ARRANGEMENT MATTERS]

Both cells have the same loading, the same spacing and two empty neighbours. In the bunched cell the D is softer (0.431 against 0.529 eV/Å²), pulls far more weakly on its neighbours (0.128 against 0.399) and is pushed harder (0.109 against 0.070 eV/Å) than in the spread cell.

[04 · A TWO-TERM LAW]

(6)K_ii = K_ii(x = 1) + c_ε (a / 4.084 Å − 1) + c_n n_vac

K_ii is the on-site stiffness, a the lattice constant and n_vac the number of empty nearest-neighbour sites. The fit gives c_ε = −14.4 eV/Å² and c_n = −0.077 eV/Å² per empty neighbour. The root-mean-square residual is 0.045 eV/Å², and the law is at most 0.108 off when a whole cell is left out and predicted. The law carries the trend with loading. It does not distinguish bunched from spread; the computed cells are the reference for that. The couplings are fitted shell by shell in the same form.

Local well of one D in a 64-atom PdD cell by loading and vacancy arrangement
Loading x	Arrangement	Lattice constant a Å	Empty nearest neighbours	On-site stiffness eV/Å²	Sum of neighbour pulls eV/Å²	Force at its site toward the vacancy eV/Å
1	every site filled	4.084	0	0.591	0.389	0
0.875	spread	4.062	0	0.703	0.485	0.002
0.875	spread	4.062	2	0.529	0.399	0.070
0.75	bunched	4.039	0	0.732	0.508	0.002
0.75	bunched	4.039	4	0.477	0.246	0.136
0.875	bunched	4.062	2	0.431	0.128	0.109

One D is moved by ±0.03 Å. The on-site stiffness comes from forces by eq. (2), and the sum of pulls from eq. (4). The force is along the direction probed, which faces a vacancy when the D has empty neighbours. 4.084 Å is the published extrapolation to full loading of measurements at 77 K over loadings 0.80 to 0.98; the reduced-loading lattice constants follow the same line. The x = 0.75 loading lies below 0.80, so its value is extrapolated along that line.

[04.8 · THE DETAIL]

Equations, method and numbers

The well

(1)V(u) − V(0) = ½ K₂ u² + ¼ K₄ u⁴

u is the displacement of the hydrogen from the centre of its site (Å), along [100] or [111]. V is the total energy (eV). K₂ is the stiffness, the curvature at the centre (eV/Å²). K₄ is the quartic coefficient (eV/Å⁴), which says how fast the walls steepen. K₂ is fixed by eq. (2); K₄ is fitted to energies out to 0.30 Å in the two-atom cell and 0.24 Å in the 64-atom cell, where it is fitted to energies and forces together. The site is symmetric, so V(−u) = V(u).

The stiffness from forces

(2)K₂ = −[F(+h) − F(−h)] / 2h

F is the force on the displaced hydrogen along the displacement (eV/Å) and h the step: 0.01 Å in the two-atom cell, 0.03 Å in the 64-atom cell.

The vibration energy

(3)ħω = ħ √(K₂ / m)

m is the proton mass (NiH) or the deuteron mass (PdD) and ħ the reduced Planck constant. This is the harmonic energy of the bottom of the well, not an anharmonic excitation energy.

Collective, local and coupling

(4)K_coll = K_ii − Σ_j≠i J_ij, J_ij = −Φ_ij

K_ii is the on-site stiffness of hydrogen i with every other nucleus held fixed. Φ_ij = −[F_j(+h) − F_j(−h)] / 2h is the force constant between i and hydrogen j, both along x, and J_ij is the pull: positive when j is drawn along with i. The sum runs over the other 31 hydrogen sites of the periodic 64-atom cell, each counted once; the metal atoms stay fixed. K_coll is then the stiffness for all hydrogens moving together.

The strain response

(5)d ln K / dε ≈ [K(+1%) − K(−1%)] / (0.02 K(0))

ε is the strain. For uniform strain the cell and all positions scale by 1 + ε. For strain along or across the motion only the x axis scales, and the hydrogen moves along x or along y. The PdD strain law fits all computed strain cases (uniform, uniaxial, shear and film) with a quadratic in the strain tensor of 12 coefficients; its root-mean-square error is 0.008 eV/Å² and its largest error 0.044 eV/Å².

Method

Code: GPAW 26.7, plane waves, PAW setups with 10 valence electrons for Ni and Pd
Functional: PBE, non-magnetic, Fermi smearing 0.1 eV
Cutoff: 500 eV. At 600 eV the stiffness changes by −5.55% (NiH) and +1.78% (PdD)
Lattice: NiH a = 3.732 Å, PdD a = 4.084 Å (published extrapolation to full loading), rock-salt, hydrogen on the octahedral site
Metal atoms: held fixed; the electrons relax self-consistently around every displaced hydrogen
Two-atom cell: k-mesh 8×8×8 (NiH) and 12×12×12 (PdD), each the first mesh within 5% of the next denser one
64-atom cell: 32 metal atoms, 32 hydrogen sites. PdD and NiH: 3×3×3 k-points, 500 eV, tight force convergence. The earlier runs (one k-point, 400 eV) are kept in Tables C and D for comparison
Compute: 164 core-minutes for the first set of wells (g2_input), 200 for the PdD shear and film cases, about 2,510 for the PdD 64-atom wells at all loadings (1,091 of it for x = 1) and about 2,870 with the NiH check. Failed and stopped attempts are included

g2_input g3_skew/partA g3_loading

A · Collective stiffness against k-mesh two-atom cell, 500 eV

Material	k-mesh	K₂ eV/Å²	ħω meV
NiH	6	1.173	69.8
NiH	8 accepted	1.050	66.0
NiH	12	1.028	65.3
NiH	16	1.017	65.0
PdD	6	0.421	29.5
PdD	8	0.221	21.4
PdD	12 accepted	0.181	19.4
PdD	16	0.175	19.1

PdD is very sensitive to k-point sampling: at 6×6×6 its stiffness is more than twice the converged value.

B · Collective stiffness under strain eV/Å², two-atom cell

Case	−1%	0	+1%	dK/dε	d ln K/dε
NiH, uniform	1.295	1.050	0.834	−23.1	−22.0
NiH, along the motion	1.314	1.050	0.812	−25.1	−23.9
NiH, across the motion	1.033	1.050	1.057	+1.2	+1.2
PdD, uniform	0.397	0.181	−0.013	−20.5	−113
PdD, along the motion	0.339	0.181	0.046	−14.6	−80.9
PdD, across the motion	0.209	0.181	0.152	−2.9	−15.9
PdD, uniform, at 4.134 Å			−0.055 at +1.2%

C · Quartic coefficients K₄ eV/Å⁴

Case	[100]	[111]
NiH, all H together	16.2	1.35
PdD, all D together	11.2	3.29
PdD, all D together, +1% uniform stretch	9.88	3.26
PdD, one D alone, 64-atom cell, 3×3×3 k	11.6	2.23
NiH, one H alone, 64-atom cell	16.2 recomputed	1.13 earlier run

NiH [100] was recomputed at 500 eV and 3×3×3 k-points from one displaced point at 0.24 Å, with K₂ held at its recomputed value (the earlier run gave 18.3). NiH [111] was not recomputed: it is still the value from the earlier run at 400 eV and one k-point.

D · The 64-atom cell eV/Å², eq. (4)

Run	K_ii	ΣJ	K_coll	ħω meV
PdD, 3×3×3 k, 500 eV	0.591 force 0.599 energy	0.389	0.202	35.0
NiH, 3×3×3 k, 500 eV	1.712 force 1.714 energy	0.676	1.036
NiH, earlier run at coarser settings (one k-point, 400 eV) replaced	1.47	0.23	1.23	78.1
PdD, earlier run at coarser settings (one k-point, 400 eV) replaced	1.09	0.064	1.03	47.6

The earlier runs used coarser numerical settings together: a 400 eV cutoff, a single k-point and looser force convergence. For PdD they gave a much smaller pull (0.064 against 0.389) and a stiffer site. At the tighter settings the NiH pull is about three times larger (0.68 against 0.23). The NiH collective value of 1.036 now agrees with the two-atom cell's 1.050 to within 1.3 percent, and the earlier 1.234 did not.

E · Where the pull comes from PdD, 64-atom cell, 3×3×3 k

Neighbour distance Å	Sites in the cell	Sum of pulls eV/Å²
2.888	12	+0.226
4.084	3	+0.269
5.002	12	−0.109
5.776	3	+0.007
7.074	1	−0.005
Total	31	+0.389

Sites are the distinct hydrogen sites of the periodic cell at each distance, so the far shells hold fewer sites than in an infinite crystal. Negative entries push rather than pull.

F · Shear PdD, lowest collective stiffness, eV/Å²

Shear	0%	2%	4%	6%	8%
along [100]	0.181	0.155	0.119	0.078	0.035
along [110]					0.030

G · Thin film with its height relaxed PdD, eV/Å²

In-plane strain	Height / reference	In the plane	Across the layer
−4%	1.103	0.868	−0.341
−2%	1.066	0.395	−0.290
0%	1.032	0.084	−0.181
+2%	1.000	−0.127	0.057
+4%	0.971	−0.281	0.421

These are small-vibration values at the site centre. A negative entry marks a direction where the well is flat or double-bottomed; zero-point motion and the anharmonic walls add stiffness back. The +2% row kept its starting height of 1.000 because that already met the stress tolerance (residual stress 1.6 × 10⁻⁴ eV/Å³ against a tolerance of 5 × 10⁻⁴; the other rows reached about 10⁻⁶). It lies on the trend of the other four rows, which interpolate to about 1.001.

[WHAT THIS MEANS FOR THE NEXT MODEL]

One hydrogen in its real well

The next model puts one hydrogen into exactly this well, with the stiffness, the steep walls and the neighbour pull computed here. It then asks when the motion locks into step. Strain is its main control: in PdD about one percent of stretch takes the well from stiff to flat.

MODEL Model 4: the anharmonic oscillator Next: a static route under stretch and a pumped route MODEL Model 2: what could make them order Previous: five candidate ingredients and the two routes that survive SOURCE The computed wells in the repository Two-atom cells, strain response, the first 64-atom pass SOURCE Strain, shear and film laws The PdD strain law and the spacing check SOURCE Loading The x = 1 local well, the reduced loadings, the fitted law and the NiH check

In short