Molecular Forces & Interactions

A comprehensive treatment of the fundamental forces that govern biomolecular structure and function: van der Waals interactions, hydrogen bonds, electrostatic forces, and the hydrophobic effect. We derive key expressions from first principles and connect them to measurable biological phenomena.

Table of Contents

1. Van der Waals Forces & the Lennard-Jones Potential

Van der Waals forces are weak, short-range interactions between all atoms and molecules. They arise from fluctuations in electron distributions and comprise three contributions: orientation forces (Keesom, between permanent dipoles), induction forces (Debye, between permanent and induced dipoles), and dispersion forces (London, between instantaneous and induced dipoles).

The Lennard-Jones 12-6 Potential

The Lennard-Jones potential is the most widely used empirical pair potential in molecular simulations. It combines a short-range repulsive term (Pauli exclusion) with a long-range attractive term (London dispersion):

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

where:

  • $\varepsilon$ is the depth of the potential well (binding energy)
  • $\sigma$ is the distance at which $V(r) = 0$
  • • The $r^{-12}$ term represents short-range Pauli repulsion
  • • The $r^{-6}$ term represents London dispersion attraction

Hydrogen Bonds in Biological Systems

In proteins, backbone hydrogen bonds stabilize secondary structures:

  • Alpha-helix: N-H of residue $i$ bonds to C=O of residue $i-4$, energy ~12-17 kJ/mol per bond
  • Beta-sheet: Interstrand N-H...O=C bonds, energy ~12-17 kJ/mol per bond
  • DNA base pairing: A-T forms 2 H-bonds (~28 kJ/mol), G-C forms 3 H-bonds (~42 kJ/mol)

The net stabilization in aqueous solution is modest because forming an intramolecular H-bond requires breaking H-bonds with water. The free energy gain per H-bond in protein folding is only ~2-6 kJ/mol, a subtle balance between enthalpy and entropy.

4. Electrostatic Interactions & Debye-Hückel Screening

Coulomb's Law in a Dielectric Medium

The electrostatic interaction between two charges $q_1$ and $q_2$in a medium of relative permittivity $\varepsilon_r$ is:

$$V(r) = \frac{q_1 q_2}{4\pi\varepsilon_0 \varepsilon_r r}$$

For water at 25°C, $\varepsilon_r \approx 78.5$, which dramatically reduces electrostatic interactions. The Bjerrum length$l_B$ is the distance at which the Coulomb energy between two elementary charges equals the thermal energy:

$$l_B = \frac{e^2}{4\pi\varepsilon_0 \varepsilon_r k_B T}$$

In water at 25°C, $l_B \approx 0.71$ nm. In vacuum, $l_B \approx 55.6$ nm. This factor of ~80 reduction is why charged biomolecules can exist in aqueous solution.

Debye-Hückel Theory: Derivation

In an ionic solution, mobile ions screen electrostatic interactions. We derive the screened potential by combining the Poisson equation with the Boltzmann distribution.

Step 1: Poisson equation. The electrostatic potential $\phi(\mathbf{r})$ satisfies:

$$\nabla^2\phi = -\frac{\rho(\mathbf{r})}{\varepsilon_0\varepsilon_r}$$

where $\rho$ is the total charge density, including the central charge and the ionic atmosphere.

Step 2: Boltzmann distribution of ions. For a $z:z$ electrolyte at number density $n_0$, the local ion concentrations are:

$$n_\pm(\mathbf{r}) = n_0 \exp\!\left(\mp\frac{ze\phi(\mathbf{r})}{k_BT}\right)$$

Step 3: Linearization. For weak potentials ($ze\phi \ll k_BT$), we expand the exponential:

$$\rho_{\text{ion}} = zen_+ - zen_- \approx -\frac{2n_0 z^2 e^2}{k_BT}\phi$$

Step 4: Poisson-Boltzmann linearized. Combining:

$$\nabla^2\phi = \kappa^2\phi, \quad \text{where} \quad \kappa^2 = \frac{2n_0 z^2 e^2}{\varepsilon_0\varepsilon_r k_BT}$$

Step 5: Spherical solution. For a point charge $Q$ at the origin, the solution with boundary conditions$\phi \to 0$ as $r \to \infty$ is:

$$\phi(r) = \frac{Q}{4\pi\varepsilon_0\varepsilon_r r}\exp(-\kappa r)$$

The Debye screening length is:

$$\lambda_D = \kappa^{-1} = \sqrt{\frac{\varepsilon_0\varepsilon_r k_BT}{2n_0 z^2 e^2}}$$

For physiological salt concentration (~150 mM NaCl), $\lambda_D \approx 0.8$ nm. This means electrostatic interactions between charged residues on a protein surface are effectively screened beyond about 1 nm.

The Screened Coulomb (Yukawa) Potential

The interaction energy between two charges $q_1$ and $q_2$ in an ionic solution is the Yukawa potential:

$$V(r) = \frac{q_1 q_2}{4\pi\varepsilon_0\varepsilon_r r}\exp(-r/\lambda_D)$$

This exponential decay is why electrostatic interactions in biology are effectively short-ranged despite being formally long-ranged in vacuum. The energy at contact ($r = 0.5$ nm) between two elementary charges at physiological ionic strength is approximately $1\text{-}2\, k_BT$, comparable to thermal fluctuations.

5. The Hydrophobic Effect

The hydrophobic effect is the tendency of nonpolar molecules to aggregate in aqueous solution. It is the dominant driving force for protein folding, membrane formation, and micelle assembly. Counterintuitively, it is primarily an entropic phenomenon at room temperature.

Thermodynamic Derivation from Entropy

Consider transferring a nonpolar solute from a nonpolar environment into water. The Gibbs free energy of transfer is:

$$\Delta G_{\text{transfer}} = \Delta H_{\text{transfer}} - T\Delta S_{\text{transfer}}$$

The entropic argument: When a nonpolar solute enters water, water molecules form an ordered "cage" (clathrate-like structure) around the hydrophobic surface. This ordering reduces the entropy of water:

$$\Delta S_{\text{water}} < 0 \quad (\text{ordering of water cage})$$

