Lab: electrostatics by hand

The setup

We stay in GROMACS units. One new constant: the Coulomb energy of two unit charges at distance r is f/r with

For comparison, k_B T at 300 K is 2.494 kJ/mol: an ion pair at contact is about 200 k_B T. Electrostatics is the largest energy in an ionic system, so small relative errors in it are large errors in everything else.

Step 1

Why Coulomb is different from Lennard-Jones

Why: In Lab 2 we cut LJ off at 2.5σ and lost about 3 % of the energy. Before doing the same with charges, compare how the two fall off with distance, and how many neighbours sit at each distance.

What the numbers say: LJ falls as 1/r⁶ while the number of neighbours grows as r², so far shells contribute almost nothing: a cutoff is safe. Coulomb falls only as 1/r, so a far shell holds more energy than a near one, 10 000 kJ/mol at 2 nm, unless its positive and negative charges cancel almost perfectly. The whole electrostatic energy is a small difference between huge positive and negative pieces. Any method that treats far shells carelessly, such as a simple cutoff, can be wrong by more than the answer itself.

Step 2

Truncation and the order of summation

Why: To see the problem cleanly, take crystals whose exact energy is known. A 1D chain + − + − with spacing a has the energy of one ion equal to −(2f/a) × (1 − ½ + ⅓ − ¼ + …) = −(2f/a) ln 2. Add the neighbours one by one, as a cutoff would.

In 3D the rock-salt (NaCl) crystal has Madelung constant 1.747565. Add up all ions around one Na⁺ inside a growing sphere, or inside a growing cube (ions on the cube faces counted ½, edges ¼, corners ⅛):

What the numbers say: In 1D, cutting off after n neighbours leaves an error of about 1/(2n), and its sign flips each time a new ion is added: with 100 neighbours per side it is still 0.7 % off. In 3D it gets worse: a sphere never converges at all (0.07, 3.16, −1.21, 1.89, …) because each new spherical shell is not neutral. The cube, which keeps every layer neutral, settles on 1.7476. The sum depends on the order of addition, which mathematicians call conditionally convergent. The lesson: never cut Coulomb off bluntly. Either make every piece neutral, or use a method that sums the infinite periodic system exactly: Ewald.

Step 3

Two ions in a periodic box: what the minimum image misses

Why: In a periodic simulation, each ion interacts with the other ion and with all of their images, an infinite lattice of ion pairs. The minimum image (Lab 2) keeps only the nearest copy. For LJ that is enough; let us see what it misses for charges.

What the numbers say: At 0.3 nm the images add −29.8 kJ/mol (6 %, 12 k_B T). The error grows as the ions separate, to −103 kJ/mol at 0.5 nm, where the Cl⁻ sits halfway between the Na⁺ and its image and is pulled by both, but the minimum image sees only one. So the minimum image gives the wrong forces too: it would push the pair together. To get the periodic energy we need all images, and Step 2 showed we cannot just add them up one by one.

Step 4

Ewald summation by hand

Why: Ewald's trick: split 1/r into a short-range part erfc(αr)/r, which dies off fast and can be summed in real space like LJ, and a smooth long-range part erf(αr)/r, which is summed as waves (k-vectors) in reciprocal space, where it converges fast. A self term removes the interaction of each charge with its own smoothed cloud. α decides where the split happens.

What the numbers say: Each piece is ordinary arithmetic: one erfc, one square root, and for each k-vector an exponential and a cosine. The three parts change enormously with α (the self term from −314 to −1254), but their sum stays at −492.9 kJ/mol whenever both sums are converged, which is the check that Ewald is right. The limits show up at the edges: with α = 2 the real-space part needs images beyond the minimum image (error +7.5); with α = 8 the k-sum needs more than 728 vectors (error −5.4). Choosing α is balancing the two costs; codes pick it from the cutoff and a tolerance (GROMACS ewald-rtol, LAMMPS kspace_style ... 1e-5).

Step 5

A charged box: the background and its energy

Why: Ewald needs the box to be neutral; otherwise the infinite lattice of identical charges has infinite energy. Codes then silently add a uniform neutralising background. This is common when one computes the solvation free energy of a single ion, and the background has an energy that depends on the box size.

What the numbers say: An artificial energy of −99 kJ/mol in a 2 nm box (40 k_B T), and it falls only as 1/L: making the box 5 times bigger still leaves −39 kJ/mol. For a charged system, any energy that is compared between box sizes or between charge states, such as an ion solvation free energy or a redox potential, must be corrected analytically for this term (and its dielectric screening). For a neutral electrolyte it cancels, which is why simulations should be neutral whenever possible.

Step 6

The cost: why PME

Why: Ewald is exact but its reciprocal sum loops over every k-vector for every charge. How does that scale with the number of atoms N, compared with adding all pairs, or with particle-mesh Ewald (PME), which spreads the charges on a grid and uses the fast Fourier transform?

What the numbers say: For 100 000 atoms, PME needs about 20 times fewer operations than optimally tuned Ewald and 3000 times fewer than all pairs (which would not even include the images). That is why every modern code uses PME or a relative (PPPM in LAMMPS). Its price is an interpolation error from the grid, set by the grid spacing (~0.12 nm) and the interpolation order; it is the same α trade-off as in step 4, now with a grid. Limitation: the FFT needs all-to-all communication, which limits scaling on many nodes.

Step 7

Induced dipoles: polarization by iteration

Why: Real atoms are not fixed charges: an electric field shifts their electron cloud and creates a dipole μ = α E, which in turn changes the field on its neighbours. Polarizable force fields include this, and it matters most next to ions, where fields are strong. We compute it, then find where it fails.

What the numbers say: One induced dipole is worth −27 kJ/mol next to an ion, about 10 k_B T: not a small correction. With two atoms, the mutual induction raises μ_B from 0.006 to 0.011, almost double, and needs about 5 iterations at every time step (or extra degrees of freedom, as in Drude oscillators). Limitation, by numbers: if two polarizable sites come closer than 0.144 nm, each iteration amplifies the dipoles instead of correcting them and the energy goes to −∞. Real polarizable force fields therefore damp the interaction at short range (Thole screening). The cheap alternative many non-polarizable models for ionic liquids use is to scale all ionic charges by about 0.7–0.9, an average version of the same screening.

Summary