grimme-lab
py-xtb and xtb report different solvation energies
I'm getting different solvation energies for water by using the xtb commandline interface and python interface.
Input:
3
O 0.000000 -0.000000 -0.066467
H 0.000000 0.758588 0.527436
H 0.000000 -0.798588 0.527436
commandline invocation:
xtb --alpb water --json --parallel 1 water.xyz
python script:
import numpy as np
from scipy import constants
from xtb.interface import Calculator, Param
from xtb.libxtb import VERBOSITY_MINIMAL , VERBOSITY_FULL
from xtb.utils import Solvent
ANGSTROM_IN_BOHR = constants.physical_constants['Bohr radius'][0] * 1.0e10
def run_xtb_python(atomic_numbers, positions, net_charge, solvent=None):
"""
Run xTB with the provided python API
There is no error check
"""
calc = Calculator(
Param.GFN2xTB, atomic_numbers, positions / ANGSTROM_IN_BOHR, net_charge
)
calc.set_verbosity(VERBOSITY_FULL)
if solvent is not None:
calc.set_solvent(solvent)
res = calc.singlepoint()
energy = res.get_energy()
# serialize and store a list object
gradient = (res.get_gradient() / ANGSTROM_IN_BOHR).tolist()
charges = res.get_charges().tolist()
print("eigenvalues", res.get_orbital_eigenvalues())
print("total charge", sum(charges))
print("dipole", res.get_dipole())
return {
"energy": energy,
"charges": charges,
"gradient": gradient
}
if __name__ == "__main__":
PERIODIC_TABLE = {"H": 1, "O": 8}
atomic_numbers, positions = [], []
net_charge = 0
with open("water.xyz") as fin:
token = next(fin)
natoms = int(token)
_ = next(fin)
for i in range(natoms):
elmnt, x, y, z = next(fin).split()
atomic_numbers.append(PERIODIC_TABLE[elmnt])
positions.append([float(x), float(y), float(z)])
solvent = Solvent.h2o
output_data = run_xtb_python(
np.array(atomic_numbers),
np.array(positions),
net_charge,
solvent
)
The xtb executable gives the following summary:
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: SUMMARY ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: total energy -5.086650340935 Eh ::
:: total w/o Gsasa/hb -5.072450640868 Eh ::
:: gradient norm 0.034523773934 Eh/a0 ::
:: HOMO-LUMO gap 14.203449851402 eV ::
::.................................................::
:: SCC energy -5.115460341784 Eh ::
:: -> isotropic ES 0.055279536996 Eh ::
:: -> anisotropic ES -0.001553606638 Eh ::
:: -> anisotropic XC 0.000001590124 Eh ::
:: -> dispersion -0.000124241951 Eh ::
:: -> Gsolv -0.024975539514 Eh ::
:: -> Gborn -0.010775839447 Eh ::
:: -> Gsasa 0.004887985565 Eh ::
:: -> Ghb -0.019499345345 Eh ::
:: -> Gshift 0.000411659713 Eh ::
:: repulsion energy 0.028810000815 Eh ::
:: add. restraining 0.000000000000 Eh ::
:: total charge 0.000000000000 e ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
the python execution gives this summary
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: SUMMARY ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: total energy -5.079191678312 Eh ::
:: total w/o Gsasa/hb -5.071742364786 Eh ::
:: gradient norm 0.037702233171 Eh/a0 ::
:: HOMO-LUMO gap 13.870116488842 eV ::
::.................................................::
:: SCC energy -5.108001654563 Eh ::
:: -> isotropic ES 0.043341913075 Eh ::
:: -> anisotropic ES -0.000309776201 Eh ::
:: -> anisotropic XC -0.000217021538 Eh ::
:: -> dispersion -0.000131539071 Eh ::
:: -> Gsolv -0.011920230479 Eh ::
:: -> Gborn -0.004470916953 Eh ::
:: -> Gsasa 0.000215386523 Eh ::
:: -> Ghb -0.009522143176 Eh ::
:: -> Gshift 0.001857443127 Eh ::
:: repulsion energy 0.028809976217 Eh ::
:: add. restraining 0.000000000000 Eh ::
:: total charge -0.000000000000 e ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
Apologies if I missed an import or other error in the python, I'm extracting this from a longer workflow.
ps. this is on OSX 10.14.16 and I'm using the following xtb version
- xtb version 6.3.3 (71d3805) compiled by 'conda@b85dec0bf610' on 2021-01-07
The API still returns GBSA solvation free energies, which are available from the binary with --gbsa. I haven't got a chance to rework the API to support both GBSA and ALPB solvation free energies, but this is on my TODO list (see https://github.com/grimme-lab/xtb/issues/336).
Thanks for the comment Sebastian. Confirmed that xtb --gbsa water --json --parallel 1 water.xyz returns the same energy as the python api. I'll just use the CLI for now.