A Dispersion-Corrected Modified Embedded-Atom Method Bond

Subscriber access provided by Kaohsiung Medical University

A: New Tools and Methods in Experiment and Theory

A Dispersion-Corrected Modified Embedded-Atom Method Bond Order Interatomic Potential for Sulfur Doyl Dickel, Steven R. Gwaltney, Sungkwang Mun, Michael I. Baskes, and Mark F. Horstemeyer J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b07410 • Publication Date (Web): 06 Nov 2018 Downloaded from http://pubs.acs.org on November 21, 2018

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 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 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.

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 19 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 Journal of Physical Chemistry

A Dispersion-Corrected Modified Embedded-Atom Method Bond Order Interatomic Potential for Sulfur Doyl Dickel,∗,† Steven R. Gwaltney,‡ Sungkwang Mun,† Michael I. Baskes,¶,§,k,⊥ and Mark F. Horstemeyer†,# †Center for Advanced Vehicular Systems, Mississippi State University, Starkville, MS 39759 ‡Department of Chemistry, Mississippi State University, Starkville, MS 39759 ¶Office of Research and Development, Mississippi State University, Starkville, MS 39759 §Los Alamos National Laboratory, Los Alamos, NM 87545 kJacobs School of Engineering, University of California at San Diego, La Jolla, CA 92093 ⊥Department of Materials Science and Engineering, University of North Texas, Denton, TX 76203 #Department of Mechanical Engineering, Mississippi State University, Starkville, MS 39759 E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

Abstract An interatomic potential for sulfur has been developed using the bond order addition to the modified embedded-atom method (MEAM-BO). In order to correctly model the interaction between molecules, dispersion forces have been included via the DFT-D3 modification. It is demonstrated that this semi-empirical classical potential correctly reproduces the behavior of the S2 dimer, various cyclic sulfur rings, the molecular solids: α-, β-, and γ-, sulfur, and a number of theoretical, high symmetry sulfur structures. This potential will serve as a useful tool in the atomistic modeling of sulfur, and, ultimately, in the modeling of sulfur containing organic compounds using this updated MEAM-BO formalism.

Introduction Sulfur, and in particular thiol functional groups, play an important role in biological chemistry. The formation of disulfide bridges is critical to the tertiary structure of protein chains. Also, sulfhydryl groups can contribute to catalytic activity in enzymes. 1 However, to date, most computational research on the atomic scale behavior of sulfur has been performed using ab initio techniques which are reliable but computationally demanding, precluding the investigation of systems of more than a few dozen atoms. A few specialized potentials, using formalisms such as ReaxFF 2,3 and REBO, 4 have been developed for sulfur as part of various systems. The recently developed Modified Embedded Atom Method with Bond Order (MEAM-BO) interatomic potential formalism 5 allows for the first time the generation of semi-empirical classical potentials consistent across metallic and covalent bonding. The success of this method in generating a potential for carbon-hydrogen systems which correctly describes the bonding character and energies of single, double, and triple covalent bonds using a formalism which has already been extensively used to describe metals 6,7 suggests that it will be capable of successfully modeling interactions among a variety of different materials in composite 2


Page 2 of 19

Page 3 of 19 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


systems, notably between organic compounds and metals. As potentials are developed and applied to other elements typical of organic molecules, such as oxygen, nitrogen, phosphorus, and sulfur, the range of molecules and interactions which can be consistently analyzed will grow significantly. Correctly modeling the atomistic behavior of sulfur also requires modeling of dispersion interactions. 8 These interactions are not captured by most Density Functional Theory (DFT) methods but are crucial for correctly describing the significant long-range interactions in many covalently bonded materials, notably including the solid phases of sulfur. We present here a MEAM-BO potential for sulfur including a DFT-D3 correction 8,9 which correctly reproduces the properties of a number of sulfur structures. These include the S2 dimer, a number of ringed structures such as the cyclic S8 molecule, sulfur chains, and a number of high coordination structures not observed in nature but for which first principles data can be generated. Due to the design of the MEAM-BO formalism, the developed potential reproduces, by construction, the experimental results for the bonding of S2 and the binding energy and bond length of S8 . Combined with the inclusion of the dispersion correction in the form of the D3 potential, the potential should offer accuracy beyond what is currently available from other classical potentials at both short and long range. The potential is usable as a reactive potential to describe the behavior of sulfur systems and can also serve as the basis for an extension to more complicated organic systems including the existing C-H MEAM-BO potential and future potentials which incorporate a wide variety of both biological and metallic elements.

