A third-generation charge optimized many body (COMB3) potential for nitrogen-containing organic molecules

A third-generation charge optimized many body (COMB3) potential for nitrogen-containing organic molecules

Computational Materials Science 139 (2017) 153–161 Contents lists available at ScienceDirect Computational Materials Science journal homepage: www.e...

966KB Sizes 3 Downloads 20 Views

Computational Materials Science 139 (2017) 153–161

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A third-generation charge optimized many body (COMB3) potential for nitrogen-containing organic molecules Jackelyn Martinez a, Tao Liang (梁涛) b, Susan B. Sinnott b, Simon R. Phillpot a,⇑ a b

Department of Materials Science and Engineering, University of Florida, Gainesville 32611, USA Department of Materials Science and Engineering, The Pennsylvania State University, University Park, PA 16802, USA

a r t i c l e

i n f o

Article history: Received 9 May 2017 Received in revised form 10 July 2017 Accepted 18 July 2017

Keywords: Molecular dynamics Reactive potential Organic Nitrogen COMB3

a b s t r a c t This work presents a new empirical, variable charge potential for NCOH systems in the third-generation charge-optimized many-body (COMB3) potential framework. The potential parameters were determined by fitting each binary system separately to the experimental enthalpy of formation and geometries for experimentally important molecules. The resulting parameter set reproduces the energetics and geometries of many nitrogen-based organic molecules with a similar average error as other reactive potentials for the same systems. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Complex polyatomic NCOH molecules have a wide range of applications, including polymer synthesis [1,2] and organic electronics [2]. Large nitrogen containing molecules are of particular interest for explosives [3,4] such as the recently synthesized pentazole [5], a five-nitrogen ring. Such molecular systems are notoriously difficult to synthesize. A variety of computational methods have been exploited to study NCOH systems. Quantummechanical and electronic-structure methods offer a high degree of materials fidelity for small molecular systems. However, in order to calculate the properties of interest for many applications, the dynamical evolution at long time and length scales is required, for which the computational load is often too high for electronicstructure methods [6]. Atomic-level simulation methods such as molecular dynamics (MD) using classical interatomic potentials offer an alternative to electronic-structure calculations by accounting for the effects of the interatomic interactions associated with the valence electrons without an explicit description of the electrons themselves. This allows for increased computational efficiency and thus allows for simulations with longer length and time scales. Reactive potentials, unlike conventional force fields, have the capacity to describe dynamical processes such as bond breaking ⇑ Corresponding author. E-mail address: [email protected] (S.R. Phillpot). http://dx.doi.org/10.1016/j.commatsci.2017.07.019 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.

and formation, essential to any realistic description of organic materials. They also have the capacity to describe multiple bonding types and are well-suited to simulate the complex applications of complex polyatomic NCOH molecules. REBO [7,8], MEAM [9], ReaxFF [10] and COMB3 (third-generation charge-optimized many body) are all reactive potentials that can describe complex bonding environments. While potentials for NCOH systems have already been developed within the REBO, MEAM, and ReaxFF frameworks, each additional potential set developed in the COMB3 framework is transferable and compatible with other parameter sets, allowing for the study of a wide variety of disparate systems including but not limited to, Ti [11], Al [12], Zr [13], and Cu [14] as well as the corresponding oxide and nitride systems. REBO potentials have been developed to specifically focus on carbon and carbon–oxygen systems with little focus on metal or metal oxide systems. The MEAM potential was developed to model alloys well; while parameterizations exist for NCOH systems, they do not model diamond and graphite systems or unsaturated carbon molecules as these were not specifically included in the training set [15]. Also, the parameter sets for the MEAM formalism are often developed for specific targeted applications, with the result that the various parameter sets are often not consistent and cannot be used with each other. The ReaxFF formalism does include charge transfer as does COMB3. However, like MEAM, ReaxFF parameter sets are developed for specific applications and are often not transferable.

154

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

The rest of the paper is organized as follows: Section 2 outlines the general form for COMB3. The fitting procedure used and the targeted properties for NCOH are given in Section 3. This builds on the previously developed COH parameterization of COMB3. Section 4 summaries the predicted properties of the potential. Our conclusions are in Section 5.

spline function specifically to model the complex bonding behavior of carbon [14]. The long-range interactions take the form of a classic LennardJones formalism where eij is the energy scale and rij is the length:

U v dW

2 XX v dW 4 ¼ 4eij i

