Classical Interatomic Potentials Research Articles

A fully automated approach to calculate the melting temperature of elemental crystals

The interface method is a well established approach for predicting melting points of materials using interatomic potentials. However, applying the interface method is tedious and involves significant human intervention. The whole procedure involves several successive tasks: estimate a rough melting point, set up the interface structure, run molecular dynamic calculations and analyze the data. Loop calculations are necessary if the predicted melting point is different from the estimated one by more than a certain convergence criterion, or if full melting/solidification occurs. In this case monitoring the solid–liquid phase transition in the interface structure becomes critical. As different initial random seeds for the molecular dynamic simulations within the interface method induce slightly different melting points, a few ten or hundred interface method calculations with different random seeds are necessary for performing a statistical analysis on these melting points. Considering all these technical details, the work load for manually executing and combining the various involved scripts and programs quickly becomes prohibitive. To simplify and automatize the whole procedure, we have implemented the interface method into pyiron (http://pyiron.org). Our fully automatized procedure allows to efficiently and precisely predict melting points of stable unaries represented by arbitrary potentials with only two user-specified parameters (interatomic potential file and element). For metastable or dynamically unstable unary phases, the crystal structure needs to be provided as an additional parameter. We have applied our automatized approach on fcc Al, Ni, dynamically unstable bcc Ti and hcp Mg and employed a large set of available interatomic potentials. Melting points for classical interatomic potentials of these metals have been obtained with a numerical precision well below 1 K.

Efficient determination of solid-state phase equilibrium with the multicell Monte Carlo method.

Building on our previously introduced multicell Monte Carlo (MC)^{2} method for modeling phase coexistence, this paper provides important improvements for efficient determination of phase equilibria in solids. The (MC)^{2} method uses multiple cells, representing possible phases. Mass transfer between cells is modeled virtually by solving the mass balance equation after the composition of each cell is changed arbitrarily. However, searching for the minimum free energy during this process poses a practical problem. The solution to the mass balance equation is not unique away from equilibrium, and consequently the algorithm is in risk of getting trapped in nonequilibrium solutions. Therefore, a proper stopping condition for (MC)^{2} is currently lacking. In this work, we introduce a consistency check via a predictor-corrector algorithm to penalize solutions that do not satisfy a necessary condition for equivalence of chemical potentials and steer the system toward finding equilibrium. The most general acceptance criteria for (MC)^{2} is derived starting from the isothermal-isobaric Gibbs ensemble for mixtures. Using this ensemble, translational MC moves are added to include vibrational excitations as well as volume MC moves to ensure the condition of constant pressure and temperature entirely with a MC approach, without relying on any other method for relaxation of these degrees of freedom. As a proof of concept the method is applied to two binary alloys with miscibility gaps and a model quaternary alloy, using classical interatomic potentials.

Prediction of the lattice thermal conductivity of zircon and the cubic and monoclinic phases of zirconia by molecular dynamics simulation

Zirconia (ZrO2) and zircon (ZrSiO4) have high strengths and stabilities at high temperatures and also have remarkable thermal properties, such as rather low thermal conductivities. In the present research, the phonon thermal conductivity of zircon and two polymorphs (cubic and monoclinic phases) of zirconia are investigated. For this purpose, equilibrium molecular dynamics simulations using classical interatomic potentials and the Green-Kubo formalism are presented. The results are given in detail over a wide temperature range, with a 100 K temperature step, from 200 K to 2400 K and 200 K to 1400 K for the cubic and monoclinic phases of zirconia, respectively, and 200 K to 1400 K for zircon. The temperature dependence of the lattice parameters and the equilibrium atomic volumes show good agreement with available experimental data. Next, the phonon thermal conductivity is calculated by analysing the raw data of the heat current autocorrelation function. The results illustrate that the above-mentioned materials have quite low thermal conductivities that are dependent on temperature. It is also shown that the lattice thermal conductivity of the different phases of zirconia and zircon can be decomposed into three contributions due to the acoustic short-range, long-range phonon and optical phonons modes. Finally, the thermal conductivity results from this study are compared with previous experimental studies and very good agreement is found.

On the elaboration of the next generation of thermodynamic models of solid solutions.

Thermodynamic models of solid solutions used in computational thermochemistry have not been modernized in recent years. With the advent of fast and cheap computers, it is nowadays possible to add, at a minimal computational cost, physical ingredients such as coordination numbers, inter-atomic distances and classical interatomic potentials to the function describing the energetics of ordered and disordered solid solutions. As we show here, the integration of these elements into a robust statistical thermodynamic model of solution establishes natural connections with other deterministic and stochastic atomistic methods such as Monte Carlo and molecular dynamics simulations. Ultimately, all these numerical approaches need to be self-consistent and generate complementary sets of numerical thermo-physical properties. The present work proposes a new formalism to define the Gibbs free energy of ordered and disordered solid solutions. It allows for a complete prediction of the thermal, volumetric and compositional dependence of the Gibbs free energy by solving a constrained minimization problem. As a proof of concept, we explore the energetic behavior of pure face-centered cubic gold as well as the AuCu L10 ordered solution as a function of both temperature and pressure. We finally compare these results with the average properties obtained from classical molecular dynamics simulations and explain the origin of the existing differences between the two approaches based on how the temperature is accounted for in each method.

Continuous and optimally complete description of chemical environments using Spherical Bessel descriptors

Recently, machine learning potentials have been advanced as candidates to combine the high-accuracy of electronic structure methods with the speed of classical interatomic potentials. A crucial component of a machine learning potential is the description of local atomic environments by some set of descriptors. These should ideally be invariant to the symmetries of the physical system, twice-differentiable with respect to atomic positions (including when an atom leaves the environment), and complete to allow the atomic environment to be reconstructed up to symmetry. The stronger condition of optimal completeness requires that the condition for completeness be satisfied with the minimum possible number of descriptors. Evidence is provided that an updated version of the recently proposed Spherical Bessel (SB) descriptors satisfies the first two properties and a necessary condition for optimal completeness. The Smooth Overlap of Atomic Position (SOAP) descriptors and the Zernike descriptors are natural counterparts of the SB descriptors and are included for comparison. The standard construction of the SOAP descriptors is shown to not satisfy the condition for optimal completeness and, moreover, is found to be an order of magnitude slower to compute than that of the SB descriptors.

Evaluation of classical interatomic potentials for molecular dynamics simulations of borosilicate glasses

In this study, we compare two sets of classical interatomic potentials for borosilicate glasses by investigating two major commercial borosilicate glasses, Borofloat®33 (Boro33) and N-BK7®, via molecular dynamics simulations. To evaluate the performance of the two potentials, structural and elastic properties are analyzed for the two glass compositions and are compared with available experimental values. We find that the potential by Wang et al. [M. Wang, N.M. Anoop Krishnan, B. Wang, M.M. Smedskjaer, J.C. Mauro, and M. Bauchy, J. Non. Cryst. Solids, 498 294–304 (2018)] provides a closer N4 value for Boro33 but underpredicts the N4 value for N-BK7. In contrast, the N4 value of N-BK7 using the potential of Deng and Du [L. Deng and J. Du, J. Am. Ceram. Soc., 102 [5] 2482–2505 (2019)] agrees well with the experimental data, but that of Boro33 sample is overpredicted. Our result also indicates that Wang's potential gives a better prediction in the short-range structure, while Du's potential provides a closer medium-range structure compared with the experimental data. Neither set of potentials is able to provide accurate predictions of elastic moduli. Wang's potential predicts lower elastic modulus due to the under predicted N4 value, whereas Du's potential yields higher elastic modulus compared with the experimental values, resulting from its overpredicted N4 value.

Machine Learning Classical Interatomic Potentials for Molecular Dynamics from First-Principles Training Data

The ever-increasing power of modern supercomputers, along with the availability of highly scalable atomistic simulation codes, has begun to revolutionize predictive modeling of materials. In partic...

A charge optimized many-body potential for iron/iron-fluoride systems.

A classical interatomic potential for iron/iron-fluoride systems is developed in the framework of the charge optimized many-body (COMB) potential. This interatomic potential takes into consideration the effects of charge transfer and many-body interactions depending on the chemical environment. The potential is fitted to a training set composed of both experimental and ab initio results of the cohesive energies of several Fe and FeF2 crystal phases, the two fluorine molecules F2 and the F2-1 dissociation energy curve, the Fe and FeF2 lattice parameters of the ground state crystalline phase, and the elastic constants of the body centered cubic Fe structure. The potential is tested in an NVT ensemble for different initial structural configurations as the crystal ground state phases, F2 molecules, iron clusters, and iron nanospheres. In particular, we model the FeF2/Fe bilayer and multilayer interfaces, as well as a system of square FeF2 nanowires immersed in an iron solid. It has been shown that there exists a reordering of the atomic positions for F and Fe atoms at the interface zone; this rearrangement leads to an increase in the charge transfer among the atoms that make the interface and put forward a possible mechanism of the exchange bias origin based on asymmetric electric charge transfer in the different spin channels.

Development of interatomic potentials for Fe-Cr-Al alloy with the particle swarm optimization method

Due to its potential application in nuclear cladding system, Fe-Cr-Al alloys have attracted significant interest recently. Thus, it is of importance to acquire the classical interatomic interaction potential in order to predict the mechanical behavior of FeCrAl structures at the microscopic scale. In this work, the particle swarm optimization method has been applied and improved to develop the interatomic potentials for Fe, Cr, Al unitary systems and Fe-Cr, Fe-Al binary systems. The potential parameters are fitted to experimental elastic constants, first principle calculated lattice information and defect formation energies with consideration of force matching. The developed Fe-Cr-Al ternary potentials are used to predict the phase stability region for FeCrAl ferrite alloys. The bulk modulus for the single crystal and the plastic behavior under grain boundary tensile stress are also modeled to find the influence of Cr and Al contents on the mechanical properties of FeCrAl alloys. The solutes Cr would improve the mechanical properties while high Al content would lead to poor elastic and plastic properties. The predictions from the newly developed potential are in good agreement with experimental results.

Geometries of edge and mixed dislocations in bcc Fe from first-principles calculations

We use DFT to compute core structures of $a_0[100](010)$ edge, $a_0[100](011)$ edge, $a_0/2[\bar{1}\bar{1}1](1\bar{1}0)$ edge, and $a_0/2[111](1\bar{1}0)$ $71^{\circ}$ mixed dislocations in bcc Fe. The calculations use flexible boundary conditions (FBC), which allow dislocations to relax as isolated defects by coupling the core to an infinite harmonic lattice through the lattice Green function (LGF). We use LGFs of dislocated geometries in contrast to previous FBC-based dislocation calculations that use the bulk crystal LGF. Dislocation LGFs account for changes in topology in the core as well as strain throughout the lattice. A bulk-like approximation for the force constants in a dislocated geometry leads to LGFs that optimize the cores of the $a_0[100](010)$ edge, $a_0[100](011)$ edge, and $a_0/2[111](1\bar{1}0)$ $71^{\circ}$ mixed dislocations. This approximation fails for the $a_0/2[\bar{1}\bar{1}1](1\bar{1}0)$ dislocation, so here we derive the LGF using accurate force constants from a Gaussian approximation potential. The standard deviations of dislocation Nye tensor distributions quantify the widths of the cores. The relaxed cores are compact, and the magnetic moments on the Fe atoms closely follow the volumetric strain distributions in the cores. We also compute the core structures of these dislocations using eight different classical interatomic potentials, and quantify symmetry differences between the cores using the Fourier coefficients of their Nye tensor distributions. Most of the core structures computed using the classical potentials agree well with DFT results. The DFT geometries provide benchmarking for classical potential studies of work-hardening, as well as substitutional and interstitial sites for computing solute-dislocation interactions that serve as inputs for mesoscale models of solute strengthening and solute diffusion near dislocations.

Elastic constants of graphene: Comparison of empirical potentials and DFT calculations

The capacity of popular classical interatomic potentials to describe elastic properties of graphene is tested. The Tersoff potential, Brenner reactive bond-order potentials REBO-1990, REBO-2000, REBO-2002 and AIREBO as well as LCBOP, PPBE-G, ReaxFF-CHO and ReaxFF-C2013 are considered. Linear and non-linear elastic response of graphene under uniaxial stretching is investigated by static energy calculations. The Young's modulus, Poisson's ratio and high-order elastic moduli are verified against the reference data available from experimental measurements and ab initio studies. The density functional theory calculations are performed to complement the reference data on the effective Young's modulus and Poisson's ratio at small but finite elongations. It is observed that for all the potentials considered, the elastic energy deviates remarkably from the simple quadratic dependence already at elongations of several percent. Nevertheless, LCBOP provides the results consistent with the reference data and thus realistically describes in-plane deformations of graphene. Reasonable agreement is also observed for the computationally cheap PPBE-G potential. REBO-2000, AIREBO and REBO-2002 give a strongly non-linear elastic response with a wrong sign of the third-order elastic modulus and the corresponding results are very far from the reference data. The ReaxFF potentials drastically overestimate the Poisson's ratio. Furthermore, ReaxFF-C2013 shows a number of numerical artefacts at finite elongations. The bending rigidity of graphene is also obtained by static energy calculations for large-diameter carbon nanotubes. The best agreement with the experimental and ab initio data in this case is achieved using the REBO-2000, REBO-2002 and ReaxFF potentials. Therefore, none of the considered potentials adequately describes both in-plane and out-of-plane deformations of graphene.

Defect Creation in Crystals: A Portal to Directional Dark Matter Searches

A large body of astrophysical observations indicate that around 85% of the matter in the universe is not made of recognized standard model particles. Understanding the nature of this so-called dark matter is of fundamental importance to cosmology, astrophysics, and high-energy particle physics. We examine the response of commonly used semiconductor materials to low-mass WIMP interactions using numerical simulations based on classical interatomic potentials in these materials. These simulations, backed up by more precise density functional theory simulations and experiments, predict a nonlinear energy loss that never produces phonons due to the nonzero energy required to form crystallographic defects. We argue that such nonlinear effects related to defect formation in electron-volt-scale resolution semiconductor detectors allows for very effective directional sensitivity and possible statistical nuclear recoil discrimination to dark matter signals for masses below 1 $$\hbox {GeV}/c^2$$ .

Developing an interatomic potential for martensitic phase transformations in zirconium by machine learning

Atomic simulations provide an effective means to understand the underlying physics of structural phase transformations. However, this remains a challenge for certain allotropic metals due to the failure of classical interatomic potentials to represent the multitude of bonding. Based on machine-learning (ML) techniques, we develop a hybrid method in which interatomic potentials describing martensitic transformations can be learned with a high degree of fidelity from ab initio molecular dynamics simulations (AIMD). Using zirconium as a model system, for which an adequate semiempirical potential describing the phase transformation process is lacking, we demonstrate the feasibility and effectiveness of our approach. Specifically, the ML-AIMD interatomic potential correctly captures the energetics and structural transformation properties of zirconium as compared to experimental and density-functional data for phonons, elastic constants, as well as stacking fault energies. Molecular dynamics simulations successfully reproduce the transformation mechanisms and reasonably map out the pressure–temperature phase diagram of zirconium.

Classical interatomic potential for quaternary Ni–Fe–Cr–Pd solid solution alloys

In this paper, we present a new quaternary interatomic potential for the NiFeCrPd system, which is an extension on the previous NiFeCr potential. Density functional theory is used to calculate the quantities to be fitted, with particular focus on the energetics of point defects with solutes, for the potential to be used towards understanding radiation damage properties. The potential thus will enable the modeling of multi-elemental solid solution alloys consisting of up to four elements. To test the potential, we have performed atomistic kinetic Monte Carlo simulations to investigate the effect of configurational entropy on the self-diffusion coefficients. The self-diffusion coefficients are found to increase with chemical complexity, contrary to the common postulation of sluggish diffusion in high entropy alloys (HEAs). In addition, we have performed molecular dynamics simulations to elucidate the effect of Pd on vacancy diffusion and clustering in pure Ni and binary alloys. In agreement with recent irradiation experiments, our simulations show that while large vacancy clusters, such as stacking fault tetrahedra, are formed in pure Ni, Ni–Fe and Ni–Cr systems, negligible vacancy clustering is observed in Ni–Pd systems, indicating a possible effect of Pd in reducing cluster sizes. We suggest that this potential will be useful for studying the defect evolution in multi-component HEAs.

Molecular Dynamics Simulations of the Elastic Anisotropy of Pd at Extreme Conditions**Supported by the National Natural Science Foundation of China (41574076), the Basic Research of Technology Program of China under Grant No. JSHS2014404B002, the Young Core Teacher Scheme of Henan Province under Grant No. 2014GGJS-108, and key project of science and technology research of Henan Provincial Education Department under Grant

It is very interesting to discover the elastic properties of engineering material palladium, especially its elastic anisotropy along Hugoniot states. We here investigate the evolution of its high pressure and temperature (PT) elastic ansotropy along Hugoniot using molecular dynamics simulations based on accurate classical interatomic potential. In order to testify the validity of the interatomic potential of Pd in describing the high PT elastic properties, we calculate its isothermal and adiabatic elastic moduli using molecular dynamics method. The obtained data are in good agreement with experimental data. From the isothermal elastic constants, we deduce the Hugoniot acoustic velocities and find that the resulting data are in good agreement with experimental acoustic velocity data. Based on the reliable elastic constants, we further investigate the spacial elastic ansotropy along Hugoniot PT states. It is found that the spacial elastic anisotropy of Pd increases along Hugoniot states.

A tungsten-rhenium interatomic potential for point defect studies

A tungsten-rhenium (W-Re) classical interatomic potential is developed within the embedded atom method interaction framework. A force-matching method is employed to fit the potential to ab initio forces, energies, and stresses. Simulated annealing is combined with the conjugate gradient technique to search for an optimum potential from over 1000 initial trial sets. The potential is designed for studying point defects in W-Re systems. It gives good predictions of the formation energies of Re defects in W and the binding energies of W self-interstitial clusters with Re. The potential is further evaluated for describing the formation energy of structures in the σ and χ intermetallic phases. The predicted convex-hulls of formation energy are in excellent agreement with ab initio data. In pure Re, the potential can reproduce the formation energies of vacancies and self-interstitial defects sufficiently accurately and gives the correct ground state self-interstitial configuration. Furthermore, by including liquid structures in the fit, the potential yields a Re melting temperature (3130 K) that is close to the experimental value (3459 K).

Evaluating variability with atomistic simulations: the effect of potential and calculation methodology on the modeling of lattice and elastic constants

Atomistic simulations using classical interatomic potentials are powerful investigative tools linking atomic structures to dynamic properties and behaviors. It is well known that different interatomic potentials produce different results, thus making it necessary to characterize potentials based on how they predict basic properties. Doing so makes it possible to compare existing interatomic models in order to select those best suited for specific use cases, and to identify any limitations of the models that may lead to unrealistic responses. While the methods for obtaining many of these properties are often thought of as simple calculations, there are many underlying aspects that can lead to variability in the reported property values. For instance, multiple methods may exist for computing the same property and values may be sensitive to certain simulation parameters. Here, we introduce a new high-throughput computational framework that encodes various simulation methodologies as Python calculation scripts. Three distinct methods for evaluating the lattice and elastic constants of bulk crystal structures are implemented and used to evaluate the properties across 120 interatomic potentials, 18 crystal prototypes, and all possible combinations of unique lattice site and elemental model pairings. Analysis of the results reveals which potentials and crystal prototypes are sensitive to the calculation methods and parameters, and it assists with the verification of potentials, methods, and molecular dynamics software. The results, calculation scripts, and computational infrastructure are self-contained and openly available to support researchers in performing meaningful simulations.

Amorphous tantala and its relationship with the molten state

The structure factors of molten $\mathrm{T}{\mathrm{a}}_{2}{\mathrm{O}}_{5}$ and $\mathrm{N}{\mathrm{b}}_{2}{\mathrm{O}}_{5}$ have been measured by high-energy x-ray and pulsed neutron diffraction. These are compared to transmission-mode x-ray diffraction through a self-supported 15-\ensuremath{\mu}m ion-beam sputtered amorphous tantala film. Atomistic models derived from the diffraction data by means of empirical potential structure refinement reveal that tantala and niobia liquids are very close to isomorphous, as confirmed by measurement of a molten mixture, $\mathrm{T}{\mathrm{a}}_{0.8}\mathrm{N}{\mathrm{b}}_{1.2}{\mathrm{O}}_{5}$. Nonetheless, peak Nb-O bond lengths are about $1%$ shorter than those for Ta-O, at temperatures, ${T}^{*}=T/{T}_{\mathrm{melt}}$, scaled to the melting points. Mean coordination numbers are ${n}_{M\mathrm{O}}\ensuremath{\simeq}5.6(1),{n}_{\mathrm{O}M}\ensuremath{\simeq}2.23(4)$ in the liquid state, and ${n}_{\mathrm{TaO}}\ensuremath{\simeq}6.6(2),{n}_{\mathrm{OTa}}\ensuremath{\simeq}2.63(8)$ in the solid. The liquids are built from five- and six-fold $M$-O polyhedra which connect principally by corner sharing, with a minority of edge sharing; a-$\mathrm{T}{\mathrm{a}}_{2}{\mathrm{O}}_{5}$ on the other hand has a local structure more akin to the crystalline polymorphs, built primarily from six- and seven-fold polyhedra, with a larger degree of edge sharing. The structural differences between liquid and amorphous $\mathrm{T}{\mathrm{a}}_{2}{\mathrm{O}}_{5}$, coupled with observations of increasing peak bond lengths upon cooling, are consistent with the interpretation that the amorphous film reaches a supercooled liquidlike metastable equilibrium during deposition. In other words, the amorphous film shares a common progenitor state with a hypothetical glass quenched from a fragile melt. In addition, we show that recent classical interatomic potentials do not fully reproduce the diffraction data, and infer that inclusion of attractive (non-Coulombic) Ta-Ta interactions is important, particularly for obtaining the correct degree of edge sharing, coordination numbers, and densities. Nanoscale inhomogeneity of the amorphous film is confirmed by the observation of small-angle x-ray scattering.

Comparing modeling predictions of aluminum edge dislocations: semidiscrete variational Peierls-Nabarro versus atomistics.

Multiple computational methods for modeling dislocations are implemented within a high-throughput calculation framework allowing for rigorous investigations comparing the methodologies. Focusing on aluminum edge dislocations, twenty-one classical aluminum interatomic potentials are used to directly model dislocation core structures using molecular dynamics, as well as provide input data for solving the semidiscrete variational Peierls-Nabarro dislocation model. The predicted dislocation core spreading obtained from both computational methods show similar trends across the potentials. Additionally, tests are done to rigorously determine if a recent correction to the Peierls-Nabarro model results in better agreement with the atomistic calculations.

Extending the accuracy of the SNAP interatomic potential form.

The Spectral Neighbor Analysis Potential (SNAP) is a classical interatomic potential that expresses the energy of each atom as a linear function of selected bispectrum components of the neighbor atoms. An extension of the SNAP form is proposed that includes quadratic terms in the bispectrum components. The extension is shown to provide a large increase in accuracy relative to the linear form, while incurring only a modest increase in computational cost. The mathematical structure of the quadratic SNAP form is similar to the embedded atom method (EAM), with the SNAP bispectrum components serving as counterparts to the two-body density functions in EAM. The effectiveness of the new form is demonstrated using an extensive set of training data for tantalum structures. Similar to artificial neural network potentials, the quadratic SNAP form requires substantially more training data in order to prevent overfitting. The quality of this new potential form is measured through a robust cross-validation analysis.