Blog

Geometry of molecules

Sulfur dioxide in 3D: a yellow sulfur nucleus bent between two red oxygens, each joined by a pair of cylinders for the double bond, with a pale blue lone-pair lobe rising from the sulfur

Pick an element, and it appears. Select an atom, pick another element, and the new one arrives bonded to it. Every atom is a small nucleus inside a transparent shell, and every bond is a cylinder between two nuclei — steel for covalent, amber for ionic, one cylinder per bond order.

Pick an element to drop the first atom in, or load an example.

covalent ionic lone pairdrag to rotate · scroll to zoom · click an atom to select it

Nothing yet0 atoms · 0 bonds

Atoms you add will be listed here. Selecting one here is the same as clicking it in the view.

Add an atomunbonded

38 of the 118 have radii and oxidation states in the table this builder works from. The rest are greyed out — real elements, no numbers to draw them with.

Saved moleculesionic when the electronegativities differ by 1.8 or more

Nothing saved yet. Saves live in this browser only.

Nothing in there consults a table of bond angles. It started with two rules and it now has four, one fitted number, and a short list of places where it gives up and looks the answer up — all of them in the d block, and all of them labelled. This is an account of what each rule bought, which is a more interesting question than it sounds, because two of the things I added bought nothing and had to come back out.

The rules

One: bonds have a length. The distance between two nuclei is the sum of the two atoms’ covalent radii, which is where those radii come from in the first place — they were reverse-engineered from measured bond lengths precisely so that adding them up works. Carbon’s is 0.77 Å and hydrogen’s is 0.37 Å, so a C–H bond is 1.14 Å, and the real thing is 1.09 Å. Close enough to be recognisable, which is all this is trying to be. Multiple bonds are shorter, by factors chosen to land on the familiar carbon–carbon series: 1.54 Å single, 1.42 aromatic, 1.32 double, 1.20 triple.

Two: things attached to the same atom push each other away, with an inverse-square force, measured between directions rather than positions.

Three: a double or aromatic bond holds the groups at its two ends in one plane.

Four: a single bond prefers the groups at its two ends staggered.

Then the program rolls downhill: nudge everything along the forces, over and over, until nothing is moving.

Rules three and four are late arrivals and they are the smallest terms in the program by a wide margin. Rule two does nearly all the work, and it is worth seeing why before adding anything to it.

Why rule two is almost enough

Load methane and look at it. The four hydrogens sit at 109.47° from each other — the tetrahedral angle — and nobody told the program that number. Here is why they have to.

Fix a central atom and ask where its kk neighbours go. Bond lengths are already settled by the first rule, so each neighbour is pinned to a sphere around the centre and the only freedom left is direction. The second rule says those directions repel each other with an inverse-square force. Spreading kk points over a sphere to minimise inverse-square repulsion is a named problem — the Thomson problem, originally about electrons in Thomson’s plum-pudding atom — and for small kk its answers are exactly the shapes chemistry teaches:

NeighboursThomson solutionWhat chemistry calls it
2opposite poleslinear, 180°
3triangle on a great circletrigonal planar, 120°
4tetrahedrontetrahedral, 109.47°
5triangular bipyramidtrigonal bipyramidal
6octahedronoctahedral, 90°

So VSEPR, the rule every chemistry course teaches for predicting molecular shape, is the Thomson problem wearing a different hat. The playground is not approximating a table of ideal angles; it is solving the problem the table was an answer to.

The measuring-between-directions part matters more than it looks. The repulsion acts between unit vectors, so it is a pure function of angles and cannot stretch a bond however hard it pushes. Every later attempt to break that property caused a bug, and there is a whole section below about one of them.

What counts as a thing

The first version of this only let atoms push. Water came out straight, because an oxygen with two bonded neighbours puts them at opposite poles, and real water is bent at 104.5°.

The reason is famous: the oxygen is carrying two lone pairs, pairs of electrons that are in no bond but still take up room and still shove. Adding them needs no new rule. A lone pair is a thing attached to an atom, so rule two already covers it. All it needs is a position and a note that it pushes slightly harder than a bonding pair does, since it is held by one nucleus rather than shared between two.

