Molecular Simulation/Printable version

Molecular Simulation

The current, editable version of this book is available in Wikibooks, the open-content textbooks collection, at
https://en.wikibooks.org/wiki/Molecular_Simulation

Permission is granted to copy, distribute, and/or modify this document under the terms of the Creative Commons Attribution-ShareAlike 3.0 License.

Charge-Charge Interactions

All intramolecular and intermolecular interactions are the result of electrostatic interactions between charged particles. All molecules are comprised the three subatomic particles: protons, neutrons, and electrons. Neutrons do not carry a charge, but protons and electrons carry charges with equal magnitude but opposite sign. The magnitude of these charges are fixed. This value is elementary charge, e. By convention, protons are defined as having positive charges and electrons are defined as having negative charges. The magnitude of these charges have a constant value known as the elementary charge, e=1.602176565(35) × 10^-19 C. ε₀ is the vacuum permittivity constant, which is equal to 8.854187817... 10^-12 F/m (farads per metre).

The potential energy surface of the interaction of two charged particles, calculated using Coulomb's law.

The force between two charged particles due to this electrostatic interactions is,

Coulomb's law (force) $F(r)={\frac {-1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}^{2}}}$

In this equation, $r_{AB}$ is the distance between particles A and B. The charge of a particle is given by the variable q. A charge is a scalar quantity with a sign and a magnitude.

It is often more convenient to discuss intermolecular forces in terms of the potential energy of the interaction. The potential energy of the interaction of two charged particles, labelled A and B, can be determined by integrating the force experienced between the particles if they were moved from infinite separation where the intermolecular interaction is zero, to the distance ($r_{AB}$) they are actually separated by,

${\mathcal {V}}(r)=\int _{\infty }^{r}{\frac {-1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}^{2}}}dr$

$={\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}}}{\bigg |}_{\infty }^{r_{AB}}$

$=\left[{\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}}}\right]-\left[{\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{\infty }}\right]={\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}}}$

Coulomb's law (potential energy) ${\mathcal {V}}(r)={\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}}}$

Ionic molecules have charge-charge Coulombic interactions. If the charges have the same sign (e.g., two positive ions), the interaction is repulsive. If the charges have the opposite sign, the interaction is attractive.

Example

What is the potential energy of the electrostatic interaction between an Na⁺ ion and a Cl^- ion separated by 5.00 Å?

Solution: These are both ions at intermediate distances, so their electrostatic interaction can be approximated using Coulomb's law. The charge-charge distance is given (5 Å). The charges, q_A and q_B, are the elementary charge of a proton/electron for the Na+ cation and the Cl- anion, respectively. These distances and charges must be converted to the SI units of m and C, respectively.

${\mathcal {V}}(r)={\frac {1}{4\pi \epsilon _{0}}}{\frac {q_{A}q_{B}}{r_{AB}}}$

$={\frac {1}{4\pi \epsilon _{0}}}{\frac {(1e)(-1e)}{5.00{\text{Å}}}}$

$={\frac {1}{4\pi (8.8541\times 10^{-12}{\text{ F m}}^{-1})}}{\frac {(1.602\times 10^{-19}C)(-1.602\times 10^{-19}C)}{5.00\times 10^{-10}m}}$

$=-4.61\times 10^{-19}{\text{ J}}$

$=-278{\text{ kJ/mol}}$

Further Reading

Wikipedia:Coulomb's law

Charge-Dipole Interactions

A cation (+) has an attractive interaction with the negative pole of a dipolar molecule, but a repulsive interaction with the positive pole of a dipolar molecule.

Charge-Dipole Interactions

Vector geometry can be used to derive an expression for the potential energy of an ion interacting with a dipolar molecule in terms of the distance and angle of orientation between the ion and the center of the dipole.

Charge-Dipole interactions occur in the presence of a atom with a formal net charge such as Na⁺ (q_ion = +1) or Cl^- (q_ion=-1) and a dipole. A dipole is a vector ($\mu$) which connects two charged species of different signs i.e (q_ion=+1 with q_ion=-1 NaCl) over a distance $r$

The dipole moment of a molecule depends on a few factors.

1. Differences in electronegativity of the atoms in the molecule. The larger the difference, the greater the dipole moment.
2. The distance between the positive and negative ends of the dipole. A larger distance corresponds to a larger dipole moment.
3. Multiple polar bonds pointing in the same direction increase the molecular dipole moment were as if they are pointing in opposite directions the dipoles will cancel out each other.

For diatomic molecules, the dipole moment can be calculated by the formula $\mu =qr$ were $\mu$ is the dipole moment q is the atomic charges and r is the distance.

polar bonds point in the same direction increase the molecular dipole moment

polar bonds in opposite directions will cancel out

from the potential energy of a Charge-Dipole interaction equation

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos \theta }{4\pi \epsilon _{0}r^{2}}}$

you can see that there are 4 factors that affect the potential energy

• the distance $(1/r^{2})$ from the charge (q_ion) to the dipole ($\mu$)
• the magnitude of the dipole ($\mu$)
• the magnitude of the charge (q_ion)
• the angle of the interaction $(\theta )$

Depending on the angle the interactions can be repulsive when the charge is close to the like side of the dipole (i.e, + +) attractive if the charge is close to the opposite charged end of the dipole (i.e., + -) and zero if the positive and negative ends of the dipole are equal distances from the charge (e.g., perpendicular).

Derivation of Charge-Dipole Interactions Energy

Express interaction as the sum of Coulombic interactions.

${\mathcal {V}}(r)=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}|{\overrightarrow {AB}}|}}+{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}|{\overrightarrow {AC}}|}}$

$=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {1}{|{\overrightarrow {AB}}|}}-{\frac {1}{|{\overrightarrow {AC}}|}}\right]$

Distance from charges to poles

$|{\overrightarrow {AB}}|={\sqrt {(r-0.5l\cos \theta )^{2}+(0.5l\sin \theta )^{2}}}\approx {\sqrt {(r-0.5l\cos \theta )^{2}}}=r-0.5l\cos \theta$

This approximation is valid of the distance between the charge and dipole (r) is much greater than the length of the dipole, making the $0.5l\sin \theta ^{2}$ term negligible.

$|{\overrightarrow {AC}}|={\sqrt {(r+0.5l\cos \theta )^{2}+(0.5l\sin \theta )^{2}}}\approx {\sqrt {(r+0.5l\cos \theta )^{2}}}=r+0.5l\cos \theta$

Substitute lengths into the potential energy equation

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {1}{r+0.5l\cos \theta }}-{\frac {1}{r-0.5l\cos \theta }}\right]$

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {r-0.5l\cos \theta }{(r+0.5l\cos \theta )(r-0.5l\cos \theta )}}-{\frac {r+0.5l\cos \theta }{(r+0.5l\cos \theta )(r-0.5l\cos \theta )}}\right]$

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {(r-0.5l\cos \theta )-(r+0.5l\cos \theta )}{r^{2}-0.25l^{2}\cos ^{2}\theta }}\right]$

$=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {(l\cos \theta )}{(r^{2}-0.25l^{2}\cos ^{2}\theta )}}\right]$

$r^{2}>>0.25l^{2}\cos ^{2}\theta$

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}q_{d}}{4\pi \epsilon _{0}}}\left[{\frac {l\cos \theta }{r^{2}}}\right]$

$\mu =q_{d}l$

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos \theta }{4\pi \epsilon _{0}r^{2}}}$

Example: Charge-Dipole Interaction

What is the interaction energy of a ammonia molecule (μ=1.42 D) with a Cl^- ion that is 2.15 Å away when the dipole moment of ammonia is aligned with the Cl^- ion at (θ=0)?

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos(\theta )}{4\pi \epsilon _{0}r^{2}}}$

${\mathcal {V}}(r,\theta )=-{\frac {(-1e)(1.42D)\cos(0)}{4\pi \epsilon _{0}(2.15)^{2}}}$

${\mathcal {V}}(r,\theta )=44\;{\text{kJ/mol}}$

Dipole-Dipole Interactions

Two molecules that possess permanent dipole moments have a dipole-dipole interaction. Attractive interactions occur between the opposite-sign poles of the molecules (+/-). Repulsive interactions occur between the like-sign poles (+/+ and -/-).

Dipole-dipole interactions occur between two molecules when each of the molecules has a permanent dipole moment. These interactions can be attractive or repulsive depending on the variable orientation of the two dipoles.

A permanent dipole moment occurs when a molecule exhibits uneven charge distribution. This is caused by atoms of different electronegativities being bonded together. The electron density of the molecule will be polarized towards the more electronegative atoms. This excess of negative charge in one part of the molecule and the deficit of charge in another area of the molecule results in the molecule having a permanent dipole moment.

The orientation of the two interacting dipoles can be described by three angles (θ₁, θ₂, Φ). Angle Φ is the dihedral angle between the two dipoles, which corresponds to the rotation. It is this angle which accounts for the greatest change in potential energy of the electrostatic interaction. Attractive dipole-dipole interactions occur when oppositely charged ends of the dipole moments of two molecules are in closer proximity than the likely charged ends. When the opposite and like charged ends of the two dipoles are equidistant, there is a net-zero electrostatic interaction. This net-zero interaction occurs when the dihedral angle Φ is at 90°, i.e. the two dipoles are perpendicular to each other.

Two dipoles arranged in an orientation with a net-zero electrostatic interaction because the opposite-signed poles are at the same distance as the like-signed poles.

Dipole-dipole interaction in repulsive-parallel interaction orientation

The interaction between two dipolar molecules (AB and CD) can be described using vector geometry.

Derivation

The potential energy of a dipole-dipole interaction can be derived by applying a variation of Coulomb's law for point charges. By treating the two dipoles as a series of point charges, the interaction can be described by a sum of all charge-charge interactions, where the distance between charges is described using vector geometry.

${\mathcal {V}}(r)={\frac {q_{d1}q_{d2}}{4\pi \epsilon _{0}}}{\frac {1}{\left\vert AC\right\vert }}-{\frac {q_{d1}q_{d2}}{4\pi \epsilon _{0}}}{\frac {1}{\left\vert AD\right\vert }}-{\frac {q_{d1}q_{d2}}{4\pi \epsilon _{0}}}{\frac {1}{\left\vert BC\right\vert }}+{\frac {q_{d1}q_{d2}}{4\pi \epsilon _{0}}}{\frac {1}{\left\vert BD\right\vert }}$

Similarly to the derivation of a charge-dipole interaction, the derivation for dipole-dipole interactions begins by defining the distance between the centre of the two dipoles as r, and the length of each dipole as l. Using a method similar to the derivation of a charge-dipole interaction, trigonometry can be used to determine the distance of each vector corresponding to an interaction: $\left\vert AC\right\vert$, $\left\vert AD\right\vert$, $\left\vert BC\right\vert$, and $\left\vert BD\right\vert$.

This equation can be simplified so that the potential energy can be determined as a function of four variables:

${\mathcal {V}}(r,\theta _{1},\theta _{2},\phi )=-{\frac {2\mu _{1}\mu _{2}}{4\pi \epsilon _{0}r^{3}}}\left(\cos \theta _{1}\cos \theta _{2}-{\frac {1}{2}}\sin \theta _{1}\sin \theta _{2}\cos \phi \right)$

