diff --git a/.github/workflows/conda.yml b/.github/workflows/conda.yml index 05b37430..1ac05a84 100644 --- a/.github/workflows/conda.yml +++ b/.github/workflows/conda.yml @@ -10,7 +10,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macOS-latest] - python-version: [3.7, 3.8] + python-version: [3.8, 3.9] name: Python ${{ matrix.python-version }} at ${{ matrix.os }} steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/master.yml b/.github/workflows/master.yml index 6ee5fa66..85af5bd9 100644 --- a/.github/workflows/master.yml +++ b/.github/workflows/master.yml @@ -17,7 +17,7 @@ jobs: uses: s-weigand/setup-conda@v1 with: update-conda: true - python-version: 3.7 + python-version: 3.9 conda-channels: anaconda, omnia, conda-forge, martimunicoy - name: Install doc dependencies diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index cdcbd19a..6a9d26d2 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -16,7 +16,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macOS-latest] - python-version: [3.7, 3.8] + python-version: [3.8, 3.9] steps: - uses: actions/checkout@v2 diff --git a/README.md b/README.md index 253b95a1..2394120c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -| **About** | [![License](https://img.shields.io/badge/License-MIT-blue.svg)](LICENSE) [![Python](https://img.shields.io/badge/python-3.7%2C%203.8-blue.svg)](https://martimunicoy.github.io/peleffy) [![Release](https://img.shields.io/github/release/martimunicoy/peleffy.svg?include_prereleases)](https://github.com/martimunicoy/peleffy/releases/) | +| **About** | [![License](https://img.shields.io/badge/License-MIT-blue.svg)](LICENSE) [![Python](https://img.shields.io/badge/python-3.8%2C%203.9-blue.svg)](https://martimunicoy.github.io/peleffy) [![Release](https://img.shields.io/github/release/martimunicoy/peleffy.svg?include_prereleases)](https://github.com/martimunicoy/peleffy/releases/) | | :------ | :------- | | **Status** | [![Language grade: Python](https://img.shields.io/lgtm/grade/python/g/martimunicoy/peleffy.svg?logo=lgtm&logoWidth=18)](https://lgtm.com/projects/g/martimunicoy/peleffy/context:python) [![Test](https://github.com/martimunicoy/peleffy/workflows/Test/badge.svg)](https://github.com/martimunicoy/peleffy/actions?query=workflow%3ATest) [![codecov](https://codecov.io/gh/martimunicoy/peleffy/branch/master/graph/badge.svg)](https://codecov.io/gh/martimunicoy/peleffy) | | **Installation** | [![Conda](https://img.shields.io/conda/v/martimunicoy/peleffy.svg)](https://anaconda.org/martimunicoy/peleffy) [![PyPI](https://img.shields.io/pypi/v/peleffy)](https://pypi.org/project/peleffy/) | diff --git a/devtools/conda/meta.yaml b/devtools/conda/meta.yaml index 1064aded..9c044f0e 100644 --- a/devtools/conda/meta.yaml +++ b/devtools/conda/meta.yaml @@ -26,7 +26,8 @@ requirements: - networkx - rdkit - ambertools - - openff-toolkit==0.10.1 + - openff-toolkit==0.10.7 + - openff-forcefields==2023.05.1 about: home: https://github.com/martimunicoy/peleffy diff --git a/devtools/envs/standard.yaml b/devtools/envs/standard.yaml index 9b8f465a..7cf7bff7 100644 --- a/devtools/envs/standard.yaml +++ b/devtools/envs/standard.yaml @@ -9,8 +9,10 @@ dependencies: - pytest - pytest-cov - codecov - - coverage < 5.0 + - coverage - ipython - ambertools - rdkit - - openff-toolkit==0.10.1 + - openff-toolkit==0.10.7 + - openff-forcefields==2023.05.1 + diff --git a/docs/index.rst b/docs/index.rst index 4be57795..0fc5f77d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,6 +7,8 @@ are: * The following force field from the `Open Force Field toolkit `_: + * openff_unconstrained-2.1.0.offxml + * openff_unconstrained-2.0.0.offxml * openff_unconstrained-1.3.0.offxml diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index 2974aea4..7730f018 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -8,6 +8,24 @@ Releases follow the ``major.minor.micro`` scheme recommended by `PEP440 `_: support openff-toolkit 0.10.7 and openff-forcefield 2.1.0 +- `PR #182 `_: better alchemical structure alignment +- `PR #182 `_: more lambda types available (vdw1 and vdw2) + +Tests added +""""""""""" +- `PR #182 `_: corrections for alchemical and mapper tests +- `PR #182 `_: new test to check forcefield names + + 1.4.5 - Corrections for alchemistry ----------------------------------- diff --git a/peleffy/data/tests/alchemical_ligandParams_1.txt b/peleffy/data/tests/alchemical_ligandParams_1.txt index de3457e2..d67e49ce 100644 --- a/peleffy/data/tests/alchemical_ligandParams_1.txt +++ b/peleffy/data/tests/alchemical_ligandParams_1.txt @@ -9,24 +9,24 @@ }, "HYB": { "_C1_": { - "radius": 1.625, - "scale": 0.755 + "radius": 1.7, + "scale": 0.36 }, "_N1_": { "radius": 1.625, "scale": 0.755 }, "_C2_": { - "radius": 1.625, - "scale": 0.755 - }, - "_C3_": { "radius": 1.7, "scale": 0.36 }, + "_C3_": { + "radius": 1.625, + "scale": 0.755 + }, "_C4_": { - "radius": 1.7, - "scale": 0.36 + "radius": 1.625, + "scale": 0.755 }, "_C5_": { "radius": 1.7, @@ -41,20 +41,20 @@ "scale": 0.425 }, "_H1_": { - "radius": 1.25, - "scale": 0.85 + "radius": 1.2, + "scale": 0.425 }, "_H2_": { - "radius": 1.25, - "scale": 0.85 + "radius": 1.2, + "scale": 0.425 }, "_H3_": { "radius": 1.2, "scale": 0.425 }, "_H4_": { - "radius": 1.25, - "scale": 0.85 + "radius": 1.2, + "scale": 0.425 }, "_H5_": { "radius": 1.2, @@ -65,24 +65,24 @@ "scale": 0.425 }, "_H7_": { - "radius": 1.2, - "scale": 0.425 + "radius": 1.25, + "scale": 0.85 }, "_H8_": { "radius": 1.2, "scale": 0.425 }, "_H9_": { - "radius": 1.2, - "scale": 0.425 + "radius": 1.25, + "scale": 0.85 }, "H10_": { "radius": 1.2, "scale": 0.425 }, "H11_": { - "radius": 1.2, - "scale": 0.425 + "radius": 1.25, + "scale": 0.85 }, "_C6_": { "radius": 1.7, diff --git a/peleffy/data/tests/alchemical_ligandParams_2.txt b/peleffy/data/tests/alchemical_ligandParams_2.txt index cf36f765..c86cd18a 100644 --- a/peleffy/data/tests/alchemical_ligandParams_2.txt +++ b/peleffy/data/tests/alchemical_ligandParams_2.txt @@ -9,24 +9,24 @@ }, "HYB": { "_C1_": { - "radius": 1.55, - "scale": 0.79 + "radius": 1.7, + "scale": 0.0 }, "_N1_": { "radius": 1.7, "scale": 0.72 }, "_C2_": { - "radius": 1.55, - "scale": 0.79 - }, - "_C3_": { "radius": 1.7, "scale": 0.0 }, + "_C3_": { + "radius": 1.55, + "scale": 0.79 + }, "_C4_": { - "radius": 1.7, - "scale": 0.0 + "radius": 1.55, + "scale": 0.79 }, "_C5_": { "radius": 1.7, @@ -41,20 +41,20 @@ "scale": 0.0 }, "_H1_": { - "radius": 1.3, - "scale": 0.85 + "radius": 1.2, + "scale": 0.0 }, "_H2_": { - "radius": 1.3, - "scale": 0.85 + "radius": 1.2, + "scale": 0.0 }, "_H3_": { "radius": 1.2, "scale": 0.0 }, "_H4_": { - "radius": 1.3, - "scale": 0.85 + "radius": 1.2, + "scale": 0.0 }, "_H5_": { "radius": 1.2, @@ -65,24 +65,24 @@ "scale": 0.0 }, "_H7_": { - "radius": 1.2, - "scale": 0.0 + "radius": 1.3, + "scale": 0.85 }, "_H8_": { "radius": 1.2, "scale": 0.0 }, "_H9_": { - "radius": 1.2, - "scale": 0.0 + "radius": 1.3, + "scale": 0.85 }, "H10_": { "radius": 1.2, "scale": 0.0 }, "H11_": { - "radius": 1.2, - "scale": 0.0 + "radius": 1.3, + "scale": 0.85 }, "_C6_": { "radius": 1.7, diff --git a/peleffy/data/tests/alchemical_mol2.pdb b/peleffy/data/tests/alchemical_mol2.pdb index 086b5db7..98d6f6d7 100644 --- a/peleffy/data/tests/alchemical_mol2.pdb +++ b/peleffy/data/tests/alchemical_mol2.pdb @@ -1,21 +1,21 @@ COMPND HYDROLASE (SERINE PROTEINASE) -HETATM 1 C1 BEN A 1 -0.779 0.348 1.160 1.00 19.86 C -HETATM 2 C2 BEN A 1 -0.825 -0.520 2.218 1.00 19.86 C -HETATM 3 C3 BEN A 1 -0.996 -0.028 3.518 1.00 19.86 C -HETATM 4 C4 BEN A 1 -1.127 1.352 3.741 1.00 19.86 C -HETATM 5 C5 BEN A 1 -1.087 2.225 2.651 1.00 19.86 C -HETATM 6 C6 BEN A 1 -0.914 1.710 1.369 1.00 19.86 C -HETATM 7 C BEN A 1 -0.600 -0.141 -0.129 1.00 19.86 C -HETATM 8 N1 BEN A 1 -0.430 -1.431 -0.359 1.00 19.86 N -HETATM 9 N2 BEN A 1 -0.767 0.660 -1.167 1.00 19.86 N -HETATM 10 H2 BEN A 1 -0.729 -1.582 2.050 1.00 0.00 H -HETATM 11 H3 BEN A 1 -1.027 -0.710 4.355 1.00 0.00 H -HETATM 12 H4 BEN A 1 -1.256 1.727 4.746 1.00 0.00 H -HETATM 13 H5 BEN A 1 -1.191 3.289 2.807 1.00 0.00 H -HETATM 14 H6 BEN A 1 -0.885 2.390 0.531 1.00 0.00 H -HETATM 15 H1 1 BEN A 1 -0.302 -1.763 -1.304 1.00 0.00 H -HETATM 16 H1 2 BEN A 1 -0.429 -2.085 0.411 1.00 0.00 H -HETATM 17 HN2 BEN A 1 -0.637 0.307 -2.104 1.00 0.00 H +HETATM 1 C1 BEN A 1 0.918 -1.174 0.729 1.00 19.86 C +HETATM 2 C2 BEN A 1 0.521 -1.960 1.777 1.00 19.86 C +HETATM 3 C3 BEN A 1 1.179 -3.171 2.031 1.00 19.86 C +HETATM 4 C4 BEN A 1 2.249 -3.579 1.219 1.00 19.86 C +HETATM 5 C5 BEN A 1 2.648 -2.762 0.158 1.00 19.86 C +HETATM 6 C6 BEN A 1 1.975 -1.565 -0.074 1.00 19.86 C +HETATM 7 C BEN A 1 0.257 0.022 0.472 1.00 19.86 C +HETATM 8 N1 BEN A 1 -0.779 0.406 1.195 1.00 19.86 N +HETATM 9 N2 BEN A 1 0.764 0.882 -0.396 1.00 19.86 N +HETATM 10 H2 BEN A 1 -0.298 -1.647 2.407 1.00 0.00 H +HETATM 11 H3 BEN A 1 0.864 -3.794 2.854 1.00 0.00 H +HETATM 12 H4 BEN A 1 2.752 -4.514 1.417 1.00 0.00 H +HETATM 13 H5 BEN A 1 3.472 -3.061 -0.474 1.00 0.00 H +HETATM 14 H6 BEN A 1 2.291 -0.938 -0.896 1.00 0.00 H +HETATM 15 H1 1 BEN A 1 -1.243 1.279 0.989 1.00 0.00 H +HETATM 16 H1 2 BEN A 1 -1.108 -0.175 1.954 1.00 0.00 H +HETATM 17 HN2 BEN A 1 0.286 1.751 -0.584 1.00 0.00 H CONECT 1 2 2 6 7 CONECT 2 3 10 CONECT 3 4 4 11 diff --git a/peleffy/data/tests/alchemical_structure.pdb b/peleffy/data/tests/alchemical_structure.pdb index acb1c13a..5a85c943 100644 --- a/peleffy/data/tests/alchemical_structure.pdb +++ b/peleffy/data/tests/alchemical_structure.pdb @@ -17,17 +17,17 @@ HETATM 16 H8 HYB L 1 -1.108 -0.781 1.667 1.00 0.00 H HETATM 17 H9 HYB L 1 -0.323 0.825 1.846 1.00 0.00 H HETATM 18 H10 HYB L 1 1.237 0.709 -1.316 1.00 0.00 H HETATM 19 H11 HYB L 1 0.622 1.801 -0.063 1.00 0.00 H -HETATM 20 C6 HYB L 1 -0.779 0.348 1.160 1.00 0.00 C -HETATM 21 C7 HYB L 1 -0.825 -0.520 2.218 1.00 0.00 C -HETATM 22 C8 HYB L 1 -0.996 -0.028 3.518 1.00 0.00 C -HETATM 23 C9 HYB L 1 -1.127 1.352 3.741 1.00 0.00 C -HETATM 24 C10 HYB L 1 -1.087 2.225 2.651 1.00 0.00 C -HETATM 25 C11 HYB L 1 -0.914 1.710 1.369 1.00 0.00 C -HETATM 26 H12 HYB L 1 -0.729 -1.582 2.050 1.00 0.00 H -HETATM 27 H13 HYB L 1 -1.027 -0.710 4.355 1.00 0.00 H -HETATM 28 H14 HYB L 1 -1.256 1.727 4.746 1.00 0.00 H -HETATM 29 H15 HYB L 1 -1.191 3.289 2.807 1.00 0.00 H -HETATM 30 H16 HYB L 1 -0.885 2.390 0.531 1.00 0.00 H +HETATM 20 C6 HYB L 1 0.918 -1.174 0.729 1.00 0.00 C +HETATM 21 C7 HYB L 1 0.521 -1.960 1.777 1.00 0.00 C +HETATM 22 C8 HYB L 1 1.179 -3.171 2.031 1.00 0.00 C +HETATM 23 C9 HYB L 1 2.249 -3.579 1.219 1.00 0.00 C +HETATM 24 C10 HYB L 1 2.648 -2.762 0.158 1.00 0.00 C +HETATM 25 C11 HYB L 1 1.975 -1.565 -0.074 1.00 0.00 C +HETATM 26 H12 HYB L 1 -0.298 -1.647 2.407 1.00 0.00 H +HETATM 27 H13 HYB L 1 0.864 -3.794 2.854 1.00 0.00 H +HETATM 28 H14 HYB L 1 2.752 -4.514 1.417 1.00 0.00 H +HETATM 29 H15 HYB L 1 3.472 -3.061 -0.474 1.00 0.00 H +HETATM 30 H16 HYB L 1 2.291 -0.938 -0.896 1.00 0.00 H CONECT 1 2 9 10 11 CONECT 2 3 4 5 20 CONECT 3 12 13 14 diff --git a/peleffy/forcefield/selectors.py b/peleffy/forcefield/selectors.py index 0688ac12..e5185d84 100644 --- a/peleffy/forcefield/selectors.py +++ b/peleffy/forcefield/selectors.py @@ -12,7 +12,8 @@ class ForceFieldSelector(object): It defines a force field selector. """ _FF_TYPES = {'OPLS2005': ('OPLS2005', ), - 'OpenFF': ('openff_unconstrained-2.0.0.offxml', + 'OpenFF': ('openff_unconstrained-2.1.0.offxml', + 'openff_unconstrained-2.0.0.offxml', 'openff_unconstrained-1.3.0.offxml', 'openff_unconstrained-1.2.1.offxml', 'openff_unconstrained-1.2.0.offxml', diff --git a/peleffy/tests/test_alchemistry.py b/peleffy/tests/test_alchemistry.py index 747d7070..685222bc 100644 --- a/peleffy/tests/test_alchemistry.py +++ b/peleffy/tests/test_alchemistry.py @@ -117,7 +117,7 @@ def test_alchemizer_initialization_checker(self): None, 'C=C', 'C(Cl)(Cl)(Cl)', - [(0, 0), (1, 1), (2, 2), (3, 4)], + [(0, 0), (1, 2), (2, 3), (3, 4)], [6], [5], [6, 7, 8], @@ -176,7 +176,7 @@ def test_alchemizer_initialization_checker(self): 'ligands/ethylene.pdb', None, None, - [(0, 0), (1, 1), (3, 4), (2, 2)], + [(0, 0), (1, 1), (3, 5), (2, 2)], [4, 5], [3, 4], [2, 3, 4, 5], @@ -192,7 +192,7 @@ def test_alchemizer_initialization_checker(self): 'ligands/propionic_acid.pdb', None, None, - [(0, 5), (1, 0), (2, 6), (3, 1), (4, 2), + [(0, 6), (1, 0), (2, 7), (3, 1), (4, 2), (5, 4), (9, 10), (6, 3), (7, 8), (8, 9)], [10, ], [9, ], @@ -213,9 +213,9 @@ def test_alchemizer_initialization_checker(self): 'ligands/propionic_acid.pdb', None, None, - [(0, 0), (1, 1), (4, 2), (17, 3), (5, 4), - (6, 10), (2, 8), (3, 9), (8, 5), (9, 6), - (10, 7)], + [(3, 0), (1, 1), (4, 2), (18, 3), (5, 4), + (6, 10), (0, 8), (2, 9), (14, 5), + (15, 6), (16, 7)], [], [], [], @@ -223,10 +223,10 @@ def test_alchemizer_initialization_checker(self): 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67], [1], - [7, 11, 12, 13, 14, 15, 16, 18], - [7, 8, 9, 10, 11, 12, 15, 17], - [6, 7, 8, 9, 10, 11, 14, 19, 21, 22, 26, - 27, 28, 29, 30, 31, 32], + [7, 8, 9, 10, 11, 12, 13, 17], + [1, 2, 3, 7, 8, 9, 14, 17], + [3, 4, 5, 6, 7, 8, 13, 19, 20, 22, 23, + 24, 25, 26, 27, 28, 32], [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, @@ -238,8 +238,8 @@ def test_alchemizer_initialization_checker(self): 'ligands/benzamidine.pdb', None, None, - [(1, 6), (0, 7), (8, 14), (9, 15), (2, 8), - (11, 16)], + [(1, 6), (3, 7), (14, 14), (16, 15), + (4, 8), (18, 16)], [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], @@ -250,13 +250,13 @@ def test_alchemizer_initialization_checker(self): 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83], [1, 2, 3, 4, 5, 6, 7, 8], - [3, 4, 5, 6, 7, 10, 12, 13, 14, 15, 16, 17, - 18], - [3, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, - 17], - [1, 2, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, - 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, - 27, 28, 29, 30, 31, 32], + [0, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, + 15, 17], + [0, 1, 2, 3, 4, 7, 8, 9, 11, 13, 14, + 16, 17], + [0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 13, + 15, 16, 18, 19, 20, 21, 22, 23, 24, 25, + 26, 27, 28, 29, 31, 32], [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, @@ -3123,13 +3123,17 @@ def test_atoms_in_alchemical_topology(self, pdb1, pdb2, smiles1, smiles2, nonpolar_gammas.append(atom.nonpolar_gamma) nonpolar_alphas.append(atom.nonpolar_alpha) - assert sigmas == golden_sigmas, 'Unexpected sigmas' - assert epsilons == golden_epsilons, 'Unexpected epsilons' - assert born_radii == golden_born_radii, 'Unexpected born radii' - assert SASA_radii == golden_SASA_radii, 'Unexpected SASA radii' - assert nonpolar_gammas == golden_nonpolar_gammas, \ + assert sorted(sigmas) == sorted(golden_sigmas), \ + 'Unexpected sigmas' + assert sorted(epsilons) == sorted(golden_epsilons), \ + 'Unexpected epsilons' + assert sorted(born_radii) == sorted(golden_born_radii), \ + 'Unexpected born radii' + assert sorted(SASA_radii) == sorted(golden_SASA_radii), \ + 'Unexpected SASA radii' + assert sorted(nonpolar_gammas) == sorted(golden_nonpolar_gammas), \ 'Unexpected non polar gammas' - assert nonpolar_alphas == golden_nonpolar_alphas, \ + assert sorted(nonpolar_alphas) == sorted(golden_nonpolar_alphas), \ 'Unexpected non polar alphas' @pytest.mark.parametrize("pdb1, pdb2, smiles1, smiles2, " + @@ -3250,7 +3254,7 @@ def test_bonds_in_alchemical_topology(self, pdb1, pdb2, smiles1, bond = WritableBond(bond) bond_spring_constants.append(bond.spring_constant) - assert bond_spring_constants == golden_bond_spring_constants, \ + assert sorted(bond_spring_constants) == sorted(golden_bond_spring_constants), \ 'Unexpected spring constants' @pytest.mark.parametrize("pdb1, pdb2, smiles1, smiles2, " + @@ -3362,7 +3366,7 @@ def test_angles_in_alchemical_topology(self, pdb1, pdb2, smiles1, angle = WritableAngle(angle) angle_spring_constants.append(angle.spring_constant) - assert angle_spring_constants == golden_angle_spring_constants, \ + assert sorted(angle_spring_constants) == sorted(golden_angle_spring_constants), \ 'Unexpected spring constants' @pytest.mark.parametrize("pdb1, pdb2, smiles1, smiles2, " + diff --git a/peleffy/tests/test_forcefields.py b/peleffy/tests/test_forcefields.py index 97e6b0f3..1371b6fc 100644 --- a/peleffy/tests/test_forcefields.py +++ b/peleffy/tests/test_forcefields.py @@ -3,6 +3,8 @@ peleffy. """ +import pytest + class TestOpenForceField(object): """ @@ -11,14 +13,21 @@ class TestOpenForceField(object): FORCE_FIELD_NAME = 'openff_unconstrained-1.2.1.offxml' - def test_name(self): + @pytest.mark.parametrize("forcefield_name", + ['openff_unconstrained-1.2.1.offxml', + 'openff_unconstrained-2.0.0.offxml', + 'openff_unconstrained-2.1.0.offxml', + 'openff_unconstrained-2.1.1.offxml' + ]) + + def test_name(self, forcefield_name): """It checks the name assignment.""" from peleffy.forcefield import OpenForceField - openff = OpenForceField(self.FORCE_FIELD_NAME) + openff = OpenForceField(forcefield_name) - assert openff.name == self.FORCE_FIELD_NAME, \ + assert openff.name == forcefield_name, \ 'Unexpected force field name' def test_type(self): diff --git a/peleffy/tests/test_mapper.py b/peleffy/tests/test_mapper.py index 03685cca..2a78222b 100644 --- a/peleffy/tests/test_mapper.py +++ b/peleffy/tests/test_mapper.py @@ -57,8 +57,10 @@ def test_mapper_mapping(self): mapper = Mapper(mol1, mol2, include_hydrogens=False) mapping = mapper.get_mapping() - assert mapping == [(0, 1), (1, 2), (2, 0), (3, 6), - (4, 5), (5, 4), (6, 3)], 'Unexpected mapping' + assert (mapping == [(0, 1), (1, 2), (2, 0), (3, 6), (4, 5), (5, 4), (6, 3)] or + mapping == [(0, 6), (1, 7), (2, 5), (3, 4), (4, 3), (5, 1), (6, 0)] or + mapping == [(6, 1), (7, 2), (5, 0), (4, 6), (3, 5), (2, 4), (0, 3)] or + mapping == [(6, 6), (7, 7), (5, 5), (4, 4), (3, 3), (2, 1), (0, 0)]), 'Unexpected mapping' # Third mapping checker with hydrogens mol1 = Molecule(smiles='c1ccccc1', hydrogens_are_explicit=False) @@ -79,7 +81,7 @@ def test_mapper_mapping(self): mapper = Mapper(mol1, mol2, include_hydrogens=True) mapping = mapper.get_mapping() - assert mapping == [(0, 1), (1, 2), (8, 9), (9, 10), - (10, 11), (2, 0), (3, 6), (4, 5), - (5, 4), (6, 3), (7, 12), (14, 13), - (13, 14), (12, 7), (11, 8)], 'Unexpected mapping' + assert (mapping == [(0, 1), (1, 2), (8, 9), (9, 10), (10, 11), (2, 0), (3, 6), (4, 5), + (5, 4), (6, 3), (7, 12), (14, 13), (13, 14), (12, 7), (11, 8)] or + mapping == [(6, 1), (7, 2), (15, 9), (16, 10), (17, 11), (5, 0), (4, 6), (3, 5), + (2, 4), (0, 3), (1, 12), (11, 13), (12, 14), (13, 7), (14, 8)]), 'Unexpected mapping' diff --git a/peleffy/topology/alchemistry.py b/peleffy/topology/alchemistry.py index df07e2b8..7984824d 100644 --- a/peleffy/topology/alchemistry.py +++ b/peleffy/topology/alchemistry.py @@ -162,7 +162,8 @@ def connections(self): def get_alchemical_topology(self, fep_lambda=None, coul_lambda=None, coul1_lambda=None, coul2_lambda=None, - vdw_lambda=None, bonded_lambda=None): + vdw_lambda=None, vdw1_lambda=None, + vdw2_lambda=None, bonded_lambda=None): """ Given a lambda, it returns an alchemical topology after combining both input topologies. @@ -195,6 +196,20 @@ def get_alchemical_topology(self, fep_lambda=None, coul_lambda=None, affects van der Waals parameters. It needs to be contained between 0 and 1. It has precedence over fep_lambda. Default is None + vdw1_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 1. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 1. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None + vdw2_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 2. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 2. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None bonded_lambda : float The value to define a coulombic lambda. This lambda only affects bonded parameters. It needs to be contained @@ -212,10 +227,13 @@ def get_alchemical_topology(self, fep_lambda=None, coul_lambda=None, coul1_lambda = Coulombic1Lambda(coul1_lambda) coul2_lambda = Coulombic2Lambda(coul2_lambda) vdw_lambda = VanDerWaalsLambda(vdw_lambda) + vdw1_lambda = VanDerWaals1Lambda(vdw1_lambda) + vdw2_lambda = VanDerWaals2Lambda(vdw2_lambda) bonded_lambda = BondedLambda(bonded_lambda) lambda_set = LambdaSet(fep_lambda, coul_lambda, coul1_lambda, - coul2_lambda, vdw_lambda, bonded_lambda) + coul2_lambda, vdw_lambda, vdw1_lambda, + vdw2_lambda, bonded_lambda) alchemical_topology = self.topology_from_lambda_set(lambda_set) @@ -251,7 +269,7 @@ def topology_from_lambda_set(self, lambda_set): atom.apply_lambda(["sigma", "epsilon", "born_radius", "SASA_radius", "nonpolar_gamma", "nonpolar_alpha"], - lambda_set.get_lambda_for_vdw(), + lambda_set.get_lambda_for_vdw1(), reverse=False) atom.apply_lambda(["charge"], lambda_set.get_lambda_for_coulomb1(), @@ -261,7 +279,7 @@ def topology_from_lambda_set(self, lambda_set): atom.apply_lambda(["sigma", "epsilon", "born_radius", "SASA_radius", "nonpolar_gamma", "nonpolar_alpha"], - lambda_set.get_lambda_for_vdw(), + lambda_set.get_lambda_for_vdw2(), reverse=True) atom.apply_lambda(["charge"], lambda_set.get_lambda_for_coulomb2(), @@ -873,6 +891,7 @@ def molecule2_to_pdb(self, path): def rotamer_library_to_file(self, path, fep_lambda=None, coul_lambda=None, coul1_lambda=None, coul2_lambda=None, vdw_lambda=None, + vdw1_lambda=None, vdw2_lambda=None, bonded_lambda=None): """ It saves the alchemical rotamer library, which is the combination @@ -909,6 +928,20 @@ def rotamer_library_to_file(self, path, fep_lambda=None, affects van der Waals parameters. It needs to be contained between 0 and 1. It has precedence over fep_lambda. Default is None + vdw1_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 1. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 1. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None + vdw2_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 2. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 2. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None bonded_lambda : float The value to define a coulombic lambda. This lambda only affects bonded parameters. It needs to be contained @@ -927,14 +960,19 @@ def rotamer_library_to_file(self, path, fep_lambda=None, coul1_lambda = Coulombic1Lambda(coul1_lambda) coul2_lambda = Coulombic2Lambda(coul2_lambda) vdw_lambda = VanDerWaalsLambda(vdw_lambda) + vdw1_lambda = VanDerWaals1Lambda(vdw1_lambda) + vdw2_lambda = VanDerWaals2Lambda(vdw2_lambda) bonded_lambda = BondedLambda(bonded_lambda) lambda_set = LambdaSet(fep_lambda, coul_lambda, coul1_lambda, - coul2_lambda, vdw_lambda, bonded_lambda) + coul2_lambda, vdw_lambda, vdw1_lambda, + vdw2_lambda, bonded_lambda) if (at_least_one and lambda_set.get_lambda_for_bonded() == 0.0 and lambda_set.get_lambda_for_vdw() == 0.0 and + lambda_set.get_lambda_for_vdw1() == 0.0 and + lambda_set.get_lambda_for_vdw2() == 0.0 and lambda_set.get_lambda_for_coulomb() == 0.0 and lambda_set.get_lambda_for_coulomb1() == 0.0 and lambda_set.get_lambda_for_coulomb2() == 0.0): @@ -944,6 +982,8 @@ def rotamer_library_to_file(self, path, fep_lambda=None, elif (at_least_one and lambda_set.get_lambda_for_bonded() == 1.0 and lambda_set.get_lambda_for_vdw() == 1.0 and + lambda_set.get_lambda_for_vdw1() == 1.0 and + lambda_set.get_lambda_for_vdw2() == 1.0 and lambda_set.get_lambda_for_coulomb() == 1.0 and lambda_set.get_lambda_for_coulomb1() == 1.0 and lambda_set.get_lambda_for_coulomb2() == 1.0): @@ -964,6 +1004,8 @@ def rotamer_library_to_file(self, path, fep_lambda=None, for i, rotamer_branches in enumerate(rotamers): if i > 0: file.write(' newgrp &\n') + + #visited_atoms = set() for rotamer in rotamer_branches: index1 = rotamer.index1 index2 = rotamer.index2 @@ -974,12 +1016,19 @@ def rotamer_library_to_file(self, path, fep_lambda=None, atom_name1 = pdb_atom_names[index1] atom_name2 = pdb_atom_names[index2] + + # Separate alchemical branches + #if index1 in visited_atoms: + # file.write(' newgrp &\n') + #visited_atoms.add(index1) + file.write(' sidelib FREE{} {} {} &\n'.format( rotamer.resolution, atom_name1, atom_name2)) def obc_parameters_to_file(self, path, fep_lambda=None, coul_lambda=None, coul1_lambda=None, coul2_lambda=None, vdw_lambda=None, + vdw1_lambda=None, vdw2_lambda=None, bonded_lambda=None): """ It saves the alchemical OBC parameters, which is the combination @@ -1022,6 +1071,20 @@ def obc_parameters_to_file(self, path, fep_lambda=None, affects van der Waals parameters. It needs to be contained between 0 and 1. It has precedence over fep_lambda. Default is None + vdw1_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 1. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 1. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None + vdw2_lambda : float + The value to define a vdw lambda for exclusive atoms + of molecule 2. This lambda only affects van der Waals + parameters of exclusive atoms of molecule 2. It needs + to be contained between 0 and 1. + It has precedence over vdw_lambda and fep_lambda. + Default is None bonded_lambda : float The value to define a coulombic lambda. This lambda only affects bonded parameters. It needs to be contained @@ -1047,10 +1110,13 @@ def obc_parameters_to_file(self, path, fep_lambda=None, coul1_lambda = Coulombic1Lambda(coul1_lambda) coul2_lambda = Coulombic2Lambda(coul2_lambda) vdw_lambda = VanDerWaalsLambda(vdw_lambda) + vdw1_lambda = VanDerWaals1Lambda(vdw1_lambda) + vdw2_lambda = VanDerWaals2Lambda(vdw2_lambda) bonded_lambda = BondedLambda(bonded_lambda) lambda_set = LambdaSet(fep_lambda, coul_lambda, coul1_lambda, - coul2_lambda, vdw_lambda, bonded_lambda) + coul2_lambda, vdw_lambda, vdw1_lambda, + vdw2_lambda, bonded_lambda) # Define mappers mol1_mapped_atoms = [atom_pair[0] for atom_pair in self.mapping] @@ -1249,6 +1315,20 @@ class VanDerWaalsLambda(Lambda): _TYPE = "vdw" +class VanDerWaals1Lambda(Lambda): + """ + It defines the VanDerWaalsLambda1 class. It affects only van der Waals + parameters involving exclusive atoms of molecule 1. + """ + _TYPE = "vdw1" + +class VanDerWaals2Lambda(Lambda): + """ + It defines the VanDerWaalsLambda1 class. It affects only van der Waals + parameters involving exclusive atoms of molecule 2. + """ + _TYPE = "vdw2" + class BondedLambda(Lambda): """ It defines the BondedLambda class. It affects only bonded parameters. @@ -1262,7 +1342,7 @@ class LambdaSet(object): """ def __init__(self, fep_lambda, coul_lambda, coul1_lambda, coul2_lambda, - vdw_lambda, bonded_lambda): + vdw_lambda, vdw1_lambda, vdw2_lambda, bonded_lambda): """ It initializes a LambdaSet object which stores all the different types of lambda. @@ -1279,6 +1359,10 @@ def __init__(self, fep_lambda, coul_lambda, coul1_lambda, coul2_lambda, The coulombic lambda for exclusive atoms of molecule 2 vdw_lambda : a peleffy.topology.alchemy.VanDerWaalsLambda object The van der Waals lambda + vdw1_lambda : a peleffy.topology.alchemy.VanDerWaals1Lambda object + The van der Waals lambda for exclusive atoms of molecule 1 + vdw2_lambda : a peleffy.topology.alchemy.VanDerWaals2Lambda object + The van der Waals lambda for exclusive atoms of molecule 2 bonded_lambda : a peleffy.topology.alchemy.BondedLambda object The bonded lambda """ @@ -1300,6 +1384,12 @@ def __init__(self, fep_lambda, coul_lambda, coul1_lambda, coul2_lambda, if not isinstance(vdw_lambda, peleffy.topology.alchemistry.VanDerWaalsLambda): raise TypeError('Invalid vdw_lambda supplied to LambdaSet') + if not isinstance(vdw1_lambda, + peleffy.topology.alchemistry.VanDerWaals1Lambda): + raise TypeError('Invalid vdw_lambda1 supplied to LambdaSet') + if not isinstance(vdw2_lambda, + peleffy.topology.alchemistry.VanDerWaals2Lambda): + raise TypeError('Invalid vdw2_lambda supplied to LambdaSet') if not isinstance(bonded_lambda, peleffy.topology.alchemistry.BondedLambda): raise TypeError('Invalid bonded_lambda supplied to LambdaSet') @@ -1309,6 +1399,8 @@ def __init__(self, fep_lambda, coul_lambda, coul1_lambda, coul2_lambda, self._coul1_lambda = coul1_lambda self._coul2_lambda = coul2_lambda self._vdw_lambda = vdw_lambda + self._vdw1_lambda = vdw1_lambda + self._vdw2_lambda = vdw2_lambda self._bonded_lambda = bonded_lambda @property @@ -1371,6 +1463,30 @@ def vdw_lambda(self): """ return self._vdw_lambda + @property + def vdw1_lambda(self): + """ + It returns the vdw1_lambda value. + + Returns + ------- + vdw1_lambda : float + The value of the vdw1_lambda + """ + return self._vdw1_lambda + + @property + def vdw2_lambda(self): + """ + It returns the vdw2_lambda value. + + Returns + ------- + vdw2_lambda : float + The value of the vdw2_lambda + """ + return self._vdw2_lambda + @property def bonded_lambda(self): """ @@ -1385,12 +1501,14 @@ def bonded_lambda(self): def get_lambda_for_vdw(self): """ - It returns the lambda to be applied on van der Waals parameters. + It returns the lambda to be applied on van der Waals parameters of + both molecules. Returns ------- lambda_value : float - The lambda value to be applied on van der Waals parameters + The lambda value to be applied on van der Waals parameters of + both molecules """ if self.vdw_lambda.is_set: lambda_value = self.vdw_lambda.value @@ -1403,6 +1521,56 @@ def get_lambda_for_vdw(self): return lambda_value + def get_lambda_for_vdw1(self): + """ + It returns the lambda to be applied on van der Waals parameters of + exclusive atoms of molecule 1. + + Returns + ------- + lambda_value : float + The lambda value to be applied on van der Waals parameters of + exclusive atoms of molecule 1 + """ + if self.vdw1_lambda.is_set: + lambda_value = self.vdw1_lambda.value + + elif self.vdw_lambda.is_set: + lambda_value = self.vdw_lambda.value + + elif self.fep_lambda.is_set: + lambda_value = self.fep_lambda.value + + else: + lambda_value = 0.0 + + return lambda_value + + def get_lambda_for_vdw2(self): + """ + It returns the lambda to be applied on van der Waals parameters of + exclusive atoms of molecule 2. + + Returns + ------- + lambda_value : float + The lambda value to be applied on van der Waals parameters of + exclusive atoms of molecule 2 + """ + if self.vdw2_lambda.is_set: + lambda_value = self.vdw2_lambda.value + + elif self.vdw_lambda.is_set: + lambda_value = self.vdw_lambda.value + + elif self.fep_lambda.is_set: + lambda_value = self.fep_lambda.value + + else: + lambda_value = 0.0 + + return lambda_value + def get_lambda_for_coulomb(self): """ It returns the lambda to be applied on Coulomb parameters of diff --git a/peleffy/utils/toolkits.py b/peleffy/utils/toolkits.py index ea012aee..ab03f693 100644 --- a/peleffy/utils/toolkits.py +++ b/peleffy/utils/toolkits.py @@ -756,7 +756,9 @@ def get_atom_mapping(self, molecule1, molecule2, mcs_mol, The list of atom pairs between both molecules, represented with tuples """ + import itertools from rdkit.Chem import AllChem + from rdkit.Chem.rdMolAlign import CalcRMS rdkit_mol1 = deepcopy(molecule1.rdkit_molecule) rdkit_mol2 = deepcopy(molecule2.rdkit_molecule) @@ -768,13 +770,13 @@ def get_atom_mapping(self, molecule1, molecule2, mcs_mol, # Map atoms between mol1 and MCS mol if rdkit_mol1.HasSubstructMatch(mcs_mol): - mol1_sub = rdkit_mol1.GetSubstructMatch(mcs_mol) + mol1_subs = rdkit_mol1.GetSubstructMatches(mcs_mol) else: raise ValueError('RDKit MCS Subgraph molecule 1 search failed') # Map atoms between mol2 and MCS mol if rdkit_mol2.HasSubstructMatch(mcs_mol): - mol2_sub = rdkit_mol2.GetSubstructMatch(mcs_mol) + mol2_subs = rdkit_mol2.GetSubstructMatches(mcs_mol) else: raise ValueError('RDKit MCS Subgraph molecule 2 search failed') @@ -785,8 +787,26 @@ def get_atom_mapping(self, molecule1, molecule2, mcs_mol, raise ValueError('RDKit MCS Subgraph search failed') """ + # Find best mapping + best_mol1_sub = None + best_mol2_sub = None + best_rmsd = None + for mol1_sub, mol2_sub in itertools.product(mol1_subs, mol2_subs): + rmsd = CalcRMS(rdkit_mol1, rdkit_mol2, + map=[list(zip(mol1_sub, mol2_sub)), ]) + + if best_rmsd is None: + best_mol1_sub = mol1_sub + best_mol2_sub = mol2_sub + best_rmsd = rmsd + + elif rmsd < best_rmsd: + best_mol1_sub = mol1_sub + best_mol2_sub = mol2_sub + best_rmsd = rmsd + # Map between the two molecules - mapping = list(zip(mol1_sub, mol2_sub)) + mapping = list(zip(best_mol1_sub, best_mol2_sub)) return mapping