protocaller / ProtoCaller

Full automation of relative protein-ligand binding free energy calculations in GROMACS
http://protocaller.readthedocs.io
GNU General Public License v3.0
43 stars 15 forks source link

How to restore ligand from parametrised file? #32

Closed kexul closed 3 years ago

kexul commented 3 years ago

Hi, this is the code I used to parametrise my ligand:

smiles = 'ClC1=CC(C2=CC=NC3=C2SC(C(CCl)(CO)O)=C3)=C4C(C[C@H](C(N5CCNCC5)=O)O4)=C1'
lig = Ligand(smiles, protonated=False, minimise=True, workdir='test_case')
lig.protonate()
lig.parametrise()

When I want to restore the ligand with already parametrised files,

smiles = 'ClC1=CC(C2=CC=NC3=C2SC(C(CCl)(CO)O)=C3)=C4C(C[C@H](C(N5CCNCC5)=O)O4)=C1'
lig = Ligand(smiles, protonated=False, minimise=True, workdir='test_case', parametrised_files=['ligand1.prmtop', 'ligand1.inpcrd'])
lig.protonate()
lig.parametrise()

I got an error:

[16:06:09] Explicit valence for atom # 30 N, 4, is greater than permitted
Traceback (most recent call last):
  File "/root/repo/prepare/teth_validate.py", line 184, in <module>
    lig = Ligand(smiles, protonated=False, minimise=True, workdir='test_case', parametrised_files=['ligand1.prmtop', 'ligand1.inpcrd'])
  File "/data/FEP/prepare/ProtoCaller/Ensemble/Ligand.py", line 56, in __init__
    template=input)
  File "/data/FEP/prepare/ProtoCaller/Wrappers/rdkitwrapper.py", line 194, in openAsRdkit
    mol = AssignBondOrdersFromTemplate(template, mol)
  File "/data/FEP/prepare/ProtoCaller/Wrappers/rdkitwrapper.py", line 296, in AssignBondOrdersFromTemplate
    _Chem.SanitizeMol(mol)
rdkit.Chem.rdchem.AtomValenceException: Explicit valence for atom # 30 N, 4, is greater than permitted

I tried set protonated=True and minimise=False, the error still exists.

kexul commented 3 years ago

I tried use the generated ligand1_protonated.pdb file and ligand1_protonated.sdf, which produce no error. But, when I used the ligand1_protonated_antechamber.mol2 as input, another error occurred:

****
Post-condition Violation
Element 'o' not found
Violation occurred on line 91 in file /home/conda/feedstock_root/build_artifacts/rdkit_1578668820522/work/Code/GraphMol/PeriodicTable.h
Failed Expression: anum > -1
****

/data/FEP/prepare/ProtoCaller/Wrappers/rdkitwrapper.py:189: UserWarning: No template containing bond information was supplied. Resulting geometry might be wrong
  _warnings.warn("No template containing bond information was supplied. "
Traceback (most recent call last):
  File "/root/repo/prepare/teth_validate.py", line 184, in <module>
    lig = Ligand('ligand1_protonated_antechamber.mol2', protonated=True, minimise=False, workdir='test_case', parametrised_files=['ligand1.prmtop', 'ligand1.inpcrd'])
  File "/data/FEP/prepare/ProtoCaller/Ensemble/Ligand.py", line 56, in __init__
    template=input)
  File "/data/FEP/prepare/ProtoCaller/Wrappers/rdkitwrapper.py", line 194, in openAsRdkit
    mol = AssignBondOrdersFromTemplate(template, mol)
  File "/data/FEP/prepare/ProtoCaller/Wrappers/rdkitwrapper.py", line 282, in AssignBondOrdersFromTemplate
    atom1 = match[b.GetBeginAtomIdx()]
KeyError: 7
msuruzhon commented 3 years ago

This must be an issue with how RDKit sanitises SMILES strings, I think that's one of the reasons why SDF files should be preferred. An AMBER mol2 file will never work with RDKit because it's not actually a valid mol2 file (it's more like a parameter file, much like GROMACS top), it's just called like that for historical reasons.