Where μ₁ and μ₂ are the magnitude of the two dipole moments, r is the distance between the two dipoles and θ₁, θ₂, and Φ describe the orientation of the two dipoles. This potential energy decays at a rate of $r^{3}$.

The ideal attractive potential occurs when the two dipoles are aligned such that Φ = θ₁ = θ₂ = 0. The interaction potential energy can be rotationally averaged to provide the interaction energy thermally-averaged over all dipole orientations.

Quadrupole-Quadrupole Interactions

The Quadrupolar Phenomenon

Contour plot of the equipotential surfaces of an electric quadrupole field.

A molecular quadrupole moment arises from an uneven distribution of charge inside a molecule. Unlike a molecular dipole moment, quadrupole moments cannot be described using two point charges separated by a distance.

Quadrupoles are the 2nd-order term in the multipole expansion of an electric field of a molecule. The general form of a molecular quadrupole is a rank-two tensor (i.e., a $3\times 3$ matrix). For linear quadrupoles, we can describe the quadrupolar interaction by only looking at the $\Theta _{zz}$ moment.

In chemistry, some special cases arise when looking at linear molecules such as carbon dioxide (CO₂) and carbon disulfide (CS₂). Each of these molecules has a quadrupole moment, but they are of opposite sign. These phenomena can be explained by invoking Pauling electronegativity arguments. Oxygen has a greater electronegativity than carbon, so there is more electron density around the terminal oxygen atoms than in the centre. The opposite is true for carbon disulfide: the sp hybridized carbon is more electronegative than the terminal sulfur atoms and so electron density is localized around the centre of the molecule. The physical manifestation of the two different electronic distributions is that carbon dioxide has a quadrupole moment of $\Theta =-1.5\times 10$^-39 C·m², while that of carbon disulfide is $\Theta$= +1.2×10^-39 C·m²

Figure 1. Molecular structures showing the regions of partial charge in a) carbon dioxide and b) carbon disulfide

Calculating Quadrupole-Quadrupole Interactions

Figure 2. CO₂ molecules have quadrupole–quadrupole interactions with each other. The central carbon carries a partial positive charge, while the terminal oxygens carry partial negative charges. The like-charge atoms have repulsive electrostatic interactions (+/+ and -/-) while the oppositely-charged atoms have attractive electrostatic interactions (+/-).

The interaction energy of two linear quadrupoles is given by,

Linear Quadrupole–Quadrupole Interaction Energy ${\mathcal {V}}(r)={\frac {-6\,\Theta _{1}\Theta _{2}}{4\pi \epsilon _{0}r^{5}}}\times \Gamma (\theta _{1},\theta _{2},\phi )$

The quadrupole–quadrupole interaction energy is proportional to the reciprocal of the distance between the two particles ($r$) to the power of 5. Due to this higher exponent, quadrupole–quadrupole interactions become weaker as the intermolecular distance is increased at a faster rate than charge-charge interactions or charge-dipole interactions. The quadrupolar interaction energy also depends on the orientation of the two molecules with respect to each other. This is expressed in the orientational term, $\Gamma (\theta _{1},\theta _{2},\phi )$

$\Gamma (\theta _{1},\theta _{2},\phi )={\frac {1}{8}}[1-5\cos ^{2}\theta _{1}-5\cos ^{2}\theta _{2}-15\cos ^{2}\theta _{1}\cos ^{2}\theta _{2}+2(4\cos \theta _{1}\cos \theta _{2}-\sin \theta _{1}\sin \theta _{2}\cos \phi )^{2}]$

Applications of Quadrupole-Quadrupole Interactions to Molecular Organization

Quadrupole-quadrupole interactions fall under the category of electric phenomena. That is, they are governed by Coulomb's Law just like charge-charge and dipole-dipole interactions. Thus, the conventions of 'like-repels-like and 'opposites attract' are permitted in describing quadrupole-quadrupole interactions. See Figure 1. for a visualization.

Benzene (C₆H₆) has a quadrupole moment resulting from its π bonding system. The most stable geometries of a complex comprised of two benzene molecules are the 'T-shaped' and 'parallel displaced', where the quadrupole-quadrupole attractive interaction is strongest.

Figure 3. Schematics of a) the regions of partial charge in benzene, b) the parallel-displaced attractive orientation of two benzene molecules and c) the T-shaped attractive orientation of two benzene molecules

Hexafluorobezene (C₆F₆) is another aromatic molecule which possesses a quadrupole moment, although it is negative. A partial negative charge on the periphery of the ring is created by the strongly electron-withdrawing fluorine atoms, which cause the centre of the ring to be electron deficient. The favoured alignments of 2 hexafluorobenzene molecules are still 'T-shaped' and 'parallel displaced'. Because of the complementary differences in electronic distribution, benzene and hexafluorobezene can now stack in a parallel geometry.

Figure 4. Schematics of a) the regions of partial charge in hexfluorobenzene and b) the interaction orientation of hexafluorobenzene and benzene molecules

Appendix

Complete Equation for Charge-Linear Quadrupole Interaction Energy ${\mathcal {V}}(r)={\frac {q\mu }{4\pi \epsilon _{0}r^{3}}}(1+\cos(2\theta ))$

Complete Equation for Dipole-Linear Quadrupole Interaction Energy ${\mathcal {V}}(r)={\frac {3\,\mu _{1}\Theta _{2}}{4\pi \epsilon _{0}r^{4}}}\cdot {\frac {1}{2}}[\cos \theta _{1}(3\cos ^{2}\theta _{2}-1)-\sin \theta _{1}\sin 2\theta _{2}\cos \phi ]$

Complete Equation for Linear Quadrupole–Quadrupole Interaction Energy ${\mathcal {V}}(r)={\frac {6\,\Theta _{1}\Theta _{2}}{4\pi \epsilon _{0}r^{5}}}\cdot {\frac {1}{8}}[1-5\cos ^{2}\theta _{1}-5\cos ^{2}\theta _{2}-15\cos ^{2}\theta _{1}\cos ^{2}\theta _{2}+2(4\cos \theta _{1}\cos \theta _{2}-\sin \theta _{1}\sin \theta _{2}\cos \phi )^{2}]$

Induced Polarization

In isolation, atoms have spherical charge distributions with no dipole moment (top). Applied electric fields ($\mathrm {E}$) push electron density in the direction of the field, creating a depletion of charge on one side of the atom ($\delta +$)and an increase in charge on the other side ($\delta -$). This creates an induced dipole moment.

When a molecule is exposed to an electric field ɛ, electron density is polarized away from the nuclei. A dipole is thus created due to the asymmetric distribution of electron density. The induced molecular dipole created is proportional to the strength of the electric field (ɛ) as well as the polarizability of the molecule, α.

Polarizability allows us to better understand the interactions that occur between non polar atoms and molecules and other electrically charged species. In general, polarizability correlates with the interaction between electrons and the nucleus.

Units of Polarizability

Dimensional analysis of induced polarization.

$\mu _{ind}=\alpha \epsilon$

$\alpha ={\frac {\mu _{ind}}{\epsilon }}$

$\alpha ={\frac {Cm}{Vm^{-1}}}$

$\alpha =(\mathrm {C} \cdot \mathrm {m} ^{2}\cdot \mathrm {V} ^{-1})$

For convenience we will let:

$\alpha ^{'}={\frac {1}{4\pi \varepsilon _{0}}}\alpha$

$\alpha ^{'}={\frac {1}{CV^{-1}m^{-1}}}(\mathrm {C} \cdot \mathrm {m} ^{2}\cdot \mathrm {V} ^{-1})=m^{3}$

Polarizability has units of volume, which roughly correspond to the molecular volume. A large polarizability corresponds to a "soft" cloud of electron density, which is easily distorted by electric field.

Polarizability in Intermolecular Interactions

When we consider the effect that a charge (Q), has on a neutral atom at distance r, an induced dipole is created:

$\mu _{ind}=\alpha \epsilon =q\Delta r$

Where the electrostatic force is defined by:

$f=(charge)\times E$

Charge-induced dipole electrostatic force can be approximated by the electrostatic force between $Q$ and $q^{-}$ minus the force between $Q$ and $q^{+}$.

$f\approx q\Delta r{\frac {d\epsilon }{dr}}$

Since the electric field at distance r from a point charge Q:

$\epsilon ={\frac {Q}{4\pi \epsilon _{0}r^{2}}}$

Take the derivative

${\frac {d\epsilon }{dr}}=-2{\frac {Q}{4\pi \epsilon _{0}r^{3}}}$

$\epsilon {\frac {d\epsilon }{dr}}=\left({\frac {Q}{4\pi \epsilon _{0}r^{3}}}\right)\left(-2{\frac {Q}{4\pi \epsilon _{0}r^{3}}}\right)$

$\epsilon {\frac {d\epsilon }{dr}}=-2\left({\frac {Q}{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{5}}}$

Therefore:

$f=\alpha \epsilon {\frac {d\epsilon }{dr}}$

$f=-2\alpha \left({\frac {Q}{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{5}}}$

We can thus get the intermolecular potential by integrating this force over r:

${\mathcal {V}}(r)=-\int _{\infty }^{r}f(r)dr$

${\mathcal {V}}(r)=\int _{\infty }^{r}-2\alpha \left({\frac {Q}{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{5}}}dr=-{\frac {\alpha }{2}}\left({\frac {Q}{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{4}}}$

A schematic of how a (i) two neutral atoms can have an attractive electrostatic interaction when (ii) the electron density of one atom (left) spontaneously polarizes to form an instantaneous dipole. (iii) This polarizes the neighbouring atom.

Repulsive Interactions

Pauli Exclusion Principle

For a system of two indistinguishable fermions (particles of half integer spin: protons, electrons,...etc.), the total wavefunction is a mixture of independent single-particle states $\psi _{a}$ and $\psi _{b}$:

$\psi (\mathbf {r_{1}} ,\mathbf {r_{2}} )=A\left[\psi _{a}(\mathbf {r_{1}} )\psi _{b}(\mathbf {r_{2}} )-\psi _{b}(\mathbf {r_{1}} )\psi _{a}(\mathbf {r_{2}} )\right]$ , where $\mathbf {r}$ a spatial variable.

Notice that the wavefunction vanishes if $\psi _{a}=\psi _{b}$. In the case of an electron, the wavefunction also contains a spinor, $\chi (s)$, which distinguishes between spin up and spin down states. For simplicity, assume that both electrons have the same spin. This results in the Pauli Exclusion Principle:

Two electrons of the same spin cannot exist in the same state.

Since there are only two allowable spins (+1/2,-1/2), this guarantees that each orbital contains a maximum of two electrons of opposing spin.

Exchange Force

The apparent force that is observed when electron orbitals overlap originates from the requirement that the wavefunction be antisymmetric under exchange,

$\psi (\mathbf {r_{1}} ,\mathbf {r_{2}} )\chi (s_{1},s_{2})=-\psi (\mathbf {r_{2}} ,\mathbf {r_{1}} )\chi (s_{2},s_{1})$.

Electron Density Repulsion

This is commonly known as the exchange force. When the spin is symmetric (ex. $\uparrow \uparrow$, $\downarrow \downarrow$) the spatial wavefunction must be antisymmetric (antibonding), resulting in a repulsive interaction. This is also the case for filled electron orbitals; helium with an electron configuration $1s^{2}$ does not bond with itself because such interacting electrons reside in the destabilized ${\textstyle \sigma ^{*}}$ antibonding orbital (verified with MO-diagram of He–He).

The Pauli repulsion energy is of the form

${\mathcal {V}}(r)=ae^{-b(r-c)}$ , where $a$, $b$, $c$ are constants.

Notice that this parallels with the radial distribution of electron density around a nucleus (Figure: Electron Density Repulsion). Intuitively, the phenomenon of Pauli repulsion should only exist in the environment of electrons.

Comparison of Pauli Repulsion Formulae

It is of practical advantage to replace the exponential function with the simpler polynomial expression as in

${\mathcal {V}}(r)={\frac {C_{12}}{r^{12}}}$ , where $C_{12}$ is a constant.

A comparison of the polynomial and exponential form of the potential is illustrated on the right (Figure: Comparison of Pauli Repulsion Formulae). While the functions appear to be very different for small $\mathbf {r}$, the difference is immaterial considering that the probability of being in such a high energy regime is so low.

Comparison to Attractive Forces

Electrostatic Potential ${\mathcal {V}}\left(r\right)\propto {\frac {1}{r^{n+m+1}}}$

The potential for electrostatic interactions, such as charge-charge, charge-dipole, dipole-dipole, and quadrupole-quadrupole, always takes on some variation of the form above, where $n$ and $m$ are integers signifying the order of the pole in the multipole expansion (see Table: Range of Electrostatic Interactions).

Table: Range of Electrostatic Interactions

	Monopole n=0	Dipole n=1	Quadrupole n=2
Monopole m=0	$V(r)\propto {\frac {1}{r}}$	$V(r)\propto {\frac {1}{r^{2}}}$	$V(r)\propto {\frac {1}{r^{3}}}$
Dipole m=1		$V(r)\propto {\frac {1}{r^{3}}}$	$V(r)\propto {\frac {1}{r^{4}}}$
Quadrupole m=2			$V(r)\propto {\frac {1}{r^{5}}}$

For the sake of completeness also consider the charge-induced interactions

Attraction and Repulsion

${\mathcal {V}}(r)=-{\frac {\alpha }{2}}\left({\frac {Q}{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{4}}}\propto {\frac {1}{r^{4}}}$ ,

and and dipole-induced interactions

${\mathcal {V}}(r)=-{\frac {C_{6}}{r^{6}}}\propto {\frac {1}{r^{6}}}$ , $C_{6}={\frac {3}{2}}\alpha _{1}^{'}\alpha _{2}^{'}{\frac {I_{1}I_{2}}{I_{1}+I_{2}}}$, where $\alpha$ is the polarizability of the molecule, $I$ is the ionization potential, and $Q$ is the charge.

For electrostatic interactions following this general form, the potential energy for the interaction between oppositely-charged particles will approach negative infinity as the distance between approaches zero. This would suggest that if electrostatic interactions are the only forces in our model, oppositely-charged ions would collapse onto each other as $r$ approaches zero. Clearly this is an inaccurate representation of close range interactions, and some repulsive term must also be considered.

Pauli repulsion is a comparatively short range interaction, with $V(r)\propto 1/r^{12}$. The slope of this curve is very large for small internuclear separations, $r$. That is to say, the repulsive force rapidly diverges to infinity as the distance between nuclei approaches zero. For this reason it is said that the Pauli repulsive term is much "harder" than the other electrostatic forces.

Van der Waals Model of Fluid (Hard Sphere)

In liquids it is very improbable to have two particles at a highly repulsive distance (i.e., at distances where these functions diverge). Therefore, it is unnecessary to provide a perfect model at these distances, as long as it is repulsive. It is insightful to consider a fluid as collection of hard spheres, interacting with only Pauli repulsive forces. Neglecting the attractive forces between particles it justified by the following:

• Attractive forces are weak and isotropic. For a uniform (constant density) fluid these forces will average to zero.
• The cohesive force attributed to the change in entropy of gas to liquid phase is small for a dense gas; there is a small structural change between gas and liquid.
• Repulsive forces are extremely strong to prevent any arrangement where atoms overlap.
• Liquids arrange to minimize volume by maximizing packing of these hard spheres.  

See Also

• Lennard-Jones Potential
• Exchange_Interaction
• Pauli_exclusion_principle
• Multipole Expansion

Hydrogen Bonds

Hydrogen bonding is a type polar binding between a hydrogen atom that is bound to another atom with a large electronegativity, these then form electrostatic interactions with a hydrogen bond acceptor. Hydrogen bond acceptors fall into two main categories: atoms with a large negative charge such as Cl^- or Br^- and another compound with a polar bond, such as a ketone or carboxylic acid. These bonds are a subclass of dipole-dipole interactions. These interactions are found in many different scenarios, they are present in some organic solvents, solid crystal structures, proteins and DNA.

Hydrogen Bonding in Liquids

Hydrogen bonds are present in organic liquids like water, ethanol, hydrogen fluoride, and N-methylacetamide. These hydrogen bonds are partially responsible for the large enthalpies of vaporization, along with higher boiling points than hydrocarbons with comparable sizes. For example, the boiling point of hexane is 68 °C but the boiling point is hexanol is 135 °C. This is due to the hydrogen bonds that can be formed between the hydrogen atoms on the alcohols interacting with the oxygen atoms on other hexanol molecules.

The potential energy surfaces of the dimerization of hydrogen fluoride (black) and chloromethane (blue) calculated using MP2/aug-cc-pVTZ. Although the dipole moments of the two molecules are very similar, HF has a much more negative dimerization energy because the molecules can reach a shorter intermolecular distance before repulsive interactions become dominant.

Hydrogen Bonding in Solids

Crystal structure of the Ih phase of ice

Hydrogen bonds are also present in solids. The bonding in the solids is observed to allow for maximum packing of molecules that will give the solid the maximum number of hydrogen bonds leading to a stronger crystal. This can be observed in solid water, and is responsible for some of waters unique properties. Each water molecule is capable of donating two hydrogen bonds and accepting two hydrogen bonds.

Strength of Interactions

Magnitudes of hydrogen bonds are dependent on a couple different factors, the first is the difference in electronegativity of the atom that hydrogen is attached. The larger the difference in electronegativity the stronger the interaction will be. The second being the radius of the acceptor, it needs to be of a size that will allow for the hydrogen to get close enough to interact. Generally the Slater radius should be between 0.65 – 1.00 Å, any larger than this will not allow for the hydrogen to become close enough to the atom to attractively interact. The figure on the right illustrates the potential energies of the dimerization of HF and CH₃Cl, They both have very similar dipole moments, and but due to the hydrogen bonding the interactions between the HF molecules a larger potential energy minimum can be obtained, illustrated by the larger well depth. Also in this figure the well appears at a lower radius for the HF, and this is due to the previously mentioned dependence on the radius of the acceptor, the F is smaller than the Cl allowing for the hydrogen to get closer for stronger interactions. The directionality of the hydrogen bond also has effects on the strength of the bond, the optimal angle of interaction is 180°, while the optimal distance for interaction is between 1.5 – 2.5 Å. Planar molecules can contain very strong hydrogen bonding capabilities. The molecule being flat allows for constructive bond dipoles and it avoids steric clashes between other atoms in the molecule.

Dependence on Distance

The strength of hydrogen bonds dissipate rather quickly, they follow the same exponential trend that dipole dipole interactions follow, which is illustrated by the equation below:

Dipole-Dipole Interactions ${\mathcal {V}}(r,\theta _{1},\theta _{2},\Phi )={\frac {-\mu _{1}\mu _{2}}{4\pi \epsilon _{0}r^{3}}}(2\cos \theta _{1}\cos \theta _{2}-\sin \theta _{1}\sin \theta _{2}\cos \Phi )$

Where θ₁ , θ₂ and Φ are the two angles of orientation and the dihedral angles respectively. This equation illustrates that there is an inverse relationship between the strength of the bond and the distance between the dipoles cubed.

Effect of Radius ${\mathcal {V}}\propto {\frac {1}{r^{3}}}$

Rotational Averaging

Rotational Averaging

An ion interacts with a polar molecule. Vector geometry can be used to describe the energy of this interaction.

Rotational averaging describes the contribution to the potential energy from the rotational orientation of a charge-dipole interaction. Expectation values are utilized to give a single optimal value for the system's potential energy due to rotation.

For example, take a charged particle and a molecule with a permanent dipole. When they interact, the potential energy of this interaction can easily be calculated. For a dipole of length $l$, with a radius of $r$ between the dipole centre and the charged particle, the energy of interaction can be described by:

Potential Energy of a Charge Dipole Interaction ${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos {\theta }}{4\pi \epsilon _{0}\ r^{2}}}$

where $q$ is the charge of the particle, $\mu$ is the dipole moment, $\theta$ is the angle between $r$ and the dipole vector, $\epsilon _{0}$ is the vacuum permittivity constant, and $r$ is the radius between the particle and dipole.

Geometrically, this interaction is dependant on the radius and the length of dipole, as well as the orientation angle. If the radius $r$ between the ion and dipole is taken to be a fixed value, the angle $\theta$ still has the ability to change. This varied orientation of $\theta$ results in rotation of the dipole about its center, relative to the interacting charged particle. The weights of various orientations are described by a Boltzmann distribution expectation value, described generally by:

Expectation value (discrete states) $\langle M\rangle ={\frac {\sum _{i}M_{i}\exp \left({\frac {-{\mathcal {V}}_{i}}{k_{B}{\text{T}}}}\right)}{\sum _{i}\exp \left({\frac {-{\mathcal {V}}_{i}}{k_{B}{\text{T}}}}\right)}}$

Expectation value (continuous states) $\langle M\rangle ={\frac {\int M(r)\exp \left({\frac {-{\mathcal {V}}(r)}{k_{B}T}}\right){\textrm {d}}r}{\int \exp \left({\frac {-{\mathcal {V}}(r)}{k_{B}T}}\right){\textrm {d}}r}}$

where $\langle M\rangle$ is the expectation value, ${\mathcal {V}}$ is the energy value for a particular configuration, $k_{B}$ is Boltzmann's constant, and T is temperature. This Boltzmann-described weighting is the sum over the quantum mechanical energy levels of the system. Therefore, the probability $P$ is directly proportional to $\exp \left({\frac {-{\mathcal {V}}}{k_{B}{\text{T}}}}\right)$, indicating that at a specific temperature lower energy configurations are more probable. An equation can then be derived from this general expression, in order to relate it to the geometry and energy of a charge-dipole interaction.

Derivation of Rotationally-Averaged Charge-Dipole Interaction Potential

The orientationally averaged potential energy is the expectation value of the charge-dipole potential energy averaged over $\theta$

Starting with the potential energy of a charge-dipole interaction

${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos \theta }{4\pi \epsilon _{0}r^{2}}}$

We let $C=-{\frac {q_{ion}\mu }{4\pi \epsilon _{0}r^{2}}}$

This makes ${\mathcal {V}}(r,\theta )=-{\frac {q_{ion}\mu \cos \theta }{4\pi \epsilon _{0}r^{2}}}=C\cos \theta$

The average over the dipole orientation using the expectation value in classical statistical mechanics is:

$\langle {\mathcal {V}}\left(r\right)\rangle ={\frac {\int _{0}^{\pi }C\cos \theta \exp \left({\frac {-C\cos \theta }{k_{B}T}}\right)\sin \theta {\textrm {d}}\theta }{\int _{0}^{\pi }\exp \left({\frac {-C\cos \theta }{k_{B}T}}\right)\sin \theta {\textrm {d}}\theta }}$

Note: When integrating over an angle, the variable of integration becomes $\sin \theta {\textrm {d}}\theta$

To solve this integral we must first use first order Taylor's series approximation because integrals of exponential of $\cos \theta$ do not have analytical solutions.

$\int \exp \left({\frac {-C\cos \theta }{k_{B}T}}\right){\textrm {d}}\theta$

The first order Taylor's series approximations is as follows:

$f(x)\approx f(0)+x\cdot f'(0)$

$\exp(-x)\approx \exp(-0)-x\cdot \exp(-0)$

$\exp(-x)\approx 1-x$

Using the Taylor series with $\int \exp \left({\frac {-C\cos \theta }{k_{B}T}}\right){\textrm {d}}\theta$ gives: $\exp \left({\frac {C\cos \theta }{k_{B}T}}\right)\approx 1-{\frac {C\cos \theta }{k_{B}T}}$

The integral now becomes:

$\langle {\mathcal {V}}\left(r\right)\rangle ={\frac {\int _{0}^{\pi }C\cos \theta \left[1-{\frac {C\cos \theta }{k_{B}T}}\right]\sin \theta {\textrm {d}}\theta }{\int _{0}^{\pi }\left[1-{\frac {C\cos \theta }{k_{B}T}}\right]\sin \theta {\textrm {d}}\theta }}$

Multiplying into the brackets will give:

$\langle {\mathcal {V}}\rangle ={\frac {\int _{0}^{\pi }C\cos \theta \sin \theta -C\cos \theta \sin \theta {\frac {C\cos \theta }{k_{B}T}}{\textrm {d}}\theta }{\int _{0}^{\pi }\sin \theta -\sin \theta {\frac {C\cos \theta }{k_{B}T}}{\textrm {d}}\theta }}$

All terms that do not depend on ${\textrm {d}}\theta$ are constants and can be factored out of the integral. The terms can be expressed as 4 integrals:

$\langle {\mathcal {V}}\left(r\right)\rangle ={\frac {C\int _{0}^{\pi }\cos \theta \sin \theta {\textrm {d}}\theta -{\frac {C^{2}}{k_{B}T}}\int _{0}^{\pi }\cos ^{2}\theta \sin \theta {\textrm {d}}\theta }{\int _{0}^{\pi }\sin \theta {\textrm {d}}\theta -{\frac {C}{k_{B}T}}\int _{0}^{\pi }\sin \theta \cos \theta {\textrm {d}}\theta }}$

We must use trigonometric integrals to solve each of these for integrals:

First:

$C\int _{0}^{\pi }\cos \theta \sin \theta {\textrm {d}}\theta \longrightarrow C\int _{0}^{\pi }\cos \theta \sin \theta {\textrm {d}}\theta =C\cdot \int _{0}^{\pi }\sin \theta d(\sin \theta )=C\cdot \left[{\frac {1}{2}}\sin ^{2}\theta \right]{\Big |}_{0}^{\pi }=C\cdot 0=0$

Second:

${\frac {C^{2}}{k_{B}T}}\int _{0}^{\pi }\cos ^{2}\theta \sin \theta {\textrm {d}}\theta \longrightarrow {\frac {C^{2}}{k_{B}T}}\int _{0}^{\pi }\cos ^{2}\theta \sin \theta {\textrm {d}}\theta ={\frac {C^{2}}{k_{B}T}}\cdot -\int _{0}^{\pi }\cos ^{2}\theta d(\cos \theta )={\frac {C^{2}}{k_{B}T}}\left[-{\frac {1}{3}}\cos ^{3}\theta \right]{\Big |}_{0}^{\pi }={\frac {C^{2}}{k_{B}T}}\cdot {\frac {2}{3}}$

Third:

$\int _{0}^{\pi }\sin \theta {\textrm {d}}\theta \longrightarrow \int _{0}^{\pi }\sin \theta {\textrm {d}}\theta =-\cos \theta {\Big |}_{0}^{\pi }=2$

Fourth:

${\frac {C}{k_{B}T}}\int _{0}^{\pi }\sin \theta \cos \theta {\textrm {d}}\theta \longrightarrow {\frac {C}{k_{B}T}}\int _{0}^{\pi }\sin \theta \cos \theta {\textrm {d}}\theta ={\frac {C}{k_{B}T}}\int _{0}^{\pi }\sin \theta d(\sin \theta )={\frac {C}{k_{B}T}}\cdot \left[{\frac {1}{2}}\sin \theta \right]{\Big |}_{0}^{\pi }={\frac {C}{k_{B}T}}\cdot 0=0$

Plugging each trigonometric solved integral back into the equation gives:

$\langle {\mathcal {V}}\left(r\right)\rangle ={\frac {C\int _{0}^{\pi }\cos \theta \sin \theta {\textrm {d}}\theta -{\frac {C^{2}}{k_{B}T}}\int _{0}^{\pi }\cos ^{2}\theta \sin \theta {\textrm {d}}\theta }{\int _{0}^{\pi }\sin \theta {\textrm {d}}\theta -{\frac {C}{k_{B}T}}\int _{0}^{\pi }\sin \theta \cos \theta {\textrm {d}}\theta }}={\frac {0-{\frac {2C^{2}}{3k_{B}T}}}{2-0}}=-{\frac {C^{2}}{3k_{B}T}}$

Finally replace $C$ with $C=-{\frac {q_{ion}\mu }{4\pi \epsilon _{0}r^{2}}}$ to give:

The Orientational Average of the Charge-Dipole Interaction Energy $\langle {\mathcal {V}}\left(r\right)\rangle =-{\frac {1}{3k_{B}T}}\left({\frac {q_{ion}\mu }{4\pi \epsilon _{0}}}\right)^{2}{\frac {1}{r^{4}}}$

The Lennard-Jones Potential

The Lennard-Jones 6-12 potential approximates the intermolecular interactions of two atoms due to Pauli repulsion and London dispersion attraction. The potential is defined in terms of the well-depth ($\epsilon$) and the intercept ($\sigma$). Other formulations use the radius where the minimum occurs, $R_{min}$, instead of $\sigma$.

The Lennard-Jones potential describes the interactions of two neutral particles using a relatively simple mathematical model. Two neutral molecules feel both attractive and repulsive forces based on their relative proximity and polarizability. The sum of these forces gives rise to the Lennard-Jones potential, as seen below:

${\mathcal {V}}\left(r\right)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]=\varepsilon \left[\left({\frac {R_{min}}{r}}\right)^{12}-2\left({\frac {R_{min}}{r}}\right)^{6}\right]$

where ε is the potential well depth, σ is the distance where the potential equals zero (also double the van der Waals radius of the atom), and R_min is the distance where the potential reaches a minimum, i.e. the equilibrium position of the two particles.

The relationship between the $\sigma$ and $R_{\min }$ is $R_{min}={\sqrt[{6}]{2}}\sigma$

Pauli Repulsion

The first half of the Lennard-Jones potential is Pauli-Repulsion. This occurs when two closed shell atoms of molecules come in close proximity to each other and their electron density distributions overlap. This causes high inter-electron repulsion and extremely short distances, inter-nuclei repulsion. This repulsion follows an exponential distribution of electron density:

${\mathcal {V}}_{rep}\left(r\right)=Ae^{-cr}$

where A and c are constants and r is the intermolecular distance. In a liquid, however, it is very unlikely that two particles will be at highly repulsive distances, and thus a simplified expression can be used by assuming the potential has a r^-12 dependence (note that this high exponent means that the energy of repulsion drops off extremely fast as molecules separate). The resulting simple polynomial is as follows:

${\mathcal {V}}\left(r\right)={\frac {C_{12}}{r^{12}}}$

Where the C₁₂ coefficient is defined as:

$C_{12}=4\varepsilon \sigma ^{12}=\varepsilon R_{min}^{12}$

$\Phi _{12}(r)=A\exp \left(-Br\right)-{\frac {C}{r^{6}}}$

Buckingham

Further Reading

• Buckingham potential
• Exchange interaction

London Dispersion

The Second half of the Lennard-Jones potential is known as London dispersion, or induced Dipole-Dipole interaction. While a particular molecule may not normally have a dipole moment, at any one instant in time its electrons may be asymmetrically distributed, giving an instantaneous dipole. The strength of these instantaneous dipoles and thus the strength of the attractive force depends on the polarizability and ionization potential of the molecules. The ionization potential measures how strongly outer electrons are held to the atoms. The more polarizable the molecule, the more its electron density can distort, creating larger instantaneous dipoles. Much like with Pauli Repulsion, this force is dependent on a coefficient, C₆, and also decays as the molecules move further apart. In this case, the dependence is r^-6

${\mathcal {V}}\left(r\right)={\frac {-C_{6}}{r^{6}}}$

$C_{6}={\frac {3}{2}}\alpha _{1}^{'}\alpha _{2}^{'}{\frac {I_{1}I_{2}}{I_{1}+I_{2}}}$

Where α' is the polarizability, usually given as a volume and I is the ionization energy, usually given as electron volts. Lastly, the C₆ constant can also be expressed by the variables as seen in the Lennard-Jones equation:

$C_{6}=4\varepsilon \sigma ^{6}=2\varepsilon R_{min}^{6}$

Combination Rules

The estimation of Lennard-Jones potential parameters for mixed pairs of atoms (i.e., AB) using the Lorentz-Berthelot combination rule.

In the case of two separate molecules interacting, a combination rule called the Lorentz-Berthelot combination rule can be imposed to create new σ and ε values. These values are arithmetic and geometric means respectively. For example, an Ar-Xe L-J plot will have intermediate σ and ε values between Ar-Ar and Xe-Xe. An example of this combination rule can be seen in the figure to the right.

$\sigma _{AB}={\frac {\sigma _{AA}+\sigma _{BB}}{2}}$

$\varepsilon _{AB}={\sqrt {\varepsilon _{AA}\varepsilon _{BB}}}$

Example

Calculate the intermolecular potential between two Argon (Ar) atoms separated by a distance of 4.0 Å (use ϵ=0.997 kJ/mol and σ=3.40 Å).

${\mathcal {V}}\left(r\right)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]$

${\mathcal {V}}\left(r\right)=4(0.997~{\text{kJ/mol}})\left[\left({\frac {3.40}{4.00}}\right)^{12}-\left({\frac {3.40}{4.00}}\right)^{6}\right]$

${\mathcal {V}}\left(r\right)=3.988(0.142242-0.377150)$

${\mathcal {V}}\left(r\right)=-0.94~{\text{kJ/mol}}$

See Also

Wikipedia:Lennard-Jones potential Chemwiki: Lennard-Jones

Intramolecular Forces

The representation of the potential energy of a torsional angle rotation using a cosine representation.

Bond Stretches

Bond stretches are the alteration of the bond length within a molecule. Using water as an example we will explain the different forms of stretches. The two modes of stretching are symmetric; to which in the case of water the two oxygen hydrogen bonds lengthen and shorten at the same rate and asymmetric; to which they lengthen and shorten at opposite times. The energy associated with this motion can be described using the following equation:

${\mathcal {V}}_{bond}(r)={\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}$

We see a Hooke's law relationship where k is a constant for the bond length and the x² term in the Hooke's law equation is replaced with the change in bond length subtracted from the equilibrium bond length. This will give an average change in bond length. This relationship will give us an accurate value for the energy of the bond length of the given molecule.

Bond Angle Bending

Bond angle energies are brought about from an alteration of the optimized bond angle to a less favorable conformation. This can be either a decrease in angle; the orbitals get closer to one another which will increase the potential energy, or increase in angle; the orbitals more away from the ideal overlap which give the molecule favorable orbital interactions to increase the stability. This will increase the potential energy by moving the molecule out of the low energy conformation. This energy relationship can be observed through the following equation:

${\mathcal {V}}_{angle}(\theta )={\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}$

We again see a hooks law relationship where k is a constant for the bond length and the x² term in the Hooks law equation is replaced with the change in bond angle subtracted from the equilibrium bond angle. This will give an average change in bond angle. This relationship will give us an accurate value for the energy of the bond length of the given molecule.

Torsional Rotations

The torsional potential energy surface of butane.

Torsional rotations, which can also be described as dihedral interactions come about when a system of 4 or more atoms form a molecule. This interaction relies upon the rotation central bond into a favorable orientation with the most stable overlap conformation. Some of these known conformations are called periplanar, anticlinal, gauche, and antiperiplaner, listed in order of highest to lowest energy. Observing Newman projections upon the energy we can see where the unfavorable confirmational interactions will lie. The methyl groups here will be considered large and have unfavorable overlap with one another. Mathematically this can be represented using the following equation:

${\mathcal {V}}_{dihedral}(\Phi )={k_{\Phi }}({1+\cos(n\Phi -\delta }))$

Where k is the rotational barrier of the system, n is the periodicity as the rotations repeat around 360 degrees, ɸ is the dihedral angle, and δ is the offset of the function

Improper Torsion

An improper torsion is a special type of torsional potential that is designed to enforce the planarity of a set of atoms (ijkl) where i is the index of the central atom. The torsional potential is defined as the angle between the planes defined by atoms ijk and jkl.

Internal Non-Bonded Interactions

The atoms of a molecule will also have intramolecular steric, dispersive, and electrostatic interactions with each other. These interactions within a molecule determine its conformation. The steric repulsion as well as Coulombic and dispersive interactions all together in the conformation of the molecule. Every atom in the system can be assigned a partial charge ($q_{i}$). These interactions are calculated using standard Lennard-Jones and Coulombic pairwise interaction terms,

${\mathcal {V}}_{NB}(r_{ij})={4\epsilon _{ij}}\left[\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right]+{\frac {q_{1}q_{2}}{4\pi \epsilon _{0}}}{\frac {1}{r_{ij}}}$

Here we see the combination of the Leonard Jones equation as well as Coulomb's law. With this relationship we can calculate the force resulting between the non bonding intermolecular forces. It should be noted at this point that interactions between bonded or atoms forming angles are excluded from this relationship. These interactions have been included though the calculations in the previous sections.

General Form of Potential Energy Function

The overall energy of a system is based on an enormous number of factors, and it is the sum of these individual attractive and repulsive forces between all the atoms in the system that results in a final "force field". The potential energy for a set of atomic coordinates (r) is the sum of bonded (ie intermolecular forces) and non-bonded (i.e., repulsion/dispersion and electrostatics) terms. This sum ultimately defines the molecular mechanical force field.

${\mathcal {V}}({\textrm {r}})=\sum _{bonds}{\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}+\sum _{angles}{\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}+\sum _{dihedrals}{k_{\phi }}({1+\cos(n\phi -\delta }))+\sum _{pairs:i,j}{4\epsilon _{ij}}\left[\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right]+{\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$

Due to this complexity, calculating the potential energy for a set of molecules using a force field would require a complete set of known parameters, which also require several different methods of determination. These methods are outlined in this section:

Interaction	Term type	Expression	Parameters
Bonds	Bonded	${\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}$	spring constant (kbond), equilibrium length (re)
Angles	Bonded	${\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}$	spring constant (kangle), equilibrium angle (θe)
Dihedrals	Bonded	${k_{\phi }}({1+\cos(n\phi -\delta }))$	barrier height (kφ), periodicity (n), offset (δ)
Repulsion/Dispersion	Non-bonded	${4\epsilon _{ij}}\left[\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right]$	atomic radii (σ), well depth (ε)
Electrostatics	Non-bonded	${\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$	partial atomic charges (q)

Parameterization

Empirical Parameters

For small molecules, parameters may be defined to fit experimental properties. Liquids will have different bulk physical properties for each set of non-bonded parameters. Each set of parameters in the force field will result in a different set of physical properties, the combinations of which can be tested until properly representing the properties of the liquid. An example would be determining the proper model for bulk water assuming no Lennard-Jones interactions for hydrogen atoms. Parameters would be tested together until the model properly represented the expected characteristics of bulk water.

Calculating Electrostatic Interactions

Starting with the non-bonded terms, we can look at electrostatics. Every intermolecular pair of atoms has this electrostatic interaction, and can be summed up in a form comparing all atoms in one molecule (A) to all atoms in another molecule (B). This quickly adds up to a massive sum even when considering systems of a few complex molecules:

${\mathcal {V}}(r_{ij})=\sum _{i}^{N_{A}}\sum _{j}^{N_{B}}{\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$

In regards to partial charges, as there is no rigorous way to translate charge distribution to point charges, certain constraints are formed: sum of partial charges should equal formal charge of molecule (seen below), and equivalent atoms should carry the same charge. Charges can be chosen to match gas phase dipole moment, be treated as free parameters and adjusted for target properties, and to fit quantum mechanical electrostatics:

${\vec {\mu }}_{exptl}=\sum _{i}^{N_{atom}}q_{i}{\vec {r}}_{i}$

Charge fitting schemes fit partial charges to the quantum mechanical electrostatic potential (the potential energy a unit point charge has at position in space r). Potential can be calculated at a set of points around the molecule to generate electrostatic potential surface (ie: can map out the electronic density around the molecule and determine the localized area of charges):

Target Electrostatic Potential Data from Quantum Chemical Calculations

First term represents electron density, second term represents nuclear charges at positions R_i:

${\mathcal {V}}_{QM}(r)=\int {\frac {\rho (r')}{|r-r'|}}dr'+\sum _{i}^{N}{\frac {Z_{i}}{|r-R_{i}|}}$

Molecular Mechanical Partial Atomic Charges

Looks at partial atomic charges at positions R_i:

${\mathcal {V}}_{MM}(r)=\sum _{i}^{N}{\frac {q_{i}}{|r-R_{i}|}}$

In terms of electrostatic charges, point charges (q) must be defined as such that $V_{QM}(r)$ ≈ $V_{MM}(r)$ as closely as possible, thereby fitting both models.

Restraining terms can be introduced to give more practical models of molecules (known as Restrained Electronic Potential (RESP) fitting). Such restraints limits charges to moderate and chemically relevant values. Restraints are also placed to ensure chemically equivalent groups retain the same charges. These groups are given an atom type distinguished by element, hybridization and bonding partners, and parameters are defined for each type. For example, all aromatic carbons part of the same benzene ring should have the same parameters.

Generalized Amber Force Field (GAFF)

GAFF is a commonly used force field designed for reasonable parameters for simulations of small organic molecules. GAFF involves the use of the RESP method for assigning partial charges, and atom type defines the Lennard-Jones parameters and internal bonding terms that will be associated with the atom. Atom type definitions can be extremely specific (e.g., H on aliphatic C with 3 electron-withdrawing groups; inner sp2 C in conjugated ring systems).

Lennard-Jones Potential

A combination of the Pauli Repulsion and London dispersion attraction terms, the Lennard-Jones Potential is defined using two parameters: van der Waals radius (σ) and well depth (ε). These parameters are responsible for packing of molecules in liquids, and must compensate for errors in all other properties. van der Waals radius essentially sets the limit at which molecules may approach one another (potential energy > 0 when closer than σ). Well depth sets the energy well of the potential energy, occurring when distance between molecules is at R_min (ie radius related to minimum potential energy). First term represents Pauli repulsion, second term represents London dispersion attraction:

${\mathcal {V}}(r)={\frac {C_{12}}{r^{12}}}-{\frac {C_{6}}{r^{6}}}$ $={4\epsilon }\left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]$ $={\epsilon }\left[\left({\frac {R_{min}}{r}}\right)^{12}-2\left({\frac {R_{min}}{r}}\right)^{6}\right]$

Bonding Interactions

Bond Stretch

Atoms within a molecule may vary their distances (stretch or compress), and this relationship can be described as a spring. This is why Hooke's Law can be applied to bond lengths to define the energy required to stretch/compress, according to spring constant (k_bond) and equilibrium length (r_e). Spring constant may be determined via the infrared vibrational frequency of the bond stretch (according to $\nu ={\frac {1}{2\pi }}\left({\frac {k}{\mu }}\right)^{\frac {1}{2}}$, while equilibrium bond lengths can be found through the use of rotational spectroscopy and crystallography.

${\mathcal {V}}_{bond}={\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}$

Angle Stretch

The potential energy of bending a bond angle can also be approximated by Hooke's law. Hooke's Law can again be applied to bond angles to define the energy required to stretch/compress, according to spring constant (k_angle) and equilibrium angle (θ_e). Again, the spring constant may be determined via the infrared vibrational frequency of the bond stretch, following the same formula, while the equilibrium angle can be found through the use of rotational spectroscopy and crystallography.

${\mathcal {V}}_{angle}={\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}$

Dihedral Parameters

Also known as torsional rotations, dihedral relationships are produced with 4 atoms. This relates to how the central bond is rotated and the overlap of the outer two atoms (designated as 1,4 interactions, with 2 and 3 being the central atoms). These overlaps may result in a higher or lower potential energy depending on confirmation (i.e. gauche, trans, periplanar, etc), and the relationship is dependent on barrier height (k_φ) (energy required to rotate the central axis), periodicity (n) and offset (δ). Determination of these parameters can be done using quantum chemistry or experimental isomer populations.

${\mathcal {V}}_{dihedral}={k_{\phi }}({1+\cos(n\phi -\delta }))$

Calculation of molecular forces using quantum chemistry

Intermolecular Interactions

Intermolecular interactions can be calculated using quantum chemical programs.

$\Delta E_{interaction}=E_{complex}-\sum _{i}E_{fragment,i}$

The fragments are the atoms or molecules that comprise the complex.

Basis Set Superposition Error

Many quantum chemical programs use atom centered basis functions to represent the atomic orbitals. For a given calculation, the number of basis functions used is fixed and finite. Inevitably, the description of the orbitals will be incomplete, in that a more accurate answer could be obtained if more basis functions were added.

When calculating intermolecular interactions, the calculation of the $\Delta E_{interaction}$ will have the basis functions from all the fragments. At close intermolecular contacts, the tails of these orbitals will extend close to the space of atoms in other fragments. In this way, the fragments will be stabilized by the presence of these additional basis functions, so the energy of the complex includes both the stabilization due to the intermolecular interaction and this "basis set superposition error."

Statistical properties

Statistical thermodynamics describes physical descriptions according to probability distributions.

Probability Distributions

A probability distribution is a function that shows the likelihood of an outcome. In statistical mechanics, this generally means the probability of of a system being in a particular state.

Gaussian Distribution

A Gaussian distribution, also called a normal distribution is a function with a bell shaped distribution. Any random variable with this normal probability density

$p(x)={\frac {1}{\sqrt {2\pi \sigma ^{2}}}}\;e^{-{\frac {(x-\alpha )^{2}}{2\sigma ^{2}}}}$

is called a normal variable with mean $\alpha$ and variance $\sigma ^{2}$

The Boltzmann Distribution is one such probability distribution that gives the probability distribution as a function of a states energy $E$ and temperature if a system $T$.

$p_{i}={\frac {e^{-{\varepsilon }_{i}/kT}}{\sum _{j=1}^{M}{e^{-{\varepsilon }_{j}/kT}}}}$

As shown below physical properties of a system can be calculated using this Boltzmann distribution.

Conformational Distribution

Within solids, liquids, and gases, atoms can take on different orientations and arrangements. Conformation changes within a molecule can occur by rotations around bonds. For example, gauche or eclipsed conformations, trans or cis and any conformation in between these would all fall under the distribution of conformation.

Calculating Physical Properties

Macroscopic properties are the average arrangement a system take son over time. Assuming all possible energy levels are continuous, we can employ classical mechanics to calculate physical properties. This assumption is only valid when particles and heavy and forces are relatively soft.

Averages

The ensemble average in statistical mechanics refers to the mean of all possible states of a system as given by its probability distribution. It is dependant on the ensemble chosen (i.e canonical, microcanonical etc...).

The expectation value $\langle M\rangle$ gives the average value when any physical property of a system is measured. depending on the type of distribution, it may not be the most probable value but is the probability weighted value which we expect to measure. In a classical system, the expectation value is an integral over all of the possible configurations. The possible configurations are integrated over the interval $[x_{i},x_{f}]$.

$\langle M\rangle ={\frac {\int _{x_{i}}^{x_{f}}M(x)\exp \left({\frac {-{\mathcal {V}}(x)}{k_{B}T}}\right)\,dx}{\int _{x_{i}}^{x_{f}}\exp \left({\frac {-{\mathcal {V}}(x)}{k_{B}T}}\right)\,dx}}$

where ${M}(x)$ represents the property being calculated, ${\mathcal {V}}(x)$ represents the potential energy, $k_{B}$ represents the Boltzmann constant, and ${T}$ represents the temperature of the system.

In classical statistical mechanics, the classical ensemble average is the normalized Boltzmann weighted integral over all phase space.

$\langle M\rangle ={\frac {{\int _{-\infty }^{\infty }}{\int _{-\infty }^{\infty }}{\int _{-\infty }^{\infty }}{\int _{0}^{L}}{\int _{0}^{L}}{\int _{0}^{L}}{\mathcal {H}}(x,y,z,p_{x},p_{y},p_{z})e^{\frac {-{\mathcal {H}}(x,y,z,p_{x},p_{y},p_{z})}{k_{B}T}}dxdydzdp_{x}dp_{y}dp_{z}}{{\int _{-\infty }^{\infty }}{\int _{-\infty }^{\infty }}{\int _{-\infty }^{\infty }}{\int _{0}^{L}}{\int _{0}^{L}}{\int _{0}^{L}}e^{\frac {-{\mathcal {H}}(x,y,z,p_{x},p_{y},p_{z})}{k_{B}T}}dxdydzdp_{x}dp_{y}dp_{z}}}$

Where ${\hat {H}}$ is the Hamiltonian for the described system. This expression can be used to find many physical properties such as the average energy of a single particle system and many-body systems by integrating over spatial coordinates ($dxdydz$ and integrating over the Maxwell distribution ($dp_{x}dp_{y}dp_{z}$

In a quantum mechanical system, the expectation value is the Boltzmann weighted sum over energy levels.

$\langle M\rangle ={\frac {\sum _{i}g_{i}M_{i}e^{\frac {-E_{j}}{k_{B}T}}}{\sum _{j}g_{j}e^{\frac {-E_{j}}{k_{B}T}}}}$

Example: Conformationally Averaged Dipole Moment

Conformational averaging using the Boltzmann distribution gives a way of finding the average dipole moment, as described in previous section.

The dipole moment of a molecule changes as its conformation changes as is it the vector sum of all bond moment and therefore the observed (or 'expected') value is a linear combination of the two conformations. For example, in 1,2-dichloroethane the trans conformation is non-polar but the cis conformation is polar. Since the molecular dipole is a vector quantity, the conformationally-averaged dipole moment is the average of the square of the individual dipole moments.

$\mu _{average}={\sqrt {\mu ^{2}}}={\sqrt {\frac {\sum _{i}g_{i}\mu ^{2}e^{\frac {-\nu }{k_{B}T}}}{\sum _{j}g_{j}e^{\frac {-\nu }{k_{B}T}}}}}$

${\sqrt {\mu ^{2}}}={\sqrt {\frac {\mu _{trans}^{2}e^{\frac {-\nu _{trans}}{k_{B}T}}+\mu _{cis+}^{2}e^{\frac {-\nu _{cis}}{k_{B}T}}+\mu _{cis-}^{2}e^{\frac {-\nu _{cis}}{k_{B}T}}}{e^{\frac {-\nu _{trans}}{k_{B}T}}+e^{\frac {-\nu _{cis}}{k_{B}T}}+e^{\frac {-\nu _{cis}}{k_{B}T}}}}}$

Using this formula, the conformationally averaged dipole moment for a molecule like 1,2 dichloroethane where the trans conformation has no dipole and the cis does can be calculated.

Variance

Variance is value that indicates how widely spread a set of data is when compared to the average. The variance is the expectation of the square of the standard deviation of a given set of data. If the standard deviation or variance is low, the data is close to the expected value. The variation of ${X}$, a random variable, is the expected value of the square of the standard deviation from the mean. ${E}$ represents the expectation and $\mu$ represents the mean. Variance can be represented by $\operatorname {Var} (X)$ or $\sigma ^{2}$.

$\operatorname {Var} (X)=\operatorname {E} \left[(X-\mu )^{2}\right]$

Thermodynamic ensembles

The entire universe is governed by the laws of thermodynamics through the transfer of energy between matter. This process is too complex to consider directly, instead we consider several simplifications where part of the universe (i.e., the system) is considered separately from the rest of its surroundings. The system houses the process or thermodynamic state under consideration, whereas the surroundings contain everything else.

The system can be described using an ensemble, also known as a statistical ensemble, which is an idealization consisting a large number of copies of a system, considered all at once, each of which represents a possible state that the real system might be in. A thermodynamic ensemble is a statistical ensemble that is in statistical equilibrium. It provides a way to derive the properties of a real thermodynamic systems from the laws of classical and quantum mechanics.

We can consider different systems with different degrees of separation from their surroundings, from completely isolated (e.g., microcanonical ensemble) to completely open to its surroundings (e.g., grand canonical ensemble). However, for the purposes of molecular simulations only three main ensembles will be considered. Namely, the microcanonical, canonical, and isothermal-isobaric ensembles

Microcanonical ensemble

The microcanonical ensemble is a system that is completely isolated from its surroundings. There is no transfer of energy or matter between the system and the surroundings, and the volume of the system remains fixed. In other words, this is a physical system where the total number of particles (N), the total volume (V), and the total energy (E) are constant. The microcanonical ensemble is typically abbreviated NVE.

Constant number of particles (N), volume (V), and total energy (E).

The total energy of the system is constant, but there is no defined temperature of the microcanonical ensemble. Temperature can only be defined for systems interacting with their surroundings. Isolated systems, such as microcanonical, are only described in terms of their composition, volume, and total energy. In statistical mechanics, a microcanonical ensemble is used to represent the possible states of a system which have a specified total energy.

Canonical ensemble

In the canonical ensemble, the volume of the system is fixed, and energy can transfer across the boundary between system and surroundings, but matter cannot. We can describe systems as being immersed in a heat bath at a temperature (T), where the heat bath is much larger than the system. No amount of heat put off by the system will significantly increase the temperature of the surroundings. The canonical ensemble applies to systems of any size; while it is necessary to assume that the relative size of the heat bath is very large, the system itself may be small or large. Now because the system and surroundings are in thermal contact, the system will transfer heat (q) to and from the surroundings until they are in thermal equilibrium. Therefore, unlike the microcanonical ensemble, the temperature of the canonical ensemble can be a defined constant (T). This ensemble is typically abbreviated in terms of constant number of particles, total volume, and temperature (NVT).

Constant number of particles (N), volume (V), and temperature (T).

The canonical ensemble is a statistical mechanics ensemble that represents the possible states of a mechanical system in thermal equilibrium with a heat bath at some fixed temperature. The system can exchange energy with the heat bath, so that the states of the system will differ only in total energy. In the microcanonical ensemble, the total energy is fixed, but in the canonical ensemble the energy is no longer constant. It can take on a range of values depending on the temperature. The canonical ensemble is important when attempting to describe the Helmholtz free energy of a system, which is the maximum amount of work a system can do at a constant volume (V) and temperature (T).

Isothermal-Isobaric ensemble

In the isothermal-isobaric ensemble, energy can transfer across the boundary, but matter cannot. The volume of the system can change such that the internal pressure of the system matches the pressure exerted on the system by its surroundings. Similar to the canonical ensemble, we can describe the isothermal-isobaric ensemble as a system immersed in a heat bath at a temperature (T), where the heat bath is much larger than the system. No amount of heat put off by the system will significantly increase the temperature of the surroundings.

The isothermal-isobaric ensemble is a statistical mechanical ensemble that maintains a constant total number of particles, and constant temperature (T) and pressure (p), typically abbreviated NpT. This ensemble plays an important role in chemistry since the majority of important chemical reactions are carried out under constant pressure conditions. The isothermal-isobaric ensemble is also important when attempting to describe the Gibbs free energy of a system, which is the maximum amount of work a system can do at constant pressure (p) and temperature (T).

Constant number of particles (N), pressure (p), and temperature (T).

The Hamiltonian

The Hamiltonian, ${\mathcal {H}}$, is a function that calculates the sum of the potential and kinetic energy of a system. For the systems considered here, the Hamiltonian is a function of the positions (r) and momenta (p) of the particles in the system. The potential energy, ${\mathcal {V}}$, depends only on the positions of the atoms $({\textbf {r}})$ and the kinetic energy, ${\mathcal {T}}$, depends only on the momentum of the particles $({\textbf {p}})$ [1]. Therefore, the Hamiltonian can be separated into a potential energy term that depends only on the positions of the particles (${\mathcal {V}}({\textbf {r}})$) and a kinetic energy term that depends only on the momenta of the particles (${\mathcal {T}}({\textbf {p}})$),

${\mathcal {H}}({\textbf {r}},{\textbf {p}})={\mathcal {V}}({\textbf {r}})+{\mathcal {T}}({\textbf {p}})$

The microcanonical ensemble is the statistical ensemble used to represent the isolated system, where the number of particles (N), the volume (V), and the total energy (E) are constant.

Phase Space

The unique positions and momenta of each particle, for a system of N particles, can be accurately described by a 6 N dimensional spaced called phase space. Each particle in a system has three unique position variables $(x,y,z)$ ranging from 0 to L and three unique momentum variables $(p_{x},p_{y},p_{z})$ ranging from $-\infty$ to $\infty$ [1]. These 6N numbers constitute the microscopic state of the system at time t [2,3]. A system of N particles corresponds to a phase space with 6 N variables, which describes every conceivable position and momentum combination of all particles in a mechanical system [1,2,3].

Phase Point

A phase point refers to any one point in an N-body system at any time t [2]. The dynamics of the phase point are fully described by the motion and trajectory of the phase point as it travels through phase space [2]. Every classical state can be considered to be equally probable as N, V, and E are held constant [3].

Two Argon atoms inside phase space. The blue and green arrow vectors correspond to the atoms position and momenta, respectively. The atoms are non-interacting, however they do change trajectory as they bounce off the walls of the box.

Example

Consider a system containing two particles of Argon (labelled 1 and 2, respectively), where N = 2. Since phase space is described by 6 N dimensions the total number of spatial variables becomes 6(2) = 12.

• r₁, position of particle 1: $x_{1},y_{1},z_{1}$
• r₂, position of particle 2: $x_{2},y_{2},z_{2}$
• p₁, momentum of particle 1: $p_{x1},p_{y1},p_{z1}$
• p₂, momentum of particle 2: $p_{x2},p_{y2},p_{z2}$

Argon is a noble gas and inert to intermolecular interactions with other particles in a system. Accordingly, the two Argon atoms in this hypothetical system are non-interacting. When atoms are assumed to be non-interacting, potential energy interactions between particles are no longer accounted for and ${\mathcal {V}}$ = 0.

The Hamiltonian becomes the sum of the kinetic energies of the two particles:

${\mathcal {H}}={\frac {p_{1}^{2}}{2m}}+{\frac {p_{2}^{2}}{2m}}$

The Classical Partition Function

In quantum mechanical systems, the system must occupy a set of discrete (quantized) states. These states correspond to the energy levels of the system, labelled $E_{1},E_{2},...,E_{n}$ A partition function can be defined a sum over these states.

$Q=\sum _{i}\exp \left({\frac {-E_{i}}{k_{B}T}}\right)$

Where $E_{i}$ is the energy of state $i$, $k_{B}$ is the Boltzmann constant and $T$ is the temperature.

Explicit sums are practical when the energy states of molecules are independent and there are no intermolecular interactions. In systems where intermolecular interactions are significant, it is generally more practical to approximate the system as existing in terms of continuous variables. The states are approximated to be continuous and can be integrated with respect to the phase space of the system while simultaneously sampling every state of the system. This leads to the classical partition function:

$Q_{classical}=\int \dots \iint \dots \int \dots \iint \exp \left({\frac {-{\mathcal {H}}({\textbf {r}},{\textbf {p}})}{k_{B}T}}\right)\,dr_{1}\,dr_{2}\dots \,dr_{N}\,dp_{1}\,dp_{2}\dots \,dp_{N}$

Many Body Systems

The expectation values of a classical system can be evaluated by integrating over phase space. The classical partition function is an integral over 6N phase space variables:

${\textrm {d}}{\textbf {r}}=dx_{1}\,dy_{1}\,dz_{1},dx_{2},dy_{2},dz_{2},...,dx_{N}\,dy_{N}\,dz_{N}$

${\textrm {d}}{\textbf {p}}=dp_{x,1}\,dp_{y,1},dz_{z,1},...,dp_{x,N}\,dp_{y,N}\,dp_{z,N}$$Q_{classical}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{0}^{L}\int _{0}^{L}\int _{0}^{L}\exp \left({\frac {-{\mathcal {H}}({\textbf {r}},{\textbf {p}})}{k_{B}T}}\right){\textrm {d}}{\textbf {r}}{\textrm {d}}{\textbf {p}}$

References

1. McQuarrie, Donald A. 1973. Statistical Thermodynamics. New York, N.Y. Harper & Row, Publishers, Inc. Pages 117-118.
2. Tuckerman, Mark E. 2010. Statistical Mechanics: Theory and Molecular Simulation. Oxford England. New York: Oxford University Press. Pages 2-5.

Dilute gases

A molecular dynamics simulation of high-density argon gas at 300 K. Argon-argon interactions are described using a Lennard-Jones potential. The gas density is 5282 mol per cubic meter. The pressure is approximately 13200 kPa.

Ideal Gas Law

Gasses are often described using ideal gas law, which relates the pressure of a gas to its density with a simple expression.

${\frac {P}{k_{B}T}}=\rho$

Where $\rho$ is the density of the gas, $P$ is the pressure of the gas, $k_{B}$ is the Boltzmann constant, equal to $1.38064852\times 10^{-23}m^{2}kgs^{-2}K^{-1}$, and $T$ is the temperature in Kelvin. The expression is derived from approximating a gas as point masses that have only kinetic energy and experience perfectly elastic collisions. Unfortunately, this theory breaks down at densities where intermolecular forces become significant, that is, when the potential energy is non-zero. Ideal gas law is suitable for only very dilute gasses.

Virial Theorem

To more accurately describe the properties of dilute gasses the Virial equation of state is used. Virial theorem accounts for the effects of intermolecular forces, through an expansion into higher order functions of the density. Mathematically, this is described using infinite power series, where $B_{2}(T)$ and $B_{3}(T)$ are the second and third Virial coefficients.

${\frac {P}{k_{B}T}}=\rho +B_{2}(T)\rho ^{2}+B_{3}(T)\rho ^{3}+...$

At low densities the deviations from ideal gas behavior can be sufficiently described in the second Virial coefficient, $B_{2}(T)$.

$B_{2}(T)=-{\frac {1}{2V}}(Z_{2}-Z_{1}^{2})$

Where $V$ is the volume of a gas and $Z$ is the configurational integral. The configurational integral for the second Virial coefficient is the contribution of every possible pair of positions weighted over its Boltzmann distribution. For the third Virial coefficient the configurational integral would be the contribution of three interacting particles. By extension, the configurational integral for any nth Virial coefficient would be the contribution of n particles interacting. For this reason higher Virial coefficients are quite complicated to derive, fortunately they are only necessary for describing gasses at pressures above 10atm. The configurational integral for two interacting particles is as follows:

$Z=\int \int e^{\frac {-u(r_{1},r_{2})}{k_{B}T}}dr_{1}dr_{2}$

Where $u$ is the potential energy of the interaction of a single pair of particles, and $r_{1}$ and $r_{2}$ are the positions of particles 1 and 2. The second virial coefficient can be written in terms of pairwise intermolecular interaction potential, $u_{r}$, if the position of particle 2 is defined relative to the position of particle 1. The distance, $r$, would then be the distance between the two interacting particles. The equation derived from this modification is as follows:

$B_{2}=-2\pi \int \limits _{0}^{\infty }\left[e^{\frac {-u(r)}{k_{B}T}}-1\right]r^{2}dr$

A different $B_{2}$ is derived for the hard sphere potential and the Lennard-Jones potential.

Hard Sphere Potential

The hard sphere model approximates particles as hard spheres that cannot overlap, if the spheres are not overlapping then the potential energy is zero and if they are overlapping then the potential energy is infinitely high. This approximation represents the very strong short range Pauli repulsion forces. The equation for the potential energy is as follows:

${\mathcal {V}}(r)={\begin{cases}{\mathcal {V}}=\infty &r<r_{cut}\\{\mathcal {V}}=0&r\geq r_{cut}\end{cases}}$

Where ${\mathcal {V}}$ is the potential energy and $r_{cut}$ is the radius of the hard sphere. Integrating the configurational integral for the hard sphere potential gives, $B_{2}(T)={\frac {2\pi }{3}}r_{cut}^{3}$, as the second Virial coefficient. This model is crude and only accounts for repulsive forces, a slightly more accurate model would be the Lennard-Jones potential model.

Lennard-Jones Potential

The Lennard-Jones potential is a combination of a polynomial repulsion term, $\upsilon (r)={\frac {C_{12}}{r^{12}}}$, and a London dispersion attractive term, $\upsilon (r)={\frac {C_{6}}{r^{6}}}$.

$\upsilon (r)={\frac {C_{12}}{r^{12}}}-{\frac {C_{6}}{r^{6}}}$

The $C_{12}$ and $C_{6}$ terms can be expanded, and internal energy, $u(r)$, can be expressed as:

$u(r)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]$

Where $\varepsilon$ is the potential well depth, $\sigma$ is the intercept, and $r$ is the distance between the particles. The second Virial coefficient derived from the Lennard-Jones potential has no analytical solution and must be solved numerically.

$B_{2}=\int \limits _{0}^{\infty }\left[e^{-{\frac {4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]}{k_{B}T}}}-1\right]r^{2}dr$

The Lennard-Jones model is a more accurate than the hard sphere model as it accounts for attractive interactions and the repulsive term is more realistic than hard sphere repulsion. That being said, it is still limited in the fact that only London dispersion attractive interactions are considered, making the model applicable to only nobel gasses.

Potential of mean force

The potential of mean force is one of the 3 major forces considered in stochastic dynamics models. The potential of mean force provides a free energy profile along a preferred coordinate, be it a geometric or energetic coordinate, such as the distance between two atoms or the torsional angle of a molecule. This free energy profile describes the average force of all possible configurations of a given system (the ensemble average of the force) on particles of interest.. The potential of mean force can be determined in both Monte Carlo Simulations as well as Molecular Dynamics.

Relation to Radial Distribution Function

For a liquid, the potential of mean force is related to the radial distribution function,

$w(r)=-k_{B}T\ln \left(g(r)\right)$

The potential of mean force can be calculated directly by using histogram analysis of a trajectory from an MD or MC simulation. This analysis calculates the probability of each possible configuration, and determines the free energy change associated with that state. In systems where parts of the coordinate of interest will not be sufficiently sampled from a naive simulation, umbrella sampling, free energy perturbation, or other enhanced sampling techniques can be used.

Derivation

$\langle F\rangle =-\langle {{\frac {d}{dr}}U(r^{N})}\rangle _{r_{1},r_{2}}={\frac {-\int dr_{3}\ldots dr_{N}{\frac {dU}{dr_{1}}}e^{-\beta U}}{\int dr_{3}\ldots dr_{N}e^{-\beta U}}}$

$=+k_{B}T{\frac {{\frac {d}{dr_{1}}}\int dr_{3}\ldots dr_{N}e^{-\beta U}}{\int dr_{3}\ldots dr_{N}e^{-\beta U}}}$

$=+k_{B}T{\frac {d}{dr_{1}}}\ln \int dr_{3}\ldots dr_{N}e^{-\beta U}$

$=+k_{B}T{\frac {d}{dr_{1}}}\ln(N(N-1){\frac {\int dr_{3}\ldots dr_{N}e^{-\beta U}}{\int dr^{N}e^{-\beta U}}}$

$=+k_{B}T{\frac {d}{dr_{1}}}\ln \left(g(r)\right)$

By integrating this function over 2 distance values, the reversible work associated with moving the two particles can be found. This work is equal to the change in Helmholtz energy associated with the change in the two states.

$W_{1\rightarrow 2}=\int _{r_{1}}^{r_{2}}\langle F_{1}\rangle dr_{1}=\Delta A_{1\rightarrow 2}$

Sources

1. Chandler, D. Introduction To Modern Statistical Mechanics, Oxford University Press, 1987, QC174.8.C47

2. Tuckerman, M. Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2010

See Also

Wikipedia:Potential of mean force

Langevin dynamics

Langevin dynamics is used to describe the acceleration of a particle in a liquid.

${\frac {dv}{dt}}=-\epsilon v(t)+A(t)$ 1

The term $-\epsilon v(t)$ term corresponds to the drag force divided by mass, m, of the particle. The drag force of the system is based on weak, long ranged intermolecular forces between atoms or molecules. This drag force is created by the liquid surrounding a particle. The drag force is based on the frictional constant of the system, $\gamma$, and the velocity of the particle, v. $\epsilon$ is ${\frac {\gamma }{m}}$. The frictional constant is proportional to the viscosity of the surrounding liquid, and the radius of the particle found from Stokes' Law.

The term $A(t)$ represents the random force. The random force of the system is based on short, strong Pauli repulsion between atoms or molecules. The random force in the liquid is experienced as a result of the collisions of a molecule with the surrounding liquid. This force changes rapidly with time, and changes much faster than it's velocity in denser liquids.

Velocity of the Particle

(1) is a first order linear differential equation and can be solved formally to give (2).

$v(t)=v_{0}e^{-\epsilon t}+e^{-\epsilon t}\int \limits _{1}^{3}e^{\epsilon t^{'}}A(t^{'})dt^{'}$ 2

The ensemble average of velocity is,

$\langle v(t)\rangle =v_{0}e^{-\epsilon t}$

This is because the acceleration of the particle is due to a random force, and so the average acceleration will be 0 (i.e., $\langle A(t)\rangle =0$).

At short time intervals ($t\longrightarrow 0$ ), the average velocity becomes,

$\langle v(t)\rangle =v_{0}e^{-\epsilon 0}=v_{0}$

At long time intervals ( $t\longrightarrow \infty$ ), the average velocity becomes,

$\langle v(t)\rangle =v_{0}e^{-\epsilon \infty }=0$

Displacement of the Particle

The Langevin Equation is used to express the rate of change of a particle's velocity. This in turn can be used to calculate the diffusion as diffusion depends on the velocity of a particle in a liquid. The Langevin equation can be used to sample the canonical ensembles of states. By integrating the velocity from time 0 to time t, the displacement of a particle can be found.

$r(t)-r(0)=\int \limits _{0}^{t}v(t^{''})dt^{''}$ 3

Substituting the definition for velocity from (2) into (3),and integrating by parts gives

$r(t)-r(0)={\frac {1}{\epsilon }}\int \limits _{0}^{t}[1-e^{\epsilon (t^{'}-t)}]A(t^{'})dt^{'}+v_{0}{\frac {1-e^{-\epsilon t}}{\epsilon }}$ 4

Squaring both sides of (4) and taking the ensemble average gives

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}(1-e^{-\epsilon t})^{2}}{\epsilon ^{2}}}+{\frac {3k_{B}T}{m\epsilon ^{2}}}(2\epsilon t-3+4e^{-\epsilon t}-e^{-2\epsilon t})$ 5

The derivation is as follows,

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}(1-e^{-\epsilon t})^{2}}{\epsilon ^{2}}}+{\frac {1}{\epsilon ^{2}}}\int \limits _{0}^{t}\int \limits _{0}^{t}[1-e^{-\epsilon (t-t^{'})}]^{2}\langle A(t^{''})A(t^{'})\rangle dt^{''}dt^{'}$

$\langle A(t^{''})A(t^{'})\rangle$ is a rapidly varying function of $\left|t^{'}-t^{''}\right|$ only and it is non zero only when $\left|t^{'}-t^{''}\right|$ is small. So $\langle A(t^{''})A(t^{'})\rangle$ can be re-written as $\phi \left|t^{'}-t^{''}\right|$ . Let $t^{'}-t^{''}=y$ and $t-t^{'}=x$, and the above equation becomes,

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}(1-e^{-\epsilon t})^{2}}{\epsilon ^{2}}}+{\frac {1}{\epsilon ^{2}}}{\frac {2}{\epsilon }}e^{-\epsilon t}-{\frac {1}{\epsilon ^{2}}}\int \limits _{0}^{2t}[1-e^{-\epsilon x}]^{2}dx\int \limits _{0}^{\infty }\phi |y|dy$

Let $\int \limits _{0}^{\infty }\phi |y|dy=\tau$, and simplify, the above equation becomes

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}(1-e^{-\epsilon t})^{2}}{\epsilon ^{2}}}+{\frac {\tau }{\epsilon ^{2}}}(2\epsilon t-3+4e^{-\epsilon t}-e^{-2\epsilon t})$

Equipartition theory applies when $t\longrightarrow \infty$, and so $\tau ={\frac {3k_{B}T}{m}}$, and the above equation becomes (5).

At short time intervals ($t\longrightarrow 0$) the part of the equation representing Brownian motion goes to zero ( i.e. ${\frac {3k_{B}T}{m\epsilon ^{2}}}(2\epsilon t-3+4e^{-\epsilon t}-e^{-2\epsilon t})=0$

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}(1-e^{-\epsilon t})^{2}}{\epsilon ^{2}}}$

Let $e^{-\epsilon t}=1-\epsilon t$. Then,

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {v_{0}^{2}}{\epsilon ^{2}}}(1-1+\epsilon t)^{2}$ 6

