Lesson 5.4 · 16 min

Polarizable models

Fixed-charge force fields give every atom one charge forever. Real electron clouds shift in an electric field. Polarizable models add this response, as induced dipoles, Drude oscillators or fluctuating charges, at two to four times the cost.

Why fixed charges miss induction, and how big the induction energy is near an ion
Why polarization is many-body (not pairwise additive), and the polarization catastrophe
The three common models (induced dipoles, Drude oscillators, fluctuating charges), their cost and when to use them

The idea

In a fixed-charge force field (Lessons 4.3 and 4.5) each atom carries one partial charge, chosen once. But electrons respond to their surroundings: an electric field E pulls the electron cloud one way and the nucleus the other, creating an induced dipole

μ_ind = α E

where α is the polarizability (about 1.45 Å³ for a water molecule, 0.03 Å³ for Li⁺, 4.5 Å³ for Cl⁻). The induced dipole lowers the energy by ½αE². Next to a Li⁺ ion at 2 Å this is about −15 kcal/mol, as large as a hydrogen bond, but it falls off as r⁻⁴ and is under 0.2 kcal/mol at 6 Å (code tab, Part 1).

Polarization is many-body. The induced dipole depends on the total field, and the energy goes as the square of it. Three ions on the same side of a molecule give −41 kcal/mol; the same three ions spread evenly around it give zero, because their fields cancel. A sum of pair terms gives −31 in both cases (Part 2). Induced dipoles also polarize each other (mutual polarization), which is why water's dipole grows from 1.85 D in the gas to about 2.5–3 D in the liquid.

The polarization catastrophe. Two point dipoles that get too close keep reinforcing each other, and below r = (4αα′)^1/6 (1.43 Å for two water-like sites in line) the induced dipoles become infinite (Part 3). Real electron clouds overlap and the field is weaker, so models damp the short-range dipole field, most often with Thole's smeared-charge functions.

Three ways to put polarization into MD:

Induced point dipoles (AMOEBA, AMBER ff02/ipolq): solve μ = α(E₀ + Tμ) every step, iteratively to self-consistency (or propagate the dipoles with an extended Lagrangian).
Drude oscillators (CHARMM Drude, SWM4-NDP water, OPLS/CL&Pol for ionic liquids): attach a charged massive particle to each atom with a spring; α = q_D²/k_D. It only needs ordinary charges and springs, so any MD code can run it. The Drude particles are kept "cold" (about 1 K) with a separate thermostat, which is close to keeping them relaxed.
Fluctuating charges (CHARMM FQ): charges flow within each molecule to equalize electronegativity. Cheap, but no out-of-plane response.

Cost and when to use it. Polarizable runs cost about 2–4× a fixed-charge run (more charges, SCF iterations or a smaller time step). They matter most where fields are strong or change a lot: ions and ion selectivity, interfaces (air–water, electrodes), ionic liquids and molten salts, and spectroscopy. A cheap mean-field shortcut for ionic liquids and ions is to scale the ionic charges by 0.7–0.8 (the "electronic continuum correction").

Key idea. Polarization is the response of the electron cloud to the local field: μ = αE, energy −½αE². It is many-body and short-ranged, needs damping at short range, and costs 2–4× more. Use it when the field around your molecules changes a lot.

The algorithm

One MD step with induced point dipoles (Drude version below):

E0_i = field at site i from all permanent charges (with PME)
μ_i = α_i E0_i                              # initial guess (or last step's dipoles)
repeat:                                      # self-consistent field (SCF)
    E_i = E0_i + Σ_j T_ij(damped) μ_j
    μ_i_new = α_i E_i
until max |μ_new - μ| < 1e-6 D (typically 5-10 iterations)
U_pol = -½ Σ_i μ_i · E0_i
forces = permanent + dipole terms; integrate as usual

Drude (extended Lagrangian, e.g. LAMMPS fix drude, GROMACS, OpenMM):
each polarizable atom = core (q - q_D) + Drude particle (q_D, mass ~0.4 u) on a spring k_D
integrate the centre of mass of each pair at T (normal thermostat)
integrate the core-Drude relative motion at T_D ~ 1 K (second, cold thermostat)
time step ~1 fs; add a hard wall (~0.2 Å) so a Drude cannot fly off

Use PME for the permanent and induced fields (codes extend PME to dipoles).
Exclude or screen interactions between dipoles on bonded neighbours (1–2, 1–3); Thole damping for the rest.
For SCF, a tight tolerance matters: a loose one leaks energy and makes NVE runs drift.
Software: AMOEBA in Tinker, Tinker-HP, OpenMM, AMBER; Drude in CHARMM, NAMD, GROMACS (drude branch), OpenMM and LAMMPS (fix drude, pair thole). Parameters: CHARMM Drude FF 2019, SWM4-NDP water, CL&Pol for ionic liquids.