ChIMES: A Force Matched Potential with Explicit ... - ACS Publications

Subscriber access provided by READING UNIV

Article

ChIMES: A Force Matched Potential with Explicit Three-body Interactions for Molten Carbon Rebecca Kristine Lindsey, Laurence E. Fried, and Nir Goldman J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00867 • Publication Date (Web): 07 Nov 2017 Downloaded from http://pubs.acs.org on November 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

ChIMES: A Force Matched Potential with Explicit Three-body Interactions for Molten Carbon Rebecca K. Lindsey,∗ Laurence E. Fried, and Nir Goldman Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, United States E-mail: [email protected] Phone: +1-925-422-0915 Abstract We present a new force field and development scheme for atomistic simulations of materials under extreme conditions. These models, which explicitly include twoand three-body interactions, are generated by fitting linear combinations of Chebyshev polynomials through force matching to trajectories from Kohn–Sham density functional theory (DFT). We apply our method to liquid carbon near the diamond/graphite/liquid triple point and at higher densities and temperatures, where metallization and manybody effects may be substantial. We show that explicit inclusion of three-body interaction terms allows our model to yield improved descriptions of both dynamic and structural properties over previous empirical potential efforts, while exhibiting transferability to nearby state points. The simplicity of our functional form and subsequent efficiency of parameter determination allows for extension of DFT to experimental time and length scales while retaining most of its accuracy.

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

Experiments targeting equation of state data at extreme conditions (i.e. 10’s to 100’s of GPa and 1000’s of K) are highly challenging. The rapidly changing chemical environments and relatively short time scales (nanoseconds) associated with these experiments make features such as temperature and evolution of chemical species extremely difficult to observe, often leaving them to be deduced from thermodynamic arguments. Alternatively, computer simulations can be leveraged for design of experiments, interpretation of data, and microscopic insights. While the aforementioned time scales and associated lengths scales are relatively small for physical experiments, they can surpass the feasibility limits of atomistic simulations. For example, quantum simulation methods such as Kohn–Sham density functional theory (DFT) 1 are an attractive choice for reactive atomistic simulations as they can yield accurate chemistry, particularly under extreme conditions, where the choice of exchange-correlation functional is less important. However, DFT is generally limited to nanometer system sizes (e.g., hundreds of atoms) and time scales on the order of picoseconds. In contrast, experiments can probe nanosecond time scales and micron length scales, or larger. Semi-empirical quantum methods such as density functional tight binding (DFTB) 2 can yield significant improvements in computational efficiency by combining approximate quantum mechanics with empirical functions. However, DFTB still requires calculation of electronic state eigenvalues, which can limit usability to nanometer spatial scales. Thus, there remains a time and length scale gap between experimental needs and current capabilities for reactive molecular dynamics (MD) simulation. Molecular mechanics force fields (MMFF) are a promising tool for extending simulations to the time and length regimes relevant to chemistry at extreme conditions. 3 However, difficulty can arise under reactive conditions, where chemical bonds break and form rapidly, and classical force fields are unable to accurately model the subsequent dynamically changing potential energy surface. Furthermore, the MMFFs that are able to capture reactivity are typically comprised of numerous complicated non-linear equations, making parameteriza2


Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tion a significant endeavor, even when combined with advanced fitting approaches such as machine learning. For example, the Gaussian Approximation Potential (GAP) and Spectral Neighbor Analysis Potential (SNAP) machine learning models describe the local environment of an atom through order parameters corresponding to its bispectrum components, 4–6 while the ReaxFF, Long-Range Carbon Bond Order Potential (LCBOP), and Reactive Empirical Bond Order (REBO) models rely on a reactive bond-order type description, 3,7,8 and Embedded Atom Method (EAM)-style potentials leverage density-functional type descriptions for many-body interactions. 9 Additionally, machine learning approaches have been leveraged for the study of solid and liquid phases of carbon. 10,11 Recently, steps have been taken towards bridging the capability gap between the predictive power of DFT and the speed and scalability of MMFFs through use of force matching. 12 This fitting scheme, which fits interaction potentials to the 3N force components of a chemical system, has been heavily utilized for development of coarse-grained potentials from all-atom models and first-principles data, 13–17 as well as for development of atomic potentials from DFT data, 18–22 with some degree of success. Although force matching has been used successfully to improve semi-empirical quantum approaches, 23 it has been used for bound systems most typically, where molecular systems are treated as unreactive. Within the past year, efforts have been made towards development of force matched pairwise reactive interatomic models through linear combinations of Chebyshev polynomials. 24 In that work, reactive water models that retained the accuracy of DFT were developed for prediction of speciation and dynamics of a high density dissociate liquid state as well as the superionic phase at 2000 K, while yielding a roughly 5–6 order of magnitude improvement in efficiency compared to DFT. It was found, however, that inclusion of a non-linear overbonding potential was required to accurately describe the higher-bodied contributions to chemistry and dynamics. In this work, we build upon previous efforts and develop force matched models using linear combinations of Chebyshev polynomials for molten carbon. In contrast to previous

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

