Charge not passed for GFN-FF API call

[pavankum]

Hi, I encountered a problem with xtb python interface for a gfnff calculation that raised EEQ charge constrain error for a charged molecule. The same works well with direct invocation of xtb from command line so it's a bit puzzling for me. I have the following versions in my env

xtb                       6.4.1                hf06ca72_1    conda-forge
xtb-python                20.2             py37h5e8e339_3    conda-forge

I attached the xyz file I used to run the calculation via command line using xtb --gfnff xtb_mol.xyz -charge -1, this works well and the output is in xtb_from_CLI.out.

The python code to reproduce the behaviour is here, and the output is in xtb_:

from xtb.interface import Calculator
from xtb.utils import get_method
import numpy as np

positions = np.array([[ 1.60661680e+00, -3.54849264e-01,  1.40573187e+00],
       [ 3.74941117e-01, -5.45532477e-01,  7.78876693e-01],
       [ 1.72053005e+00, -6.24272886e-01,  2.74744734e+00],
       [-6.64451602e-01, -1.00863329e+00,  1.56734114e+00],
       [ 6.04983698e-01, -1.07927167e+00,  3.49672596e+00],
       [-1.90987134e+00, -1.64422697e+00,  6.90357388e+00],
       [-8.84291779e-01, -2.16110423e+00,  7.62236059e+00],
       [-1.77323096e-01, -1.54235367e+00,  5.64622607e+00],
       [-5.82876394e-01, -1.27257602e+00,  2.87036935e+00],
       [-1.46803934e+00, -1.25750181e+00,  5.66181125e+00],
       [ 8.19764156e-01, -1.34396241e+00,  4.79452364e+00],
       [ 2.32357934e-01, -2.11295396e+00,  6.84982814e+00],
       [ 2.46506874e+00, -1.69232048e-03,  8.48547597e-01],
       [ 2.27953406e-01, -3.49412789e-01, -2.73016559e-01],
       [ 2.65602892e+00, -4.95349741e-01,  3.27334037e+00],
       [-1.63897466e+00, -1.18578591e+00,  1.12017610e+00],
       [-2.93673078e+00, -1.52735109e+00,  7.20417879e+00],
       [-7.84421837e-01, -2.57471112e+00,  8.60801640e+00]])
numbers = np.array([6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 8, 1, 1, 1, 1, 1, 1])

calc = Calculator(get_method("GFNFF"), numbers, positions, charge = -1)
res = calc.singlepoint()
res.get_energy()

The traceback for the error:

Traceback (most recent call last):
  File "xtb_interface.py", line 26, in <module>
    res = calc.singlepoint()
  File "/home/maverick/anaconda3/envs/qcarchive-worker-openff-xtb/lib/python3.7/site-packages/xtb/interface.py", line 814, in singlepoint
    raise XTBException(self.get_error("Single point calculation failed"))
xtb.interface.XTBException: Single point calculation failed:
-2- gfnff_calculator_singlepoint: Force-field method terminated
-1- gfnff_eg: EEQ charge constrain error

xtb_inspection.zip

Please let me know if there is a mistake in using the python interface. Thank you for your time.

[awvwgk]

Currently, GFN-FF is ignoring the molecular charge set in the geometry input, by passing a global variable, ichrg, rather than the actual total system charge, nint(mol%chrg), at:

xtb/src/gfnff/gfnff_setup.f90

Line 80 in 79825ab

call gfnff_ini(env,verbose,ini,mol,ichrg,gen,param,topo,accuracy)

This is not a problem in the standalone, since the global variable ichrg is always in sync with the total system charge:

xtb/src/prog/main.F90

Line 336 in 79825ab

ichrg = nint(mol%chrg)

But breaks in the API since the provided input data is not processed correctly in the gfnff_setup routine.

[awvwgk]

cc @Thomas3R

[pavankum]

Thank you so much for the clarification @awvwgk!!

[awvwgk]

I proposed a patch, which will change the behaviour of GFN-FF to respect the total charge provided with the input molecule. Open for comments.

[pavankum]

Thanks for putting in the fix for this. I tried to test this by pulling down your branch, building the xtb executable and xtb-python extensions. xtb executable shows the correct commit * xtb version 6.4.1 (a0a5b9d). I followed the instructions here to build xtb-python, I copied all the _libxtb files created in the xtb-python/build directory to xtb-python/xtb. I ran the example above and still see EEQ charge constrain error. I am now a bit doubtful whether I built the packages the right way or not, would appreciate any pointers if I am doing this wrong.