Theory We first briefly summarize both the MEAM and MEAM-BO potential formalisms. More in depth discussion of both can be found in 7 and 5 respectively. In MEAM, the total energy of a system is given by

3



E=

X

Fτi (¯ ρi ) +

i

Page 4 of 19

1X Sij Φτi τj (Rij ) 2 i,j

(1)

where the embedding function Fτ , for an atom type τ , represents the energy cost to insert an atom i at a site with background electron density ρ¯i , which is a combination of spherically symmetric and angular components. Φ is a pairwise interaction whose form is determined by requiring that the energy versus volume curve of a particular reference structure follow the Universal Equation of State (UEOS) of Rose et al. 10 with respect to nearest neighbor distance, Rij . Sij is a screening function which considers three body interactions. In order to correctly distinguish bond order as is necessary in the evaluation of covalent bonds and organic molecules, MEAM has been modified to include a bond order term (MEAM-BO). In this formulation, each bond I, is assigned a fraction of double and triple bond character, denoted fI2 and fI3 respectively. For a perfect double bond, fI2 = 1 and fI3 = 0, while for a perfect triple bond, fI2 = 0 and fI3 = 1. This modifies the standard total energy from MEAM as follows:

EM EAM −BO = EM EAM +

X

fI3 · E3 (RI ) +

I

X

fI2 · E2 (RI , BOI )

(2)

I

where EBC (BC = 2 for a double bond, BC = 3 for a triple bond) is an additional energy contribution that depends on the bond length and the bond order. For sulfur, triple bonds are not observed in nature, so the relevant parameters are set to zero. In addition, the original MEAM-BO equation (initially developed for hydrocarbons) has been slightly modified to handle more general cases including sulfur. To show the change, we consider the following energy fraction term for double bond in Eq. 2.

4


Page 5 of 19 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


h h i2 i2 (0) (0) (0) (0) (0) fI2 = Si1 i2 · Zi1 − Z2 · D2 Zi2 − Z2 h h h i2 i2 i2 (1) (3) (3) (1) (1) (1) (1) (1) − Z2 · D2 ZI − Z2 · D2 Zi2 ·D2 Z i1 h i2 (3) (3) +fI3 · 1 − D3 ZI (0) D2

(3)

(0)

(1)

where Si1 i2 is the angular screening function for the two atoms in the bond. Zi , Zi

and

(3)

ZI are partial density functions for an atom i or a bond I to satisfy the geometric condition (h)

of a reference structure. Dbond is a Gaussian type function where the peak of the curve is one and the position of the center of the peak is at zero. The last part of the equation captures the transition from triple to double bond which is irrelevant in this study. The change has h i2 (1) (1) (1) (1) − Z2 where a parameter Z2 is added as Zi been made to all occurrences of D2 the geometric condition is satisfied if Z (0) ∼ 2, Z (1) ∼ 1, and Z (3) ∼ 0 in sulfur case. See Mun et. al. 5 for more detail. Next, we consider the DFT-D3 correction, accounting for dispersion forces, necessary to properly calculate the long-range interaction between molecules. In previous work using the MEAM-BO formalism, a van der Waals type interaction was used to model the long range interaction. We have here instead used the DFT-D3 correction with BJ damping. 9 Most currently used density functionals, being semilocal, cannot correctly describe these long-range interactions. However, a number of corrections to account for this interaction have been developed. 8,9 The DFT-D3 method adds a pair interaction and a three-body interaction term to the total energy with coefficients which have been fit to correctly model the long-range dispersion interaction between various atomic species. We use only the pair potential for the DFT-D3 parameters fit for the PBE pseudopotential set 9 and the dispersion interaction is smoothly cut off at distances greater than 10 ˚ A. This insured the important region of the interaction was considered but that the number of relevant neighbors was not so large as to become computationally inefficient. This method was applied directly as an 5



Page 6 of 19

energy correction to the MEAM results. In order for MEAM to still produce the correct short-range interaction, the UEOS used by the MEAM potential is modified to account for the short-range contribution of the DFT-D3 interaction:

EEOS = −Ec0 (1 + a∗ + δ Ec0 = Ec +

re ∗3 −a∗ a )e r

Z 0 E 2