work, over-bonding terms are not included. Instead, explicit three-body interactions are introduced. These force fields are not confined to any particular shape, and thus, can fully adapt to changes in dominant physics (i.e. from insulating to metallic, or from non-reactive to reactive) that can occur over the broad range of temperature and pressure encompassed by the extreme condition regime. Hot compressed carbon is an ideal test case for this model due to its monoatomic nature, tendency to become metallic at sufficiently extreme conditions (yielding significant three-body interactions), and relevance to the geochemical and alternative energy communities. 25–28 Furthermore, the complete phase diagram of the material remains unknown despite its ubiquity, owing to the significant experimental challenges faced when exceeding the temperatures and pressures associated with melting of the diamond phase. 29–31 The models developed in this work are compared to results from DFT and two popular force fields for carbon under extreme conditions: LCBOP 7 and the 2002 parameterization of REBO. 8 Our approach is simpler than the GAP model 10 in terms of functional form complexity and training methodology, while achieving good accuracy in the liquid structure and dynamics.

2

Method and Computational Details

The model developed in the present work is referred to as the Chebyshev Interaction Model for Efficient Simulation (“ChIMES”), and is given by

EChIMES =

XNX i

Eij +

j>i

N X XX i

Eijk ,

(1)

j>i k>j

where N is the total number of atoms in the system and EChIMES , Eij , and Eijk are the total energy of the system, the energy between a pair of atoms i and j, and the energy between a set of three atoms, i, j, and k, respectively. Though the general form contains explicit terms for both two- and three-body interactions, a standalone pairwise additive model can be obtained by removing the latter terms. For clarity, we refer to models containing only 4


Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


two-body interactions as “ChIMES-2B,” and models containing both two-and three-body interactions as “ChIMES-3B.” The two-body energy is expressed in terms of a series of Chebyshev polynomials of the first kind, which has been shown elsewhere to be particularly well suited for molecular model development 24

Eij = fp rij + fs rij

O2 X n=1

(2)

cn Tn sij .

In the above summation, Tn sij is a Chebyshev polynomial of order O2 = n, which depends on a transformed pair distance, sij , defined in [−1,1] and the corresponding coefficient for each Tn sij is cn .

While the expansion in Chebyshev polynomials in sij is formally equivalent to an ex-

pansion in simple polynomials of sij , the use of Chebyshev polynomials has advantages when performing force matching. The orthogonality and oscillatory nature of the Chebyshev polynomials leads to improved convergence properties of the Chebyshev coefficients as the polynomial order is increased. Since the Chebyshev polynomial values are bounded in the range [−1, 1], convergence of the expansion can be determined by inspection of the polynomial coefficients. 32 Note that for simplicity, the above equation has been given for a monoatomic system; for a system with more than one atom type, the values of cn and sij would also depend on the ee

ee

elements, e, contained in the atom pair ij, i.e. cni j and siji j . Of the remaining functions, fp rij introduces a penalty to prevent rij values below the minimum distance sampled in the training trajectory (rc,in ) and is defined by

fp rij =

   A

p

rc,in + dp − rij

  0,

3

, if rij < rc,in + dp

,

(3)

otherwise

where Ap is the penalty prefactor, rc,in is the inner cutoff radius, and dp is the penalty initiation distance. 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

The remaining function, fs rij , ensures the potential goes smoothly to zero at the outer

cutoff and is given by

fs rij =

1−

rij rc,out

!3

,

(4)

where rc,out is the outer cutoff distance.

Mapping of the pair distance from rij to sij is achieved through the following

s rij =

2 exp

−rc,in −rc,out − exp − exp λ λ . −rc,in −rc,out − exp exp λ λ

−rij λ

(5)

The above equation assumes rij exists between the rc,in and rc,out , and through a Morseinspired transformation, 33–35 maps the inner and outer cutoff distances to the sij interval of [−1,1]. The Morse variables lead to a natural decrease in the interaction strength as distance is increased. Testing has shown that the quality of the force matching fit is insensitive to the Morse variable when a λ of between 0.75 and 1.5 Å is used. Building from two-body descriptions, the three-body energy is given as the product of Chebyshev polynomials for constituent atom pairs:

Eijk = fs rij fs (rik ) fs rjk

O3 O3 X O3 X X m=0 p=0 q=0

′

cmpq Tm sij Tp (sik ) Tq sjk .

(6)