$\langle |r(t)-r(0)|^{2}\rangle =v_{0}^{2}t^{2}$

This corresponds to ballistic motion of the particle. At short time scales, the effect of friction and collisions have not affected the movement of the particle.

At long time intervals ($t\longrightarrow \infty$), the velocity of the particle goes to zero.

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {3k_{B}T}{m\epsilon ^{2}}}(2\epsilon t-3+4e^{-\epsilon t}-e^{-2\epsilon t})$

At large values of t $2\epsilon t>>4e^{-\epsilon t}{\text{and}}-e^{-2\epsilon t}$ and so,

$\langle |r(t)-r(0)|^{2}\rangle ={\frac {3k_{B}T}{m\epsilon ^{2}}}(2\epsilon t)={\frac {6k_{B}T}{m\epsilon }}t$

The diffusion constant D, is equal to ${\frac {k_{B}T}{m\epsilon }}$.

$\langle |r(t)-r(0)|^{2}\rangle =6Dt$ 7

This corresponds to the particle undergoing a random walk.

Periodic Boundary Conditions

A schematic representation of periodic boundary conditions. The unit cell of the system is in the center cell (pink). The potential energy of the system includes interactions between this center cell and its periodic images.

A schematic representation of the non-bonded interactions treated by the minimum image convention.

