(373f) Toward a New Generation of Reactive Force Fields: Polarizable Reactive Force Field (ReaxFP)

Modern quantum mechanics (QM) methods can accurately describe intermolecular forces in solids. However, QM calculations are computationally expensive and impractical for the description of the dynamics of systems with hundreds of atoms or more and also for more than nanoseconds of dynamics. Thus, to describe the dynamics of much larger systems for much longer time scales such as dynamical interactions between and within macromolecules, proteins, ceramics, semiconductors, and metallic alloys, it is useful to have an analytic potential or force field (FF) that describes how the forces respond to changing interatomic distances. Development of accurate classical molecular dynamics (MD) FFs, reactive and non-reactive, relies upon accurate expressions for non-bond potential energy terms including electrostatic, van der Waals (VdW), and H-bond interactions. The bond energy contributions are well known from experimental studies and ab-initio calculations.

In particular, the description of the electrostatic interactions is of critical importance because reactive FFs change connectivity during the MD runs and thus require updates of assigned atomic charges to match the new electrostatic environment. Therefore, many of the current approaches which use fixed charges fail to accurately describe the charge distribution within molecules.  Since QM leads directly to the electron density throughout space, this might seem as a solution, but most FFs describe the electron density as atom centered point charges. Therefore, other ways (e.g. Mulliken charges, Electrostatic Potential charges, Bader charges, etc) must be used to convert electron densities to atom-centered charges. This can add additional uncertainties due to limitations of these methods. In addition, for large systems, such as virtual ligand screening of ligand-protein interactions, we would prefer to avoid QM calculations for millions of molecules in the library for nanoseconds of the simulations.  

For this reason, we propose a polarizable reactive FF (ReaxFP) where a Polarizable Charge Equilibration scheme (PQEq) is used for calculation of electrostatic energy. In PQEq model, the charge on each atom is partitioned into a Gaussian-shaped core with positive fixed charge (nucleus and inner electron) and a Gaussian-shaped shell with negative variable charge (valence electrons). Shell charges can flow from atom to atom based on our earlier Charge Equilibration (QEq) scheme. The shells can be displaced from the cores and the primary restoring force between a core and its shell given by the electrostatic interaction between the two charge distributions. The model uses three atomic parameters for each element (bond radius, electronegativity, and hardness) and 2nd and 4th order spring forces to adjust to the proper atomic polarizability for each element. The charge on each atom is obtained by requiring that the electron chemical potential be the same on all atoms. The optimal charges are obtained using pseudodynamics, which reduces the computational cost to O(N), a large improvement over the O(N³) required for matrix inversion.

To optimize the electrostatic energy parameters (i.e. PQEq parameters) a constrained optimization was performed on a training set containing a variety of structures for the desired atoms and molecules. Optimization of the three atomic parameters is carried out with respect to the reference charges that are obtained from the electrostatic potential (ESP) of accurate ab-initio calculations and also the Mulliken population method. The polarization energies of electric dipole probes computed by quantum mechanics (QM) are used in the next step to optimize the spring force constants of the model. 

ReaxFP uses similar functional forms for valence energy terms as in the reactive FF (ReaxFF), which our group at Caltech also has developed. However, all parameters of the ReaxFP including bonded and non-bonded energy terms must be developed for similar training sets. This is easily done in our Genetic Algorithm based Reactive Force Field optimizer (GARFfield) which is a multi-platform, multi-objective parallel hybrid genetic algorithm (GA) / conjugate-gradient (CG) based force field optimization framework.  It enables first-principles based force fields prepared from large quantum mechanical data sets.

We have already optimized ReaxFP for several materials of interest, including Cu_2-xSe and GaN. Cu_2-xSe based materials are in the forefront of the research for thermoelectric materials with both superionic and electronic conducting behavior. GaN is a very hard, mechanically stable wide bandgap semiconductor material with high heat capacity and thermal conductivity suitable for applications such as in solar cells. We use ReaxFP to address the issues related to the highly mobile cations and the polarization effects in Cu_2-xSe and to understand the electrical polarization effects at the GaN interfaces which is a key to proper device simulation. By accounting for polarization, ReaxFP provides an improved description of the dynamics of these materials.