Here, the single sum given for the two-body energy is now replaced with a triple sum over the orders, O3 , for the ij, ik, and jk polynomials and yields a single coefficient for each set of powers, cmpq . While the form of Equation 6 is simple and numerically efficient to evaluate, it is general enough to represent an arbitrary finite-ranged three-body interaction. Coefficients exhibit permutational invariance, such that for a set of three identical atoms, c112 = c121 = c211 , etc. The primed sum indicates that only terms for which two or more of the m, p, q indices are greater than zero are included, which guarantees that three bodies i, j, k enter into the Chebyshev polynomials. The expression for Eijk also contains fs smoothing functions for each constituent pair distance, but penalty functions are not included for three-body 6


Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interactions; instead, all penalties for rij < rc,in are handled by the two-body Ap terms. Comparing Equation 2 with Equation 6, it is clear that the method used in construction Eijk could easily be extended towards construction of any n-body interaction; for example, a four-body interaction would simply be described as a product over polynomials for the 6 constituent pairs, which is the subject of future work. Aside from the Chebyshev coefficients, ChIMES models contain 6 user-defined parameters, O2 , O3 , rc,in , rc,out , λ, Ap , and dp . For the models discussed in the present work, O2 , O3 , rc,out , and λ were chosen by scanning each parameter to find the smallest value that yielded a converged result. The value of rc,in is not considered adjustable, and is always set to the smallest pair distance observed for each interaction type, during the training trajectory. For the remaining parameters, Ap , and dp , values are chosen to be as small as possible to avoid contamination of the conserved quantity, while still being large enough to prevent pair distances smaller than rc,in . The training trajectories for our force matching procedure contain M frames of N atoms from DFT-MD simulations, separated by enough time to decorrelate the forces. Coefficients for the ChIMES force fields (i.e. cn and cmpq ) are obtained by minimizing the following objective function v u u RMSE = t

N M 3 2 1 XXX Fijk,DFT − Fijk,ChIMES{c} , 3M N i=1 j=1 k=1

(7)

where RMSE, and {c} are the root-mean-squared error and model coefficients, respectively. Fijk indicates the k th Cartesian component of the force on atom j in configuration i. The subscripts “ChIMES” and “DFT” indicate forces predicted from the present force-matched model and the DFT molecular dynamics (DFT-MD) training trajectory. A distinguishing feature of the ChIMES force field is its purely-linear dependence on fitted coefficients, which allows parameters to be directly obtained through least squares solution of the overdeter-

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

mined matrix equation: (8)

M c = F DFT .

In the above, force component indices {ijk} form the rows of the matrix M and the elements of the vector F DFT , and two and/or three body permutationally unique parameter indices {n, mpq} form the columns of the matrix M and the elements of the vector c. The matrix M is defined by: Mtu =

∂Ft,ChIMES{c} , ∂cu

(9)

where t represents a combined index over ijk, and u represents a combined index over the permutationally invariant coefficients n, mpq. There are P parameters in total, and Q force components in total, with P < Q. Equation 8 is solved through the standard technique of the reduced (skinny) singular value decomposition. 32 The singular value decomposition expresses M as a product of three matrices: M = U wV T , where U is a Q×P column-orthogonal matrix, w is a P ×P diagonal matrix, and V is a P × P orthogonal matrix. The least squares solution of Equation 8 is given by: 32 c=

X U (i) · F DFT i

wi

V(i) ,

(10)

where U (i) is the ith column of U , and V (i) is the ith column of V . In Equation 10, only terms with |wi | > ǫ|wmax | are considered. This is a standard regularization technique which controls the magnitude of the fitting parameters while preserving the accuracy of the fit. 32 We use ǫ = 10−5 in the work presented here. The singular value decomposition (SVD) typically takes 101 to 103 seconds on a single core to yield a unique set of parameters. The fast parameterization makes the present model ideal for a self-consistent refinement scheme, where configurations from ChIMES molecular dynamics simulations are used to refine the force matching fit. This is in contrast to models containing non-linear terms that must be solved through relatively slow iterative approaches that can become trapped in local minima. 8


Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The force matching procedure can easily be extended to include components of the DFT stress tensor as well as DFT forces, where only the potential energy contribution to the stress is considered. Weights must be applied to the stresses to compensate for the many more force components per frame (3N ) than stress components (3 diagonal for hydrostatic, and optionally, 3 off-diagonal for deviatoric). Unless otherwise stated, training trajectories consist of 26 frames based on 3–5 ps DFT-MD simulations of 256 carbon atoms in the canonical ensemble at 5000 K and 2.43 g/cm3 , and are generated in a self-consistent fashion. Further details on training trajectory generation and specification of user-defined parameters are provided in the supporting information. All DFT simulations are run with VASP 36–39 utilizing a planewave cutoff of 1000 eV, the Perdew-Burke-Ernzerhof generalized gradient approximation functional 40,41 (PBE), and projector-augmented wave pseudopotentials 42,43 with a 0.5 fs timestep. The models developed in the present work are validated and compared with other reactive carbon force fields via molecular dynamics (MD) simulations. All MD results presented herein are from N V T ensemble simulations using a global Nose-Hoover thermostat 44,45 (i.e. on all atoms), and periodic boundary conditions are enforced through use of the explicit image method (i.e. ghost atoms). Simulations involving ChIMES potentials are run using locally developed MD software, while simulations for remaining models are run with the LAMMPS suite. 46 A 0.25 fs time step was used for all ChIMES simulations, whereas a value of 0.5 fs was used for LCBOP and REBO, and unless otherwise specified, all simulations/models are fit and/or run at 5000 K and 2.43 g/cm3 .

3

