Parameterization issue for protein receptor + Zn ion with espaloma charges

[itaneja2]

When using Meeko to create a Polymer consisting of a protein receptor and Zn ion, I receive the following error when trying to parameterize the Polymer with espaloma charges:

mk_config = { 
    "merge_these_atom_types": [], 
    "load_atom_params": ["openff"],
    "charge_model": "espaloma",
    "dihedral_model": "openff",
    "rigid_macrocycles": True, 
}   
mk_prep = MoleculePreparation.from_config(mk_config)
polymer.parameterize(mk_prep)
File "/gpfs/home/itaneja/Meeko/meeko/polymer.py", line 1131, in parameterize
    self.monomers[residue_id].parameterize(mk_prep, residue_id)
File "/gpfs/home/itaneja/Meeko/meeko/polymer.py", line 2080, in parameterize
    molsetups = mk_prep(self.padded_mol)
File "/gpfs/home/itaneja/Meeko/meeko/preparation.py", line 443, in __call__
    return self.prepare(*args)
File "/gpfs/home/itaneja/Meeko/meeko/preparation.py", line 510, in prepare
    molgraph = self.espaloma_model.get_espaloma_graph(setup)
File "/gpfs/home/itaneja/Meeko/meeko/espalomatyper.py", line 49, in get_espaloma_graph
    molgraph = self.EspalomaGraph(openffmol)
File "/gpfs/home/itaneja/micromamba/envs/autodock/lib/python3.9/site-packages/espaloma/graphs/graph.py", line 63, in __init__
    heterograph = self.get_heterograph_from_graph_and_mol(
File "/gpfs/home/itaneja/micromamba/envs/autodock/lib/python3.9/site-packages/espaloma/graphs/graph.py", line 137, in get_heterograph_from_graph_and_mol
    heterograph = esp.graphs.utils.read_heterogeneous_graph.from_homogeneous_and_mol(
File "/gpfs/home/itaneja/micromamba/envs/autodock/lib/python3.9/site-packages/espaloma/graphs/utils/read_heterogeneous_graph.py", line 150, in from_homogeneous_and_mol
    for x in idxs["n%s" % big_idx][
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

input file: 5msb_A_receptor_water.json

[rwxayheee]

Hi @itaneja2
Could you please share the input file and command/script you used? It looks like an error raised from espaloma, but happy to help and see if we have workaround with meeko.

[rwxayheee]

Hi @itaneja2

I can reproduce the same error on an oxide ion ([O2-]). See details at the end of comment. Therefore, the origin of this issue might not be limited to Zn being an unsupported element (see choderalab/espaloma#206) but also any lone atom.

I will look at Meeko's codes, the get_espaloma_graph part, to see if we have a workaround (like direct assignment of formal charge), and if we can implement guards and better error handling for this situation.

Thanks again for reporting.

espaloma version

'0.3.2+3.g907d183'

input

from meeko import MoleculePreparation
from meeko import Monomer
from rdkit import Chem
from rdkit.Chem import AllChem

mk_config = { 
    "merge_these_atom_types": [], 
    "load_atom_params": ["openff"],
    "charge_model": "espaloma",
    "dihedral_model": "openff",
    "rigid_macrocycles": True, 
} 
mk_prep = MoleculePreparation.from_config(mk_config)

oa_smi = "[O-2]"
oa_mol = Chem.MolFromSmiles(oa_smi)
AllChem.EmbedMolecule(oa_mol)

oa_monomer = Monomer(
    raw_input_mol=oa_mol, 
    rdkit_mol=oa_mol, 
    mapidx_to_raw = {0:0}, 
)
oa_monomer.padded_mol = oa_mol
oa_monomer.parameterize(mk_prep, ["A:1"])

output

Traceback (most recent call last):
  File "/Users/amyhe/Downloads/meeko_esp.py", line 26, in <module>
    oa_monomer.parameterize(mk_prep, ["A:1"])
  File "/Users/amyhe/Desktop/0_forks/Meeko/meeko/polymer.py", line 2456, in parameterize
    molsetups = mk_prep(self.padded_mol)
  File "/Users/amyhe/Desktop/0_forks/Meeko/meeko/preparation.py", line 474, in __call__
    return self.prepare(*args, **kwargs)
  File "/Users/amyhe/Desktop/0_forks/Meeko/meeko/preparation.py", line 543, in prepare
    molgraph = self.espaloma_model.get_espaloma_graph(setup)
  File "/Users/amyhe/Desktop/0_forks/Meeko/meeko/espalomatyper.py", line 49, in get_espaloma_graph
    molgraph = self.EspalomaGraph(openffmol)
  File "/Users/amyhe/micromamba/envs/has_espaloma/lib/python3.10/site-packages/espaloma-0.3.2+3.g907d183-py3.10.egg/espaloma/graphs/graph.py", line 63, in __init__
    heterograph = self.get_heterograph_from_graph_and_mol(
  File "/Users/amyhe/micromamba/envs/has_espaloma/lib/python3.10/site-packages/espaloma-0.3.2+3.g907d183-py3.10.egg/espaloma/graphs/graph.py", line 137, in get_heterograph_from_graph_and_mol
    heterograph = esp.graphs.utils.read_heterogeneous_graph.from_homogeneous_and_mol(
  File "/Users/amyhe/micromamba/envs/has_espaloma/lib/python3.10/site-packages/espaloma-0.3.2+3.g907d183-py3.10.egg/espaloma/graphs/utils/read_heterogeneous_graph.py", line 150, in from_homogeneous_and_mol
    for x in idxs["n%s" % big_idx][
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

[itaneja2]

Simple workaround if ion charge parameters can be treated independent of ligand:

mk_config = { 
    "merge_these_atom_types": [], 
    "load_atom_params": ["openff"],
    "charge_model": "espaloma",
    "dihedral_model": "openff",
    "rigid_macrocycles": True,
}   

#espaloma doesnt like ions 
mk_config_ion = { 
    "merge_these_atom_types": [], 
    "load_atom_params": ["openff"],
    "charge_model": "gasteiger",
    "dihedral_model": "openff",
    "rigid_macrocycles": True,  
}   

mk_prep = MoleculePreparation.from_config(mk_config)
mk_prep_ion = MoleculePreparation.from_config(mk_config_ion)

for residue_id, residue in polymer.get_valid_monomers().items():
    if residue.input_resname in ['CL', 'MG', 'CA', 'MN', 'FE', 'ZN']:
        residue.parameterize(mk_prep_ion, residue_id)
    else:
        residue.parameterize(mk_prep, residue_id)

[rwxayheee]

According to choderalab/espaloma#191
It's a current limitation within espaloma that molecules with less than 3 atoms can't get a heterograph without the same error.

For single-atom molecules, we can pass the formal charge directly.

For possible diatomic molecules in biological systems (molecular oxygen, CO, ...), we still need a bypass to get charges and bond parameters. Not very sure how to do that at the moment..

For unsupported elements, it seems like espaloma usually raises UnassignedBondError with indications of element (atom) types. The error message is informative, and we can let it be, or if we want we can flag earlier and abort the parameterize attempt for unsupported elements.