
==== Front
Nat Commun
Nat Commun
Nature Communications
2041-1723
Nature Publishing Group UK London

39261455
50793
10.1038/s41467-024-50793-0
Article
Implementing reactivity in molecular dynamics simulations with harmonic force fields
Winetrout Jordan J. 12
http://orcid.org/0000-0002-3622-2655
Kanhaiya Krishan 123
Kemppainen Joshua 4
in ‘t Veld Pieter J. 5
Sachdeva Geeta 6
Pandey Ravindra 6
Damirchi Behzad 7
http://orcid.org/0000-0002-3478-4945
van Duin Adri 7
http://orcid.org/0000-0001-7577-6565
Odegard Gregory M. 4
http://orcid.org/0000-0002-6776-7404
Heinz Hendrik hendrik.heinz@colorado.edu

12
1 https://ror.org/02ttsq026 grid.266190.a 0000 0000 9621 4564 Department of Chemical and Biological Engineering, University of Colorado at Boulder, Boulder, CO USA
2 https://ror.org/02ttsq026 grid.266190.a 0000 0000 9621 4564 Materials Science and Engineering Program, University of Colorado at Boulder, Boulder, CO USA
3 https://ror.org/04tsk2644 grid.5570.7 0000 0004 0490 981X Insitute of Physics, Ruhr University Bochum, Universitätstrasse 150, Bochum, Germany
4 https://ror.org/0036rpn28 grid.259979.9 0000 0001 0663 5937 Department of Mechanical Engineering - Engineering Mechanics, Michigan Technological University, Houghton, MI USA
5 grid.3319.8 0000 0001 1551 0781 BASF SE, Molecular Modeling & Drug Discovery, Ludwigshafen, Germany
6 https://ror.org/0036rpn28 grid.259979.9 0000 0001 0663 5937 Department of Physics, Michigan Technological University, Houghton, MI USA
7 https://ror.org/04p491231 grid.29857.31 0000 0001 2097 4281 Department of Mechanical Engineering, Pennsylvania State University, University Park, PA USA
11 9 2024
11 9 2024
2024
15 794516 3 2023
17 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ 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/.
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.

Molecular dynamics is a common tool to study microscopic physicochemical systems, however, it is limited by the inhability to form and break chemical bonds. Here the authors present a method to modify traditional force-fields implementing bond dissociation and bond forming.

Subject terms

Structural properties
Reaction kinetics and dynamics
Computational biophysics
Structural properties
https://doi.org/10.13039/100000001 National Science Foundation (NSF) https://doi.org/10.13039/100000104 National Aeronautics and Space Administration (NASA) issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

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–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–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–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–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–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).Fig. 1 Concept of IFF-R parameterization.

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–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–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–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.Fig. 2 Examples of bond scans and experimental bond dissociation energies for the parameterization of the reactive INTERFACE Force Field, IFF-R.

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–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–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–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–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–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–73. The computed mechanical properties from IFF-R agree with most reference data within 10%.Fig. 3 Stress-strain curves and snapshots of the materials under uniaxial tensile strain using the reactive INTERFACE force field (IFF-R), the reactive polymer consistent force field (PCFF-R), and ReaxFF parameter sets.