Counting them is plain bookkeeping. Take the atom’s valence electrons, subtract its formal charge, subtract one for each bond order it is carrying, and halve what is left:

lone pairs=⌊valence−charge−∑bond orders2⌋\text{lone pairs} = \left\lfloor \frac{\text{valence} - \text{charge} - \sum \text{bond orders}}{2} \right\rfloor

Select any atom in the playground and it shows you that sum, so you can check it rather than take it on faith.

The bond-order term is not a detail. It is why bond order had to be built first. Carbon dioxide’s carbon has four valence electrons and two double bonds, so 4 − 4 = 0 lone pairs and the molecule stays straight. Count those same bonds as single and the carbon is left holding a pair it does not have, acquires a third thing to make room for, and bends to 120° — the model would have gone backwards on a molecule it already had right. You can do this to it deliberately: load carbon dioxide, select the carbon, and drop either bond back to single.

One number, fitted once

There is exactly one fitted constant in the whole program: how much harder a lone pair pushes than a bonding pair. I fitted it to a single measurement — water’s 104.5° — and got 1.2442.

Fitting one parameter to one number is not impressive on its own. What the parameter then does elsewhere is the whole test, because nothing was adjusted afterwards:

DomainsThe model saysMeasured
Water, H₂OAX₂E₂104.5°104.5° (fitted)
Ammonia, NH₃AX₃E₁106.7°106.7°
Ammonium, NH₄⁺AX₄109.5°109.5°
Xenon difluoride, XeF₂AX₂E₃linear, 180°linear, 180°
Xenon tetrafluoride, XeF₄AX₄E₂square planar, 90°square planar, 90°
Chlorine trifluoride, ClF₃AX₃E₂T-shaped, 85.5°T-shaped, 87.5°
Bromine pentafluoride, BrF₅AX₅E₁square pyramidal, 87.7°square pyramidal, 84.8°
Sulfur dioxide, SO₂AX₂E₁114.2°119.5°
Sulfur tetrafluoride, SF₄AX₄E₁see-saw, 123.5°see-saw, 102°

Ammonia is the one I would point at. Its angle was not fitted, its molecule was not consulted, and it lands on the measured value exactly.

There is a small bonus in there too. Textbook VSEPR has to state as a separate rule that lone-pair/lone-pair repulsion beats lone-pair/bond, which beats bond/bond. Here each pair’s shove is the product of the two weights, so 1.2442² > 1.2442 > 1, and that ordering comes free rather than being asserted.

The rule I tried to add and could not

The bottom two rows of that table are wrong, and the standard explanation for why is that a lone pair is not a point. It occupies a wider cone near the nucleus — it is fat — and in sulfur tetrafluoride it is that width, pressing sideways within the equatorial plane, that closes those fluorines from 123.5° down to 102°.

That is a satisfying story and it suggests an obvious fix: give lone pairs a width as well as a weight. So I built a standalone study of the problem — just points on a unit sphere, every way of seating nn bonds and mm lone pairs on the relevant Thomson solution, relaxed and scored — and fitted three versions of the model against the whole measured table.

ModelFree parametersBest fitError across the tableSF₄ equatorial
Weight only, fitted to water alone1λ = 1.24437.35°123.5°
Weight only, fitted to everything1λ = 1.07656.18°121.0°
Weight and a free exponent21.0869, 1.7446.19°121.0°
Weight and an angular width21.0764, 0.0005 rad6.18°121.0°

The width parameter fits to half a milliradian. Given the freedom to be fat, the model asks to be a point. And the best any of them manages is about a degree of improvement across the table, bought by pushing water off its measured 104.5° to 107.8° and dragging ammonia off the value it currently hits exactly — trading two right answers for a marginally better average.

Then I stopped fitting and checked whether SF₄ was reachable at all. Set the lone pair’s advantage to nothing, λ = 1, and its equatorial angle is 120.0° — a textbook trigonal bipyramid, which is what it should be with five identical things around an atom. Now turn the advantage up:

λ1.01.21.41.61.7
SF₄ equatorial120.0°122.8°126.5°132.4°140.3°

