Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Expanding for multistates #12

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open

Expanding for multistates #12

wants to merge 16 commits into from

Conversation

RiesBen
Copy link
Contributor

@RiesBen RiesBen commented Jul 19, 2023

This pull request is all about Multistate Free Energy Methods.

Features:

  • implementation of multistate mapping
  • Visualization of multistate mapping

Motivation:

The goal is to enable methods like (RE-)EDS, LadyBugs or other multistate methods to get quick access to hybrid topology approach atom-mappings.

Theory:

The algorithm concept:

  • generate all NxN mappings
  • filter the mappings
  • merge down the mappings to a global mapping, that contains linked atom mappings for all states.
    • get the largest connected set of the global mapping.

Usage

Here example result, all ligands have a mapped core-region, such that all 14 ligands can be used in one hybrid topology simulation.

Code:

# Get input Data
from openfe_benchmarks import benzenes
components = benzenes.get_system().ligand_components

# Exclude cycle breakers! as not feasible for Hybrid topology approaches.
not_lig = ["lig_4", "lig_7",  "lig_10"] #"lig_3", "lig_2",
components = [c for c in components if(c.name not in not_lig)][:15]

#for visualization
#Chem.Draw.MolsToGridImage([c.to_rdkit() for c in components])

#mapp:
from kartograf import KartografAtomMapper

atomMapper = KartografAtomMapper()
multi_state_mapping = atomMapper.suggest_multistate_mapping(components, greedy=True, map_hydrogens=False)

multi_state_mapping

Result Visualization:

  • 3D- Representation
from kartograf.utils.multistate_visualization_3D import visualize_multistate_mappings_3D
visualize_multistate_mappings_3D(components, multi_state_mapping)

image

  • 2D- Representation
from kartograf.utils.multistate_visualization import visualize_multistate_mappings_2D
visualize_multistate_mappings_2D(components, multi_state_mapping)

image

returned value:
[{'lig_1': 4,
'lig_10': 6,
'lig_11': 3,
'lig_12': 3,
'lig_13': 3,
'lig_14': 4,
'lig_15': 4,
'lig_16': 4,
'lig_2': 10,
'lig_3': 10,
'lig_5': 6,
'lig_6': 4,
'lig_8': 8,
'lig_9': 8},
{'lig_1': 6,
'lig_10': 3,
'lig_11': 6,
'lig_12': 6,
'lig_13': 5,
'lig_14': 6,
'lig_15': 6,
'lig_16': 6,
'lig_2': 4,
'lig_3': 4,
'lig_5': 3,
'lig_6': 6,
'lig_8': 10,
'lig_9': 5},
{'lig_1': 8,
'lig_10': 15,
'lig_11': 11,
...
}]

@RiesBen RiesBen changed the title {WIP] Expanding for multistates [WIP] Expanding for multistates Jul 19, 2023
@codecov
Copy link

codecov bot commented Jul 19, 2023

Codecov Report

Attention: Patch coverage is 3.77358% with 153 lines in your changes missing coverage. Please review.

Project coverage is 77.06%. Comparing base (8ef8f88) to head (2d35bf3).
Report is 11 commits behind head on main.

Files with missing lines Patch % Lines
src/kartograf/atom_mapper.py 3.77% 153 Missing ⚠️
Additional details and impacted files
@@             Coverage Diff             @@
##             main      #12       +/-   ##
===========================================
- Coverage   96.35%   77.06%   -19.30%     
===========================================
  Files          13       13               
  Lines         604      763      +159     
===========================================
+ Hits          582      588        +6     
- Misses         22      175      +153     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@pep8speaks
Copy link

pep8speaks commented Oct 19, 2023

