luancarvalhomartins / PyAutoFEP

PyAutoFEP: an automated FEP workflow for GROMACS integrating enhanced sampling methods
161 stars 74 forks source link

The molecules didn't align together as expected #64

Closed 1200Caixia closed 1 year ago

1200Caixia commented 2 years ago

I don't know why the molecules in the red box don't align together, they only differ on the NH and O atom image

luancarvalhomartins commented 2 years ago

This is not because of the alignment, but rather because of the current dual-topology model implemented in PyAutoFEP. From #61:

Currently, PyAutoFEP cannot (a) perturb bonds and (b) concomitantly perturb charges and VdW parameters of an atom; so that it is not possible to perturb a C to a N in ring right now. That's why the MCS code won't allow for such core (it would give spurious, unphysical results). I am finishing a code for doing charge-changing perturbations and that code will allow for such perturbations to be done (as a side effect of the implementation). Should you be able to wait for a few weeks, it will make into the tree.

So with the current version of the code, it is not possible to build a dual-topology containing the benzofuran/indole carbons in the common core region. Therefore, PyAutoFEP will add the benzofuran and indole to the parts A and B. In general, using this dual-topology implementation should work similarly to perturbing only the O to N.

1200Caixia commented 2 years ago

Thanks for your detailed reply. PyAutoFEP is a convenient program. Looking forward the new code.

1200Caixia commented 2 years ago

Is FXR_78 for a similar reason? On FXR_12→FXR_78, generate_perturbation_map.py and prepare_dual_topology.py run for a long time without results image

luancarvalhomartins commented 1 year ago

No, that's the find_mcs code falling back to rdkit.Chem.rdFMCS.FindMCS with explicit hydrogens because of the symmetry of the 2,6-DiClPh ring. find_mcs does a lot of stuff to avoid calling rdkit.Chem.rdFMCS.FindMCS with explicit Hs because it takes forever to run. You hitted a bad case where that fallback was happening. I am committing a patch fixing it right now.

luancarvalhomartins commented 1 year ago

Note: fixes only the FXR_78 issue. Perturbing atoms in rings is yet to come.

1200Caixia commented 1 year ago

Yes, it works

1200Caixia commented 1 year ago

image On FXR_12→FXR_78,FXR78 changed to the same molecule as FXR12 in the fullsystem_step6.pdb.

luancarvalhomartins commented 1 year ago

It most likely isn't. This is likely an artifact of the representation you are using. There should be an hydrogen overlapped by the C-Cl bond in the left side of the image. Try using a ball and stick representation, eg:

fxr12fxr78

1200Caixia commented 1 year ago

when "pose_loader = superimpose", FXR_12→FXR_78 have the right structure. however, when "pose_loader = generic", the wrong structure will appear. image

luancarvalhomartins commented 1 year ago

Sorry, I cannot reproduce your issue. Using pose_loader = generic worked for me (apart, sure, for the poses being loaded from the input structure instead of superimposed to a ref.).

Would you mind sharing your complete input file and command line args? Also, would you mind retrying without a progress file (eg, remove or rename your progress.pkl)?

1200Caixia commented 1 year ago

workdir.zip These are the input files. "prepare_dual_topology.py --config_file=step2-superimpose.ini" will get the right structure, However, "prepare_dual_topology.py --config_file=step2-generic.ini" will get the wrong structure. progress.pkl has been removed every time.