It opens. It only ever opens, monotonically, for every admissible value. The measured 102° is not merely unreached by my fitted parameter; it is on the wrong side of a bound that holds for the entire family of point-on-a-sphere models, including all the ones with extra parameters. A lone pair that pushes harder pushes those fluorines apart. Closing them needs a repulsion that is anisotropic — stronger within the equatorial plane than perpendicular to it — and a point charge on a sphere has no way to express a preferred plane.

So the fat-lone-pair rule is not in the program. The shipped model is the one-parameter version, which gets water and ammonia exactly right and SF₄ wrong by 21°, and I would rather have that than a two-parameter version that is wrong everywhere by less.

The third rule, and a piece of algebra that saved it

Rules one and two say nothing at all about twisting. Every angle around every individual atom can be perfect while the two halves of a molecule are rotated arbitrarily with respect to each other, because nothing connects them. Load ethene with rule three switched off and the two CH₂ groups settle fifty degrees out of plane, each locally perfect, jointly nonsense. Real ethene is flat and rigidly so; that is what the π bond is for.

The textbook way to write this down is a dihedral angle — the angle between the plane through the first three atoms and the plane through the last three — with an energy that is lowest when the two planes coincide. I wrote that, and it destroyed the program. Benzene came apart into a three-dimensional tangle, the relaxation stopped converging, and bond lengths drifted by a seventh of an ångström.

The cause is a genuine trap rather than a typo. The derivative of a dihedral angle carries a factor of 1/∣b1×b2∣21/|b_1 \times b_2|^2, which blows up when three of the four atoms line up — and during a relaxation they transiently do. Clamping the denominator only trades an explosion for a discontinuity in the energy, which the search then gets stuck against; I have the sweep showing the answer jumping around unpredictably as the clamp moves.

The escape is to notice what the rule is actually trying to say. Not “this dihedral angle is zero” but “these four atoms are coplanar” — and coplanarity has a clean algebraic measure. The volume of the box spanned by the three bond vectors is zero exactly when they lie in a plane, so with b1,b2,b3b_1, b_2, b_3 the three bonds along the chain:

Eflat=V(det⁡[b1,b2,b3]∣b1∣ ∣b2∣ ∣b3∣)2E_{\text{flat}} = V \left( \frac{\det[b_1, b_2, b_3]}{|b_1|\,|b_2|\,|b_3|} \right)^{2}

No cross products in a denominator, nothing to blow up, nothing to tune. And the two formulations are the same thing, because (b1×b2)×(b2×b3)=det⁡[b1,b2,b3] b2(b_1 \times b_2) \times (b_2 \times b_3) = \det[b_1,b_2,b_3]\, b_2, which makes that normalised volume equal to sin⁡(dihedral)sin⁡θ1sin⁡θ2\sin(\text{dihedral}) \sin\theta_1 \sin\theta_2 — a dihedral term that has already faded itself out near collinearity, which is precisely where a dihedral stops meaning anything.

The fourth rule, which cannot be dodged

Staggered ethane is the other thing nothing was enforcing. With the first three rules, ethane settles eclipsed, and cyclohexane’s chair is half flattened out.

The usual hand-waving explanation of the rotation barrier is that the hydrogens on opposite carbons get too close when eclipsed, so I tried the cheap version of that: let atoms three bonds apart push each other apart and let staggered win on its own. It does win — ethane snapped to a perfect 60° — and I threw it out anyway, because a repulsion between two atoms three bonds apart acts partly along the bonds between them. It put a fifth of an ångström into ethanol’s bond lengths and made cyclohexane’s chair worse. This is the same mistake as measuring rule two between positions instead of directions, and it is worth stating as a rule of thumb: a term that changes an answer it was not introduced to change is doing harm, however physical its story sounds.

There is also no shortcut available. Ethane has nine of these terms across its central bond, and because the three substituents at each end are 120° apart, a term with onefold or twofold periodicity in the dihedral contributes exactly the same total at every rotation — it sums to a constant and cannot prefer anything. Only a genuine threefold harmonic survives the sum. So the fourth rule has to be a real threefold torsion, which brings back the singularity from the third one.

