Unhandled xTB errors such as NaN energy, incorrect net charge, and huge partial charges

[jminuse]

Describe the bug

For some inputs, xtb-python produces impossible outputs which are not handled as errors.

To Reproduce

This is a script to demonstrate unhandled xTB errors for simple cubes.
The command line arguments are atomic number and grid spacing (in Angstroms) for a cube of the given element.
Metals tend to be more unstable.

import sys
import numpy as np
from scipy import constants
import xtb
from xtb.interface import Calculator, XTBException
from xtb.libxtb import VERBOSITY_MUTED
from xtb.utils import get_method
ANGSTROM_IN_BOHR = constants.physical_constants['Bohr radius'][0] * 1.0e10


def run_xtb(atomic_numbers, xyz, net_charge, n_unpaired_electrons, max_iter=250, conv=1.0, method="GFN2-xTB"):
    '''
    Run a single xtb job and return the result.
    Throws XTBException for errors or unrealistic results.

    Arguments:
        atomic_numbers (np.array): integer element labels for all atoms
        xyz (np.array): Nx3 array of coordinates for all atoms, in Angstroms
        n_unpaired_electrons (int): number of unpaired electrons in the system
        max_iter (int): max number of SCF iterations before giving up
        conv (float): SCF cutoffs, smaller means tighter convergence
    Returns:
        E_xtb: total system energy in Hartree
        gradient_xtb: dE/dxyz for each atom, in Hartree/Angstrom
        dipole_xtb: system dipole in electron-Angstrom (not available with PBC)
        qs_xtb: Mulliken atomic charges, in elementary charge units
    '''
    calc = Calculator(get_method(method), atomic_numbers, xyz / ANGSTROM_IN_BOHR,
                      charge=net_charge, uhf=n_unpaired_electrons)
    calc.set_verbosity(VERBOSITY_MUTED)
    calc.set_accuracy(conv)  # default = 1.0, smaller = tighter convergence
    calc.set_max_iterations(max_iter)  # default = 250
    result = calc.singlepoint()
    if result.check() != 0:
        raise XTBException(result.get_error())
    E_xtb = result.get_energy()
    if np.isnan(E_xtb):  # yet another convergence issue
        raise XTBException('Energy is NaN')
    qs_xtb = result.get_charges()
    if abs(net_charge - np.sum(qs_xtb)) > 1e-5:
        raise XTBException('Bad net charge: %d vs %f' % (net_charge, np.sum(qs_xtb)))
    if np.max(np.abs(qs_xtb)) > 4.0:
        raise XTBException('Very implausible partial charge: %f' % np.max(np.abs(qs_xtb)))
    dipole_xtb = result.get_dipole() * ANGSTROM_IN_BOHR  # from electron-Bohr to electron-Angstrom
    gradient_xtb = result.get_gradient() / ANGSTROM_IN_BOHR  # from Hartree/Bohr to Hartree/Angstrom
    # units: Hartree, Hartree/Angstrom, electron-Angstrom, electrons
    return E_xtb, gradient_xtb, dipole_xtb, qs_xtb

def test_metal_cubes(atomic_number, default_scale=3.0, cube_size=3):
    print('Cube of atomic number', atomic_number, 'with', cube_size**3, 'atoms, default interatomic r =', default_scale)
    for scale in np.arange(0.2, 1.2, 0.01):
        scale *= default_scale
        xyz = []
        for i in range(cube_size):
            for j in range(cube_size):
                for k in range(cube_size):
                    xyz.append((i, j, k))
        xyz = np.array(xyz, dtype=np.float64) * scale
        atomic_numbers = np.array([atomic_number] * len(xyz))
        q = np.sum(atomic_numbers) % 2  # if odd n electrons, remove 1
        try:
            E_xtb, gradient_xtb, dipole_xtb, qs_xtb = run_xtb(atomic_numbers, xyz, q, 0)
            # print('r = %.2f' % scale, 'E_xtb =', E_xtb)
        except XTBException as ex:
            print('r = %.2f' % scale, 'Failed,', ex)