Hello @RiesBen! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 887:9: E128 continuation line under-indented for visual indent
Line 888:9: E128 continuation line under-indented for visual indent
Line 889:9: E128 continuation line under-indented for visual indent
Line 890:9: E128 continuation line under-indented for visual indent
Line 891:9: E128 continuation line under-indented for visual indent
Line 892:9: E128 continuation line under-indented for visual indent
Line 893:9: E125 continuation line with same indent as next logical line
Line 893:9: E128 continuation line under-indented for visual indent
Line 903:9: E303 too many blank lines (2)
Line 931:80: E501 line too long (98 > 79 characters)
Line 934:80: E501 line too long (85 > 79 characters)
Line 955:5: E303 too many blank lines (2)
Line 955:80: E501 line too long (85 > 79 characters)
Line 958:42: E251 unexpected spaces around keyword / parameter equals
Line 958:44: E251 unexpected spaces around keyword / parameter equals
Line 961:9: E265 block comment should start with '# '
Line 963:18: E225 missing whitespace around operator
Line 963:18: E712 comparison to True should be 'if cond is True:' or 'if cond:'
Line 977:80: E501 line too long (128 > 79 characters)
Line 984:80: E501 line too long (81 > 79 characters)
Line 985:80: E501 line too long (93 > 79 characters)
Line 986:80: E501 line too long (87 > 79 characters)
Line 987:80: E501 line too long (92 > 79 characters)
Line 989:13: E265 block comment should start with '# '
Line 990:80: E501 line too long (95 > 79 characters)
Line 991:80: E501 line too long (103 > 79 characters)
Line 994:9: E265 block comment should start with '# '
Line 1003:25: E117 over-indented
Line 1005:46: E231 missing whitespace after ':'
Line 1021:9: E303 too many blank lines (2)
Line 1026:80: E501 line too long (96 > 79 characters)
Line 1046:80: E501 line too long (80 > 79 characters)
Line 1058:9: E303 too many blank lines (2)
Line 1061:5: E303 too many blank lines (2)
Line 1063:9: E731 do not assign a lambda expression, use a def
Line 1070:61: E261 at least two spaces before inline comment
Line 1070:62: E262 inline comment should start with '# '
Line 1072:40: E261 at least two spaces before inline comment
Line 1072:41: E262 inline comment should start with '# '
Line 1074:63: E261 at least two spaces before inline comment
Line 1074:64: E262 inline comment should start with '# '
Line 1078:80: E501 line too long (85 > 79 characters)
Line 1079:80: E501 line too long (84 > 79 characters)
Line 1095:80: E501 line too long (83 > 79 characters)
Line 1096:80: E501 line too long (127 > 79 characters)
Line 1100:80: E501 line too long (95 > 79 characters)
Line 1101:80: E501 line too long (81 > 79 characters)
Line 1104:80: E501 line too long (86 > 79 characters)
Line 1105:80: E501 line too long (91 > 79 characters)
Line 1109:80: E501 line too long (87 > 79 characters)
Line 1113:80: E501 line too long (94 > 79 characters)
Line 1116:80: E501 line too long (83 > 79 characters)
Line 1121:80: E501 line too long (96 > 79 characters)
Line 1122:80: E501 line too long (82 > 79 characters)
Line 1126:80: E501 line too long (107 > 79 characters)
Line 1131:80: E501 line too long (113 > 79 characters)
Line 1132:9: E125 continuation line with same indent as next logical line
Line 1139:80: E501 line too long (85 > 79 characters)
Line 1147:80: E501 line too long (98 > 79 characters)
Line 1157:80: E501 line too long (108 > 79 characters)

Line 60:1: W391 blank line at end of file

Line 27:1: E303 too many blank lines (3)
Line 27:26: E231 missing whitespace after ':'
Line 28:5: E125 continuation line with same indent as next logical line
Line 54:1: E303 too many blank lines (3)
Line 56:22: E128 continuation line under-indented for visual indent
Line 56:44: E252 missing whitespace around parameter equals
Line 57:22: E128 continuation line under-indented for visual indent
Line 58:22: E128 continuation line under-indented for visual indent
Line 59:22: E124 closing bracket does not match visual indentation
Line 61:80: E501 line too long (86 > 79 characters)
Line 73:80: E501 line too long (98 > 79 characters)
Line 99:1: E302 expected 2 blank lines, found 1

Comment last updated at 2024-08-28 23:24:13 UTC