The same trick works twice. Writing S=sin⁡θ1sin⁡θ2S = \sin\theta_1 \sin\theta_2 for the product of the two bond-angle sines, the two combinations Ssin⁡φS\sin\varphi and Scos⁡φS\cos\varphi are both free of small denominators — the first is the normalised volume above, and the second is just cos⁡θ1cos⁡θ2−cos⁡ψ\cos\theta_1\cos\theta_2 - \cos\psi, three dot products. And

S3(1+cos⁡3φ)=S3+4w3−3wS2,w=Scos⁡φS^3\left(1 + \cos 3\varphi\right) = S^3 + 4w^3 - 3wS^2, \qquad w = S\cos\varphi

which is the whole term, built from quantities that cannot explode. S2S^2 is a plain rational function of the positions; S3S^3 is its three-halves power, whose derivative carries a factor of SS and therefore dies quietly at collinearity instead of detonating there.

What the two new rules cost, measured: nothing. Every bond length, every bond angle, and every row of the AXₙEₘ table above is unchanged to the last digit printed. Ethene goes from 50° twisted to flat, ethane from eclipsed to 60°, cyclohexane’s chair from half-collapsed to a proper 60° everywhere, and benzene’s residual 2° ring twist to zero. That is the entire effect, which is what a well-behaved new term should look like.

Their strengths are not fitted, because there is no single measurement to fit them to, so both are set as low as will still do the job. The planarity one also has a ceiling that has nothing to do with chemistry: past about twice its current value the energy landscape gets stiff enough that the relaxer stops reliably finding benzene’s regular hexagon and starts settling for a planar ring with bonds alternating by a tenth of an ångström. A rule strong enough to beat the search is worse than a rule that is merely too gentle.

Ions are a different size

Table salt used to be drawn wrong, and wrong in a way that undercut the whole picture. Sodium’s van der Waals radius is 2.27 Å and chlorine’s is 1.75 Å, so the program drew a big sodium next to a small chlorine — which is correct for the two neutral atoms and backwards for the compound, where the chloride is unmistakably the larger of the two.

The fix is a second radius per element. Sodium gives an electron away and loses its outermost shell completely, ending up at 1.02 Å; chlorine takes one on and swells to 1.81 Å. Those are Shannon’s effective ionic radii, the numbers crystallographers use.

What makes this more than a table lookup is that the charges are not typed in. Marking a bond ionic is the statement that one electron has moved from the less electronegative atom to the more electronegative one, so the charges are derived from the bonds rather than stored alongside them, and the two can never disagree. Switch table salt’s bond to covalent and the charges vanish, the shells return to neutral size, and the bond shortens; switch it back and everything returns. The bookkeeping follows too — an ionic bond shares nothing, so it spends no electrons, and chloride correctly ends up holding a full four lone pairs where a covalent chlorine holds three.

The length changes as well. An ionic bond is two ions in contact, so it is measured with ionic radii: 2.83 Å for sodium chloride, which is the crystal figure. Covalent radii would give 2.53 Å, the length of a shared pair that in this case is not being shared.

The one thing it looks up

Four things around an atom make a tetrahedron. That is rule two, it is the four-point Thomson problem, and it is right nearly everywhere. It is wrong for copper(II), where four ligands sit in a square, for reasons that live in the d orbitals and are invisible to anything counting electron pairs.

They are invisible to this program too, and I want to be precise about how badly, because the obvious diagnosis is wrong. The square is not unstable here. Put two extra bodies above and below the plane — a square is an octahedron with its poles filled in — and the square becomes a perfectly respectable minimum: push a ligand a fifth of an ångström out of plane and it springs straight back. The problem is that it is not the lowest minimum. A lopsided alternative sits one part in six thousand below it, and the search, doing exactly what it is built to do, goes and finds it.

So copper’s square is placed, not predicted. The four ligands are put where copper is known to put them, and the two poles that hold them there are pinned — which stops the search looking for something better. The axis is an input, not a result, and pinning is the honest way to say so. Load copper(II) chloride and you are looking at the one shape in the playground that is a fact rather than a consequence.