The enthalpy change is near zero or slightly favorable at room temperature (van der Waals contacts between water and solute are replaced by water-water hydrogen bonds). Therefore:

$$\Delta G_{\text{transfer}} \approx -T\Delta S_{\text{water}} > 0 \quad (\text{unfavorable transfer to water})$$

The transfer free energy per unit area of nonpolar surface is approximately:

$$\frac{\Delta G}{\Delta A} \approx 25\text{-}50 \text{ mJ/m}^2 \approx 0.1\text{-}0.2 \; k_BT/\text{\AA}^2$$

Temperature Dependence and the Heat Capacity Signature

A hallmark of the hydrophobic effect is the large, negative heat capacity change:

$$\Delta C_p = \left(\frac{\partial \Delta H}{\partial T}\right)_p < 0$$

This means the enthalpy and entropy are both strongly temperature-dependent:

$$\Delta H(T) = \Delta H(T_0) + \Delta C_p(T - T_0)$$
$$\Delta S(T) = \Delta S(T_0) + \Delta C_p \ln\!\left(\frac{T}{T_0}\right)$$

There exists a temperature $T_S$ where $\Delta S = 0$:

$$T_S = T_0 \exp\!\left(-\frac{\Delta S(T_0)}{\Delta C_p}\right)$$

For many proteins, $T_S \approx 385$ K, above which $\Delta S > 0$ and the hydrophobic driving force reverses, contributing to thermal denaturation. Similarly, at low temperatures, entropy-driven stability weakens, contributing to cold denaturation of proteins.

Quantitative Model: Cavity Formation

The free energy of solvating a nonpolar solute can be decomposed:

$$\Delta G_{\text{solv}} = \Delta G_{\text{cavity}} + \Delta G_{\text{interaction}}$$

The cavity free energy is related to the probability $p_0$ of finding a solute-sized cavity in the solvent:

$$\Delta G_{\text{cavity}} = -k_BT \ln p_0$$

For small solutes (radius $< 1$ nm), $\Delta G_{\text{cavity}}$scales with volume. For large solutes, it crosses over to scaling with surface area, which is why the surface-area-based model works well for protein folding.

6. Boltzmann Distribution & Thermal Energy Scale

The Boltzmann Factor

The probability of a system being in a microstate with energy $E_i$at temperature $T$ is given by the Boltzmann distribution:

$$P(E_i) = \frac{e^{-E_i/k_BT}}{Z}, \quad Z = \sum_i e^{-E_i/k_BT}$$

where $Z$ is the partition function. The thermal energy scale $k_BT$ at room temperature (298 K) is:

$$k_BT \approx 4.11 \times 10^{-21}\text{ J} \approx 2.48\text{ kJ/mol} \approx 0.593\text{ kcal/mol} \approx 25.7\text{ meV}$$

This single number sets the scale for all molecular biophysics: processes with energy barriers $\gg k_BT$ are essentially frozen, while those with barriers $\sim k_BT$ occur spontaneously by thermal fluctuations.

Applications to Molecular Biophysics

Two-state equilibrium: For a molecule with two states separated by energy $\Delta E$:

$$\frac{P_2}{P_1} = e^{-\Delta E/k_BT}$$

Transition rates (Arrhenius): The rate of crossing an energy barrier of height $\Delta G^\ddagger$ is:

$$k = A\exp\!\left(-\frac{\Delta G^\ddagger}{k_BT}\right)$$

Force-energy relationship: A force$F$ applied over a distance $\delta$ does work $F\delta$. To be biologically relevant, this must be comparable to $k_BT$:

$$F \sim \frac{k_BT}{\delta} \sim \frac{4 \times 10^{-21}\text{ J}}{1\text{ nm}} \approx 4\text{ pN}$$

This piconewton force scale is precisely the range measured in single-molecule experiments on molecular motors, protein unfolding, and DNA stretching.

Summary: Energy and Force Scales in Biology

Hierarchy of Interaction Energies

The relative strengths of molecular interactions at room temperature (in units of $k_BT \approx 2.5$ kJ/mol):

  • Covalent bonds: ~150 $k_BT$ (essentially permanent)
  • Salt bridges (buried): ~5-15 $k_BT$
  • Hydrogen bonds: ~2-8 $k_BT$ (in vacuum), ~0-2 $k_BT$ (net in water)
  • Hydrophobic effect: ~1-2 $k_BT$ per methylene group
  • Van der Waals contacts: ~0.2-1 $k_BT$ per atom pair
  • Screened Coulomb (physiological salt): ~0.5-2 $k_BT$ at contact

Protein stability arises from the sum of thousands of weak interactions, with a typical net stability of only 5-15 $k_BT$ (12-40 kJ/mol). This marginal stability is essential for biological function: it allows proteins to undergo conformational changes, be regulated by ligand binding, and be degraded when necessary.

Related Video Lectures

Molecular Interactions are Due to Electrons

Induced Dipole-Dipole Interactions

Van der Waals: Buckingham & Lennard-Jones

Hydrogen Bonds Stabilize Biomolecules

7. Derivation: Lennard-Jones Potential & Second Virial Coefficient

We now present a rigorous derivation of the Lennard-Jones potential from the underlying physics of Pauli repulsion and London dispersion, and show how it connects to measurable thermodynamic quantities through the second virial coefficient.

Physical Origin of the Two Terms

Repulsive term ($r^{-12}$): At short range, the electron clouds of two atoms overlap. By the Pauli exclusion principle, electrons cannot occupy the same quantum state. The kinetic energy cost of orthogonalizing the wavefunctions scales exponentially with overlap. The exact repulsive energy from Hartree-Fock calculations is:

$$V_{\text{rep}}(r) = A' e^{-\beta r}$$

where $\beta$ is related to the orbital decay length. The Born-Mayer exponential form is more physically accurate, but the $r^{-12}$ power law is chosen as a computationally efficient approximation. We require that $r^{-12}$matches the exponential repulsion in the physically relevant range $0.9\sigma < r < 1.1\sigma$.

