Correspondence between dxtb and xtb-python of xtb-gfn1

[hannesvdc]

Dear Grimme lab

We're evaluating dxtb as a GPU-accelerated drop-in for the CPU xtb-python calculator. We found correspondence is excellent overall on a range of molecules between 2 and 182 atoms. GFN1 energies agree to ~1e-6 Ha and forces to <= 2e-5 eV/Å in float64 on CPU.

But at a couple of peculiar cases, the forces disagree, and I wanted to ask whether you've seen this before.

Test case

At the equilibrium (ASE G2) geometries of methanol ( minimal code example below), dxtb and xtb-python GFN1 forces differ by approximately 0.22 eV/Å on the x-coordinate of the carbon atom. Energies still agree within a relative error of 1e-5 eV.

dxtb is correct because the negative central-finite-difference gradient of xtb-python's own energy matches dxtb to ~1e-6 eV/Å, and
disagrees with xtb-python's analytical force by the full 0.225 eV/Å. This totally makes sense because dxtb forces are computed through automatic differentiation. Apparently, xtb-python's reported force is not the gradient of its own energy there, while
dxtb's autograd force is.

The discrepancy is invariant to SCF accuracy (1.0→1e-4), electronic temperature, and finite-difference step size. The issue does not appear to be a convergence artifact. It's also geometry-localized: a non-equilibrium methanol conformer agrees perfectly.

Numbers (methanol equilibrium, worst atom = C), eV/Å

  xtb-python analytical : [-0.0652 -0.3899  0.    ]
  dxtb autograd         : [-0.2904 -0.3899  0.    ]
  -FD(xtb-python energy): [-0.2904 -0.3899  0.    ]   (neutral truth)

  max |xtb-python analytical - (-FD(xtb energy))| = 2.25e-01 eV/A
  max |dxtb autograd        - (-FD(xtb energy))| = 1.6e-06 eV/A

Is this a known incompatibility between dxtb and the xtb/tblite analytical gradients at (near-)symmetric geometries? Have you encountered it in other situations or molecules, and is there a known cause (e.g. a gradient term that xtb/tblite treats differently, or a documented edge case)? We'd like to know how broadly it can occur before relying on either code's forces near symmetric
stationary points.

Have you also seen this issue with xtb-gfn2?

Minimal-reproducing script:

import numpy as np, torch
from ase import Atoms
from xtb.ase.calculator import XTB
from dxtb_potential import DxtbPotential   # thin dxtb wrapper, GFN1, f64

SYMBOLS = ["C","O","H","H","H","H"]; Z = np.array([6,8,1,1,1,1])
XYZ = np.array([[-0.047131, 0.664389, 0.0],[-0.047131,-0.758551,0.0],
                [-1.092995, 0.969785, 0.0],[ 0.878534,-1.048458,0.0],
                [ 0.437145, 1.080376, 0.891772],[0.437145,1.080376,-0.891772]])

def E(x):
    a=Atoms(symbols=SYMBOLS,positions=x)
    a.calc=XTB(method="GFN1-xTB",electronic_temperature=300.0); return a.get_potential_energy()

a=Atoms(symbols=SYMBOLS,positions=XYZ); a.calc=XTB(method="GFN1-xTB",electronic_temperature=300.0)
F_xtb = a.get_forces()
_, F_dxtb = DxtbPotential(Z=Z, method="GFN1-xTB", dtype=torch.double).energy_forces(XYZ)

h=2e-4; g=np.zeros_like(XYZ)
for i in range(6):
    for j in range(3):
        xp=XYZ.copy(); xp[i,j]+=h; xm=XYZ.copy(); xm[i,j]-=h
        g[i,j]=(E(xp)-E(xm))/(2*h)
F_fd = -g
print("xtb analytical  - FD:", np.abs(F_xtb  - F_fd).max())   # ~0.225
print("dxtb autograd   - FD:", np.abs(F_dxtb - F_fd).max())   # ~1e-6

The full code for this example is in our repo https://github.com/hannesvdc/ChemDM/blob/main/examples/xtb_gpu/reference_gradient_bug.py

Thank you,

Hannes Vandecasteele

[marvinfriede]

Hi, I don't know of any similar issues. From your results, this seems to be an issue with xtb-python and not with dxtb. Please note that xtb-python is not supported anymore. I would highly encourage you to switch over to tblite. It is under active development, bugs get fixed, and it already has a dedicated ASE calculator.

Just a word of caution regarding dxtb and GPU usage: While the code works on GPU, we did not spend much time optimizing for GPU performance, because we did not have access to state-of-the-art GPUs back then. If you are planning on going to very large systems (~10.000 atoms), you will probably run out of GPU memory sooner or later (especially if you need forces).

[hannesvdc]

Hello Marvin

Thank you for your quick reply! Yes, the force field discrepancy indeed seems like an issue with xtb-python.

I'm mainly interested in replacing xtb-python with dxtb wholesale, especially for (reduced) Hessian calculations. But I will also check out tblite as you suggested.

Regarding GPU acceleration, to what extent is dxtb fully batched? From a look at the source, the SCF itself seems genuinely batched on CUDA (batched Fock build and a single eighb/torch.linalg.eigh over the [B, #atoms, #atoms] stack). But the integral construction via libcint appears to be evaluated sequentially per molecule and forced onto the CPU (force_cpu_for_libcint).

Is this a correct reading of the code? Does the CPU step (transfer and compute) become a bottleneck in your opinion? I know you haven't done a full performance check yet, but do you have an intuition for which batch size dxtb overtakes tblite?

Thank you for the help, Hannes

[marvinfriede]

Regarding GPU acceleration, to what extent is dxtb fully batched? From a look at the source, the SCF itself seems genuinely batched on CUDA (batched Fock build and a single eighb/torch.linalg.eigh over the [B, #atoms, #atoms] stack). But the integral construction via libcint appears to be evaluated sequentially per molecule and forced onto the CPU (force_cpu_for_libcint).

That is correct. Integral evaluation is very difficult (if not impossible) to vectorize, and hence, difficult to make performant within PyTorch (i.e., without writing custom CUDA kernels). The libcint evaluation is significantly faster, also for gradients.

Does the CPU step (transfer and compute) become a bottleneck in your opinion?

I don't think so. Compared to MLIPs and FFs, the diagonalization in the xTB methods is much slower, and will probably always be the bottleneck.

I know you haven't done a full performance check yet, but do you have an intuition for which batch size dxtb overtakes tblite?

Hard to say. This also depends on many parameters.
What are you comparing (energy, gradients, ...)? Analytical gradients will always be faster with tblite than autodiff gradients. That is just the price to pay for not having to implement analytical derivatives.
What molecules are batched? If you have different sizes, there will be some overhead due to padding.
How large are the systems being batched? In general, batching is more beneficial for small molecules (let's say up to 50 atoms), because you save setup overhead and parallelization is more effective for larger systems.
There are probably even more details influencing batch performance. For example, if there is one system in the batch with poor SCF convergence, it could make batching less efficient because the other systems are needlessly dragged through the SCF too.