Liquids like water can be simulated using Periodic Boundary Conditions (PBC). A unit cell of the liquid is replicated periodically.

Molecular dynamics simulation of hexane in a periodic 30 Å × 30 Å × 30 Å simulation cell. Hydrogen atoms are omitted for clarity.

Macroscopic systems are extremely large and are, therefore, expensive to compute by molecular simulations. For instance, one gram of water has about 3 x 10²² molecules, a number too large to be calculated, even on computers. Fortunately, periodic boundary conditions enable us to mimic an infinite system by treating a relatively small part of a system to achieve a reasonable representation of the infinite system. The particles of this small subsystem are controlled by a set of boundary conditions called a unit cell (e.g., a three-dimensional box). During the simulation, particles are free to move in the central (original) cell; therefore, their periodic images of the adjacent cells move in an identical way. This means any particle that crosses one boundary of the cell, will reappear on the opposite side.

Non-Bonded Interactions Under Periodic Boundary Conditions

A particle in a periodic cell has non-bonded interactions with the other particles both in the original cell, and also in the surrounding cells (the periodic images of the original cell). The non-bonded potential energy (${\mathcal {V}}^{NB}$) can be written as:

${\mathcal {V}}^{NB}=\sum \limits _{s}^{}{\sum \limits _{i>j}^{}{v_{ij}^{NB}({r_{ij,s}})}}$