print('xTB version', xtb.__version__)
atomic_number, default_scale = int(sys.argv[1]), float(sys.argv[2])
test_metal_cubes(atomic_number, default_scale)

Example input:
python xtb_huge_error.py 3 2.7
Output:

xTB version 20.1
Cube of atomic number 3 with 27 atoms, default interatomic r = 2.7
r = 0.92 Failed, Very implausible partial charge: 209.337257
r = 0.95 Failed, Very implausible partial charge: 96.060772
r = 0.97 Failed, Very implausible partial charge: 66.404307
r = 1.00 Failed, Very implausible partial charge: 163.766930
r = 1.05 Failed, Very implausible partial charge: 145.753860
r = 1.08 Failed, Bad net charge: 1 vs 1.000011
r = 1.13 Failed, Very implausible partial charge: 29.071110
r = 1.19 Failed, Very implausible partial charge: 257.358086
r = 1.22 Failed, Very implausible partial charge: 273.118857
r = 1.24 Failed, Bad net charge: 1 vs 0.412483
r = 1.27 Failed, Very implausible partial charge: 115.909595
r = 1.30 Failed, Very implausible partial charge: 19332.830109
r = 1.35 Failed, Very implausible partial charge: 5.019581
r = 1.38 Failed, Very implausible partial charge: 33.592007
r = 1.40 Failed, Very implausible partial charge: 33.003970
r = 1.43 Failed, Very implausible partial charge: 66.802335
r = 1.46 Failed, Bad net charge: 1 vs 1.000011
r = 1.49 Failed, Very implausible partial charge: 91.379054
r = 1.51 Failed, Bad net charge: 1 vs 1471240208384.000000
r = 1.76 Failed, Energy is NaN
r = 2.11 Failed, Energy is NaN

Expected behaviour

These error conditions should be handled within xtb or xtb-python, rather than requiring user checks. Probably the most alarming error here is Bad net charge: 1 vs 0.412483, where a significant fraction of an electron has appeared from nowhere.

Additional context

GFN1-xTB and GFN0-xTB don't have the same problems with bad net charges or implausible partial charges, but they can still give NaN for energy.

[jminuse]

Update: with a clean install of the latest xtb-python, this script gives the output "Failed, Energy is NaN" only for GFN1-xTB and GFN0-xTB. GFN2-xTB now raises the appropriate exception, "scf: Self consistent charge iterator did not converge", and none of the methods give an incorrect net charge. However, it would still be nice to have appropriate handling of a NaN energy (possibly related to #80).

[jminuse]

Another example: with the latest xtb-python, the following system returns non-finite values for the gradient in GFN0-xTB & GFN2-xTB. With GFN0-xTB, the gradient is an array of mixed np.inf and -np.inf, while with GFN2-xTB, the gradient is all NaN.

atomic_numbers = np.array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3])
xyz = np.array([[ 0.305,  1.818,  2.417],
       [-1.227,  0.658, -2.75 ],
       [ 0.021,  2.861,  0.379],
       [-1.928,  1.384,  0.584],
       [-1.211,  2.915, -1.727],
       [ 0.72 ,  1.346, -1.534],
       [-1.089, -0.193,  2.436],
       [ 1.151, -0.511, -2.935],
       [-0.115, -0.055,  0.266],
       [-1.341, -2.106,  1.163],
       [-2.026, -0.473, -0.915],
       [-0.122, -2.054, -1.216],
       [ 1.956,  1.279,  0.876],
       [ 3.207,  1.115, -1.338],
       [ 1.125, -0.637,  2.324],
       [ 3.241, -0.997,  1.191],
       [ 1.39 , -2.447,  0.827],
       [ 1.987, -0.837, -0.745],
       [-3.279, -0.566,  1.313],
       [-2.764, -2.499, -0.615]])
charge = 0
uhf = 0

Energy, charge, and dipole all appear perfectly normal. With GFN1-xTB, or the DFT functional wB97X-D3BJ/def2-TZVPD in Psi4, this system is fine, so I don't think there's any physical reason for it to give invalid gradients. In any case, I would think it should throw an exception rather than silently return non-finite values.