E 0 = E D3 (re ) r Ec0 r a∗ = α ( − 1) Ec re

(4) (5) (6) (7)

where Ec and re are the cohesive energy and nearest neighbor distance of the reference state of the element, E 0 is the dispersion energy at re , Ec0 is the corrected cohesive energy, Z is the number of nearest neighbors in the MEAM reference structure, α and δ are the standard MEAM parameters used in the UEOS (see, for example 11 where δ is divided into δ + and δ − depending on the sign of

r re

− 1) and a∗ is the scale parameter corrected to account for the

D3 interaction. Here, the energy correction used was E 0 = 0.07368 eV. In the absence of this correction and without the DFT-D3 interaction, the original MEAM-BO form is recovered. As will be shown, the addition of this interaction allows the sulfur potential to correctly predict the long range interaction between sulfur molecules.

Sulfur Potential The parameters used in the MEAM potential for sulfur and its bond order modification are given in Tables 1 and 2 respectively. As required for the MEAM potential formalism, a reference structure must be selected to determine the pairwise interaction. Typically, a low-energy structure, similar to those observed in nature, is selected for a potential. Because of this, high coordination structures such as fcc, bcc, or hcp, typical for metals, would not be suitable for sulfur. The two simplest singly bonded, low-energy, naturally occurring sul6


Page 7 of 19 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


fur structures are the ringed sulfur molecules and sulfur chains. Because of its low energy, octasulfur, an eight membered sulfur ring was chosen to match the reference state. However, if only first nearest neighbors are considered, all ring and chain structures can be described with only two parameters: bond length and angle. As such, these are all equivalent, considering only nearest neighbor interactions, to a zigzag chain. As such, a zigzag structure with specified bond length 2.048 ˚ A and angle 108◦ was chosen as the reference state. It is important to note, however, that for chain structures, the interaction of lone pairs of electrons on each sulfur atom cause the chains to develop a chirality. Since the MEAM implementation cannot account directly for the presence of lone pairs, it will be insensitive to this effect directly and can only approximate it through interactions between atoms. As mentioned, previously, triple bond terms may be ignored for sulfur as, due to the six atoms in its valence shell, it will form, at most, double bonds. In addition, no hybridized bonding, such as is observed in carbon for benzene or graphene, need be considered. The parameters for the single covalent bond were determined using the properties of the S8 ring structure. The parameters for the sulfur-sulfur double bond were determined from experimental data on the S2 dimer. By construction, the bond length (189 pm), dissociation energy (2.11 eV), and curvature at the equilibrium bond length (715 cm−1 vibrational band) all agree with experiment for S2 . Table 1: MEAM potential parameters for the sulfur single bond potential. The reference structure is the zig-zag structure described in the text. MEAM

Ec 2.75

alat 2.059

α 6.0

A 1.9

β (0) 5.75

β (1) 7.5

β (2) 7.5

t(1) 0.01

t(2) 1.90

t(3) -2.70

Cmin 0.20

Cmax 2.00

δ+ 0.098

δ− 0.2

7


β (3) 7.5


Page 8 of 19

Table 2: MEAM-BO potential parameters for the sulfur double bond. Bond Order 2

Rbond 1.89

β (0) 0.1

β (1) 5.5

β (2) 8.0

β (3) 0.0

p(0) 0.0

p(1) 1.0

p(2) 1.0

p(3) 1.0

g0,0 0.585

g0,1 44.5

g0,2 -90.0

Results Short-range interactions Since only sulfur is considered for this potential, the only relevant doubly bonded molecule is S2 , which is correct by construction as discussed in Section 3. This is a significant improvement over existing sulfur potentials which cannot fit both single and double bond characters exactly and therefore must approximate one or both of these bond energies. See for example. 3 For the remainder of the validation, we will only consider single covalent bonds and dispersion interactions among sulfur atoms. In principle, this behavior can be reproduced using only the standard MEAM contribution (Eq. 1) combined with the D3 correction. We first consider pure sulfur structures with a variety of high coordination numbers. Table 3 shows the nearest neighbor distances and energy of each of these structures for both the new MEAM potential and DFT results. DFT results were calculated using the potential parametrization of Perdew, Burke and Ernzerhof (PBE) 12 in the Vienna Ab-initio Software Package (VASP) 13,14 except where a different method is explicitly given. As can be seen, all high coordination structures have energy and bond lengths much higher than S8 , and bond lengths increase with increased coordination. The notable exception is the low energy in DFT of the simple cubic sulfur structure. This is due to the ability for sulfur to break its octet in this structure as it does, for example, in sulfur hexafluoride (SF6 ). 15 The MEAM formalism has no way of predicting this result. However, as the simple cubic crystal solid does not occur in nature, (indeed, DFT predicts it to have an energy per atom almost 0.5 eV higher than S8 ) and this ability in sulfur only rarely affect its chemistry in organic molecules,