@RiesBen RiesBen changed the title [WIP] Expanding for multistates Expanding for multistates Oct 24, 2023
@eruijsena
Copy link

Hey @RiesBen !

So I tried using this to generate mappings for the Schindler et al. FEP benchmark, and it generally works very well. There are just some cases were the behavior isn't exactly what I would expect (just from an alchemical pathway perspective, I didn't look at your code so I don't know if these are bugs or features).

Here's the code to generate these problematic mappings:

from rdkit import Chem
from kartograf import KartografAtomMapper
from gufe.components.smallmoleculecomponent import SmallMoleculeComponent
import requests

# get ligands
url = "https://raw.githubusercontent.com/MCompChem/fep-benchmark/master/cmet/ligands.sdf"

response = requests.get(url)
with open("ligands_cmet.sdf", 'wb') as file:
    file.write(response.content)

mols = [x for x in Chem.SDMolSupplier('ligands_cmet.sdf', removeHs=False)]
components = [SmallMoleculeComponent.from_rdkit(mol) for mol in mols]

atomMapper = KartografAtomMapper()

# Get the second subset of 10 ligands from the original set
subset_components = components[10:20]
# Mapping
multi_state_mapping = atomMapper.suggest_multistate_mapping(subset_components, alternative=False, map_hydrogens=False)

This code yields the following mapping:
image

The alchemical transformations suggested by this mapping involve a weird ring closure/opening, which I imagine should not be allowed? I may be wrong though.

If I use the same code on another ligand set from the benchmark, url = "https://raw.githubusercontent.com/MCompChem/fep-benchmark/master/hif2a/ligands.sdf", I get

image

This one is clearly suspicious, as the nitrile and the sulfoxide oxygen definitely don't have the same geometry.

Wonder what you think about theses cases.

@RiesBen
Copy link
Contributor Author

RiesBen commented Nov 13, 2023

Hej @eruijsena thanks for your cool feedback! :)
I will look into the strange cases you got, I assume it might due to how the states are finally merged... let me see :)

@RiesBen RiesBen mentioned this pull request Nov 20, 2023
@RiesBen RiesBen linked an issue Nov 20, 2023 that may be closed by this pull request
@RiesBen RiesBen self-assigned this Dec 12, 2023
@RiesBen RiesBen added the enhancement New feature or request label Dec 12, 2023
@RiesBen
Copy link
Contributor Author

RiesBen commented Jan 29, 2024

Hi @eruijsena,
I finally had some time to improve our filter strategy and implemented a greedy mapping strategy, making everything faster! (I adapted the code example above to use the greedy approach and I added new visualization possibilities)

the new mappings give me a better impression, what do you think?

HERE we have a subset of Cmet:

2D - representation:

from kartograf.utils.multistate_visualization import visualize_multistate_mappings_2D
visualize_multistate_mappings_2D(components, multi_state_mapping)

image

3D - representation:

from kartograf.utils.multistate_visualization_3D import visualize_multistate_mappings_3D
visualize_multistate_mappings_3D(components, multi_state_mapping)

image

@RiesBen
Copy link
Contributor Author

RiesBen commented Jan 29, 2024

and here a subset of hif2a
image

I see again the nitrile-Sulfoxide oxygen mapping. I think for this we will need a new filter (see #35) to not map this, as the distances of our alignment make it possible to map those atoms (see the 3D):

image

@eruijsena
Copy link

Hey Ben.

Yes these do look better!

Would it be possible to have an explanation of what the greedy version does? I'm not sure I understand how it assigns a mapping.

For instance, if you turn the flag map_hydrogens to True in your code example with the benzenes set, I get this rather funky mapping:
Screenshot from 2024-03-12 17-00-28

Is this intended?

@RiesBen
Copy link
Contributor Author

RiesBen commented Apr 17, 2024

@eruijsena
aehm, no that is not a feature ;)
I will try to write an explanation on what the greedy algorithm does. (pls. keep in mind, this version of kartograf is still experimental, but I try to get it out of this phase :) )

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

MultiState Enablement
3 participants