OPLS-AA Usage Notes

Background

In recent years, all-atom molecular dynamics simulations have been widely used in studies of polymers and condensed matter materials. Compared with coarse-grained models, all-atom models retain specific chemical structures and can therefore describe local conformations, functional-group interactions, and chain-segment motions.

To compare the applicability of different classical force fields, we first examine the parameter sources, potential energy functions, and automatic parameterization methods of OPLS-AA.

Introduction

OPLS-AA (Optimized Potentials for Liquid Simulations - All Atom) is a classical all-atom force field for organic molecules and condensed-phase systems. It originated from OPLS-UA, developed by William L. Jorgensen and co-workers in the 1980s, and was initially used to simulate intermolecular interactions in organic liquids. The early OPLS-UA (United Atom) model merged some hydrogen atoms into the corresponding heavy-atom interaction centers.

In 1996, Jorgensen, Maxwell, and Tirado-Rives systematically established the OPLS-AA force field. This model explicitly describes hydrogen atoms and uses quantum chemical conformational energies and liquid experimental data to determine dihedral and nonbonded interaction parameters. Since then, OPLS-AA has been continuously extended and revised for parameters of different functional groups and chemical environments.

In 2015, OPLS-AA/M was optimized for the conformational energies of peptides and proteins, with a focus on revising backbone and side-chain dihedral parameters. In 2019, this series was further extended and improved for the description of RNA.

In 2017, LigParGen provided a method for automatically generating OPLS-AA parameters and molecular topologies, and it supports CM1A-LBCC atomic charges. Public implementations of the LigParGen Python wrapper already exist, but the core parameter assignment relies on the non-open-source BOSS program and has many limitations.

In 2019, PolyParGen adopted a molecular fragmentation method to generate force field parameters for larger molecules and polymers.

In addition, Schrödinger published the commercial force field OPLS4 in 2021, which improved parameters and prediction accuracy for a broader chemical space.

OPLS-AA Form

Jorgensen, Maxwell, and Tirado-Rives published Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids in Journal of the American Chemical Society in 1996.

This work mainly borrowed the bond stretching and bond angle parameters of the AMBER all-atom force field, used rotational energy profiles calculated at RHF/6-31G*//RHF/6-31G* to fit dihedral parameters, and optimized nonbonded interactions using experimental liquid density and heat of vaporization data.

The total potential energy of OPLS-AA can be written as:

$$U=U_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}+U_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}+U_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}+U_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}+U_{\mathrm{L}\mathrm{J}}+U_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}.$$

The first four terms describe bonded interactions, and the last two terms describe nonbonded interactions.

Bond stretching

Bond stretching uses the harmonic oscillator approximation:

$$U_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}=\sum _{(ij)\in B}\frac{k_{ij}}{2}(r_{ij}-r_{ij}^{(0)}{)}^2.$$

Here, $r_{ij}^{(0)}$ is the equilibrium bond length and $k_{ij}$ is the force constant. Only the second-order expansion of the potential energy around the equilibrium bond length is considered. This is a very simple harmonic approximation.

The $1/2$ coefficient here is a parameter-definition convention. Different force field files may absorb this coefficient into the force constant.

Angle bending

Bond angles use a similar harmonic approximation:

$$U_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}=\sum _{(ijk)\in A}\frac{k_{ijk}}{2}({\theta }_{ijk}-{\theta }_{ijk}^{(0)}{)}^2.$$

Here, ${\theta }_{ijk}^{(0)}$ is the equilibrium bond angle. This ignores higher-order anharmonic terms and explicit coupling between bond lengths and bond angles.

Proper torsion

The classical OPLS dihedral potential uses a Fourier series truncated at the fourth harmonic:

$$\begin{array}{rl}{U}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}(\varphi )=& \frac{V_1}{2}(1+\cos \varphi )\\ & +\frac{V_2}{2}(1-\cos 2\varphi )\\ & +\frac{V_3}{2}(1+\cos 3\varphi )\\ & +\frac{V_4}{2}(1-\cos 4\varphi ).\end{array}$$

Here, $V_1,\dots ,V_4$ control the relative energies and rotational barriers of different dihedral conformations. Compared with bond stretching and angle bending, the energy required for rotation about many single bonds is lower, so multiple thermodynamically accessible conformations can be generated. Therefore, for flexible polymers, the energy distribution of backbone dihedrals significantly affects the local flexibility of the chain, conformational statistics, and conformational transition rates.

Improper torsion

Improper dihedrals are mainly used to constrain local geometry, for example, to maintain the planarity of aromatic rings or other planar functional groups. Depending on the specific force field parameters and program implementation, a harmonic potential or a periodic potential may be used. For example, the harmonic form is:

$$U_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}}=\sum _{\xi }\frac{k_{\xi }}{2}(\xi -{\xi }_0{)}^2.$$

Here, $\xi$ is the improper dihedral angle and ${\xi }_0$ is the reference angle.

Nonbonded interactions

The nonbonded interactions of OPLS-AA consist of the Lennard-Jones potential and the Coulomb potential.

The Lennard-Jones potential is:

$$U_{\mathrm{L}\mathrm{J}}=\sum _{i<j}{f}_{ij}^{\mathrm{L}\mathrm{J}}4{ϵ}_{ij}[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6].$$

Here, ${ϵ}_{ij}$ is the potential well depth, and ${\sigma }_{ij}$ is the distance parameter corresponding to the zero-crossing position of the LJ potential.

Traditional OPLS-AA describes only 1–4 (i.e., separated by three covalent bonds) and higher Lennard-Jones interactions for intramolecular interactions; nonbonded parameters use geometric mixing rules:

$$\boxed{{\displaystyle {\sigma }_{ij}=\sqrt{{\sigma }_i{\sigma }_j},{ϵ}_{ij}=\sqrt{{ϵ}_i{ϵ}_j}.}}$$

Electrostatic interactions use a fixed point-charge model:

$$U_{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}=\sum _{i<j}{f}_{ij}^{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}\frac{q_i{q}_j}{4\pi {ϵ}_0{r}_{ij}}.$$

Standard OPLS-AA applies a special $1/2$ scaling to 1–4 nonbonded interactions (i.e., connected through three covalent bonds):

$$U_{14}=0.5U_{14}^{\mathrm{L}\mathrm{J}}+0.5U_{14}^{\mathrm{C}\mathrm{o}\mathrm{u}\mathrm{l}}.$$

That is, both the LJ and Coulomb 1–4 nonbonded interactions must be multiplied by 0.5; together with the proper torsion term, they define the dihedral angle. Therefore, when parameterizing a certain dihedral angle, one should compare the complete potential energy change and cannot consider only the explicit Fourier dihedral term. In particular, note that the scaling coefficients of 1–4 interactions are part of the force field definition. After modifying these coefficients, it is usually necessary to recheck or refit the related dihedral parameters.