Abstract
The simulation of chemical reactions and mechanical properties including failure from atoms to the micrometer scale remains a longstanding challenge in chemistry and materials science. Bottlenecks include computational feasibility, reliability, and cost. We introduce a method for reactive molecular dynamics simulations using a clean replacement of non-reactive classical harmonic bond potentials with reactive, energy-conserving Morse potentials, called the Reactive INTERFACE Force Field (IFF-R). IFF-R is compatible with force fields for organic and inorganic compounds such as IFF, CHARMM, PCFF, OPLS-AA, and AMBER. Bond dissociation is enabled by three interpretable Morse parameters per bond type and zero energy upon disconnect. Use cases for bond breaking in molecules, failure of polymers, carbon nanostructures, proteins, composite materials, and metals are shown. The simulation of bond forming reactions is included via template-based methods. IFF-R maintains the accuracy of the corresponding non-reactive force fields and is about 30 times faster than prior reactive simulation methods.
Similar content being viewed by others
Introduction
The prediction of chemical reactions and mechanisms for bond breaking, failure, and toughness of advanced materials from the atomistic scale to the microstructural and macroscopic length scale is a grand challenge in materials science1,2,3,4. The interplay of chemistry, defects, cross-linking in polymers, hierarchical assembly from the molecular scale to the macroscale and multi-scale dynamics of materials generates a wide range of chemical, optical, conductive, thermal, elastic, plastic, and failure properties. The development of bioinspired and engineered functional materials is limited by difficulties in probing associated patterns of chemical reactions and stress response mechanisms using experimental and computational methods4,5,6,7,8. Simulation methods, such as molecular dynamics (MD) simulations, can accelerate understanding of the underlying mechanisms from the scale of atoms to micrometers, provide guidance in materials selection and processing. For example, the role of compositions, defects, reactions, and processing conditions has been explored utilizing feedback loops from experiment to simulation and back to experiment to shorten discovery and development times9,10,11,12,13,14,15,16,17,18. Current limitations include the modeling chemical reactions, restrictions in size ( < μm) and computation time ( < ms). In this contribution, we introduce a method for modeling chemical reactions, especially bond breaking and failure of materials, without sacrificing the high accuracy, relatively short computation time and large system size in classical MD simulations.
Many current force fields for MD simulations in the materials and biomolecular community use harmonic bond energy terms, in addition to other bonded and nonbonded terms, and are limited in predictions of bond-breaking and bond-forming reactions (IFF5, PCFF19, AMBER20, CVFF21, CHARMM22, OPLS-AA23, COMPASS24, DREIDING25)14,26. Elements and specific chemistries are represented by atom types and often include multiple atom types for the same element. In some instances, such as metals, oxides, and ionic solids, only nonbonded terms are necessary and reactivity can be implemented with relative ease15,17,18. For example, alloying10,27, cis-trans photoisomerization reactions14, thiol-metal surface linking16,28, and hydration reactions18,29,30,31 can be sufficiently described without the need for explicit changes in bonding. The dissociation and formation of covalent bonds, however, requires new approaches. To introduce reactivity for this class of potentials, we consider the INTERFACE force field (IFF) which, among all-atom bonded force fields, covers the most diverse chemical space across the periodic table in highest accuracy and compatibility5. IFF has been developed to exclusively use interpretable parameters, accurately represent chemical bonding, and reproduce the structural as well as the energetic properties of included compounds under standard conditions relative to experimental data and theory, exceeding DFT methods up to an order of magnitude in accuracy17,32,33. IFF was conceived to be transferable and compatible with existing harmonic force fields for organic molecules and biopolymers by using the same potential energy terms, including 12-6 and 9-6 options for the Lennard-Jones potential5. IFF covers, for example, metals, minerals, oxides, 2D materials, common water models, solvents, gases, organic compounds, and performs exceptionally well for mixtures, solid-electrolyte interfaces, and hybrid materials. Lattice parameters and densities of pure phases usually deviate <0.5%, surface and hydration energies <5%, and mechanical properties <10% from reference data, which is about one to two orders of magnitude more accurate than other force fields for the same compounds5,17,18,33. The reliability exceeds common electron density functionals34,35,36,37,38 especially with regard to simulating interfaces17,33. IFF can be merged with biomolecular force fields such as CHARMM22, AMBER20, OPLS-AA23, as well as PCFF19,39 and COMPASS24, extending access to protein, DNA, lipid, and carbohydrate based multifunctional materials.
To-date, bond order potentials such as the Reactive Empirical Bond Order (REBO) Potential (REBO)40 and the Reactive Force Field (ReaxFF)41 have been developed to address bond dissociation in large-scale MD simulations, requiring a high number of additional energy terms. Specifically, ReaxFF uses Pauling’s bond order concept to describe the strength of a bond based on its chemical environment and interatomic distances41. ReaxFF can predict bond breakage and accurately reproduce reaction barriers and reaction energies derived from quantum mechanics for certain reactions42. However, the potential function contains more than triple the amount of energy terms compared to IFF or CHARMM, including extensive bond order terms, over-coordination, angle, dihedral, conjugation, lone pair, H-bond effects, and correction terms43. These terms add the ability to describe complex, concerted reaction paths and, at the same time, come at the cost of significant complexity, which can be difficult to interpret and which increase the computational expense. In contrast to IFF and biomolecular force fields, ReaxFF employs only one atom type for each element. This choice greatly simplifies the simulation setup by not requiring atom typing, however, it is also challenging for the force field developer and users. For example, the simulation of multiphase materials such as oxides in contact with water and oxygen-containing polymers requires multiple oxygen types. Therefore, different ReaxFF ‘branches’ have been developed, such as the combustion and the aqueous phase branches42. ReaxFF also does not cover complex polymers such as proteins, DNA, and their interfaces with various engineering materials with relatively high accuracy44, and it is desirable to combine the simplicity of a harmonic potential with the capabilities of ReaxFF.
As chemical reactions are highly sensitive to bond polarity, molecular geometry, and interaction energies, IFF is uniquely prepared to incorporate reactive processes. In this contribution, we augment IFF into the reactive Interface force field, IFF-R, to facilitate simulations of bond breaking and formation reactions. IFF-R uses a Morse potential to quantitatively represent the bond energy between pairs of atoms i and j as a function of distance, as opposed to a harmonic potential in IFF, and density functional theory (DFT) or Møller-Plesset perturbation theory (MP2) which feature unphysical energy maxima upon dissociation. The Morse representation aligns with experimentally measured energy functions and has a quantum mechanical justification. (Fig. 1a)45,46. Accordingly, we substituted the harmonic bond energy term in the energy expression with a Morse bond energy term, including either all types of bonds or a subset thereof, e.g., only bonds of lowest energy known to dissociate before others (Fig. 1b, c). We demonstrate that this clean replacement requires no changes in other force field parameters to add bond breaking capabilities, maintains the full benefits of the non-reactive IFF, and speeds up reactive simulation ~30 times relative to existing methods. We also demonstrate the formation of bonds, or reversibility of bond breaking with IFF-R using template-based reaction simulations, specifically, the REACTER toolkit47. The formation and breaking of bonds is thus possible in a continuous simulation without the need for the more complex bond-order potentials. We demonstrate the application of IFF-R with PCFF and CHARMM functional forms (Fig. 1d).
a The non-reactive harmonic bond potential in IFF (gray) is substituted for a dissociative Morse bond potential (pink) to yield IFF-R, focusing on chemical bonds at the lower spectrum of bond strength, or on all atoms in the system. The graph shows the bond energy vs. bond distance for the resulting Morse potential and the equivalent harmonic bond potential, as well as a bead-spring model for the resting state (at r0), the stretched state (upper middle), and dissociated bond states, respectively (far right). The Morse potential in IFF-R is tuned to fit the harmonic bond potential near equilibrium (α) and the bond dissociation energy (D) without changes in other, not bond-related force field parameters. The same modifications can be applied to other harmonic force fields (CHARMM, AMBER, OPLS-AA, DREIDING, CVFF, PCFF), given a certain level of thermodynamic consistency and interpretability of the parameters. b Changes in mathematical form and parameters for the harmonic and Morse bond potentials. Two parameters r0,ij and Kr, ij from the harmonic parent force field are extended into three parameters in the Morse potential, of which all can be obtained from experimental data. If r0,ij and Kr,ij are known from a parent force field, only Dij needs to be determined using additional information from experimental databases51,52,53,54 or from quantum mechanical calculations, preferably at the MP2, CCSD, or higher level55. The Morse bond equation 1 is used for parameterization purposes only and the Morse bond equation 2 is strictly to be used in simulations to ensure a transition to zero bond energy during dissociation at larger atomic displacements. c List of parameters, interpretation, and data sources. d Workflow diagram to convert IFF into IFF-R using Morse potentials for specific chemistries.
Results
In the following, we describe the underlying theory, the workflow of parameter additions, examples of bond-breaking and stress-strain curves up to failure for validation, and the incorporation of bond forming reactions. The examples cover bond dissociation curves for small molecules, their representation in IFF-R, the prediction of stress-strain curves up to failure for multiple materials classes including carbon nanotubes (SWCNT), syndiotactic polyacrylonitrile (PAN), cellulose Iβ, spider silk protein, polymer/ceramic composites, and γ-iron. The computational speed is at least 30x higher than current state-of-the-art bond-order potentials. We also verify that IFF-R maintains the high-level performance for bulk and interfacial properties characteristic of IFF parameters, including validation of mass densities, vaporization energies, interfacial energies, and elastic moduli. We compare these critical computed pre-reaction and post-reaction properties between IFF-R and ReaxFF. We share an initial graphical user interface to aid in input preparation and simulation setup. We discuss limitations and potential extensions of IFF-R. Hereby, our contribution focuses on the introduction of the method and providing sufficient evidence that the described IFF to IFF-R conversion can be applied to any carefully parameterized class I or class II force field to facilitate accurate simulations of bond breaking and formation. A large database of use cases, exhaustive lists of parameters, and optimizations of code are beyond the scope of this work.
Description of the approach
The simulation of chemical reactions, especially in heterogeneous materials, has remained a major challenge in molecular dynamics simulations. The replacement of harmonic chemical bonds by Morse bonds (Fig. 1a–c) in a simulation offers a simple and interpretable description of bond dissociation without the need for complex fit parameters, and has not been comprehensively realized in reactive force fields to-date48,49,50. To obtain the Morse parameters for a specific reactive bond, i.e., a pair of covalently bonded atom types ij, we can first plot the bond energy vs atomic distance using the harmonic bond potential, as known from experiment or the nonreactive potential (Fig. 1a, d). Next, using experimental data51,52,53,54 or high-level quantum mechanical estimates (CCSD(T), MP2, possibly DFT)55 for the respective bond dissociation energy Dij (Fig. 1c, d), we plot an initial Morse potential curve (Fig. 1b) for atoms i and j. Thereby, the equilibrium bond length ro,ij remains the same as in the harmonic potential, or as known from experiment. The parameter αij can be used to fit the Morse bond curve to the harmonic bond curve near the resting state (Fig. 1c, d) and typically in a narrow range of 2.1 ± 0.3 Å−1 21. The parameter αij can be refined to match the wavenumber of bond vibrations in the simulations to experimental data from Infrared and Raman spectroscopy56.
Typical parameter ranges are 0.5 to 2.5 Å for equilibrium bond lengths, 1.8 to 2.4 Å−1 for the width α (whereby larger values are equal to higher vibrational wavenumbers and a smaller width), and 0 to 250 kcal/mol for the bond dissociation energy. Examples of Morse bond parameters are shown in Fig. 2 and in Table S1 in the Supplementary Information. The remaining force field parameters in IFF or other force fields (e.g., PCFF used here for some organic molecules, AMBER, CHARMM) can be retained without changes, including the atomic charges, Lennard-Jones parameters, angles, torsions, out-of-plane terms, and other non-reactive harmonic bonds. For simplicity of parameterization, Python code is provided in the Supplementary Software (see Figure S1 in the Supplementary Information for an overview) so that any available bond in IFF (or other class I and class II force fields) can be replaced with a Morse bond potential. The provided code performs an empirical fitting of the Morse bond potential to a given harmonic potential and may not represent the same level of IFF accuracy that will be discussed in the following subsections.
We consider the dissociation of several bonds and bond pairs in two amine monomers, 3,5-diethyl-toluene-1,6-diamine (DETDA) a, b and 4,4’-diaminodiphenylsulfone (4,4’-DDS) c–f, indicated by red dashed lines. The experimentally determined bond dissociation energy D for similar chemistry is shown as a horizontal line and labeled in the bottom right corner of the panels (refs. 51,52,53,54,63). These values, however, are somewhat approximate and do not consider the exact substitution pattern. Chemically specific bond scans were obtained using simulations at the level of 2nd order Møller-Plesset perturbation theory (MP2), density functional theory (DFT) with the functional B3LYP, and simulations with ReaxFF (ref. 60). The drop in bond energies between 2.5 and 3 Å when using MP2 and DFT calculations shows an unphysical limitation that is related to electron delocalization errors and affects all reactive simulations at MP2 and DFT levels of theory (refs. 61,62). IFF-R parameters were informed by experimental, MP2 and DFT data, and two options with different values for the bond vibration constant α are shown (solid and dashed blue lines). (a, b) Bond dissociation curves for Car – Car in DETDA in two distinct molecular planes. (c, d) Bond dissociation curve of Car – Car bonds in 4-4’-DDS in two distinct planes. (e, f) Bond dissociation curve of Car – S and Car – N bonds in 4,4’-DDS. Source data are provided as a Source Data file.
Further validation and quantum mechanical simulations indicate that the dissociation of reactive bonds can be considered complete between approximately 140% and 170% of the equilibrium bond length (Figures S2 to S6 and Section S3 in the Supplementary Information). We recommend removing topological connections between atoms i and j at approximately 200% of the equilibrium bond length to minimize artificial correlations in the simulation near and past failure. Following bond disconnection, different atomic charges and atom types may be assigned at the reaction site to maintain charge neutrality and accurately represent the reaction products.
Morse parameterization and customization of IFF-R
Bond dissociation curves with accurate values for the bond dissociation energy D are essential for reliable predictions of bond breaking mechanisms and mechanical properties in complex materials. Hereby, tabulated “bond dissociation energies” for specific bonds from the literature and available databases52 are most suitable for IFF-R. “Average bond dissociation energies” are less suitable as they represent an average bond dissociation energy for all equivalent bonds in a molecule, which typically differ from the energy required to break the first bond, second bond, and so forth (e.g. multiple C-H bonds in CH4).
As an example, we compared bond dissociation energies obtained using different methods for diethyltoluenediamine (DETDA) (Fig. 2a, b) and 4,4’-diaminodiphenylsulfone (4,4’-DDS) (Fig. 2c–f). The molecules are used in the synthesis of polyurethanes, polyimides and epoxy resins to manufacture high-strength CNT/epoxy composites57,58,59. We compare data from experiments (horizontal lines), 2nd order Møller-Plesset perturbation theory (MP2), density functional theory (DFT) with B3LYP, molecular dynamics simulation with tuned curves in IFF-R and ReaxFF (parameters from ref. 60). The bond scans generally look similar but differ in detail. While experimental data are often an average for a given type of chemical bond, MP2 calculations resolve substituent effects and are more accurate than plane wave DFT methods. As a limitation, energy profiles at the DFT and MP2 levels invert toward lower energies at bond distances of 2.5 to 3 Å, as opposed to convergence to zero energy (Fig. 2a–d). The unphysical curve shape is related to electron delocalization error (see further examples in ref. 61,62) and impacts the reliability of DFT and MP2 methods in all instances where bonds form or dissociate, especially in catalysis such as CO2 reduction, water splitting, fuel cells, and nitrogen fixation. The Morse potential in IFF-R as well as higher-level quantum mechanical calculations such as CCSD(T), if computationally feasible62, overcome these errors. Bond scans by ReaxFF appear somewhat irregular and less aligned with experimental or quantum mechanical results. IFF-R curves were tuned to a combination of experimental, MP2 and DFT values and yield smooth curve shapes (Fig. 2a–f).
Highly resolved quantum mechanical calculations such as MP2, or preferably CCSD(T), can be used to distinguish differences in bond dissociation energies relative to average values for Car-Car, C-N, and C-S bonds from experiments51,52,53,54,63. Substituent effects, such as inductive and mesomeric effects, can be incorporated into IFF-R by choosing α and D parameters that match MP2 or other high-level results. For example, in DETDA, the bond dissociation energy of the 1, 6 and 3, 4 bonds (MP2 in Fig. 2a) is ~10 kcal mol-1 lower than that of the 2, 3 and 5, 6 bonds (MP2 in Fig. 2b). Hereby, the 1, 6 bonds and 3, 4 bonds have less double bond character compared to the 2, 3 and 5, 6 bonds according to an analysis of the resonance structures (Figure S7 and Section S4 in the Supplementary Information). In 4,4’-DDS, a push-pull effect by the NH2 and SO2 groups yields nearly full double bond character for the 2, 3 and 5, 6 bonds (MP2 in Fig. 2c), which have a ~ 30 kcal mol-1 stronger bond dissociation energy than the 1, 2 and 4, 5 bonds (MP2 in Fig. 2d)64. We did not incorporate specific adjustments in atom types and bond energies for substituent effects here, however, the data illustrate how IFF-R can include fine detail of the electronic structure when necessary. The MP2 bond scans also reproduce different bond strengths for Car-S and Car-N bonds in comparison to Car-Car bonds, which are 75-85 kcal mol-1, 110 kcal/mol, and 120-170 kcal mol-1, respectively, according to experiments (Fig. 2e, f).
In summary, IFF-R can consider all chemical bonds as breakable or focus on a sub-set of bonds, such as bonds of lowest strength only, and represent the remaining bonds as harmonic. All three Morse parameters for each bond type can be conveniently obtained using extensive databases with thousands of experimental measurements (Table S1 in the Supplementary Information), as well as using chemical analogy for similar bonds51,52,53,54,63. IFF-R accuracy can be further improved by tuning parameters based on observations from high-level quantum mechanical calculations. Thereby, IFF-R can account for specific chemical environments and electronic structure effects as needed.
Applications of IFF-R to compute stress-strain curves and analyze failure
Stress-strain curves are critical to characterize mechanical properties and were simulated for different materials classes using IFF-R, including single wall carbon nanotubes (SWCNTs), polyacrylonitrile (PAN), cellulose Iβ, and γ-iron (Fig. 3). Atomic-level deformations and the occurrence of bond scission upon failure can be clearly seen, which non-reactive simulations cannot predict. The comparison to experimentally measured Young’s modulus and tensile strength65,66,67,68,69, or high-level quantum mechanical results if not available, shows good agreement (Fig. 3a, b, d, f, h and Table 1)3,66,67,70,71,72,73. The computed mechanical properties from IFF-R agree with most reference data within 10%.
When available, the experimentally determined tensile strengths, or high-level theory estimates, are shown as horizontal lines a, b, d, f, h (refs. 65,66,67,68,69). a, b Single wall carbon nanotubes (SWCNT) of 1 nm diameter and comparison to simulations with ReaxFF (refs. 43,74,76,77,78). Models of the fractured CNTs are shown in the top right corner of a and top left corner of b (using ref. 76. in b). Virtual π-electrons are used for the IFF-R CNT model. (c) Crack propagation demonstrated in a SWCNT using IFF-R. d Syndiotactic, crystalline polyacrylonitrile (PAN). e Bond scission, progression of deformation, and failure in PAN. f Cellulose Iβ. g A broken chain in a strained cellulose Iβ crystal and evolution of cellulose Iβ failure. h Crystalline γ-iron. Strain was applied in the plane perpendicular to the (111) surface. i Plastic deformation and ductile failure in γ-iron. Source data are provided as a Source Data file.
The reliability is, on average, notably better and more consistent compared to ReaxFF, where typical deviations are 10% to 50% and exceed 1000% for iron (Fig. 3 and Table 1). ReaxFF required the use of multiple parameter sets74,75,76. The choice of the ReaxFF parameter set is somewhat arbitrary as the parameter sets were not explicitly parameterized for the purpose of uniaxial tensile simulations, and in some cases not even for the specific materials of interest. Therefore, the accuracy varied widely (Fig. 3b, d, f, h). Given the difficulty of developing ReaxFF parameters, it was necessary to use discretion to choose the best available ReaxFF parameters as a comparison. For convenience to the user, IFF-R and IFF, use one well-defined, interpretable parameter set that avoids the ambiguity of multiple parameter sets, and the concept of thermodynamic consistency does not require specific fitting to compute mechanical properties43,60,74,75,76,77,78,79,80. Distinct functional groups in IFF-R and IFF utilize different atom types of the same element as necessary.
Regarding failure mechanisms, single wall carbon nanotubes (SWCNTs) often showed smooth cleavage planes (Fig. 3a–c). The SWCNT fails by strain-induced bond dissociation propagation and reaches a strain of ~20% at break in the absence of defects (Fig. 3a, c). Mechanical failure of the polymers PAN and cellulose Iβ arises from breaking of load-bearing bonds in single polymer chains (Fig. 3d–g). Dissociation initially affects single chains, involves chain sliding, and then propagates further through the crystal (Fig. 3e, g). The IFF-R model of γ-iron shows the evolution of microcracks and plastic deformation before failure (Fig. 3h, i). Slip and crack propagation along partially preserved crystal planes is seen, leading to moderately ductile failure.
IFF-R produces continuous, non-linear stress-strain curves and reliable mechanical properties. The Young’s modulus and tensile strength values for the SWCNT are comparable to experiments of a CNT with similar dimensions (Fig. 3a, b)67. For 100% crystalline PAN (Fig. 3d, e), we found no experimental reference data for the tensile strength, and the elastic modulus was estimated to be 172 GPa from quantum mechanical DFT (GGA PBE) simulations using the Cambridge Serial Total Energy Package (CASTEP) program, which agrees with IFF-R data (Table 1)81. ReaxFF with both the Liu and Damirchi parameters (refs. 74,76) suggests similar stress-strain properties, although greater deviation from the DFT (GGA PBE) values are seen. The modulus of cellulose was computed about 20% lower than in experiment and the strength was overestimated by 20%, which is in better agreement and less ambiguous than the corresponding ReaxFF predictions (Fig. 3f, g and Table 1). Cellulose parameters were used as given in PCFF and modified to utilize the Morse potential following the IFF-R protocol (Fig. 1). Updating PCFF parameters for cellulose for higher IFF and IFF-R quality may require the representation of stereo-electronic effects and updated nonbond parameters and was not further considered here.
ReaxFF parameters for the γ-iron model led to rather divergent results, including more than twice the experimental strength and extreme, unphysical ductility of 40%. For γ-iron and other metals, IFF-R does not require the definition of Morse bonds due to the simplicity of the Lennard-Jones (LJ) potential and is equal to IFF without modifications. For most chemistries, including metals and oxides10,17, Young’s Modulus tends to be better predicted using IFF and IFF-R with 12-6 LJ options (same energy expression as CVFF, AMBER, CHARMM, OPLS-AA, DREIDING) than with equivalent 9-6 LJ parameters (same energy expression as CFF, PCFF, and COMPASS).
We further tested the simulations of stress-strain curves up to failure with IFF-R for less idealized, more realistic nanostructures and biopolymers (Figures S8 and S9, and Section S5 in the Supplementary Information). These include SWCNTs with various levels of vacancy defects, which successively reduce the tensile strength and ultimate strain by ~20% (Figure S8 in the Supplementary Information). We simulated the failure of the spider silk protein (spidroin) at two different temperatures, demonstrating good agreement of the computed Young’s modulus (3.6 GPa) and tensile strength (1.6 GPa) with experimental data (Figure S9 in the Supplementary Information). For detailed quantitative studies of the mechanics and failure mechanisms of complex materials, usually a large set of simulations using multiple morphologies, relaxation times, and strain rates is necessary.
Computational efficiency
A major benefit of simulations with IFF-R, besides accuracy and interpretability, is a high computational speed and low memory usage relative to alternative methods. The advantage over quantum mechanical simulations and DFT-machine learned potentials is obvious, about 106 or 104 times faster, respectively82. A comparison of the computational speed of IFF-R and ReaxFF using the open-source code LAMMPS shows on average about 30 times faster speed for systems sized between 1000 and 10,000 atoms and needs only about 10% of the memory (RAM) (Fig. 4). Therefore, IFF-R allows access to systems of larger size and more extensive dynamics. To arrive at this benchmark, we measured the computation time for the above stress-strain curves using molecular dynamics simulation in the NPT ensemble, normalized to 1 CPU with a typical time steps of 1 fs for IFF-R and 0.5 fs for ReaxFF. IFF-R is approximately 60x faster for simulations of metals, about 26x faster for small polymeric systems, and 8x faster for carbon nanotubes, which include virtual π electrons, when compared to simulations with ReaxFF (Fig. 4)74. IFF-R is faster by about the same magnitude for larger systems with over 100,000 atoms (Figure S10 and Section S6 in the Supplementary Information). The reason for the speed up lies in the simpler energy expression with an order of magnitude fewer terms, including constant atomic charges, predefined bonded terms, and approximately two orders of magnitude fewer algorithmic operations per unit time overall.
The models were 600 to 14,000 atoms in size and include γ-iron, crystalline cellulose Iβ (ref. 119), crystalline syndiotactic poly(acrylonitrile) (PAN), and a bundle of single-walled carbon nanotubes (SWCNTs) of 1 nm diameter (which includes virtual π electrons in IFF-R). Snapshots at failure are shown and the speed-up in multiples of CPU time by IFF-R is indicated. The benchmark simulations were carried out on 1 CPU, using 1000 steps of molecular dynamics simulations with a constant strain rate in the tensile direction. The isobaric-isothermal ensemble (NPT) was employed with pressure control in x and y directions, and no pressure control in the tensile (z) direction. IFF-R completes 1 ns molecular dynamics simulation on 1 CPU in less than 12 h for all 4 materials while ReaxFF requires 2 to 4 days for γ-iron and SWCNTs, and between 6 days and 1½ weeks for (bio)polymers. The relative CPU times remain similar for model systems over 100,000 atoms (Figure S10 in the Supplementary Information). Source data are provided as a Source Data file.
For the comparisons, exact numbers depend on the simulation code and settings of both IFF-R and ReaxFF. For example, 1 fs is a conservative time step for IFF-R and time steps of 1.5 fs or 2 fs can often be used along with the “shake” algorithm for H atoms and virtual π electrons83. The average speedup can then double to approximately 50x. Simulations can further be accelerated by using simulation codes other than LAMMPS, which is one of the slower open-source codes for both IFF/IFF-R and ReaxFF. IFF/IFF-R is most efficiently used with the open-source code GROMACS, which runs 3 to 6 times faster than LAMMPS84,85, offering approximately 100x speed-up. ReaxFF runs multiple times faster with the Amsterdam Modeling Suite (AMS), which is available commercially. The efficiency of ReaxFF in LAMMPS is strongly dependent on its implementation42. In this study, the “reax/c” LAMMPS implementation was used which is faster than “reax”86. A major bottleneck in computational speed of ReaxFF is the handling of atomic charges76. This study used the “fix qeq/reax” method, however, some faster alternatives may exist and are not tested here (e.g., “fix qeq/fire”)87,88. The relative speedup can also be affected by use of GPUs rather than CPUs.
Compatibility and transferability
We determined that the replacement of the harmonic bond potential in IFF with the Morse bond potential in IFF-R maintains the accuracy of all physical, chemical, and interfacial properties of IFF5, while adding the desired bond dissociation capabilities (Table 2). The key functionality consists in the unique compatibility to examine simulate mixed inorganic, organic, biomolecular, composite, and electrolyte model systems without additional fit parameters. Consistent with IFF theory5, validation of such functionality includes (1) the evaluation of correct structural properties via lattice parameters and the mass density of individual compounds, as well as (2) the evaluation of relative energy differences via surface or vaporization energies, both at 298 K. (3) In addition, compatible models require physically correct, interpretable representations of chemical bonding, which achieved by using the same atomic charges in IFF and IFF-R, and physically realistic Morse potentials. Altogether, keeping the same accuracy of predicted properties before and after the replacement of the harmonic potential with a Morse potential ensures that the key functions of the Hamiltonian are maintained, which is to predict structures and associated changes in energies5.
Specifically, the predicted lattice parameters, mass density, and surface energy of graphite using IFF-R remain the same as in IFF within statistical errors of ±0.5% (Table 2). The calculated lattice parameters and the surface energy are close to experimental data, and nearly equivalent for IFF and IFF-R. The same is true for the mass density and vaporization enthalpy of propionitrile, a small molecule most closely related to PAN, which is liquid at 298 K as opposed to solid graphite. The computed mass densities are within 1% of experimental data and comparable between IFF and IFF-R (Table 2). The simulated vaporization enthalpies of propionitrile are a few percent higher than experiment yet remain the same using IFF and IFF-R within the statistical uncertainty. α-D-glucose, the monomer of cellulose, shows nearly identical properties using IFF and IFF-R (Table 2). γ-iron shows identical performance for IFF and IFF-R, where covalent bonds were not necessary to facilitate reactivity (Table 2).
For α-D-glucose, we used “PCFF-R” parameters that correspond to PCFF parameters with addition of a Morse potential for the weakest C-O bonds, analogous to the extension of IFF parameters to IFF-R (Figure S11 in the Supplementary Information). PCFF and PCFF-R are of much lower accuracy than IFF, even though we included some improvements in the distribution of atomic charges (Table 2). The lattice parameters and vaporization (sublimation) energies computed by PCFF-R and PCFF match, and systematic errors relative to experimental data remain the same89. In terms of absolute values, the lattice parameters a and b exhibit 10% to 15% shrinkage and the c parameter expands by 36%, which is far from IFF standards ( < 0.5%)5, and only the mass density and vaporization energy are comparable to experiment. Better force field parameters for glucose require extensive analysis and changes, such as stereo-electronic effects at C1 carbon atoms and more interpretable bonded and nonbonded parameters (Figure S11 in the Supplementary Information).
The agreement in 3D structures and relative energies of diverse compounds in IFF and IFF-R, as well as in PCFF and PCFF-R, confirms that a modification of harmonic bonds to Morse bonds maintains original accuracy and is efficient to introduce reactivity. The IFF-R protocol can also introduce reactivity into other force fields that use the same energy expression as IFF, e.g., CHARMM, AMBER, OPLS-AA, DREIDING, and PCFF5.
Comparison of the accuracy of IFF-R and ReaxFF
IFF-R typically computes structural parameters within 0% to 1% deviation relative to experiments for validated compounds, and ReaxFF with deviations of 0% to 5%17. In addition, a Hamiltonian foremost aims to reproduce energy differences associated with structural differences under standard conditions, such as surface energies (cleavage energies) or vaporization energies under standard conditions5. Mechanical, thermal, and reactive properties can only be meaningfully validated after this condition is met as they are defined as higher order derivatives of the energy with respect to interatomic distances and temperature. These concepts of thermodynamic consistency are integral to IFF-R and are often neglected in other interatomic potentials, making them less suitable and less extensible for the simulation of reactions. For example, some force fields reproduce structural and mechanical properties while energy differences exceed 100% deviation from experimental observations7.
IFF-R reproduces the surface energies and vaporization energies for various materials with less than 5% deviation with experiment (Table 3). For graphite, virtual π-electrons are essential to simultaneously reproduce surface, wetting, and adsorption properties90,91. Deviations using ReaxFF exceed 30% for graphite and propionitrile (Table 3). Mechanical properties such as modulus and tensile strength deviate in the 10-20% range from experiment with IFF-R, compared to a 50% range and maximum deviations over 100% with ReaxFF (Table 1) (see detailed comparisons in Section S7 in the Supplementary Information).
IFF-R performs better overall than ReaxFF for validated compounds, works reliably with a single parameter set, and can be applied to any mixture of materials (Tables 1 and 3). Further validation of IFF-R for predicting the thermomechanical properties of amorphous epoxy, benzoxazine, PEEK, and polyamide resins has recently been reported by Odegard et al. 92,93,94,95.
Disconnection of Morse bonds after dissociation
IFF-R is suited to examine bond breaking mechanisms and analyze stress-strain curves up to failure. The Morse potential alone, however, does not disconnect the specified bonds upon elongation. Therefore, we systematically tested the disconnection of Morse bonds upon elongation for a range of cutoff distances, as well as the role of the shifted Morse potential using the program LAMMPS (Figures S3 to S5 and Section S3 in the Supplementary Information)57. As a result of our tests, we recommend a cutoff at a bond elongation of about 200% of the equilibrium bond length. Depending on the bond type, a minimum elongation of 140% to 170% was required for physically meaningful results, and cutoffs much beyond 200% of the equilibrium bond length (e.g., 300%) introduce artificial energy contributions from the remaining bonded and nonbonded terms. At approximately 200% of the equilibrium bond length, typical Morse bond parameters overcome most of the attraction toward the minimum bond force (Fig. 2). An additive shift of bond-specific Morse potentials to zero at the bond cutoff distance supports energy continuity when Morse bonds are disconnected and reduces discontinuities in other energy contributions (Fig. 1). Upon disconnection of bonds, atomic charges and atom types can be reassigned to those of the reaction products, whereby keeping charge neutrality at the reaction site is essential (Figure S2 to S6 and Section S3 in the Supplementary Information). When simulations extend beyond bond disconnection, it is also possible to adjust chemistries for follow-on reactions, e.g., using reaction templates (Figure S6 in the Supplementary Information).
In comparison, we recall that quantum calculations at DFT and MP2 levels lead to unphysical shapes of bond dissociation curves and have potentially serious errors at around 200% of the equilibrium bond length (Fig. 2a-d)61,62. Therefore, DFT and MP2 data on bond dissociation are not suitable to train IFF-R beyond estimates of bond dissociation energies. IFF-R, in addition to good interpretability, is also promising to outperform reactive simulations based on machine learning with DFT and MP2 data, which incorporate the flawed bond dissociation curves (Fig. 2a-d). Trustworthy guidance for bond formation or dissociation from ab-initio methods requires a theory level of CCSD(T) or higher.
Simulation of complete chemical reactions including bond formation, bond dissociation, and specific stoichiometries
In a similar way, the reversible process of forming bonds was implemented as reactive groups dynamically approach each other during the simulation. The physical basis to decide on bond connections or the occurrence of multi-step chemical reactions can involve local screening functions below a threshold distance, followed by changes in force field types, bond connectivity, and energy minimization as described by Barr et al.96. Similar concepts were implemented in the REACTER framework in LAMMPS, for example, to polymerize nylon-6,6’ and polystyrene with up to 99% cross-linking47. These methods can also be used to simulate bond cleavage, however, IFF-R with a shifted Morse potential better reproduces the characteristic nonlinearity of the bond energy before failure (Fig. 1), eliminates discontinuities in bond energy, requires less parameters and less user input. IFF-R is compatible with REACTER to handle bond formation and bond breaking in the same simulation (Figs. 5 and 6).
Tensile stress was applied to a model of aligned polyacrylonitrile (PAN) chains up to failure, followed by compression and allowing bond forming interactions using the reactive INTERFACE force field (IFF-R) and REACTER47. Bond breaking and bond formation occurred in a single, continuous simulation. a An initial, energy-minimized structure of 9 aligned PAN chains. 4 chains are shown in color to highlight the formation of 4 corresponding 4 PAN oligomers after chains failure and allowing reconnection of new bonds. b Bond failure during application of tensile stress in molecular dynamics simulation with IFF-R. c Compression of the disconnected structure from (b) and allowance of bond formation between broken chain ends using REACTER. Red, green, blue, and purple chains reacted with uncolored chains and formed 4 new polymer strands. A simplistic assumption of reversible bond formation was made for demonstration purposes. Specific mechanisms and statistical reactivity can be incorporated in detailed simulations.
a The two monomers bis[4-(3-aminophenoxy)phenyl] sulfone (m-BAPS), a difunctional amine, and bisphenol A diglycidyl ether (DGEBA), a difunctional epoxide, were polymerized using the reactive INTERFACE force field (IFF-R) and REACTER47. The primary and secondary amine reactions were included. b A cluster analysis to visualize the number of independent molecules of the pre-reacted model (upper left), cross-linked neat polymer structures (upper right), and cross-linked polymer structures in the presence of a CNT (bottom left and right). c Snapshots of the cured epoxy resin in b during the strain simulation. (d) Snapshot of failure of the CNT/epoxy composite in b. e Stress-strain curve (on logarithmic scale) for the cross-linked neat m-BAPS-DGEBA network, the CNT/epoxy composite, and a single CNT model. Source data are provided as a Source Data file.
As a demonstration, we carried out reactive MD simulations to dissociate carbon-carbon bonds in polyacrylonitrile (PAN) and subsequently allow bond formation via a hypothetical self-repair mechanism (Figs. 3d, e and 5). First, nine aligned PAN chains were exposed to strain during IFF-R MD simulation until failure (Fig. 5a, b). Subsequently, the model system was compressed and the fragments were allowed to simultaneously bond using the REACTER toolkit (Fig. 5c)47. While the mechanism is oversimplified here (and can be customized for specificity), the simulation demonstrates that IFF-R and REACTER are compatible to simulate bond cleavage and bond formation in a continuous simulation. The methods can be used, for example, to simulate self-healing of polymers.
As a second example, we simulated the polymerization and stress-induced failure of aerospace-grade epoxy polymers and a CNT/epoxy matrix composite (Fig. 6). First, a stoichiometric mixture of epoxy and amine monomer precursors (Fig. 6a) containing approximately 34,000 atoms was subjected to polymerization using IFF-R and REACTER47 using the LAMMPS program (Fig. 6b). Cross-linking simulations were carried out for the pristine monomers and for the monomers in the presence of a carbon nanotube. The polymerization reached ~84% cross-linking in both cases, whereby less than full conversion agrees with experimental observations due to limitations in diffusion of monomers and oligomers. Every molecule was color-coded and the progress of the reaction towards the cross-linked matrix can be seen as a color change (dark red color for cross-linked matrix molecules, due to superposition it may not appear quantitatively in the 2D view) (Fig. 6b). Discontinuous matrix coverage near the top of the CNT reinforcement is seen in the side view. The models can be upscaled beyond 105 atoms using IFF-R, which would be difficult using quantum calculations and bond order potentials. We used predefined reaction templates under standard pressure to form the cross-links (see Methods).
After the reactions reached completion, the matrix and composite models were subjected to MD simulations with an applied strain program until failure (Fig. 6c–e). Stress-strain curves were recorded in single continuous simulations, showing clear differences in modulus, tensile strength, and strain at break. Failure mechanisms involved void formation, polymer-CNT bonding, and stress concentration in areas with minimal thickness of the matrix (Fig. 6c, d). For the neat polymer matrix, the stress-strain curve shows an elastic region at low strain (0.0 to ~0.3), followed by yielding up to a strain of approximately 1.5, plastic strain hardening up to approximately 4.3, and failure with a strain at break of ~5.5 (Fig. 6c, e). At a nominal strain of 3.0 (Fig. 6c, top), the chains in the polymer network aligned with the strain direction (strain hardening) and at a nominal strain of 5.5, the polymer chains fully dissociated. The tensile strength was up to approximately 0.2 GPa, in the range of experimental data (Fig. 6e). In the case of the CNT/epoxy composite, CNT reinforcement led to bond dissociation in the area with minimal matrix coverage and highest stress concentration (Fig. 6d). The CNT/epoxy composite was stronger and less ductile than the polymer. Relative to the pristine CNT, the ductility was enhanced (strain at break 0.23) and the tensile strength decreased to ~6 GPa (Fig. 6d, e). The pristine CNT had a tensile strength of 115 GPa and a strain at break of 0.20, consistent with laboratory data (Table 1).
Accordingly, IFF-R and other harmonic force fields can simulate bond formation and bond dissociation without the computational complexity of a bond order force field or extensive reparameterization, respectively. During the reaction phase, we recommend using the harmonic bond parameters to avoid unintentional dissociation of some atoms, which can induce failure of the simulation. Morse bond parameters should be available from the outset to run simulations continuously and applied once the reaction has reached the desired degree of conversion. To aid in this setup, we developed a tool for Morse bond parameterization and a sample script “in.xlink” (Section 2 in the Supplementary Information and Supplementary Software).
Opportunities, potential pitfalls, and limitations
The proposed Morse bond modification can be applied to all force fields using a harmonic, or otherwise exponential bond potential, such as IFF, CHARMM, AMBER, PCFF, COMPASS, OPLS-AA, and DREIDING. When the parent force fields do not reproduce experimentally known properties such as lattice parameters or vaporization energies as seen for α-D-glucose and cellulose in PCFF (Table 2),we recommend a derivation of updated parameters following IFF protocol5 to improve the accuracy and transferability, which is of equal importance to adding Morse potentials.
In the derivation of the three Morse parameters for new reactive bonds, we recommend to strictly follow the physical interpretation, which includes the equilibrium bond length (r0ij), bond stiffness with associated vibration spectra (α), and bond dissociation energies (D) (Fig. 1). Experimental data as well as simulations with high-level quantum mechanical methods such as CCSD(T) and MP2 are suitable. Quantum mechanical simulations at simplified DFT level often have large errors and should be avoided for parameterization purposes97,98. Alternatively, Morse parameters could be empirically adjusted to match computed stress-strain properties to experimental reference data. We would give low priority to this approach as it may not correctly weigh contributions from other energy terms to the mechanical properties (angles, torsions, charges, LJ potential) and insert additional uncertainties from stress-strain measurements related to chemical composition and crystallinity of the samples into the force field.
Chemical reactions after bond dissociation can involve a variety of intermediates and products, including rearrangements with new chemical bonds, addition, elimination, and other multistep reactions (Figure S6 in the Supplementary Information). Morse bonds in IFF-R do not consider secondary reactions unless reaction templates with defined changes in bond connectivity are added. The creation of reaction templates remains specific for a given type of chemistry. Likely reaction mechanisms and their representation in molecular simulations need to be known a priori from chemistry knowledge, for example, simple one-step stoichiometric conversions to products, inclusion of intermediates, and multiple products. Parameters for reaction intermediates and products meeting IFF standards may be needed. Predictions of reaction mechanisms by computational methods such as quantum mechanics/DFT, ReaxFF, or IFF-R are hardly feasible since the details often depend on solvents, pH value, temperature, catalysts, byproducts, and time scale (e.g., seconds to hours)64,99. Synthetic chemists continue to rely on experimentally derived knowledge, enhanced by insights from computational methods.
Discussion
IFF-R is well suited for the analysis of deformation and failure mechanisms of complex biological and material structures from atomic to micrometer scales. IFF-R builds on an effective replacement of nonreactive harmonic bond potentials with reactive Morse bond potentials in the Interface Force Field (IFF) and in other harmonic force fields to allow the simulation of bond scission. The modifications are easy to implement via substitution of the harmonic bond expressions with 2 interpretable parameters by a Morse-bond expression with 3 interpretable parameters, which are shifted to zero energy at the bond cutoff. We call this force field IFF-R and other harmonic reactive force fields equivalently, e.g., CHARMM-R and PCFF-R. Morse bonds can be assigned to selected bonded atom pairs or to all bonded atom pairs, and the remaining force field parameters can be used as is without loss in accuracy. Bond formation during chemical reactions can be incorporated by template-based methods such as the REACTER framework, enabling the dissociation and formation of bonds in single, continuous simulations using IFF-R.
IFF-R is suitable to simulate non-linear stress-strain curves up to failure and predicting elastic moduli and tensile strengths for all material classes and multiphase materials across the periodic table with high accuracy, at low computational cost, and with high interpretability. Relative to IFF, IFF-R improves the accuracy of the representation of chemical bonding and retains the accuracy of computed physical, chemical, and interfacial properties, including lattice parameters, mass densities, surface energies, and elastic moduli. The parameters for Morse bonds can be derived from experimental and high-accuracy ab-initio data and, as an option, also include local electronic structure effects via customized atom types with specific bond dissociation energies.
IFF-R uses a single parameter set, simplifying the user-experience and avoiding the challenge of selecting from a multitude of differently performing parameter sets, as is common with electronic density functionals and ReaxFF. The accuracy is up to several times higher than other reactive force fields, such as ReaxFF. The computational speed is about 30 times faster when using LAMMPS and up to 100 times faster using more efficient MD codes like GROMACS. IFF-R provides dependable predictions for deformations and failure of bulk materials, heterogeneous interfaces, and composite materials up to a billion atoms and length scales of a micrometer, supporting the generation of reliable large data sets for machine learning and accelerated multiscale simulations.
The representation of chemical reactions involving bond formation requires additional assumptions, e.g., reaction templates, and can be integrated using the REACTER framework. Hereby, reaction pathways and parameter sets for products need to be provided, or determined if necessary, before the simulation. The derivation of force field parameters remains the same as for IFF, the only difference lies in using a Morse potential for covalent-type bonds. A general on-the-fly parameterization of IFF-R would be desirable to further simplify reactive simulations for the user, however, such a method requires further work due to the complexity of physical principles and chemical analogy involved, case-specific search and curation of reliable reference data, and extensive simulations for validation.
Methods
Molecular models were built using the Materials Studio 7.0 Graphical User Interface100, and molecular dynamics simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, version Sep2021)57, as well as using the Discover program. OVITO101 2.8, VMD102, the Multiwfn103 software, and Materials Studio were employed to visualize MD trajectories and prepare images of the simulation models. REACTER map files were created manually, and the AutoMapper tool developed by Bone et al.104. was used for conversion of REACTER templates to LAMMPS molecule files.
Bond scan simulations
Ab-initio simulations of bond scans were performed using MP2 calculations based on Hartree-Fock theory with added electron correlations and DFT with the B3LYP gradient-corrected hybrid functional using Gaussian 2016. The amine monomers were aligned in the x-y plane of the simulation cell and the bond length was increased in increments of 0.1 Å or less. Thereby, atoms not participating in the bond were constrained in their relative positions. Energy minimization was performed at each bond length with 10-6 eV total energy convergence to construct the zero-point bond energy profile. Molecular mechanics simulations of bond scans with IFF-R and ReaxFF involved energy minimization using the conjugate gradient algorithm in LAMMPS for 100 to 500 steps (up to 10,000 steps for verification) with an energy tolerance of 10-4 and a force tolerance of 10-6 kcal mol-1 Å-1) (Section S1.1. in the Supplementary Information).
Simulations of stress–strain curves
For the simulations of stress-strain properties up to failure, models were initially subjected to energy minimization. Subsequent molecular dynamics simulations with IFF or IFF-R employed a temperature of 298.15 K, setting initial atomic velocities using a Gaussian distribution, and the Nose-Hoover thermostat for temperature control including dampening within 100 timesteps. We employed the NPT ensemble, except for CNTs that include some vacuum in the simulation box where the NVT ensemble is more suitable, and the Nose-Hoover barostat for pressure control at 1 atm. The tensile strain was increased continuously during molecular dynamics simulations in a strain program from zero until after failure. The maximum strain varied depending on the material from 0.4 for CNTs to 6.0 for polymers. The strain rate was between 0.04 and 0.10 per 100 ps, equal to between 0.4 ns-1 and 1.0 ns-1. All components of the stress tensor were measured using the NPT ensemble for polymers and metals105,106,107, and using the NVT ensemble for CNTs. In the case of CNTs, the engineering stress in the tensile direction was calculated using the cross-sectional area of the CNTs in relation to the total cross-sectional area of the box that included some vacuum (Section S1.2 in the Supplementary Information).
The summation of Coulomb interactions used the PPPM K-space solver, including a cutoff for the direct summation at 8.0 Å and relative energy tolerance of 10-4 for the long range, corresponding to high accuracy. The summation of Lennard-Jones interactions, which represent van-der-Waals interactions, was carried out with a spherical cutoff at 12.0 Å which is standard for IFF. The time step was typically 1 fs and lower at 0.5 fs for models with hydrogen or virtual π-electrons. Tensile simulations were used to strain periodic models along the axial direction of the carbon nanotube or polymer chains, and perpendicular to the [111] facet for the iron crystal at a strain rate of 10-6 every 10 timesteps with pressure control at 1 atm in the lateral directions. The simulations were run until 100% strain or material failure. Stress was monitored relative to the initial cross-sectional area of the models. The tensile modulus and the tensile strength were determined from the computed stress-strain curves, e.g., modulus as a ratio of stress over strain at low deformation (strain <0.01). Standard deviations for the tensile modulus and strength were obtained from three independent simulated models with different initial velocities.
Additional simulations at the DFT level for PAN are described in Section S1.3 in the Supplementary Information.
Evaluation of structural and surface properties
The evaluation of equilibrium crystal structures and liquid densities also utilized molecular dynamics simulations in the NPT ensemble. The computation of surface energies (using equilibrium lateral dimensions of the surface from simulations in the NPT ensemble), simulations of stress-strain curves of CNTs, and of the average energy of gas molecules as part of the calculation of the vaporization energies utilized the NVT ensemble. The reason to use the NVT ensemble is that these simulation setups include major parts of vacuum or gas phase.
Equilibrium crystal structures and liquid densities were obtained by subjecting 3D periodic bulk structures to 1 ns molecular dynamics at 1 atm pressure and 298 K in the NPT ensemble. The average lattice parameters from this simulation, in case of crystalline materials graphite and iron, were further used to create models with cleaved surfaces by expanding the simulation cell 4.0 nm in the z-direction, creating an implicit surface because of the periodic z-boundary. Hereby, we used the (0001) cleavage plane for graphite and the (111) cleavage plane for iron. The cleaved structure and the bulk structure with average lattice dimensions were subjected to 1 ns of molecular dynamics at 298 K in the NVT ensemble. The average energies for the cleaved and bulk models were used to calculate surface energy, normalized relative to the surface area (Section S1.4 in the Supplementary Information).
Evaluation of vaporization and sublimation enthalpies
Calculations of the molar enthalpy of vaporization and sublimation (ΔH) were performed by obtaining average energies for the crystalline (α-D-glucose) or liquid (propionitrile) structures, respectively, and for the same structures in the gaseous state, using the energy difference (ΔU) and a correction for volume work (ΔH = ΔU + PV). Simulations of the crystalline and liquid models were run for 1 ns at 298 K and 1 atm pressure in the NPT ensemble. Simulations of the gases were carried out for 1 ns at 298 K in the NVT ensemble. The average energy per molecule was used to calculate the enthalpy of sublimation or vaporization as a difference between the gaseous and the solid or liquid state, respectively. Standard deviations for computed surface energies, enthalpies of vaporization, and enthalpies of sublimation were obtained from the block averages of the total energies in equilibrium (Section S1.5 in the Supplementary Information).
Examples for deriving IFF-R parameters and improving nonreactive base parameters (e.g., PCFF) are given in Sections S1.6 and S1.7 in Supplementary Information.
Simulation of bond-forming reactions
We used the REACTER framework to facilitate bond formation47. The required reaction templates were created in Materials Studio and the molecules (PAN, amine/epoxy) were allowed to react during MD simulation at 1 atm pressure over a period of 150 ps using a reaction distance of 6.0 Å, which included equilibration and a reaction phase of 50 ps. Two thermostats were used during the reaction phase: (1) a Nose-Hoover thermostat at 600 K and 100 fs dampening for all atoms not participating in the reaction and (2) a velocity rescale thermostat at 200 K and 100 fs dampening for all atoms participating in the reaction. This setup uses time-temperature equivalence to accelerate the reaction time from milliseconds in real time to less than nanoseconds in the simulation. Atoms not participating in the reaction received enough kinetic energy at 600 K to move more rapidly and explore the range of potential reaction constellations in a short timeframe of 50 ps. During the reaction, atoms are relatively near to each other and may experience large forces upon connection, as well as strongly repulsive interactions when nonbonded at short distance or overlapping. Limiting the temperature to 200 K in the reacting regions mitigates such interactions and prevents failure of the simulations. Lower reaction distances, e.g., 4 Å to 5 Å, along with longer simulation times can also be explored and reduce the need for two dissimilar thermostats (Section S1.8 in the Supplementary Information).
The concept of IFF-R was originally described on the Heinz group website in January 2020 along with a tutorial108 and an ArXiv preprint109. Supplementary Data 1 include sample force field files, simulation input scripts, 3D molecular models to reproduce simulations, as well as a tutorial describing the workflow and examples for IFF-R. Supplementary Software contains code for automated Morse bond conversion and a new LAMMPS option for the shifted Morse potential. The Supplementary Information includes additional figures and discussion on the benefit of the shifted Morse potential versus the standard Morse potential, stress-strain curves and total energies up to failure for different bond cutoffs, individual energy contributions as a function of strain, reaction pathways after bond dissociation, the relationship between resonance structures and bond strength for selected examples, further mechanical property analysis for select materials, the comparison of the computational speed of IFF-R for large systems >100,000 atoms relative to ReaxFF, as well as details of the IFF and IFF-R parameterization of organic molecules. Supplementary tables contain example parameters for the Morse potential and further comparisons of computed properties using IFF-R and ReaxFF to experimental data.
Data availability
All data to reproduce the results are available in the manuscript, Supplementary Information, Supplementary Data 1, and can also be obtained upon request from the corresponding author. Specifically, Supplementary Data 1 contains 3D atomic models, force field files, and simulation scripts to reproduce the results. Supplementary Software contains code to deploy the newly developed IFF-R tools, including examples. The same data and code are also shared on Zenodo110. Source data are provided with this paper.
Code availability
All code is available in the Supplementary Software, including the GUI and the custom LAMMPS option for the shifted Morse potential. The same data and code are also shared on Zenodo110.
References
Liu, Z. Q., Zhang, Z. F. & Ritchie, R. O. Structural orientation and anisotropy in biological materials: Functional designs and mechanics. Adv. Funct. Mater. 30, 1908121 (2020).
Ling, S. J., Kaplan, D. L. & Buehler, M. J. Nanofibrils in nature and materials engineering. Nat. Rev. Mater. 3, 1–15 (2018).
Chang, H. B., Luo, J., Gulgunje, P. V. & Kumar, S. Structural and functional fibers. Annu. Rev. Mater. Res. 47, 331–359 (2017).
Nepal, D. et al. Hierarchically structured bioinspired nanocomposites. Nat. Mater. 22, 18–35 (2023).
Heinz, H., Lin, T. J., Mishra, R. K. & Emami, F. S. Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The interface force field. Langmuir 29, 1754–1765 (2013).
Lindorff-Larsen, K., Piana, S., Dror, R. O. & Shaw, D. E. How fast-folding proteins fold. Science 334, 517–520 (2011).
Heinz, H. & Ramezani-Dakhel, H. Simulations of inorganic–bioorganic interfaces to discover new materials: Insights, comparisons to experiment, challenges, and opportunities. Chem. Soc. Rev. 45, 412–448 (2016).
Pisani, W. A. et al. Multiscale modeling of PEEK using reactive molecular dynamics modeling and micromechanics. Polymer 163, 96–105 (2019).
Ruan, L. Y. et al. A rational biomimetic approach to structure defect generation in colloidal nanocrystals. ACS Nano 8, 6934–6944 (2014).
Kanhaiya, K., Kim, S., Im, W. & Heinz, H. Accurate simulation of surfaces and interfaces of ten FCC metals and steel using lennard–jones potentials. npj Comput. Mater. 7, 17 (2021).
Jamil, T., Javadi, A. & Heinz, H. Mechanism of molecular interaction of acrylate-polyethylene glycol acrylate copolymers with calcium silicate hydrate surfaces. Green. Chem. 22, 1577–1593 (2020).
Gissinger, J. R., Pramanik, C., Newcomb, B., Kumar, S. & Heinz, H. Nanoscale structure-property relationships of polyacrylonitrile/CNT composites as a function of polymer crystallinity and CNT diameter. ACS Appl. Mater. Interfaces 10, 1017–1027 (2018).
Hoff, S. E., Di Silvio, D., Ziolo, R. F., Moya, S. E. & Heinz, H. Patterning of self-assembled monolayers of amphiphilic multisegment ligands on nanoparticles and design parameters for protein interactions. ACS Nano 16, 8766–8783 (2022).
Heinz, H., Vaia, R. A., Koerner, H. & Farmer, B. L. Photoisomerization of azobenzene grafted to layered silicates: Simulation and experimental challenges. Chem. Mater. 20, 6444–6456 (2008).
Ramezani-Dakhel, H., Mirau, P. A., Naik, R. R., Knecht, M. R. & Heinz, H. Stability, surface features, and atom leaching of palladium nanoparticles: Toward prediction of catalytic functionality. Phys. Chem. Chem. Phys. 15, 5488–5492 (2013).
Coppage, R. et al. Exploiting localized surface binding effects to enhance the catalytic reactivity of peptide-capped nanoparticles. J. Am. Chem. Soc. 135, 11048–11054 (2013).
Kanhaiya, K. et al. Accurate force fields for atomistic simulations of oxides, hydroxides, and organic hybrid materials up to the micrometer scale. J. Chem. Theor. Comput. 19, 8293–8322 (2023).
Mishra, R. K., Kanhaiya, K., Winetrout, J. J., Flatt, R. J. & Heinz, H. Force field for calcium sulfate minerals to predict structural, hydration, and interfacial properties. Cem. Concr. Res. 139, 106262 (2021).
Sun, H., Mumby, S. J., Maple, J. R. & Hagler, A. T. An ab-initio cff93 all-atom force-field for polycarbonates. J. Am. Chem. Soc. 116, 2978–2987 (1994).
Weiner, P. K. & Kollman, P. A. Amber - assisted model-building with energy refinement - a general program for modeling molecules and their interactions. J. Comput. Chem. 2, 287–303 (1981).
Dauber-Osguthorpe, P. et al. Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system. Proteins: Struct., Funct., Genet. 4, 31–47 (1988).
Brooks, B. R. et al. Charmm - a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4, 187–217 (1983).
Jorgensen, W. L., Maxwell, D. S. & TiradoRives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225–11236 (1996).
Sun, H. Compass: An ab initio force-field optimized for condensed-phase applications - overview with details on alkane and benzene compounds. J. Phys. Chem. B 102, 7338–7364 (1998).
Mayo, S. L., Olafson, B. D. & Goddard, W. A. Dreiding - a generic force-field for molecular simulations. J. Phys. Chem. 94, 8897–8909 (1990).
Gissinger, J. R., Jensen, B. D. & Wise, K. E. Modeling chemical reactions in classical molecular dynamics simulations. Polymer 128, 211–217 (2017).
Liu, J. et al. Understanding chemical bonding in alloys and the representation in atomistic simulations. J. Phys. Chem. C. 122, 14996–15009 (2018).
Bedford, N. M. et al. Elucidation of peptide-directed palladium surface structure for biologically tunable nanocatalysts. ACS Nano 9, 5082–5092 (2015).
Pustovgar, E. et al. Influence of aluminates on the hydration kinetics of tricalcium silicate. Cem. Concr. Res. 100, 245–262 (2017).
Mishra, R. K., Fernández-Carrasco, L., Flatt, R. J. & Heinz, H. A force field for tricalcium aluminate to characterize surface properties, initial hydration, and organically modified interfaces in atomic resolution. Dalton Trans. 43, 10602–10616 (2014).
Mishra, R. K., Flatt, R. J. & Heinz, H. Force field for tricalcium silicate and insight into nanoscale properties: Cleavage, initial hydration, and adsorption of organic molecules. J. Phys. Chem. C. 117, 10417–10432 (2013).
Heinz, H. & Suter, U. W. Atomic charges for classical simulations of polar systems. J. Phys. Chem. B 108, 18341–18352 (2004).
Zhu, C., Hoff, S. E., Hémadi, M. & Heinz, H. Accurate and ultrafast simulation of molecular recognition and assembly on metal surfaces in four dimensions. ACS Nano 17, 9938–9952 (2023).
Sauer, J. Ab initio calculations for molecule-surface interactions with chemical accuracy. Acc. Chem. Res. 52, 3502–3510 (2019).
Ruiz, V. G., Liu, W. & Tkatchenko, A. Density-functional theory with screened van der waals interactions applied to atomic and molecular adsorbates on close-packed and non-close-packed surfaces. Phys. Rev. B 93, 035118 (2016).
Goerigk, L. & Grimme, S. A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions. Phys. Chem. Chem. Phys. 13, 6670–6688 (2011).
Dharmawardhana, C. C. et al. Reliable computational design of biological-inorganic materials to the large nanometer scale using interface-ff. Mol. Simul. 43, 1394–1405 (2017).
Pramanik, C. et al. Molecular engineering of interphases in polymer/carbon nanotube composites to reach the limits of mechanical performance. Compos. Sci. Technol. 166, 86–94 (2018).
Sun, H. Ab-initio calculations and force-field development for computer-simulation of polysilanes. Macromolecules 28, 701–712 (1995).
Brenner, D. W. et al. A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons. J. Phys. Condens. Mat. 14, 783–802 (2002).
Pauling, L. The nature of the chemical-bond - 1992. J. Chem. Educ. 69, 519–521 (1992).
Senftle, T. P. et al. The reaxff reactive force-field: Development, applications and future directions. npj Comput. Mater. 2, 15011 (2016).
Chenoweth, K., van Duin, A. C. T. & Goddard, W. A. Reaxff reactive force field for molecular dynamics simulations of hydrocarbon oxidation. J. Phys. Chem. A 112, 1040–1053 (2008).
Moerman, E., Furman, D. & Wales, D. J. Systematic evaluation of reaxff reactive force fields for biochemical applications. J. Chem. Theor. Comput. 17, 497–514 (2021).
Morse, P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev. 34, 57 (1929).
Costa Filho, R. N., Alencar, G., Skagerstam, B. S. & Andrade, J. S. Morse potential derived from first principles. EPL 101, 10009 (2013).
Gissinger, J. R., Jensen, B. D. & Wise, K. E. Reacter: A heuristic method for reactive molecular dynamics. Macromolecules 53, 9953–9961 (2020).
Tersoff, J. Empirical interatomic potential for silicon with improved elastic properties. Phys. Rev. B 38, 9902–9905 (1988).
van Duin, A. C. T., Dasgupta, S., Lorant, F. & Goddard, W. A. Reaxff: A reactive force field for hydrocarbons. J. Phys. Chem. A 105, 9396–9409 (2001).
Stuart, S. J., Tutein, A. B. & Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 112, 6472–6486 (2000).
Lide, D. R. CRC handbook of chemistry and physics 96th edn (CRC Press, Boca Raton, FL, 2015).
Blanksby, S. J. & Ellison, G. B. Bond dissociation energies of organic molecules. Acc. Chem. Res. 36, 255–263 (2003).
Darwent, B. D. In National Standard Reference Data Series, National Bureau of Standards (U. S) Vol. 31 (ed. U.S. Department of Commerce) 9–47 (Washington, D.C, 1970).
NIST chemistry webbook. (US Department of Commerce, https://doi.org/10.18434/T4D303, 2018).
St John, P. C., Guan, Y. F., Kim, Y., Kim, S. & Paton, R. S. Prediction of organic homolytic bond dissociation enthalpies at near chemical accuracy with sub-second computational cost. Nat. Commun. 11, 2328 (2020).
Heinz, H., Koerner, H., Anderson, K. L., Vaia, R. A. & Farmer, B. L. Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite. Chem. Mater. 17, 5658–5669 (2005).
Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19 (1995).
Wang, S., Liang, R., Wang, B. & Zhang, C. Covalent addition of diethyltoluenediamines onto carbon nanotubes for composite application. Polym. Compos. 30, 1050–1057 (2009).
Xie, H. et al. Thermal characterization of carbon-nanofiber-reinforced tetraglycidyl-4,4′-diaminodiphenylmethane/4,4′-diaminodiphenylsulfone epoxy composites. J. Appl. Polym. Sci. 100, 295–298 (2006).
Singh, S. K. et al. Thermal properties of fluorinated graphene. Phys. Rev. B 87, 104114 (2013).
Geada, I. L., Ramezani-Dakhel, H., Jamil, T., Sulpizi, M. & Heinz, H. Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential. Nat. Comm. 9, 716 (2018).
Bryenton, K. R., Adeleke, A. A., Dale, S. G. & Johnson, E. R. Delocalization error: The greatest outstanding challenge in density-functional theory. WIREs Computational Mol. Sci. 13, e1631 (2023).
Zavitsas, A. A. The relation between bond lengths and dissociation energies of carbon-carbon bonds. J. Phys. Chem. A 107, 897–898 (2003).
Smith, M. B. March’s advanced organic chemistry: Reactions, mechanism, and structure. 7th edn (Wiley, Hoboken, NJ, 2013).
Moon, R. J., Martini, A., Nairn, J., Simonsen, J. & Youngblood, J. Cellulose nanomaterials review: Structure, properties and nanocomposites. Chem. Soc. Rev. 40, 3941–3994 (2011).
Brenner, S. S. Tensile strength of whiskers. J. Appl. Phys. 27, 1484–1491 (1956).
Bai, Y. X. et al. Carbon nanotube bundles with tensile strength over 80 GPa. Nat. Nanotechnol. 13, 589–595 (2018).
Yu, M. F., Files, B. S., Arepalli, S. & Ruoff, R. S. Tensile loading of ropes of single wall carbon nanotubes and their mechanical properties. Phys. Rev. Lett. 84, 5552–5555 (2000).
Zhao, Q. Z., Nardelli, M. B. & Bernholc, J. Ultimate strength of carbon nanotubes: A theoretical study. Phys. Rev. B 65, 144105 (2002).
Benito, J. A., Manero, J. M., Jorba, J. & Roca, A. Change of young’s modulus of cold-deformed pure iron in a tensile test. Metall. Trans. A 36A, 3317–3324 (2005).
Clatterbuck, D. M., Chrzan, D. C. & Morris, J. W. The inherent tensile strength of iron. Philos. Mag. Lett. 82, 141–147 (2002).
Diddens, I., Murphy, B., Krisch, M. & Muller, M. Anisotropic elastic properties of cellulose measured using inelastic X-ray scattering. Macromolecules 41, 9755–9759 (2008).
Dri, F. L., Hector, L. G., Moon, R. J. & Zavattieri, P. D. Anisotropy of the elastic properties of crystalline cellulose i-beta from first principles density functional theory with van der waals interactions. Cellulose 20, 2703–2718 (2013).
Damirchi, B. et al. Reaxff reactive force field study of polymerization of a polymer matrix in a carbon nanotube-composite system. J. Phys. Chem. C. 124, 20488–20497 (2020).
Islam, M. M., Zou, C., van Duin, A. C. T. & Raman, S. Interactions of hydrogen with the iron and iron carbide interfaces: A reaxff molecular dynamics study. Phys. Chem. Chem. Phys. 18, 761–771 (2016).
Liu, L. C. et al. Reaxff-/g: Correction of the ReaxFF reactive force field for London dispersion, with applications to the equations of state for energetic materials. J. Phys. Chem. A 115, 11016–11022 (2011).
Kanski, M. et al. Development of a charge-implicit reaxff potential for hydrocarbon systems. J. Phys. Chem. Lett. 9, 359–363 (2018).
Strachan, A., van Duin, A. C. T., Chakraborty, D., Dasgupta, S. & Goddard, W. A. Shock waves in high-energy materials: The initial chemical events in nitramine RDX. Phys. Rev. Lett. 91, 098301 (2003).
Shin, Y. K., Gai, L. L., Raman, S. & van Duin, A. C. T. Development of a reaxff reactive force field for the Pt-Ni. alloy Catal. J. Phys. Chem. A 120, 8044–8055 (2016).
Aryanpour, M., van Duin, A. C. T. & Kubicki, J. D. Development of a reactive force field for iron-oxyhydroxide systems. J. Phys. Chem. A 114, 6298–6307 (2010).
Clark, S. J. et al. First principles methods using CASTEP. Z. Kristallog 220, 567–570 (2005).
Tadmor, E. B., Elliott, R. S., Sethna, J. P., Miller, R. E. & Becker, C. A. Open knowledgebase of interatomic models (openkim). (https://openkim.org, 2011–2023).
Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. Numerical-integration of cartesian equations of motion of a system with constraints - molecular-dynamics of n-alkanes. J. Comput. Phys. 23, 327–341 (1977).
Loeffler, H. H. & Winn, M. D. Large biomolecular simulation on HPC platforms III. Amber, CHARMM, GROMACS, LAMMPS and NAMD. (Technical Report, STFC Daresbury Laboratory, 2013).
Loeffler, H. H. & Winn, M. D. Large biomolecular simulation on HPC platforms ii. Dl poly, gromacs, LAMMPS and NAMD (Technical Report, STFC Daresbury Laboratory, 2010).
Li, X. X., Mo, Z., Liu, J. & Guo, L. Revealing chemical reactions of coal pyrolysis with GPU-enabled reaxff molecular dynamics and cheminformatics analysis. Mol. Simula. 41, 13–27 (2015).
Thompson, A. P. et al. LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 271, 108171 (2022).
Han, Y. et al. Development, applications and challenges of reaxff reactive force field in molecular simulations. Front. Chem. Sci. Eng. 10, 16–38 (2016).
Mcdonald, T. R. R. & Beevers, C. A. The crystal and molecular structure of alpha-glucose. Acta Crystallogr 5, 654–659 (1952).
Pramanik, C., Gissinger, J. R., Kumar, S. & Heinz, H. Carbon nanotube dispersion in solvents and polymer solutions: Mechanisms, assembly, and preferences. ACS Nano 11, 12805–12816 (2017).
Akkineni, S. et al. Amyloid-like amelogenin nanoribbons template mineralization via a low-energy interface of ion binding sites. Proc. Natl Acad. Sci. Usa. 119, e2106965119 (2022).
Odegard, G. M. et al. Molecular dynamics modeling of epoxy resins using the reactive interface force field. Macromolecules 54, 9815–9824 (2021).
Gaikwad, P. S. et al. Understanding the origin of the low cure shrinkage of polybenzoxazine resin by computational simulation. Acs Appl. Polym. Mater. 3, 6407–6415 (2021).
Odegard, G. M. et al. Accurate predictions of thermoset resin glass transition temperatures from all-atom molecular dynamics simulation. Soft Matter 18, 7550–7558 (2022).
Pisani, W. A., Newman, J. K. & Shukla, M. K. Multiscale modeling of polyamide 6 using molecular dynamics and micromechanics. Ind. Eng. Chem. Res 60, 13604–13613 (2021).
Barr, S. A. et al. Bond breaking in epoxy systems: A combined qm/mm approach. J. Chem. Phys. 144 (2016).
Check, C. E. & Gilbert, T. M. Progressive systematic underestimation of reaction energies by the B3LYP model as the number of C-C bonds increases: Why organic chemists should use multiple DFT models for calculations involving polycarbon hydrocarbons. J. Org. Chem. 70, 9828–9834 (2005).
Wodrich, M. D., Corminboeuf, C., Schreiner, P. R., Fokin, A. A. & Schleyer, P. V. How accurate are DFT treatments of organic energies? Org. Lett. 9, 1851–1854 (2007).
Huheey, J. E., Keiter, E. A. & Keiter, R. L. Inorganic chemistry - principles of structure and reactivity. 4th edn (Harper Collins College Publishers, New York, 1993).
Materials studio 2019 program suite and user guide. (Biovia/Dassault Systemes, Cambridge, UK, 2019).
Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool. Modell. Simul. Mater. Sci. Eng. 18, 015012 (2010).
Humphrey, W., Dalke, A. & Schulten, K. Vmd: Visual molecular dynamics. J. Mol. Graph. Modell. 14, 33–38 (1996).
Lu, T. & Chen, F. W. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 33, 580–592 (2012).
Bone, M. A., Howlin, B. J., Hamerton, I. & Macquart, T. Automapper: A Python tool for accelerating the polymer bonding workflow in LAMMPS. Comp. Mater. Sci. 205, 111204 (2022).
Zartman, G. D., Liu, H., Akdim, B., Pachter, R. & Heinz, H. Nanoscale tensile, shear, and failure properties of layered silicates as a function of cation density and stress. J. Phys. Chem. C. 114, 1763–1772 (2010).
Wu, X. W., Moon, R. J. & Martini, A. Tensile strength of iβ crystalline cellulose predicted by molecular dynamics simulation. Cellulose 21, 2233–2245 (2014).
Pan, Z. L., Li, Y. L. & Wei, Q. Tensile properties of nanocrystalline tantalum from molecular dynamics simulations. Acta Mater. 56, 3470–3480 (2008).
Interface Force Field (IFF) and a surface model database. https://bionanostructures.com/interface-md/ (2013−2023).
Winetrout, J. J. et al. Implementing reactivity in molecular dynamics simulations with the interface force field (IFF-R) and other harmonic force fields. Preprint at https://arxiv.org/abs/2107.14418 (2021).
Winetrout, J. J. et al. Implementing reactivity in molecular dynamics simulations with harmonic force fields: Associated software to use the Reactive Interface Force Field (IFF-R). Zenodo, https://doi.org/10.5281/zenodo.12518681 (2024).
Iwamoto, S., Kai, W. H., Isogai, A. & Iwata, T. Elastic modulus of single cellulose microfibrils from tunicate measured by atomic force microscopy. Biomacromolecules 10, 2571–2576 (2009).
Lamb, S. & Bringas, J. E. Casti handbook of stainless steels & nickel alloys. 2nd edn (CASTI Pub., Edmonton, 2002).
Peckner, D. & Bernstein, I. M. Handbook of stainless steels. (McGraw-Hill, New York, 1977).
Trucano, P. & Chen, R. Structure of graphite by neutron-diffraction. Nature 258, 136–137 (1975).
Liu, Z. et al. Interlayer binding energy of graphite: A mesoscopic determination from deformation. Phys. Rev. B 85, 205418 (2012).
Zacharia, R., Ulbricht, H. & Hertel, T. Interlayer cohesive energy of graphite from thermal desorption of polyaromatic hydrocarbons. Phys. Rev. B 69, 155406 (2004).
Oja, V. & Suuberg, E. M. Vapor pressures and enthalpies of sublimation of d-glucose, d-xylose, cellobiose, and levoglucosan. J. Chem. Eng. Data 44, 26–29 (1999).
Basinski, Z. S., Humerothery, W. & Sutton, A. L. The lattice expansion of iron. Proc. R. Soc. Lond., Ser. A 229, 459–467 (1955).
Nishiyama, Y., Langan, P. & Chanzy, H. Crystal structure and hydrogen-bonding system in cellulose 1 beta from synchrotron X-ray and neutron fiber diffraction. J. Am. Chem. Soc. 124, 9074–9082 (2002).
Acknowledgements
This work was supported by the NASA Space Technology Research Institute (STRI) for Ultra-Strong Composites by Computational Design (NNX17AJ32G), the National Science Foundation (CMMI 1940335, OAC 1931587, DMREF 2323546), BASF SE (Ludwigshafen, Germany), and the University of Colorado at Boulder. We acknowledge the allocation of computational resources at the Summit supercomputer, a joint effort of the University of Colorado Boulder and Colorado State University, which is supported by the National Science Foundation (ACI−1532235 and ACI−1532236), and resources at the Argonne Leadership Computing Facility, a DoE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Publication of this article was partially funded by the University of Colorado Boulder Libraries Open Access Fund.
Author information
Authors and Affiliations
Contributions
J.J.W. developed the workflow, carried out nonreactive and reactive MD simulations with IFF, IFF-R and ReaxFF, analyzed the data, and wrote the manuscript. K.K. contributed to the conception and workflow of the study, simulations, and proofreading. J.K. developed workflows, coded and tested the graphical user interface to assign Morse bonds in IFF-R, carried out simulations, and wrote the manuscript. P.J.i.V. coded the custom LAMMPS option for the IFF-R Morse potential, tested calculations, and proofread the manuscript. G.S. and R.P. carried out MP2 and DFT calculations of bond scans. B.D. and A.v.D. carried out reactive simulations with ReaxFF and wrote the manuscript. G.M.O. advised J.K., wrote and revised the manuscript. H.H. conceived the study, developed the workflow, coordinated the effort, analyzed the data, wrote, revised, and proofread the manuscript and supporting files.
Corresponding author
Ethics declarations
Competing interests
The authors declare no conflicts of interest.
Peer review
Peer review information
Nature Communications thanks Jacob F. N. Dethan and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Source data
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Winetrout, J.J., Kanhaiya, K., Kemppainen, J. et al. Implementing reactivity in molecular dynamics simulations with harmonic force fields. Nat Commun 15, 7945 (2024). https://doi.org/10.1038/s41467-024-50793-0
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467-024-50793-0