Attractive term ($r^{-6}$): The London dispersion interaction (derived in Section 8) gives an exact $-C_6/r^6$ dependence at long range. This is not an approximation but follows rigorously from second-order perturbation theory for the dipole-dipole interaction.

Combining: We write the total pair potential as:

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

Define $\sigma$ as the zero-crossing ($V(\sigma) = 0$) and $\varepsilon$ as the well depth ($V(r_{\min}) = -\varepsilon$). From $V(\sigma) = 0$:

$$\frac{C_{12}}{\sigma^{12}} = \frac{C_6}{\sigma^6} \quad\Rightarrow\quad C_{12} = C_6 \sigma^6$$

Substituting back and finding the minimum by setting $dV/dr = 0$:

$$\frac{dV}{dr} = -\frac{12 C_{12}}{r^{13}} + \frac{6 C_6}{r^{7}} = 0 \quad\Rightarrow\quad r_{\min}^6 = \frac{2C_{12}}{C_6} = 2\sigma^6$$
$$\boxed{r_{\min} = 2^{1/6}\sigma \approx 1.1225\,\sigma}$$

The well depth:

$$V(r_{\min}) = \frac{C_6\sigma^6}{(2\sigma^6)} - \frac{C_6}{(2\sigma^6)^{1}} = C_6\sigma^6 \cdot \frac{1}{4\sigma^{12}} - \frac{C_6}{2\sigma^6} \cdot \frac{1}{1}$$

More directly, evaluating at $r_{\min} = 2^{1/6}\sigma$:

$$V(r_{\min}) = C_6\sigma^6 \cdot \frac{1}{(2^{1/6}\sigma)^{12}} - \frac{C_6}{(2^{1/6}\sigma)^6} = \frac{C_6}{4\sigma^6} - \frac{C_6}{2\sigma^6} = -\frac{C_6}{4\sigma^6}$$

Setting $\varepsilon = C_6/(4\sigma^6)$, we obtain $C_6 = 4\varepsilon\sigma^6$ and$C_{12} = 4\varepsilon\sigma^{12}$, giving the standard form:

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

Force from the LJ Potential

The force between two atoms is $F(r) = -dV/dr$. Computing the derivative:

$$\frac{dV}{dr} = 4\varepsilon\left[-12\frac{\sigma^{12}}{r^{13}} + 6\frac{\sigma^6}{r^7}\right] = -\frac{24\varepsilon}{r}\left[2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

Therefore:

$$\boxed{F(r) = \frac{24\varepsilon}{r}\left[2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]}$$

The force changes sign at $r = r_{\min} = 2^{1/6}\sigma$: repulsive for $r < r_{\min}$and attractive for $r > r_{\min}$. The maximum attractive force occurs at:

$$\frac{dF}{dr} = 0 \quad\Rightarrow\quad r^* = \left(\frac{26}{7}\right)^{1/6}\sigma \approx 1.2445\,\sigma$$

with magnitude $F_{\max} \approx 2.397\,\varepsilon/\sigma$. This maximum attractive force sets the scale for adhesion between atoms.

Second Virial Coefficient from the LJ Potential

The second virial coefficient $B_2(T)$ connects the pair potential to the equation of state of a real gas. The virial expansion for pressure is:

$$\frac{P}{k_BT} = \rho + B_2(T)\rho^2 + B_3(T)\rho^3 + \cdots$$

where $\rho = N/V$ is the number density. The second virial coefficient is given by the Mayer cluster integral:

$$B_2(T) = -2\pi\int_0^{\infty}\left[e^{-V(r)/k_BT} - 1\right]r^2\,dr$$

For the LJ potential, this integral cannot be evaluated analytically in closed form. However, by substituting $u = r/\sigma$ and $T^* = k_BT/\varepsilon$(reduced temperature), we obtain:

$$B_2^*(T^*) = \frac{B_2(T)}{\frac{2}{3}\pi\sigma^3} = -3\int_0^{\infty}\left[\exp\!\left(-\frac{4}{T^*}\left(u^{-12} - u^{-6}\right)\right) - 1\right]u^2\,du$$

This integral is evaluated numerically. Key limiting behaviors:

  • High $T^*$: $B_2^* \to \frac{2}{3}\pi\sigma^3$ (hard-sphere limit, repulsion dominates)
  • Low $T^*$: $B_2^* \to -\infty$ (attraction dominates, condensation)
  • Boyle temperature $T_B^*$: $B_2(T_B) = 0$, numerically $T_B^* \approx 3.42$

The connection to van der Waals equation parameters is:

$$B_2(T) \approx b - \frac{a}{k_BT}$$

where $b = \frac{2}{3}\pi\sigma^3$ is the excluded volume and$a = \frac{2}{3}\pi\int_\sigma^\infty |V(r)|r^2\,dr$ measures the integrated attraction strength. This directly links the microscopic LJ parameters to the macroscopic van der Waals equation $(P + a/V^2)(V - b) = Nk_BT$.

8. Derivation: London Dispersion Forces from Perturbation Theory

We derive the London dispersion interaction rigorously from quantum mechanical second-order perturbation theory. This derivation reveals why every atom and molecule, regardless of its permanent dipole moment, experiences an attractive $r^{-6}$ interaction.

Step 1: The Dipole-Dipole Interaction Hamiltonian

Consider two neutral atoms A and B separated by a distance $R$ along the$z$-axis. Each atom consists of a nucleus and an electron cloud. The instantaneous displacement of electron $i$ in atom A from its nucleus is $\mathbf{r}_{Ai}$, creating an instantaneous dipole $\hat{\mathbf{p}}_A = -e\sum_i \mathbf{r}_{Ai}$.

The interaction between two dipoles separated by $R \gg |\mathbf{r}_{Ai}|, |\mathbf{r}_{Bj}|$(the multipole expansion) gives:

$$\hat{H}' = \frac{e^2}{4\pi\varepsilon_0 R^3}\sum_{i,j}\left(x_{Ai}x_{Bj} + y_{Ai}y_{Bj} - 2z_{Ai}z_{Bj}\right)$$

where $x, y, z$ are the Cartesian components of the electron displacements. Note: higher-order multipole terms ($R^{-4}$ dipole-quadrupole, $R^{-5}$quadrupole-quadrupole) are neglected.

Step 2: First-Order Correction Vanishes

The first-order energy correction is:

$$E^{(1)} = \langle 0_A 0_B | \hat{H}' | 0_A 0_B \rangle$$

For atoms in spherically symmetric ground states ($S$-states), the expectation value of any component of the position operator vanishes:

$$\langle 0_A | x_{Ai} | 0_A \rangle = \langle 0_A | y_{Ai} | 0_A \rangle = \langle 0_A | z_{Ai} | 0_A \rangle = 0$$

Since $\hat{H}'$ is bilinear in the coordinates of A and B, and the expectation values factor, $E^{(1)} = 0$. This is why atoms without permanent dipoles have no first-order electrostatic interaction.

Step 3: Second-Order Correction Gives the $R^{-6}$ Law

The second-order perturbation theory gives:

$$E^{(2)} = -\sum_{n \neq 0}\sum_{m \neq 0} \frac{|\langle n_A m_B | \hat{H}' | 0_A 0_B \rangle|^2}{(E_n^A - E_0^A) + (E_m^B - E_0^B)}$$

Since $\hat{H}' \propto R^{-3}$, the matrix element squared contributes $R^{-6}$. We can write $E^{(2)} = -C_6/R^6$ where:

$$C_6 = \frac{e^4}{(4\pi\varepsilon_0)^2}\sum_{n,m \neq 0} \frac{|\langle n_A | \mathbf{r}_A | 0_A \rangle|^2 \cdot |\langle m_B | \mathbf{r}_B | 0_B \rangle|^2 \cdot f(\text{angular})}{(E_n^A - E_0^A) + (E_m^B - E_0^B)}$$

The angular factor, after averaging over orientations of the two transition dipoles, contributes a factor of 6 (from $x_Ax_B + y_Ay_B - 2z_Az_B$ giving terms like$1 + 1 + 4 = 6$).

Step 4: London's Approximation

The exact sum over excited states is intractable. London's brilliant approximation replaces the excitation energies by the ionization energies:

$$E_n^A - E_0^A \approx I_A, \quad E_m^B - E_0^B \approx I_B$$

This allows us to factor the sum. Using the oscillator strength sum rule and the static polarizability:

$$\alpha_A = \frac{2e^2}{4\pi\varepsilon_0}\sum_{n \neq 0}\frac{|\langle n_A | z_A | 0_A \rangle|^2}{E_n^A - E_0^A} \approx \frac{2e^2}{4\pi\varepsilon_0}\frac{1}{I_A}\sum_{n \neq 0}|\langle n_A | z_A | 0_A \rangle|^2$$

we can express the sum over matrix elements in terms of polarizability. After careful evaluation, the $C_6$ coefficient becomes:

$$\boxed{C_6 = \frac{3}{4}\frac{\alpha_A \alpha_B}{(4\pi\varepsilon_0)^2} \cdot \frac{I_A I_B}{I_A + I_B}}$$

The factor of $3/4$ arises from combining the angular factor of 6 with the denominator structure $1/(I_A + I_B) = 1/I_A \cdot 1/I_B \cdot I_AI_B/(I_A + I_B)$and normalization. The full London dispersion energy is:

$$\boxed{V_{\text{London}}(R) = -\frac{C_6}{R^6} = -\frac{3}{4}\frac{\alpha_A \alpha_B}{(4\pi\varepsilon_0)^2} \cdot \frac{I_A I_B}{I_A + I_B} \cdot \frac{1}{R^6}}$$

Universality of Dispersion Forces

A remarkable feature of the London formula is its universality:

  • Every atom and molecule has dispersion forces. The interaction requires only a nonzero polarizability $\alpha > 0$, which is true for all matter (even noble gases). No permanent dipole is needed.
  • Always attractive. Since $\alpha > 0$ and$I > 0$ for all species, $C_6 > 0$ and $V_{\text{London}} < 0$ always.
  • Scales with size. Larger atoms/molecules have greater polarizability, giving stronger dispersion. For example,$C_6(\text{Xe-Xe}) \approx 10 \times C_6(\text{He-He})$.
  • Dominates for nonpolar species. For nonpolar molecules like alkanes, dispersion is the sole source of intermolecular attraction. This explains why boiling points increase with chain length.

Representative $C_6$ values (in units of $10^{-79}$ J·m$^6$):

  • • He-He: $C_6 \approx 1.5$
  • • Ne-Ne: $C_6 \approx 6.3$
  • • Ar-Ar: $C_6 \approx 64$
  • • Kr-Kr: $C_6 \approx 130$
  • • Xe-Xe: $C_6 \approx 286$

9. Derivation: Debye-Hückel Screening

We derive the screened Coulomb potential from the linearized Poisson-Boltzmann equation in full detail, establishing the Debye length and its physical consequences for electrostatic interactions in ionic solutions.

Step 1: The Poisson Equation

We begin with Gauss's law in a linear dielectric medium. The electrostatic potential$\phi(\mathbf{r})$ generated by a total charge density $\rho_{\text{total}}$satisfies Poisson's equation:

$$\nabla^2\phi(\mathbf{r}) = -\frac{\rho_{\text{total}}(\mathbf{r})}{\varepsilon_0 \varepsilon_r}$$

The total charge density consists of the fixed charge (the ion of interest) plus the mobile ionic atmosphere: $\rho_{\text{total}} = \rho_{\text{fixed}} + \rho_{\text{ions}}$.

Step 2: Boltzmann Distribution of Ions

In thermal equilibrium at temperature $T$, the local concentration of ions of species $i$ with valence $z_i$ follows the Boltzmann distribution:

$$n_i(\mathbf{r}) = n_i^0 \exp\!\left(-\frac{z_i e \phi(\mathbf{r})}{k_BT}\right)$$

where $n_i^0$ is the bulk number concentration of species $i$ far from the central charge. For a symmetric $z{:}z$ electrolyte (e.g., NaCl where $z = 1$), we have cations ($+z$) and anions ($-z$) both at bulk concentration $n_0$:

$$n_+(\mathbf{r}) = n_0 e^{-ze\phi/k_BT}, \quad n_-(\mathbf{r}) = n_0 e^{+ze\phi/k_BT}$$

The net ionic charge density is:

$$\rho_{\text{ions}} = ze\,n_+ - ze\,n_- = zen_0\left(e^{-ze\phi/k_BT} - e^{+ze\phi/k_BT}\right) = -2zen_0\sinh\!\left(\frac{ze\phi}{k_BT}\right)$$

Step 3: The Full Poisson-Boltzmann Equation

Substituting the ionic charge density into Poisson's equation (away from the central charge):

$$\nabla^2\phi = \frac{2zen_0}{\varepsilon_0\varepsilon_r}\sinh\!\left(\frac{ze\phi}{k_BT}\right)$$

This is the nonlinear Poisson-Boltzmann (PB) equation, which is generally difficult to solve analytically.

Step 4: Debye-Hückel Linearization

The key approximation: when the electrostatic energy is small compared to thermal energy ($ze\phi \ll k_BT$, i.e., $\phi \ll 25.7$ mV at room temperature for monovalent ions), we can linearize $\sinh(x) \approx x$:

$$\nabla^2\phi = \frac{2n_0 z^2 e^2}{\varepsilon_0\varepsilon_r k_BT}\phi \equiv \kappa^2\phi$$

where we define the inverse Debye length:

$$\boxed{\kappa = \frac{1}{\lambda_D} = \sqrt{\frac{2n_0 z^2 e^2}{\varepsilon_0\varepsilon_r k_BT}}}$$

Equivalently, the Debye screening length is:

$$\boxed{\lambda_D = \kappa^{-1} = \sqrt{\frac{\varepsilon_0\varepsilon_r k_BT}{2n_0 z^2 e^2}}}$$

Step 5: Solving the Linearized PB Equation in Spherical Coordinates

For a point charge $q$ at the origin, we seek a spherically symmetric solution. In spherical coordinates, the Laplacian is:

$$\nabla^2\phi = \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d\phi}{dr}\right) = \frac{1}{r}\frac{d^2}{dr^2}(r\phi) = \kappa^2\phi$$

