Ab initio-based force fields

Table of contents:

Why use analytical forcefields?

In classical Molecular Dynamics (MD) and Monte Carlo (MC) simulations, the particle interactions are described by analytical potentials. These potentials can be divided into two main groups: empirical potentials and non-empirical ab initio-derived potentials. Empirical potentials are derived by fitting the potential expression to certain properties
observable by experiments, such as vaporization, lattice or binding enthalpies, densities, and various structural and dynamical properties. Such potentials are thus more or less guaranteed to reproduce those very properties, but their ability to reproduce other properties can be substantially worse. A further disadvantage with empirical potentials is that simple functional expressions with few parameters are needed to restrict the time taken in the development of the potential, and to keep the computational expense low. Thus, it is in practice impossible to find the "correct" function, especially for complex systems.

In contrast, in non-empirical potential expressions, the parameters are derived from (high-level) ab initio calculations. The properties of the model are thus not dependent on any experimental values, so if some calculated properties compare reasonably well with the corresponding experimental values, other properties derived from the simulation can be expected to be reasonably described as well, including the detailed molecular-level information.

We are working on several different forcefields, each too complex to rely entirely on empirical data. However, in certain cases it is still necessary to refer to experimental data for some properties. From ab initio calculations we make use of the following data:

interaction energies for various geometries
binding/dissociation energies
electrostatic potential
multipole moments
polarizability
lattice enthalpy
equation of state

These data is fitted to the analytical potential function (forcefield). Two examples of how to use this data are shown below.

In the best of cases, an ab initio-based force field can be used as a very efficient - and therefore very powerful - alternative to high-level ab initio calculations. Combined with MD simulations, one can then perform (a type of) "ab initio MD" simulations for very large-scale systems and long time spans. We are particularly interested in making such force-fields which are based on high-level ab initio calculations; this strategy is discussed more in the last section below.

Example 1: Fitting to the electrostatic potential of the DMSO molecule.

The first image shows the electrostatic potential around a DMSO molecule obtained from an ab initio calculation at the MP2/aug-cc-pVTZ level. Yellow to red correspond to negative values, i.e. where a positive ion would be attracted to the DMSO molecule. Blue areas correspond to positive values. The green "sheet" is the zero-level. Notice that the dipole moment (blue to red direction) is not parallel to the S-O bond (yellow and red balls). Also notice that, while the negative values are highest close to the oxygen atom, there is a clearly negative region close to the sulphur atom as well. This causes the zero level to show a "bent" behaviour around the molecule.

The second image shows the electrostatic potential resulting from point charges assigned to the locations of all nuclei. This is fairly similar to the standard CHELPG charges. Notice the failure of these charges to reproduce the behaviour in the negative value area close to the sulphur atom. The dipole moment vector has the correct orientation however; this results in negatively charged methyl groups.

The third image shows the electrostatic potential resulting from point charges assigned to the locations of all nuclei and two extra positions close to the sulphur atom. The values of the point charges as well as the positions of the extra positions were fitted to the electrostatic potential in the first picture above. Notice that the curvature of the zero level and the negative area close to the sulphur is more accurately reproduced. The total charge of the methyl groups are fairly close to zero, although it may not be apparent from the image.

Example 2: Fitting of interaction energies to a potential expression

Choosing or creating new potential expressions (forcefields) can be very time consuming, as a lot of care need to be taken to select reasonable interactions. Between which atoms should which potentials be used? Should all parameters be included for all pairs? Are there any significant many-body interactions present? Sometimes the only way to obtain the answer to these questions is to start performing a fit and learn as one goes. Often important information of the chemical properties can be learned by this procedure as well. Fitting potentials (in general) takes some time to learn and has even been described as an art, not science!

This image shows a set of so-called potential curves where each point corresponds to an ab initio interaction energy calculation at a certain geometry. Each curve corresponds to some geometry where only one internal coordinate is changed. The x-axis shows the value of the internal coordinate, which is some arbitrary distance between two molecules, and the y-axis the interaction energy of all the molecules. The red curves correspond to the ab initio interaction energies and the green curves to the model interaction energies. In this particular example, so-called Buckingham potentials are used. A drawback with these potentials is that the potential energy collapses towards -infinity at short distances. This is why the green curves bend quickly towards large negative values at short distances (they go outside the plot at a distance of 3). With care in the fitting procedure this has no effect in the resulting properties of the model, as the activation barrier for reaching this part of the potential can be kept high. In the figure above, the barrier is about 50 kcal/mol (about 200 kJ/mol), which is much too high to be crossed at room temperature (or higher, reasonable temperatures) at the timescales possible to cover during a simulation.

A procedure for increasing the accuracy of the calculated ab initio energies used in the potential fit using extrapolation techniques

The forcefield cannot be any better than the quality of the ab initio data used to "train" the forcefield. In many cases it is essential that the quality of the ab initio data is very high, for instance in chemical systems where the energy differences are small, or where cheap(er) ab initio methods are unable to describe the kind of interactions which are important for some systems. A notable example is weakly bound systems where van der Waals interactions are important. A method which partly helps overcome such limitations have been developed by us (see Spangberg, Acta Universitatis Upsaliensis Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 892. ISBN 91-554-5751-7 for details). In this procedure, a large number of energies at points related by a simple geometrical change are calculated using a relatively cheap method. A few of these points are then recalculated using a more accurate, but expensive, method. The difference in energy at each of these points can be seen as the "error" in the cheap method. Assuming that a simple function, f(r), can be fitted to the "error", we have:

E(r)-C(r)=f(r),

where E is the expensive energy and C the cheap energy. The simplest possible relation between the geometry points chosen is the distance between them, i.e. the expression above is applied to a series of points differing only in some distance coordinate r.

Potential energy curves at the HF (dotted line) and MP2 (solid line and five crosses) levels are shown. The correction function 447334.5exp(-8.989798r) +24.26814/r**4+5.208182 was obtained by fitting to the error for the five selected geometry points, and this function was added to the uncorrected potential energy curve, yielding the final "corrected" curve (dashed line).

Reactive force fields (ReaxFF)

The ReaxFF force field [van Duin et al, J. Phys. Chem. A, 2001, 105(4) 9396] is a reactive force field. As such, it allows for chemical bonds to be broken and formed dynamically during simulation, without the need for a predefined reaction path. The force field uses a bond length/bond order relation, that allows for a continuous description of the dissociation of covalent bonds. Terms for the bond, valence angle, and torsion angle energies are bond order dependent, i.e. dependent on the local environment of each atom. The force field also includes non-bonded interaction terms (van der Waals and Coulomb) that are calculated between every pair of atoms. Atomic charges for the Coulomb interaction are determined by the Electronegativity Equalization Method.