Voltage equilibration for reactive atomistic simulations of electrochemical processes
We introduce EChemDID, a model to describe electrochemical driving force in reactive molecular dynamics simulations. The method describes the equilibration of external electrochemical potentials (voltage) within metallic structures and their effect on the self consistent partial atomic charges used in reactive molecular dynamics. An additional variable assigned to each atom denotes the local potential in its vicinity and we use fictitious, but computationally convenient, dynamics to describe its equilibration within not-simply connected metallic structures on-the-fly during the molecular dynamics simulation. This local electrostatic potential is used to dynamically modify the atomic electronegativities used to compute partial atomic changes via charge equilibration. Validation tests show that the method provides an accurate description of the electric fields generated by the applied voltage and the driving force for electrochemical reactions. We demonstrate EChemDID via simulations of the operation of electrochemical metallization cells. The simulations predict the switching of the device between a high-resistance to a low-resistance state as a conductive metallic bridge is formed and resistive currents that can be compared with experimental measurements. In addition to applications in nanoelectronics, EChemDID could be useful to model electrochemical energy conversion devices.
School of Materials Engineering and Birck Nanotechnology Center Purdue University, West Lafayette, IN 47906 USA
Keywords: molecular dynamics, electrochemistry, reactive potentials, electrochemical metallization cell
Molecular dynamics  (MD) modeling has become a key tool in many areas of science and engineering including chemistry, [2, 3] material science [4, 5] and biology . Cloud computing has significantly simplified access to these MD simulations  to the point where they can now be used in undergraduate education . While ab initio MD  (where atomic forces are computed from a quantum mechanical electronic structure calculation) plays an important role in many applications [10, 11], its computational intensity limits its use to relatively small systems and short times. Therefore, dynamical simulations with empirical interatomic potentials represent an indispensable tool to connect electronic processes with the response of materials or biological systems and in determining the atomistic mechanisms that govern their behavior. The last decades witnessed the development of accurate interatomic potentials for a wide range of materials and applications including metals [12, 13], semiconductors [14, 15], oxides [16, 17], and molecular systems [18, 19, 20]. More recently, reactive potentials [21, 22, 23] have demonstrated the ability to describe complex chemistry including the reaction and decomposition of high-energy density materials , the description of solid/liquid interfaces [25, 26] and the mechanisms of complex catalytic chemical reaction .
A wide range of applications of technological interest require the description of electrochemical processes; examples include batteries , capacitors , and electrochemical metallization (ECM) cells for nanoelectronics . Despite recent progress [31, 32] we lack a self-consistent approach to describe how externally applied electrochemical potentials affect partial atomic charges and drive electrochemical reactions. In this paper we introduce the electrochemical dynamics with implicit degrees of freedom (EChemDID) method that describes the electrochemical driving force resulting from the application of an external voltage to metallic electrodes in reactive MD simulations. A significant challenge addressed by our method is describing the equilibration of the external potential over metallic electrodes with changing composition and morphology as ions electrochemically dissolve into the electrolyte or are deposited. We demonstrate the power of EChemDID via simulations of nanoscale ECM cells of interest for resistive switching non-volatile memory applications. The original version of the approach used a cluster analysis to determine the extent of the electrodes and enabled the first atomistic description of the operation of ECM cells . This paper extends and validates the method, provides a detailed description of the underlying physics and further exemplifies its use.
The remainder of the paper is organized as follows. Section 2 describes the EChemDID method and its implementation within reactive MD simulations. Section 3 verifies our implementation and validates the model while Section 4 exemplifies its use to describe ECM cells of interest in nanoelectronics. Finally, conclusions are drawn in Section 5.
2 The EChemDID method
This section starts with a brief review of reactive MD simulations focusing on the self-consistent calculation of partial atomic charges using charge equilibration (QEq ). We then introduce the EChemDID approach to describe the application of an electrochemical potential to metallic electrodes.
2.1 Reactive MD and charge equilibration
Reactive interatomic potentials like ReaxFF , REBO [35, 36] or COMB , use the concept of bond order to describe covalent interactions that include bond stretch, angle bending and torsional potentials as well as terms to penalize over and under coordination. In addition to covalent interactions, van der Waals interactions describe London dispersion and Pauli repulsion and electrostatics are described in terms of partial atomic charges. Theses charges need to be computed self-consistently based on the atomic structure of the system to capture charge transfer during chemical reactions. ReaxFF uses the charge equilibration method , also known as electronegativity equalization (EEM ), where the atomic charges are obtained equating the electronic chemical potential of every atoms in the system. The atomic electronegativity is obtained from a simple expression for the total electrostatic energy of the system:
where the second term in the right-hand side (RHS) represents the electrostatic interaction between ions and ; takes the form of shielded Coulomb interaction. The first term represents the energy cost associated with changing the charge of individual atoms and depends on their electronegativity () and hardness (). Given an atomic structure, the equilibrium partial atomic charges are obtained by equating the chemical potential of each atom:
leading to equations with unknowns (charges for the atoms and the chemical potential ). The system of equations is closed by constraining the total charge in the system: .
Important recent advances to charge equilibration methods include methods based on charge transfer between bonded atoms like the Split Charge Equilibration (SQE ), the addition of integer charge changes to describe reactions with ions dissolved in an electrolyte  and more recently, the extension of SQE and its derivation from density functional theory (DFT) as proposed in the atom-condensed Kohn-Sham DFT approximated to second order method (ACKS2 ). In this paper we develop EChemDID starting from the original formulation of QEq that remains widely used, the method is easily extensible to SQE.
We are interested in simulating electrochemical reactions that originate from the application of an external electrochemical potential difference () between two electrodes; this has the effect of changing the relative energy of electrons in the two electrodes by . From the energy expression in Eq. 1 we see that this can be accomplished by changing the atomic electronegativity of the atoms in one electrode to and of those in the other electrode to . Changes in atomistic electronegativity has been used in the past to create external fields and electrostatic potentials [40, 41, 32]. The EChemDID method introduced in the next subsection, describes the process of equilibration of the applied electrochemical potential in non-simply-connected structure of each contact accounting for the loss and gain of atoms during electro-dissolution and metallization.
2.2 EChemDID equations for electrochemical potential equilibration
The problem we need to address now is how the external potential applied to a group of atoms in the electrodes propagates and equilibrates throughout the metallic structure given that their constituent atoms change as electrochemical reactions occur at the electrodes. We start by assigning an additional dynamical variable to each metallic atom that represents the local electrochemical potential in their vicinity: . This is done in the same spirit as the local electronic temperature in two-temperature models designed to capture the thermal transport role of conduction electrons in metals [42, 43] or to represent the temperature of internal degrees of freedom in coarse grain simulations . We consider that the electrochemical potential is externally applied to a pre-determined set of atoms located away from the chemically active regions, see Figure 1 (at 0 ps), and the electrochemical potential of these atoms remain constant throughout the simulation. We now propose a dynamical equation for the temporal evolution of the electrochemical potential of all atoms in the system.
The voltage applied to a group of atoms within the metallic electrode propagates through it at the speed of light following Maxwell’s wave equations  until the entire metal is at the same electrochemical potential. We are not interested in the actual dynamics of this voltage equilibration but on the much slower electrochemical reactions; thus we use fictitious diffusive dynamics to describe the equilibration process. We assume that the voltage within the metal follows:
where the dot denotes time derivative, is the Laplacian operator and is an effective diffusivity. If a voltage difference is applied to electrodes that are not connected to each other, the solution of Eq. 3 will result in the voltage being equilibrated within each electrode at a timescale that depends on the value of . Furthermore, we will assume that if an atom becomes detached from the electrode its electronegativity should go back to its atomic value ; that is, the external electrochemical potential should go to zero. The numerical solution of Eq. 3 is performed on-the-fly during the MD simulation using atoms as a grid and a local weighting function as in the eleDID method for thermal transport :
The first term on the RHS solves the diffusion equation (Eq. 3); denotes the distance between atoms and and is a local weighting function defined below. The second term in the RHS relaxes the electronegativity of atoms detached from the electrodes to their atomic value. is the relaxation rate, is a switching function that turns on the relaxation as an atom becomes detached from other metallic atoms and is the total metallic coordination of atom . This external electrochemical potential is then added to the atomic electronegativity and is used for charge equilibration.
The local weighting function takes the form:
with a normalization constant. The range of this weighting function () represents a critical separation distance below which two atoms are considered as part of the same metallic cluster and, thus, be at the same electrochemical potential. We define the normalization factor from a reference lattice (in dimensions, with neighbors at distance ) as:
is the total number of atoms and represents the total metallic coordination of atom in the reference lattice . The last term in Eq. 4 relaxes the external potential of isolated atoms to zero. The switching function increases from zero to one as the total metallic coordination decreases to zero as:
represents a critical metallic coordination below which an atom will evolve its electronegativity towards the equilibrium value. We take this value as that of the weighting function when two atoms are separated by a distance close to but smaller than the cutoff (Eq. 5). In practice, we take this value as .
Figure 1 exemplifies the equilibration of the chemical potential in a two electrodes setup. Each electrode consists of Cu atoms forming a perfect fcc lattice and an applied external electrochemical potential of V is applied to the regions shown at time ps. Solving the EChemDID equations (with atoms fixed in space in this first example) shows the equilibration of the voltage in each of the contacts.
The following Sections demonstrate the coupling of EChemDID with charge equilibration that generates electrostatics and with ionic dynamics to simulate electrochemical reactions.
2.3 EChemDID and electronic transport: estimating current densities
While the main purpose of the method is to equilibrate the electrochemical potential in metallic structures for electrochemical simulations, an interesting by-product of the approach is that Eqs. 3 and 4 () can be used to estimate electrical currents when a potential difference is applied across a metallic system. Assuming diffusive transport, a combination of Ohm’s law (where is the electrical current and the electrical resistivity) with the continuity equation (where is the electrical charge density) shows that the time derivative in the electrochemical potential in EChemDID Eq. 3 is proportional to the time derivative of the free charge density as . Let’s consider a continuous region of a metallic wire of cross sectional area sandwiched between the regions and , the time derivative of the total charge in can be expressed as:
with an atomic volume. From Eq. 4 (assuming ) can be written as the sum of two terms representing the current flowing to/from and :
Finally, under steady state, the total current density from the EChemDID calculation can be obtained as:
The equality of the currents coming in and out of the region of interest offers a simple verification of the numerical implementation.
2.4 EChemDID implementation
The EChemDID method has been implemented as an external LAMMPS  user-package  including a LAMMPS “compute” that solves the Laplacian in Eq. 4 and a “fix” that integrates the voltage diffusion in time (Eq. 3). The whole implementation is consistent with the parallel scheme employed in LAMMPS. The modifications are independent to the actual LAMMPS program, only minor modifications to the reax/c package  allow switching from regular QEq (with constant electronegativity) to dynamical electronegativity in EChemDID. In the following MD simulations, the interactions between atoms are described by ReaxFF for Cu-SiO from the Ref. . While EChemDID has been fully implemented following the model presented in the previous section, a simplified version will be applied in the following ignoring the relaxation term and we will refer to the numerical solution of the Laplacian as Eq. 4 (). This last approximation assumes that the particles retain memory of their previous electronegativity.
3 Implementation verification and EChemDID validation
We first verify our implementation of the EChemDID equilibration using atoms as a grid by analyzing the 1D propagation of the electrochemical potential through a metallic electrode. We then validate the approach by computing the electric fields generated by the atomic charges induced by the difference in potential between two electrodes.
3.1 Electrochemical potential equilibration and electrostatics
We start by verifying the EChemDID equilibration of voltage by comparing the numerical solution with an analytical one in a 1D transport over a contact of length with one end set to +4 V. The initial voltage through the sample is zero and the analytical solution for the diffusion equation leads to voltage profiles as a function of time given by:
Figure 2 compares the numerical solution (red circles) with the analytical one (lines) at different times; the EChemDID numerical integration uses a timestep of 0.05 fs. The excellent agreement verifies our numerical solution of the equation Eq. 4. For this example, the Cu atoms are fixed in space; coupling Eq. 4 with molecular dynamics will be demonstrated in Section 4.
Having demonstrated that EChemDID equilibrates the electrochemical potential within metallic electrodes we now show how changing the atomic electronegativity affects the QEq-derived atomic charges and, thus, creates electric fields. We do this first for a simple parallel plates capacitor setup separated by vacuum. The application of an electrochemical potential between infinite, parallel electrodes separated by vacuum of thickness should generate a constant electric field between them equal to and a surface charge density equal to (with the relative permittivity). After solving the EChemDID equations and performing charge equilibration we compute the magnitude of the electric field between the electrodes using a probe charge that scans space. The resulting force on the probe atom of charge is proportional to the electric field as . Figure 3 shows 2D maps of the resulting electric field generated by two types of electrodes. The electric field created by the flat electrodes represented in 3a is spatially constant whereas the roughness induced by the triangular electrode showed in 3b results in electric field enhancement at the tip of the electrode. In both cases depicted in Figure 3, the amplitude of the electric field agrees with the applied voltage with approximately of error.
3.2 Driving force for electrochemistry
In order to quantify the driving force for electrochemistry we compare the energy of a series of atomic configurations under different voltages. We move a probe consisting of a Cu atom through a 1nm-thick layer of crystalline silicon dioxide sandwiched between two copper electrodes. At each location of the probe atom, we relax the structure via energy minimization and plot in Figure 4 the resulting energy landscape. The voltage applied between the parallel metallic electrodes acts on the ion dissolved in the dielectric and the average slope of the energy landscape is given by: where =1 nm is the separation between electrodes and =0.35 e the (averaged) charge of the dissolved ion. Therefore, we can extract an effective voltage from the total energies as a function of ion position obtained from the simulations. Figure 4 compares the extracted effective (in parenthesis) to the applied voltages. We find a reasonable agreement between the two considering the large fluctuations originating from the rugged underlying atomistic energy landscape.
4 Simulating electrochemical reactions using EChemDID
We now focus on full EChemDID simulations where Eq. 4 () is solved together with a reactive MD simulation. To demonstrate the approach we simulate the operation of nanoscale ECM cells of interest in nanoelectronics . These resistance-switching devices consist of an electro-active electrode (typically Cu or Ag) and an electro-inactive one (often W or Pt) separated by a solid electrolyte (amorphous SiO in our case). When a voltage of the appropriate polarity is applied, the atoms in the active electrode dissolve into the electrolyte and are field-driven towards the inactive electrode. Eventually, this process leads to the formation of a metallic conductive filament between the two electrodes, with the subsequent decrease in cell resistance.
4.1 ECM simulation details
The simulation cell is composed by two metallic electrodes separated by amorphous SiO generated following an annealing procedure described in Ref.  and similar to the setup presented in Ref. . The active electrode is made of Cu and so is the inactive one, which is described as a rigid body to capture its electrochemical inertness. In order to mimic the roughness of real interfaces we randomly deleted 50 % of the Cu atoms part of the first layer at the interfaces. The structure is then relaxed via energy minimization and equilibrated using isobaric/isothermal MD for 50 ps. The ECM cell is 44 nm in cross section where we apply periodic boundary conditions and the separation between electrodes is approximately 1.0 nm.
To simulate the operation of the cell, we impose a constant potential to a thin layer of atoms (0.2 nm) at each boundary of the electrodes (in the direction perpendicular to the interfaces) resulting in a voltage of 8 V, characteristic during the forming phase of this type of device . The MD equations of motion are solved with a timestep of 0.5 fs within the canonical ensemble at 300 K. The voltage diffusion equation (Eq. 3) is integrated 10 times per MD timestep (for numerical stability) and the diffusion constant is chosen to be equal to 4 Å/fs with a cutoff distance 0.4 nm (Eq. 4). For the purpose of demonstration, we use EChemDID with we the coefficient set to 0; a movie of the ECM simulation is included in the Supplementary Material.
4.2 Resistance Switching
Figure 5 shows atomistic snapshots at various times during the switching of the ECM cell, for clarity only the copper atoms are shown; the bottom panel shows the time evolution of the corresponding switching state of the cell. The color on the top row of snapshots denotes the local electronegativity (blue is -4 V and red +4 V) and the bottom row uses color to denote partial atomic charges. We see that QEq correctly localizes charges at the interfaces between the copper electrodes and the amorphous silica with positive ions on the active electrode and negative at the inactive one. The low electronegativity in the active electrode leads to positive partial charges which facilitates its dissolution into the dielectric; this can be formally thought of as an oxidation reaction: Cu Cu. Similarly, the electronegativity of the dissolved copper ions increases as they become closer to the inactive electrode (within the cutoff R), lowering their partial charge until they bond to the electrode and become metallic following chemical reduction, formally: Cu Cu. Throughout the simulation, EChemDID maintains constant electrochemical potential within each electrode even as their morphology changes dramatically during the operation of the device. As shown above, the charged electrodes generate an electric field throughout the dielectric that acts as the driving force for the diffusion of copper ions inside the solid electrolyte from the active toward the inactive electrode during the forming or programing phase.
EChemDID also provides insightful information regarding the formation of a conducting filament. To study this process we use two approaches: i) a distance-based cluster analysis to determine whether the two electrodes are connected by a Cu filament as was done in Ref.  and ii) using the current computed from EChemDID Eq. 10. The distance-based cluster analysis uses a cutoff of 0.4 nm and we assign a value of 1 to the signal when the two electrodes are bridged (low resistivity state: LR), 0 otherwise (high resistivity state: HR); the red line in Figure 5 (bottom scale) shows a running average (with a time-window of 25 ps) of the resulting signal for clarity. Therefore, a signal between LR and HR corresponds to switching states where the cell oscillates between ON and OFF within the time-window of the average. The current represented on the bottom scale of Figure 5 (blue) has also been averaged with a time-window of 5 ps.
From the cluster analysis, we see that a filament was established for a short duration just before 0.5 ns and a more stable one forms at time 1.4 ns; this second filament remains for the duration of the run. The atomic configuration of the first bridge consists of a linear single-atom chain (inset ) and bridges the two electrodes for approximately 0.1 ns. As in our prior work  single-atom chains are metastable and do not survive long timescales. The second, more stable, filament (inset ) is slightly thicker with a zig-zag shape. We show in the insets of Figure 5 ( and ) the time-averaged coordination number of the atoms in the filament. While the atoms in the first bridge exhibit an averaged coordination number between 2 and 3, typical of a linear configuration, the atoms in the zig-zag chain have higher coordination number suggesting a stronger connection.
The current computed via Eq. 10 correlates with the distance-based cluster analysis however additional informations can be extracted from its amplitude. Contrary to the binary state of the cluster analysis, the amplitude of the current corroborates the strength of the connection, hence the quantity of current that can flow through the filament. We see on Figure 5 a peak in the current approximately 3 times larger than in in agreement with the inset showing a zig-zag configuration of the bridge. We demonstrate a direct connection between the atomic structure of the filament and the current leaving and entering the electrodes computed from the dynamical equilibration of the EChemDID method. This last analysis is the qualitative in silico equivalent to atomically-controlled quantum conductance experiments proposed recently .
5 Discussion and conclusions
In summary, we introduced a novel approach to describe electrochemical reaction using reactive molecular dynamics simulations. EChemDID describes the equilibration of electrochemical potential throughout non-simply connected metallic structures and this potential shifts the energetics used to obtain atomic charges in charge equilibration. Thus, the application of an external potential leads to atomic charges that can drive electrochemical process. We demonstrated that the driving force for electrochemistry is accurately described by EChemDID. In addition, the electrochemical potential equilibration is performed on-the-fly during the MD simulation and thus adjusts to changes of composition and topology of the electrodes. As atoms dissolve into the electrolyte their electronegativity evolves back to its original value or keeps memory of previous voltage. The implementation in LAMMPS is particularly straightforward and the method can be easily extended to more sophisticated charge equilibration methods.
EChemDID is not without limitations. The ions dissolved in the electrolyte could trap electrons and change their charge state; such processes are currently ignored in EChemDID. Recently, Müser and Dapp proposed a method to handle integer charge states and used it to explore redox reactions in a model battery; such an approach could be incorporated into EChemDID. In addition, while EChemDID enables an estimation of electronic currents, they do not affect atomic dynamics; processes like Joule heating and electromigration are ignored in this first version of the model.
The demonstration of the use of EChemDID to simulate the resistance switching in ECM cells shows the power of the approach which is generally applicable with reactive and non-reactive force fields (as long as dynamical charges are computed with QEq) and can be used to study electrochemical processes with atomistic detail. There is growing interest in ab initio descriptions of electrochemistry , and EChemDID can help bridge such methods toward full device simulations.
This work was supported by the FAME Center, one of six centers of STARnet, a Semiconductor Research Corporation program sponsored by MARCO and DARPA. We would like to express our gratitude to Steve Plimpton and Aidan Thompson for the interesting discussions related to the software LAMMPS. Finally, we thank nanoHUB.org and Purdue for the computational resources.
- Corresponding author: firstname.lastname@example.org
- BJ ALDER and TE WAINWRIGHT. PHASE TRANSITION FOR A HARD SPHERE SYSTEM. JOURNAL OF CHEMICAL PHYSICS, 27(5):1208–1209, 1957.
- DR HERSCHBACH. MOLECULAR-DYNAMICS OF ELEMENTARY CHEMICAL-REACTIONS (NOBEL LECTURE). ANGEWANDTE CHEMIE-INTERNATIONAL EDITION IN ENGLISH, 26(12):1221–1243, DEC 1987.
- M KARPLUS. Development of Multiscale Models for Complex Chemical Systems: From H+H2 to Biomolecules. Nobel lecture, 8 Dec., 2013.
- B.L Holian and P.S Lomdahl. Plasticity induced by shock waves in nonequilibrium molecular-dynamics simulations. Science, 280(5372):2085, 1998.
- Jun Song and W. A. Curtin. Atomic mechanism and prediction of hydrogen embrittlement in iron. Nature Materials, 12(2):145–151, 2013.
- M Karplus and JA McCammon. Molecular dynamics simulations of biomolecules. NATURE STRUCTURAL BIOLOGY, 9(9):646–652, SEP 2002.
- Alejandro Strachan, Gerhard Klimeck, and Mark Lundstrom. Cyber-enabled simulations in nanoscale science and engineering introduction. Computing in Science & Engineering, 12(2):12–17, 2010.
- SEAN P Brophy, Alejandra J Magana, and ALEJANDRO Strachan. Lectures and simulation laboratories to improve learners conceptual understanding. Advances in Engineering Education, 3(3):1–27, 2013.
- R CAR and M PARRINELLO. UNIFIED APPROACH FOR MOLECULAR-DYNAMICS AND DENSITY-FUNCTIONAL THEORY. PHYSICAL REVIEW LETTERS, 55(22):2471–2474, 1985.
- M Sprik, J Hutter, and M Parrinello. Ab initio molecular dynamics simulation of liquid water: Comparison three gradient-corrected density functionals. JOURNAL OF CHEMICAL PHYSICS, 105(3):1142–1152, JUL 15 1996.
- ME Tuckerman, D Marx, ML Klein, and M Parrinello. On the quantum nature of the shared proton in hydrogen bonds. SCIENCE, 275(5301):817–820, FEB 7 1997.
- MS DAW and MI BASKES. EMBEDDED-ATOM METHOD - DERIVATION AND APPLICATION TO IMPURITIES, SURFACES, AND OTHER DEFECTS IN METALS. PHYSICAL REVIEW B, 29(12):6443–6453, 1984.
- MI BASKES. MODIFIED EMBEDDED-ATOM POTENTIALS FOR CUBIC MATERIALS AND IMPURITIES. PHYSICAL REVIEW B, 46(5):2727–2742, AUG 1 1992.
- FH STILLINGER and TA WEBER. COMPUTER-SIMULATION OF LOCAL ORDER IN CONDENSED PHASES OF SILICON. PHYSICAL REVIEW B, 31(8):5262–5271, 1985.
- J TERSOFF. NEW EMPIRICAL-APPROACH FOR THE STRUCTURE AND ENERGY OF COVALENT SYSTEMS. PHYSICAL REVIEW B, 37(12):6991–7000, APR 15 1988.
- GV LEWIS and CRA CATLOW. POTENTIAL MODELS FOR IONIC OXIDES. JOURNAL OF PHYSICS C-SOLID STATE PHYSICS, 18(6):1149–1161, 1985.
- F. H. Streitz and J. W. Mintmire. Electrostatic potentials for metal-oxide surfaces and interfaces. Physical Review B, 50(16):11996, 1994.
- SL MAYO, BD OLAFSON, and WA GODDARD. DREIDING - A GENERIC FORCE-FIELD FOR MOLECULAR SIMULATIONS. JOURNAL OF PHYSICAL CHEMISTRY, 94(26):8897–8909, DEC 27 1990.
- WD Cornell, P Cieplak, CI Bayly, IR Gould, KM Merz, DM Ferguson, DC Spellmeyer, T Fox, JW Caldwell, and PA Kollman. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules (vol 117, pg 5179, 1995). JOURNAL OF THE AMERICAN CHEMICAL SOCIETY, 118(9):2309, MAR 6 1996.
- B. R. Brooks, C. L. Brooks, III, A. D. Mackerell, Jr., L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York, and M. Karplus. CHARMM: The Biomolecular Simulation Program. JOURNAL OF COMPUTATIONAL CHEMISTRY, 30(10, SI):1545–1614, JUL 30 2009.
- ACT van Duin, S Dasgupta, F Lorant, and WA Goddard. ReaxFF: A reactive force field for hydrocarbons. JOURNAL OF PHYSICAL CHEMISTRY A, 105(41):9396–9409, OCT 18 2001.
- HS JOHNSTON and C PARR. ACTIVATION ENERGIES FROM BOND ENERGIES .1. HYDROGEN TRANSFER REACTIONS. JOURNAL OF THE AMERICAN CHEMICAL SOCIETY, 85(17):2544–&, 1963.
- SW RICK, SJ STUART, and BJ BERNE. DYNAMICAL FLUCTUATING CHARGE FORCE-FIELDS - APPLICATION TO LIQUID WATER. JOURNAL OF CHEMICAL PHYSICS, 101(7):6141–6156, OCT 1 1994.
- A Strachan, ACT van Duin, D Chakraborty, S Dasgupta, and WA Goddard. Shock waves in high-energy materials: The initial chemical events in nitramine RDX. PHYSICAL REVIEW LETTERS, 91(9), AUG 29 2003.
- Joseph C. Fogarty, Hasan Metin Aktulga, Ananth Y. Grama, Adri C. T. van Duin, and Sagar A. Pandit. A reactive molecular dynamics simulation of the silica-water interface. JOURNAL OF CHEMICAL PHYSICS, 132(17), MAY 7 2010.
- Adri C. T. van Duin, Vyacheslav S. Bryantsev, Mamadou S. Diallo, William A. Goddard, Obaidur Rahaman, Douglas J. Doren, David Raymand, and Kersti Hermansson. Development and Validation of a ReaxFF Reactive Force Field for Cu Cation/Water Interactions and Copper Metal/Metal Oxide/Metal Hydroxide Condensed Phases. JOURNAL OF PHYSICAL CHEMISTRY A, 114(35):9507–9514, SEP 9 2010.
- KD Nielson, ACT van Duin, J Oxgaard, WQ Deng, and WA Goddard. Development of the ReaxFF reactive force field for describing transition metal catalyzed reactions, with application to the initial stages of the catalytic formation of carbon nanotubes. JOURNAL OF PHYSICAL CHEMISTRY A, 109(3):493–499, JAN 27 2005.
- JM Tarascon and M Armand. Issues and challenges facing rechargeable lithium batteries. NATURE, 414(6861):359–367, NOV 15 2001.
- Patrice Simon and Yury Gogotsi. Materials for electrochemical capacitors. NATURE MATERIALS, 7(11):845–854, NOV 2008.
- RaineR Waser and Masakazu Aono. Nanoionics-based resistive switching memories. NATURE MATERIALS, 6(11):833–840, NOV 2007.
- Wolf B. Dapp and Martin H. Mueser. Redox reactions with empirical potentials: Atomistic battery discharge simulations. JOURNAL OF CHEMICAL PHYSICS, 139(6), AUG 14 2013.
- Celine Merlet, Benjamin Rotenberg, Paul A. Madden, Pierre-Louis Taberna, Patrice Simon, Yury Gogotsi, and Mathieu Salanne. On the molecular origin of supercapacitance in nanoporous carbon electrodes. NATURE MATERIALS, 11(4):306–310, APR 2012.
- Nicolas Onofrio, David Guzman, and Alejandro Strachan. Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells. Nature Materials, 14:440–446, 2015.
- AK RAPPE and WA GODDARD. CHARGE EQUILIBRATION FOR MOLECULAR-DYNAMICS SIMULATIONS. JOURNAL OF PHYSICAL CHEMISTRY, 95(8):3358–3363, APR 18 1991.
- SJ Stuart, AB Tutein, and JA Harrison. A reactive potential for hydrocarbons with intermolecular interactions. JOURNAL OF CHEMICAL PHYSICS, 112(14):6472–6486, APR 8 2000.
- DW Brenner, OA Shenderova, JA Harrison, SJ Stuart, B Ni, and SB Sinnott. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. JOURNAL OF PHYSICS-CONDENSED MATTER, 14(4):783–802, FEB 4 2002.
- WJ MORTIER, SK GHOSH, and S SHANKAR. ELECTRONEGATIVITY EQUALIZATION METHOD FOR THE CALCULATION OF ATOMIC CHARGES IN MOLECULES. JOURNAL OF THE AMERICAN CHEMICAL SOCIETY, 108(15):4315–4320, JUL 23 1986.
- Razvan A. Nistor, Jeliazko G. Polihronov, Martin H. Muser, and Nicholas J. Mosey. A generalization of the charge equilibration method for nonmetallic materials. JOURNAL OF CHEMICAL PHYSICS, 125(9), SEP 7 2006.
- T. Verstraelen, P. W. Ayers, V. Van Speybroeck, and M. Waroquier. ACKS2: Atom-condensed Kohn-Sham DFT approximated to second order. JOURNAL OF CHEMICAL PHYSICS, 138(7), FEB 21 2013.
- Jiahao Chen and Todd J. Martinez. Charge conservation in electronegativity equalization and its implications for the electrostatic properties of fluctuating-charge models. JOURNAL OF CHEMICAL PHYSICS, 131(4), JUL 28 2009.
- O. Assowe, O. Politano, V. Vignal, P. Arnoux, B. Diawara, O. Verners, and A. C. T. van Duin. Reactive Molecular Dynamics of the Initial Oxidation Stages of Ni(111) in Pure Water: Effect of an Applied Electric Field. JOURNAL OF PHYSICAL CHEMISTRY A, 116(48):11796–11805, DEC 6 2012.
- H HAKKINEN and U LANDMAN. SUPERHEATING, MELTING, AND ANNEALING OF COPPER SURFACES. PHYSICAL REVIEW LETTERS, 71(7):1023–1026, AUG 16 1993.
- Keng-Hua Lin, Brad Lee Holian, Timothy C. Germann, and Alejandro Strachan. Mesodynamics with implicit degrees of freedom. JOURNAL OF CHEMICAL PHYSICS, 141(6), AUG 14 2014.
- A Strachan and BL Holian. Energy exchange between mesoparticles and their internal degrees of freedom. PHYSICAL REVIEW LETTERS, 94(1), JAN 14 2005.
- John David Jackson. Classical Electrodynamics. John Wiley & Sons, Inc., New Jersey, USA, 3 edition, 1991.
- Plimpton Steve. Fast parallel algorithms for Short-Range molecular dynamics. Journal of Computational Physics, 117(1):1–19, March 1995.
- Nicolas Onofrio and Alejandro Strachan. https://nanohub.org/groups/strachangroup/lammpsmodules.
- H. M. Aktulga, J. C. Fogarty, S. A. Pandit, and A. Y. Grama. Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques. PARALLEL COMPUTING, 38(4-5):245–259, APR-MAY 2012.
- Rainer Waser, Regina Dittmann, Georgi Staikov, and Kristof Szot. Redox-Based Resistive Switching Memories - Nanoionic Mechanisms, Prospects, and Challenges. ADVANCED MATERIALS, 21(25-26):2632+, JUL 13 2009.
- Nathan L. Anderson, Ravi Pramod Vedula, Peter A. Schultz, R. M. Van Ginhoven, and Alejandro Strachan. First-Principles Investigation of Low Energy E ‘ Center Precursors in Amorphous Silica. PHYSICAL REVIEW LETTERS, 106(20), MAY 17 2011.
- Tohru Tsuruoka, Kazuya Terabe, Tsuyoshi Hasegawa, Ilia Valov, Rainer Waser, and Masakazu Aono. Effects of Moisture on the Switching Characteristics of Oxide-Based, Gapless-Type Atomic Switches. ADVANCED FUNCTIONAL MATERIALS, 22(1):70–77, JAN 11 2012.
- Ilia Valov, Ina Sapezanskaia, Alpana Nayak, Tohru Tsuruoka, Thomas Bredow, Tsuyoshi Hasegawa, Georgi Staikov, Masakazu Aono, and Rainer Waser. Atomically controlled electrochemical nucleation at superionic solid electrolyte surfaces. NATURE MATERIALS, 11(6):530–535, JUN 2012.
- Nicephore Bonnet and Nicola Marzari. First-Principles Prediction of the Equilibrium Shape of Nanoparticles Under Realistic Electrochemical Conditions. PHYSICAL REVIEW LETTERS, 110(8), FEB 20 2013.