Letting $u(r) = r\phi(r)$, the equation becomes:

$$\frac{d^2u}{dr^2} = \kappa^2 u$$

The general solution is $u(r) = Ae^{-\kappa r} + Be^{+\kappa r}$. Applying boundary conditions:

  • $\phi \to 0$ as $r \to \infty$: requires $B = 0$
  • $\phi \to q/(4\pi\varepsilon_0\varepsilon_r r)$ as $r \to 0$ (bare Coulomb near the charge): gives $A = q/(4\pi\varepsilon_0\varepsilon_r)$

The screened Coulomb (Yukawa) potential is therefore:

$$\boxed{\phi(r) = \frac{q}{4\pi\varepsilon_0\varepsilon_r r}\exp(-\kappa r)}$$

Physical Interpretation: Effective Short-Range Interactions

The exponential factor $e^{-\kappa r}$ screens the bare Coulomb potential. The physical picture: mobile ions of opposite sign cluster near the central charge, forming an ionic atmosphere that neutralizes its field beyond $\lambda_D$. Key consequences:

  • At $r = \lambda_D$: the potential is reduced to $1/e \approx 37\%$ of the bare Coulomb value
  • At $r = 3\lambda_D$: the potential is reduced to $e^{-3} \approx 5\%$
  • At $r = 5\lambda_D$: the potential is reduced to $e^{-5} \approx 0.7\%$

Debye lengths at 25°C in water for various NaCl concentrations:

  • • 1 mM: $\lambda_D \approx 9.6$ nm
  • • 10 mM: $\lambda_D \approx 3.0$ nm
  • • 100 mM: $\lambda_D \approx 0.96$ nm
  • • 150 mM (physiological): $\lambda_D \approx 0.78$ nm

At physiological ionic strength, electrostatic interactions are screened within about 1 nm. This means that two charged amino acids on a protein surface effectively interact only when they are nearly in contact. Electrostatic steering of protein-protein encounters occurs at longer range only because large patches of charge create collective electrostatic fields that extend beyond individual Debye lengths.

10. Derivation: Hydrophobic Effect & Scaled Particle Theory

We derive the free energy of cavity formation in water using scaled particle theory (SPT), revealing the fundamental crossover between entropy-dominated (small solutes) and enthalpy-dominated (large solutes) regimes of hydrophobic solvation.

Scaled Particle Theory: Setup

The free energy of inserting a hard sphere of radius $R$ into a solvent of molecular radius $r_s$ and number density $\rho_s$ is related to the probability $p_0(R)$ of finding a cavity of radius $R + r_s$ (the closest approach distance) that is devoid of solvent centers:

$$\Delta G_{\text{cavity}} = -k_BT \ln p_0(R)$$

For water, $r_s \approx 0.14$ nm and $\rho_s \approx 33.3$ nm$^{-3}$(55.5 M).

Small Solute Regime ($R \lesssim 1$ nm): Volume Scaling

For small cavities (comparable to the size of a water molecule), we can use the exact result from the theory of hard-sphere fluids. The probability of finding a cavity of radius $R_c = R + r_s$ free of solvent centers is, in the limit of small $R_c$:

$$p_0(R_c) = \exp\left(-\rho_s \frac{4}{3}\pi R_c^3\right) \quad \text{(Poisson, for small } R_c\text{)}$$

This gives:

$$\Delta G_{\text{cavity}}^{\text{small}} = k_BT \cdot \rho_s \cdot \frac{4}{3}\pi R_c^3 \propto V$$

The free energy scales with the excluded volume. This is entropy-dominated: inserting a small solute reduces the number of configurations available to solvent molecules. The entropy cost is:

$$\Delta S_{\text{cavity}}^{\text{small}} = -k_B \rho_s \frac{4}{3}\pi R_c^3 < 0$$

with $\Delta H \approx 0$ for small cavities (no hydrogen bonds need to be broken, the water network simply rearranges around the solute). This is the molecular basis of the iceberg model: water forms structured cages around small nonpolar solutes without breaking H-bonds.

