Keith Refson

First principles methods using CASTEP

Abstract The CASTEP code for first principles electronic structure calculations will be described. A brief, non-technical overview will be given and some of the features and capabilities highlighted. Some features which are unique to CASTEP will be described and near-future development plans outlined.

Moldy: a portable molecular dynamics simulation program for serial and parallel computers

Abstract Moldy is a highly portable C program for performing molecular-dynamics simulations of solids and liquids using periodic boundary conditions. It runs in serial mode on a conventional workstation or on a parallel system using an interface to a parallel communications library such as MPI or BSP. The “replicated data” parallelization strategy is used to achieve reasonable performance with a minimal difference between serial and parallel code. The code has been optimized for high performance in both serial and parallel cases. The model system is completely specified in a run-time input file and may contain atoms, molecules or ions in any mixture. Molecules or molecular ions are treated in the rigid-molecule approximation and their rotational motion is modeled using quaternion methods. The equations of motion are integrated using a modified form of the Beeman algorithm. Simulations may be performed in the usual NVE ensemble or in isobaric and/or isothermal ensembles. Potential functions of the Lennard–Jones, 6-exp and MCY forms are supported and the code is structured to give an straightforward interface to add a new functional form. The Ewald method is used to calculate long-ranged electrostatic forces.

Reproducibility in density functional theory calculations of solids

A comparison of DFT methods Density functional theory (DFT) is now routinely used for simulating material properties. Many software packages are available, which makes it challenging to know which are the best to use for a specific calculation. Lejaeghere et al. compared the calculated values for the equation of states for 71 elemental crystals from 15 different widely used DFT codes employing 40 different potentials (see the Perspective by Skylaris). Although there were variations in the calculated values, most recent codes and methods converged toward a single value, with errors comparable to those of experiment. Science, this issue p. 10.1126/science.aad3000; see also p. 1394 A survey of recent density functional theory methods shows a convergence to more accurate property calculations. [Also see Perspective by Skylaris] INTRODUCTION The reproducibility of results is one of the underlying principles of science. An observation can only be accepted by the scientific community when it can be confirmed by independent studies. However, reproducibility does not come easily. Recent works have painfully exposed cases where previous conclusions were not upheld. The scrutiny of the scientific community has also turned to research involving computer programs, finding that reproducibility depends more strongly on implementation than commonly thought. These problems are especially relevant for property predictions of crystals and molecules, which hinge on precise computer implementations of the governing equation of quantum physics. RATIONALE This work focuses on density functional theory (DFT), a particularly popular quantum method for both academic and industrial applications. More than 15,000 DFT papers are published each year, and DFT is now increasingly used in an automated fashion to build large databases or apply multiscale techniques with limited human supervision. Therefore, the reproducibility of DFT results underlies the scientific credibility of a substantial fraction of current work in the natural and engineering sciences. A plethora of DFT computer codes are available, many of them differing considerably in their details of implementation, and each yielding a certain “precision” relative to other codes. How is one to decide for more than a few simple cases which code predicts the correct result, and which does not? We devised a procedure to assess the precision of DFT methods and used this to demonstrate reproducibility among many of the most widely used DFT codes. The essential part of this assessment is a pairwise comparison of a wide range of methods with respect to their predictions of the equations of state of the elemental crystals. This effort required the combined expertise of a large group of code developers and expert users. RESULTS We calculated equation-of-state data for four classes of DFT implementations, totaling 40 methods. Most codes agree very well, with pairwise differences that are comparable to those between different high-precision experiments. Even in the case of pseudization approaches, which largely depend on the atomic potentials used, a similar precision can be obtained as when using the full potential. The remaining deviations are due to subtle effects, such as specific numerical implementations or the treatment of relativistic terms. CONCLUSION Our work demonstrates that the precision of DFT implementations can be determined, even in the absence of one absolute reference code. Although this was not the case 5 to 10 years ago, most of the commonly used codes and methods are now found to predict essentially identical results. The established precision of DFT codes not only ensures the reproducibility of DFT predictions but also puts several past and future developments on a firmer footing. Any newly developed methodology can now be tested against the benchmark to verify whether it reaches the same level of precision. New DFT applications can be shown to have used a sufficiently precise method. Moreover, high-precision DFT calculations are essential for developing improvements to DFT methodology, such as new density functionals, which may further increase the predictive power of the simulations. Recent DFT methods yield reproducible results. Whereas older DFT implementations predict different values (red darts), codes have now evolved to mutual agreement (green darts). The scoreboard illustrates the good pairwise agreement of four classes of DFT implementations (horizontal direction) with all-electron results (vertical direction). Each number reflects the average difference between the equations of state for a given pair of methods, with the green-to-red color scheme showing the range from the best to the poorest agreement. The widespread popularity of density functional theory has given rise to an extensive range of dedicated codes for predicting molecular and crystalline properties. However, each code implements the formalism in a different way, raising questions about the reproducibility of such predictions. We report the results of a community-wide effort that compared 15 solid-state codes, using 40 different potentials or basis set types, to assess the quality of the Perdew-Burke-Ernzerhof equations of state for 71 elemental crystals. We conclude that predictions from recent codes and pseudopotentials agree very well, with pairwise differences that are comparable to those between different high-precision experiments. Older methods, however, have less precise agreement. Our benchmark provides a framework for users and developers to document the precision of new applications and methodological improvements.

