py-xtb and xtb report different solvation energies

[leifjacobson]

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.

[leifjacobson]

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

[awvwgk]

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 grimme-lab/xtb#336).

[leifjacobson]

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.