Large Solute Regime ($R \gtrsim 1$ nm): Area Scaling

For large cavities, the interface between the cavity and water behaves like a macroscopic liquid-vapor interface. The free energy becomes:

$$\Delta G_{\text{cavity}}^{\text{large}} = \gamma \cdot 4\pi R_c^2 + \kappa_{\text{curv}} \cdot \frac{4}{3}\pi R_c^3 \propto A$$

where $\gamma \approx 72$ mJ/m$^2$ is the liquid-vapor surface tension of water at 25°C, and $\kappa_{\text{curv}}$ is a curvature correction (the Tolman length correction). The dominant $\gamma A$ term is:

$$\boxed{\Delta G_{\text{cavity}}^{\text{large}} \approx \gamma \cdot A}$$

This regime is enthalpy-dominated: creating a large cavity requires breaking hydrogen bonds at the interface. Each water molecule at the surface loses approximately one H-bond (~10 kJ/mol), and the surface density of such molecules gives the macroscopic surface tension.

Crossover Length and the Full SPT Expression

The crossover between volume and area scaling occurs at a characteristic length$R^*$ where the two contributions are comparable. Setting$k_BT\rho_s V \sim \gamma A$:

$$k_BT\rho_s \cdot \frac{4}{3}\pi R^{*3} \sim \gamma \cdot 4\pi R^{*2}$$
$$\boxed{R^* \sim \frac{3\gamma}{k_BT\rho_s} \approx \frac{3 \times 0.072}{4.1 \times 10^{-21} \times 3.33 \times 10^{28}} \approx 1.6 \text{ nm}}$$

The full scaled particle theory expression interpolates smoothly between these regimes. For a hard-sphere solute of radius $R$ in a hard-sphere solvent:

$$\frac{\Delta G_{\text{cavity}}}{k_BT} = -\ln(1-\eta) + \frac{3\eta}{1-\eta}\frac{R}{r_s} + \frac{3\eta}{1-\eta}\left(\frac{R}{r_s}\right)^2 + \frac{P}{\rho_s k_BT}\frac{4}{3}\pi(R+r_s)^3$$

where $\eta = (4/3)\pi r_s^3 \rho_s$ is the packing fraction ($\eta \approx 0.38$for water) and $P$ is the pressure. The $R^2$ term is the dominant surface-area contribution for large $R$, while the $R^3$ pressure-volume work term dominates at very large scales (relevant for macroscopic phase separation).

Implications for Biophysics

The crossover explains several biological phenomena:

  • Small hydrophobes (methane, noble gases): $R \sim 0.2$ nm, entropy-dominated.$\Delta G \propto V$. Water can maintain its H-bond network by forming clathrate cages. Large negative $\Delta S$, small $\Delta H$.
  • Medium hydrophobes (amino acid side chains): $R \sim 0.3\text{-}0.5$ nm. Intermediate regime. The cavity free energy per unit area is ~25 cal/(mol·A$^2$), the empirical value used in protein folding models.
  • Large hydrophobes (protein cores, membranes): $R > 1$ nm, enthalpy-dominated.$\Delta G \propto A$. Water cannot maintain H-bonds at the interface, leading to a drying transition. This drives protein folding, membrane assembly, and aggregation of hydrophobic particles.

11. Derivation: DLVO Theory

DLVO theory (Derjaguin-Landau-Verwey-Overbeek) describes the interaction between charged colloidal particles by combining van der Waals attraction with electrostatic double-layer repulsion. This framework is foundational for understanding colloidal stability, protein crystallization, and nanoparticle self-assembly.

Van der Waals Attraction Between Two Flat Surfaces

To obtain the interaction between macroscopic bodies, we integrate the $-C_6/r^6$pair interaction over all pairs of atoms. For two infinite flat surfaces separated by distance $d$, with atomic number densities $\rho_1$ and $\rho_2$:

$$\frac{V_{\text{vdW}}(d)}{A_{\text{surface}}} = -\rho_1\rho_2\int_0^\infty dz_1\int_0^\infty dz_2\int_0^\infty 2\pi r\,dr\;\frac{C_6}{(r^2 + (d+z_1+z_2)^2)^3}$$

Evaluating the integrals step by step. The radial integral gives:

$$\int_0^\infty \frac{2\pi r\,dr}{(r^2 + D^2)^3} = \frac{\pi}{2D^4}$$

where $D = d + z_1 + z_2$. The $z_1$ and $z_2$ integrals then give:

$$\int_0^\infty\int_0^\infty \frac{dz_1\,dz_2}{(d+z_1+z_2)^4} = \frac{1}{6d^2}$$

Combining and defining the Hamaker constant$A_H = \pi^2 C_6 \rho_1 \rho_2$:

$$\boxed{\frac{V_{\text{vdW}}(d)}{A_{\text{surface}}} = -\frac{A_H}{12\pi d^2}}$$

This is the non-retarded van der Waals interaction per unit area between two half-spaces. The $d^{-2}$ dependence (much slower decay than the$r^{-6}$ pair interaction) means that van der Waals forces between macroscopic bodies are significant at much larger separations than between individual atoms.

The Hamaker Constant from Lifshitz Theory