Computer simulation of interlayer water in 2:1 clays

Monte Carlo computer simulation has been used to study water confined between the layers of 2:1 clay minerals. The model systems are based on natural Mg and Na smectites. The simulation cells contain one clay layer, 64 water molecules and four magnesium or eight sodium interlayer cations. These atoms and molecules interact with each other through a new set of effective pair potentials, which we discuss. The calculations are conducted in constant (N,p,T) ensembles, at T=300 K and with a uniaxial pressure, p, of 1 M Pa applied normal to the clay sheets. All the molecules, including the clay sheets, are therefore allowed to move during the simulations. The calculated equilibrium layer spacing is 14.7±0.1 A with interlayer Mg2+ and 14.2±0.1 A with interlayer Na+. These spacings compare with experimental values of 15.1 A and 14.5 A, measured for Mg and Na saturated Chambers montmorillonite, at 79% relative humidity. The corresponding densities and average potential energies of the interlayer water molecules ar...

WATER CHEMISORPTION AND RECONSTRUCTION OF THE MGO SURFACE

Cavendish Laboratory, University of Cambridge, Madingley Road, Cambridge CB3 0HE(December 6, 1994)The observed reactivity of MgO with water is in apparent conflict with theoretical calculationswhich show that molecular dissociation does not occur on a perfect (001) surface. We have per-formed ab-initio total energy calculations which show that a chemisorption reaction involving areconstruction to form a (111) hydroxyl surface is strongly preferred with ∆E = −90.2 kJ mol

An ab initio study of hydrogen in forsterite and a possible mechanism for hydrolytic weakening

Even small amounts of water can profoundly effect the physical properties of minerals. In olivine 2 eV more favourable to put the first proton into the Si vacancy than the magnesium site, the presence of water will certainly act to increase the population of silicon vacancies. In fact, in the presence of water the energy required to form a Si vacancy is perhaps less than that to form an Mg vacancy. This is in stark contrast to dry olivine where Si vacancies are many eV less favourable. If creep is rate limited by the diffusion of the slowest species, silicon in olivine, then increasing the Si vacancy concentration could provide a mechanism for hydrolytic weakening.

The kinetics and mechanism of MgO dissolution