When available, the experimentally determined tensile strengths, or high-level theory estimates, are shown as horizontal lines a, b, d, f, h (refs. 65–69). a, b Single wall carbon nanotubes (SWCNT) of 1 nm diameter and comparison to simulations with ReaxFF (refs. 43,74,76–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.

Table 1 Young’s Modulus and the tensile strength for SWCNTs, crystalline syndiotactic PAN, Iβ cellulose, and γ-iron (representative of steels) according to experimental measurements and molecular dynamics simulations with IFF-R, PCFF and ReaxFF (refs. 74–76)

	SWCNT	PAN	Cellulose Iβ (PCFF-R)	γ-Iron, [111] directionc	
Modulus (GPa)	Strength (GPa)	Modulus (GPa)	Strength (GPa)	Modulus (GPa)	Strength (GPa)	Modulus (GPa)	Strength (GPa)	
Experiment	1007 ± 118a	121.6 ± 22a	172 ± 10 (DFT)	N/A	150 ± 5.0111	4–1065	210 ± 1070	13.1 ± 1.566	
IFF-R (or PCFF-R)	1001 ± 15	115 ± 1.5	168 ± 9.0	18.5 ± 2.0	120 ± 1.0	12.7 ± 1.5	237 ± 20	14.5 ± 1.1	
IFF-R dev to expt (%)	-1 ± 12	-5 ± 19	-2 ± 6	NA	-20 ± 3	+81 ± 43	+13 ± 5	+11 ± 11	
ReaxFFb	752 ± 8.0	173 ± 2.0	219 ± 6.0	29 ± 2.0	131 ± 4.0	14.1 ± 2.0	4627 ± 12.0	35 ± 5.0	
ReaxFF dev to expt (%)	-25 ± 12	+42 ± 19	+27 ± 6	NA	-13 ± 3	+101 ± 43	+2100 ± 5	+169 ± 11	
a From tensile testing of a 2.1 nm CNT using a micro/nanoscale material testing system (ref. 67).

b Using the ReaxFF parameter set from Liu et al. (ref. 76) for SWCNTs, PAN, cellulose and from Islam et al. (ref. 75) for iron. Values may change for a different ReaxFF parameter set.

c Experimental data for γ-iron are similar to the mechanical properties of steel112,113. The [111] crystallographic direction was used, i.e., vertical to the (111) plane.

The mechanical properties are reported in the alignment direction of CNTs and polymers, and the models are representative of different materials classes, including ceramics, polymers, biopolymers, and metals

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–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–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.Fig. 4 Comparison of computing time (CPU time) needed for the simulation of stress-strain curves of different materials using the reactive INTERFACE force field (IFF-R), or reactive polymer consistent force field (PCFF-R) in case of cellulose, and ReaxFF.

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.Table 2 Preservation of key properties upon changing the force field from IFF to IFF-R

Material		Lattice parameters	Mass density (g cm-3)	Surface energy (S) (mJ/m2) or vaporization enthalpy (V) (kJ mol-1)	
a, b, c (nm)	α, β, γ (°)	
Graphite	Expt114–116	0.246, 0.426, 0.680	90, 90, 90	2.26	186 ± 2.0 (S)	
IFF	0.247, 0.428, 0.675	90, 90, 90	2.24	188 ± 2.0	
IFF-R	0.247, 0.428, 0.676	90, 90, 90	2.23	188 ± 1.9	
Propionitrile	Expt51,54	N/A (liquid at 298 K)	0.782	36.1 ± 1.6 (V)	
IFF	0.786	37.7 ± 0.15	
IFF-R	0.783	37.0 ± 0.35	
α-D-glucose	Expt89,117	1.036, 1.483, 0.493	90, 90, 90	1.56	181 ± 5.0 (V)	
PCFF	0.912, 1.268, 0.670	90, 90, 90	1.54	176 ± 0.26	
PCFF-R	0.912, 1.267, 0.670	90, 90, 90	1.55	176 ± 0.33	
γ-irona	Expt118	3.566, 3.566, 3.566	90, 90, 90	8.18	2250 ± 100 (S)	
IFF10	3.587, 3.587, 3.587	90, 90, 90	8.03	2270 ± 7.0	
IFF-R	3.587, 3.587, 3.587	90, 90, 90	8.03	2270 ± 7.0	
a Experimental lattice parameters and surface energies were extrapolated for 298 K from experimental data in ref. 118, see also ref. 10.

Lattice parameters, mass densities, and surface energies (or vaporization enthalpies) from experimental measurements are compared to results from molecular dynamics simulations with IFF and IFF-R under standard conditions. For α-D-glucose, PCFF/PCFF-R parameters were utilized.

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).Table 3 Comparison of surface energies and vaporization energies of representative materials (ceramic, metal, organic) from experimental measurements with computed values using IFF-R and ReaxFF. Examples include the surface energies of graphite, iron), the vaporization energy of a molecular liquid (propionitrile) and the sublimation energy of a molecular solid (α-D-glucose)

Material	Experiment	IFF-R	Dev. of IFF-R (%)	ReaxFFa	Dev. of ReaxFF (%)	
Surface energy (mJ m-2)	
Graphite	186 ± 2.0115,116	188 ± 1.9	+1.1	242 ± 2.0	+30.1	
γ-Iron [111]	2250 ± 10098	2270 ± 7.0	+0.9	2376 ± 8.0	+5.3	
Vaporization energy (kJ mol-1)	
Propionitrile	36.1 ± 1.654	37.0 ± 0.35	+2.5	47 ± 0.9	+30	
α-D-glucose	181.7 ± 5.0b, 117	176 ± 0.33	-3.1	169 ± 1.6	-7.0	
The reliability to compute these energy-related properties scales with the reliability to predict adsorption, stability of mixtures, interfaces, and reactions.

a ReaxFF parameters from Liu et al. for graphite, propionitrile, and α-D-glucose (ref. 76). Alternative Damirchi parameters (ref. 74) are likely less suited. ReaxFF parameters by Shin et al. were used to calculate the surface energy of iron (ref. 79). The alternative Aryanpour et al. parameter set for Fe (ref. 80) would incur a deviation of 200% from experimental reference data and was not chosen here. b The experimental value of the vaporization enthalpy was extrapolated for 298 K from an exponential fit of the integrated form of the Clausius-Clapeyron equation (see ref. 117).

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–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).Fig. 5 Simulation of micro-reversibility of bond breaking interactions.

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.

Fig. 6 Molecular dynamics simulation of cross-linking an amine-epoxy network and a carbon nanotube (CNT)/epoxy polymer matrix composite, followed by the simulation of tensile strain ɛ and failure mechanisms.

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

Supplementary information

Supplementary Information

Peer Review File

Description of Additional Supplementary Files

Supplementary Dataset 1

Supplementary Software

Source data

Source Data

Supplementary information

The online version contains supplementary material available at 10.1038/s41467-024-50793-0.

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

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.

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.

Competing interests