Results and discussion

We begin by comparing reproduction of DFT forces by the 17-parameter ChIMES-2B and 48-parameter ChIMES-3B potentials, for the same set of 26 configurations. The 26 frame training set was found sufficient to converge the root-mean-squared errors (RMSE), and

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

inclusion of additional frames was not found to lead to substantial improvements in the resulting model. In the present work, all ChIMES models use a 2-body order of 12 and when included, a 3-body order of 4. Remaining parameters for each of the models developed in this work, including user-specified values, are provided in the supporting information. The results in Figure 1 show that both models exhibit slight over prediction of small-magnitude DFT forces. The reduction in peak height observed for ChIMES-3B, however, indicates clear enhancement over ChIMES-2B model. A quantitative description of force recovery can be obtained by examining the RMSEs between forces stemming from DFT and from each force field. Here, we present these values as well as reduced RMSEs, which have been divided by the average of the absolute value of DFT forces, i.e. RMSEr = RMSE/h| FDFT |i, to account for the inherent relationship between RMSE and the magnitude of forces observed in the training set (i.e. an equivalent percent deviation for two sets of data of different magnitudes will yield much different RMSE values). We find that RMSE (RMSEr ) is reduced from 40.7 kcal/mol/Å (0.678) to 26.4 (0.439) upon introduction of three-body terms, further demonstrating improved performance of the ChIMES-3B model relative to ChIMES-2B. As an alternative metric of force reproduction, Figure 2 provides a direct comparison of the forces from DFT for each training configuration with those coming from either the ChIMES-2B or ChIMES-3B models. Once again, both models are found to exhibit good reproduction of DFT forces (R2 of 0.85 and 0.95 for ChIMES-2B and ChIMES-3B, respectively). Regarding the efficiency of the ChIMES models, we find that, for the present system, DFT takes roughly 50,000 CPU hours to yield 5 ps of N V T -dynamics, while the ChIMES-2B and -3B models require approximately 0.15 and 5 CPU hours, respectively. Figure 3 shows that the enhancement in model performance upon introduction of threebody terms is also reflected by improved accuracy of the radial distribution functions (RDFs) predicted by the ChIMES models, relative to DFT. Specifically, there is a clear reduction in first peak height and an improvement in location of the second, when moving from ChIMES2B to ChIMES-3B. The height of the first peak in an RDF can be used to compute the

10


Page 10 of 32

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


potential of mean force associated with removing an atom from its first solvation shell, ∆W 1

st solv

= −kB T ln[∆g(r = r1 st

st peak

)]. Thus, the difference in this value between DFT and st

st

1 solv 1 solv 1 solv the various models, ∆WDFT→FF = ∆WDFT − ∆WFF provides a numeric measure of how

well the force field methods reproduce the DFT RDFs. Reduction of the first peak height in the RDF from 2.92 to 2.23 is observed when moving from the ChIMES-2B to ChIMES-3B st

1 solv model, improving ∆WDFT→FF from −2.68 to 0.13 kcal/mol, respectively. The significant

reduction in height of the first RDF peak observed upon introduction of three-body terms can be understood through examination of Figure 4, which shows the likelihood of finding a set of 3-atoms with given ij, ik, and jk distances. The metallic nature of the system and subsequent repulsive effects give rise to disallowed three-body configurations (i.e. in the vicinity of rij ≤ 1.75; rik = rjk = 1.5 Å). As can be seen in the ChIMES-2B histogram, this model is unable to detect these configurations, leading to a large well feature in the histogram (indicated by a the blue arrow) and consequently increased first RDF peak height. Though there is no expectation that the present parameterization should yield good predictions for graphite and diamond structures, 0 K lattice constants were computed. Given the small value for our cutoff radius (3.15 Å), for starting structures with graphite layer spacing near the experimental value (i.e., 3.354 Å), the stopping criteria for that dimension is immediately satisfied. Hence, we choose to report the value of the a lattice constant for graphite, only. Results for the ChIMES models are provided in Table 1 along with experimental values. Good agreement with experimental diamond values is found for both ChIMES models, and for graphite, ChIMES-2B yields a slightly closer match to experiment than ChIMES-3B. Structural prediction by the ChIMES models is also benchmarked against two reactive carbon potentials, LCBOP 7 and the 2002 parameterization of REBO. 8 Both of these models use a bond order description for greater than two-bodied interactions, and are parameterized to reproduce features of graphite and diamond such as binding energies and shear elastic constants. Though these force fields have been successfully utilized in applications including prediction of the graphite/liquid and graphite/diamond melt lines, 47 they are designed to

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

capture the three- and four-fold coordination associated with the structures of graphite and diamond, while liquid carbon has been shown to also exhibit linear configurations. 48 Figure 3 st

1 solv indicates that both models yield substantially larger magnitude ∆WDFT→FF than ChIMES-

