Molecular Mechanics' Cornerstone: Force Fields - Part 2

In molecular modeling, the factors that determine the energy of a system are not limited to bonded interactions such as bonds, angles, and torsions. Even when atoms are not directly connected through covalent bonds, they can still influence each other via non-bonded interactions. These interactions form the basis of many physical and biological processes, including molecular folding, stability, binding affinity, and phase behavior. In force fields, non-bonded interactions are typically described by two main components: electrostatic interactions (Coulomb interactions) and van der Waals interactions. In this article, we will explore the fundamental principles of non-bonded interactions and examine the different force fields used in molecular simulations.

Coulomb Potential

Coulomb’s Law

Ionic interactions between groups carrying full or partial charges can be approximately expressed for each atom pair (i, j) using Coulomb’s law.

Here, F denotes the force, while qᵢ represents the effective charge of atom i. If the particles carry the same charge, the force is positive (repulsion), whereas if they carry opposite charges, the force is negative (attraction). This gives rise to the following Coulomb potential.

Here, ε denotes the dielectric constant, while  K_coul represents the conversion factor required to obtain the energy in kcal/mol with the chosen charge units.

Electrostatic interactions play a crucial role in the stability of biomolecular conformations in the solvent environment and in bringing together residues that are distant in the primary sequence within the folded structure (1).

Figure 1. Coulomb potentials (K_coul qᵢqⱼ / rᵢⱼ), calculated for C...C and O...H interactions using CHARMM partial charge parameters (in esu units). −0.1800 (arginine carbon), −0.8340 (water oxygen), and 0.4170 (water hydrogen) (1).

Van der Waals Potential

For a non-bonded distance (rᵢⱼ), the van der Waals potential takes the commonly used 6/12 Lennard-Jones form in force fields:

The attraction and repulsion coefficients (A and B) depend on the types of the two interacting atoms. The r⁶ term accounts for the attractive component arising from quantum mechanics, while the r¹² term is mostly chosen for computational convenience. Together, these terms explain why atoms repel each other at very short distances but attract each other when approaching an optimal interatomic separation (1).

Figure 2. Change in the energy of a molecule as a function of the distance between non-bonded atoms. Atoms A and B may be within the same molecule, or the interactions may be intermolecular. The minimum energy occurs due to van der Waals interactions. For small apolar atoms, the minimum energy point corresponds to a decrease of a few kJ·mol⁻¹ (e.g., for CH₄/CH₄, Eₘᵢₙ=−1.2 kJ·mol⁻¹), but very short distances can destabilize a molecule by tens of kJ·mol⁻¹ due to non-bonded interactions (2).

The parameters of the van der Waals potential can be obtained either by fitting them to lattice energies and crystal structures or through liquid simulations designed to reproduce liquid properties (3).

Van der Waals energy can be interpreted as the non-bonded, apolar component of the interaction, separate from electrostatic energy arising from atomic charges. For example, it may involve the interaction between two methane molecules or between two methyl groups at opposite ends of the same molecule. The van der Waals energy approaches zero at large interatomic distances and becomes highly repulsive at short distances. From a quantum mechanical perspective, this repulsion arises from the overlap of electron clouds of two atoms; oppositely charged electrons repel each other. At intermediate distances, there is a slight attraction between the electron clouds due to induced dipole–dipole interactions, which physically arises from electron correlation. Even if a molecule (or part of it) does not have a permanent dipole moment, the motion of electrons can create a transient, slightly uneven charge distribution at a given moment. This dipole induces a polarization in a neighboring molecule (or another part of the same molecule), leading to attraction (4).

Types of Force Fields

There are many different force fields used in molecular simulations. These force fields differ in terms of the functional form of each energy term, the inclusion of cross terms, and the data used to parameterize them.