$v_{ij}^{NB}({r_{ij,s}})=4{\varepsilon _{ij}}\left[{{{\left({\frac {\sigma _{ij}}{r_{ij,s}}}\right)}^{12}}-{{\left({\frac {\sigma _{ij}}{r_{ij,s}}}\right)}^{6}}}\right]+{\frac {{q_{i}}{q_{j}}}{r_{ij,s}}}$

where ${r_{ij,s}}=\left|{{{\rm {r}}_{i}}-{{\rm {r}}_{j}}+s}\right|$, and ${s}$ is a translational vector to different cell images. For a periodic cubic cell of length ${L}$, ${\mathcal {s}}$ becomes

$s=\left[{\begin{array}{l}{n_{x}}L\\{n_{y}}L\\{n_{z}}L\end{array}}\right]$

where ${n_{x}}$, ${n_{y}}$, and ${n_{z}}$ are vectors of integers. The number of possible translational vectors, ${s}$, are infinite. This will lead us to an infinite number of non-bonded interactions. Thus, we need to perform some approximations to deal with this problem.

The first term of the $v_{ij}^{NB}({r_{ij,s}})$ is the Lennard-Jones potential. This potential has an attractive component (${C_{6}}{r^{-6}}$) for London dispersion, which is a short range interaction. Accordingly, if the distance between interacting particles increases (i.e., r > 10 Å), the interaction becomes very weak. Thus, the interactions can be truncated for distances greater than 10 A. However, the potential will lose its continuity because of the unexpected truncation. This problem can be fixed by employing a switching function, which will smoothly scale down van der Waals interactions to zero at the cutoff distance. For faster computation, pair lists are designed to search for particles that could be inside the cutoff area within a specified range ($r_{pairlistdist}$ is ~2 Å larger than $r_{cutoff}$), and, during the simulation, lists are updated.

The second term of the $v_{ij}^{NB}({r_{ij,s}})$ represents the long rang ($r^{-<!--\; -\; -->1}$