Xenon tetrafluoride, in that same table above, is the uncomfortable neighbour of this. It is the same four-plus-two arrangement, it is genuinely square in reality, and the program draws it square — but from the same near-degeneracy, which means it is square because the search happens to start symmetric and never stumbles into the alternative. It is in the table as a success and it is really a coin landing the right way up. Measuring that gap was more useful than not knowing about it.

How many cobalts are there?

For a while the element list stopped at twenty-three, and somebody reasonably asked why — why no cobalt, given that nothing in the program cares what an element is called. And nothing does: the element lookup is a map, and the only place in the entire codebase that mentions a specific symbol is a cosmetic sort that puts carbon first in a formula. The list was short because every row is hand-entered, and the rows are not all the same kind of fact.

Cobalt turned out to be the best possible example of the problem, because the answer to “how big is a cobalt, and how many things does it gather around itself” is that there is no single answer:

Ionic radiusLigands it gathers
Cobalt(II)0.745 Åfour or six, both perfectly ordinary
Cobalt(III)0.545 Å, low spinsix, essentially always

Those are not two estimates of one number. They are a quarter apart, they are both right, and the table had one slot. So it had to either pick one and be wrong about the other, or not have cobalt in it. It had been quietly doing the first thing for iron and copper all along.

The fix is that an element now carries a list of oxidation states, each with its own radius, its own coordination number and its own opinion about whether that arrangement is square. What makes this cheap rather than a rewrite is that an oxidation state is not a new concept: it is a formal charge, atoms have carried one since charges were added, and it was already saved to disk and already derived from ionic bonds. Nothing needed a new field. The picker in the panel is only offering the charges that mean something for that element instead of any integer, and every lookup that used to ask the element now asks the charge first.

The payoff is a thing the program could not previously express at all. Load tetrachloroplatinate(II) and then hexachloroplatinate(IV). Same metal, same ligand, nothing different but the charge — and one is a square of four and the other an octahedron of six. No amount of counting electron pairs distinguishes those, because the electrons being counted are not the ones making the decision.

It also fixes a smaller dishonesty. An ionic bond to iron(III) is shorter than one to iron(II) by 0.135 Å, entirely because the ion is smaller, and the program now gets that right for the same reason it gets sodium chloride right.

The two limitations are worth stating rather than burying. Every radius here is Shannon’s at six-fold coordination, including for the states that are not six-coordinate — silver(I) is 1.15 Å at six and 0.67 Å at two, and quoting one column is a choice. It is the column everyone quotes and the only one that is complete, but a two-coordinate cation is drawn noticeably too large. And a coordination number is the commonest case rather than the only one. Nickel(II) is six-coordinate in most of its salts and square with four in the cyanide, so it carries both; beyond that pair there is no mechanism here for “it depends”, because there is no mechanism here for ligand field strength, which is the thing it actually depends on.

Why every metal preset wears bare chlorides

There is a reason diamminesilver(I) is the only complex here with a polyatomic ligand, and it is not taste.

I built hexaamminecobalt(III) first — Werner’s compound, the one that settled how coordination works, six ammonias on a cobalt. It came out 44° away from octahedral. So I tried harder: ten times the relaxation steps, five times the basin hops. The numbers came back identical to the digit. That is the tell that you are not looking at a search that gave up, you are looking at the model’s honest answer, and the model’s honest answer is wrong.

The cause is that the ammonia hydrogens on neighbouring ligands crowd each other, and the term that keeps atoms from overlapping is measured between positions, not directions. Rule two is careful to be a pure function of angles, which is exactly why it can set a shape without disturbing a bond length. The clash term has no such protection, and given six ligands with three hydrogens each it bends the L–M–L angles wholesale. It is the same failure that killed an earlier experiment of mine for the same reason, and I should have recognised it faster.

With bare halides — nothing to crowd — the same code lands on 90.0° and 180.0° exactly, for both the square and the octahedron. So the presets use bare halides, and cisplatin, whose whole point is that the shape is the drug, is described in a sentence rather than shipped as a wonky picture. Tetrachloroplatinate(II) is that same square with two of the chlorides still on it.