3B, at −3.55 and 1.95 kcal/mol for LCBOP and REBO, respectively. Furthermore, the force fields predict much higher structural order relative to both DFT and ChIMES-3B, evidenced by taller, narrower peaks in conjunction with the nearly-zero value of the first minima of both RDFs. We note that the agreement between the ChIMES-3B model and DFT at 5000 K, 2.43 g/cm3 is similar to that achieved by the GAP approach 10 with two-, three-, and many-body interactions at 5000 K and 2.5 g/cm3 . Predicted pressures for each model are presented in Table 2. Of the models, ChIMES3B provides the closest prediction, over-estimating the DFT result by 43%. The LCBOP model yields the next best result, with a 47% over-estimation, followed by REBO with a 81% under-estimation, and ChIMES-2B with a 97% over-prediction, respectively. None of the models predict an accurate equation of state at the conditions studied. In the case of ChIMES-3B, this indicates that many-body forces play an important role in determining the DFT pressure, as it has been shown necessary for other monoatomic systems at similar densities. 49 It is unsurprising that the ChIMES-3B model exhibits better predictions compared to ChIMES-2B, as the latter model is unable to describe potential three-body contributions to the pressure. It is important to recall, however, that matching to just the forces (rather than forces and stress tensors) does not constrain the energy dependence on the system volume, which can also contribute to inaccurate pressure predictions. 23,24 In applications where better pressure prediction is desired, weighted inclusion of stress tensor components in the training data has been shown to be effective. 23 Doing so for an additional two-body only fit, ChIMES2B+σ, results in a pressure of 19.2 GPa, consistent with the ChIMES-3B result. Furthermore, including stress tensor data in the fit leads to only small changes in the structure and dynamics of ChIMES-2B+σ compared to ChIMES-2B. In contrast, including tensor data in

12


Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a fit with three-body interactions does not yield changes in pressure for modest weighting (i.e. 19.7 and 20.0 GPa for w = 1.0 and 100.0, respectively), whereas models exhibit singularities at higher weight (i.e. w > 500.0). The inability to recover the DFT pressure, even with inclusion of stress tensor data, may be due to a fundamental issue; the present models are very short-ranged, while pressure in condensed-phase systems is typically determined by long-ranged interactions. 50 Further comparison of properties predicted between the two models can be found in the supporting information. We now present results for time-dependent properties of the carbon system. Figure 5 provides the vibrational power spectra along with predicted diffusion constants for ChIMES3B compared to DFT, ChIMES-2B, LCBOP, and REBO. Beginning with comparison of the ChIMES models, both capture the overall shape of the DFT power spectrum, however ChIMES-3B yields a substantially improved result. In particular, the relative peak intensity of the spectrum near 100–400 cm−1 is reduced from approximately 21.5 to 17.3, for ChIMES2B and ChIMES-3B, respectively, compared to DFT at roughly 15.0. Furthermore, the values at 0 cm−1 , which are directly related to system diffusivity, are reduced from 23.8 to 12.8 upon addition of three-body interactions, compared to 9.8 for DFT. In contrast, neither LCBOP or REBO capture the shape of the DFT spectrum, with both exhibiting a value of nearly zero at 0 cm−1 , indicative of little to no atomic diffusion. Notably, both LCBOP and REBO have a similar shape at frequencies below 750 cm−1 , however REBO exhibits a distinct second peak just above 1500 cm−1 . The trends in power spectra at 0 cm−1 for each of the models relative to DFT are also reflected in predicted self-diffusion constants, ds . Relative to DFT, only the ChIMES-2B and ChIMES-3B models yield values of the same order of magnitude, with ChIMES-3B exhibiting the closest result. The REBO model predicts ds of only one order of magnitude less than DFT while LCBOP underestimates the DFT result by several orders. Finally, we investigate transferability of our ChIMES-3B model. Since these models are fit to a single state point, there is no a priori expectation that they should exhibit

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

transferability over large temperature or pressure ranges; nevertheless there is utility in understanding model limitations. The ChIMES force fields discussed above are fit close to the carbon triple point. Because we want to investigate isothermal transferability for liquid carbon rather than across a phase boundary, the present studies are conducted at 6000 K, safely above the predicted diamond/liquid melt line for the temperature and density conditions examined. 30,31,47 Structure along an isotherm at 6000 K, and densities of 2.25, 2.75, and 3.00 g/cm3 was determined from DFT simulations, and a new ChIMES-3B model including both two- and three-body interactions, fit only to 3.00 g/cm3 , was run at each state point. The results, shown in Figure 6, demonstrate excellent agreement between DFT and the force field at each simulation condition, both in terms of overall RDF shape and st

1 solv in ∆WDFT→FF , with values ranging from −0.25 to 0.14 kcal/mol across the three sets of

simulations. Furthermore, Table 3, which provides ds and pressure values predicted by DFT and the present ChIMES-3B model at each of the three state points, shows that the model yields pressures within the range of deviations observed for the system at 5000 K, for each state point, and that accuracy in dynamical properties is also maintained. Thus, the results from the present study demonstrates that, without any additional tuning, models can indeed provide similar performance at nearby state points. We note that the ChIMES-2B force field yields poor transferability relative to ChIMES-3B; while the model performs reasonably well at the state point for which it was fit (i.e. 6000 K and 3.0 g/mL), performance at other state points is questionable.