8


Page 9 of 19 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


this is considered to be an unimportant discrepancy in cases for which this potential in intended. Next, we consider the octasulfur (S8 ) molecule, other low index sulfur rings and linear and polymer chain sulfur structures, all with coordination 2. Sulfur is unique in that it has been readily observed to form a large number of ring structures containing between five and twenty members. 16 S8 is particularly important for validation as it is the most common allotrope of sulfur, and its solid crystal structure has been well studied. We specifically consider: S5 , which has been detected in sulfur vapor; S6 , hexasulfur, which adopts a chair confirmation, similar to cyclohexane; S7 , a less stable structure; S8 , octasulfur which adopts a crown conformation consisting of a pair of concentric squares offset along the axis normal to them; and S12 the most stable cyclo-allotrope of sulfur. Figure 1 shows the structure of each of these rings. Each of these individual molecules can be characterized by relative stability, bond distance and bond angle. The comparison of these structures between DFT and the MEAM potential as well as a comparison of the linear molecule and zigzag and relaxed spiral chains are presented in Table 4. Due to the complex bonding in the linear and zigzag structures and issues related to lone electron pairs discussed above, we do not find good agreement with DFT for those configurations. However, for a spiral structure and for the cyclic rings considered, we see good agreement between the MEAM potential and DFT. In particular, we see the correct ordering of the relative energies. We also consider an alternate “boat” configuration of hexasulfur, to approximate the internal dynamic of these cyclic rings. MEAM predicts the “boat” configuration to be higher in energy than the “chair” configuration (see Figure 1), in agreement with first principles. We also consider the elastic behavior of S8 . Using MEAM, the force constant matrix for an isolated equilibrium S8 molecule is determined. The eigenvalues of these matrices are then compared as a quantitative measure of the low temperature elastic properties of S8 using the MEAM potential. The results for the force constants of the 18 non-trivial vibrational modes (24 degrees of freedom minus 3 translational and 3 rotational) of S8 are presented in

9



Page 10 of 19

Figure 2 and compared to experiment. For the 10 low frequency modes (bending modes) the RMS error between MEAM and experiment is 0.92 eV/˚ A2 , while for the 8 higher frequency stretching modes the RMS error is 11.0 ev/˚ A2 . As the frequency of the stretching modes will be highly dependent on the interactions among third and fourth neighbors and on the variation in the chiral angles, which are not explicitly considered in the reference configuration, this larger deviation in stretching modes is to be expected. We note that the force constants were not used in the fitting of the potential. Table 3: Comparison of per atom energies Ef , D3 energy contributions ED3 , and equilibrium bond lengths re for a variety of high coordination sulfur structures. ∆E/ES8 shows the energy difference as a fractional deviation from the S8 energy. For all listed structures, the formation energy is relative to S8 (-2.78 eV). Structure Graphene

Coordination 3

Diamond Cubic

4

Simple Cubic

6

BCC

8

HCP (Ideal ac )

12

FCC

12

MEAM-BO PBE MEAM-BO PBE MEAM-BO PBE MEAM-BO PBE MEAM-BO PBE MEAM-BO PBE

Ef (eV) 0.58 0.76 0.61 0.82 0.78 0.45 0.75 0.94 0.78 1.11 0.78 1.11

∆E/ES8 0.21 0.27 0.22 0.29 0.28 0.16 0.27 0.34 0.28 0.40 0.28 0.40

ED3 (eV) -0.13 -0.13 -0.3 -0.26 -0.25 -0.24 -0.28 -0.26 -0.28 -0.26 -0.28 -0.26

re (˚ A) 2.28 2.31 2.33 2.49 2.84 2.58 3.07 2.74 3.16 2.82 3.17 2.82

Figure 1: Schematics of the various sulfur rings considered in this paper. (a) S5 (b) “chair” configuration of S6 (c) “boat” configuration of S6 (d) S7 (e) S8 (f) S12

