Calculations on charged molecules

[vymetal]

Dear developers, I tried to calculate charge distribution of positively charged molecules using espaloma charge v 0.0.8 from pip distribution, but the resulting partial atomic charges always sum to 0.0. The following piece of code demonstrates the behavior

from espaloma_charge import charge
from rdkit import Chem
mol=Chem.MolFromSmiles("[NH4+]")
mol=Chem.AddHs(mol)
tc=Chem.GetFormalCharge(mol)
print(tc);
print([atom.GetFormalCharge() for atom in mol.GetAtoms()])
c=charge(mol,total_charge=tc)
print(c)
print(sum(c))`

The corresponding output is

1
[1, 0, 0, 0, 0]
/opt/conda/lib/python3.10/site-packages/dgl/heterograph.py:92: DGLWarning: Recommend creating graphs by `dgl.graph(data)` instead of `dgl.DGLGraph(data)`.
  dgl_warning(
[-0.93732476  0.2343312   0.2343312   0.2343312   0.23433116]
1.4901161193847656e-08

Is it a bug or I just missinterpreted the output?

Thank you

Jiri

[mikemhenry]

@ijpulidos @yuanqing-wang Thoughts?

[hannahbruce]

This doesn't seem to be an issue when using the openff-toolkit to charge, if that is a work-around or indicates where the issue might be coming from

from espaloma_charge import charge
from rdkit import Chem
from openff.toolkit.topology import Molecule
from rdkit import Chem
from rdkit.Chem import AllChem
from espaloma_charge.openff_wrapper import EspalomaChargeToolkitWrapper
toolkit_registry = EspalomaChargeToolkitWrapper()

smiles = "[NH4+]"

molecule = Molecule.from_smiles(smiles)
molecule.assign_partial_charges('espaloma-am1bcc', toolkit_registry=toolkit_registry)
espaloma_charges = molecule.partial_charges
print(espaloma_charges)
print(sum(molecule.partial_charges))

mol=Chem.MolFromSmiles(smiles)
mol=Chem.AddHs(mol)
tc=Chem.GetFormalCharge(mol)
print(tc);
print([atom.GetFormalCharge() for atom in mol.GetAtoms()])
c=charge(mol,total_charge=tc)
print(c)
print(sum(c))

[yuanqing-wang]

This is an issue that we've fixed a while ago. Could you try pulling from the source again?