The simple pairwise summation (Hamaker's original approach) neglects many-body effects and retardation. Lifshitz theory provides a more rigorous treatment by treating the van der Waals interaction as arising from fluctuating electromagnetic fields. The non-retarded Hamaker constant for two bodies 1 and 2 interacting across medium 3 is:

$$A_H = \frac{3}{4}k_BT\sum_{n=0}^{\infty}{}' \sum_{s=1}^{\infty}\frac{1}{s^3}\left(\frac{\varepsilon_1(i\xi_n) - \varepsilon_3(i\xi_n)}{\varepsilon_1(i\xi_n) + \varepsilon_3(i\xi_n)}\right)^s\left(\frac{\varepsilon_2(i\xi_n) - \varepsilon_3(i\xi_n)}{\varepsilon_2(i\xi_n) + \varepsilon_3(i\xi_n)}\right)^s$$

where $\xi_n = 2\pi n k_BT/\hbar$ are the Matsubara frequencies and the prime on the sum means the $n = 0$ term is weighted by 1/2. For the leading term ($s = 1$):

$$A_H \approx \frac{3}{4}k_BT\left(\frac{\varepsilon_1 - \varepsilon_3}{\varepsilon_1 + \varepsilon_3}\right)\left(\frac{\varepsilon_2 - \varepsilon_3}{\varepsilon_2 + \varepsilon_3}\right) + \frac{3\hbar}{4\pi}\int_{\xi_1}^{\infty}\frac{\Delta_{13}(i\xi)\Delta_{23}(i\xi)}{1}\,d\xi$$

In the simplified Tabor-Winterton approximation for identical particles ($\varepsilon_1 = \varepsilon_2$) in vacuum ($\varepsilon_3 = 1$):

$$A_H \approx \frac{3}{4}k_BT\left(\frac{\varepsilon - 1}{\varepsilon + 1}\right)^2 + \frac{3h\nu_e}{16\sqrt{2}}\frac{(n^2 - 1)^2}{(n^2 + 1)^{3/2}}$$

where $n$ is the refractive index and $\nu_e$ is the main electronic absorption frequency. Typical Hamaker constants: proteins in water $A_H \sim 1\text{-}5 \times 10^{-21}$ J ($\sim 0.25\text{-}1.2\,k_BT$); metals in vacuum $A_H \sim 4 \times 10^{-19}$ J ($\sim 100\,k_BT$).

Electrostatic Double-Layer Repulsion

Two like-charged surfaces in an electrolyte repel each other because their diffuse double layers overlap. For two flat surfaces with surface potential $\psi_0$(or equivalently, surface charge density $\sigma_s$), the Debye-Hückel theory gives the interaction energy per unit area:

$$\frac{V_{\text{DH}}(d)}{A_{\text{surface}}} = \frac{64 n_0 k_BT}{\kappa}\tanh^2\!\left(\frac{ze\psi_0}{4k_BT}\right)\exp(-\kappa d)$$

In the Debye-Hückel limit ($ze\psi_0 \ll 4k_BT$), $\tanh(x) \approx x$and:

$$\boxed{\frac{V_{\text{DH}}(d)}{A_{\text{surface}}} = \frac{64 n_0 k_BT}{\kappa}\left(\frac{ze\psi_0}{4k_BT}\right)^2\exp(-\kappa d) = \frac{n_0 z^2 e^2 \psi_0^2}{\kappa k_BT}\exp(-\kappa d)}$$

The repulsion decays exponentially with decay length $\kappa^{-1} = \lambda_D$(the Debye length). Adding salt increases $\kappa$, compresses the double layer, and reduces the repulsion.

The Total DLVO Interaction

The DLVO theory sums the van der Waals attraction and the double-layer repulsion:

$$\boxed{V_{\text{DLVO}}(d) = V_{\text{vdW}}(d) + V_{\text{DH}}(d) = -\frac{A_H}{12\pi d^2} + \frac{64 n_0 k_BT}{\kappa}\Gamma^2 e^{-\kappa d}}$$

where $\Gamma = \tanh(ze\psi_0/4k_BT)$ (per unit area). The competition between the power-law attraction and exponential repulsion creates a characteristic profile with up to three key features:

  • Primary minimum: At very short separations ($d \sim 0.2$ nm), the van der Waals attraction dominates. This is a deep, often irreversible minimum (coagulation/fusion).
  • Energy barrier: At intermediate separations ($d \sim 1\text{-}5$ nm), the exponential repulsion exceeds the power-law attraction, creating a barrier $V_{\max}$. If $V_{\max} \gg k_BT$, the system is kinetically stable (the particles cannot overcome the barrier by thermal fluctuations).
  • Secondary minimum: At larger separations ($d \sim 5\text{-}20$ nm), there may be a shallow minimum of depth ~$1\text{-}5\,k_BT$. Particles can reversibly aggregate in this minimum (flocculation).
DLVO theory diagram showing the total interaction energy as a function of particle separation, with contributions from van der Waals attraction and double-layer repulsion
DLVO interaction energy: the total potential (bold curve) results from summing van der Waals attraction and electrostatic double-layer repulsion, producing a characteristic energy barrier and secondary minimum. — Source: Wikimedia Commons (CC-BY-SA)

Stability Criterion

The critical condition for colloidal stability is that the energy barrier$V_{\max} \gg k_BT$. The barrier height depends on the ratio of the Hamaker constant to the double-layer parameters. The critical coagulation concentration (CCC) is found by simultaneously requiring $V_{\text{DLVO}} = 0$ and$dV_{\text{DLVO}}/dd = 0$:

$$-\frac{A_H}{12\pi d_c^2} + \frac{64 n_0 k_BT}{\kappa}\Gamma^2 e^{-\kappa d_c} = 0$$
$$\frac{A_H}{6\pi d_c^3} + \frac{64 n_0 k_BT \kappa}{\kappa}\Gamma^2 e^{-\kappa d_c} \cdot (-\kappa) = 0$$

Dividing the second equation by the first gives $\kappa d_c = 2$, i.e., the barrier occurs at $d_c = 2/\kappa = 2\lambda_D$. Substituting back and solving for$n_0$, the Schulze-Hardy rule emerges:

$$\boxed{n_{\text{CCC}} \propto \frac{1}{z^6}} \quad \text{(for high surface potentials)}$$

This dramatic $z^{-6}$ dependence on valence explains why divalent ions (Ca$^{2+}$) are ~64 times more effective at destabilizing colloids than monovalent ions (Na$^+$), and trivalent ions (Al$^{3+}$) are ~729 times more effective. This is why aluminum salts are used as flocculants in water treatment.

12. Applications in Molecular Biophysics

Protein-Protein Interactions

The second virial coefficient $B_{22}$ measured by static light scattering captures the net protein-protein interaction. DLVO-like analysis reveals that protein interactions are controlled by a delicate balance: van der Waals attraction ($A_H \sim 3\text{-}10\,k_BT$for proteins), electrostatic repulsion (tunable via pH and ionic strength), and hydrophobic patches. Crystallization occurs when $B_{22}$ falls within a narrow "crystallization slot" ($-8 \times 10^{-4}$ to $-2 \times 10^{-4}$ mL·mol/g$^2$).

DNA Condensation

DNA is a highly charged polyelectrolyte ($-2e$ per base pair, linear charge density$\sim 6\,e/\text{nm}$). In the presence of multivalent cations (spermidine$^{3+}$, spermine$^{4+}$, Co(NH$_3$)$_6^{3+}$), DNA condenses into compact toroids. Manning counterion condensation theory predicts condensation when the fraction of charge neutralized by counterions exceeds ~89% for $z = 3$. The balance between electrostatic attraction (from correlated counterion fluctuations) and bending rigidity ($l_p \approx 50$ nm) determines the toroid radius.

Colloidal Stability

Colloidal particles (1 nm to 1 $\mu$m) are stabilized against aggregation by electrostatic repulsion (charge stabilization) or steric repulsion (polymer coating). DLVO theory quantitatively predicts the critical salt concentration for destabilization. Biological examples: casein micelles in milk (stabilized by $\kappa$-casein), blood cells (stabilized by the glycocalyx), and liposomes for drug delivery.

Self-Assembly

Molecular self-assembly is driven by the interplay of all molecular forces. Amphiphilic molecules (lipids, surfactants, block copolymers) self-assemble when the critical micelle concentration (CMC) is reached. The CMC depends on the hydrophobic free energy per molecule:$\ln(\text{CMC}) \approx -\Delta G_{\text{hydrophobic}}/(k_BT)$. The assembled structure (micelles, vesicles, bilayers) is determined by the packing parameter $p = v/(a_0 l_c)$, where $v$ is the hydrophobic volume, $a_0$ is the optimal head-group area, and $l_c$ is the chain length.

Drug-Receptor Binding & Molecular Recognition

The binding affinity between a drug and its receptor is determined by the net free energy change $\Delta G_{\text{bind}} = -k_BT\ln K_d$, where $K_d$ is the dissociation constant. This free energy is a sum of contributions from all molecular forces: van der Waals contacts at the binding interface (~0.1-0.5 $k_BT$ per atom pair), hydrogen bonds (~1-3 $k_BT$ net per H-bond in water), salt bridges (~2-5 $k_BT$ if desolvation is favorable), and the hydrophobic effect (~0.04 $k_BT$ per A$^2$ of buried nonpolar surface). The entropic penalty of immobilizing the ligand ($-T\Delta S_{\text{trans+rot}} \sim +5\text{-}10\,k_BT$) must be overcome by these favorable interactions. Typical drug affinities range from$K_d \sim$ nM ($\Delta G \sim -20\,k_BT$) to $\mu$M ($\Delta G \sim -14\,k_BT$).

13. Historical Context

Timeline of Key Developments

  • 1873 — Johannes Diderik van der Waals proposed his equation of state$(P + a/V^2)(V - b) = RT$, introducing the concepts of intermolecular attraction ($a$) and finite molecular volume ($b$). This was the first quantitative recognition that molecules interact beyond hard-core repulsion. Van der Waals received the Nobel Prize in Physics (1910) for this work.
  • 1912 — Peter Debye derived the interaction between a permanent dipole and an induced dipole (induction or Debye forces, $\propto r^{-6}$).
  • 1921 — Willem Hendrik Keesom derived the thermally averaged interaction between two permanent dipoles (orientation or Keesom forces, $\propto r^{-6}$).
  • 1923 — Peter Debye and Erich Hückel published their theory of electrolyte solutions, deriving the screened Coulomb potential and the concept of the ionic atmosphere. This explained the concentration dependence of activity coefficients and conductivities in electrolyte solutions. The linearization of the Poisson-Boltzmann equation remains one of the most useful approximations in physical chemistry.
  • 1930 — Fritz London applied quantum mechanical second-order perturbation theory to derive the dispersion interaction between two neutral atoms. London showed that the $-C_6/r^6$ attraction is universal, explaining for the first time why noble gases condense. This resolved a long-standing puzzle: van der Waals had postulated attractive forces in 1873, but their origin was unknown until London's quantum mechanical derivation.
  • 1937 — Hugo Christiaan Hamaker developed the pairwise summation method for calculating van der Waals forces between macroscopic bodies, introducing the Hamaker constant.
  • 1941-1948 — DLVO Theory was developed independently by Boris Derjaguin and Lev Landau (1941, Soviet Union) and by Evert Verwey and Theo Overbeek (1948, Netherlands). They combined van der Waals attraction with electrostatic double-layer repulsion to explain the stability of lyophobic colloids. The Schulze-Hardy rule ($\text{CCC} \propto z^{-6}$) was explained quantitatively for the first time.
  • 1956 — Evgeny Lifshitz reformulated the theory of van der Waals forces between macroscopic bodies using quantum electrodynamics, eliminating the need for pairwise additivity and naturally incorporating retardation effects. The Lifshitz theory expresses the Hamaker constant in terms of measurable dielectric spectra.
  • 1969 — Walter Kauzmann published his seminal review on the hydrophobic effect in protein folding, establishing the entropic nature of hydrophobic interactions at room temperature. The concept of the "hydrophobic bond" (now called the hydrophobic effect) became central to protein biophysics.
  • 1976-present — Jacob Israelachvili and colleagues developed the Surface Forces Apparatus (SFA) for directly measuring forces between surfaces with sub-angstrom distance resolution. SFA measurements confirmed DLVO theory, revealed non-DLVO forces (hydration forces, hydrophobic forces, oscillatory solvation forces), and measured the range and magnitude of the hydrophobic attraction. Israelachvili's textbook Intermolecular and Surface Forces (1985, 3rd ed. 2011) remains the definitive reference in the field.

14. Interactive Simulations

Lennard-Jones Potential & Debye-Huckel Screening

Python

Compute and visualize the LJ potential, its components, and Debye-Huckel screening lengths at various ionic strengths.

script.py87 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Boltzmann Distribution & Hydrophobic Effect

Python

Analyze the thermodynamics of the hydrophobic effect, hydrogen bond energetics, and van der Waals interactions between macroscopic bodies.

script.py101 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server