10


Page 11 of 19 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 4: Per atom formation energies and bond lengths for coordination 2 structures including several isolated cyclic sulfur rings using both DFT and MEAMBO results. DFT results were calculated using ωB97M-V/def2-TZPPD density functional. With the exception of S8 , per atom formation energy is given relative to S8 . For S8 the absolute formation energy is given. ∆E/ES8 shows the energy difference as a fractional deviation from the S8 energy. We note that both predict S12 to be the most stable structure followed by S8 . Multiple bond lengths are listed for S7 and the “boat” configuration of S6 due to the broken symmetry of these molecules. Structure zigzag spiral Line S5 S6 (“chair”) S6 (“boat”) S7 S8 S12

MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO ωB97M-V/def2-TZPPD MEAM-BO Experiment MEAM-BO ωB97M-V/def2-TZPPD

Ef 0.013 0.38 0.001 .01 0.90 0.72 0.01 0.136 0.005 0.026 0.024 0.127 0.012 0.025 -2.756 -2.78 -0.005 -0.027

∆E/ES8 0.005 0.137 0.0004 0.004 0.323 0.259 0.004 0.049 0.002 0.009 0.009 0.046 0.004 0.009

-0.002 -0.010

re (˚ A) 2.06 2.06 2.06 2.07 2.35 2.08 2.05 2.04 2.05 2.06 2.05, 2.06 2.03, 2.09 2.05, 2.06 2.00, 2.04 2.04 2.055 17 2.05 2.05

Figure 2: Spring constants for the 18 nontrivial modes of an isolated S8 molecule for both experiment 18 and MEAM-BO calculated from the equilibrium force constant matrix.

11



Long-range interactions (Dispersion Effects) We here explore the long-range interaction between sulfur atoms. This includes the bonding in molecular solids and the interaction among sulfur molecules in a liquid or gas phase. We note that the MEAM results in this section are pure predictions and were not used to fit the sulfur potential. It is due to the importance of the dispersion and the ability of the D3 correction to describe it that these results match closely with experimental results. To demonstrate the necessity of including dispersion effects in the potential, we first calculate the interaction between two S8 molecules. The molecules are translated relative to one another along an axis which passes through the center of their rings. The relative energy of these molecules is shown in Figure 3 with and without the dispersion term included, as well as a DFT calculation using the RI-MP2 method with the aug-cc-pVTZ basis set to account for dispersion forces. While, MEAM-BO + D3 gives a slightly longer minimum energy distance relative to RI-MP2, the binding energies agree within 15 meV. A more refined D3 parametrization, designed specifically for the MEAM potential, might decrease this discrepancy, however, the small differences between methods should not significantly alter the observed behavior. It is evident, however, that inclusion of the D3 term is essential for correctly describing the interaction of sulfur molecules. Furthermore, the dispersion interaction also plays a role in the stability of the molecular solid. The three low energy molecular solid structures of S8 rings, α-, β-, and γ-sulfur are used to test the interaction here. Results for these solids are summarized in Table 5. We see that for the MEAM potential, all three phases are stable, with the favorable formation energy compared to isolated S8 and lattice parameters in agreement with experiment. Simple thermodynamic quantities were also considered. The α and β phases were stable up to 400 K (β-sulfur is reported to melt at approximately 390 K experimentally), and the γ phase was stable up to 200 K. The heat capacities and thermal expansion of the solids is in agreement with experiment.

12


Page 12 of 19

Page 13 of 19 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


Figure 3: Energy of two parallel concentric S8 molecules relative to isolated molecules compared to DFT. The D3 term correctly models the dispersion interaction giving the approximately correct energy and distance for equilibrium compared with the electronic structure calculation. α-S

β-S

Experiment

MEAM-BO

Ef (eV)

-0.13

-0.14

a (˚ A)

10.46a

10.93

10.8125b

11.35

8.442c

8.99

b (˚ A)

12.87a

14.34

10.7232b

11.47

13.025c

13.38

c (˚ A)

24.49a

24.19

10.6883b

11.06

9.356c

7.88

95.746b

99.5

124.98c

113.0

0.75a

0.85

0.82

β (◦ ) J CP ( gK ) ˚ ∆V A3 ∆T ( K )

0.73a

0.84

0.0037d

0.0042

Experiment

γ-S MEAM-BO

