This project stands on software other people wrote and give away for free. Almost all of our quantum chemistry runs on PySCF, a Python library kept up by a group of researchers and volunteers. Next to it we use a handful of smaller packages: one for anharmonic corrections, one for geometry optimisation, one for a fast approximate method.

Software that many people use still has corners that nobody has walked into yet. Our calculations are large and unusual enough that we walked into several. Each time we had two choices: work around the problem in our own code, or fix it at the source so that the next person does not hit it. Where we could, we fixed it at the source. Each fix goes to the maintainers as a pull request: a proposed change, usually with a test that fails without the change and passes with it. The maintainers read it, ask questions, and decide whether to take it.

One thing should be said plainly. The code in these pull requests was written and tested by an AI coding assistant, Claude Code, which works on this project with me. I decide what we take on and hand in the pull requests; tracking down the cause, writing the change and proving it with tests is its work.

So far that has been nineteen pull requests to six projects. Eleven have been accepted, six are waiting for review (four of them submitted on 27 October), and two we withdrew because a better fix came along. Below is each of them, in a sentence or two.

PySCF

PySCF is a free Python library for quantum chemistry, written and kept up by researchers and volunteers; almost every calculation in this project runs on it. Its documentation is at pyscf.org.

  • #3469, accepted 6 October. When you ask PySCF for the gradient of a CCSD(T) energy without first computing one particular intermediate, it quietly used the intermediate of the simpler CCSD method. The result looked fine but was the slope of a different energy. We lost six days of computing before we noticed; now PySCF computes the right intermediate by itself.
  • #3470, accepted 6 October. Two of the slowest steps of the CCSD(T) gradient were written in Python. We rewrote them in C, which made them 16 and 14 times faster on benzene.
  • #3477, accepted 12 October. On its way to the gradient, the calculation builds the electron density, and from that density the dipole moment follows almost for free. PySCF threw the density away. PySCF now keeps it, so every gradient also gives the dipole moment, which decides how bright an infrared line is.
  • #3507, accepted 16 October. A computed Hessian should be symmetric, but small numerical errors make its two halves differ slightly. PySCF’s frequency routine read only one half, so the answer depended on which half happened to carry the error. It now averages the two halves first.
  • #3531, accepted 22 October. Inside a Hessian calculation a large set of equations is solved step by step (CPHF). That solver never wrote its progress to the log, so a run that was stuck looked exactly like a run that was working. It now reports every step. We had first reported the problem as #3530.
  • #3508, withdrawn. The same CPHF solver sometimes stops before it has reached the accuracy it was asked for, a problem another user had reported in #2702. We added measurements to that report and proposed a fix, then withdrew it in favour of a fix by another contributor that repairs the solver itself, for every part of PySCF that uses it.

PySCF-forge (PySCF’s add-ons)

PySCF-forge holds new methods for PySCF until they are ready to move into PySCF itself.

  • #212, accepted 30 September. The add-on for local correlation methods stopped working with the newest versions of PySCF and of NumPy, the number-crunching library underneath. Two small changes made it work again.
  • #213, accepted 9 October. When you install the add-ons, the build looks up where PySCF lives. It could ask a different Python than the one doing the installing, and then build against the wrong copy of PySCF. It now asks the right one.

PySCF properties

PySCF properties computes how molecules and crystals respond to electric and magnetic fields, among them the brightness of infrared lines.

  • #16, waiting for review since 11 October. The infrared routine crashed when you handed it your own orbitals instead of letting it compute them; only the default route worked.
  • #17, waiting for review since 11 October. The existing tests of the infrared intensities compare with a commercial program. This PR adds a test that needs no other program: it checks the routine against finite differences, small steps of the atoms with the dipole moment computed at each step.

pyVPT2 (anharmonic corrections)

pyVPT2 computes the anharmonic corrections to vibrational frequencies by second-order perturbation theory, on top of the quantum-chemistry package Psi4 (documentation).

  • #59, accepted 19 October. pyVPT2 has a built-in check that is meant to warn when its numbers are noisy. The check compared two numbers that were equal by construction, so it could never fire. It now compares two independent routes to the same number. This PR replaced our first attempt, #58, which only reported the problem; the report itself is #57.
  • #61, accepted 19 October. When two vibrations share almost the same energy they mix (a Fermi resonance), and such groups have to be solved together. The solver could crash, or treat one connected group as two separate ones. It now merges them. Report: #60.

optking (geometry optimisation)

optking searches for the shape in which a molecule is most stable, a geometry optimisation. It started as part of Psi4 and is now a separate Python program.

  • #116, accepted 13 October. An optimisation could hang forever, using a full processor core and printing nothing, when an angle close to zero was mistaken for a straight line. optking now recognises that case and starts again cleanly. Report: #115.

dxtb (a fast approximate method)

dxtb, from Stefan Grimme’s group in Bonn, is a fast approximate quantum-chemistry method (extended tight binding) written in PyTorch, so that it can be combined with machine learning (documentation).

  • #278, accepted 23 October. The forces on the atoms came out wrong when the calculator reused results it had stored for an earlier geometry. It now clears those stored results when it gets a new geometry. Report: #277.

Submitted to PySCF, waiting for review

We submitted these four on the morning of 27 October, after reworking all four to follow the style rules that PySCF publishes for contributors. None of them has been reviewed yet, so they may still change or be turned down. The first two fix problems that can hit any molecule. The other two are about cations, molecules that have lost one electron. Many PAHs in space are cations, and the most accurate method we use, CCSD(T), ran out of memory on every cation beyond the smallest.

  • #3545, empty block. The gradient code works through the basis functions in blocks. When an atom’s first group of functions was wider than one block, the first block came out empty and the calculation crashed. We met it on benzene with a large basis set.
  • #3546, CCSD(T) tables on disk. PySCF has a route that is meant to keep the biggest tables on disk instead of in memory. For the CCSD(T) gradient it still built three of them in memory, 28 GB each for benzene with a large basis set. They now go to disk.
  • #3547, the same for cations. The cation version of the same calculation had the same problem. On a small test molecule the memory peak drops from 4.1 to 3.4 GB.
  • #3548, CCSD(T) for cations in pieces. The cation CCSD(T) gradient built several tables of about 22 GB each for the benzene cation, and 830 GB each for a three-ring cation. The work is now done in blocks small enough to fit. The whole gradient of the benzene cation ran on a laptop with a memory peak of 6.1 GB.