Molecular dynamics of carbon and classical and quantum hydrogen in iron

In the HEmS project, we got interested in how hydrogen behaves in the crystal lattice of α-iron. Part of our approach was to make molecular dynamics (MD) simulations, and on this page I want to show off some movies that we made.

In order to do MD you need a way to calculate interatomic forces. In all our work we use the quantum mechanical magnetic tight binding approximation [1]. Details of the models can be found in Ref [2].

First we made some classical MD simulations to look at the dynamics of hydrogen, H, and carbon, C, when they are trapped at a vacancy (that is a lattice site with an atom missing). These two movies show simulations at 400K.

The left hand, or upper if your display is small, movie (simulation time, 2.9 picoseconds) has two C atoms and one H atom trapped in the vacancy. Note these points:

The right hand, or lower, movie (simulation time, 3.6 picoseconds) shows two hydrogen atoms trapped at the vacancy. We know that the vacancy will trap up to five H atoms and that they like to occupy positions close to the centres of the faces of the cube enclosing the vacancy [3].

Hydrogen is the lightest element and generally speaking it is not sufficient to treat is as a classical object—you would never treat electrons as classical, which is why we use quantum mechanics to describe interatomic forces. In the first two videos it's almost legitimate to use classical MD since the temperature is 400K; but even at room temperature hydrogen behaves non classically—its diffusivity as a function of temperature is non-Arrhenius as it can tunnel between tetrahedral sites of the perfect iron crystal lattice.

If I showed a video of hydrogen in perfect iron at room temperature then it would not last long as the H atom would quickly wander away out of shot (it's not trapped as in the vacancy movies above). So the next movie (simulation time 5 picoseconds) shows hydrogen in perfect iron at 100K. It does jump between tetrahedral sites but no too frequently in the total simulation time of 5 picoseconds.

We do not treat the quantum mechanics of the electrons and the H nucleus—the proton—on the same footing. The electrons are represented as Bloch states of the tight binding hamiltonian; the proton dynamics are simulated using path integral molecular dynamics, based on Richard Feynman's path integral representation of the quantum partition function. In this theory, a quantum particle is simulated as a "necklace" of beads, each bead is an image of the proton. The dynamics of this object is represented in the movie exactly as such a necklace. Even the iron atoms look a bit fuzzy as their nuclei are treated in the same way, and indeed in that picture their position is indeed uncertain as expressed in Heisenberg's uncertainty principle.

To show that the proton is quite noticeably delocalised even at room temperature, the final movie shows the proton, again trapped at a vacancy so that it won't wander away. Note the perspective is the same as before and now the cube in front has its body centred atom missing.

Acknowledgments. The MD simulations were made using the TBE branch of the questaal code (see, which has been interfaced to the i-ipi suite by Michele Ceriotti. The movies were made using C++/vtk software built by Dimitar L. Pashov, using files generated for Michael Methfessel's xbs program.

  1. Anthony T. Paxton and Michael W. Finnis, Magnetic tight binding and the iron-chromium enthalpy anomaly, Physical Review B, 77, art. 024428 (2008)
  2. Anthony T. Paxton and Christian Elsässer, Analysis of a carbon dimer bound to a vacancy in iron using density functional theory and a tight binding model, Physical Review B, 87, art. 224110 (2013)
  3. Anthony T. Paxton From quantum mechanics to the physical metallurgy of steels, Materials Science and Technology, 30, 1063-1070 (2013)
  4. Anthony T. Paxton and Ivaylo H. Katzarov, Quantum and isotope effects on hydrogen diffusion, trapping and escape in iron, Acta Materialia, 103, 71-76 (2016)
  5. Ivaylo H. Katzarov, Dimitar L. Pashov and Anthony T. Paxton, Fully quantum mechanical calculation of the diffusivity of hydrogen in iron using the tight-binding approximation and path integral theory, Physical Review B, 88, art. 054107 (2013)
  6. Bingqing Cheng, Anthony T. Paxton and Michele Ceriotti, Hydrogen diffusion and trapping in α-Fe: the role of quantum and anharmonic fluctuations, Physical Review Letters, 120, art. 225901 (2018)

Last modified: 27 Aug 2020 (Tony Paxton)