4

Conclusions

New force matched models have been developed for hot compressed carbon using a purely two-bodied description as well as with explicit inclusion of three-body terms. The ChIMES3B force field was shown to improve upon ChIMES-2B for both structural and dynamic properties, and to outperform existing reactive carbon models LCBOP and REBO. In par-

14


Page 14 of 32

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ticular, it was shown that inclusion of three-body terms was necessary to recover the correct structure of the carbon system, and leads to significant improvements in the predicted distribution of forces, vibrational power spectrum, and diffusion coefficient. Furthermore, ChIMES-3B has been shown to exhibit some degree of transferability for a given temperature over a finite range of densities, without any additional tuning. The ChIMES models are entirely linear in fitted coefficients, allowing for rapid generation of parameter sets through simple linear least squares regression. Furthermore, these models are completely flexible in shape, allowing the tunability required to accurately recover system physics described by DFT. The ChIMES model may require additional higher-body interaction terms or longer outer cutoff values to generate pressures closer to those predicted by DFT. However, overall, our ChIMES potentials have been demonstrated capable of maintaining DFT-level accuracy for prediction of various properties, and thus have potential use for studying complex long time and length scale reactivity at extreme conditions, where there is a significant need for computationally efficient atomistic simulations methods to aid in the interpretation and design of experiments. Current efforts are underway to apply our methodology to the complex reactivity of C-H-O-N containing materials under extreme pressures and temperatures.

Acknowledgement This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The project 16-LW-020 was funded by the Laboratory Directed Research and Development Program at LLNL with N.G. as principal investigator, and is assigned release number LLNL-JRNL-736850.

Supporting Information Available The following files are available free of charge.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

• ChIMES-SI.pdf: Parameters and additional details for the four ChIMES models developed in this work This material is available free of charge via the Internet at http://pubs.acs.org/.

16


Page 16 of 32

Page 17 of 32

-3

8×10

Normalized Distribution

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-3

6×10

4×10

2×10

DFT ChIMES-3B ChIMES-2B

-3

-3

0 0

200 400 |F| [kcal/mol/Å]

600

Figure 1: Normalized distribution of forces acting on carbon atoms at 5000 K and 2.43 g/cm3 .

17



600

2

ChIMES-2B, R = 0.85, m = 0.73 2 ChIMES-3B, R = 0.95, m = 0.89

400 FMM [kcal/mol/Å]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

200 0 -200 -400 -600 -600 -400 -200 0 200 400 FDFT [kcal/mol/Å]

600

Figure 2: Comparison of forces for carbon at 5000 K and 2.43 g/cm3 predicted by DFT, ChIMES-2B, and ChIMES-3B. The black line gives y = x and serves as a guide to the eye, and the legend provides R2 values and slopes (m) for linear regressions to each data set.

18


Page 18 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DFT ChIMES-3B ChIMES-2B LCBOP REBO

3.0

g(r)

Page 19 of 32

2.0

1.0

0.0 1.0

2.0

3.0 rCC [Å]

4.0

25

Figure 3: Radial distribution function for carbon at 5000 K and 2.43 g/cm3 .

19


rjk [Å] rjk [Å]

Page 20 of 32

0.010 0.008 0.006 0.004 0.002 0.000

rjk [Å]

ChIMES-2B ChIMES-3B

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

DFT


Figure 4: Histogram slices of atom triplet distances for DFT (top) ChIMES-2B (middle) and ChIMES-3B (bottom). Contour lines are drawn at intervals of 0.001, and data are normalized for r ≤ 3.15 Å. The blue arrow indicates regions in three-body space disallowed by DFT and ChIMES-3B, where ChIMES-2B predicts significant population.

20


Page 21 of 32

rCC [Å] 25 DFT ChIMES-3B ChIMES-2B LCBOP REBO

Intensity (Arb. units)

20

15

10

5

ChIMES-3B: 2.30E-4

10–3

3000

ChIMES-2B: 4.96E-4

REBO : 1.25E-5

100

1000 2000 -1 Frequency [cm ] LCBOP: 4.10E-7

Self-diffusion constant [cm2/s]

0 0

DFT: 1.67E-4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10–5

Figure 5: Vibrational power spectra (top) and self-diffusion coefficients (bottom) for carbon at 5000 K and 2.43 g/cm3 .

21



DFT 3 ChIMES-3B fit at 3.00 g/cm

g(r), ρ = 2.25 g/cm

3

2

1

0

g(r), ρ = 2.75 g/cm

3

2

1

0

3

2 g(r); ρ = 3.00 g/cm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

0

1

2

3 r [Å]

4

5

Figure 6: Radial distribution functions for simulations at 6000 K and 2.25 (top), 2.75 (middle), and 3.00 g/cm3 (bottom), using a force field fit at 6000 K and 3.00 g/cm3 . 22