2. COMB3 potential The third-generation COMB3 formalism builds upon the second generation COMB potential (COMB2) developed by Liang et al. [14]. The third generation includes minor changes to the functional form that increase the flexibility of the formalism to describe many differing bonding environments, and a major extension to the coordination function to include a spline function that is more capable of describing the multiple bonding environments present in carbon molecules and compounds. A complete description of the COMB3 formalism can be found elsewhere [14]; only a brief description is provided below, together with details for the targeting and fitting strategies used in this paper. The COMB3 formalism is a bond-order type potential with variable charge. The general form for COMB3 is:

U Tot ¼ U Ele ðq; rÞ þ U Short ðq; rÞ þ U v dW ðrÞ þ U Corr ði; j; kÞ

ð1Þ

where the total potential energy (UTot) is given as the sum of the electrostatic energy (UEle), the short-range energy (UShort), the long-range or van der Waals energy (UvdW) and a correction term (UCorr). The electrostatic energy is:

U Ele ðq; rÞ ¼ U q ðq; rÞ þ U qq ðq; rÞ þ U qZ ðq; rÞ þ U Polar ðq; rÞ

ð2Þ

The electrostatic energy is the sum of the energy to form a charge from a neutral atom (Uq), the charge-charge interactions (Uqq), the charge-nuclear interactions (UqZ) and the energy associated with the polarizability (UPolar). The energy to form a charge (Uq) takes into account the electron affinity and ionization energies as well as field affects such as the change in electronegativity and chemical hardness as a function of the bonding environment [16]. The charge-charge interaction is an integration over the Streitz and Mintmire style charge density of two atoms [17], while the chargenuclear interaction is the sum of the interactions between the point charge of one atom and the shielded nuclear interaction of a neighboring atom. These two terms are handled via the Wolf charge-neutralized direct-summation method [18]. The energy associated with the polarizability represents the relative tendency of a charge distribution to distort in the presence of an external field and uses the fluctuation dipole model of Stern et al. [19]. The short range interactions are captured by an extended Tersoff-like potential, in which the repulsive and attractive interactions are modified by a charge dependent correction term:

U Short ðq; rÞ ¼

XX i

  F c ðr ij Þ V R ðr ij ; qi ; qj Þ  Bij V A ðrij ; qi ; qj Þ

ð3Þ

j>i

In Eq. (3), Fc is a Tersoff cutoff function and smoothly terminates the function between two defined cutoff radii, VR is the repulsive term, VA is the attractive term, and Bij is the bond order function. The bond order function:

( Bij ¼

" #gi )1=2gi N X 1þ ½F c ðr ik Þ  nðr ij ; r ik Þ  g ij ðcosðhijk ÞÞ þ Pij

ð4Þ

k–ij

was the subject of the largest change from COMB2 to COMB3. Here, n is an asymmetric fitting parameter that weakens the longer bond length if rij is not equal to rik; gij is the angular term and takes the form of a sixth-order polynomial; Pij is the coordination function. The coordination function was extended in COMB3 to include a

j>i

rijv dW r ij

!12 

rvij dW r ij

!6 3 5

ð5Þ

The final contribution to the energy shown is the correction term, which allows for the distinction between similar phases such as hexagonal closed packed and face center cubic phases as well as increased accuracy in ternary potentials.

U Corr ðcosðhijk ÞÞ ¼

6 X ang n bijk cos ðhijk Þn

ð6Þ

n¼1

3. Parameterization of COMB3 NCOH The COMB3 parameters for COH systems were previously developed by Liang et al. [14] and were used for the development of the NCOH parameter sets described here. The fitting procedure follows the basic procedure outlined by Martinez et al. [20], using the POSMat fitting software [21]. Elemental parameters for all species were previously developed for the COMB3 formalism by Liang et al. [14] and Cheng et al. [11]. These were used here as the basis of the development of the bond specific or binary parameters. For the development of all of the binary parameters, training sets were constructed for molecules of interest, which included experimentally stable molecules such as NH3 for NAH parameterization and molecules that exhibit a range of bonding environments, such as CH2NH and CH3NH2 (double versus single bond CAN bond) for the CAN parameterization. A complete list of targeted molecules is given in subsequent sections. Training sets for the parameterization included the molecular geometry, such as bond angles and lengths, and the equilibrium and non-equilibrium energetics. The binary parameter set was then optimized using the Simplex algorithm [22] and simulated annealing [23] where the error in the target properties is given by:

FðaÞ ¼

2 X  W j Pj ðaÞ  Poj

ð7Þ

j

where a denotes the current parameter set for a given optimization iteration, Poj is the target value for a property j, Pj(a) is the calculated property based on the current parameter set a, and WJ is the inputted weight for that property. Target energies were calculated by converting the experimental enthalpy of formation into an energy per atom (eV/atom) through the use of the general chemical reaction: DH f a b d N2 þ O2 þ cC graphite þ H2 ! Na Ob C c Hd 2 2 2

ð8Þ

where the energies used on the left side of the equation are COMB3 values from our previous parameterizations for these systems, DHf is an experimental value, and the energy of the molecule on the right is the target for the fitting (cohesive energy). Throughout the remainder of the work, the cohesive energies (NaObCCHd) are reported in eV/atom, while the enthalpies of formation (DHf) are reported in kJ/mol. Both are the typical units used for reporting these properties; for reference: 1 eV/atom = 96.5 kJ/mol. In order to target the bond strength and the bond bending characteristics, curves of energy versus bond length and/or bond angle were calculated for key molecules using B3LYP/6-311G⁄⁄. While such quantum chemical calculations can offer a high degree of accuracy, there generally still remains a discrepancy between the B3LYP/6-311G⁄⁄ values and the experimental values. In order to

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

155

Fig. 1. Target energy curve calculation for NH3 illustrating the shifting of target quantum chemical values to experiment.

maintain consistency within the training database, the calculated values from the quantum chemical calculations were adjusted to better match the experimental values. Specifically, the B3LYP curve was shifted up or down by the difference in energy from the electronic-structure calculations and experimentally derived cohesive energy value for the equilibrium structure from Eq. (8), as illustrated in Fig. 1 for ammonia. Experimental values are more suitable than the results of quantum chemical calculations as target values for obvious reasons. However, experimental values do not exist for all targeted properties, such as bond strength (first derivative of energy). Such quantities must, therefore, be calculated from quantum chemical methods which, of necessity, introduces inconsistency into the fitting database when combined with experimental values. Shifting the entire energy curve is one way to account, at least in part, for this inconsistency. For the molecules used in this study for which there are experimental values for the energetics available, this shift

was between 0.0084 eV/atom and 1.21 eV/atom, with an average of 0.138 eV/atom (13.3 kJ/mol). In the context of cohesive energies, a shift of 0.2 eV is not considered large when the expected precision in COMB3 for the property may be within 1%, as given in Martinez et al. [20]. However, within the context of enthalpies of the equilibrium structure, such a shift can be very impactful. In the case of NH3, using the B3LYP cohesive energy to calculate the enthalpy of formation in tandem with the COMB3 values for N2, O2, H2, and graphite, as shown in Eq. (8), results in an enthalpy of formation of 107.4 kJ/mol as compared to the experimental value of 34.09 kJ/mol. From this we can conclude that the enthalpy of formation is a more robust measure of accuracy than the B3LYP cohesive energy. There are experimental values for the enthalpies of formation of all but two (NH2ONH2 and HNCNH) of the molecules in the training database. In these cases, the cohesive energies calculated from B3LYP were targeted. This may have introduced an inconstancy

Fig. 2. Target energy curve calculation for NH3 illustrating the targeting of bond strength (curvature) where the COMB3 values have been shifted to match quantum chemical bond lengths in order to illustrate the error in curvature.

156

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

Table 1 COMB3 binary parameters for NCOH. Parameters

NAH

HAN

NAO

OAN

NAC

CAN

A B1 B2 B3 K

280.617 3.44934 7.61271 5.32756 5.7607 4.13002 0.541527 5.86348 1.39802 0.0001135 0.000850371 7.933 2.1359 0.283076 0.249446 0.0950972 0.000246534 0.000524908 0.00056752 0.00014139 1.6 1.8 0.585232 0.316466

280.617 3.44934 7.61271 5.32756 5.7607 4.13002 0.541527 5.86348 7.41726 0.697904 0.251519 1.47748 0.229852 0.954102 0.382336 6.75946 0 0 0 0 1.6 1.8 0.001 0.001

13110.9 327.321 17.6785 37.065 6.89799 3.73064 2.59741 1.91364 0.268183 0.0549322 0.0120463 3.38298 1.85513 1.19956 0.301961 0.761188 6.9895E05 7.9784E05 5.1698E05 0.00010862 2.1 2.4 0.234888 0.169673

13110.9 327.321 17.6785 37.065 6.89799 3.73064 2.59741 1.91364 0.0439581 1.80325 1.80508 0.32263 0.0499692 2.00816 0.521228 3.3922 4.3052E06 0.00207904 7.2895E06 0.00015165 2.1 2.4 0.00010149 0.0010149

9142.53 51.5103 0.0123541 0.0053851 6.75861 1.23296 0.159932 8.20817 0.29099 2.1165 1.8033 2.5096 0.1825 2.2362 2.4748 0.985207 2.53728 0.450555 1.04495 0.113829 1.8 2 4.44248 2.69147

9142.53 51.5103 0.0123541 0.0053851 6.75861 1.23296 0.159932 8.20817 0 0.0209992 0.141332 0.197857 0.0086649 0.41238 0.299886 0.007702 0.280786 0.0381319 1.97195 0.111791 1.8 2 9.80336 4.90356

a1 a2 a3 B b6 b5 b4 b3 b2 b1 b0 c0 c1 c2 c3 Rs(Å) Ss(Å) PX PJ

into the training database. However, because these molecules were primarily included to target energy vs bond angles, the energies of the equilibrium structures were not targeted very heavily, while the variations of the energy with bond angle were. The bond strength (the first derivative of the energy) is targeted by specifically determining the dependence of the energy on the bond length. In the case of NH3, the error in the bond length for the equilibrium structure is 5.4% which is within the expected accuracy for COMB3 [20]. When this error is taken into account as shown in the right side of Fig. 2, it can be seen that the error in curvature near the equilibrium structure (±0.1 Å) is less than 5%, which is also within the expected accuracy of COMB3 [20], found in potentials prior COMB3 parameterizations. Such levels of accuracy are consistent with other methods [24].

The energy and curvature (bond strength) close to the equilibrium bond length are reproduced within the expected accuracy previously mentioned. The equilibrium structures were allowed to structurally relax during the fitting process and the bond lengths and angles were targeted directly. Charges are not directly targetable; however, reasonable bounds were set on the charges for each element in the fitting and the charges were monitored. Nitrogen-hydrogen parameters were developed first, followed by nitrogen–oxygen parameters and finally nitrogen-carbon parameters. Once the short range interactions were fit and all predicted properties for the targeted molecules were optimized from the COMB3 formalism, the parameterization was tested on a much larger group of molecules and the correction term was

Table 2 COMB3 binary parameters for COH. Parameters

CAO

OAC

CAH

HAC

OAH

HAO

A B1 B2 B3 K

1534.88 889.569 14.0948 0.0 3.41466 2.77391 1.39992 0.0 0.0 0.0 0.0 0.032709 0.104852 0.06711 0.115834 0.446323 0.0 0.0 0.0 0.0 1.6 1.9 2.10127 1.09580

1534.88 889.569 14.0948 0.0 3.41466 2.77391 1.39992 0.0 0.0 0.0 0.0 0.00278771 0.00037601 0.0007348 0.0105115 0.170033 0.0 0.0 0.0 0.0 1.6 1.9 0.157626 0.186831

273.134 5.97783 0.0 0.0 6.26017 0.418769 0.0 0.0 0.0 0.155586 0.0110468 0.207023 0.0206675 0.173997 0.44355 0.0462052 0.0 0.0 0.0 0.0 1.5 1.8 1.23224 0.939023

273.134 5.97783 0.0 0.0 6.26017 0.418769 0.0 0.0 3.82336 0.697905 0.251519 1.47748 0.229851 0.954102 0.382336 6.75945 0.0 0.0 0.0 0.0 1.5 1.8 0.001 0.001

245.338 98.4247 11.6990 0.0 3.36932 2.22378 2.06777 0.0 4.42858 0.505441 1.82534 1.87431 1.91884 2.75901 0.757369 0.334137 0.502193 0.0169267 2.36744 0.00006836 1.5 1.8 1.96828 0.98415

245.338 98.4247 11.6990 0.0 3.36932 2.22378 2.06777 0.0 1.58139 0.697905 0.251519 0.147748 0.229851 0.954192 0.382336 4.06954 0.0 0.0 0.0 0.0 1.5 1.8 0.001 0.001

a1 a2 a3 B b6 b5 b4 b3 b2 b1 b0 c0 c1 c2 c3 Rs(Å) Ss(Å) PX PJ

157

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161 Table 3 COMB3 correction parameters for NCOH. BOND

P6P0

P6P1

P6P2

P6P3

P6P4

P6P5

P6P6

NANAC CANAC CANAH OANAC OANAH HANAC NACAN NACAO NACAH CACAN NAOAC HAOAN CAHAC* HAHAC*

1.61566 2.34738 0.867325 1.34847 0.916703 0.228459 1.41381 3.14066 2.03915 0.533823 1.03263 0.330134 0.636633 0.217

0.0 0.0 0.0997167 0.0 0.0 0.0 0.0 0.0 0.103805 0.0 0.0 0.0 1.27327 0.435

0.0 0.0 0.154539 0.0 0.0 0.0 0.0 0.0 0.0618152 0.0 0.0 0.0 0.636633 0.0

0.0 0.0 0.0636347 0.0 0.0 0.0 0.0 0.0 0.0017006 0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

Table 4 Summary of targeted properties and their experimental, B3LYP and COMB3 values for NAH molecules in the fitting database. NAH Properties Property

Molecule NH3

Enthalpy of Formation (kJ/mol) (0 K) Binding Energy (eV/atom) Bond Length (Å) Bond Angle (°) Charge on N (eV)

NH2

Exp [27]

COMB3

38.95 3.096 1.012 106.7

35.34 3.011 1.070 105.7 0.410

B3LYP [27] 3.204 1.019 106.5 0.568

Fig. 3. Charge on Nitrogen atom as a function of bond length for NH.

used to more accurately reproduce energies for systems containing more than two elements. The COMB3 NCOH parameters are given in Table 1. The COH parameters are given in Table 2; the CH parameters were previously published by Liang et al. [14]. The correction parameters is given in Table 3 [11,25,26]. The correction parameters denoted with an asterisk were developed outside the scope of this work. All of these parameters are now available in LAMMPS.

NH

Exp [27]

COMB3

189.10 2.553 1.020 103.4

140.84 2.653 1.078 105.7 0.282

B3LYP [27] 2.651 1.031 102.1 0.356

Exp [27]

COMB3

B3LYP [27]

358.43 1.702 1.036

310.49 1.970 1.090

1.419 1.065

0.142

0.176

are those with less than three bonds. Due to this limited range of environments, only three molecules were targeted during this step of the parameterization; NH, NH2, and NH3. Ammonia (NH3) was deemed the most important molecule as it is experimentally stable and was therefore targeted more heavily than the other two molecules. The results of the targeted energy and geometries can be seen in Table 4 and can be compared to the B3LYP values in the same table. While COMB3 does not perfectly replicate the experimental values for the energetics, it is not qualitatively worse than the quantum mechanical values. Bond-length versus energy curves were targeted for all three molecules as well as the bond angle versus energy curve for NH2. The COMB3 reproduces the enthalpies of formation for all cases well and does a particularly good job in the case of NH3. COMB3 slightly overestimates the bond lengths, but reasonably reproduces the bond length trend for all three molecules. NH3 was targeted the most heavily, and as the experimental bond angles for NH2 and NH3 are close, the COMB3 bond angle formulation could not distinguish between the two. This resulted in a compromise in the bond angle that slightly favored NH3. In tests in which NH2 was not included in the targeting database, the NH3 angle was reproduced with much greater accuracy. As mentioned previously, the charges were not directly targeted during parameterization, but were monitored; it was found that the charge on nitrogen was quite small and became less negative as the nitrogen–hydrogen bond length increased for all structures, which is consistent with electronic-structure calculations, as shown in Fig. 3 for NH.

3.1. Hydrides 3.2. Oxides NAH bonding exhibits a limited range of bonding environments. Nitrogen generally bonds with sp3 hybridization. The only NAH molecules with radically different bonding environments

For the parameterization of pair-wise nitrogen–oxygen interactions, bond scans were targeted for nitric oxide (NO) and nitrogen

158

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

Table 5 Summary of targeted properties and their experimental, B3LYP and COMB3 values for NAO containing molecules in the fitting database. NAO Properties Property

Molecule NO2

Enthalpy of Formation (kJ/mol) (0 K) Binding Energy (eV/atom) Bond Length (Å) Bond Angle (°) Charge on N (eV)

NO

Exp [27]

COMB3

33.97 3.228 1.193 134.1

42.20 3.223 1.188 134.5 0.290

B3LYP 3.296 1.194 134.3 0.333

Fig. 4. Energy versus bond length for NO and NO2.

dioxide (NO2) as representative cases. While N tends towards sp3 hybridization, the stabilities of NAO systems are limited by the sp hybridization of oxygen. The bond-order terms for OANAO were developed by targeting the bond angle versus energy for nitrogen dioxide. The bond-order terms for NAOAN interactions presented a greater challenge. While Nitrous Oxide (N2O) is an experimentally viable molecule, it has an NANAO bond angle instead of the NAOAN bond angle desired for the development of the bond order terms, for which symmetric molecules were required. It was decided to use the previously developed nitrogen–hydrogen parameters to target the bond angle versus energy of NH2AOANH2. The results for the targeted energy and properties can be seen in Table 5. After a parameter set was obtained, the equilibrium energies for subsequent molecules were tested to ascertain the viability of the parameter set. It was found that the target binding energy for nitric

NH2ONH2

Exp

COMB3

B3LYP

91.04 3.273 1.154

243.42 2.524 1.190

3.303 1.148

0.115

0.159

Exp

COMB3

B3LYP

3.131 1.363 114.8 0.466

2.957 1.437 106.4 0.524

oxide (NO) skewed the accuracy of more experimentally viable molecules. Therefore, the energy for nitric oxide was weighted less than the energy of other molecules; however, its energy versus bond length curvature was still weighed heavily. COMB3 was found to reproduce the binding energy extremely well for NO2, and has a greater accuracy than the quantum chemical methods. The enthalpy of formation was also reproduced reasonably well and the geometry is also adequately described. The bond strength was reproduced extremely well as seen in Fig. 4. As in the previous case, the charges were not directly targeted but were monitored during the fitting process. It was found that in NO2 and NO the nitrogen exhibits a positive charge. In the case of NH2ONH2 however, the nitrogen exhibited a negative charge due to the presence of hydrogen. In all cases the charge on the nitrogen decreased as the nitrogen–oxygen bond length increased. This is consistent with electronic-structure calculations for NO2 and NH2ONH2. COMB3 predicted that the nitrogen charge decreases with bond length for nitric oxide as well; however, B3LYP predicted that it increases. This discrepancy is attributable to the presence of a nitrogen radical, which COMB3 does not replicate. 3.3. Cyanides The carbon–nitrogen parameterization offered the greatest challenge because single, double and triple bonds are all possible. Each of the three bond types has a distinctive bond length and bond strength (first derivative of energy as a function of bond length). Three molecules were chosen, each to target one of these bond types; HC„N, [email protected], and H3CANH2. Only the bond scan for the double bond was targeted heavily because it was found that targeting the energy, bond length and bond strength of three different bond types did not result in viable parameter fits. The double bond was chosen in order to maintain a balance between the energetics of the single bond and triple bond. However, this

Table 6 Summary of targeted properties and their experimental, B3LYP and COMB3 values for CAN containing molecules in the fitting database. CAN Properties Property

Molecule HCN

CH2NH

Exp

COMB3

Enthalpy of Formation (kJ/mol) Binding Energy (eV/atom) Bond Length (Å) Charge on N (eV)

132.38 4.208 1.168

162.76 4.290 1.195 0.216

Enthalpy of Formation (kJ/mol) Binding Energy (eV/atom) Bond Angle (o) Charge on N (eV)

18.6 3.635 112.2

B3LYP 3.920 1.149 0.252

CH3NHCH3

CH3NH2

Exp

COMB3

91.2 3.682 1.273

82.42 3.644 1.259 0.0153

B3LYP 3.794 1.266 0.092

184.9 3.939 180 0.453

4.196 133 0.422

HNCNH 2.87 3.533 111 0.560

3.613 113.2 0.385

Exp

COMB3

B3LYP

7.8 3.450 1.471

68.89 3.472 1.339 0.282

3.589 1.465 0.199

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

159

Fig. 5. COMB3 (abscissa) vs experimental (ordinate) values for enthalpy of formation for NACHO molecules. The diagonal line corresponds to a perfect fit.

4. Results and discussion

Table 7 RMS error in the enthalpy of formation for various NCOH molecules. System

CAH NAH NAO NAOH NAC NACH NACHO All N

RMS Error (kcal/mol) COMB3

ReaxFF [3]

5.7 [14] 8.2 25.9 11.7 15.9 40.4 36.7 34.8

5.1 [14] 12.9 20.3 9.6 15.0 39.2 44.1 36.9

resulted in a much shorter predicted bond length in the case of the single bond than was found experimentally. While COMB3 can accurately distinguish between the double and single bond, the shorter bond lengths were found to cause higher strains in molecules containing rings. This will be discussed in greater detail in a subsequent section. The equilibrium energies and bond lengths for all three molecules are shown in Table 6. The angular function was fit to bond angle versus energy curves for CH3NHCH3 and HNCNH. As in previous cases, the enthalpy of formation was targeted except in the case of HNCNH where it was unavailable. COMB3 accurately reproduced the trend in enthalpies of formation for HCN, CH2NH, and CH3NH2. The enthalpy of formation for the double bond was produced the most accurately. COMB3 also reproduced the binding energy with greater accuracy than the quantum chemical methods for all three cases. The bond angle for CH3NHCH3 was also accurately reproduced. However, it was found that targeting the 180° bond angle greatly destabilized all of the ring molecules in the testing suite. Therefore, the weight was greatly reduced on this molecule and the minimum of the angular function was shifted in order to accommodate the ring molecules of the testing suite. The charge trend seen in the triple to single CN bond was correctly reproduced as well as the trend for CANAC bonds versus NACAN bonds.

Having parameterized the potential parameters with 12 fitting molecules, more than eighty molecules were used to test the entire NCOH potential. The COMB3 vs. experimental values for enthalpy of formation values for the fitting molecules and nearly all of the test molecules are shown in Fig. 5. A small subset of the test molecules were found to be unstable in COMB3 and are not included in Fig. 5. These special cases will be addressed in greater detail below. The RMS error was calculated for all of the stable molecules; Table 7 shows that these errors are comparable to the corresponding RMS errors for the RDX/High Energy ReaxFF potential which includes NACAOAH and is available in LAMMPS [3,28–30]. The most significant differences are for NAH and NACHO, for which COMB3 shows a significantly lower RMS error; that is, the COMB3 provides a better description than the corresponding ReaxFF potential. COMB3 does a particularly good job describing the larger linear molecules in the testing suite. A small subset of molecules from the testing suite are shown in Table 8. The full list of the results of the testing suite can be found in the supplementary materials. These results show promise in describing polymeric materials such as Nylon 6,6. However, the COMB3 potential did not describe all of the molecules well. These outlier molecules had one of two physical characteristics: an internal nitrogen-nitrogen bond which COMB3 predicted dissociated, or a 5-membered ring with smaller than equilibrium bond angles, which COMB3 predicted higher energy than experimental. The equilibrium bond angle for NAC interactions is 111° while there are at least one or two bond angles in these structures that exhibit a bond angle less than 90°. Examples of both types of structures are shown in Table 9. 5. Conclusions This paper has presented a potential for NCOH systems which was developed in order to fill an existing gap within the COMB3

160

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

Table 8 The enthalpy of formation and binding energies for a selection of large complex molecules from the testing suite which COMB3 does a good job of reproducing. Energetics and Geometries of well reproduced molecules Molecule Name/Formula

Molecule Geometry

Enthalpy of Formation (kJ/mol)

Binding Energy (eV/atom)

COMB3

Exp

COMB3

Exp

C4H2N2 Fumaronitrile

340.00

365.69

4.991

5.062

C5H7N (E)-3-Pentenenitrile

125.56

133.34

4.309

4.383

CH3CHNOH Acetaldoxime

22.55

17.88

3.729

3.807

H2NCH2COOH Glycine

390.50

406.26

4.015

4.066

C4H5NO 3-Methylisoxazole

35.65

59.01

4.319

4.400

NH2COOC2H5 Urethane

446.30

392.20

3.989

4.102

CH3NHCH2COOH Sarcosine

367.24

325.35

3.939

4.039

NH2CH2CH2COOH beta-Alanine

424.00

415.53

4.008

4.084

Table 9 Properties of molecules that are not well described by the COMB3 NCOH Potential. Energy and Structures of Outliers Structure Type

Inner NAN Bond Dissociation

5-membered Ring

Molecule

Binding Energy (eV/atom)

Structure

Exp

B3LYP

COMB3

CH2N4

4.045

4.113

3.650

N2H4

3.042

3.137

2.930

C2H6N2O2

3.599

3.664

3.693

C3H4N2

4.389

4.454

3.862

C2H3N2

4.317

4.382

3.433

C4H7NO

4.219

4.274

4.028

framework. It was shown that the potential can accurately reproduce the energetic and geometries of a wide variety of nitrogen containing organic molecules with an accuracy similar to that of other reactive potentials. However, unlike other reactive potentials, the NCOH potential developed here is compatible with other potentials available within the COMB3 framework offering a wide variety of possible applications including but, not limited to, oxidation mechanics of intermetallics, material behavior of oxides and

Predicted

COMB3

intermetallics in water or organic containing environments, organics, polymer containing systems, and high energy materials such as the recently synthesized pentazole [5]. Acknowledgements The authors acknowledge the support of the National Science Foundation (DMR-1005779).

J. Martinez et al. / Computational Materials Science 139 (2017) 153–161

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.commatsci.2017. 07.019. References [1] V.F. Chemie, Technischen Univeritat Darmstadt, 2012. [2] H. Sitter, C. Draxl, M. Ramsey, Small Organic Molecules on Surfaces: Fundamentals and Applications, Springer, Berlin Heidelberg, 2013. [3] A. Strachan, A.C.T. van Duin, D. Chakraborty, S. Dasgupta, W.A. Goddard, Phys. Rev. Lett. 91 (2003). [4] Z.A. Dreger, Y.A. Gruzdkov, Y.M. Gupta, J. Phys. Chem. B 106 (2002) 247. [5] C. Zhang, C. Sun, B. Hu, C. Yu, L. Ming, Science 355 (2017). [6] D.S. Sholl, J.A. Steckel, Density Functional Theory: A Practical Introduction, John Wiley & Sons, Hoboken, NJ, 2009. [7] J. Tersoff, Phys. Rev. Lett. 56 (1986) 632. [8] D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, S.B. Sinnott, J. Phys.: Condens. Matter 14 (2002) 783. [9] M.I. Baskes, Phys Rev B 46 (1991). [10] A.C.T. Van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, J Phys Chem A 105 (2001) 9396. [11] Y.-T. Cheng, T. Liang, J.A. Martinez, S.R. Phillpot, S.B. Sinnott, J. Phys.: Condens. Matter 26 (2014). [12] K. Choudhary, T. Liang, A. Chernatynskiy, S.R. Phillpot, S.B. Sinnott, J. Phys.: Condens. Matter 27 (2015). [13] M.J. Noordhoek, T. Liang, Z. Lu, T.-R. Shan, S.B. Sinnott, S.R. Phillpot, J. Nucl. Mater. 441 (2013) 274.

161

[14] T. Liang, B. Devine, S.R. Phillpot, S.B. Sinnott, J. Phys. Chem. A 116 (2012) 7976. [15] S. Nouranian, M.A. Tschopp, S.R. Gwaltney, M.I. Baskes, M.F. Horstemeyer, Phys. Chem. Chem. Phys. 16 (2014). [16] H. Toufar, K. Nulens, G.O.A. Janssens, W.J. Mortier, R.A. Schoonheydt, F. DeProft, P. Geerlings, J. Phys. Chem. 100 (1996). [17] F.H. Streitz, J.W. Mintmire, Phys. Rev. B 50 (1994) 11996. [18] D. Wolf, P. Keblinski, S.R. Phillpot, J. Eggebrecht, J. Chem. Phys. 110 (1999) 8254. [19] H.A. Stern, G.A. Kaminski, J.L. Banks, R.H. Zhou, B.J. Berne, R.A. Friesner, J. Phys. Chem. B 103 (1999) 4730. [20] J.A. Martinez, D.E. Yilmaz, T. Liang, S.B. Sinnott, S.R. Phillpot, Curr. Opin. Solid State Mater. Sci. 17 (2013) 263. [21] J.A. Martinez, A. Chernatynskiy, D.E. Yilmaz, T. Liang, S.B. Sinnott, S.R. Phillpot, Comput. Phys. Commun. 203 (2016) 201. [22] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambride University Press, New York, 1986. [23] J.A. Nelder, R. Mead, Comput. J. 7 (1965) 308. [24] F. Ercolessi, J.B. Adams, Europhys. Lett. 26 (1994). [25] Y.-T. Cheng, T.-R. Shan, T. Liang, R.K. Behera, S.R. Phillpot, S.B. Sinnott, J. Phys. Condens. Matter 26 (2014). [26] T. Liang, M. Ashton, K. Choudhary, D. Zhang, A.F. Fonseca, B.C. Revard, R.G. Hennig, S.R. Phillpot, S.B. Sinnot, J. Phys. Chem. C 120 (2016) 12530. [27] R. D. I. Johnson, NIST Computational Chemistry Comparison and Benchmark Database, Vol. Release 16a. [28] A. Strachan, E. Kober, A.C.T. van Duin, J. Oxgaard, W.A. Goddard, J. Chem. Phys. 122 (2005). [29] L. Zhang, A.C.T. van Duin, S. Zybin, W.A. Goddard, J. Phys. Chem. 113 (2009) 10770. [30] L. Zhang, S. Zybin, A.C.T. van Duin, S. Dasgupta, W.A. Goddard, E. Kober, J. Phys. Chem. 113 (2009) 10619.