Rings, and the most useful thing I learned

Load benzene: a flat, regular hexagon, each ring bond drawn as the one-and-a-half it is counted as. Each carbon has three neighbours and no lone pair, and a flat hexagon is the one arrangement that gives all three of them 120° at every carbon simultaneously.

Then load cyclohexane, the same six-carbon ring with two hydrogens on each carbon instead of one. Now every carbon has four neighbours and wants 109.47°, but a flat hexagon forces 120° on the ring bonds. Something has to give, and what gives is the plane: the ring buckles into the shape chemists call a chair. It arrives there with no idea that it is famous for doing so.

Rings are also where rolling downhill stops being enough. Closing a ring yanks together two ends that were nowhere near each other, and the strain of that can settle into a twist where every atom is genuinely in equilibrium and the molecule is genuinely the wrong shape — benzene stuck at 115° in a shallow saddle. The fix is to shove the whole thing and roll downhill again, keeping the result only if it came out calmer. That is basin hopping, and it is what the tidy up button does. If something you have built looks wrong, press it.

Getting that working taught me the single most useful thing in this whole exercise. Basin hopping compares candidate shapes by scoring them, so the force the relaxer follows must be the exact gradient of the score it is judged by. My first attempt at the repulsion projected the force sideways by hand, to stop it stretching bonds, which worked and was not the derivative of anything. The relaxer then cheerfully walked away from arrangements its own scoring called better, and benzene never flattened.

I now have a test that differentiates the score numerically and compares it, component by component, against every force the program reports — across a hundred awkward arrangements with lone pairs, mixed bond orders, crowding, rings and near-collinear chains. When I added the third rule, that test caught the error on its first run: my hand-written dihedral derivative had two sign conventions mixed together, and the failures pointed straight at the two middle atoms. Working out what the derivative should be took twenty minutes of measuring it numerically and reading off the coefficients. Without the test I would have spent the afternoon looking at benzene and wondering.

Ionic or covalent

The colour of a bond is decided by electronegativity — how hard an atom pulls on shared electrons. When the two atoms differ by less than 1.8 on the Pauling scale, the electrons are shared and the bond is drawn in steel. Above that, one atom has effectively taken them, and the bond is drawn thinner and in amber. Sodium and chlorine differ by 2.23, so table salt comes out ionic; oxygen and hydrogen differ by 1.24, so water does not.

That cutoff is a convention rather than a fact — real bonds shade continuously from one kind to the other, and 1.8 is just where textbooks like to draw the line. Select any atom and its bonds are listed with the actual Δχ, so you can see how close to the line each one sits, and override any of them by hand. Overriding is not paint any more, though: because the charges and the radii follow from the bond type, forcing a bond ionic really does turn both atoms into ions and change the picture.

Where it is still wrong

Sulfur tetrafluoride and sulfur dioxide, for the reason measured above: a lone pair here is a point with a weight, and what those two molecules need is a lone pair with a preferred plane. This is a limit of the whole approach rather than a missing constant, which is a more comfortable thing to know than to suspect.

Copper is told, not derived, and xenon tetrafluoride is right by luck. Both are the same near-degeneracy. Every transition metal’s coordination number is a lookup too, and now so is every oxidation state’s.

Crowded polyatomic ligands come out bent, measurably and irreparably, for the reason above: the clash term works on positions and rule two works on directions. Anything with six ammonias on it is beyond this program, and more computing time does not touch it.

One weight cannot be right at every coordination number. A lone pair’s advantage over a bonding pair is not the same fraction when there are three things around an atom as when there are four. The weight was fitted at four, and sulfur dioxide, with three, pays for it.

Nothing here is an energy. The score has no units and no temperature; it is a shape-finding device, not a force field. It will not tell you a barrier height or which of two conformers you would actually find in a bottle.

There is no solvent, no crystal and no neighbours. Every molecule here is alone in a vacuum, which is why the ionic bond length is quoted against the crystal figure and not the gas-phase one.

Things worth trying

If you like watching simple rules produce structure, the game behind TREE(3) is two rules about coloured trees that produce a number too large to write down.

Comments