Experiment

-0.13

MEAM-BO -0.14

C11 (GPa)

24

e

19.

18.

16.

C22 (GPa)

20.5

e

19.

22.

24.

C33 (GPa)

48.3

e

22.

17.

18.

C12 (GPa)

13.3

e

8.

11.

9.

C13 (GPa)

17.1

e

9.

8.

9.

C23 (GPa)

15.9

e

14.

10.

11.

C44 (GPa)

4.3

e

5.

5.

3.

C55 (GPa)

8.7

e

7.

4.

6.

C66 (GPa)

7.6

e

6.

3.

2.

a Reference 19 b Reference 20

13



High Temperature Behavior As sulfur is heated, sulfur vapor, predominately in the form of S2 dimers is formed. The S3 trimer, with a structure similar to ozone, has also been observed in sulfur vapor, comprising 10% of species at 713 K. 19 S3 has also been observed in high temperature liquid sulfur. In order to confirm that the above developed potential will adequately describe the thermodynamics of sulfur, it is important to confirm the correct behavior at high temperature. In order to test the composition of sulfur opon heating, a simulation was performed wherein a collection of S8 rings were allowed to interact and evolve dynamically with gradually increasing temperature at zero pressure. At a temperature of approximately 400 K, the rings began to open, forming S8 chains. In principle, under high pressure, these chains could link, forming amorphous sulfur consisting of long chains. While such linking was rarely observed in simulation, the frequency with which possible linkages could occur was too low and the timescale for polymerization too long to observed using molecular dynamics methods. Given the low energy of the periodic sulfur spiral used as the reference configuration, these long chains are expected to form for the developed potentials on longer timescales than are computationally feasible. Upon further heating, the S8 begin to break, forming smaller S2 , S3 , and S4 molecules. The S4 molecules then quickly break into pairs of S2 dimers such that, at a temperature of 800K and 2 ns of simulation time, no S4 remains. The remaining sulfur is predominately in the form of S2 dimers while 10% is S3 , in agreement with experiment. No isolated sulfur atoms are observed. Based on these observations, the newly developed potential will be appropriate for high temperature simulations on the behavior of liquid and gaseous sulfur. As S3 is an important factor in the reactivity of liquid sulfur, the potential should also have applications, combined with other elemental potential, studying chemical reactions.

14


Page 14 of 19

Page 15 of 19 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


Conclusions We present an interatomic potential for sulfur based on the new MEAM-BO formalism. 5 To account for dispersion effects which must be included to correctly model the long range interactions important for organic molecules, the DFT-D3 8 energy term was directly added to the MEAM-BO formalism, with the MEAM-BO energy being modified in a systematic way to correct for the short range effects of this inclusion. The optimized sulfur potential correctly models the S2 dimer by construction, due to the bond-order terms added to MEAM to produce MEAM-BO. In addition, using a zigzag reference structure, which should also be appropriate for other covalently bonded elements, particularly oxygen, the properties of a wide variety of singly bonded sulfur structures are correctly predicted by the new potential, including a number of cyclic sulfur rings, sulfur molecular solids, and a number of high coordination structures not observed in nature. In addition, the dispersion correction insures that the long range interaction between sulfur rings is correctly predicted. By combining the MEAM-BO formalism and the independent D3 dispersion correction, we were able to correctly predict the long-range interaction observed between sulfur molecules without explicitly fitting to these values. In addition, the potential serves as a component towards the modeling of more complicated organic molecules using MEAM-BO, supplementing the already developed C-H potential. This potential can be used in conjunction that potential, 5 and future potentials including oxygen, nitrogen, phosphorous, and other elements common in organic chemistry to create a comprehensive, reactive model. In addition, as the MEAM-BO formalism is fully compatible with existing MEAM potentials, commonly used for metals, the ability to reliably model atomic systems containing both organic and metallic components only requires a small number of additional parameters once potentials for the individual elements have been developed. In the absence of bond-order terms, the original MEAM formalism is recovered. Finally, the modification to the MEAM-BO formalism to allow the inclusion of D3 dispersion like effects (detailed in Equations 4-7) allow 15



for a significant expansion of the types of interactions MEAM can account for. Due to the correction to the UEOS the addition of the D3 term is fully consistent with both MEAM and MEAM-BO. While the D3 addition used here was fit for PBE density functionals, a reparametrization for MEAM could be possible as the database of MEAM potentials for covalently bonded molecules increases. This will significantly improve the predictive power of such potentials at intermolecular length scales which preserving the high accuracy of the MEAM-BO formalism.

