Hello, I’ve been working on merging itp files for my ligands using the ‘ligandHybrid’ script, but I’m encountering issues with specific ligands. While the tool works fine for some, others result in errors during the hybrid topology generation process. Here’s a brief overview of the error messages I’m receiving:
ligandHybridTop__log_> ERROR: fitting mol2 on mol1 failed
ligandHybridTop__log_> Initializing the main object for LigandHybridTopology
ligandHybridTop__log_> Starting hybrid topology generation
ligandHybridTop__log_> Making atom pairs…
…
Violation occurred on line 62 in file /home/conda/feedstock_root/build_artifacts/rdkit_1707396357915/work/Code/GraphMol/Conformer.cpp
Failed Expression: 19 < 7
…
self._make_atom_pairs( )
File “/usr/people/daniel/yehon/anaconda3/envs/rdkit-env/lib/python3.8/site-packages/pmx/ligand_alchemy.py”, line 2172, in _make_atom_pairs
a4 = self.m4.fetch_atoms(n2,how=‘byid’)[0]
IndexError: list index out of range
The process seems to fail when trying to fit molecule 2 on molecule 1, and eventually throws an IndexError during the creation of atom pairs. I think it might be related to the ‘rdkit’ package, which I installed via Conda following the standard procedure. Furthermore, in the ‘Failed Expression: 19 < 7’ the number 19 is not constant, it change with the increase of the ligand number of atoms.
I’ve attached the complete error log and the files relevant to the process (PDB files for the ligands, modified itp files, and the pairs.dat file) for reference. error.txt (17.1 KB)
It seems like rdkit isn’t finding an appropriate alignment between these two molecules, which could be because of the difference in the two molecules. NC1 contains just two heavy atoms whereas NC18 is a linear chain of 19 heavy atoms.
You could play with the -n1 and -n2 flags of “pmx ligandHybrid” and see.
Also, do you see a similar issue when working with two more similar molecules, say NC1 and NC2?
I tried using -n1 and -n2 with ndx files created by ‘gmx make_ndx’ without success; the problem remains the same. I believe this problem is directly related to the number of atoms. For example, in NC18, I get “Failed Expression: 19 < 7,” while for NC15, I get “Failed Expression: 14 < 7.” Indeed, it works for the transition of smaller molecules to NC1.
It seems to be an rdkit issue, and I don’t know if there is any workaround.
Even if you succeed in aligning two largely different molecules like NC1 and NC18, it would be hard to get a converged free energy estimate, because of sampling issues. NC18 can have a variety of conformations, and it would not be easy to sample those with nonequilibrium transitions from NC1.
Is it possible to divide the transformation into smaller steps, like NC1-NC3, NC3-NC5, …, NC15-NC17?
You have mixed the order of your molecules in the command line. The pairs.dat file has atoms of NC1 in the first column and atoms of NC18 in the second column. So you need to consistently provide -i1 NC1.pdb and -i2 NC18.pdb (the same with -itp options). The following command works fine: