ASE support for net charge

[bieniekmateusz]

trafficstars

Describe the bug

Directly using XTB-python interface the charge= parameter can be passed, but employing the same approach with ASE/XTB, charge doesn't get propagated

Ie, this works for charge=0 and charge=1 and changes energy accordingly

calc = Calculator(Param.GFN2xTB, mol.get_atomic_numbers(), mol.get_scaled_positions(), charge=1)
res = calc.singlepoint()
print('energy', res.get_energy())

However, this doesn't remember the charge and always sets charge to 0:

mol2 = read('qca_seed10.xyz') 
mol2.set_calculator(XTB(method="GFN2-xTB", charge=1))
print('pe', mol2.get_potential_energy())

To Reproduce

unzip the file and run the script f.zip

1. xtb.version is 20.1
2. output:

direct calc energy 71.23435304058115
direct calc energy 71.54337595108338
energy before conversion -63.89323100198211
ase/xtb pe -1738.623373139004
energy before conversion -63.89323100198211
ase/xtb pe -1738.623373139004

Expected behaviour

charge affecting the final energy when ase/xtb is used

The issue seems to be in the XTB class : _charge = self.atoms.get_initial_charges().sum()

Thanks!

[awvwgk]

ASE is using an array of initial_charges which is part of the Atoms object rather than the Calculator object, therefore the XTB ASE calculator does not support setting a charge but determines it from the initial_charges property of the Atoms object it is working with.

[bieniekmateusz]

Thanks for having a look @awvwgk . I know what you mean, I am still worried about the inconsistency. Also, as I am using a .xyz file without the charges, this leaves me with no room to correct that with ase/xtb-python interface. Thanks!

[awvwgk]

My point here is similar to the argument in https://github.com/grimme-lab/xtb-python/issues/42 and https://github.com/grimme-lab/xtb-python/issues/47 for the multiplicity. Since ASE provides a canonical way to solve this issue, the XTB calculator will not implement something to bypass this mechanism.

Keep in mind that in ASE you can actually change the charge of your system after creating the calculator, the calculator must adjust to this change correctly, by adding this property to the calculator it becomes more difficult to keep the information in sync with the geometry used.

[bieniekmateusz]

I must have missed the canonical way that you're mentioning, would you mind pointing me towards? I naively assumed that the XTB calculator would a similar charge parameter to the example case.

Thanks

[awvwgk]

Have a look at the Atoms object constructor which supports:

charges: list of float
    Initial atomic charges.

And set_initial_charges as way to modify the system charge.

[aizvorski]

I have been working with other ase calculators (eg psi4) and can confirm that it is possible to have a net charge for the whole system without having the initial_charges array set. I don't know too much about how this works under the hood, but this is a minimal example to reproduce

atoms = ase.io.read('test.sdf')
calc = ase.calculators.psi4.Psi4(atoms = atoms,  charge = -1, multiplicity = 1)
atoms.calc = calc
print(atoms.get_initial_charges())  # this is all zeros
print(atoms.get_potential_energy()) # this works correctly and uses the charge passed in

It would be cool is the xtb calculator worked the same way

[awvwgk]

I asked in the ASE matrix channel and got the following answer:

Total charge and multiplicity are not normalized in any way. So if your code accepts charge=5 in the input file, then the calculator would presumably accept the same thing. (atoms.get_initial_charges() exists, but its purpose is more like a help for wavefunction/orbital initialization, mostly for the benefit of GPAW, rather than a way to actually control the total charge.)

I'm a bit unhappy about the lack of standardization here, since this means we have to implement our own caching logic for the total charge and multiplicity for the XTB ASE calculator.

[awvwgk]

For backwards compatibility we should keep the possibility to read from initial_charges but I'm open to accept a patch to implement an overwrite with charge and multiplicity parameters for the XTB ASE calculator (no hacking in Fortran or C involved, only Python needed).

[Shunyang2018]

Thank you all for the discussion! I also try to get this charge function recently.

1. I pass the ASE Atoms object to xtb calculator to bypass the API:

from xtb.interface import Molecule,Calculator, Param, Environment
calc = Calculator(Param.IPEAxTB, mol.arrays['numbers'],mol.arrays['positions'],charge=charge,uhf=spin)
calc.set_output('log.txt')
res =  calc.singlepoint().get_energy()

2. can we just read the charge from ase-xtb api rather than try to modify the initial charges array? instead of: https://github.com/grimme-lab/xtb-python/blob/0a6fd5c40b26750e18324beffe3e83e71a21d6de/xtb/ase/calculator.py#L190

we can use similar codes in the psi4 API:

# https://wiki.fysik.dtu.dk/ase/_modules/ase/calculators/psi4.html#Psi4

charge = self.parameters.get('charge')
mult = self.parameters.get('multiplicity')

3. we can use Atoms.info to save the spin and charge information here is an example output:

xyz.info = {'charge':1.0,'spin':0.0}

xyz.info
Out[37]: {'charge': 1.0, 'spin': 0.0}

now, you can use

_charge = self.atoms.info['charge']

Hope those can help! It would be cool to have this function!

[njzjz]

A simple hack to assign net charge for XTB:

import numpy as np
from xtb.ase.calculator import XTB as OldXTB


class XTB(OldXTB):
    """ASE calculator for XTB with net charge."""

    def __init__(self, *args, charge=0, **kwargs):
        self.charge = charge
        super().__init__(*args, **kwargs)

    def _create_api_calculator(self):
        initial_charges = np.zeros(len(self.atoms))
        initial_charges[0] = self.charge
        self.atoms.set_initial_charges(initial_charges)
        return super()._create_api_calculator()