The authors declare no conflicts of interest.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Liu ZQ Zhang ZF Ritchie RO Structural orientation and anisotropy in biological materials: Functional designs and mechanics Adv. Funct. Mater. 2020 30 1908121 10.1002/adfm.201908121
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).10.1002/adfm.201908121
2. Ling SJ Kaplan DL Buehler MJ Nanofibrils in nature and materials engineering Nat. Rev. Mater. 2018 3 1 15 10.1038/natrevmats.2018.16
Ling, S. J., Kaplan, D. L. & Buehler, M. J. Nanofibrils in nature and materials engineering. Nat. Rev. Mater. 3, 1–15 (2018).10.1038/natrevmats.2018.16
3. Chang HB Luo J Gulgunje PV Kumar S Structural and functional fibers Annu. Rev. Mater. Res. 2017 47 331 359 10.1146/annurev-matsci-120116-114326
Chang, H. B., Luo, J., Gulgunje, P. V. & Kumar, S. Structural and functional fibers. Annu. Rev. Mater. Res. 47, 331–359 (2017).10.1146/annurev-matsci-120116-114326
4. Nepal D Hierarchically structured bioinspired nanocomposites Nat. Mater. 2023 22 18 35 10.1038/s41563-022-01384-1 36446962
Nepal, D. et al. Hierarchically structured bioinspired nanocomposites. Nat. Mater. 22, 18–35 (2023).36446962 10.1038/s41563-022-01384-1
5. Heinz H Lin TJ Mishra RK Emami FS Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The interface force field Langmuir 2013 29 1754 1765 10.1021/la3038846 23276161
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).23276161 10.1021/la3038846
6. Lindorff-Larsen K Piana S Dror RO Shaw DE How fast-folding proteins fold Science 2011 334 517 520 10.1126/science.1208351 22034434
Lindorff-Larsen, K., Piana, S., Dror, R. O. & Shaw, D. E. How fast-folding proteins fold. Science 334, 517–520 (2011).22034434 10.1126/science.1208351
7. Heinz H Ramezani-Dakhel H Simulations of inorganic–bioorganic interfaces to discover new materials: Insights, comparisons to experiment, challenges, and opportunities Chem. Soc. Rev. 2016 45 412 448 10.1039/C5CS00890E 26750724
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).26750724 10.1039/C5CS00890E
8. Pisani WA Multiscale modeling of PEEK using reactive molecular dynamics modeling and micromechanics Polymer 2019 163 96 105 10.1016/j.polymer.2018.12.052
Pisani, W. A. et al. Multiscale modeling of PEEK using reactive molecular dynamics modeling and micromechanics. Polymer 163, 96–105 (2019).10.1016/j.polymer.2018.12.052
9. Ruan LY A rational biomimetic approach to structure defect generation in colloidal nanocrystals ACS Nano 2014 8 6934 6944 10.1021/nn501704k 24937767
Ruan, L. Y. et al. A rational biomimetic approach to structure defect generation in colloidal nanocrystals. ACS Nano 8, 6934–6944 (2014).24937767 10.1021/nn501704k
10. 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. 2021 7 17 10.1038/s41524-020-00478-1
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).10.1038/s41524-020-00478-1
11. Jamil T Javadi A Heinz H Mechanism of molecular interaction of acrylate-polyethylene glycol acrylate copolymers with calcium silicate hydrate surfaces Green. Chem. 2020 22 1577 1593 10.1039/C9GC03287H
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).10.1039/C9GC03287H
12. Gissinger JR 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 2018 10 1017 1027 10.1021/acsami.7b09739 29231715
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).29231715 10.1021/acsami.7b09739
13. Hoff SE Di Silvio D Ziolo RF Moya SE Heinz H Patterning of self-assembled monolayers of amphiphilic multisegment ligands on nanoparticles and design parameters for protein interactions ACS Nano 2022 16 8766 8783 10.1021/acsnano.1c08695 35603431
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).35603431 10.1021/acsnano.1c08695
14. Heinz H Vaia RA Koerner H Farmer BL Photoisomerization of azobenzene grafted to layered silicates: Simulation and experimental challenges Chem. Mater. 2008 20 6444 6456 10.1021/cm801287d
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).10.1021/cm801287d
15. Ramezani-Dakhel H Mirau PA Naik RR Knecht MR Heinz H Stability, surface features, and atom leaching of palladium nanoparticles: Toward prediction of catalytic functionality Phys. Chem. Chem. Phys. 2013 15 5488 5492 10.1039/c3cp00135k 23474536
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).23474536 10.1039/c3cp00135k
16. Coppage R Exploiting localized surface binding effects to enhance the catalytic reactivity of peptide-capped nanoparticles J. Am. Chem. Soc. 2013 135 11048 11054 10.1021/ja402215t 23865951
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).23865951 10.1021/ja402215t
17. Kanhaiya K Accurate force fields for atomistic simulations of oxides, hydroxides, and organic hybrid materials up to the micrometer scale J. Chem. Theor. Comput. 2023 19 8293 8322 10.1021/acs.jctc.3c00750
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).10.1021/acs.jctc.3c00750
18. Mishra RK Kanhaiya K Winetrout JJ Flatt RJ Heinz H Force field for calcium sulfate minerals to predict structural, hydration, and interfacial properties Cem. Concr. Res. 2021 139 106262 10.1016/j.cemconres.2020.106262
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).10.1016/j.cemconres.2020.106262
19. Sun H Mumby SJ Maple JR Hagler AT An ab-initio cff93 all-atom force-field for polycarbonates J. Am. Chem. Soc. 1994 116 2978 2987 10.1021/ja00086a030
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).10.1021/ja00086a030
20. Weiner PK Kollman PA Amber - assisted model-building with energy refinement - a general program for modeling molecules and their interactions J. Comput. Chem. 1981 2 287 303 10.1002/jcc.540020311
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).10.1002/jcc.540020311
21. Dauber-Osguthorpe P Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system Proteins: Struct., Funct., Genet. 1988 4 31 47 10.1002/prot.340040106 3054871
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).3054871 10.1002/prot.340040106
22. Brooks BR Charmm - a program for macromolecular energy, minimization, and dynamics calculations J. Comput. Chem. 1983 4 187 217 10.1002/jcc.540040211
Brooks, B. R. et al. Charmm - a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4, 187–217 (1983).10.1002/jcc.540040211
23. Jorgensen WL Maxwell DS TiradoRives J Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids J. Am. Chem. Soc. 1996 118 11225 11236 10.1021/ja9621760
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).10.1021/ja9621760
24. 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 1998 102 7338 7364 10.1021/jp980939v
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).10.1021/jp980939v
25. Mayo SL Olafson BD Goddard WA Dreiding - a generic force-field for molecular simulations J. Phys. Chem. 1990 94 8897 8909 10.1021/j100389a010
Mayo, S. L., Olafson, B. D. & Goddard, W. A. Dreiding - a generic force-field for molecular simulations. J. Phys. Chem. 94, 8897–8909 (1990).10.1021/j100389a010
26. Gissinger JR Jensen BD Wise KE Modeling chemical reactions in classical molecular dynamics simulations Polymer 2017 128 211 217 10.1016/j.polymer.2017.09.038 33149370
Gissinger, J. R., Jensen, B. D. & Wise, K. E. Modeling chemical reactions in classical molecular dynamics simulations. Polymer 128, 211–217 (2017).33149370 10.1016/j.polymer.2017.09.038
27. Liu J Understanding chemical bonding in alloys and the representation in atomistic simulations J. Phys. Chem. C. 2018 122 14996 15009 10.1021/acs.jpcc.8b01891
Liu, J. et al. Understanding chemical bonding in alloys and the representation in atomistic simulations. J. Phys. Chem. C. 122, 14996–15009 (2018).10.1021/acs.jpcc.8b01891
28. Bedford NM Elucidation of peptide-directed palladium surface structure for biologically tunable nanocatalysts ACS Nano 2015 9 5082 5092 10.1021/acsnano.5b00168 25905675
Bedford, N. M. et al. Elucidation of peptide-directed palladium surface structure for biologically tunable nanocatalysts. ACS Nano 9, 5082–5092 (2015).25905675 10.1021/acsnano.5b00168
29. Pustovgar E Influence of aluminates on the hydration kinetics of tricalcium silicate Cem. Concr. Res. 2017 100 245 262 10.1016/j.cemconres.2017.06.006
Pustovgar, E. et al. Influence of aluminates on the hydration kinetics of tricalcium silicate. Cem. Concr. Res. 100, 245–262 (2017).10.1016/j.cemconres.2017.06.006
30. Mishra RK Fernández-Carrasco L Flatt RJ Heinz H A force field for tricalcium aluminate to characterize surface properties, initial hydration, and organically modified interfaces in atomic resolution Dalton Trans. 2014 43 10602 10616 10.1039/C4DT00438H 24828263
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).24828263 10.1039/C4DT00438H
31. Mishra RK Flatt RJ Heinz H Force field for tricalcium silicate and insight into nanoscale properties: Cleavage, initial hydration, and adsorption of organic molecules J. Phys. Chem. C. 2013 117 10417 10432 10.1021/jp312815g
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).10.1021/jp312815g
32. Heinz H Suter UW Atomic charges for classical simulations of polar systems J. Phys. Chem. B 2004 108 18341 18352 10.1021/jp048142t
Heinz, H. & Suter, U. W. Atomic charges for classical simulations of polar systems. J. Phys. Chem. B 108, 18341–18352 (2004).10.1021/jp048142t
33. Zhu C Hoff SE Hémadi M Heinz H Accurate and ultrafast simulation of molecular recognition and assembly on metal surfaces in four dimensions ACS Nano 2023 17 9938 9952 10.1021/acsnano.2c10953 37260141
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).37260141 10.1021/acsnano.2c10953
34. Sauer J Ab initio calculations for molecule-surface interactions with chemical accuracy Acc. Chem. Res. 2019 52 3502 3510 10.1021/acs.accounts.9b00506 31765121
Sauer, J. Ab initio calculations for molecule-surface interactions with chemical accuracy. Acc. Chem. Res. 52, 3502–3510 (2019).31765121 10.1021/acs.accounts.9b00506
35. Ruiz VG 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 2016 93 035118 10.1103/PhysRevB.93.035118
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).10.1103/PhysRevB.93.035118
36. Goerigk L Grimme S A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions Phys. Chem. Chem. Phys. 2011 13 6670 6688 10.1039/c0cp02984j 21384027
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).21384027 10.1039/c0cp02984j
37. Dharmawardhana CC Reliable computational design of biological-inorganic materials to the large nanometer scale using interface-ff Mol. Simul. 2017 43 1394 1405 10.1080/08927022.2017.1332414
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).10.1080/08927022.2017.1332414
38. Pramanik C Molecular engineering of interphases in polymer/carbon nanotube composites to reach the limits of mechanical performance Compos. Sci. Technol. 2018 166 86 94 10.1016/j.compscitech.2018.04.013
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).10.1016/j.compscitech.2018.04.013
39. Sun H Ab-initio calculations and force-field development for computer-simulation of polysilanes Macromolecules 1995 28 701 712 10.1021/ma00107a006
Sun, H. Ab-initio calculations and force-field development for computer-simulation of polysilanes. Macromolecules 28, 701–712 (1995).10.1021/ma00107a006
40. Brenner DW A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons J. Phys. Condens. Mat. 2002 14 783 802 10.1088/0953-8984/14/4/312
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).10.1088/0953-8984/14/4/312
41. Pauling L The nature of the chemical-bond - 1992 J. Chem. Educ. 1992 69 519 521 10.1021/ed069p519
Pauling, L. The nature of the chemical-bond - 1992. J. Chem. Educ. 69, 519–521 (1992).10.1021/ed069p519
42. Senftle TP The reaxff reactive force-field: Development, applications and future directions npj Comput. Mater. 2016 2 15011 10.1038/npjcompumats.2015.11
Senftle, T. P. et al. The reaxff reactive force-field: Development, applications and future directions. npj Comput. Mater. 2, 15011 (2016).10.1038/npjcompumats.2015.11
43. Chenoweth K van Duin ACT Goddard WA Reaxff reactive force field for molecular dynamics simulations of hydrocarbon oxidation J. Phys. Chem. A 2008 112 1040 1053 10.1021/jp709896w 18197648
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).18197648 10.1021/jp709896w
44. Moerman E Furman D Wales DJ Systematic evaluation of reaxff reactive force fields for biochemical applications J. Chem. Theor. Comput. 2021 17 497 514 10.1021/acs.jctc.0c01043
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).10.1021/acs.jctc.0c01043
45. Morse PM Diatomic molecules according to the wave mechanics. II. Vibrational levels Phys. Rev. 1929 34 57 10.1103/PhysRev.34.57
Morse, P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev. 34, 57 (1929).10.1103/PhysRev.34.57
46. Costa Filho RN Alencar G Skagerstam BS Andrade JS Morse potential derived from first principles EPL 2013 101 10009 10.1209/0295-5075/101/10009
Costa Filho, R. N., Alencar, G., Skagerstam, B. S. & Andrade, J. S. Morse potential derived from first principles. EPL 101, 10009 (2013).10.1209/0295-5075/101/10009
47. Gissinger JR Jensen BD Wise KE Reacter: A heuristic method for reactive molecular dynamics Macromolecules 2020 53 9953 9961 10.1021/acs.macromol.0c02012
Gissinger, J. R., Jensen, B. D. & Wise, K. E. Reacter: A heuristic method for reactive molecular dynamics. Macromolecules 53, 9953–9961 (2020).10.1021/acs.macromol.0c02012
48. Tersoff J Empirical interatomic potential for silicon with improved elastic properties Phys. Rev. B 1988 38 9902 9905 10.1103/PhysRevB.38.9902
Tersoff, J. Empirical interatomic potential for silicon with improved elastic properties. Phys. Rev. B 38, 9902–9905 (1988).10.1103/PhysRevB.38.9902
49. van Duin ACT Dasgupta S Lorant F Goddard WA Reaxff: A reactive force field for hydrocarbons J. Phys. Chem. A 2001 105 9396 9409 10.1021/jp004368u
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).10.1021/jp004368u
50. Stuart SJ Tutein AB Harrison JA A reactive potential for hydrocarbons with intermolecular interactions J. Chem. Phys. 2000 112 6472 6486 10.1063/1.481208
Stuart, S. J., Tutein, A. B. & Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 112, 6472–6486 (2000).10.1063/1.481208
51. Lide, D. R. CRC handbook of chemistry and physics 96th edn (CRC Press, Boca Raton, FL, 2015).
52. Blanksby SJ Ellison GB Bond dissociation energies of organic molecules Acc. Chem. Res. 2003 36 255 263 10.1021/ar020230d 12693923
Blanksby, S. J. & Ellison, G. B. Bond dissociation energies of organic molecules. Acc. Chem. Res. 36, 255–263 (2003).12693923 10.1021/ar020230d
53. 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).
54. NIST chemistry webbook. (US Department of Commerce, 10.18434/T4D303, 2018).
55. St John PC Guan YF Kim Y Kim S Paton RS Prediction of organic homolytic bond dissociation enthalpies at near chemical accuracy with sub-second computational cost Nat. Commun. 2020 11 2328 10.1038/s41467-020-16201-z 32393773
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).32393773 10.1038/s41467-020-16201-z
56. Heinz H Koerner H Anderson KL Vaia RA Farmer BL Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite Chem. Mater. 2005 17 5658 5669 10.1021/cm0509328
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).10.1021/cm0509328
57. Plimpton S Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 1995 117 1 19 10.1006/jcph.1995.1039
Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19 (1995).10.1006/jcph.1995.1039
58. Wang S Liang R Wang B Zhang C Covalent addition of diethyltoluenediamines onto carbon nanotubes for composite application Polym. Compos. 2009 30 1050 1057 10.1002/pc.20654
Wang, S., Liang, R., Wang, B. & Zhang, C. Covalent addition of diethyltoluenediamines onto carbon nanotubes for composite application. Polym. Compos. 30, 1050–1057 (2009).10.1002/pc.20654
59. Xie H Thermal characterization of carbon-nanofiber-reinforced tetraglycidyl-4,4′-diaminodiphenylmethane/4,4′-diaminodiphenylsulfone epoxy composites J. Appl. Polym. Sci. 2006 100 295 298 10.1002/app.23106
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).10.1002/app.23106
60. Singh SK Thermal properties of fluorinated graphene Phys. Rev. B 2013 87 104114 10.1103/PhysRevB.87.104114
Singh, S. K. et al. Thermal properties of fluorinated graphene. Phys. Rev. B 87, 104114 (2013).10.1103/PhysRevB.87.104114
61. Geada IL 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. 2018 9 716 10.1038/s41467-018-03137-8
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).10.1038/s41467-018-03137-8
62. Bryenton KR Adeleke AA Dale SG Johnson ER Delocalization error: The greatest outstanding challenge in density-functional theory WIREs Computational Mol. Sci. 2023 13 e1631 10.1002/wcms.1631
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).10.1002/wcms.1631
63. Zavitsas AA The relation between bond lengths and dissociation energies of carbon-carbon bonds J. Phys. Chem. A 2003 107 897 898 10.1021/jp0269367
Zavitsas, A. A. The relation between bond lengths and dissociation energies of carbon-carbon bonds. J. Phys. Chem. A 107, 897–898 (2003).10.1021/jp0269367
64. Smith, M. B. March’s advanced organic chemistry: Reactions, mechanism, and structure. 7th edn (Wiley, Hoboken, NJ, 2013).
65. Moon RJ Martini A Nairn J Simonsen J Youngblood J Cellulose nanomaterials review: Structure, properties and nanocomposites Chem. Soc. Rev. 2011 40 3941 3994 10.1039/c0cs00108b 21566801
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).21566801 10.1039/c0cs00108b
66. Brenner SS Tensile strength of whiskers J. Appl. Phys. 1956 27 1484 1491 10.1063/1.1722294
Brenner, S. S. Tensile strength of whiskers. J. Appl. Phys. 27, 1484–1491 (1956).10.1063/1.1722294
67. Bai YX Carbon nanotube bundles with tensile strength over 80 GPa Nat. Nanotechnol. 2018 13 589 595 10.1038/s41565-018-0141-z 29760522
Bai, Y. X. et al. Carbon nanotube bundles with tensile strength over 80 GPa. Nat. Nanotechnol. 13, 589–595 (2018).29760522 10.1038/s41565-018-0141-z
68. Yu MF Files BS Arepalli S Ruoff RS Tensile loading of ropes of single wall carbon nanotubes and their mechanical properties Phys. Rev. Lett. 2000 84 5552 5555 10.1103/PhysRevLett.84.5552 10990992
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).10990992 10.1103/PhysRevLett.84.5552
69. Zhao QZ Nardelli MB Bernholc J Ultimate strength of carbon nanotubes: A theoretical study Phys. Rev. B 2002 65 144105 10.1103/PhysRevB.65.144105
Zhao, Q. Z., Nardelli, M. B. & Bernholc, J. Ultimate strength of carbon nanotubes: A theoretical study. Phys. Rev. B 65, 144105 (2002).10.1103/PhysRevB.65.144105
70. Benito JA Manero JM Jorba J Roca A Change of young’s modulus of cold-deformed pure iron in a tensile test Metall. Trans. A 2005 36A 3317 3324 10.1007/s11661-005-0006-6
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).10.1007/s11661-005-0006-6
71. Clatterbuck DM Chrzan DC Morris JW The inherent tensile strength of iron Philos. Mag. Lett. 2002 82 141 147 10.1080/095008302317262642
Clatterbuck, D. M., Chrzan, D. C. & Morris, J. W. The inherent tensile strength of iron. Philos. Mag. Lett. 82, 141–147 (2002).10.1080/095008302317262642
72. Diddens I Murphy B Krisch M Muller M Anisotropic elastic properties of cellulose measured using inelastic X-ray scattering Macromolecules 2008 41 9755 9759 10.1021/ma801796u
Diddens, I., Murphy, B., Krisch, M. & Muller, M. Anisotropic elastic properties of cellulose measured using inelastic X-ray scattering. Macromolecules 41, 9755–9759 (2008).10.1021/ma801796u
73. Dri FL Hector LG Moon RJ Zavattieri PD Anisotropy of the elastic properties of crystalline cellulose i-beta from first principles density functional theory with van der waals interactions Cellulose 2013 20 2703 2718 10.1007/s10570-013-0071-8
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).10.1007/s10570-013-0071-8
74. Damirchi B Reaxff reactive force field study of polymerization of a polymer matrix in a carbon nanotube-composite system J. Phys. Chem. C. 2020 124 20488 20497 10.1021/acs.jpcc.0c03509
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).10.1021/acs.jpcc.0c03509
75. Islam MM Zou C van Duin ACT Raman S Interactions of hydrogen with the iron and iron carbide interfaces: A reaxff molecular dynamics study Phys. Chem. Chem. Phys. 2016 18 761 771 10.1039/C5CP06108C 26626108
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).26626108 10.1039/C5CP06108C
76. Liu LC 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 2011 115 11016 11022 10.1021/jp201599t 21888351
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).21888351 10.1021/jp201599t
77. Kanski M Development of a charge-implicit reaxff potential for hydrocarbon systems J. Phys. Chem. Lett. 2018 9 359 363 10.1021/acs.jpclett.7b03155 29291618
Kanski, M. et al. Development of a charge-implicit reaxff potential for hydrocarbon systems. J. Phys. Chem. Lett. 9, 359–363 (2018).29291618 10.1021/acs.jpclett.7b03155
78. Strachan A van Duin ACT Chakraborty D Dasgupta S Goddard WA Shock waves in high-energy materials: The initial chemical events in nitramine RDX Phys. Rev. Lett. 2003 91 098301 10.1103/PhysRevLett.91.098301 14525217
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).14525217 10.1103/PhysRevLett.91.098301
79. Shin YK Gai LL Raman S van Duin ACT Development of a reaxff reactive force field for the Pt-Ni alloy Catal. J. Phys. Chem. A 2016 120 8044 8055 10.1021/acs.jpca.6b06770
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).10.1021/acs.jpca.6b06770
80. Aryanpour M van Duin ACT Kubicki JD Development of a reactive force field for iron-oxyhydroxide systems J. Phys. Chem. A 2010 114 6298 6307 10.1021/jp101332k 20455552
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).20455552 10.1021/jp101332k
81. Clark SJ First principles methods using CASTEP Z. Kristallog 2005 220 567 570 10.1524/zkri.220.5.567.65075
Clark, S. J. et al. First principles methods using CASTEP. Z. Kristallog 220, 567–570 (2005).10.1524/zkri.220.5.567.65075
82. 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).
83. Ryckaert JP Ciccotti G Berendsen HJC Numerical-integration of cartesian equations of motion of a system with constraints - molecular-dynamics of n-alkanes J. Comput. Phys. 1977 23 327 341 10.1016/0021-9991(77)90098-5
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).10.1016/0021-9991(77)90098-5
84. 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).
85. 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).
86. Li XX Mo Z Liu J Guo L Revealing chemical reactions of coal pyrolysis with GPU-enabled reaxff molecular dynamics and cheminformatics analysis Mol. Simula. 2015 41 13 27 10.1080/08927022.2014.913789
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).10.1080/08927022.2014.913789
87. Thompson AP LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales Comput. Phys. Commun. 2022 271 108171 10.1016/j.cpc.2021.108171
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).10.1016/j.cpc.2021.108171
88. Han Y Development, applications and challenges of reaxff reactive force field in molecular simulations Front. Chem. Sci. Eng. 2016 10 16 38 10.1007/s11705-015-1545-z
Han, Y. et al. Development, applications and challenges of reaxff reactive force field in molecular simulations. Front. Chem. Sci. Eng. 10, 16–38 (2016).10.1007/s11705-015-1545-z
89. Mcdonald TRR Beevers CA The crystal and molecular structure of alpha-glucose Acta Crystallogr 1952 5 654 659 10.1107/S0365110X52001787
Mcdonald, T. R. R. & Beevers, C. A. The crystal and molecular structure of alpha-glucose. Acta Crystallogr 5, 654–659 (1952).10.1107/S0365110X52001787
90. Pramanik C Gissinger JR Kumar S Heinz H Carbon nanotube dispersion in solvents and polymer solutions: Mechanisms, assembly, and preferences ACS Nano 2017 11 12805 12816 10.1021/acsnano.7b07684 29179536
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).29179536 10.1021/acsnano.7b07684
91. Akkineni S Amyloid-like amelogenin nanoribbons template mineralization via a low-energy interface of ion binding sites Proc. Natl Acad. Sci. Usa. 2022 119 e2106965119 10.1073/pnas.2106965119 35522709
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).35522709 10.1073/pnas.2106965119
92. Odegard GM Molecular dynamics modeling of epoxy resins using the reactive interface force field Macromolecules 2021 54 9815 9824 10.1021/acs.macromol.1c01813
Odegard, G. M. et al. Molecular dynamics modeling of epoxy resins using the reactive interface force field. Macromolecules 54, 9815–9824 (2021).10.1021/acs.macromol.1c01813
93. Gaikwad PS Understanding the origin of the low cure shrinkage of polybenzoxazine resin by computational simulation Acs Appl. Polym. Mater. 2021 3 6407 6415 10.1021/acsapm.1c01164
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).10.1021/acsapm.1c01164
94. Odegard GM Accurate predictions of thermoset resin glass transition temperatures from all-atom molecular dynamics simulation Soft Matter 2022 18 7550 7558 10.1039/D2SM00851C 36149371
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).36149371 10.1039/D2SM00851C
95. Pisani WA Newman JK Shukla MK Multiscale modeling of polyamide 6 using molecular dynamics and micromechanics Ind. Eng. Chem. Res 2021 60 13604 13613 10.1021/acs.iecr.1c02440
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).10.1021/acs.iecr.1c02440
96. Barr, S. A. et al. Bond breaking in epoxy systems: A combined qm/mm approach. J. Chem. Phys. 144 (2016).
97. Check CE Gilbert TM 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. 2005 70 9828 9834 10.1021/jo051545k 16292812
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).16292812 10.1021/jo051545k
98. Wodrich MD Corminboeuf C Schreiner PR Fokin AA Schleyer PV How accurate are DFT treatments of organic energies? Org. Lett. 2007 9 1851 1854 10.1021/ol070354w 17417862
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).17417862 10.1021/ol070354w
99. 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).
100. Materials studio 2019 program suite and user guide. (Biovia/Dassault Systemes, Cambridge, UK, 2019).
101. Stukowski A Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool Modell. Simul. Mater. Sci. Eng. 2010 18 015012 10.1088/0965-0393/18/1/015012
Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool. Modell. Simul. Mater. Sci. Eng. 18, 015012 (2010).10.1088/0965-0393/18/1/015012
102. Humphrey W Dalke A Schulten K Vmd: Visual molecular dynamics J. Mol. Graph. Modell. 1996 14 33 38 10.1016/0263-7855(96)00018-5
Humphrey, W., Dalke, A. & Schulten, K. Vmd: Visual molecular dynamics. J. Mol. Graph. Modell. 14, 33–38 (1996).10.1016/0263-7855(96)00018-5
103. Lu T Chen FW Multiwfn: A multifunctional wavefunction analyzer J. Comput. Chem. 2012 33 580 592 10.1002/jcc.22885 22162017
Lu, T. & Chen, F. W. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 33, 580–592 (2012).22162017 10.1002/jcc.22885
104. Bone MA Howlin BJ Hamerton I Macquart T Automapper: A Python tool for accelerating the polymer bonding workflow in LAMMPS. Comp Mater. Sci. 2022 205 111204
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).
105. Zartman GD 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. 2010 114 1763 1772 10.1021/jp907012w
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).10.1021/jp907012w
106. Wu XW Moon RJ Martini A Tensile strength of iβ crystalline cellulose predicted by molecular dynamics simulation Cellulose 2014 21 2233 2245 10.1007/s10570-014-0325-0
Wu, X. W., Moon, R. J. & Martini, A. Tensile strength of iβ crystalline cellulose predicted by molecular dynamics simulation. Cellulose 21, 2233–2245 (2014).10.1007/s10570-014-0325-0
107. Pan ZL Li YL Wei Q Tensile properties of nanocrystalline tantalum from molecular dynamics simulations Acta Mater. 2008 56 3470 3480 10.1016/j.actamat.2008.03.025
Pan, Z. L., Li, Y. L. & Wei, Q. Tensile properties of nanocrystalline tantalum from molecular dynamics simulations. Acta Mater. 56, 3470–3480 (2008).10.1016/j.actamat.2008.03.025
108. Interface Force Field (IFF) and a surface model database. https://bionanostructures.com/interface-md/ (2013−2023).
109. 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).
110. 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, 10.5281/zenodo.12518681 (2024).
111. Iwamoto S Kai WH Isogai A Iwata T Elastic modulus of single cellulose microfibrils from tunicate measured by atomic force microscopy Biomacromolecules 2009 10 2571 2576 10.1021/bm900520n 19645441
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).19645441 10.1021/bm900520n
112. Lamb, S. & Bringas, J. E. Casti handbook of stainless steels & nickel alloys. 2nd edn (CASTI Pub., Edmonton, 2002).
113. Peckner, D. & Bernstein, I. M. Handbook of stainless steels. (McGraw-Hill, New York, 1977).
114. Trucano P Chen R Structure of graphite by neutron-diffraction Nature 1975 258 136 137 10.1038/258136a0
Trucano, P. & Chen, R. Structure of graphite by neutron-diffraction. Nature 258, 136–137 (1975).10.1038/258136a0
115. Liu Z Interlayer binding energy of graphite: A mesoscopic determination from deformation Phys. Rev. B 2012 85 205418 10.1103/PhysRevB.85.205418
Liu, Z. et al. Interlayer binding energy of graphite: A mesoscopic determination from deformation. Phys. Rev. B 85, 205418 (2012).10.1103/PhysRevB.85.205418
116. Zacharia R Ulbricht H Hertel T Interlayer cohesive energy of graphite from thermal desorption of polyaromatic hydrocarbons Phys. Rev. B 2004 69 155406 10.1103/PhysRevB.69.155406
Zacharia, R., Ulbricht, H. & Hertel, T. Interlayer cohesive energy of graphite from thermal desorption of polyaromatic hydrocarbons. Phys. Rev. B 69, 155406 (2004).10.1103/PhysRevB.69.155406
117. Oja V Suuberg EM Vapor pressures and enthalpies of sublimation of d-glucose, d-xylose, cellobiose, and levoglucosan J. Chem. Eng. Data 1999 44 26 29 10.1021/je980119b
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).10.1021/je980119b
118. Basinski ZS Humerothery W Sutton AL The lattice expansion of iron Proc. R. Soc. Lond., Ser. A 1955 229 459 467 10.1098/rspa.1955.0102
Basinski, Z. S., Humerothery, W. & Sutton, A. L. The lattice expansion of iron. Proc. R. Soc. Lond., Ser. A 229, 459–467 (1955).10.1098/rspa.1955.0102
119. 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. 2002 124 9074 9082 10.1021/ja0257319 12149011
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).12149011 10.1021/ja0257319