Page 22 of 32

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Lattice constants at zero temperature [Å]. Experimental values for graphite and diamond are taken from references 51 and 52, respectively. diamond graphite

ChIMES-2B a 3.567 a 2.550

23

ChIMES-3B 3.565 2.500


Experiment 3.57 2.562


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Predicted pressures for carbon at 5000 K [GPa]. DFT ChIMES-3B ChIMES-2B ChIMES-2B+σ LCBOP REBO

24

13.8 19.7 27.2 19.2 20.3 2.6


Page 24 of 32

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Self-diffusion coefficients and predicted pressures for carbon at 6000 K. ρ [g/cm3 ] 2.25 2.75 3.00

ds [cm2 /s] DFT ChIMES-3B −4 3.5 × 10 5.8 × 10−4 2.7 × 10−4 4.4 × 10−4 1.5 × 10−4 3.8 × 10−4

25

DFT 12.3 30.1 40.1


p [GPa] ChIMES-3B 20.9 37.6 50.0


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. Physical Review 1965, 140, A1133. (2) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Fraunheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex material properties. Physical Review B 1998, 58, 7260. (3) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: a reactive force field for hydrocarbons. The Journal of Physical Chemistry A 2001, 105, 939. (4) Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Physical review letters 2010, 104, 136403. (5) Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. Physical Review B 2013, 87, 184115. (6) Thompson, A. P.; Swiler, L. P.; Trott, C. R.; Foiles, S. M.; Tucker, G. J. Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials. Journal of Computational Physics 2015, 285, 316. (7) Los, J.; Fasolino, A. Intrinsic long-range bond-order potential for carbon: performance in Monte Carlo simulations of graphitization. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 024107. (8) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J. Phys. Condens. Matter 2002, 14, 783. (9) Lipnitskii, A.; Saveliev, V. Development of n-body expansion interatomic potentials and its application for V. Computational Materials Science 2016, 121, 67–78. 26


Page 26 of 32

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Deringer, V. L.; Csányi, G. Machine learning based interatomic potential for amorphous carbon. Physical Review B 2017, 95, 094203. (11) Khaliullin, R. Z.; Eshet, H.; Kühne, T. D.; Behler, J.; Parrinello, M. Graphite-diamond phase coexistence study employing a neural-network mapping of the ab initio potential energy surface. Physical Review B 2010, 81, 100103. (12) Ercolessi, F.; Adams, J. B. Interatomic potentials from first-principles calculations: the force-matching method. Europhys. Lett. 1994, 26, 583. (13) Izvekov, S.; Voth, G. A. Multiscale coarse graining of liquid-state systems. J. Chem. Phys. 2005, 123, 134105. (14) Izvekov, S.; Violi, A.; Voth, G. A. Systematic coarse-graining of nanoparticle interactions in molecular dynamics simulation. J. Phys. Chem. B 2005, 109, 17019–17024. (15) Wang, Y.; Izvekov, S.; Yan, T.; Voth, G. A. Multiscale coarse-graining of ionic liquids. J. Phys. Chem. B 2006, 110, 3564–3575. (16) Zhou, J.; Thorpe, I. F.; Izvekov, S.; Voth, G. A. Coarse-grained peptide modeling using a systematic multiscale approach. Biophys. J. 2007, 92, 4289–4303. (17) Liu, P.; Izvekov, S.; Voth, G. A. Multiscale coarse-graining of monosaccharides. J. Phys. Chem. B 2007, 111, 11566–11575. (18) Lenosky, T. J.; Sadigh, B.; Alonso, E.; Bulatov, V. V.; de la Rubia, T. D.; Kim, J.; Voter, A. F.; Kress, J. D. Highly optimized empirical potential model of silicon. Model. Simul. Mater. Sci. Eng. 2000, 8, 825. (19) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: a new method for force-matching. J. Chem. Phys. 2004, 120, 10896–10913.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(20) Izvekov, S.; Voth, G. A. Effective force field for liquid hydrogen fluoride from ab initio molecular dynamics simulation using the force-matching method. J. Phys. Chem. B 2005, 109, 6573–6586. (21) Paesani, F.; Zhang, W.; Case, D. A.; Cheatham III, T. E.; Voth, G. A. An accurate and simple quantum model for liquid water. J. Chem. Phys. 2006, 125, 184507. (22) Spiegel, K.; Magistrato, A.; Maurer, P.; Ruggerone, P.; Rothlisberger, U.; Carloni, P.; Reedijk, J.; Klein, M. L. Parameterization of azole-bridged dinuclear platinum anticancer drugs via a QM/MM force matching procedure. J. Comput. Chem. 2008, 29, 38–49. (23) Goldman, N.; Fried, L. E.; Koziol, L. Using Force-matched potentials to improve the accuracy of density functional tight binding for reactive conditions. J. Chem. Theory Comput. 2015, 11, 4530. (24) Koziol, L.; Fried, L. E.; Goldman, N. Using force-matching to determine reactive force fields for bulk water under extreme thermodynamic conditions. J. Chem. Theory Comput. 2017, 13, 135. (25) Bundy, F.; Bassett, W.; Weathers, M.; Hemley, R.; Mao, H.; Goncharov, A. The pressure-temperature phase and transformation diagram for carbon; updated through 1994. Carbon 1996, 34, 141–153. (26) Ross, M. The ice layer in Uranus and Neptune–diamonds in the sky? Nature 1981, 292, 435–436. (27) MacKinnon, A. J.; Meezan, N. B.; Ross, J. S.; Pape, S. L.; Hopkins, L. B.; Divol, L.; Ho, D.; Milovich, J.; Pak, A.; Ralph, J.; Döppner, T.; Patel, P. K.; Thomas, C.; Tommasini, R.; Haan, S.; MacPhee, A. G.; McNaney, J.; Caggiano, J.; Hatarik, R.; Bionta, R.; Ma, T.; Spears, B.; Rygg, J. R.; Benedetti, L. R.; Town, R. P. J.;