Abstract Reactions and atomic rearrangements at fluid–crystal interfaces play an important role in catalysis and in controlling the kinetics and mechanisms of dissolution. We have studied the attachment and reactions of water molecules at the MgO–water interface by combining measurements of 1 H and 2 D surface penetration and etch pit morphology with ab initio calculations. These studies show that the most common MgO cleavage surface, (001), is thermodynamically unstable when hydrated. Proton rearrangement on such surfaces precedes proton–cation exchange and provides a general mechanism for the detachment of ions during dissolution. The kinetics of dissolution are strongly influenced by the concentration of surface defects and a simple model based on the ab initio results predicts a dissolution rate of 10−10 mol cm−2 s−1 for a typical surface defect concentration of 0.1.

Density functional theory in the solid state.

Density functional theory (DFT) has been used in many fields of the physical sciences, but none so successfully as in the solid state. From its origins in condensed matter physics, it has expanded into materials science, high-pressure physics and mineralogy, solid-state chemistry and more, powering entire computational subdisciplines. Modern DFT simulation codes can calculate a vast range of structural, chemical, optical, spectroscopic, elastic, vibrational and thermodynamic phenomena. The ability to predict structure–property relationships has revolutionized experimental fields, such as vibrational and solid-state NMR spectroscopy, where it is the primary method to analyse and interpret experimental spectra. In semiconductor physics, great progress has been made in the electronic structure of bulk and defect states despite the severe challenges presented by the description of excited states. Studies are no longer restricted to known crystallographic structures. DFT is increasingly used as an exploratory tool for materials discovery and computational experiments, culminating in ex nihilo crystal structure prediction, which addresses the long-standing difficult problem of how to predict crystal structure polymorphs from nothing but a specified chemical composition. We present an overview of the capabilities of solid-state DFT simulations in all of these topics, illustrated with recent examples using the CASTEP computer program.

Spin singlet formation in MgTi2O4: Evidence of a helical dimerization pattern

The transition-metal spinel MgTi2O4 undergoes a metal-insulator (M-I) transition on cooling below T(M-I)=260 K. A sharp reduction of the magnetic susceptibility below T(M-I) suggests the onset of a magnetic singlet state. Using high-resolution synchrotron and neutron powder diffraction, we have solved the low-temperature crystal structure of MgTi2O4, which is found to contain dimers with short Ti-Ti distances (the locations of the spin singlets) alternating with long bonds to form helices. Band structure calculations based on hybrid exchange density functional theory show that, at low temperatures, MgTi2O4 is an orbitally ordered band insulator.

Periclase surface hydroxylation during dissolution

Abstract Periclase (001) surfaces were etched in dilute acid at pH 2 and 4. X-ray reflectivity measurements on a reference crystal constrained the initial roughness of these surfaces to be approximately 30 A. The reference crystal and the crystal reacted at pH 2 were analyzed by Elastic Recoil Detection Analysis (ERDA) for proton penetration. After reaction the etched sample showed proton penetration to a depth of at least 5000 A while the reference crystal showed no significant proton inventory. Down to 900 A, the H Mg ratio in the etched sample was approximately 2, consistent with near-surface protonation of the MgO to form hydroxylated brucite-like layers. Protonation is a far more likely mechanism to explain the proton profile than precipitation because this reaction was completed over 16 orders of magnitude below saturation with brucite. Formation of a hydroxylated near-surface layer on periclase during dissolution explains why the dissolution rates of periclase and brucite are identical in the pH range 2–5; the detachment rates are the same because the surface structures are the same. This suggests that even for this ionic solid in acid, the dissolution reaction involves a two-step mechanism with a rapid single protonation step of near-surface oxygen atoms and a slower, rate determining second protonation step. In general, product phases such as brucite are likely to be better developed under natural weathering conditions of near-neutral pH because the second step of protonation (and thus full hydration of the detaching cation, e.g., Mg+2) is much slower than in acid. Our proposed protonation mechanism relates field observations of the periclase weathering reaction to laboratory dissolution, hydration, and dehydration experiments.