There are two general trends for force fields. If a force field is primarily designed to study large systems such as proteins or DNA, the functional forms are kept as simple as possible. This typically involves using only harmonic functions for bond length and bond angle potentials, excluding cross terms, and employing the Lennard-Jones potential for van der Waals interactions. Such force fields are commonly referred to as “harmonic,” “diagonal,” or “Class I” (4). Examples include AMBER (5), CHARMM (6), GROMOS (7), and DREIDING (8).

The other approach focuses on reproducing small- and medium-sized molecules with high accuracy. These force fields include a series of cross terms, use at least cubic or quartic expansions for bond length and bond angle potentials, and, if necessary, apply an exponential-type potential for van der Waals interactions. The development of small molecule force fields aims not only to reproduce geometries and relative energies but also vibrational frequencies, and they are typically referred to as “Class II” force fields. Examples include COMPASS (9) and CVFF (10).

Advanced improvements that incorporate, for example, hyperconjugation effects and electronic polarization using parameters dependent on neighboring atom types are defined as “Class III” force fields. An example is AMOEBA (11).

Force fields designed for studying macromolecules can be simplified by not explicitly considering hydrogens – this is known as the united atom approach. For instance, this option is available in the AMBER, CHARMM, GROMOS, and DREIDING force fields. The advantage of united atoms is that they reduce the number of variables by approximately 2–3 times, allowing the simulation of larger systems. Naturally, the coarser the atomic representation, the less detailed the final results. The choice of representation, and therefore the type of force field to be used, depends on the type of information sought. For example, if the goal is to study the geometries and relative energies of different conformations of a hexose, a detailed force field is required. However, if the objective is to investigate the dynamics of a protein composed of hundreds of amino acids, a coarse-grained model using each amino acid as a single unit is the only feasible option given the system size (4).

In conclusion, non-bonded interactions play a critical role in understanding the behavior of molecular systems. Van der Waals and electrostatic forces determine the stability of a system depending on interatomic distances and charge distributions. Accurate modeling of these interactions is essential for both structural accuracy and the reliability of dynamic behavior in molecular simulations. The parameterization of force fields and the choice of an appropriate approach depend on the size of the system being simulated and the type of information sought; simplified models are sufficient for large biomolecules, whereas small- and medium-sized molecules require more detailed force fields. Thus, the study of non-bonded interactions forms a cornerstone of molecular modeling, enabling the understanding of mechanisms at the atomic level as well as applications in biomolecular design and engineering.

References

1. Schlick, T. (2010). Molecular Modeling and Simulation: An Interdisciplinary Guide. 2nd Edition, Springer, Berlin. https://doi.org/10.1007/978-1-4419-6351-2

2. Errol G. Lewars. (2004). Computational Chemistry. Introduction to the Theory and Applications of Molecular and Quantum Mechanics. 3rd Edition.

3. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A Kollman. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Amer. Chem. Soc., 117:5179–5197, 1995.

4. Jensen, F. (2017). Introduction to computational chemistry (Third edition). John Wiley & Sons, Ltd.

5. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmayer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197.

6. MacKerell, A. D.; Bashford, D.; Bellott; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616.

7. Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676.

8. S. L. Mayo, B. D. Olafson, and W. A. Goddard III, “DREIDING: A generic force field for molecular simulations,” J. Phys. Chem. 94, 8897–8909 (1990).

9. Sun, H., Ren, P., & Fried, J. R. (1998). Condensed-phased optimized molecular potential for atomistic simulation studies. Comput. Theor. Polym. Sci, 8(1/2), 229.

10. Ritschl, F., Fait, M., Fiedler, K., Köhler, J. E., Kubias, B., & Meisel, M. (2002). An Extension of the Consistent Valence Force Field (CVFF) with the Aim to Simulate the Structures of Vanadium Phosphorus Oxides and the Adsorption of n‐Butane and of 1‐Butene on their Crystal Planes. Zeitschrift für anorganische und allgemeine Chemie, 628(6), 1385-1396.

11. Ren, P., & Ponder, J. W. (2003). Polarizable atomic multipole water model for molecular mechanics simulation. The Journal of Physical Chemistry B, 107(24), 5933-5947.