error when using `psi4/resp`: `Fatal Error: RHF: RHF reference is only for singlets.`

[linminhtoo]

hello esther,

(I found a fix, please see my comment, but I don't know the root cause)
i've been trying out the espsim library on some molecules. the demo notebooks are very helpful and the analyses in the paper are also nicely done. great work!

i was trying to compare the ESP similarity of ATP against another known inhibitor of CDK2.
Specifically, I used the ATP co-crystal ligand from PDB ID 1B38,
(https://models.rcsb.org/v1/1b38/ligand?auth_seq_id=381&label_asym_id=C&encoding=sdf&filename=1b38_C_ATP.sdf)
and the Dinaciclib co-crystal ligand from PDB ID 5L2W
(https://models.rcsb.org/v1/5l2w/ligand?auth_seq_id=900&label_asym_id=C&encoding=sdf&filename=5l2w_C_1QK.sdf)

I am trying EmbedAlignScore() on these 2 mols, and the calculation works when I use gasteiger, mmff, ml, but not resp.
psi4 complains about this error:

RuntimeError: 
Fatal Error: RHF: RHF reference is only for singlets.
Error occurred in file: /build/source/psi4/src/psi4/libscf_solver/rhf.cc on line: 92
The most recent 5 function calls were:
psi::PsiException::PsiException(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, char const*, int)

I've also attached the detailed run log for your reference: https://gist.github.com/linminhtoo/84182da4bf727361b23905f34a429a5d

Here was how I run the code (I did rename some variables for my convenience but I didn't change any code logic)

# load co-crystal ATP
mol_atp = Chem.SDMolSupplier(str(RELATIVE / "data/cdk2_knownhits/1b38_C_ATP.sdf"), removeHs=False)[0]
mol_atp = Chem.AddHs(mol_atp, addCoords=True)

# load co-crystal dinaciclib
mol_dina = Chem.SDMolSupplier(str(RELATIVE / "data/cdk2_knownhits/5l2w_C_1QK.sdf"), removeHs=False)[0]
mol_dina = Chem.AddHs(mol_dina, addCoords=True)

# run ESPSim
shape_sims, esp_sims = EmbedAlignScore(
    probe_mol=deepcopy(mol_dina),
    ref_mols=deepcopy(mol_atp),
    probe_num_confs=10,
    ref_num_confs=10,
    partial_charge_mode="resp",
    renormalize=True,  # to [0, 1]
    getBestESP=True,  # more accurate but slower
    randomseed=2342,
)

I tried swapping probe_mol to mol_atp and ref_mols to mol_dina and psi4 did calculate the charges for one molecule before crashing again with the same error, which means it is not happy with Dinaciclib for some reason, but is fine with ATP. (assuming it calculates probe_mol first.

Do you happen to have any idea what's the issue?

Best,
Min Htoo

[linminhtoo]

Okay, so as per my suspicion, I loaded Dinaciclib directly from SMILES, and it runs fine now... so I guess it wasn't happy with something in the sdf file provided by RCSB, but I have no clue what's the issue haha

# get from SMILES
mol_dina = Chem.MolFromSmiles("CCC1=C2N=C(C=C(N2N=C1)NCC3=C[N+](=CC=C3)[O-])N4CCCC[C@H]4CCO")

This seems to be OK since ESPSim embeds & enumerates a number of conformers, but in certain cases we might want to keep the "reference molecule"'s conformation fixed (if it's a co-crystallised conformation).
Does it make sense to want to keep the reference conformer fixed? I was thinking this would also cut down the runtime significantly.
edit: I just saw the re-scoring benchmark notebooks where you indeed use the co-crystal conformer as the ref_mol. However, I had some difficulty interpreting the re-scoring metrics at the end of your notebook. (https://github.com/hesther/espsim/blob/master/benchmarks/benchmark_3_d4-rescore.ipynb)

What would you recommend for re-scoring AutoDock Vina poses against a co-crystal ligand? (preferably as fast as possible) In my case, the docked ligands can be quite different structurally from the co-crystal (may not even share a common substructure), so I wonder if sampling random poses or using the docked pose is better.

[hesther]

Hi @linminhtoo , I will look into this and get back to you as soon as I can. Unfortunately, I am currently on maternal leave, so it might take a while, sorry! My first thought would be that the Dinaciclib SDF file has a problem, causing it to be read in wrongly in RDKit. I.e. if you just create that molecule and check all the valences, I would expect at least one of them to be off. RDKit sometimes has trouble reading SDF files from other sources.

Using a docked pose is much faster, since embedding the molecules takes a very long time with RDKit. So if you are looking for speed and already have poses that somewhat make sense, I would recommend that. Ideally, you would have a benchmark system that you could try both methods on, and decide what works best for your case.

[hesther]

Also, if the SDF is read in wrong, it might be better to not trust any of the charges (not only resp)