28


Page 28 of 32

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Bradley, D. K.; Dewald, E. L.; Fittinghoff, D.; Jones, O. S.; Robey, H. R.; Moody, J. D.; Khan, S.; Callahan, D. A.; Hamza, A.; Biener, J.; Celliers, P. M.; Braun, D. G.; Erskine, D. J.; Prisbrey, S. T.; Wallace, R. J.; Kozioziemski, B.; Dylla-Spears, R.; Sater, J.; Collins, G.; Storm, E.; Hsing, W.; Landen, O.; Atherton, J. L.; Lindl, J. D.; Edwards, M. J.; Frenje, J. A.; Gatu-Johnson, M.; Li, C. K.; Petrasso, R.; Rinderknecht, H.; Rosenberg, M.; guin, F. H. S.; Zylstra, A.; Knauer, J. P.; Grim, G.; Guler, N.; Merrill, F.; Olson, R.; Kyrala, G. A.; Kilkenny, J. D.; Nikroo, A.; Moreno, K.; Hoover, D. E.; Wild, C.; Werner, E. High-density carbon ablator experiments on the National Ignition Facility. Phys. Plasmas 2014, 21, 056318. (28) Eggert, J.; Hicks, D.; Celliers, P.; Bradley, D.; McWilliams, R.; Jeanloz, R.; Miller, J.; Boehly, T.; Collins, G. Melting temperature of diamond at ultrahigh pressure. Nat. Phys. 2010, 6, 40–43. (29) Gold, J. S.; Bassett, W. A.; Weathers, M. S.; Bird, J. M. Melting of diamond. Science 1984, 225, 921–922. (30) Wang, X.; Scandolo, S.; Car, R. Carbon phase diagram from ab initio molecular dynamics. Phys. Rev. Lett. 2005, 95, 185701. (31) Sun, J.; Klug, D. D.; Martoňák, R. Structural transformations in carbon under extreme pressure: beyond diamond. J. Chem. Phys. 2009, 130, 194512. (32) Press, W. H.; Tuekolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: Cambridge, U.K., 2007. (33) Wang, Y.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Full-dimensional, ab initio potential energy and dipole moment surfaces for water. J. Chem. Phys. 2009, 131, 054511.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(34) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 094509. (35) Braams, B. J.; Bowman, J. M. Permutationally invariant potential energy surfaces in high dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577–606. (36) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558. (37) Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal– amorphous-semiconductor transition in germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251. (38) Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15–50. (39) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169. (40) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (41) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple [Erratum to Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396– 1396. (42) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (43) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmentedwave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758. 30


Page 30 of 32

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511–519. (45) Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695. (46) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Phys. Chem. C 1995, 117, 1–19. (47) Ghiringhelli, L. M.; Meijer, E. J. Computer-Based Modeling of Novel Carbon Systems and Their Properties; Springer, 2010; pp 1–36. (48) Cannella, C. B.; Goldman, N. Carbyne fiber synthesis within evaporating metallic liquid carbon. J. Phys. Chem. C 2015, 119, 21605–21611. (49) Schwerdtfeger, P.; Steenbergen, K. G.; Pahl, E. Relativistic coupled-cluster and densityfunctional studies of argon at high pressure. Physical Review B 2017, 95, 214116. (50) Rosenberger, D.; Hanke, M.; van der Vegt, N. F. Comparison of iterative inverse coarsegraining methods. The European Physical Journal Special Topics 2016, 225, 1323–1345. (51) Zhao, Y. X.; Spain, I. L. X-ray diffraction data for graphite to 20 GPa. Physical Review B 1989, 40, 993. (52) Yamanaka, T.; Morimoto, S.; Kanda, H. Influence of the isotope ratio on the lattice constant of diamond. Physical Review B 1994, 49, 9341.

31



Graphical TOC Entry E(rijk) = ! ! ! cmpq Tm(sij) Tp(sik) Tq(sjk) c1

E(rijk) with fixed rij

!

c3

M

c4

rjk

c2

T (s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

. . .

rik

s

32


Page 32 of 32

ChIMES: A Force Matched Potential with Explicit ... - ACS Publications

Recommend Documents