ValueError: Molecule doesn't match the core, when running prepare_dual_topology.py

[zangsk]

Thanks for sharing your great work. I successfully ran tutorial01 without any issues, But when I try to read the pdbqt files, something goes wrong.

I used AutoDock Vina to perform molecular docking on the FXR ligands from the tutorial. All input pdbqt files were generated by meeko from mol files in tutorial01 , therefore contains RMARK SMILES. All other files were from tutorial01.

Here is my step2.ini.
[prepare_dual_topology]
input_ligands = lig_dock
structure = receptor_data/5q17_processed.pdb
extradirs = oplsaam.ff
pose_loader = vina
poses_reference_structure = receptor_data/5q17.pdbqt
poses_input = list.name
perturbations_dir = tutorial01
gmx_bin_local = /XXXX/bin/gmx
gmx_bin_run = /XXXX/bin/gmx_mpi
output_scripttype = slurm
output_resources = all_cpus:16; all_gpus:1; all_time: 24
output_collecttype = python

The error message is as follows
=============== mdp and run steps=================
Complex Water
min01.mdp min01.mdp
min02.mdp min02.mdp
min03.mdp nve.mdp
nve.mdp nvt.mdp
nvt.mdp npt.mdp
npt.mdp md.mdp
md.mdp -

================= Input ligands ==================
Name Molecule Topology
FXR_12 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Cl)c1[H] lig_dock/FXR_12.itp
FXR_74 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Br)c1[H] lig_dock/FXR_74.itp
FXR_76 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c([H])c1[H] lig_dock/FXR_76.itp
FXR_84 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(F)c1[H] lig_dock/FXR_84.itp
FXR_85 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C([H])([H])[H])c1[H] lig_dock/FXR_85.itp
FXR_88 [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C(F)(F)F)c1[H] lig_dock/FXR_88.itp

================= Perturbations ==================
State A State B
FXR_12 FXR_74
FXR_12 FXR_76
FXR_12 FXR_84
FXR_12 FXR_85
FXR_12 FXR_88

Reading receptor data from receptor_data/5q17.pdbqt as a pdb file
=================== Poses read ===================
Name File Details
FXR_12 FXR_12.pdbqt cluster 1
================== Align poses ===================
Molecule FXR_12 aligned
================ Working on pairs ================
Perturbation Pose Coordinates
============ Working on FXR_12-FXR_74 ============

=================== STACK INFO ===================
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/prepare_dual_topology.py", line 4480, in
savestate=progress_data, verbosity=arguments.verbose)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 657, in constrained_embed_dualmol
save_state=savestate, verbosity=verbosity, mcs_type=mcs_type)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 545, in constrained_embed_shapeselect
randomseed=randomseed, verbosity=verbosity, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 98, in constrained_embed_forcefield
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 313, in local_print
formatted_string = '\n{:=^50}\n{}{:=^50}\n'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] Failed to match molecule [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Br)c1[H] to core [#8]#6[#6]1~[#6][#6]#6[#6][#6]~1 in constrained_embed_forcefield. Cannot continue
Traceback (most recent call last):
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/prepare_dual_topology.py", line 4480, in
savestate=progress_data, verbosity=arguments.verbose)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 657, in constrained_embed_dualmol
save_state=savestate, verbosity=verbosity, mcs_type=mcs_type)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 545, in constrained_embed_shapeselect
randomseed=randomseed, verbosity=verbosity, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/public/home/public/software/PyAutoFEP/PyAutoFEP/merge_topologies.py", line 99, in constrained_embed_forcefield
raise ValueError("Molecule doesn't match the core")
ValueError: Molecule doesn't match the core

I also referred to Issues 14, but I have no idea. Could you please tell me how to solve this problem and how to use pdbqt correctly for input?
Many thanks

[caopeifang]

Hello, I encountered the same problem. Have you solved it? Below is my error message:

=================== STACK INFO ===================
File "/data/run01/scz1441/PyAutoFEP/prepare_dual_topology.py", line 4166, in
merged_data = merge_topologies.constrained_embed_dualmol(merged_data, poses_mol_data[state_a_name],
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 652, in constrained_embed_dualmol
temp_mol_b = constrained_embed_shapeselect(temp_mol_b, target, core_conf_id=core_conf_id, coord_map=coord_map,
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 544, in constrained_embed_shapeselect
constrained_embed_forcefield(molecule, core=core_mol, atom_map=atom_map, num_conformations=num_conformers,
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 95, in constrained_embed_forcefield
os_util.local_print('Failed to match molecule {} to core {} in constrained_embed_forcefield. Cannot '
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 313, in local_print
formatted_string = '\n{:=^50}\n{}{:=^50}\n'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] Failed to match molecule [H]OC([H])([H])C@@(C(=O)N([H])[H])N([H])C(=O)c1c(C([H])([H])[H])oc2c([H])c([H])c(Oc3c([H])c([H])c([H])c4nc(C([H])([H])[H])c([H])c([H])c([H])c34)c([H])c12 to core [#6]1([#6@])[#6](#6[#6]2~[#6]#6[#6][#6][#6]2[#8]~1 in constrained_embed_forcefield. Cannot continue
Traceback (most recent call last):
File "/data/run01/scz1441/PyAutoFEP/prepare_dual_topology.py", line 4166, in
merged_data = merge_topologies.constrained_embed_dualmol(merged_data, poses_mol_data[state_a_name],
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 652, in constrained_embed_dualmol
temp_mol_b = constrained_embed_shapeselect(temp_mol_b, target, core_conf_id=core_conf_id, coord_map=coord_map,
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 544, in constrained_embed_shapeselect
constrained_embed_forcefield(molecule, core=core_mol, atom_map=atom_map, num_conformations=num_conformers,
File "/data/run01/scz1441/PyAutoFEP/os_util.py", line 558, in wrap_trace
return f(*args, **kwargs)
File "/data/run01/scz1441/PyAutoFEP/merge_topologies.py", line 99, in constrained_embed_forcefield
raise ValueError("Molecule doesn't match the core")
ValueError: Molecule doesn't match the core

And below is the content of my step2.ini file:

This is the configuration section for prepare_dual_topology.py (which is the only section in this file, but the [ section ] is mandatory)

[prepare_dual_topology]

Read ligands and topologies from this folder

input_ligands = lig_data

Read the macromolecule structure from this file

structure = receptor_data/protein_processed.pdb

This is the force field directory, it will be copied to each perturbation dir and used to prepare the MD systems

extradirs = amber14sb_parmbsc1.ff

Options controlling the core-constrained superimposition

First select the use of it instead of reading all ligand poses

pose_loader = superimpose

Use this pose as the reference for the superimposition

poses_reference_pose_superimpose = lig_data/ND01203.mol2

Name of the output. This will be a self-extracting bash file

perturbations_dir = tutorial_1

Sets the path to GROMACS executable in the run node, uncomment and modify if needed.

gmx_bin_run = /data/run01/scz1441/new/gromacs-2022.5/install-gromacs/bin

Options controlling the output, see manual for more info. Uncomment as needed

FEP legs are to be submitted to a slurm scheduler

output_scripttype = slurm

Run these commands at the beginnig of the jog (useful to load modules, importing libs)

output_runbefore = module load python3; module load cuda

Fine-tune job resources

output_resources = all_cpus:6; all_gpus:1; all_time: 48

Use a python file instead of a binary during the collect step

output_collecttype = python