Acknowledgement The authors would like to thank the Center for Advanced Vehicular Systems (CAVS) at Mississippi State University for supporting this work.

References (1) Nelson, D. L.; Cox, M. M. Lehninger Principles of Biochemistry 3rd ed.; Worth Publishers, New York, 2000. (2) Van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. ReaxFF: a reactive force field for hydrocarbons. Journal of Physical Chemistry A 2001, 105, 9396–9409. (3) J¨arvi, T. T.; van Duin, A. C. T.; Nordlund, K.; Goddard, W. A. Development of Interatomic ReaxFF Potentials for Au–S–C–H Systems. Journal of Physical Chemistry A 2011, 115, 10315–10322. (4) Liang, T.; Phillpot, S. R.; Sinnott, S. B. Parametrization of a reactive many-body potential for Mo–S systems. Physical Review B 2009, 79, 245110. (5) Mun, S.; Bowman, A. L.; Nouranian, S.; Gwaltney, S. R.; Baskes, M. I.; Horstemeyer, M. F. Interatomic Potential for Hydrocarbons on the Basis of the Modified 16


Page 16 of 19

Page 17 of 19 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


Embedded-Atom Method with Bond Order (MEAM-BO). Journal of Physical Chemistry A 2017, 121, 1502–1524. (6) Baskes, M. I. Modified embedded-atom potentials for cubic materials and impurities. Physical Review B 1992, 46, 2727–2742. (7) Lee, B.-J.; Baskes, M. I. Second nearest-neighbor modified embedded atom method potential. Physical Review B 2000, 62, 8564–8567. (8) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. Journal of Chemical Physics 2010, 132, 154104. (9) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J Comput Chem 2011, 32, 1456–1465. (10) Rose, J. H.; Smith, J. R.; Guinea, F.; Ferrante, J. Universal features of the equation of state of metals. Physical Review B 1984, 29, 2963–2969. (11) Lee, T.; Baskes, M. I.; Valone, S. M.; Doll, J. D. Atomistic modeling of thermodynamic equilibrium and polymorphism of iron. J. Phys.: Condens. Matter 2012, 24, 225404. (12) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Physical Review Letters 1996, 77, 3865–3868. (13) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Physical Review B 1993, 47, 558–561. (14) Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical review B 1996, 54, 11169–11186. (15) Reed, A. E.; Weinhold, F. On the role of d orbitals in sulfur hexaflouride. Journal of the American Chemical Society 1986, 108, 3586–3593. 17



(16) Rauchfuss, T. B. Elemental Sulfur and Sulfur-Rich Compounds I. Topics in Current Chemistry; Springer-Verlag, 2005. (17) Rettig, S. J.; Trotter, J. Refinement of the structure of orthorhombic sulfur, alpha-S8. Acta Crystallographica Section C: Crystal Structure Communications 1987, 43, 2260– 2262. (18) Jackson, A. J.; Tiana, D.; Walsh, A. A Universal Chemical Potential for Sulfur Vapours. Chemical Science 2016, 7, 1082–1092. (19) Meyer, B. Elemental Sulfur. Chemical Reviews 1976, 76, 367–388. (20) David, W. I. F.; Ibberson, R. M.; Cox, S. F. J.; Wood, P. T. Order–disorder transition in monoclinic sulfur: a precise structural study by high-resolution neutron powder diffraction. Acta Crystallographica Section B: Structural Science 2006, 62, 953–959. (21) Gallacher, A. C.; Pinkerton, A. A. A Redetermination of Monoclinic γ-Sulfur. Acta Cryst. 1993, 49, 125–126. (22) George, J.; Deringer, V. L.; Wang, A.; M¨ uller, P.; Englert, U.; Dronskowski, R. Lattice thermal expansion and anisotropic displacements in alpha-sulfur from diffraction experiments and first-principles theory. The Journal of chemical physics 2016, 145, 234512. (23) Hearmon, R. F. S. The elastic constants of anisotropic materials—II. Advances in Physics 1956, 5, 323–382.

18


Page 18 of 19

Page 19 of 19 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


Graphical TOC Entry

19


A Dispersion-Corrected Modified Embedded-Atom Method Bond

Recommend Documents